DILIKernel.h

Go to the documentation of this file.

#ifndef DILIKERNEL_H_
#define DILIKERNEL_H_
 
#include "MUQ/SamplingAlgorithms/TransitionKernel.h"
 
#include "MUQ/SamplingAlgorithms/MCMCProposal.h"
#include "MUQ/Modeling/Distributions/GaussianBase.h"
#include "MUQ/Modeling/LinearAlgebra/LinearOperator.h"
 
namespace muq {
  namespace SamplingAlgorithms {
 
    class AverageHessian : public muq::Modeling::LinearOperator
    {
    public:
      AverageHessian(unsigned int                            numOldSamps,
                     std::shared_ptr<Eigen::ColPivHouseholderQR<Eigen::MatrixXd>> const& uQRIn,
                     std::shared_ptr<Eigen::MatrixXd> const& oldWIn,
                     std::shared_ptr<Eigen::VectorXd> const& oldValsIn,
                     std::shared_ptr<muq::Modeling::LinearOperator> const& newHess);
 
      virtual Eigen::MatrixXd Apply(Eigen::Ref<const Eigen::MatrixXd> const& x) override;
 
      virtual Eigen::MatrixXd ApplyTranspose(Eigen::Ref<const Eigen::MatrixXd> const& x) override;
 
    private:
      const double numSamps;
      std::shared_ptr<Eigen::MatrixXd> oldW;
      std::shared_ptr<Eigen::ColPivHouseholderQR<Eigen::MatrixXd>> uQR; // holds the QR decomp of U
 
      std::shared_ptr<Eigen::VectorXd> oldEigVals;
      std::shared_ptr<muq::Modeling::LinearOperator> newHess;
      Eigen::MatrixXd thinQ;
 
    };
 
    class CSProjector : public muq::Modeling::LinearOperator
    {
    public:
 
      CSProjector(std::shared_ptr<Eigen::MatrixXd> const& Uin,
                  std::shared_ptr<Eigen::MatrixXd> const& Win,
                  unsigned int                            lisDimIn) : LinearOperator(Uin->rows(), Win->rows()),
                                                                 U(Uin), W(Win), lisDim(lisDimIn){};
 
      virtual ~CSProjector() = default;
 
      virtual Eigen::MatrixXd Apply(Eigen::Ref<const Eigen::MatrixXd> const& x) override;
 
      virtual Eigen::MatrixXd ApplyTranspose(Eigen::Ref<const Eigen::MatrixXd> const& x) override;
 
    private:
      std::shared_ptr<Eigen::MatrixXd> U, W;
      unsigned int lisDim;
    };
 
 
    class LIS2Full : public muq::Modeling::LinearOperator
    {
    public:
      LIS2Full(std::shared_ptr<Eigen::MatrixXd> const& Uin,
               std::shared_ptr<Eigen::VectorXd> const& Lin) : LinearOperator(Uin->rows(), Lin->rows()),
                                                              U(Uin), L(Lin), lisDim(Lin->rows()){};
 
      virtual ~LIS2Full() = default;
 
      virtual Eigen::MatrixXd Apply(Eigen::Ref<const Eigen::MatrixXd> const& x) override;
 
      virtual Eigen::MatrixXd ApplyTranspose(Eigen::Ref<const Eigen::MatrixXd> const& x) override;
 
    private:
      std::shared_ptr<Eigen::MatrixXd> U;
      std::shared_ptr<Eigen::VectorXd> L;
      const unsigned int lisDim;
    };
 
 
    class DILIKernel : public TransitionKernel {
    public:
 
      DILIKernel(boost::property_tree::ptree       const& pt,
                 std::shared_ptr<AbstractSamplingProblem> problem);
 
      DILIKernel(boost::property_tree::ptree       const& pt,
                 std::shared_ptr<AbstractSamplingProblem> problem,
                 Eigen::VectorXd                   const& genEigVals,
                 Eigen::MatrixXd                   const& genEigVecs);
 
      DILIKernel(boost::property_tree::ptree                  const& pt,
                 std::shared_ptr<AbstractSamplingProblem>            problem,
                 std::shared_ptr<muq::Modeling::GaussianBase> const& prior,
                 std::shared_ptr<muq::Modeling::ModPiece>     const& noiseModel,
                 std::shared_ptr<muq::Modeling::ModPiece>     const& forwardModelIn);
 
      DILIKernel(boost::property_tree::ptree                  const& pt,
                 std::shared_ptr<AbstractSamplingProblem>            problem,
                 std::shared_ptr<muq::Modeling::GaussianBase> const& prior,
                 std::shared_ptr<muq::Modeling::ModPiece>     const& noiseModel,
                 std::shared_ptr<muq::Modeling::ModPiece>     const& forwardModelIn,
                 Eigen::VectorXd                              const& genEigVals,
                 Eigen::MatrixXd                              const& genEigVecs);
 
 
      DILIKernel(boost::property_tree::ptree                  const& pt,
                 std::shared_ptr<AbstractSamplingProblem>            problem,
                 std::shared_ptr<muq::Modeling::GaussianBase> const& prior,
                 std::shared_ptr<muq::Modeling::ModPiece>     const& likelihood);
 
      DILIKernel(boost::property_tree::ptree                  const& pt,
                 std::shared_ptr<AbstractSamplingProblem>            problem,
                 std::shared_ptr<muq::Modeling::GaussianBase> const& prior,
                 std::shared_ptr<muq::Modeling::ModPiece>     const& likelihood,
                 Eigen::VectorXd                              const& genEigVals,
                 Eigen::MatrixXd                              const& genEigVecs);
 
      virtual ~DILIKernel() = default;
 
      virtual inline std::shared_ptr<TransitionKernel> LISKernel() {return lisKernel;};
      virtual inline std::shared_ptr<TransitionKernel> CSKernel() {return csKernel;};
 
      virtual void PostStep(unsigned int const t,
                            std::vector<std::shared_ptr<SamplingState>> const& state) override;
 
      virtual std::vector<std::shared_ptr<SamplingState>> Step(unsigned int const t,
                                                               std::shared_ptr<SamplingState> prevState) override;
 
      virtual void PrintStatus(std::string prefix) const override;
 
      static std::shared_ptr<muq::Modeling::ModPiece> ExtractLikelihood(std::shared_ptr<AbstractSamplingProblem> const& problem,
                                                                        std::string                              const& nodeName);
 
 
      static std::shared_ptr<muq::Modeling::GaussianBase> ExtractPrior(std::shared_ptr<AbstractSamplingProblem> const& problem,
                                                                       std::string                              const& nodeName);
 
      static std::shared_ptr<muq::Modeling::ModPiece> ExtractNoiseModel(std::shared_ptr<muq::Modeling::ModPiece> const& likelihood);
 
      // Extract the forward model from the likelihood.  This function assumes the graph has only one node between the output of the likelihood function and the output of the forward model
      static std::shared_ptr<muq::Modeling::ModPiece> ExtractForwardModel(std::shared_ptr<muq::Modeling::ModPiece> const& likelihoodIn);
 
      static std::shared_ptr<muq::Modeling::ModPiece> CreateLikelihood(std::shared_ptr<muq::Modeling::ModPiece> const& forwardModel,
                                                                       std::shared_ptr<muq::Modeling::ModPiece> const& noiseDensity);
 
      Eigen::MatrixXd const& LISVecs() const{return *hessU;};
      Eigen::VectorXd const& LISVals() const{return *hessEigVals;};
      Eigen::MatrixXd const& LISW() const{return *hessW;};
      Eigen::VectorXd const& LISL() const{return *lisL;};
 
      unsigned int LISDim() const{return lisDim;};
 
      void CreateLIS(std::vector<Eigen::VectorXd> const& currState);
 
      Eigen::VectorXd ToLIS(Eigen::VectorXd const& x) const;
 
      Eigen::VectorXd FromLIS(Eigen::VectorXd const& r) const;
 
      Eigen::VectorXd ToCS(Eigen::VectorXd const& x) const;
 
 
    protected:
 
      std::pair<Eigen::MatrixXd, Eigen::VectorXd> ComputeLocalLIS(std::vector<Eigen::VectorXd> const& currState);
 
 
      void SetLIS(Eigen::VectorXd const& eigVals, Eigen::MatrixXd const& eigVecs);
 
 
 
      void UpdateLIS(unsigned int                        numSamps,
                     std::vector<Eigen::VectorXd> const& currState);
 
      // Used to reconstruct the LIS and CS kernels when lisU, lisW, or lisL changes.
      void UpdateKernels();
 
      boost::property_tree::ptree lisKernelOpts;
      boost::property_tree::ptree csKernelOpts;
 
      std::shared_ptr<muq::Modeling::ModPiece> logLikelihood;
      std::shared_ptr<muq::Modeling::GaussianBase> prior;
 
      std::shared_ptr<muq::Modeling::ModPiece> forwardModel;
      std::shared_ptr<muq::Modeling::ModPiece> noiseDensity;
 
      // A matrix containing the eigenvectors of the generalized eigenvalue problem Hv = lam*\Gamma^{-1}v
      std::shared_ptr<Eigen::MatrixXd> hessU;
 
      // The following matrices hold a QR decomposition of hessU
      std::shared_ptr<Eigen::ColPivHouseholderQR<Eigen::MatrixXd>> hessUQR;
 
      // A vector of eigenvalues of the generalized eigenvalue problem.  Only values above hessValTol are kept.
      std::shared_ptr<Eigen::VectorXd> hessEigVals;
 
      // W = \Gamma^{-1} v
      std::shared_ptr<Eigen::MatrixXd> hessW;
 
      // L is the Cholesky factor of the approximate posterior covariance projected onto the LIS
      std::shared_ptr<Eigen::VectorXd> lisL;
 
      // Defines the operation from the LIS to the full space
      std::shared_ptr<muq::Modeling::LinearOperator> lisToFull;
 
      // Defines the projection onto the complementary space
      std::shared_ptr<muq::Modeling::LinearOperator> fullToCS;
 
      // Transition kernel for updating the LIS
      std::shared_ptr<TransitionKernel> lisKernel;
 
      // Transition kernel for updating the CS
      std::shared_ptr<TransitionKernel> csKernel;
 
      // The type of Hessian approximation to use.  "Exact" or "GaussNewton"
      std::string hessType;
 
      // Options for the Generalized eigenvalue solver
      boost::property_tree::ptree eigOpts;
 
      const int updateInterval;
      const int adaptStart;
      const int adaptEnd;
      const int initialHessSamps;
 
      /* The size of the LIS.  Note that the LIS can have a smaller dimension
       than the number of eigenvalue/eigenvector pairs kept to approximate the
       global hessian.  The number of eigenvalues greater than lisValTol dictates
       the LIS dimension
      */
      unsigned int lisDim=0;
 
      /* Threshold on eigenvalues defining the low rank approximation of the global
        hessian.  Note that this toleran is typically smaller than the tolerance
        used to define the LIS.  In Cui et al., this was set to 1e-4.
      */
      const double hessValTol;
 
      /* Threshold on eigenvalues used to define local and global LIS
         approximations.  Modes with eigenavalues greater than this threshold
         will be used.  In Cui et al., this was set to 0.1
      */
      const double lisValTol;
 
      unsigned int numLisUpdates;
    };
  } // namespace SamplingAlgorithms
} // namespace muq
 
#endif