Eigen-unsupported  3.4.90 (git rev a4098ac676528a83cfb73d4d26ce1b42ec05f47c)
MatrixSquareRoot.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2011, 2013 Jitse Niesen <jitse@maths.leeds.ac.uk>
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10#ifndef EIGEN_MATRIX_SQUARE_ROOT
11#define EIGEN_MATRIX_SQUARE_ROOT
12
14
15namespace Eigen {
16
17namespace internal {
18
19// pre: T.block(i,i,2,2) has complex conjugate eigenvalues
20// post: sqrtT.block(i,i,2,2) is square root of T.block(i,i,2,2)
21template <typename MatrixType, typename ResultType>
22void matrix_sqrt_quasi_triangular_2x2_diagonal_block(const MatrixType& T, Index i, ResultType& sqrtT)
23{
24 // TODO: This case (2-by-2 blocks with complex conjugate eigenvalues) is probably hidden somewhere
25 // in EigenSolver. If we expose it, we could call it directly from here.
26 typedef typename traits<MatrixType>::Scalar Scalar;
27 Matrix<Scalar,2,2> block = T.template block<2,2>(i,i);
28 EigenSolver<Matrix<Scalar,2,2> > es(block);
29 sqrtT.template block<2,2>(i,i)
30 = (es.eigenvectors() * es.eigenvalues().cwiseSqrt().asDiagonal() * es.eigenvectors().inverse()).real();
31}
32
33// pre: block structure of T is such that (i,j) is a 1x1 block,
34// all blocks of sqrtT to left of and below (i,j) are correct
35// post: sqrtT(i,j) has the correct value
36template <typename MatrixType, typename ResultType>
37void matrix_sqrt_quasi_triangular_1x1_off_diagonal_block(const MatrixType& T, Index i, Index j, ResultType& sqrtT)
38{
39 typedef typename traits<MatrixType>::Scalar Scalar;
40 Scalar tmp = (sqrtT.row(i).segment(i+1,j-i-1) * sqrtT.col(j).segment(i+1,j-i-1)).value();
41 sqrtT.coeffRef(i,j) = (T.coeff(i,j) - tmp) / (sqrtT.coeff(i,i) + sqrtT.coeff(j,j));
42}
43
44// similar to compute1x1offDiagonalBlock()
45template <typename MatrixType, typename ResultType>
46void matrix_sqrt_quasi_triangular_1x2_off_diagonal_block(const MatrixType& T, Index i, Index j, ResultType& sqrtT)
47{
48 typedef typename traits<MatrixType>::Scalar Scalar;
49 Matrix<Scalar,1,2> rhs = T.template block<1,2>(i,j);
50 if (j-i > 1)
51 rhs -= sqrtT.block(i, i+1, 1, j-i-1) * sqrtT.block(i+1, j, j-i-1, 2);
52 Matrix<Scalar,2,2> A = sqrtT.coeff(i,i) * Matrix<Scalar,2,2>::Identity();
53 A += sqrtT.template block<2,2>(j,j).transpose();
54 sqrtT.template block<1,2>(i,j).transpose() = A.fullPivLu().solve(rhs.transpose());
55}
56
57// similar to compute1x1offDiagonalBlock()
58template <typename MatrixType, typename ResultType>
59void matrix_sqrt_quasi_triangular_2x1_off_diagonal_block(const MatrixType& T, Index i, Index j, ResultType& sqrtT)
60{
61 typedef typename traits<MatrixType>::Scalar Scalar;
62 Matrix<Scalar,2,1> rhs = T.template block<2,1>(i,j);
63 if (j-i > 2)
64 rhs -= sqrtT.block(i, i+2, 2, j-i-2) * sqrtT.block(i+2, j, j-i-2, 1);
65 Matrix<Scalar,2,2> A = sqrtT.coeff(j,j) * Matrix<Scalar,2,2>::Identity();
66 A += sqrtT.template block<2,2>(i,i);
67 sqrtT.template block<2,1>(i,j) = A.fullPivLu().solve(rhs);
68}
69
70// solves the equation A X + X B = C where all matrices are 2-by-2
71template <typename MatrixType>
72void matrix_sqrt_quasi_triangular_solve_auxiliary_equation(MatrixType& X, const MatrixType& A, const MatrixType& B, const MatrixType& C)
73{
74 typedef typename traits<MatrixType>::Scalar Scalar;
75 Matrix<Scalar,4,4> coeffMatrix = Matrix<Scalar,4,4>::Zero();
76 coeffMatrix.coeffRef(0,0) = A.coeff(0,0) + B.coeff(0,0);
77 coeffMatrix.coeffRef(1,1) = A.coeff(0,0) + B.coeff(1,1);
78 coeffMatrix.coeffRef(2,2) = A.coeff(1,1) + B.coeff(0,0);
79 coeffMatrix.coeffRef(3,3) = A.coeff(1,1) + B.coeff(1,1);
80 coeffMatrix.coeffRef(0,1) = B.coeff(1,0);
81 coeffMatrix.coeffRef(0,2) = A.coeff(0,1);
82 coeffMatrix.coeffRef(1,0) = B.coeff(0,1);
83 coeffMatrix.coeffRef(1,3) = A.coeff(0,1);
84 coeffMatrix.coeffRef(2,0) = A.coeff(1,0);
85 coeffMatrix.coeffRef(2,3) = B.coeff(1,0);
86 coeffMatrix.coeffRef(3,1) = A.coeff(1,0);
87 coeffMatrix.coeffRef(3,2) = B.coeff(0,1);
88
89 Matrix<Scalar,4,1> rhs;
90 rhs.coeffRef(0) = C.coeff(0,0);
91 rhs.coeffRef(1) = C.coeff(0,1);
92 rhs.coeffRef(2) = C.coeff(1,0);
93 rhs.coeffRef(3) = C.coeff(1,1);
94
95 Matrix<Scalar,4,1> result;
96 result = coeffMatrix.fullPivLu().solve(rhs);
97
98 X.coeffRef(0,0) = result.coeff(0);
99 X.coeffRef(0,1) = result.coeff(1);
100 X.coeffRef(1,0) = result.coeff(2);
101 X.coeffRef(1,1) = result.coeff(3);
102}
103
104// similar to compute1x1offDiagonalBlock()
105template <typename MatrixType, typename ResultType>
106void matrix_sqrt_quasi_triangular_2x2_off_diagonal_block(const MatrixType& T, Index i, Index j, ResultType& sqrtT)
107{
108 typedef typename traits<MatrixType>::Scalar Scalar;
109 Matrix<Scalar,2,2> A = sqrtT.template block<2,2>(i,i);
110 Matrix<Scalar,2,2> B = sqrtT.template block<2,2>(j,j);
111 Matrix<Scalar,2,2> C = T.template block<2,2>(i,j);
112 if (j-i > 2)
113 C -= sqrtT.block(i, i+2, 2, j-i-2) * sqrtT.block(i+2, j, j-i-2, 2);
114 Matrix<Scalar,2,2> X;
115 matrix_sqrt_quasi_triangular_solve_auxiliary_equation(X, A, B, C);
116 sqrtT.template block<2,2>(i,j) = X;
117}
118
119// pre: T is quasi-upper-triangular and sqrtT is a zero matrix of the same size
120// post: the diagonal blocks of sqrtT are the square roots of the diagonal blocks of T
121template <typename MatrixType, typename ResultType>
122void matrix_sqrt_quasi_triangular_diagonal(const MatrixType& T, ResultType& sqrtT)
123{
124 using std::sqrt;
125 const Index size = T.rows();
126 for (Index i = 0; i < size; i++) {
127 if (i == size - 1 || T.coeff(i+1, i) == 0) {
128 eigen_assert(T(i,i) >= 0);
129 sqrtT.coeffRef(i,i) = sqrt(T.coeff(i,i));
130 }
131 else {
132 matrix_sqrt_quasi_triangular_2x2_diagonal_block(T, i, sqrtT);
133 ++i;
134 }
135 }
136}
137
138// pre: T is quasi-upper-triangular and diagonal blocks of sqrtT are square root of diagonal blocks of T.
139// post: sqrtT is the square root of T.
140template <typename MatrixType, typename ResultType>
141void matrix_sqrt_quasi_triangular_off_diagonal(const MatrixType& T, ResultType& sqrtT)
142{
143 const Index size = T.rows();
144 for (Index j = 1; j < size; j++) {
145 if (T.coeff(j, j-1) != 0) // if T(j-1:j, j-1:j) is a 2-by-2 block
146 continue;
147 for (Index i = j-1; i >= 0; i--) {
148 if (i > 0 && T.coeff(i, i-1) != 0) // if T(i-1:i, i-1:i) is a 2-by-2 block
149 continue;
150 bool iBlockIs2x2 = (i < size - 1) && (T.coeff(i+1, i) != 0);
151 bool jBlockIs2x2 = (j < size - 1) && (T.coeff(j+1, j) != 0);
152 if (iBlockIs2x2 && jBlockIs2x2)
153 matrix_sqrt_quasi_triangular_2x2_off_diagonal_block(T, i, j, sqrtT);
154 else if (iBlockIs2x2 && !jBlockIs2x2)
155 matrix_sqrt_quasi_triangular_2x1_off_diagonal_block(T, i, j, sqrtT);
156 else if (!iBlockIs2x2 && jBlockIs2x2)
157 matrix_sqrt_quasi_triangular_1x2_off_diagonal_block(T, i, j, sqrtT);
158 else if (!iBlockIs2x2 && !jBlockIs2x2)
159 matrix_sqrt_quasi_triangular_1x1_off_diagonal_block(T, i, j, sqrtT);
160 }
161 }
162}
163
164} // end of namespace internal
165
181template <typename MatrixType, typename ResultType>
182void matrix_sqrt_quasi_triangular(const MatrixType &arg, ResultType &result)
183{
184 eigen_assert(arg.rows() == arg.cols());
185 result.resize(arg.rows(), arg.cols());
186 internal::matrix_sqrt_quasi_triangular_diagonal(arg, result);
187 internal::matrix_sqrt_quasi_triangular_off_diagonal(arg, result);
188}
189
190
205template <typename MatrixType, typename ResultType>
206void matrix_sqrt_triangular(const MatrixType &arg, ResultType &result)
207{
208 using std::sqrt;
209 typedef typename MatrixType::Scalar Scalar;
210
211 eigen_assert(arg.rows() == arg.cols());
212
213 // Compute square root of arg and store it in upper triangular part of result
214 // This uses that the square root of triangular matrices can be computed directly.
215 result.resize(arg.rows(), arg.cols());
216 for (Index i = 0; i < arg.rows(); i++) {
217 result.coeffRef(i,i) = sqrt(arg.coeff(i,i));
218 }
219 for (Index j = 1; j < arg.cols(); j++) {
220 for (Index i = j-1; i >= 0; i--) {
221 // if i = j-1, then segment has length 0 so tmp = 0
222 Scalar tmp = (result.row(i).segment(i+1,j-i-1) * result.col(j).segment(i+1,j-i-1)).value();
223 // denominator may be zero if original matrix is singular
224 result.coeffRef(i,j) = (arg.coeff(i,j) - tmp) / (result.coeff(i,i) + result.coeff(j,j));
225 }
226 }
227}
228
229
230namespace internal {
231
239template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
240struct matrix_sqrt_compute
241{
249 template <typename ResultType> static void run(const MatrixType &arg, ResultType &result);
250};
251
252
253// ********** Partial specialization for real matrices **********
254
255template <typename MatrixType>
256struct matrix_sqrt_compute<MatrixType, 0>
257{
258 typedef typename MatrixType::PlainObject PlainType;
259 template <typename ResultType>
260 static void run(const MatrixType &arg, ResultType &result)
261 {
262 eigen_assert(arg.rows() == arg.cols());
263
264 // Compute Schur decomposition of arg
265 const RealSchur<PlainType> schurOfA(arg);
266 const PlainType& T = schurOfA.matrixT();
267 const PlainType& U = schurOfA.matrixU();
268
269 // Compute square root of T
270 PlainType sqrtT = PlainType::Zero(arg.rows(), arg.cols());
272
273 // Compute square root of arg
274 result = U * sqrtT * U.adjoint();
275 }
276};
277
278
279// ********** Partial specialization for complex matrices **********
280
281template <typename MatrixType>
282struct matrix_sqrt_compute<MatrixType, 1>
283{
284 typedef typename MatrixType::PlainObject PlainType;
285 template <typename ResultType>
286 static void run(const MatrixType &arg, ResultType &result)
287 {
288 eigen_assert(arg.rows() == arg.cols());
289
290 // Compute Schur decomposition of arg
291 const ComplexSchur<PlainType> schurOfA(arg);
292 const PlainType& T = schurOfA.matrixT();
293 const PlainType& U = schurOfA.matrixU();
294
295 // Compute square root of T
296 PlainType sqrtT;
297 matrix_sqrt_triangular(T, sqrtT);
298
299 // Compute square root of arg
300 result = U * (sqrtT.template triangularView<Upper>() * U.adjoint());
301 }
302};
303
304} // end namespace internal
305
318template<typename Derived> class MatrixSquareRootReturnValue
319: public ReturnByValue<MatrixSquareRootReturnValue<Derived> >
320{
321 protected:
322 typedef typename internal::ref_selector<Derived>::type DerivedNested;
323
324 public:
330 explicit MatrixSquareRootReturnValue(const Derived& src) : m_src(src) { }
331
337 template <typename ResultType>
338 inline void evalTo(ResultType& result) const
339 {
340 typedef typename internal::nested_eval<Derived, 10>::type DerivedEvalType;
341 typedef typename internal::remove_all<DerivedEvalType>::type DerivedEvalTypeClean;
342 DerivedEvalType tmp(m_src);
343 internal::matrix_sqrt_compute<DerivedEvalTypeClean>::run(tmp, result);
344 }
345
346 Index rows() const { return m_src.rows(); }
347 Index cols() const { return m_src.cols(); }
348
349 protected:
350 const DerivedNested m_src;
351};
352
353namespace internal {
354template<typename Derived>
355struct traits<MatrixSquareRootReturnValue<Derived> >
356{
357 typedef typename Derived::PlainObject ReturnType;
358};
359}
360
361template <typename Derived>
362const MatrixSquareRootReturnValue<Derived> MatrixBase<Derived>::sqrt() const
363{
364 eigen_assert(rows() == cols());
365 return MatrixSquareRootReturnValue<Derived>(derived());
366}
367
368} // end namespace Eigen
369
370#endif // EIGEN_MATRIX_FUNCTION
const MatrixSquareRootReturnValue< Derived > sqrt() const
Definition: MatrixSquareRoot.h:362
Proxy for the matrix square root of some matrix (expression).
Definition: MatrixSquareRoot.h:320
void evalTo(ResultType &result) const
Compute the matrix square root.
Definition: MatrixSquareRoot.h:338
MatrixSquareRootReturnValue(const Derived &src)
Constructor.
Definition: MatrixSquareRoot.h:330
void matrix_sqrt_quasi_triangular(const MatrixType &arg, ResultType &result)
Compute matrix square root of quasi-triangular matrix.
Definition: MatrixSquareRoot.h:182
void matrix_sqrt_triangular(const MatrixType &arg, ResultType &result)
Compute matrix square root of triangular matrix.
Definition: MatrixSquareRoot.h:206
Namespace containing all symbols from the Eigen library.
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_real_op< typename Derived::Scalar >, const Derived > real(const Eigen::ArrayBase< Derived > &x)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_arg_op< typename Derived::Scalar >, const Derived > arg(const Eigen::ArrayBase< Derived > &x)