LSSTApplications  17.0+11,17.0+34,17.0+56,17.0+57,17.0+59,17.0+7,17.0-1-g377950a+33,17.0.1-1-g114240f+2,17.0.1-1-g4d4fbc4+28,17.0.1-1-g55520dc+49,17.0.1-1-g5f4ed7e+52,17.0.1-1-g6dd7d69+17,17.0.1-1-g8de6c91+11,17.0.1-1-gb9095d2+7,17.0.1-1-ge9fec5e+5,17.0.1-1-gf4e0155+55,17.0.1-1-gfc65f5f+50,17.0.1-1-gfc6fb1f+20,17.0.1-10-g87f9f3f+1,17.0.1-11-ge9de802+16,17.0.1-16-ga14f7d5c+4,17.0.1-17-gc79d625+1,17.0.1-17-gdae4c4a+8,17.0.1-2-g26618f5+29,17.0.1-2-g54f2ebc+9,17.0.1-2-gf403422+1,17.0.1-20-g2ca2f74+6,17.0.1-23-gf3eadeb7+1,17.0.1-3-g7e86b59+39,17.0.1-3-gb5ca14a,17.0.1-3-gd08d533+40,17.0.1-30-g596af8797,17.0.1-4-g59d126d+4,17.0.1-4-gc69c472+5,17.0.1-6-g5afd9b9+4,17.0.1-7-g35889ee+1,17.0.1-7-gc7c8782+18,17.0.1-9-gc4bbfb2+3,w.2019.22
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 {
32 namespace afw {
33 namespace math {
34 
67 class LeastSquares final {
68 public:
69  class Impl;
70 
96  };
97 
99  template <typename T1, typename T2, int C1, int C2>
100  static LeastSquares fromDesignMatrix(ndarray::Array<T1, 2, C1> const& design,
101  ndarray::Array<T2, 1, C2> const& data,
102  Factorization factorization = NORMAL_EIGENSYSTEM) {
103  LeastSquares r(factorization, design.template getSize<1>());
104  r.setDesignMatrix(design, data);
105  return r;
106  }
107 
109  template <typename D1, typename D2>
110  static LeastSquares fromDesignMatrix(Eigen::MatrixBase<D1> const& design,
111  Eigen::MatrixBase<D2> const& data,
112  Factorization factorization = NORMAL_EIGENSYSTEM) {
113  LeastSquares r(factorization, design.cols());
114  r.setDesignMatrix(design, data);
115  return r;
116  }
117 
119  template <typename T1, typename T2, int C1, int C2>
120  void setDesignMatrix(ndarray::Array<T1, 2, C1> const& design, ndarray::Array<T2, 1, C2> const& data) {
121  // n.b. "template cast<T>" below is not a special kind of cast; it's just
122  // the weird C++ syntax required for calling a templated member function
123  // called "cast" in this context; see
124  // http://eigen.tuxfamily.org/dox-devel/TopicTemplateKeyword.html
125  _getDesignMatrix() = ndarray::asEigenMatrix(design).template cast<double>();
126  _getDataVector() = ndarray::asEigenMatrix(data).template cast<double>();
127  _factor(false);
128  }
129 
131  template <typename D1, typename D2>
132  void setDesignMatrix(Eigen::MatrixBase<D1> const& design, Eigen::MatrixBase<D2> const& data) {
133  _getDesignMatrix() = design.template cast<double>();
134  _getDataVector() = data.template cast<double>();
135  _factor(false);
136  }
137 
139  template <typename T1, int C1>
140  void setDesignMatrix(ndarray::Array<T1, 2, C1> const& design) {
141  _getDesignMatrix() = ndarray::asEigenMatrix(design).template cast<double>();
142  _factor(false);
143  }
144 
146  template <typename D1, typename D2>
147  void setDesignMatrix(Eigen::MatrixBase<D1> const& design) {
148  _getDesignMatrix() = design.template cast<double>();
149  _factor(false);
150  }
151 
153  template <typename T1, typename T2, int C1, int C2>
154  static LeastSquares fromNormalEquations(ndarray::Array<T1, 2, C1> const& fisher,
155  ndarray::Array<T2, 1, C2> const& rhs,
156  Factorization factorization = NORMAL_EIGENSYSTEM) {
157  LeastSquares r(factorization, fisher.template getSize<0>());
158  r.setNormalEquations(fisher, rhs);
159  return r;
160  }
161 
163  template <typename D1, typename D2>
164  static LeastSquares fromNormalEquations(Eigen::MatrixBase<D1> const& fisher,
165  Eigen::MatrixBase<D2> const& rhs,
166  Factorization factorization = NORMAL_EIGENSYSTEM) {
167  LeastSquares r(factorization, fisher.rows());
168  r.setNormalEquations(fisher, rhs);
169  return r;
170  }
171 
173  template <typename T1, typename T2, int C1, int C2>
174  void setNormalEquations(ndarray::Array<T1, 2, C1> const& fisher, ndarray::Array<T2, 1, C2> const& rhs) {
175  if ((C1 > 0) == bool(Eigen::MatrixXd::IsRowMajor)) {
176  _getFisherMatrix() = ndarray::asEigenMatrix(fisher).template cast<double>();
177  } else {
178  _getFisherMatrix() = ndarray::asEigenMatrix(fisher).transpose().template cast<double>();
179  }
180  _getRhsVector() = ndarray::asEigenMatrix(rhs).template cast<double>();
181  _factor(true);
182  }
183 
185  template <typename D1, typename D2>
186  void setNormalEquations(Eigen::MatrixBase<D1> const& fisher, Eigen::MatrixBase<D2> const& rhs) {
187  if (bool(Eigen::MatrixBase<D1>::IsRowMajor) == bool(Eigen::MatrixXd::IsRowMajor)) {
188  _getFisherMatrix() = fisher.template cast<double>();
189  } else {
190  _getFisherMatrix() = fisher.transpose().template cast<double>();
191  }
192  _getRhsVector() = rhs.template cast<double>();
193  _factor(true);
194  }
195 
214  void setThreshold(double threshold);
215 
217  double getThreshold() const;
218 
231  ndarray::Array<double const, 1, 1> getSolution();
232 
245  ndarray::Array<double const, 2, 2> getCovariance();
246 
259  ndarray::Array<double const, 2, 2> getFisherMatrix();
260 
286  ndarray::Array<double const, 1, 1> getDiagnostic(Factorization factorization);
287 
289  int getDimension() const;
290 
297  int getRank() const;
298 
301 
308  LeastSquares(Factorization factorization, int dimension);
309 
310  LeastSquares(LeastSquares const&);
314 
315  // Need to define dtor in source file so it can see Impl declaration.
316  ~LeastSquares();
317 
318 private:
319  // We want a column-major design matrix so the self-adjoint product is cache-friendly, so we always
320  // copy a (possibly row-major) design matrix into a col-major one. This is an unnecessarily and
321  // cache-unfriendly operation when solver is DIRECT_SVD, but right now it doesn't seem to be worth
322  // special-casing the design for that case. In other cases it's a cache-unfriendly op that avoids
323  // an even worse one, and it can always be avoided by using a column-major design matrix.
324  Eigen::MatrixXd& _getDesignMatrix();
325  Eigen::VectorXd& _getDataVector();
326 
327  // Storage order matters less here, so we just go with what Eigen is most accustomed to.
328  // Because the Fisher matrix is symmetric, we can just transpose it before copying to avoid doing
329  // any expensive copies between different storage orders.
330  Eigen::MatrixXd& _getFisherMatrix();
331  Eigen::VectorXd& _getRhsVector();
332 
333  // Do the actual high-level linear algrebra factorization; the appropriate input matrices and
334  // vectors must already be set before this is called. The solution and other quantities
335  // will not be computed until they are first requested.
336  void _factor(bool haveNormalEquations);
337 
338  std::shared_ptr<Impl> _impl;
339 };
340 } // namespace math
341 } // namespace afw
342 } // namespace lsst
343 
344 #endif // !LSST_AFW_MATH_LeastSquares_h_INCLUDED
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:154
void setThreshold(double threshold)
Set the threshold used to determine when to truncate Eigenvalues.
Basic LSST definitions.
int getRank() const
Return the rank of the problem (number of nonzero Eigenvalues).
Solver for linear least-squares problems.
Definition: LeastSquares.h:67
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:140
LeastSquares(Factorization factorization, int dimension)
Construct a least-squares object for the given factorization and dimensionality.
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:100
Use a thin singular value decomposition of the design matrix.
Definition: LeastSquares.h:88
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:174
Use the normal equations with a Cholesky decomposition.
Definition: LeastSquares.h:80
Use the normal equations with a symmetric Eigensystem decomposition.
Definition: LeastSquares.h:72
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(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:132
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:120
A base class for image defects.
char * data
Definition: BaseRecord.cc:62
int getDimension() const
Return the number of parameters.
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:147
double getThreshold() const
Get the threshold used to determine when to truncate Eigenvalues.
Factorization getFactorization() const
Retun the type of factorization used by the solver.
ndarray::Array< double const, 2, 2 > getCovariance()
Return the covariance matrix of the least squares problem.
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:110
Factorization
Private implementation; forward-declared publicly so we can inherit from it in .cc.
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:186
ndarray::Array< double const, 1, 1 > getSolution()
Return the vector solution to the least squares problem.
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:164
ndarray::Array< double const, 2, 2 > getFisherMatrix()
Return the Fisher matrix (inverse of the covariance) of the parameters.
LeastSquares & operator=(LeastSquares const &)