LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
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
31namespace lsst {
32namespace afw {
33namespace math {
34
67class LeastSquares final {
68public:
69 class Impl;
70
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
314
315 // Need to define dtor in source file so it can see Impl declaration.
317
318private:
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
339};
340} // namespace math
341} // namespace afw
342} // namespace lsst
343
344#endif // !LSST_AFW_MATH_LeastSquares_h_INCLUDED
char * data
Definition BaseRecord.cc:61
Basic LSST definitions.
Solver for linear least-squares problems.
LeastSquares & operator=(LeastSquares &&)
void setDesignMatrix(ndarray::Array< T1, 2, C1 > const &design)
Reset the design matrix given as an ndarray; dimension and data are not changed.
ndarray::Array< double const, 1, 1 > getSolution()
Return the vector solution to the least squares problem.
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.
ndarray::Array< double const, 2, 2 > getFisherMatrix()
Return the Fisher matrix (inverse of the covariance) of the parameters.
int getRank() const
Return the rank of the problem (number of nonzero Eigenvalues).
ndarray::Array< double const, 2, 2 > getCovariance()
Return the covariance matrix of the least squares problem.
LeastSquares & operator=(LeastSquares const &)
Factorization
Private implementation; forward-declared publicly so we can inherit from it in .cc.
@ NORMAL_EIGENSYSTEM
Use the normal equations with a symmetric Eigensystem decomposition.
@ NORMAL_CHOLESKY
Use the normal equations with a Cholesky decomposition.
@ DIRECT_SVD
Use a thin singular value decomposition of the design matrix.
LeastSquares(Factorization factorization, int dimension)
Construct a least-squares object for the given factorization and dimensionality.
int getDimension() const
Return the number of parameters.
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.
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.
Factorization getFactorization() const
Retun the type of factorization used by the solver.
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 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.
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.
double getThreshold() const
Get the threshold used to determine when to truncate Eigenvalues.
LeastSquares(LeastSquares &&)
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.
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.
void setThreshold(double threshold)
Set the threshold used to determine when to truncate Eigenvalues.
void setDesignMatrix(Eigen::MatrixBase< D1 > const &design)
Reset the design matrix given as an Eigen object; dimension and data are not changed.
LeastSquares(LeastSquares const &)
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.