lsst::ip::diffim::RegularizedKernelSolution< InputT > Class Template Reference

#include <KernelSolution.h>

Inheritance diagram for lsst::ip::diffim::RegularizedKernelSolution< InputT >:

Public Types
typedef boost::shared_ptr < RegularizedKernelSolution < InputT > >	Ptr
Public Types inherited from lsst::ip::diffim::StaticKernelSolution< InputT >
typedef boost::shared_ptr < StaticKernelSolution< InputT > >	Ptr
Public Types inherited from lsst::ip::diffim::KernelSolution
enum	KernelSolvedBy {   NONE = 0, CHOLESKY_LDLT = 1, CHOLESKY_LLT = 2, LU = 3,   EIGENVECTOR = 4 }
enum	ConditionNumberType { EIGENVALUE = 0, SVD = 1 }
typedef boost::shared_ptr < KernelSolution >	Ptr
typedef lsst::afw::math::Kernel::Pixel	PixelT
typedef lsst::afw::image::Image < lsst::afw::math::Kernel::Pixel >	ImageT

Public Member Functions
	RegularizedKernelSolution (lsst::afw::math::KernelList const &basisList, bool fitForBackground, boost::shared_ptr< Eigen::MatrixXd > hMat, lsst::pex::policy::Policy policy)
virtual	~RegularizedKernelSolution ()
void	solve ()
double	getLambda ()
double	estimateRisk (double maxCond)
boost::shared_ptr < Eigen::MatrixXd >	getM (bool includeHmat=true)
Public Member Functions inherited from lsst::ip::diffim::StaticKernelSolution< InputT >
	StaticKernelSolution (lsst::afw::math::KernelList const &basisList, bool fitForBackground)
virtual	~StaticKernelSolution ()
virtual void	build (lsst::afw::image::Image< InputT > const &templateImage, lsst::afw::image::Image< InputT > const &scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &varianceEstimate)
virtual lsst::afw::math::Kernel::Ptr	getKernel ()
virtual lsst::afw::image::Image < lsst::afw::math::Kernel::Pixel > ::Ptr	makeKernelImage ()
virtual double	getBackground ()
virtual double	getKsum ()
virtual std::pair < boost::shared_ptr < lsst::afw::math::Kernel > , double >	getSolutionPair ()
Public Member Functions inherited from lsst::ip::diffim::KernelSolution
	KernelSolution (boost::shared_ptr< Eigen::MatrixXd > mMat, boost::shared_ptr< Eigen::VectorXd > bVec, bool fitForBackground)
	KernelSolution (bool fitForBackground)
	KernelSolution ()
virtual	~KernelSolution ()
virtual void	solve (Eigen::MatrixXd mMat, Eigen::VectorXd bVec)
KernelSolvedBy	getSolvedBy ()
virtual double	getConditionNumber (ConditionNumberType conditionType)
virtual double	getConditionNumber (Eigen::MatrixXd mMat, ConditionNumberType conditionType)
boost::shared_ptr < Eigen::MatrixXd >	getM ()
boost::shared_ptr < Eigen::VectorXd >	getB ()
void	printM ()
void	printB ()
void	printA ()
int	getId () const

Private Member Functions
std::vector< double >	_createLambdaSteps ()

Private Attributes
boost::shared_ptr < Eigen::MatrixXd >	_hMat
	Regularization weights. More...
double	_lambda
	Overall regularization strength. More...
lsst::pex::policy::Policy	_policy

Additional Inherited Members
Protected Member Functions inherited from lsst::ip::diffim::StaticKernelSolution< InputT >
void	_setKernel ()
	Set kernel after solution. More...
void	_setKernelUncertainty ()
	Not implemented. More...
Protected Attributes inherited from lsst::ip::diffim::StaticKernelSolution< InputT >
boost::shared_ptr < Eigen::MatrixXd >	_cMat
	K_i x R. More...
boost::shared_ptr < Eigen::VectorXd >	_iVec
	Vectorized I. More...
boost::shared_ptr < Eigen::VectorXd >	_ivVec
	Inverse variance. More...
lsst::afw::math::Kernel::Ptr	_kernel
	Derived single-object convolution kernel. More...
double	_background
	Derived differential background estimate. More...
double	_kSum
	Derived kernel sum. More...
Protected Attributes inherited from lsst::ip::diffim::KernelSolution
int	_id
	Unique ID for object. More...
boost::shared_ptr < Eigen::MatrixXd >	_mMat
	Derived least squares M matrix. More...
boost::shared_ptr < Eigen::VectorXd >	_bVec
	Derived least squares B vector. More...
boost::shared_ptr < Eigen::VectorXd >	_aVec
	Derived least squares solution matrix. More...
KernelSolvedBy	_solvedBy
	Type of algorithm used to make solution. More...
bool	_fitForBackground
	Background terms included in fit. More...
Static Protected Attributes inherited from lsst::ip::diffim::KernelSolution
static int	_SolutionId = 0
	Unique identifier for solution. More...

Detailed Description

template<typename InputT>
class lsst::ip::diffim::RegularizedKernelSolution< InputT >

Definition at line 147 of file KernelSolution.h.

Member Typedef Documentation

template<typename InputT >

typedef boost::shared_ptr<RegularizedKernelSolution<InputT> > lsst::ip::diffim::RegularizedKernelSolution< InputT >::Ptr

Definition at line 149 of file KernelSolution.h.

Constructor & Destructor Documentation

template<typename InputT >

lsst::ip::diffim::RegularizedKernelSolution< InputT >::RegularizedKernelSolution	(	lsst::afw::math::KernelList const &	basisList,
		bool	fitForBackground,
		boost::shared_ptr< Eigen::MatrixXd >	hMat,
		lsst::pex::policy::Policy	policy
	)

Definition at line 1103 of file KernelSolution.cc.

        :
        StaticKernelSolution<InputT>(basisList, fitForBackground),
        _hMat(hMat),
        _policy(policy)
    {};

template<typename InputT >

virtual lsst::ip::diffim::RegularizedKernelSolution< InputT >::~RegularizedKernelSolution ( )

inlinevirtual

Definition at line 156 of file KernelSolution.h.

{};

Member Function Documentation

template<typename InputT >

std::vector< double > lsst::ip::diffim::RegularizedKernelSolution< InputT >::_createLambdaSteps ( )

private

Definition at line 1305 of file KernelSolution.cc.

                                                                            {
        std::vector<double> lambdas;

        std::string lambdaStepType = _policy.getString("lambdaStepType");    
        if (lambdaStepType == "linear") {
            double lambdaLinMin   = _policy.getDouble("lambdaLinMin");        
            double lambdaLinMax   = _policy.getDouble("lambdaLinMax");
            double lambdaLinStep  = _policy.getDouble("lambdaLinStep");
            for (double l = lambdaLinMin; l <= lambdaLinMax; l += lambdaLinStep) {
                lambdas.push_back(l);
            }
        }
        else if (lambdaStepType == "log") {
            double lambdaLogMin   = _policy.getDouble("lambdaLogMin");        
            double lambdaLogMax   = _policy.getDouble("lambdaLogMax");
            double lambdaLogStep  = _policy.getDouble("lambdaLogStep");
            for (double l = lambdaLogMin; l <= lambdaLogMax; l += lambdaLogStep) {
                lambdas.push_back(pow(10, l));
            }
        }
        else {
            throw LSST_EXCEPT(pexExcept::Exception, "lambdaStepType in Policy not recognized");
        }
        return lambdas;
    }

template<typename InputT >

double lsst::ip::diffim::RegularizedKernelSolution< InputT >::estimateRisk

(

double

maxCond

)

Definition at line 1116 of file KernelSolution.cc.

                                                                         {
        Eigen::MatrixXd vMat      = (this->_cMat)->jacobiSvd().matrixV();
        Eigen::MatrixXd vMatvMatT = vMat * vMat.transpose();

        /* Find pseudo inverse of mMat, which may be ill conditioned */
        Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(*(this->_mMat));
        Eigen::MatrixXd const& rMat = eVecValues.eigenvectors();
        Eigen::VectorXd eValues = eVecValues.eigenvalues();
        double eMax = eValues.maxCoeff();
        for (int i = 0; i != eValues.rows(); ++i) {
            if (eValues(i) == 0.0) {
                eValues(i) = 0.0;
            }
            else if ((eMax / eValues(i)) > maxCond) {
                pexLog::TTrace<5>("lsst.ip.diffim.RegularizedKernelSolution.estimateRisk", 
                                  "Truncating eValue %d; %.5e / %.5e = %.5e vs. %.5e",
                                  i, eMax, eValues(i), eMax / eValues(i), maxCond);
                eValues(i) = 0.0;
            }
            else {
                eValues(i) = 1.0 / eValues(i);
            }
        }
        Eigen::MatrixXd mInv    = rMat * eValues.asDiagonal() * rMat.transpose();

        std::vector<double> lambdas = _createLambdaSteps();
        std::vector<double> risks;
        for (unsigned int i = 0; i < lambdas.size(); i++) {
            double l = lambdas[i];
            Eigen::MatrixXd mLambda = *(this->_mMat) + l * (*_hMat);
            
            try {
                KernelSolution::solve(mLambda, *(this->_bVec));
            } catch (pexExcept::Exception &e) {
                LSST_EXCEPT_ADD(e, "Unable to solve regularized kernel matrix");
                throw e;
            }
            Eigen::VectorXd term1 = (this->_aVec->transpose() * vMatvMatT * *(this->_aVec));
            if (term1.size() != 1)
                throw LSST_EXCEPT(pexExcept::Exception, "Matrix size mismatch");

            double term2a = (vMatvMatT * mLambda.inverse()).trace();

            Eigen::VectorXd term2b = (this->_aVec->transpose() * (mInv * *(this->_bVec)));
            if (term2b.size() != 1)
                throw LSST_EXCEPT(pexExcept::Exception, "Matrix size mismatch");

            double risk   = term1(0) + 2 * (term2a - term2b(0));
            pexLog::TTrace<6>("lsst.ip.diffim.RegularizedKernelSolution.estimateRisk", 
                              "Lambda = %.3f, Risk = %.5e", 
                              l, risk);
            pexLog::TTrace<7>("lsst.ip.diffim.RegularizedKernelSolution.estimateRisk", 
                              "%.5e + 2 * (%.5e - %.5e)", 
                              term1(0), term2a, term2b(0));
            risks.push_back(risk);
        }
        std::vector<double>::iterator it = min_element(risks.begin(), risks.end());
        int index = distance(risks.begin(), it);
        pexLog::TTrace<5>("lsst.ip.diffim.RegularizedKernelSolution.estimateRisk", 
                          "Minimum Risk = %.3e at lambda = %.3e", risks[index], lambdas[index]);

        return lambdas[index];
        
    }

template<typename InputT >

double lsst::ip::diffim::RegularizedKernelSolution< InputT >::getLambda ( )

inline

Definition at line 158 of file KernelSolution.h.

{return _lambda;}

template<typename InputT >

boost::shared_ptr< Eigen::MatrixXd > lsst::ip::diffim::RegularizedKernelSolution< InputT >::getM

(

bool

includeHmat = true

)

Definition at line 1182 of file KernelSolution.cc.

                                                                                           {
        if (includeHmat == true) {
            return (boost::shared_ptr<Eigen::MatrixXd>(
                        new Eigen::MatrixXd(*(this->_mMat) + _lambda * (*_hMat))
                        ));
        }
        else {
            return this->_mMat;
        }
    }

template<typename InputT >

void lsst::ip::diffim::RegularizedKernelSolution< InputT >::solve ( )

virtual

Reimplemented from lsst::ip::diffim::StaticKernelSolution< InputT >.

Definition at line 1194 of file KernelSolution.cc.

                                                  {

        pexLog::TTrace<5>("lsst.ip.diffim.RegularizedKernelSolution.solve", 
                          "cMat is %d x %d; vVec is %d; iVec is %d; hMat is %d x %d", 
                          (*this->_cMat).rows(), (*this->_cMat).cols(), (*this->_ivVec).size(), 
                          (*this->_iVec).size(), (*_hMat).rows(), (*_hMat).cols());

        if (DEBUG_MATRIX2) {
            std::cout << "ID: " << (this->_id) << std::endl;
            std::cout << "C:" << std::endl;
            std::cout << (*this->_cMat) << std::endl;
            std::cout << "Sigma^{-1}:" << std::endl;
            std::cout << Eigen::MatrixXd((*this->_ivVec).asDiagonal()) << std::endl;
            std::cout << "Y:" << std::endl;
            std::cout << (*this->_iVec) << std::endl;
            std::cout << "H:" << std::endl;
            std::cout << (*_hMat) << std::endl;
        }


        this->_mMat.reset(
            new Eigen::MatrixXd(this->_cMat->transpose() * this->_ivVec->asDiagonal() * *(this->_cMat))
            );
        this->_bVec.reset(
            new Eigen::VectorXd(this->_cMat->transpose() * this->_ivVec->asDiagonal() * *(this->_iVec))
            );
        
        
        /* See N.R. 18.5
           
           Matrix equation to solve is Y = C a + N
           Y   = vectorized version of I (I = image to not convolve)
           C_i = K_i (x) R (R = image to convolve)
           a   = kernel coefficients
           N   = noise
           
           If we reweight everything by the inverse square root of the noise
           covariance, we get a linear model with the identity matrix for
           the noise.  The problem can then be solved using least squares,
           with the normal equations
           
              C^T Y = C^T C a

           or 
           
              b = M a

           with 

              b = C^T Y
              M = C^T C
              a = (C^T C)^{-1} C^T Y


           We can regularize the least square problem
  
              Y = C a + lambda a^T H a       (+ N, which can be ignored)

           or the normal equations

              (C^T C + lambda H) a = C^T Y
 
              
           The solution to the regularization of the least squares problem is

              a = (C^T C + lambda H)^{-1} C^T Y

           The approximation to Y is 

              C a = C (C^T C + lambda H)^{-1} C^T Y
 
           with smoothing matrix 

              S = C (C^T C + lambda H)^{-1} C^T 

        */
        
        std::string lambdaType = _policy.getString("lambdaType");        
        if (lambdaType == "absolute") {
            _lambda = _policy.getDouble("lambdaValue");
        }
        else if (lambdaType ==  "relative") {
            _lambda  = this->_mMat->trace() / this->_hMat->trace();
            _lambda *= _policy.getDouble("lambdaScaling");
        }
        else if (lambdaType ==  "minimizeBiasedRisk") {
            double tol = _policy.getDouble("maxConditionNumber");
            _lambda = estimateRisk(tol);
        }
        else if (lambdaType ==  "minimizeUnbiasedRisk") {
            _lambda = estimateRisk(std::numeric_limits<double>::max());
        }
        else {
            throw LSST_EXCEPT(pexExcept::Exception, "lambdaType in Policy not recognized");
        }
        
        pexLog::TTrace<5>("lsst.ip.diffim.KernelCandidate.build", 
                          "Applying kernel regularization with lambda = %.2e", _lambda);
        
        
        try {
            KernelSolution::solve(*(this->_mMat) + _lambda * (*_hMat), *(this->_bVec));
        } catch (pexExcept::Exception &e) {
            LSST_EXCEPT_ADD(e, "Unable to solve static kernel matrix");
            throw e;
        }
        /* Turn matrices into _kernel and _background */
        StaticKernelSolution<InputT>::_setKernel();
    }

Member Data Documentation

template<typename InputT >

boost::shared_ptr<Eigen::MatrixXd> lsst::ip::diffim::RegularizedKernelSolution< InputT >::_hMat

private

Regularization weights.

Definition at line 165 of file KernelSolution.h.

template<typename InputT >

double lsst::ip::diffim::RegularizedKernelSolution< InputT >::_lambda

private

Overall regularization strength.

Definition at line 166 of file KernelSolution.h.

template<typename InputT >

lsst::pex::policy::Policy lsst::ip::diffim::RegularizedKernelSolution< InputT >::_policy

private

Definition at line 167 of file KernelSolution.h.

---

The documentation for this class was generated from the following files:

• /lsst/home/lsstsw/stack/Linux64/ip_diffim/master-gcec1da163f+63/include/lsst/ip/diffim/KernelSolution.h
• /lsst/home/lsstsw/stack/Linux64/ip_diffim/master-gcec1da163f+63/src/KernelSolution.cc