GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-09-13 01:43:13 +08:00
Code Issues Packages Projects Releases Wiki Activity

Implement matrix logarithm + test + docs.

Browse Source 
Currently, test matrix_function_1 fails due to bug #288.
This commit is contained in:
Jitse Niesen 2011-06-07 14:44:43 +01:00
parent a6d42e28fe
commit 8c8ab9ae10
6 changed files with 434 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
Eigen/src/Core/MatrixBase.h
View File

@ -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>

1
Eigen/src/Core/util/ForwardDeclarations.h
View File

@ -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>

56
unsupported/Eigen/MatrixFunctions
View File

@ -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()

340
unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h Normal file
View File

@ -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

15
unsupported/doc/examples/MatrixLogarithm.cpp Normal file
View File

@ -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";
}

21
unsupported/test/matrix_function.cpp
View File

@ -120,6 +120,26 @@ void testMatrixExponential(const MatrixType& A)
VERIFY_IS_APPROX(A.exp(), A.matrixFunction(StdStemFunctions<ComplexScalar>::exp));
}
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)
{
@ -157,6 +177,7 @@ template<typename MatrixType>
void testMatrix(const MatrixType& A)
{
testMatrixExponential(A);
testMatrixLogarithm(A);
testHyperbolicFunctions(A);
testGonioFunctions(A);
}