LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
LeastSquares.h
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008-2013 LSST Corporation.
6  *
7  * This product includes software developed by the
8  * LSST Project (http://www.lsst.org/).
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the LSST License Statement and
21  * the GNU General Public License along with this program. If not,
22  * see <http://www.lsstcorp.org/LegalNotices/>.
23  */
24 
25 #ifndef LSST_AFW_MATH_LeastSquares_h_INCLUDED
26 #define LSST_AFW_MATH_LeastSquares_h_INCLUDED
27 
28 #include "lsst/base.h"
29 #include "ndarray/eigen.h"
30 
31 namespace lsst { namespace afw { namespace math {
32 
65 class LeastSquares {
66 public:
67 
68  class Impl;
69 
95  };
96 
98  template <typename T1, typename T2, int C1, int C2>
100  ndarray::Array<T1,2,C1> const & design,
101  ndarray::Array<T2,1,C2> const & data,
102  Factorization factorization = NORMAL_EIGENSYSTEM
103  ) {
104  LeastSquares r(factorization, design.template getSize<1>());
105  r.setDesignMatrix(design, data);
106  return r;
107  }
108 
110  template <typename D1, typename D2>
112  Eigen::MatrixBase<D1> const & design,
113  Eigen::MatrixBase<D2> const & data,
114  Factorization factorization = NORMAL_EIGENSYSTEM
115  ) {
116  LeastSquares r(factorization, design.cols());
117  r.setDesignMatrix(design, data);
118  return r;
119  }
120 
122  template <typename T1, typename T2, int C1, int C2>
124  ndarray::Array<T1,2,C1> const & design,
125  ndarray::Array<T2,1,C2> const & data
126  ) {
127  // n.b. "template cast<T>" below is not a special kind of cast; it's just
128  // the weird C++ syntax required for calling a templated member function
129  // called "cast" in this context; see
130  // http://eigen.tuxfamily.org/dox-devel/TopicTemplateKeyword.html
131  _getDesignMatrix() = design.asEigen().template cast<double>();
132  _getDataVector() = data.asEigen().template cast<double>();
133  _factor(false);
134  }
135 
137  template <typename D1, typename D2>
139  Eigen::MatrixBase<D1> const & design,
140  Eigen::MatrixBase<D2> const & data
141  ) {
142  _getDesignMatrix() = design.template cast<double>();
143  _getDataVector() = data.template cast<double>();
144  _factor(false);
145  }
146 
148  template <typename T1, int C1>
150  _getDesignMatrix() = design.asEigen().template cast<double>();
151  _factor(false);
152  }
153 
155  template <typename D1, typename D2>
156  void setDesignMatrix(Eigen::MatrixBase<D1> const & design) {
157  _getDesignMatrix() = design.template cast<double>();
158  _factor(false);
159  }
160 
162  template <typename T1, typename T2, int C1, int C2>
164  ndarray::Array<T1,2,C1> const & fisher,
165  ndarray::Array<T2,1,C2> const & rhs,
166  Factorization factorization = NORMAL_EIGENSYSTEM
167  ) {
168  LeastSquares r(factorization, fisher.template getSize<0>());
169  r.setNormalEquations(fisher, rhs);
170  return r;
171  }
172 
174  template <typename D1, typename D2>
176  Eigen::MatrixBase<D1> const & fisher,
177  Eigen::MatrixBase<D2> const & rhs,
178  Factorization factorization = NORMAL_EIGENSYSTEM
179  ) {
180  LeastSquares r(factorization, fisher.rows());
181  r.setNormalEquations(fisher, rhs);
182  return r;
183  }
184 
186  template <typename T1, typename T2, int C1, int C2>
188  ndarray::Array<T1,2,C1> const & fisher,
189  ndarray::Array<T2,1,C2> const & rhs
190  ) {
191  if ((C1 > 0) == bool(Eigen::MatrixXd::IsRowMajor)) {
192  _getFisherMatrix() = fisher.asEigen().template cast<double>();
193  } else {
194  _getFisherMatrix() = fisher.asEigen().transpose().template cast<double>();
195  }
196  _getRhsVector() = rhs.asEigen().template cast<double>();
197  _factor(true);
198  }
199 
201  template <typename D1, typename D2>
203  Eigen::MatrixBase<D1> const & fisher,
204  Eigen::MatrixBase<D2> const & rhs
205  ) {
206  if (bool(Eigen::MatrixBase<D1>::IsRowMajor) == bool(Eigen::MatrixXd::IsRowMajor)) {
207  _getFisherMatrix() = fisher.template cast<double>();
208  } else {
209  _getFisherMatrix() = fisher.transpose().template cast<double>();
210  }
211  _getRhsVector() = rhs.template cast<double>();
212  _factor(true);
213  }
214 
233  void setThreshold(double threshold);
234 
236  double getThreshold() const;
237 
251 
265 
279 
306 
308  int getDimension() const;
309 
316  int getRank() const;
317 
320 
327  LeastSquares(Factorization factorization, int dimension);
328 
329  // Need to define dtor in source file so it can see Impl declaration.
330  ~LeastSquares();
331 
332 private:
333 
334  // We want a column-major design matrix so the self-adjoint product is cache-friendly, so we always
335  // copy a (possibly row-major) design matrix into a col-major one. This is an unnecessarily and
336  // cache-unfriendly operation when solver is DIRECT_SVD, but right now it doesn't seem to be worth
337  // special-casing the design for that case. In other cases it's a cache-unfriendly op that avoids
338  // an even worse one, and it can always be avoided by using a column-major design matrix.
339  Eigen::MatrixXd & _getDesignMatrix();
340  Eigen::VectorXd & _getDataVector();
341 
342  // Storage order matters less here, so we just go with what Eigen is most accustomed to.
343  // Because the Fisher matrix is symmetric, we can just transpose it before copying to avoid doing
344  // any expensive copies between different storage orders.
345  Eigen::MatrixXd & _getFisherMatrix();
346  Eigen::VectorXd & _getRhsVector();
347 
348  // Do the actual high-level linear algrebra factorization; the appropriate input matrices and
349  // vectors must already be set before this is called. The solution and other quantities
350  // will not be computed until they are first requested.
351  void _factor(bool haveNormalEquations);
352 
354 };
355 
356 }}} // namespace lsst::afw::math
357 
358 #endif // !LSST_AFW_MATH_LeastSquares_h_INCLUDED
int getDimension() const
Return the number of parameters.
ndarray::Array< double const, 2, 2 > getCovariance()
Return the covariance matrix of the least squares problem.
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...
Definition: LeastSquares.h:138
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.
Definition: LeastSquares.h:163
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.
Definition: LeastSquares.h:99
#define PTR(...)
Definition: base.h:41
Eigen matrix objects that present a view into an ndarray::Array.
Eigen::VectorXd & _getDataVector()
Eigen::VectorXd & _getRhsVector()
void setDesignMatrix(Eigen::MatrixBase< D1 > const &design)
Reset the design matrix given as an Eigen object; dimension and data are not changed.
Definition: LeastSquares.h:156
EigenView< Element, ND::value, RMC::value, XprKind, Rows, Cols > asEigen() const
Solver for linear least-squares problems.
Definition: LeastSquares.h:65
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.
Definition: LeastSquares.h:123
int getRank() const
Return the rank of the problem (number of nonzero Eigenvalues).
Use a thin singular value decomposition of the design matrix.
Definition: LeastSquares.h:87
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.
Definition: LeastSquares.h:187
Use the normal equations with a Cholesky decomposition.
Definition: LeastSquares.h:79
Eigen::MatrixXd & _getDesignMatrix()
void setThreshold(double threshold)
Set the threshold used to determine when to truncate Eigenvalues.
Factorization getFactorization() const
Retun the type of factorization used by the solver.
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.
Definition: LeastSquares.h:111
double getThreshold() const
Get the threshold used to determine when to truncate Eigenvalues.
A multidimensional strided array.
Definition: Array.h:47
Factorization
Private implementation; forward-declared publicly so we can inherit from it in .cc.
Definition: LeastSquares.h:70
void _factor(bool haveNormalEquations)
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.
Definition: LeastSquares.h:175
Use the normal equations with a symmetric Eigensystem decomposition.
Definition: LeastSquares.h:71
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...
Definition: LeastSquares.h:202
ndarray::Array< double const, 2, 2 > getFisherMatrix()
Return the Fisher matrix (inverse of the covariance) of the parameters.
boost::shared_ptr< Impl > _impl
Definition: LeastSquares.h:353
Eigen::MatrixXd & _getFisherMatrix()
ndarray::Array< double const, 1, 1 > getSolution()
Return the vector solution to the least squares problem.
ndarray::Array< double const, 1, 1 > getDiagnostic(Factorization factorization)
Return a factorization-dependent vector that can be used to characterize the stability of the solutio...
void setDesignMatrix(ndarray::Array< T1, 2, C1 > const &design)
Reset the design matrix given as an ndarray; dimension and data are not changed.
Definition: LeastSquares.h:149
LeastSquares(Factorization factorization, int dimension)
Construct a least-squares object for the given factorization and dimensionality.