Implement matrix logarithm + test + docs.

Currently, test matrix_function_1 fails due to bug #288.

6 files changed, 434 insertions, 0 deletions

diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index bc5766f04..e2d34d257 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -466,6 +466,7 @@ template<typename Derived> class MatrixBase
     const MatrixFunctionReturnValue<Derived> cos() const;
     const MatrixFunctionReturnValue<Derived> sin() const;
     const MatrixSquareRootReturnValue<Derived> sqrt() const;
+    const MatrixLogarithmReturnValue<Derived> log() const;
 
 #ifdef EIGEN2_SUPPORT
     template<typename ProductDerived, typename Lhs, typename Rhs>
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index 74be25216..9f5fbc9f6 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -284,6 +284,7 @@ template<typename MatrixType,int Direction> class Homogeneous;
 template<typename Derived> struct MatrixExponentialReturnValue;
 template<typename Derived> class MatrixFunctionReturnValue;
 template<typename Derived> class MatrixSquareRootReturnValue;
+template<typename Derived> class MatrixLogarithmReturnValue;
 
 namespace internal {
 template <typename Scalar>
diff --git a/unsupported/Eigen/MatrixFunctions b/unsupported/Eigen/MatrixFunctions
index ecdf23cef..43955d352 100644
--- a/unsupported/Eigen/MatrixFunctions
+++ b/unsupported/Eigen/MatrixFunctions
@@ -50,6 +50,7 @@ namespace Eigen {
   *  - \ref matrixbase_cos "MatrixBase::cos()", for computing the matrix cosine
   *  - \ref matrixbase_cosh "MatrixBase::cosh()", for computing the matrix hyperbolic cosine
   *  - \ref matrixbase_exp "MatrixBase::exp()", for computing the matrix exponential
+  *  - \ref matrixbase_log "MatrixBase::log()", for computing the matrix logarithm
   *  - \ref matrixbase_matrixfunction "MatrixBase::matrixFunction()", for computing general matrix functions
   *  - \ref matrixbase_sin "MatrixBase::sin()", for computing the matrix sine
   *  - \ref matrixbase_sinh "MatrixBase::sinh()", for computing the matrix hyperbolic sine
@@ -71,6 +72,7 @@ namespace Eigen {
 #include "src/MatrixFunctions/MatrixExponential.h"
 #include "src/MatrixFunctions/MatrixFunction.h"
 #include "src/MatrixFunctions/MatrixSquareRoot.h"
+#include "src/MatrixFunctions/MatrixLogarithm.h"
 
 
 
@@ -172,6 +174,60 @@ Output: \verbinclude MatrixExponential.out
 \c complex<float> or \c complex<double> .
 
 
+\section matrixbase_log MatrixBase::log()
+
+Compute the matrix logarithm.
+
+\code
+const MatrixLogarithmReturnValue<Derived> MatrixBase<Derived>::log() const
+\endcode
+
+\param[in]  M  invertible matrix whose logarithm is to be computed.
+\returns    expression representing the matrix logarithm root of \p M.
+
+The matrix logarithm of \f$ M \f$ is a matrix \f$ X \f$ such that 
+\f$ \exp(X) = M \f$ where exp denotes the matrix exponential. As for
+the scalar logarithm, the equation \f$ \exp(X) = M \f$ may have
+multiple solutions; this function returns a matrix whose eigenvalues
+have imaginary part in the interval \f$ (-\pi,\pi] \f$.
+
+In the real case, the matrix \f$ M \f$ should be invertible and
+it should have no eigenvalues which are real and negative (pairs of
+complex conjugate eigenvalues are allowed). In the complex case, it
+only needs to be invertible.
+
+This function computes the matrix logarithm using the Schur-Parlett
+algorithm as implemented by MatrixBase::matrixFunction(). The
+logarithm of an atomic block is computed by MatrixLogarithmAtomic,
+which uses direct computation for 1-by-1 and 2-by-2 blocks and an
+inverse scaling-and-squaring algorithm for bigger blocks, with the
+square roots computed by MatrixBase::sqrt().
+
+Details of the algorithm can be found in Section 11.6.2 of:
+Nicholas J. Higham,
+<em>Functions of Matrices: Theory and Computation</em>,
+SIAM 2008. ISBN 978-0-898716-46-7.
+
+Example: The following program checks that
+\f[ \log \left[ \begin{array}{ccc} 
+      \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
+      \frac12\sqrt2 & \frac12\sqrt2 & 0 \\
+      0 & 0 & 1
+    \end{array} \right] = \left[ \begin{array}{ccc}
+      0 & \frac14\pi & 0 \\ 
+      -\frac14\pi & 0 & 0 \\
+      0 & 0 & 0 
+    \end{array} \right]. \f]
+This corresponds to a rotation of \f$ \frac14\pi \f$ radians around
+the z-axis. This is the inverse of the example used in the
+documentation of \ref matrixbase_exp "exp()".
+
+\include MatrixLogarithm.cpp
+Output: \verbinclude MatrixLogarithm.out
+
+\sa MatrixBase::exp(), MatrixBase::matrixFunction(), 
+    class MatrixLogarithmAtomic, MatrixBase::sqrt().
+
 
 \section matrixbase_matrixfunction MatrixBase::matrixFunction()
 
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h
new file mode 100644
index 000000000..8cdcca659
--- /dev/null
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h
@@ -0,0 +1,340 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2011 Jitse Niesen <jitse@maths.leeds.ac.uk>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_MATRIX_LOGARITHM
+#define EIGEN_MATRIX_LOGARITHM
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+/** \ingroup MatrixFunctions_Module
+  * \class MatrixLogarithmAtomic
+  * \brief Helper class for computing matrix logarithm of atomic matrices.
+  *
+  * \internal
+  * Here, an atomic matrix is a triangular matrix whose diagonal
+  * entries are close to each other.
+  *
+  * \sa class MatrixFunctionAtomic, MatrixBase::log()
+  */
+template <typename MatrixType>
+class MatrixLogarithmAtomic
+{
+public:
+
+  typedef typename MatrixType::Scalar Scalar;
+  // typedef typename MatrixType::Index Index;
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+  // typedef typename internal::stem_function<Scalar>::type StemFunction;
+  // typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
+
+  /** \brief Constructor. */
+  MatrixLogarithmAtomic() { }
+
+  /** \brief Compute matrix logarithm of atomic matrix
+    * \param[in]  A  argument of matrix logarithm, should be upper triangular and atomic
+    * \returns  The logarithm of \p A.
+    */
+  MatrixType compute(const MatrixType& A);
+
+private:
+
+  void compute2x2(const MatrixType& A, MatrixType& result);
+  void computeBig(const MatrixType& A, MatrixType& result);
+  static Scalar atanh(Scalar x);
+  int getPadeDegree(typename MatrixType::RealScalar normTminusI);
+  void computePade(MatrixType& result, const MatrixType& T, int degree);
+  void computePade3(MatrixType& result, const MatrixType& T);
+  void computePade4(MatrixType& result, const MatrixType& T);
+  void computePade5(MatrixType& result, const MatrixType& T);
+  void computePade6(MatrixType& result, const MatrixType& T);
+  void computePade7(MatrixType& result, const MatrixType& T);
+
+  static const double maxNormForPade[];
+  static const int minPadeDegree = 3;
+  static const int maxPadeDegree = 7;
+
+  // Prevent copying
+  MatrixLogarithmAtomic(const MatrixLogarithmAtomic&);
+  MatrixLogarithmAtomic& operator=(const MatrixLogarithmAtomic&);
+};
+
+template <typename MatrixType>
+const double MatrixLogarithmAtomic<MatrixType>::maxNormForPade[] = { 0.0162 /* degree = 3 */, 0.0539, 0.114, 0.187, 0.264 };
+
+/** \brief Compute logarithm of triangular matrix with clustered eigenvalues. */
+template <typename MatrixType>
+MatrixType MatrixLogarithmAtomic<MatrixType>::compute(const MatrixType& A)
+{
+  using std::log;
+  MatrixType result(A.rows(), A.rows());
+  if (A.rows() == 1)
+    result(0,0) = log(A(0,0));
+  else if (A.rows() == 2)
+    compute2x2(A, result);
+  else
+    computeBig(A, result);
+  return result;
+}
+
+/** \brief Compute atanh (inverse hyperbolic tangent). */
+template <typename MatrixType>
+typename MatrixType::Scalar MatrixLogarithmAtomic<MatrixType>::atanh(typename MatrixType::Scalar x)
+{
+  using std::abs;
+  using std::sqrt;
+  if (abs(x) > sqrt(NumTraits<Scalar>::epsilon()))
+    return Scalar(0.5) * log((Scalar(1) + x) / (Scalar(1) - x));
+  else
+    return x + x*x*x / Scalar(3);
+}
+
+/** \brief Compute logarithm of 2x2 triangular matrix. */
+template <typename MatrixType>
+void MatrixLogarithmAtomic<MatrixType>::compute2x2(const MatrixType& A, MatrixType& result)
+{
+  using std::abs;
+  using std::ceil;
+  using std::imag;
+  using std::log;
+
+  Scalar logA00 = log(A(0,0));
+  Scalar logA11 = log(A(1,1));
+
+  result(0,0) = logA00;
+  result(1,0) = Scalar(0);
+  result(1,1) = logA11;
+
+  if (A(0,0) == A(1,1)) {
+    result(0,1) = A(0,1) / A(0,0);
+  } else if ((abs(A(0,0)) < 0.5*abs(A(1,1))) || (abs(A(0,0)) > 2*abs(A(1,1)))) {
+    result(0,1) = A(0,1) * (logA11 - logA00) / (A(1,1) - A(0,0));
+  } else {
+    // computation in previous branch is inaccurate if A(1,1) \approx A(0,0)
+    int unwindingNumber = ceil((imag(logA11 - logA00) - M_PI) / (2*M_PI));
+    Scalar z = (A(1,1) - A(0,0)) / (A(1,1) + A(0,0));
+    result(0,1) = A(0,1) * (Scalar(2) * atanh(z) + Scalar(0,2*M_PI*unwindingNumber)) / (A(1,1) - A(0,0));
+  }
+}
+
+/** \brief Compute logarithm of triangular matrices with size > 2. 
+  * \details This uses a inverse scale-and-square algorithm. */
+template <typename MatrixType>
+void MatrixLogarithmAtomic<MatrixType>::computeBig(const MatrixType& A, MatrixType& result)
+{
+  int numberOfSquareRoots = 0;
+  int numberOfExtraSquareRoots = 0;
+  int degree;
+  MatrixType T = A;
+
+  while (true) {
+    RealScalar normTminusI = (T - MatrixType::Identity(T.rows(), T.rows())).cwiseAbs().colwise().sum().maxCoeff();
+    if (normTminusI < maxNormForPade[maxPadeDegree - minPadeDegree]) {
+      degree = getPadeDegree(normTminusI);
+      int degree2 = getPadeDegree(normTminusI / RealScalar(2));
+      if ((degree - degree2 <= 1) || (numberOfExtraSquareRoots == 1)) 
+	break;
+      ++numberOfExtraSquareRoots;
+    }
+    T = T.sqrt();
+    ++numberOfSquareRoots;
+  }
+
+  computePade(result, T, degree);
+  result *= pow(RealScalar(2), numberOfSquareRoots);
+}
+
+/* \brief Get suitable degree for Pade approximation. */
+template <typename MatrixType>
+int MatrixLogarithmAtomic<MatrixType>::getPadeDegree(typename MatrixType::RealScalar normTminusI)
+{
+  for (int degree = 3; degree <= maxPadeDegree; ++degree) 
+    if (normTminusI <= maxNormForPade[degree - minPadeDegree])
+      return degree;
+  assert(false); // this line should never be reached
+}
+
+/* \brief Compute Pade approximation to matrix logarithm */
+template <typename MatrixType>
+void MatrixLogarithmAtomic<MatrixType>::computePade(MatrixType& result, const MatrixType& T, int degree)
+{
+  switch (degree) {
+    case 3:  computePade3(result, T); break;
+    case 4:  computePade4(result, T); break;
+    case 5:  computePade5(result, T); break;
+    case 6:  computePade6(result, T); break;
+    case 7:  computePade7(result, T); break;
+    default: assert(false); // should never happen
+  }
+} 
+
+template <typename MatrixType>
+void MatrixLogarithmAtomic<MatrixType>::computePade3(MatrixType& result, const MatrixType& T)
+{
+  const int degree = 3;
+  double nodes[]   = { 0.112701665379258, 0.500000000000000, 0.887298334620742 };
+  double weights[] = { 0.277777777777778, 0.444444444444444, 0.277777777777778 };
+  MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows());
+  result.setZero(T.rows(), T.rows());
+  for (int k = 0; k < degree; ++k) 
+    result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI)
+                           .template triangularView<Upper>().solve(TminusI);
+}
+
+template <typename MatrixType>
+void MatrixLogarithmAtomic<MatrixType>::computePade4(MatrixType& result, const MatrixType& T)
+{
+  const int degree = 4;
+  double nodes[]   = { 0.069431844202974, 0.330009478207572, 0.669990521792428, 0.930568155797026 };
+  double weights[] = { 0.173927422568727, 0.326072577431273, 0.326072577431273, 0.173927422568727 };
+  MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows());
+  result.setZero(T.rows(), T.rows());
+  for (int k = 0; k < degree; ++k) 
+    result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI)
+                           .template triangularView<Upper>().solve(TminusI);
+}
+
+template <typename MatrixType>
+void MatrixLogarithmAtomic<MatrixType>::computePade5(MatrixType& result, const MatrixType& T)
+{
+  const int degree = 5;
+  double nodes[]   = { 0.046910077030668, 0.230765344947158, 0.500000000000000,
+		       0.769234655052841, 0.953089922969332 };
+  double weights[] = { 0.118463442528095, 0.239314335249683, 0.284444444444444,
+		       0.239314335249683, 0.118463442528094 };
+  MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows());
+  result.setZero(T.rows(), T.rows());
+  for (int k = 0; k < degree; ++k) 
+    result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI)
+                           .template triangularView<Upper>().solve(TminusI);
+}
+
+template <typename MatrixType>
+void MatrixLogarithmAtomic<MatrixType>::computePade6(MatrixType& result, const MatrixType& T)
+{
+  const int degree = 6;
+  double nodes[]   = { 0.033765242898424, 0.169395306766868, 0.380690406958402,
+		       0.619309593041598, 0.830604693233132, 0.966234757101576 };
+  double weights[] = { 0.085662246189585, 0.180380786524069, 0.233956967286345,
+		       0.233956967286346, 0.180380786524069, 0.085662246189585 };
+  MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows());
+  result.setZero(T.rows(), T.rows());
+  for (int k = 0; k < degree; ++k) 
+    result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI)
+                           .template triangularView<Upper>().solve(TminusI);
+}
+
+template <typename MatrixType>
+void MatrixLogarithmAtomic<MatrixType>::computePade7(MatrixType& result, const MatrixType& T)
+{
+  const int degree = 7;
+  double nodes[]   = { 0.025446043828621, 0.129234407200303, 0.297077424311301, 0.500000000000000,
+		       0.702922575688699, 0.870765592799697, 0.974553956171379 };
+  double weights[] = { 0.064742483084435, 0.139852695744638, 0.190915025252559, 0.208979591836734,
+		       0.190915025252560, 0.139852695744638, 0.064742483084435 };
+  MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows());
+  result.setZero(T.rows(), T.rows());
+  for (int k = 0; k < degree; ++k) 
+    result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI)
+                           .template triangularView<Upper>().solve(TminusI);
+}
+
+/** \ingroup MatrixFunctions_Module
+  *
+  * \brief Proxy for the matrix logarithm of some matrix (expression).
+  *
+  * \tparam Derived  Type of the argument to the matrix function.
+  *
+  * This class holds the argument to the matrix function until it is
+  * assigned or evaluated for some other reason (so the argument
+  * should not be changed in the meantime). It is the return type of
+  * matrixBase::matrixLogarithm() and most of the time this is the
+  * only way it is used.
+  */
+template<typename Derived> class MatrixLogarithmReturnValue
+: public ReturnByValue<MatrixLogarithmReturnValue<Derived> >
+{
+public:
+
+  typedef typename Derived::Scalar Scalar;
+  typedef typename Derived::Index Index;
+
+  /** \brief Constructor.
+    *
+    * \param[in]  A  %Matrix (expression) forming the argument of the matrix logarithm.
+    */
+  MatrixLogarithmReturnValue(const Derived& A) : m_A(A) { }
+  
+  /** \brief Compute the matrix logarithm.
+    *
+    * \param[out]  result  Logarithm of \p A, where \A is as specified in the constructor.
+    */
+  template <typename ResultType>
+  inline void evalTo(ResultType& result) const
+  {
+    typedef typename Derived::PlainObject PlainObject;
+    typedef internal::traits<PlainObject> Traits;
+    static const int RowsAtCompileTime = Traits::RowsAtCompileTime;
+    static const int ColsAtCompileTime = Traits::ColsAtCompileTime;
+    static const int Options = PlainObject::Options;
+    typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
+    typedef Matrix<ComplexScalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType;
+    typedef MatrixLogarithmAtomic<DynMatrixType> AtomicType;
+    AtomicType atomic;
+    
+    const PlainObject Aevaluated = m_A.eval();
+    MatrixFunction<PlainObject, AtomicType> mf(Aevaluated, atomic);
+    mf.compute(result);
+  }
+
+  Index rows() const { return m_A.rows(); }
+  Index cols() const { return m_A.cols(); }
+  
+private:
+  typename internal::nested<Derived>::type m_A;
+  
+  MatrixLogarithmReturnValue& operator=(const MatrixLogarithmReturnValue&);
+};
+
+namespace internal {
+  template<typename Derived>
+  struct traits<MatrixLogarithmReturnValue<Derived> >
+  {
+    typedef typename Derived::PlainObject ReturnType;
+  };
+}
+
+
+/********** MatrixBase method **********/
+
+
+template <typename Derived>
+const MatrixLogarithmReturnValue<Derived> MatrixBase<Derived>::log() const
+{
+  eigen_assert(rows() == cols());
+  return MatrixLogarithmReturnValue<Derived>(derived());
+}
+
+#endif // EIGEN_MATRIX_LOGARITHM
diff --git a/unsupported/doc/examples/MatrixLogarithm.cpp b/unsupported/doc/examples/MatrixLogarithm.cpp
new file mode 100644
index 000000000..8c5d97054
--- /dev/null
+++ b/unsupported/doc/examples/MatrixLogarithm.cpp
@@ -0,0 +1,15 @@
+#include <unsupported/Eigen/MatrixFunctions>
+#include <iostream>
+
+using namespace Eigen;
+
+int main()
+{
+  using std::sqrt;
+  MatrixXd A(3,3);
+  A << 0.5*sqrt(2), -0.5*sqrt(2), 0,
+       0.5*sqrt(2),  0.5*sqrt(2), 0,
+       0,            0,           1;
+  std::cout << "The matrix A is:\n" << A << "\n\n";
+  std::cout << "The matrix logarithm of A is:\n" << A.log() << "\n";
+}
diff --git a/unsupported/test/matrix_function.cpp b/unsupported/test/matrix_function.cpp
index 04167abfb..c2ca5d5f1 100644
--- a/unsupported/test/matrix_function.cpp
+++ b/unsupported/test/matrix_function.cpp
@@ -121,6 +121,26 @@ void testMatrixExponential(const MatrixType& A)
 }
 
 template<typename MatrixType>
+void testMatrixLogarithm(const MatrixType& A)
+{
+  typedef typename internal::traits<MatrixType>::Scalar Scalar;
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+  typedef std::complex<RealScalar> ComplexScalar;
+
+  MatrixType scaledA;
+  RealScalar maxImagPartOfSpectrum = A.eigenvalues().imag().cwiseAbs().maxCoeff();
+  if (maxImagPartOfSpectrum >= 0.9 * M_PI)
+    scaledA = A * 0.9 * M_PI / maxImagPartOfSpectrum;
+  else
+    scaledA = A;
+
+  // identity X.exp().log() = X only holds if Im(lambda) < pi for all eigenvalues of X
+  MatrixType expA = scaledA.exp();
+  MatrixType logExpA = expA.log();
+  VERIFY_IS_APPROX(logExpA, scaledA);
+}
+
+template<typename MatrixType>
 void testHyperbolicFunctions(const MatrixType& A)
 {
   // Need to use absolute error because of possible cancellation when
@@ -157,6 +177,7 @@ template<typename MatrixType>
 void testMatrix(const MatrixType& A)
 {
   testMatrixExponential(A);
+  testMatrixLogarithm(A);
   testHyperbolicFunctions(A);
   testGonioFunctions(A);
 }