MatrixFunction.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2009-2011, 2013 Jitse Niesen <jitse@maths.leeds.ac.uk>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_MATRIX_FUNCTION
00011 #define EIGEN_MATRIX_FUNCTION
00012 
00013 #include "StemFunction.h"
00014 
00015 
00016 namespace Eigen { 
00017 
00018 namespace internal {
00019 
00021 static const float matrix_function_separation = 0.1f;
00022 
00029 template <typename MatrixType>
00030 class MatrixFunctionAtomic 
00031 {
00032   public:
00033 
00034     typedef typename MatrixType::Scalar Scalar;
00035     typedef typename stem_function<Scalar>::type StemFunction;
00036 
00040     MatrixFunctionAtomic(StemFunction f) : m_f(f) { }
00041 
00046     MatrixType compute(const MatrixType& A);
00047 
00048   private:
00049     StemFunction* m_f;
00050 };
00051 
00052 template <typename MatrixType>
00053 typename NumTraits<typename MatrixType::Scalar>::Real matrix_function_compute_mu(const MatrixType& A)
00054 {
00055   typedef typename plain_col_type<MatrixType>::type VectorType;
00056   typename MatrixType::Index rows = A.rows();
00057   const MatrixType N = MatrixType::Identity(rows, rows) - A;
00058   VectorType e = VectorType::Ones(rows);
00059   N.template triangularView<Upper>().solveInPlace(e);
00060   return e.cwiseAbs().maxCoeff();
00061 }
00062 
00063 template <typename MatrixType>
00064 MatrixType MatrixFunctionAtomic<MatrixType>::compute(const MatrixType& A)
00065 {
00066   // TODO: Use that A is upper triangular
00067   typedef typename NumTraits<Scalar>::Real RealScalar;
00068   typedef typename MatrixType::Index Index;
00069   Index rows = A.rows();
00070   Scalar avgEival = A.trace() / Scalar(RealScalar(rows));
00071   MatrixType Ashifted = A - avgEival * MatrixType::Identity(rows, rows);
00072   RealScalar mu = matrix_function_compute_mu(Ashifted);
00073   MatrixType F = m_f(avgEival, 0) * MatrixType::Identity(rows, rows);
00074   MatrixType P = Ashifted;
00075   MatrixType Fincr;
00076   for (Index s = 1; s < 1.1 * rows + 10; s++) { // upper limit is fairly arbitrary
00077     Fincr = m_f(avgEival, static_cast<int>(s)) * P;
00078     F += Fincr;
00079     P = Scalar(RealScalar(1.0/(s + 1))) * P * Ashifted;
00080 
00081     // test whether Taylor series converged
00082     const RealScalar F_norm = F.cwiseAbs().rowwise().sum().maxCoeff();
00083     const RealScalar Fincr_norm = Fincr.cwiseAbs().rowwise().sum().maxCoeff();
00084     if (Fincr_norm < NumTraits<Scalar>::epsilon() * F_norm) {
00085       RealScalar delta = 0;
00086       RealScalar rfactorial = 1;
00087       for (Index r = 0; r < rows; r++) {
00088         RealScalar mx = 0;
00089         for (Index i = 0; i < rows; i++)
00090           mx = (std::max)(mx, std::abs(m_f(Ashifted(i, i) + avgEival, static_cast<int>(s+r))));
00091         if (r != 0)
00092           rfactorial *= RealScalar(r);
00093         delta = (std::max)(delta, mx / rfactorial);
00094       }
00095       const RealScalar P_norm = P.cwiseAbs().rowwise().sum().maxCoeff();
00096       if (mu * delta * P_norm < NumTraits<Scalar>::epsilon() * F_norm) // series converged
00097         break;
00098     }
00099   }
00100   return F;
00101 }
00102 
00108 template <typename Index, typename ListOfClusters>
00109 typename ListOfClusters::iterator matrix_function_find_cluster(Index key, ListOfClusters& clusters)
00110 {
00111   typename std::list<Index>::iterator j;
00112   for (typename ListOfClusters::iterator i = clusters.begin(); i != clusters.end(); ++i) {
00113     j = std::find(i->begin(), i->end(), key);
00114     if (j != i->end())
00115       return i;
00116   }
00117   return clusters.end();
00118 }
00119 
00131 template <typename EivalsType, typename Cluster>
00132 void matrix_function_partition_eigenvalues(const EivalsType& eivals, std::list<Cluster>& clusters)
00133 {
00134   typedef typename EivalsType::Index Index;
00135   typedef typename EivalsType::RealScalar RealScalar;
00136   for (Index i=0; i<eivals.rows(); ++i) {
00137     // Find cluster containing i-th ei'val, adding a new cluster if necessary
00138     typename std::list<Cluster>::iterator qi = matrix_function_find_cluster(i, clusters);
00139     if (qi == clusters.end()) {
00140       Cluster l;
00141       l.push_back(i);
00142       clusters.push_back(l);
00143       qi = clusters.end();
00144       --qi;
00145     }
00146 
00147     // Look for other element to add to the set
00148     for (Index j=i+1; j<eivals.rows(); ++j) {
00149       if (abs(eivals(j) - eivals(i)) <= RealScalar(matrix_function_separation)
00150           && std::find(qi->begin(), qi->end(), j) == qi->end()) {
00151         typename std::list<Cluster>::iterator qj = matrix_function_find_cluster(j, clusters);
00152         if (qj == clusters.end()) {
00153           qi->push_back(j);
00154         } else {
00155           qi->insert(qi->end(), qj->begin(), qj->end());
00156           clusters.erase(qj);
00157         }
00158       }
00159     }
00160   }
00161 }
00162 
00164 template <typename ListOfClusters, typename Index>
00165 void matrix_function_compute_cluster_size(const ListOfClusters& clusters, Matrix<Index, Dynamic, 1>& clusterSize)
00166 {
00167   const Index numClusters = static_cast<Index>(clusters.size());
00168   clusterSize.setZero(numClusters);
00169   Index clusterIndex = 0;
00170   for (typename ListOfClusters::const_iterator cluster = clusters.begin(); cluster != clusters.end(); ++cluster) {
00171     clusterSize[clusterIndex] = cluster->size();
00172     ++clusterIndex;
00173   }
00174 }
00175 
00177 template <typename VectorType>
00178 void matrix_function_compute_block_start(const VectorType& clusterSize, VectorType& blockStart)
00179 {
00180   blockStart.resize(clusterSize.rows());
00181   blockStart(0) = 0;
00182   for (typename VectorType::Index i = 1; i < clusterSize.rows(); i++) {
00183     blockStart(i) = blockStart(i-1) + clusterSize(i-1);
00184   }
00185 }
00186 
00188 template <typename EivalsType, typename ListOfClusters, typename VectorType>
00189 void matrix_function_compute_map(const EivalsType& eivals, const ListOfClusters& clusters, VectorType& eivalToCluster)
00190 {
00191   typedef typename EivalsType::Index Index;
00192   eivalToCluster.resize(eivals.rows());
00193   Index clusterIndex = 0;
00194   for (typename ListOfClusters::const_iterator cluster = clusters.begin(); cluster != clusters.end(); ++cluster) {
00195     for (Index i = 0; i < eivals.rows(); ++i) {
00196       if (std::find(cluster->begin(), cluster->end(), i) != cluster->end()) {
00197         eivalToCluster[i] = clusterIndex;
00198       }
00199     }
00200     ++clusterIndex;
00201   }
00202 }
00203 
00205 template <typename DynVectorType, typename VectorType>
00206 void matrix_function_compute_permutation(const DynVectorType& blockStart, const DynVectorType& eivalToCluster, VectorType& permutation)
00207 {
00208   typedef typename VectorType::Index Index;
00209   DynVectorType indexNextEntry = blockStart;
00210   permutation.resize(eivalToCluster.rows());
00211   for (Index i = 0; i < eivalToCluster.rows(); i++) {
00212     Index cluster = eivalToCluster[i];
00213     permutation[i] = indexNextEntry[cluster];
00214     ++indexNextEntry[cluster];
00215   }
00216 }  
00217 
00219 template <typename VectorType, typename MatrixType>
00220 void matrix_function_permute_schur(VectorType& permutation, MatrixType& U, MatrixType& T)
00221 {
00222   typedef typename VectorType::Index Index;
00223   for (Index i = 0; i < permutation.rows() - 1; i++) {
00224     Index j;
00225     for (j = i; j < permutation.rows(); j++) {
00226       if (permutation(j) == i) break;
00227     }
00228     eigen_assert(permutation(j) == i);
00229     for (Index k = j-1; k >= i; k--) {
00230       JacobiRotation<typename MatrixType::Scalar> rotation;
00231       rotation.makeGivens(T(k, k+1), T(k+1, k+1) - T(k, k));
00232       T.applyOnTheLeft(k, k+1, rotation.adjoint());
00233       T.applyOnTheRight(k, k+1, rotation);
00234       U.applyOnTheRight(k, k+1, rotation);
00235       std::swap(permutation.coeffRef(k), permutation.coeffRef(k+1));
00236     }
00237   }
00238 }
00239 
00246 template <typename MatrixType, typename AtomicType, typename VectorType>
00247 void matrix_function_compute_block_atomic(const MatrixType& T, AtomicType& atomic, const VectorType& blockStart, const VectorType& clusterSize, MatrixType& fT)
00248 { 
00249   fT.setZero(T.rows(), T.cols());
00250   for (typename VectorType::Index i = 0; i < clusterSize.rows(); ++i) {
00251     fT.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i))
00252       = atomic.compute(T.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i)));
00253   }
00254 }
00255 
00278 template <typename MatrixType>
00279 MatrixType matrix_function_solve_triangular_sylvester(const MatrixType& A, const MatrixType& B, const MatrixType& C)
00280 {
00281   eigen_assert(A.rows() == A.cols());
00282   eigen_assert(A.isUpperTriangular());
00283   eigen_assert(B.rows() == B.cols());
00284   eigen_assert(B.isUpperTriangular());
00285   eigen_assert(C.rows() == A.rows());
00286   eigen_assert(C.cols() == B.rows());
00287 
00288   typedef typename MatrixType::Index Index;
00289   typedef typename MatrixType::Scalar Scalar;
00290 
00291   Index m = A.rows();
00292   Index n = B.rows();
00293   MatrixType X(m, n);
00294 
00295   for (Index i = m - 1; i >= 0; --i) {
00296     for (Index j = 0; j < n; ++j) {
00297 
00298       // Compute AX = \sum_{k=i+1}^m A_{ik} X_{kj}
00299       Scalar AX;
00300       if (i == m - 1) {
00301         AX = 0; 
00302       } else {
00303         Matrix<Scalar,1,1> AXmatrix = A.row(i).tail(m-1-i) * X.col(j).tail(m-1-i);
00304         AX = AXmatrix(0,0);
00305       }
00306 
00307       // Compute XB = \sum_{k=1}^{j-1} X_{ik} B_{kj}
00308       Scalar XB;
00309       if (j == 0) {
00310         XB = 0; 
00311       } else {
00312         Matrix<Scalar,1,1> XBmatrix = X.row(i).head(j) * B.col(j).head(j);
00313         XB = XBmatrix(0,0);
00314       }
00315 
00316       X(i,j) = (C(i,j) - AX - XB) / (A(i,i) + B(j,j));
00317     }
00318   }
00319   return X;
00320 }
00321 
00328 template <typename MatrixType, typename VectorType>
00329 void matrix_function_compute_above_diagonal(const MatrixType& T, const VectorType& blockStart, const VectorType& clusterSize, MatrixType& fT)
00330 { 
00331   typedef internal::traits<MatrixType> Traits;
00332   typedef typename MatrixType::Scalar Scalar;
00333   typedef typename MatrixType::Index Index;
00334   static const int RowsAtCompileTime = Traits::RowsAtCompileTime;
00335   static const int ColsAtCompileTime = Traits::ColsAtCompileTime;
00336   static const int Options = MatrixType::Options;
00337   typedef Matrix<Scalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType;
00338 
00339   for (Index k = 1; k < clusterSize.rows(); k++) {
00340     for (Index i = 0; i < clusterSize.rows() - k; i++) {
00341       // compute (i, i+k) block
00342       DynMatrixType A = T.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i));
00343       DynMatrixType B = -T.block(blockStart(i+k), blockStart(i+k), clusterSize(i+k), clusterSize(i+k));
00344       DynMatrixType C = fT.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i))
00345         * T.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k));
00346       C -= T.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k))
00347         * fT.block(blockStart(i+k), blockStart(i+k), clusterSize(i+k), clusterSize(i+k));
00348       for (Index m = i + 1; m < i + k; m++) {
00349         C += fT.block(blockStart(i), blockStart(m), clusterSize(i), clusterSize(m))
00350           * T.block(blockStart(m), blockStart(i+k), clusterSize(m), clusterSize(i+k));
00351         C -= T.block(blockStart(i), blockStart(m), clusterSize(i), clusterSize(m))
00352           * fT.block(blockStart(m), blockStart(i+k), clusterSize(m), clusterSize(i+k));
00353       }
00354       fT.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k))
00355         = matrix_function_solve_triangular_sylvester(A, B, C);
00356     }
00357   }
00358 }
00359 
00375 template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
00376 struct matrix_function_compute
00377 {  
00388     template <typename AtomicType, typename ResultType> 
00389     static void run(const MatrixType& A, AtomicType& atomic, ResultType &result);    
00390 };
00391 
00398 template <typename MatrixType>
00399 struct matrix_function_compute<MatrixType, 0>
00400 {  
00401   template <typename AtomicType, typename ResultType> 
00402   static void run(const MatrixType& A, AtomicType& atomic, ResultType &result)
00403   {
00404     typedef internal::traits<MatrixType> Traits;
00405     typedef typename Traits::Scalar Scalar;
00406     static const int Rows = Traits::RowsAtCompileTime, Cols = Traits::ColsAtCompileTime;
00407     static const int MaxRows = Traits::MaxRowsAtCompileTime, MaxCols = Traits::MaxColsAtCompileTime;
00408 
00409     typedef std::complex<Scalar> ComplexScalar;
00410     typedef Matrix<ComplexScalar, Rows, Cols, 0, MaxRows, MaxCols> ComplexMatrix;
00411 
00412     ComplexMatrix CA = A.template cast<ComplexScalar>();
00413     ComplexMatrix Cresult;
00414     matrix_function_compute<ComplexMatrix>::run(CA, atomic, Cresult);
00415     result = Cresult.real();
00416   }
00417 };
00418 
00422 template <typename MatrixType>
00423 struct matrix_function_compute<MatrixType, 1>
00424 {
00425   template <typename AtomicType, typename ResultType> 
00426   static void run(const MatrixType& A, AtomicType& atomic, ResultType &result)
00427   {
00428     typedef internal::traits<MatrixType> Traits;
00429     typedef typename MatrixType::Index Index;
00430     
00431     // compute Schur decomposition of A
00432     const ComplexSchur<MatrixType> schurOfA(A);  
00433     MatrixType T = schurOfA.matrixT();
00434     MatrixType U = schurOfA.matrixU();
00435 
00436     // partition eigenvalues into clusters of ei'vals "close" to each other
00437     std::list<std::list<Index> > clusters; 
00438     matrix_function_partition_eigenvalues(T.diagonal(), clusters);
00439 
00440     // compute size of each cluster
00441     Matrix<Index, Dynamic, 1> clusterSize;
00442     matrix_function_compute_cluster_size(clusters, clusterSize);
00443 
00444     // blockStart[i] is row index at which block corresponding to i-th cluster starts 
00445     Matrix<Index, Dynamic, 1> blockStart; 
00446     matrix_function_compute_block_start(clusterSize, blockStart);
00447 
00448     // compute map so that eivalToCluster[i] = j means that i-th ei'val is in j-th cluster 
00449     Matrix<Index, Dynamic, 1> eivalToCluster;
00450     matrix_function_compute_map(T.diagonal(), clusters, eivalToCluster);
00451 
00452     // compute permutation which groups ei'vals in same cluster together 
00453     Matrix<Index, Traits::RowsAtCompileTime, 1> permutation;
00454     matrix_function_compute_permutation(blockStart, eivalToCluster, permutation);
00455 
00456     // permute Schur decomposition
00457     matrix_function_permute_schur(permutation, U, T);
00458 
00459     // compute result
00460     MatrixType fT; // matrix function applied to T
00461     matrix_function_compute_block_atomic(T, atomic, blockStart, clusterSize, fT);
00462     matrix_function_compute_above_diagonal(T, blockStart, clusterSize, fT);
00463     result = U * (fT.template triangularView<Upper>() * U.adjoint());
00464   }
00465 };
00466 
00467 } // end of namespace internal
00468 
00479 template<typename Derived> class MatrixFunctionReturnValue
00480 : public ReturnByValue<MatrixFunctionReturnValue<Derived> >
00481 {
00482   public:
00483     typedef typename Derived::Scalar Scalar;
00484     typedef typename Derived::Index Index;
00485     typedef typename internal::stem_function<Scalar>::type StemFunction;
00486 
00487   protected:
00488     typedef typename internal::ref_selector<Derived>::type DerivedNested;
00489 
00490   public:
00491 
00497     MatrixFunctionReturnValue(const Derived& A, StemFunction f) : m_A(A), m_f(f) { }
00498 
00503     template <typename ResultType>
00504     inline void evalTo(ResultType& result) const
00505     {
00506       typedef typename internal::nested_eval<Derived, 10>::type NestedEvalType;
00507       typedef typename internal::remove_all<NestedEvalType>::type NestedEvalTypeClean;
00508       typedef internal::traits<NestedEvalTypeClean> Traits;
00509       static const int RowsAtCompileTime = Traits::RowsAtCompileTime;
00510       static const int ColsAtCompileTime = Traits::ColsAtCompileTime;
00511       typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
00512       typedef Matrix<ComplexScalar, Dynamic, Dynamic, 0, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType;
00513 
00514       typedef internal::MatrixFunctionAtomic<DynMatrixType> AtomicType;
00515       AtomicType atomic(m_f);
00516 
00517       internal::matrix_function_compute<NestedEvalTypeClean>::run(m_A, atomic, result);
00518     }
00519 
00520     Index rows() const { return m_A.rows(); }
00521     Index cols() const { return m_A.cols(); }
00522 
00523   private:
00524     const DerivedNested m_A;
00525     StemFunction *m_f;
00526 };
00527 
00528 namespace internal {
00529 template<typename Derived>
00530 struct traits<MatrixFunctionReturnValue<Derived> >
00531 {
00532   typedef typename Derived::PlainObject ReturnType;
00533 };
00534 }
00535 
00536 
00537 /********** MatrixBase methods **********/
00538 
00539 
00540 template <typename Derived>
00541 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename internal::stem_function<typename internal::traits<Derived>::Scalar>::type f) const
00542 {
00543   eigen_assert(rows() == cols());
00544   return MatrixFunctionReturnValue<Derived>(derived(), f);
00545 }
00546 
00547 template <typename Derived>
00548 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const
00549 {
00550   eigen_assert(rows() == cols());
00551   typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
00552   return MatrixFunctionReturnValue<Derived>(derived(), internal::stem_function_sin<ComplexScalar>);
00553 }
00554 
00555 template <typename Derived>
00556 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const
00557 {
00558   eigen_assert(rows() == cols());
00559   typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
00560   return MatrixFunctionReturnValue<Derived>(derived(), internal::stem_function_cos<ComplexScalar>);
00561 }
00562 
00563 template <typename Derived>
00564 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const
00565 {
00566   eigen_assert(rows() == cols());
00567   typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
00568   return MatrixFunctionReturnValue<Derived>(derived(), internal::stem_function_sinh<ComplexScalar>);
00569 }
00570 
00571 template <typename Derived>
00572 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const
00573 {
00574   eigen_assert(rows() == cols());
00575   typedef typename internal::stem_function<Scalar>::ComplexScalar ComplexScalar;
00576   return MatrixFunctionReturnValue<Derived>(derived(), internal::stem_function_cosh<ComplexScalar>);
00577 }
00578 
00579 } // end namespace Eigen
00580 
00581 #endif // EIGEN_MATRIX_FUNCTION
 All Classes Functions Variables Typedefs Enumerator