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.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008, 2009, 2010, 2011 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 #include "Eigen/Eigenvalues"
26 #include "Eigen/SVD"
27 #include "Eigen/Cholesky"
28 #include "boost/format.hpp"
29 #include "boost/make_shared.hpp"
30 
32 #include "lsst/pex/exceptions.h"
33 #include "lsst/pex/logging.h"
34 
35 namespace lsst { namespace afw { namespace math {
36 
38 public:
39 
40  enum StateFlags {
43  RHS_VECTOR = 0x004,
44  SOLUTION_ARRAY = 0x008,
48  };
49 
50  int state;
51  int dimension;
52  int rank;
55  double threshold;
56 
57  Eigen::MatrixXd design;
58  Eigen::VectorXd data;
59  Eigen::MatrixXd fisher;
60  Eigen::VectorXd rhs;
61 
65 
67 
68  template <typename D>
69  void setRank(Eigen::MatrixBase<D> const & values) {
70  double cond = threshold * values[0];
71  if (cond <= 0.0) {
72  rank = 0;
73  } else {
74  for (rank = dimension; (rank > 1) && (values[rank-1] < cond); --rank);
75  }
76  }
77 
78  void ensure(int desired) {
80  if (desired & FULL_FISHER_MATRIX) desired |= LOWER_FISHER_MATRIX;
81  int toAdd = ~state & desired;
82  if (toAdd & LOWER_FISHER_MATRIX) {
83  assert(state & DESIGN_AND_DATA);
84  fisher = Eigen::MatrixXd::Zero(design.cols(), design.cols());
85  fisher.selfadjointView<Eigen::Lower>().rankUpdate(design.adjoint());
86  }
87  if (toAdd & FULL_FISHER_MATRIX) {
88  fisher.triangularView<Eigen::StrictlyUpper>() = fisher.adjoint();
89  }
90  if (toAdd & RHS_VECTOR) {
91  assert(state & DESIGN_AND_DATA);
92  rhs = design.adjoint() * data;
93  }
94  if (toAdd & SOLUTION_ARRAY) {
97  }
98  if (toAdd & COVARIANCE_ARRAY) {
101  }
102  if (toAdd & DIAGNOSTIC_ARRAY) {
105  }
106  state |= toAdd;
107  }
108 
109  virtual void factor() = 0;
110 
111  virtual void updateRank() = 0;
112 
113  virtual void updateSolution() = 0;
114  virtual void updateCovariance() = 0;
115 
116  virtual void updateDiagnostic() = 0;
117 
118  Impl(int dimension_, double threshold_=std::numeric_limits<double>::epsilon()) :
119  state(0), dimension(dimension_), rank(dimension_), threshold(threshold_),
120  log("afw.math.LeastSquares")
121  {}
122 
123  virtual ~Impl() {}
124 };
125 
126 namespace {
127 
128 class EigensystemSolver : public LeastSquares::Impl {
129 public:
130 
131  explicit EigensystemSolver(int dimension) :
132  Impl(dimension), _eig(dimension), _svd(), _tmp(dimension)
133  {}
134 
135  virtual void factor() {
136  ensure(LOWER_FISHER_MATRIX | RHS_VECTOR);
137  _eig.compute(fisher);
138  if (_eig.info() == Eigen::Success) {
139  setRank(_eig.eigenvalues().reverse());
140  log.debug<5>("SelfAdjointEigenSolver succeeded: dimension=%d, rank=%d", dimension, rank);
141  } else {
142  // Note that the fallback is using SVD of the Fisher to compute the Eigensystem, because those
143  // are the same for a symmetric matrix; this is very different from doing a direct SVD of
144  // the design matrix.
145  ensure(FULL_FISHER_MATRIX);
146  _svd.compute(fisher, Eigen::ComputeFullU); // Matrix is symmetric, so V == U == eigenvectors
147  setRank(_svd.singularValues());
148  log.debug<5>(
149  "SelfAdjointEigenSolver failed; falling back to equivalent SVD: dimension=%d, rank=%d",
150  dimension, rank
151  );
152  }
153  }
154 
155  virtual void updateRank() {
156  if (_eig.info() == Eigen::Success) {
157  setRank(_eig.eigenvalues().reverse());
158  } else {
159  setRank(_svd.singularValues());
160  }
161  }
162 
163  virtual void updateDiagnostic() {
164  if (whichDiagnostic == LeastSquares::NORMAL_CHOLESKY) {
165  throw LSST_EXCEPT(
166  pex::exceptions::LogicError,
167  "Cannot compute NORMAL_CHOLESKY diagnostic from NORMAL_EIGENSYSTEM factorization."
168  );
169  }
170  if (_eig.info() == Eigen::Success) {
171  diagnostic.asEigen() = _eig.eigenvalues().reverse();
172  } else {
173  diagnostic.asEigen() = _svd.singularValues();
174  }
175  if (whichDiagnostic == LeastSquares::DIRECT_SVD) {
176  diagnostic.asEigen<Eigen::ArrayXpr>() = diagnostic.asEigen<Eigen::ArrayXpr>().sqrt();
177  }
178  }
179 
180  virtual void updateSolution() {
181  if (rank == 0) {
182  solution.asEigen().setZero();
183  return;
184  }
185  if (_eig.info() == Eigen::Success) {
186  _tmp.head(rank) = _eig.eigenvectors().rightCols(rank).adjoint() * rhs;
187  _tmp.head(rank).array() /= _eig.eigenvalues().tail(rank).array();
188  solution.asEigen() = _eig.eigenvectors().rightCols(rank) * _tmp.head(rank);
189  } else {
190  _tmp.head(rank) = _svd.matrixU().leftCols(rank).adjoint() * rhs;
191  _tmp.head(rank).array() /= _svd.singularValues().head(rank).array();
192  solution.asEigen() = _svd.matrixU().leftCols(rank) * _tmp.head(rank);
193  }
194  }
195 
196  virtual void updateCovariance() {
197  if (rank == 0) {
198  covariance.asEigen().setZero();
199  return;
200  }
201  if (_eig.info() == Eigen::Success) {
202  covariance.asEigen() =
203  _eig.eigenvectors().rightCols(rank)
204  * _eig.eigenvalues().tail(rank).array().inverse().matrix().asDiagonal()
205  * _eig.eigenvectors().rightCols(rank).adjoint();
206  } else {
207  covariance.asEigen() =
208  _svd.matrixU().leftCols(rank)
209  * _svd.singularValues().head(rank).array().inverse().matrix().asDiagonal()
210  * _svd.matrixU().leftCols(rank).adjoint();
211  }
212  }
213 
214 private:
215  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> _eig;
216  Eigen::JacobiSVD<Eigen::MatrixXd> _svd; // only used if Eigendecomposition fails, should be very rare
217  Eigen::VectorXd _tmp;
218 };
219 
220 class CholeskySolver : public LeastSquares::Impl {
221 public:
222 
223  explicit CholeskySolver(int dimension) : Impl(dimension, 0.0), _ldlt(dimension) {}
224 
225  virtual void factor() {
226  ensure(LOWER_FISHER_MATRIX | RHS_VECTOR);
227  _ldlt.compute(fisher);
228  }
229 
230  virtual void updateRank() {}
231 
232  virtual void updateDiagnostic() {
233  if (whichDiagnostic != LeastSquares::NORMAL_CHOLESKY) {
234  throw LSST_EXCEPT(
235  pex::exceptions::LogicError,
236  "Can only compute NORMAL_CHOLESKY diagnostic from NORMAL_CHOLESKY factorization."
237  );
238  }
239  diagnostic.asEigen() = _ldlt.vectorD();
240  }
241 
242  virtual void updateSolution() { solution.asEigen() = _ldlt.solve(rhs); }
243 
244  virtual void updateCovariance() {
245  ndarray::EigenView<double,2,2> cov(covariance);
246  cov.setIdentity();
247  cov = _ldlt.solve(cov);
248  }
249 
250 private:
251  Eigen::LDLT<Eigen::MatrixXd> _ldlt;
252 };
253 
254 class SvdSolver : public LeastSquares::Impl {
255 public:
256 
257  explicit SvdSolver(int dimension) : Impl(dimension), _svd(), _tmp(dimension) {}
258 
259  virtual void factor() {
260  if (!(state & DESIGN_AND_DATA)) {
261  throw LSST_EXCEPT(
262  pex::exceptions::InvalidParameterError,
263  "Cannot initialize DIRECT_SVD solver with normal equations."
264  );
265  }
266  _svd.compute(design, Eigen::ComputeThinU | Eigen::ComputeThinV);
267  setRank(_svd.singularValues());
268  log.debug<5>("Using direct SVD method; dimension=%d, rank=%d", dimension, rank);
269  }
270 
271  virtual void updateRank() { setRank(_svd.singularValues()); }
272 
273  virtual void updateDiagnostic() {
274  switch (whichDiagnostic) {
276  diagnostic.asEigen<Eigen::ArrayXpr>() = _svd.singularValues().array().square();
277  break;
279  throw LSST_EXCEPT(
280  pex::exceptions::LogicError,
281  "Can only compute NORMAL_CHOLESKY diagnostic from DIRECT_SVD factorization."
282  );
284  diagnostic.asEigen() = _svd.singularValues();
285  break;
286  }
287  }
288 
289  virtual void updateSolution() {
290  if (rank == 0) {
291  solution.asEigen().setZero();
292  return;
293  }
294  _tmp.head(rank) = _svd.matrixU().leftCols(rank).adjoint() * data;
295  _tmp.head(rank).array() /= _svd.singularValues().head(rank).array();
296  solution.asEigen() = _svd.matrixV().leftCols(rank) * _tmp.head(rank);
297  }
298 
299  virtual void updateCovariance() {
300  if (rank == 0) {
301  covariance.asEigen().setZero();
302  return;
303  }
304  covariance.asEigen() =
305  _svd.matrixV().leftCols(rank)
306  * _svd.singularValues().head(rank).array().inverse().square().matrix().asDiagonal()
307  * _svd.matrixV().leftCols(rank).adjoint();
308  }
309 
310 private:
311  Eigen::JacobiSVD<Eigen::MatrixXd> _svd;
312  Eigen::VectorXd _tmp;
313 };
314 
315 } // anonymous
316 
317 void LeastSquares::setThreshold(double threshold) {
318  _impl->threshold = threshold;
319  _impl->state &= ~Impl::SOLUTION_ARRAY;
320  _impl->state &= ~Impl::COVARIANCE_ARRAY;
321  _impl->updateRank();
322 }
323 
324 double LeastSquares::getThreshold() const { return _impl->threshold; }
325 
327  _impl->ensure(Impl::SOLUTION_ARRAY);
328  return _impl->solution;
329 }
330 
332  _impl->ensure(Impl::COVARIANCE_ARRAY);
333  return _impl->covariance;
334 }
335 
338  // Wrap the Eigen::MatrixXd in an ndarray::Array, using _impl as the reference-counted owner.
339  // Doesn't matter if we swap strides, because it's symmetric.
340  return ndarray::external(
341  _impl->fisher.data(),
342  ndarray::makeVector(_impl->dimension, _impl->dimension),
343  ndarray::makeVector(_impl->dimension, 1),
344  _impl
345  );
346 }
347 
349  if (_impl->whichDiagnostic != factorization) {
350  _impl->state &= ~Impl::DIAGNOSTIC_ARRAY;
351  _impl->whichDiagnostic = factorization;
352  }
353  _impl->ensure(Impl::DIAGNOSTIC_ARRAY);
354  return _impl->diagnostic;
355 }
356 
357 int LeastSquares::getDimension() const { return _impl->dimension; }
358 
359 int LeastSquares::getRank() const { return _impl->rank; }
360 
362 
363 LeastSquares::LeastSquares(Factorization factorization, int dimension) {
364  switch (factorization) {
365  case NORMAL_EIGENSYSTEM:
366  _impl = boost::make_shared<EigensystemSolver>(dimension);
367  break;
368  case NORMAL_CHOLESKY:
369  _impl = boost::make_shared<CholeskySolver>(dimension);
370  break;
371  case DIRECT_SVD:
372  _impl = boost::make_shared<SvdSolver>(dimension);
373  break;
374  }
375  _impl->factorization = factorization;
376 }
377 
379 
380 Eigen::MatrixXd & LeastSquares::_getDesignMatrix() { return _impl->design; }
381 Eigen::VectorXd & LeastSquares::_getDataVector() { return _impl->data; }
382 
383 Eigen::MatrixXd & LeastSquares::_getFisherMatrix() { return _impl->fisher; }
384 Eigen::VectorXd & LeastSquares::_getRhsVector() { return _impl->rhs; }
385 
386 void LeastSquares::_factor(bool haveNormalEquations) {
387  if (haveNormalEquations) {
388  if (_getFisherMatrix().rows() != _impl->dimension) {
389  throw LSST_EXCEPT(
390  pex::exceptions::InvalidParameterError,
391  (boost::format("Number of rows of Fisher matrix (%d) does not match"
392  " dimension of LeastSquares solver.")
393  % _getFisherMatrix().rows() % _impl->dimension).str()
394  );
395  }
396  if (_getFisherMatrix().cols() != _impl->dimension) {
397  throw LSST_EXCEPT(
398  pex::exceptions::InvalidParameterError,
399  (boost::format("Number of columns of Fisher matrix (%d) does not match"
400  " dimension of LeastSquares solver.")
401  % _getFisherMatrix().cols() % _impl->dimension).str()
402  );
403  }
404  if (_getRhsVector().size() != _impl->dimension) {
405  throw LSST_EXCEPT(
406  pex::exceptions::InvalidParameterError,
407  (boost::format("Number of elements in RHS vector (%d) does not match"
408  " dimension of LeastSquares solver.")
409  % _getRhsVector().size() % _impl->dimension).str()
410  );
411  }
413  } else {
414  if (_getDesignMatrix().cols() != _impl->dimension) {
415  throw LSST_EXCEPT(
416  pex::exceptions::InvalidParameterError,
417  "Number of columns of design matrix does not match dimension of LeastSquares solver."
418  );
419  }
420  if (_getDesignMatrix().rows() != _getDataVector().size()) {
421  throw LSST_EXCEPT(
422  pex::exceptions::InvalidParameterError,
423  (boost::format("Number of rows of design matrix (%d) does not match number of "
424  "data points (%d)") % _getDesignMatrix().rows() % _getDataVector().size()
425  ).str()
426  );
427  }
428  if (_getDesignMatrix().cols() > _getDataVector().size()) {
429  throw LSST_EXCEPT(
430  pex::exceptions::InvalidParameterError,
431  (boost::format("Number of columns of design matrix (%d) must be smaller than number of "
432  "data points (%d)") % _getDesignMatrix().cols() % _getDataVector().size()
433  ).str()
434  );
435  }
436  _impl->state = Impl::DESIGN_AND_DATA;
437  }
438  _impl->factor();
439 }
440 
441 }}} // namespace lsst::afw::math
void setRank(Eigen::MatrixBase< D > const &values)
Definition: LeastSquares.cc:69
int getDimension() const
Return the number of parameters.
ndarray::Array< double const, 2, 2 > getCovariance()
Return the covariance matrix of the least squares problem.
Eigen::VectorXd _tmp
Eigen::LDLT< Eigen::MatrixXd > _ldlt
Eigen::VectorXd & _getDataVector()
Eigen::SelfAdjointEigenSolver< Eigen::MatrixXd > _eig
Eigen3 view into an ndarray::Array.
Definition: eigen.h:232
Eigen::VectorXd & _getRhsVector()
detail::ExternalInitializer< T, N, Owner > external(T *data, Vector< int, N > const &shape, Vector< int, N > const &strides, Owner const &owner)
Create an expression that initializes an Array with externally allocated memory.
ndarray::Array< double, 1, 1 > diagnostic
Definition: LeastSquares.cc:64
int getRank() const
Return the rank of the problem (number of nonzero Eigenvalues).
def log
Definition: log.py:85
ndarray::Array< double, 1, 1 > solution
Definition: LeastSquares.cc:62
Use a thin singular value decomposition of the design matrix.
Definition: LeastSquares.h:87
Impl(int dimension_, double threshold_=std::numeric_limits< double >::epsilon())
Use the normal equations with a Cholesky decomposition.
Definition: LeastSquares.h:79
Vector< T, N > makeVector(T v1, T v2,..., T vN)
Variadic constructor for Vector.
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.
double getThreshold() const
Get the threshold used to determine when to truncate Eigenvalues.
detail::SimpleInitializer< N > allocate(Vector< int, N > const &shape)
Create an expression that allocates uninitialized memory for an array.
bool isEmpty() const
Return true if the array has a null data point.
Definition: ArrayBase.h:120
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
Factorization
Private implementation; forward-declared publicly so we can inherit from it in .cc.
Definition: LeastSquares.h:70
void _factor(bool haveNormalEquations)
Use the normal equations with a symmetric Eigensystem decomposition.
Definition: LeastSquares.h:71
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()
Eigen::JacobiSVD< Eigen::MatrixXd > _svd
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...
ndarray::Array< double, 2, 2 > covariance
Definition: LeastSquares.cc:63
Include files required for standard LSST Exception handling.
LeastSquares(Factorization factorization, int dimension)
Construct a least-squares object for the given factorization and dimensionality.