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