logdet.cpp

#include <iostream>
#include <cmath>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

int main()
{
    Matrix2f A;
    A << 2, 1, 1, 3;
    cout << "Here is the matrix A:\n" << A << endl;
    SelfAdjointEigenSolver<Matrix2f> eigensolver(A);
    
    if (eigensolver.info() != Success) abort();
    
    Vector2f l=eigensolver.eigenvalues();
    Matrix2f P=eigensolver.eigenvectors();
    Matrix2f Pinv=P.inverse();
    cout << "The eigenvalues of A are:\n" << l << endl;
    cout << "Here's a matrix whose columns are eigenvectors of A \n"
      << "corresponding to these eigenvalues:\n"
      << P << endl;
    // computing log(A)
    Matrix2f L=Matrix2f::Zero();
    for (int i = 0; i < l.rows(); ++i)
    {
      l(i)=log((double)l(i));
    }
    L.diagonal()=l;
    cout << "L=" << endl << L << endl;
    Matrix2f logA=P*L*Pinv;
    
    // direct log(det(A))
    cout << "log(det(A)) = " << log(A.determinant()) << endl;
    
    // direct tr(log(A))
    cout << "tr(log(A)) = " << logA.diagonal().sum() << endl;
    
    // using sampling, unit vectors
    Vector2f sample;
    float sum=0.0f;
    for (int i = 0; i < 2; ++i)
    {
      sample=Vector2f::Zero();
      sample(i)=1;
      sum += sample.transpose() * logA * sample;
    }
    cout << "sampled = " << sum << endl;
    
    return 0;
}