1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
|
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra. Eigen itself is part of the KDE project.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_CHOLESKY_WITHOUT_SQUARE_ROOT_H
#define EIGEN_CHOLESKY_WITHOUT_SQUARE_ROOT_H
/** \class CholeskyWithoutSquareRoot
*
* \brief Robust Cholesky decomposition of a matrix and associated features
*
* \param MatrixType the type of the matrix of which we are computing the Cholesky decomposition
*
* This class performs a Cholesky decomposition without square root of a symmetric, positive definite
* matrix A such that A = L D L^* = U^* D U, where L is lower triangular with a unit diagonal and D is a diagonal
* matrix.
*
* Compared to a standard Cholesky decomposition, avoiding the square roots allows for faster and more
* stable computation.
*
* Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
* the strict lower part does not have to store correct values.
*
* \sa class Cholesky
*/
template<typename MatrixType> class CholeskyWithoutSquareRoot
{
public:
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
CholeskyWithoutSquareRoot(const MatrixType& matrix)
: m_matrix(matrix.rows(), matrix.cols())
{
compute(matrix);
}
/** \returns the lower triangular matrix L */
Extract<MatrixType, UnitLower> matrixL(void) const
{
return m_matrix;
}
/** \returns the coefficients of the diagonal matrix D */
DiagonalCoeffs<MatrixType> vectorD(void) const
{
return m_matrix.diagonal();
}
/** \returns whether the matrix is positive definite */
bool isPositiveDefinite(void) const
{
// FIXME is it really correct ?
return m_matrix.diagonal().real().minCoeff() > RealScalar(0);
}
template<typename Derived>
typename Derived::Eval solve(MatrixBase<Derived> &b);
void compute(const MatrixType& matrix);
protected:
/** \internal
* Used to compute and store the cholesky decomposition A = L D L^* = U^* D U.
* The strict upper part is used during the decomposition, the strict lower
* part correspond to the coefficients of L (its diagonal is equal to 1 and
* is not stored), and the diagonal entries correspond to D.
*/
MatrixType m_matrix;
};
/** Compute / recompute the Cholesky decomposition A = L D L^* = U^* D U of \a matrix
*/
template<typename MatrixType>
void CholeskyWithoutSquareRoot<MatrixType>::compute(const MatrixType& a)
{
assert(a.rows()==a.cols());
const int size = a.rows();
m_matrix.resize(size, size);
// Note that, in this algorithm the rows of the strict upper part of m_matrix is used to store
// column vector, thus the strange .conjugate() and .transpose()...
m_matrix.row(0) = a.row(0).conjugate();
m_matrix.col(0).end(size-1) = m_matrix.row(0).end(size-1) / m_matrix.coeff(0,0);
for (int j = 1; j < size; ++j)
{
RealScalar tmp = ei_real(a.coeff(j,j) - (m_matrix.row(j).start(j) * m_matrix.col(j).start(j).conjugate()).coeff(0,0));
m_matrix.coeffRef(j,j) = tmp;
int endSize = size-j-1;
if (endSize>0)
{
m_matrix.row(j).end(endSize) = a.row(j).end(endSize).conjugate()
- (m_matrix.block(j+1,0, endSize, j) * m_matrix.col(j).start(j).conjugate()).transpose();
m_matrix.col(j).end(endSize) = m_matrix.row(j).end(endSize) / tmp;
}
}
}
/** \returns the solution of \f$ A x = b \f$ using the current decomposition of A.
* In other words, it returns \f$ A^{-1} b \f$ computing
* \f$ {L^{*}}^{-1} D^{-1} L^{-1} b \f$ from right to left.
* \param vecB the vector \f$ b \f$ (or an array of vectors)
*/
template<typename MatrixType>
template<typename Derived>
typename Derived::Eval CholeskyWithoutSquareRoot<MatrixType>::solve(MatrixBase<Derived> &vecB)
{
const int size = m_matrix.rows();
ei_assert(size==vecB.size());
return m_matrix.adjoint().template extract<UnitUpper>()
.inverseProduct(
(matrixL()
.inverseProduct(vecB))
.cwiseQuotient(m_matrix.diagonal())
);
}
#endif // EIGEN_CHOLESKY_WITHOUT_SQUARE_ROOT_H
|
