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