lsst::afw::math::LeastSquares Class Reference

Solver for linear least-squares problems. More...

#include <LeastSquares.h>

Classes
class	Impl

Public Types
enum	Factorization { NORMAL_EIGENSYSTEM, NORMAL_CHOLESKY, DIRECT_SVD }
	Private implementation; forward-declared publicly so we can inherit from it in .cc. More...

Public Member Functions
template<typename T1 , typename T2 , int C1, int C2>
void	setDesignMatrix (ndarray::Array< T1, 2, C1 > const &design, ndarray::Array< T2, 1, C2 > const &data)
	Reset the design matrix and data vector given as ndarrays; dimension must not change. More...
template<typename D1 , typename D2 >
void	setDesignMatrix (Eigen::MatrixBase< D1 > const &design, Eigen::MatrixBase< D2 > const &data)
	Reset the design matrix and data vector given as Eigen objects; dimension must not change. More...
template<typename T1 , int C1>
void	setDesignMatrix (ndarray::Array< T1, 2, C1 > const &design)
	Reset the design matrix given as an ndarray; dimension and data are not changed. More...
template<typename D1 , typename D2 >
void	setDesignMatrix (Eigen::MatrixBase< D1 > const &design)
	Reset the design matrix given as an Eigen object; dimension and data are not changed. More...
template<typename T1 , typename T2 , int C1, int C2>
void	setNormalEquations (ndarray::Array< T1, 2, C1 > const &fisher, ndarray::Array< T2, 1, C2 > const &rhs)
	Reset the terms in the normal equations given as ndarrays; dimension must not change. More...
template<typename D1 , typename D2 >
void	setNormalEquations (Eigen::MatrixBase< D1 > const &fisher, Eigen::MatrixBase< D2 > const &rhs)
	Reset the terms in the normal equations given as Eigen objects; dimension must not change. More...
void	setThreshold (double threshold)
	Set the threshold used to determine when to truncate Eigenvalues. More...
double	getThreshold () const
	Get the threshold used to determine when to truncate Eigenvalues. More...
ndarray::Array< double const, 1, 1 >	getSolution ()
	Return the vector solution to the least squares problem. More...
ndarray::Array< double const, 2, 2 >	getCovariance ()
	Return the covariance matrix of the least squares problem. More...
ndarray::Array< double const, 2, 2 >	getFisherMatrix ()
	Return the Fisher matrix (inverse of the covariance) of the parameters. More...
ndarray::Array< double const, 1, 1 >	getDiagnostic (Factorization factorization)
	Return a factorization-dependent vector that can be used to characterize the stability of the solution. More...
int	getDimension () const
	Return the number of parameters. More...
int	getRank () const
	Return the rank of the problem (number of nonzero Eigenvalues). More...
Factorization	getFactorization () const
	Retun the type of factorization used by the solver. More...
	LeastSquares (Factorization factorization, int dimension)
	Construct a least-squares object for the given factorization and dimensionality. More...
	~LeastSquares ()

Static Public Member Functions
template<typename T1 , typename T2 , int C1, int C2>
static LeastSquares	fromDesignMatrix (ndarray::Array< T1, 2, C1 > const &design, ndarray::Array< T2, 1, C2 > const &data, Factorization factorization=NORMAL_EIGENSYSTEM)
	Initialize from the design matrix and data vector given as ndarrays. More...
template<typename D1 , typename D2 >
static LeastSquares	fromDesignMatrix (Eigen::MatrixBase< D1 > const &design, Eigen::MatrixBase< D2 > const &data, Factorization factorization=NORMAL_EIGENSYSTEM)
	Initialize from the design matrix and data vector given as an Eigen objects. More...
template<typename T1 , typename T2 , int C1, int C2>
static LeastSquares	fromNormalEquations (ndarray::Array< T1, 2, C1 > const &fisher, ndarray::Array< T2, 1, C2 > const &rhs, Factorization factorization=NORMAL_EIGENSYSTEM)
	Initialize from the terms in the normal equations, given as ndarrays. More...
template<typename D1 , typename D2 >
static LeastSquares	fromNormalEquations (Eigen::MatrixBase< D1 > const &fisher, Eigen::MatrixBase< D2 > const &rhs, Factorization factorization=NORMAL_EIGENSYSTEM)
	Initialize from the terms in the normal equations, given as Eigen objects. More...

Private Member Functions
Eigen::MatrixXd &	_getDesignMatrix ()
Eigen::VectorXd &	_getDataVector ()
Eigen::MatrixXd &	_getFisherMatrix ()
Eigen::VectorXd &	_getRhsVector ()
void	_factor (bool haveNormalEquations)

Private Attributes
boost::shared_ptr< Impl >	_impl

Detailed Description

Solver for linear least-squares problems.

Linear least-squares problems are defined as finding the vector \(x\) that minimizes \(\left|A x - b\right|_2\), with the number of rows of \(A\) generally greater than the number of columns. We call \(A\) the design matrix, \(b\) the data vector, and \(x\) the solution vector. When the rank of \(A\) is equal to the number of columns, we can obtain using the solution using the normal equations:

\[ A^T A x = A^T b \]

If \(A\) is not full-rank, the problem is underconstrained, and we usually wish to solve the minimum-norm least-squares problem, which also minimizes \(|x|_2\). This can be done by computing a pseudo-inverse of \(A\) using an SVD or complete orthogonal factorization of \(A\), or by performing an Eigen decomposition of \(A^T A\).

This class can be constructed from the design matrix and data vector, or from the two terms in the normal equations (below, we call the matrix \(A^TA\) the Fisher matrix, and the vector \(A^T b\) simply the "right-hand side" (RHS) vector). If initialized with the design matrix and data vector, we can still use the normal equations to solve it. The solution via the normal equations is more susceptible to round-off error, but it is also faster, and if the normal equation terms can be computed directly it can be significantly less expensive in terms of memory. The Fisher matrix is a symmetric matrix, and it should be exactly symmetric when provided as input, because which triangle will be used is an implementation detail that is subject to change and may depend on the storage order.

The solver always operates in double precision, and returns all results in double precision. However, it can be initialized from single precision inputs. It isn't usually a good idea to construct the normal equations in single precision, however, even when the data are single precision.

Definition at line 65 of file LeastSquares.h.

Member Enumeration Documentation

enum lsst::afw::math::LeastSquares::Factorization

Private implementation; forward-declared publicly so we can inherit from it in .cc.

Enumerator
NORMAL_EIGENSYSTEM	Use the normal equations with a symmetric Eigensystem decomposition. This method is fully robust and computes the minimum-norm solution when the problem does not have full rank. It is affected slightly more by round-off error than the SVD method, but because it can handle singular matrices this usually isn't a problem.
NORMAL_CHOLESKY	Use the normal equations with a Cholesky decomposition. While this method uses a robust LDL^T decomposition that does not require square roots, it is not appropriate for problems that do not have full rank, and cannot be used to determine whether a problem has full rank. It is the fastest decomposition.
DIRECT_SVD	Use a thin singular value decomposition of the design matrix. This method is the most robust and computes the minimum-norm solution the problem does not have full rank. However, it is also the slowest and is not available when the solver is initialized with the normal equations.

Definition at line 70 of file LeastSquares.h.

                       {
        NORMAL_EIGENSYSTEM,  
        NORMAL_CHOLESKY,     
        DIRECT_SVD           
    };

Constructor & Destructor Documentation

lsst::afw::math::LeastSquares::LeastSquares	(	Factorization	factorization,
		int	dimension
	)

Construct a least-squares object for the given factorization and dimensionality.

One of the set* member functions must be called before any other operations can be performed on a LeastSquares object initialized this way.

Definition at line 363 of file LeastSquares.cc.

                                                                     {
    switch (factorization) {
    case NORMAL_EIGENSYSTEM:
        _impl = boost::make_shared<EigensystemSolver>(dimension);
        break;
    case NORMAL_CHOLESKY:
        _impl = boost::make_shared<CholeskySolver>(dimension);
        break;
    case DIRECT_SVD:
        _impl = boost::make_shared<SvdSolver>(dimension);
        break;
    }
    _impl->factorization = factorization;
}

lsst::afw::math::LeastSquares::~LeastSquares

(

)

Definition at line 378 of file LeastSquares.cc.

{}

Member Function Documentation

void lsst::afw::math::LeastSquares::_factor ( bool  haveNormalEquations )

private

Definition at line 386 of file LeastSquares.cc.

                                                   {
    if (haveNormalEquations) {
        if (_getFisherMatrix().rows() != _impl->dimension) {
            throw LSST_EXCEPT(
                pex::exceptions::InvalidParameterError,
                (boost::format("Number of rows of Fisher matrix (%d) does not match"
                               " dimension of LeastSquares solver.")
                 % _getFisherMatrix().rows() % _impl->dimension).str()
            );
        }
        if (_getFisherMatrix().cols() != _impl->dimension) {
            throw LSST_EXCEPT(
                pex::exceptions::InvalidParameterError,
                (boost::format("Number of columns of Fisher matrix (%d) does not match"
                               " dimension of LeastSquares solver.")
                 % _getFisherMatrix().cols() % _impl->dimension).str()
            );
        }
        if (_getRhsVector().size() != _impl->dimension) {
            throw LSST_EXCEPT(
                pex::exceptions::InvalidParameterError,
                (boost::format("Number of elements in RHS vector (%d) does not match"
                               " dimension of LeastSquares solver.")
                 % _getRhsVector().size() % _impl->dimension).str()
            );
        }
        _impl->state = Impl::RHS_VECTOR | Impl::FULL_FISHER_MATRIX | Impl::LOWER_FISHER_MATRIX;
    } else {
        if (_getDesignMatrix().cols() != _impl->dimension) {
            throw LSST_EXCEPT(
                pex::exceptions::InvalidParameterError,
                "Number of columns of design matrix does not match dimension of LeastSquares solver."
            );
        }
        if (_getDesignMatrix().rows() != _getDataVector().size()) {
            throw LSST_EXCEPT(
                pex::exceptions::InvalidParameterError,
                (boost::format("Number of rows of design matrix (%d) does not match number of "
                               "data points (%d)") % _getDesignMatrix().rows() % _getDataVector().size()
                ).str()
            );
        }
        if (_getDesignMatrix().cols() > _getDataVector().size()) {
            throw LSST_EXCEPT(
                pex::exceptions::InvalidParameterError,
                (boost::format("Number of columns of design matrix (%d) must be smaller than number of "
                               "data points (%d)") % _getDesignMatrix().cols() % _getDataVector().size()
                ).str()
            );
        }
        _impl->state = Impl::DESIGN_AND_DATA;
    }
    _impl->factor();
}

Eigen::VectorXd & lsst::afw::math::LeastSquares::_getDataVector ( )

private

Definition at line 381 of file LeastSquares.cc.

{ return _impl->data; }

Eigen::MatrixXd & lsst::afw::math::LeastSquares::_getDesignMatrix ( )

private

Definition at line 380 of file LeastSquares.cc.

{ return _impl->design; }

Eigen::MatrixXd & lsst::afw::math::LeastSquares::_getFisherMatrix ( )

private

Definition at line 383 of file LeastSquares.cc.

{ return _impl->fisher; }

Eigen::VectorXd & lsst::afw::math::LeastSquares::_getRhsVector ( )

private

Definition at line 384 of file LeastSquares.cc.

{ return _impl->rhs; }

template<typename T1 , typename T2 , int C1, int C2>

static LeastSquares lsst::afw::math::LeastSquares::fromDesignMatrix ( ndarray::Array< T1, 2, C1 > const &  design, ndarray::Array< T2, 1, C2 > const &  data, Factorization  factorization = NORMAL_EIGENSYSTEM  )

inlinestatic

Initialize from the design matrix and data vector given as ndarrays.

Definition at line 99 of file LeastSquares.h.

      {
        LeastSquares r(factorization, design.template getSize<1>());
        r.setDesignMatrix(design, data);
        return r;
    }

template<typename D1 , typename D2 >

static LeastSquares lsst::afw::math::LeastSquares::fromDesignMatrix ( Eigen::MatrixBase< D1 > const &  design, Eigen::MatrixBase< D2 > const &  data, Factorization  factorization = NORMAL_EIGENSYSTEM  )

inlinestatic

Initialize from the design matrix and data vector given as an Eigen objects.

Definition at line 111 of file LeastSquares.h.

      {
        LeastSquares r(factorization, design.cols());
        r.setDesignMatrix(design, data);
        return r;
    }

template<typename T1 , typename T2 , int C1, int C2>

static LeastSquares lsst::afw::math::LeastSquares::fromNormalEquations ( ndarray::Array< T1, 2, C1 > const &  fisher, ndarray::Array< T2, 1, C2 > const &  rhs, Factorization  factorization = NORMAL_EIGENSYSTEM  )

inlinestatic

Initialize from the terms in the normal equations, given as ndarrays.

Definition at line 163 of file LeastSquares.h.

      {
        LeastSquares r(factorization, fisher.template getSize<0>());
        r.setNormalEquations(fisher, rhs);
        return r;
    }

template<typename D1 , typename D2 >

static LeastSquares lsst::afw::math::LeastSquares::fromNormalEquations ( Eigen::MatrixBase< D1 > const &  fisher, Eigen::MatrixBase< D2 > const &  rhs, Factorization  factorization = NORMAL_EIGENSYSTEM  )

inlinestatic

Initialize from the terms in the normal equations, given as Eigen objects.

Definition at line 175 of file LeastSquares.h.

      {
        LeastSquares r(factorization, fisher.rows());
        r.setNormalEquations(fisher, rhs);
        return r;
    }

ndarray::Array< double const, 2, 2 > lsst::afw::math::LeastSquares::getCovariance

(

)

Return the covariance matrix of the least squares problem.

The returned array is owned by the LeastSquares object and may be modified in-place by future calls to LeastSquares member functions, so it's best to promptly copy the result elsewhere.

If you want an Eigen object instead, just use getCovariance().asEigen().

The covariance is cached the first time this member function is called, and will be reused unless the matrices are reset or the threshold is changed.

Definition at line 331 of file LeastSquares.cc.

                                                           {
    _impl->ensure(Impl::COVARIANCE_ARRAY);
    return _impl->covariance;
}

ndarray::Array< double const, 1, 1 > lsst::afw::math::LeastSquares::getDiagnostic

(

Factorization

factorization

)

Return a factorization-dependent vector that can be used to characterize the stability of the solution.

The returned array's size is always equal to getDimension(). When getRank() is less than getDimension(), some elements of the array were considered effectively zero (see setThreshold).

For the NORMAL_EIGENSYSTEM method, this is the vector of Eigenvalues of the Fisher matrix, including those rejected as being below the threshold, sorted in descending order (all values are nonnegative). This is the square of the singular values of the design matrix, and is only available when the factorization is either NORMAL_EIGENSYSTEM or DIRECT_SVD.

For the DIRECT_SVD method, this is the vector of singular values of the design matrix, including those rejected as being below the threshold, sorted in descending order (all values are nonnegative). This is the square root of the Eigenvalues of the Fisher matrix, and is only available when the factorization is either NORMAL_EIGENSYSTEM or DIRECT_SVD.

For the NORMAL_CHOLESKY method, this is \(D\) in the pivoted Cholesky factorization \(P L D L^T P^T\) of the Fisher matrix. This does not provide a reliable way to test the stability of the problem, but it does provide a way to compute the determinant of the Fisher matrix. It is only available when the factorization is NORMAL_CHOLESKY.

Definition at line 348 of file LeastSquares.cc.

                                                                                      {
    if (_impl->whichDiagnostic != factorization) {
        _impl->state &= ~Impl::DIAGNOSTIC_ARRAY;
        _impl->whichDiagnostic = factorization;
    }
    _impl->ensure(Impl::DIAGNOSTIC_ARRAY);
    return _impl->diagnostic;
}

int lsst::afw::math::LeastSquares::getDimension

(

)

const

Return the number of parameters.

Definition at line 357 of file LeastSquares.cc.

{ return _impl->dimension; }

LeastSquares::Factorization lsst::afw::math::LeastSquares::getFactorization

(

)

const

Retun the type of factorization used by the solver.

Definition at line 361 of file LeastSquares.cc.

{ return _impl->factorization; }

ndarray::Array< double const, 2, 2 > lsst::afw::math::LeastSquares::getFisherMatrix

(

)

Return the Fisher matrix (inverse of the covariance) of the parameters.

Note that the Fisher matrix is exactly the same as the matrix on the lhs of the normal equations.

The returned array is owned by the LeastSquares object and may be modified in-place by future calls to LeastSquares member functions, so it's best to promptly copy the result elsewhere.

If you want an Eigen object instead, just use getFisherMatrix().asEigen().

Definition at line 336 of file LeastSquares.cc.

                                                             {
    _impl->ensure(Impl::FULL_FISHER_MATRIX);
    // Wrap the Eigen::MatrixXd in an ndarray::Array, using _impl as the reference-counted owner.
    // Doesn't matter if we swap strides, because it's symmetric.
    return ndarray::external(
        _impl->fisher.data(),
        ndarray::makeVector(_impl->dimension, _impl->dimension),
        ndarray::makeVector(_impl->dimension, 1),
        _impl
    );
}

int lsst::afw::math::LeastSquares::getRank

(

)

const

Return the rank of the problem (number of nonzero Eigenvalues).

The returned value is always the same as getDimension() when the factorization is NORMAL_CHOLESKY (which may be incorrect, because a Cholesky decomposition is not rank-revealing).

Definition at line 359 of file LeastSquares.cc.

{ return _impl->rank; }

ndarray::Array< double const, 1, 1 > lsst::afw::math::LeastSquares::getSolution

(

)

Return the vector solution to the least squares problem.

The returned array is owned by the LeastSquares object and may be modified in-place by future calls to LeastSquares member functions, so it's best to promptly copy the result elsewhere.

If you want an Eigen object instead, just use getSolution().asEigen().

The solution is cached the first time this member function is called, and will be reused unless the matrices are reset or the threshold is changed.

Definition at line 326 of file LeastSquares.cc.

                                                         {
    _impl->ensure(Impl::SOLUTION_ARRAY);
    return _impl->solution;
}

double lsst::afw::math::LeastSquares::getThreshold

(

)

const

Get the threshold used to determine when to truncate Eigenvalues.

Definition at line 324 of file LeastSquares.cc.

{ return _impl->threshold; }

template<typename T1 , typename T2 , int C1, int C2>

void lsst::afw::math::LeastSquares::setDesignMatrix ( ndarray::Array< T1, 2, C1 > const &  design, ndarray::Array< T2, 1, C2 > const &  data  )

inline

Reset the design matrix and data vector given as ndarrays; dimension must not change.

Definition at line 123 of file LeastSquares.h.

      {
        // n.b. "template cast<T>" below is not a special kind of cast; it's just
        // the weird C++ syntax required for calling a templated member function
        // called "cast" in this context; see
        // http://eigen.tuxfamily.org/dox-devel/TopicTemplateKeyword.html
        _getDesignMatrix() = design.asEigen().template cast<double>();
        _getDataVector() = data.asEigen().template cast<double>();
        _factor(false);
    }

template<typename D1 , typename D2 >

void lsst::afw::math::LeastSquares::setDesignMatrix ( Eigen::MatrixBase< D1 > const &  design, Eigen::MatrixBase< D2 > const &  data  )

inline

Reset the design matrix and data vector given as Eigen objects; dimension must not change.

Definition at line 138 of file LeastSquares.h.

      {
        _getDesignMatrix() = design.template cast<double>();
        _getDataVector() = data.template cast<double>();
        _factor(false);
    }

template<typename T1 , int C1>

void lsst::afw::math::LeastSquares::setDesignMatrix ( ndarray::Array< T1, 2, C1 > const &  design )

inline

Reset the design matrix given as an ndarray; dimension and data are not changed.

Definition at line 149 of file LeastSquares.h.

                                                               {
        _getDesignMatrix() = design.asEigen().template cast<double>();
        _factor(false);
    }

template<typename D1 , typename D2 >

void lsst::afw::math::LeastSquares::setDesignMatrix ( Eigen::MatrixBase< D1 > const &  design )

inline

Reset the design matrix given as an Eigen object; dimension and data are not changed.

Definition at line 156 of file LeastSquares.h.

                                                             {
        _getDesignMatrix() = design.template cast<double>();
        _factor(false);
    }

template<typename T1 , typename T2 , int C1, int C2>

void lsst::afw::math::LeastSquares::setNormalEquations ( ndarray::Array< T1, 2, C1 > const &  fisher, ndarray::Array< T2, 1, C2 > const &  rhs  )

inline

Reset the terms in the normal equations given as ndarrays; dimension must not change.

Definition at line 187 of file LeastSquares.h.

      {
        if ((C1 > 0) == bool(Eigen::MatrixXd::IsRowMajor)) {
            _getFisherMatrix() = fisher.asEigen().template cast<double>();
        } else {
            _getFisherMatrix() = fisher.asEigen().transpose().template cast<double>();
        }
        _getRhsVector() = rhs.asEigen().template cast<double>();
        _factor(true);
    }

template<typename D1 , typename D2 >

void lsst::afw::math::LeastSquares::setNormalEquations ( Eigen::MatrixBase< D1 > const &  fisher, Eigen::MatrixBase< D2 > const &  rhs  )

inline

Reset the terms in the normal equations given as Eigen objects; dimension must not change.

Definition at line 202 of file LeastSquares.h.

      {
        if (bool(Eigen::MatrixBase<D1>::IsRowMajor) == bool(Eigen::MatrixXd::IsRowMajor)) {
            _getFisherMatrix() = fisher.template cast<double>();
        } else {
            _getFisherMatrix() = fisher.transpose().template cast<double>();
        }
        _getRhsVector() = rhs.template cast<double>();
        _factor(true);
    }

void lsst::afw::math::LeastSquares::setThreshold

(

double

threshold

)

Set the threshold used to determine when to truncate Eigenvalues.

The rank of the matrix is determined by comparing the product of this threshold and the first (largest) element of the array returned by getDiagnostic() to all other elements of that array. Elements smaller than this product are ignored and reduce the rank.

In the DIRECT_SVD case, the diagnostic array contains the singular values of the design matrix, while in the NORMAL_EIGENSYSTEM case the diagnostic array holds the Eigenvalues of the Fisher matrix, and the latter are the square of the former. The default threshold is the same (std::numeric_limits<double>::epsilon()) in both cases, reflecting the fact that using the normal equations squares the condition number of the problem.

The NORMAL_CHOLESKY method does not use the threshold and assumes the problem is full-rank.

Definition at line 317 of file LeastSquares.cc.

                                                {
    _impl->threshold = threshold;
    _impl->state &= ~Impl::SOLUTION_ARRAY;
    _impl->state &= ~Impl::COVARIANCE_ARRAY;
    _impl->updateRank();
}

Member Data Documentation

boost::shared_ptr< Impl > lsst::afw::math::LeastSquares::_impl

private

Definition at line 353 of file LeastSquares.h.

---

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

• /home/lsstsw/stack/Linux64/afw/11.0-2-g04d2804/include/lsst/afw/math/LeastSquares.h
• /home/lsstsw/stack/Linux64/afw/11.0-2-g04d2804/src/math/LeastSquares.cc