MatrixFunction.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2009-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_FUNCTION
11 #define EIGEN_MATRIX_FUNCTION
12 
13 #include "StemFunction.h"
14 
15 
16 namespace Eigen {
17 
18 namespace internal {
19 
21 static const float matrix_function_separation = 0.1f;
22 
29 template <typename MatrixType>
30 class MatrixFunctionAtomic
31 {
32  public:
33 
34  typedef typename MatrixType::Scalar Scalar;
35  typedef typename stem_function<Scalar>::type StemFunction;
36 
40  MatrixFunctionAtomic(StemFunction f) : m_f(f) { }
41 
46  MatrixType compute(const MatrixType& A);
47 
48  private:
49  StemFunction* m_f;
50 };
51 
52 template <typename MatrixType>
53 typename NumTraits<typename MatrixType::Scalar>::Real matrix_function_compute_mu(const MatrixType& A)
54 {
55  typedef typename plain_col_type<MatrixType>::type VectorType;
56  typename MatrixType::Index rows = A.rows();
57  const MatrixType N = MatrixType::Identity(rows, rows) - A;
58  VectorType e = VectorType::Ones(rows);
59  N.template triangularView<Upper>().solveInPlace(e);
60  return e.cwiseAbs().maxCoeff();
61 }
62 
63 template <typename MatrixType>
64 MatrixType MatrixFunctionAtomic<MatrixType>::compute(const MatrixType& A)
65 {
66  // TODO: Use that A is upper triangular
67  typedef typename NumTraits<Scalar>::Real RealScalar;
68  typedef typename MatrixType::Index Index;
69  Index rows = A.rows();
70  Scalar avgEival = A.trace() / Scalar(RealScalar(rows));
71  MatrixType Ashifted = A - avgEival * MatrixType::Identity(rows, rows);
72  RealScalar mu = matrix_function_compute_mu(Ashifted);
73  MatrixType F = m_f(avgEival, 0) * MatrixType::Identity(rows, rows);
74  MatrixType P = Ashifted;
75  MatrixType Fincr;
76  for (Index s = 1; s < 1.1 * rows + 10; s++) { // upper limit is fairly arbitrary
77  Fincr = m_f(avgEival, static_cast<int>(s)) * P;
78  F += Fincr;
79  P = Scalar(RealScalar(1.0/(s + 1))) * P * Ashifted;
80 
81  // test whether Taylor series converged
82  const RealScalar F_norm = F.cwiseAbs().rowwise().sum().maxCoeff();
83  const RealScalar Fincr_norm = Fincr.cwiseAbs().rowwise().sum().maxCoeff();
84  if (Fincr_norm < NumTraits<Scalar>::epsilon() * F_norm) {
85  RealScalar delta = 0;
86  RealScalar rfactorial = 1;
87  for (Index r = 0; r < rows; r++) {
88  RealScalar mx = 0;
89  for (Index i = 0; i < rows; i++)
90  mx = (std::max)(mx, std::abs(m_f(Ashifted(i, i) + avgEival, static_cast<int>(s+r))));
91  if (r != 0)
92  rfactorial *= RealScalar(r);
93  delta = (std::max)(delta, mx / rfactorial);
94  }
95  const RealScalar P_norm = P.cwiseAbs().rowwise().sum().maxCoeff();
96  if (mu * delta * P_norm < NumTraits<Scalar>::epsilon() * F_norm) // series converged
97  break;
98  }
99  }
100  return F;
101 }
102 
108 template <typename Index, typename ListOfClusters>
109 typename ListOfClusters::iterator matrix_function_find_cluster(Index key, ListOfClusters& clusters)
110 {
111  typename std::list<Index>::iterator j;
112  for (typename ListOfClusters::iterator i = clusters.begin(); i != clusters.end(); ++i) {
113  j = std::find(i->begin(), i->end(), key);
114  if (j != i->end())
115  return i;
116  }
117  return clusters.end();
118 }
119 
131 template <typename EivalsType, typename Cluster>
132 void matrix_function_partition_eigenvalues(const EivalsType& eivals, std::list<Cluster>& clusters)
133 {
134  typedef typename EivalsType::Index Index;
135  for (Index i=0; i<eivals.rows(); ++i) {
136  // Find cluster containing i-th ei'val, adding a new cluster if necessary
137  typename std::list<Cluster>::iterator qi = matrix_function_find_cluster(i, clusters);
138  if (qi == clusters.end()) {
139  Cluster l;
140  l.push_back(i);
141  clusters.push_back(l);
142  qi = clusters.end();
143  --qi;
144  }
145 
146  // Look for other element to add to the set
147  for (Index j=i+1; j<eivals.rows(); ++j) {
148  if (abs(eivals(j) - eivals(i)) <= matrix_function_separation
149  && std::find(qi->begin(), qi->end(), j) == qi->end()) {
150  typename std::list<Cluster>::iterator qj = matrix_function_find_cluster(j, clusters);
151  if (qj == clusters.end()) {
152  qi->push_back(j);
153  } else {
154  qi->insert(qi->end(), qj->begin(), qj->end());
155  clusters.erase(qj);
156  }
157  }
158  }
159  }
160 }
161 
163 template <typename ListOfClusters, typename Index>
164 void matrix_function_compute_cluster_size(const ListOfClusters& clusters, Matrix<Index, Dynamic, 1>& clusterSize)
165 {
166  const Index numClusters = static_cast<Index>(clusters.size());
167  clusterSize.setZero(numClusters);
168  Index clusterIndex = 0;
169  for (typename ListOfClusters::const_iterator cluster = clusters.begin(); cluster != clusters.end(); ++cluster) {
170  clusterSize[clusterIndex] = cluster->size();
171  ++clusterIndex;
172  }
173 }
174 
176 template <typename VectorType>
177 void matrix_function_compute_block_start(const VectorType& clusterSize, VectorType& blockStart)
178 {
179  blockStart.resize(clusterSize.rows());
180  blockStart(0) = 0;
181  for (typename VectorType::Index i = 1; i < clusterSize.rows(); i++) {
182  blockStart(i) = blockStart(i-1) + clusterSize(i-1);
183  }
184 }
185 
187 template <typename EivalsType, typename ListOfClusters, typename VectorType>
188 void matrix_function_compute_map(const EivalsType& eivals, const ListOfClusters& clusters, VectorType& eivalToCluster)
189 {
190  typedef typename EivalsType::Index Index;
191  eivalToCluster.resize(eivals.rows());
192  Index clusterIndex = 0;
193  for (typename ListOfClusters::const_iterator cluster = clusters.begin(); cluster != clusters.end(); ++cluster) {
194  for (Index i = 0; i < eivals.rows(); ++i) {
195  if (std::find(cluster->begin(), cluster->end(), i) != cluster->end()) {
196  eivalToCluster[i] = clusterIndex;
197  }
198  }
199  ++clusterIndex;
200  }
201 }
202 
204 template <typename DynVectorType, typename VectorType>
205 void matrix_function_compute_permutation(const DynVectorType& blockStart, const DynVectorType& eivalToCluster, VectorType& permutation)
206 {
207  typedef typename VectorType::Index Index;
208  DynVectorType indexNextEntry = blockStart;
209  permutation.resize(eivalToCluster.rows());
210  for (Index i = 0; i < eivalToCluster.rows(); i++) {
211  Index cluster = eivalToCluster[i];
212  permutation[i] = indexNextEntry[cluster];
213  ++indexNextEntry[cluster];
214  }
215 }
216 
218 template <typename VectorType, typename MatrixType>
219 void matrix_function_permute_schur(VectorType& permutation, MatrixType& U, MatrixType& T)
220 {
221  typedef typename VectorType::Index Index;
222  for (Index i = 0; i < permutation.rows() - 1; i++) {
223  Index j;
224  for (j = i; j < permutation.rows(); j++) {
225  if (permutation(j) == i) break;
226  }
227  eigen_assert(permutation(j) == i);
228  for (Index k = j-1; k >= i; k--) {
229  JacobiRotation<typename MatrixType::Scalar> rotation;
230  rotation.makeGivens(T(k, k+1), T(k+1, k+1) - T(k, k));
231  T.applyOnTheLeft(k, k+1, rotation.adjoint());
232  T.applyOnTheRight(k, k+1, rotation);
233  U.applyOnTheRight(k, k+1, rotation);
234  std::swap(permutation.coeffRef(k), permutation.coeffRef(k+1));
235  }
236  }
237 }
238 
245 template <typename MatrixType, typename AtomicType, typename VectorType>
246 void matrix_function_compute_block_atomic(const MatrixType& T, AtomicType& atomic, const VectorType& blockStart, const VectorType& clusterSize, MatrixType& fT)
247 {
248  fT.setZero(T.rows(), T.cols());
249  for (typename VectorType::Index i = 0; i < clusterSize.rows(); ++i) {
250  fT.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i))
251  = atomic.compute(T.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i)));
252  }
253 }
254 
277 template <typename MatrixType>
278 MatrixType matrix_function_solve_triangular_sylvester(const MatrixType& A, const MatrixType& B, const MatrixType& C)
279 {
280  eigen_assert(A.rows() == A.cols());
281  eigen_assert(A.isUpperTriangular());
282  eigen_assert(B.rows() == B.cols());
283  eigen_assert(B.isUpperTriangular());
284  eigen_assert(C.rows() == A.rows());
285  eigen_assert(C.cols() == B.rows());
286 
287  typedef typename MatrixType::Index Index;
288  typedef typename MatrixType::Scalar Scalar;
289 
290  Index m = A.rows();
291  Index n = B.rows();
292  MatrixType X(m, n);
293 
294  for (Index i = m - 1; i >= 0; --i) {
295  for (Index j = 0; j < n; ++j) {
296 
297  // Compute AX = \sum_{k=i+1}^m A_{ik} X_{kj}
298  Scalar AX;
299  if (i == m - 1) {
300  AX = 0;
301  } else {
302  Matrix<Scalar,1,1> AXmatrix = A.row(i).tail(m-1-i) * X.col(j).tail(m-1-i);
303  AX = AXmatrix(0,0);
304  }
305 
306  // Compute XB = \sum_{k=1}^{j-1} X_{ik} B_{kj}
307  Scalar XB;
308  if (j == 0) {
309  XB = 0;
310  } else {
311  Matrix<Scalar,1,1> XBmatrix = X.row(i).head(j) * B.col(j).head(j);
312  XB = XBmatrix(0,0);
313  }
314 
315  X(i,j) = (C(i,j) - AX - XB) / (A(i,i) + B(j,j));
316  }
317  }
318  return X;
319 }
320 
327 template <typename MatrixType, typename VectorType>
328 void matrix_function_compute_above_diagonal(const MatrixType& T, const VectorType& blockStart, const VectorType& clusterSize, MatrixType& fT)
329 {
330  typedef internal::traits<MatrixType> Traits;
331  typedef typename MatrixType::Scalar Scalar;
332  typedef typename MatrixType::Index Index;
333  static const int RowsAtCompileTime = Traits::RowsAtCompileTime;
334  static const int ColsAtCompileTime = Traits::ColsAtCompileTime;
335  static const int Options = MatrixType::Options;
336  typedef Matrix<Scalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType;
337 
338  for (Index k = 1; k < clusterSize.rows(); k++) {
339  for (Index i = 0; i < clusterSize.rows() - k; i++) {
340  // compute (i, i+k) block
341  DynMatrixType A = T.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i));
342  DynMatrixType B = -T.block(blockStart(i+k), blockStart(i+k), clusterSize(i+k), clusterSize(i+k));
343  DynMatrixType C = fT.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i))
344  * T.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k));
345  C -= T.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k))
346  * fT.block(blockStart(i+k), blockStart(i+k), clusterSize(i+k), clusterSize(i+k));
347  for (Index m = i + 1; m < i + k; m++) {
348  C += fT.block(blockStart(i), blockStart(m), clusterSize(i), clusterSize(m))
349  * T.block(blockStart(m), blockStart(i+k), clusterSize(m), clusterSize(i+k));
350  C -= T.block(blockStart(i), blockStart(m), clusterSize(i), clusterSize(m))
351  * fT.block(blockStart(m), blockStart(i+k), clusterSize(m), clusterSize(i+k));
352  }
353  fT.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k))
354  = matrix_function_solve_triangular_sylvester(A, B, C);
355  }
356  }
357 }
358 
374 template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
375 struct matrix_function_compute
376 {
387  template <typename AtomicType, typename ResultType>
388  static void run(const MatrixType& A, AtomicType& atomic, ResultType &result);
389 };
390 
397 template <typename MatrixType>
398 struct matrix_function_compute<MatrixType, 0>
399 {
400  template <typename AtomicType, typename ResultType>
401  static void run(const MatrixType& A, AtomicType& atomic, ResultType &result)
402  {
403  typedef internal::traits<MatrixType> Traits;
404  typedef typename Traits::Scalar Scalar;
405  static const int Rows = Traits::RowsAtCompileTime, Cols = Traits::ColsAtCompileTime;
406  static const int Options = MatrixType::Options;
407  static const int MaxRows = Traits::MaxRowsAtCompileTime, MaxCols = Traits::MaxColsAtCompileTime;
408 
409  typedef std::complex<Scalar> ComplexScalar;
410  typedef Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols> ComplexMatrix;
411 
412  ComplexMatrix CA = A.template cast<ComplexScalar>();
413  ComplexMatrix Cresult;
414  matrix_function_compute<ComplexMatrix>::run(CA, atomic, Cresult);
415  result = Cresult.real();
416  }
417 };
418 
422 template <typename MatrixType>
423 struct matrix_function_compute<MatrixType, 1>
424 {
425  template <typename AtomicType, typename ResultType>
426  static void run(const MatrixType& A, AtomicType& atomic, ResultType &result)
427  {
428  typedef internal::traits<MatrixType> Traits;
429  typedef typename MatrixType::Index Index;
430 
431  // compute Schur decomposition of A
432  const ComplexSchur<MatrixType> schurOfA(A);
433  MatrixType T = schurOfA.matrixT();
434  MatrixType U = schurOfA.matrixU();
435 
436  // partition eigenvalues into clusters of ei'vals "close" to each other
437  std::list<std::list<Index> > clusters;
438  matrix_function_partition_eigenvalues(T.diagonal(), clusters);
439 
440  // compute size of each cluster
441  Matrix<Index, Dynamic, 1> clusterSize;
442  matrix_function_compute_cluster_size(clusters, clusterSize);
443 
444  // blockStart[i] is row index at which block corresponding to i-th cluster starts
445  Matrix<Index, Dynamic, 1> blockStart;
446  matrix_function_compute_block_start(clusterSize, blockStart);
447 
448  // compute map so that eivalToCluster[i] = j means that i-th ei'val is in j-th cluster
449  Matrix<Index, Dynamic, 1> eivalToCluster;
450  matrix_function_compute_map(T.diagonal(), clusters, eivalToCluster);
451 
452  // compute permutation which groups ei'vals in same cluster together
453  Matrix<Index, Traits::RowsAtCompileTime, 1> permutation;
454  matrix_function_compute_permutation(blockStart, eivalToCluster, permutation);
455 
456  // permute Schur decomposition
457  matrix_function_permute_schur(permutation, U, T);
458 
459  // compute result
460  MatrixType fT; // matrix function applied to T
461  matrix_function_compute_block_atomic(T, atomic, blockStart, clusterSize, fT);
462  matrix_function_compute_above_diagonal(T, blockStart, clusterSize, fT);
463  result = U * (fT.template triangularView<Upper>() * U.adjoint());
464  }
465 };
466 
467 } // end of namespace internal
468 
479 template<typename Derived> class MatrixFunctionReturnValue
480 : public ReturnByValue<MatrixFunctionReturnValue<Derived> >
481 {
482  public:
483  typedef typename Derived::Scalar Scalar;
484  typedef typename Derived::Index Index;
485  typedef typename internal::stem_function<Scalar>::type StemFunction;
486 
487  protected:
488  typedef typename internal::ref_selector<Derived>::type DerivedNested;
489 
490  public:
491 
497  MatrixFunctionReturnValue(const Derived& A, StemFunction f) : m_A(A), m_f(f) { }
498 
503  template <typename ResultType>
504  inline void evalTo(ResultType& result) const
505  {
506  typedef typename internal::nested_eval<Derived, 10>::type NestedEvalType;
507  typedef typename internal::remove_all<NestedEvalType>::type NestedEvalTypeClean;
508  typedef internal::traits<NestedEvalTypeClean> Traits;
509  static const int RowsAtCompileTime = Traits::RowsAtCompileTime;
510  static const int ColsAtCompileTime = Traits::ColsAtCompileTime;
511  static const int Options = NestedEvalTypeClean::Options;
512  typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
513  typedef Matrix<ComplexScalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType;
514 
515  typedef internal::MatrixFunctionAtomic<DynMatrixType> AtomicType;
516  AtomicType atomic(m_f);
517 
518  internal::matrix_function_compute<NestedEvalTypeClean>::run(m_A, atomic, result);
519  }
520 
521  Index rows() const { return m_A.rows(); }
522  Index cols() const { return m_A.cols(); }
523 
524  private:
525  const DerivedNested m_A;
526  StemFunction *m_f;
527 };
528 
529 namespace internal {
530 template<typename Derived>
531 struct traits<MatrixFunctionReturnValue<Derived> >
532 {
533  typedef typename Derived::PlainObject ReturnType;
534 };
535 }
536 
537 
538 /********** MatrixBase methods **********/
539 
540 
541 template <typename Derived>
542 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename internal::stem_function<typename internal::traits<Derived>::Scalar>::type f) const
543 {
544  eigen_assert(rows() == cols());
545  return MatrixFunctionReturnValue<Derived>(derived(), f);
546 }
547 
548 template <typename Derived>
549 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const
550 {
551  eigen_assert(rows() == cols());
552  typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
553  return MatrixFunctionReturnValue<Derived>(derived(), internal::stem_function_sin<ComplexScalar>);
554 }
555 
556 template <typename Derived>
557 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const
558 {
559  eigen_assert(rows() == cols());
560  typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
561  return MatrixFunctionReturnValue<Derived>(derived(), internal::stem_function_cos<ComplexScalar>);
562 }
563 
564 template <typename Derived>
565 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const
566 {
567  eigen_assert(rows() == cols());
568  typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
569  return MatrixFunctionReturnValue<Derived>(derived(), internal::stem_function_sinh<ComplexScalar>);
570 }
571 
572 template <typename Derived>
573 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const
574 {
575  eigen_assert(rows() == cols());
576  typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
577  return MatrixFunctionReturnValue<Derived>(derived(), internal::stem_function_cosh<ComplexScalar>);
578 }
579 
580 } // end namespace Eigen
581 
582 #endif // EIGEN_MATRIX_FUNCTION
Namespace containing all symbols from the Eigen library.
Definition: CXX11Meta.h:13
Proxy for the matrix function of some matrix (expression).
Definition: MatrixFunction.h:479
MatrixFunctionReturnValue(const Derived &A, StemFunction f)
Constructor.
Definition: MatrixFunction.h:497
void evalTo(ResultType &result) const
Compute the matrix function.
Definition: MatrixFunction.h:504