[eigen] Problem with LU and Cholesky inversion

[Timothy Hunter <tjhunter@xxxxxxxxxxxx>]

Hello list,

after running into numerical trouble when implementing a convex optimizer, I wrote this small test (see attachment). It inverses a positive semidefinite matrix using all the methods available in eigen. Can you tell me if I am doing a mistake here or if it exposes a bug in LU and Cholesky modules?

Thanks in advance

Tim Hunter

--

Timothy Hunter

Student (Stanford University)

T. 404 421 3075

#include <Eigen/Core>
#include <Eigen/LU>
#include <Eigen/QR>
#include <Eigen/Cholesky>
#include <Eigen/Array>
#include <Eigen/SVD>
#include <iostream>
using namespace Eigen;
using namespace std;

template<typename Scalar, int N, int M>
void testLU1(int n, int m)
{
  typedef Matrix<Scalar, N, 1> VectorType;
  typedef Matrix<Scalar, N, N> MatrixType;
  typedef Matrix<Scalar, M, 1> CVectorType;
  typedef Matrix<Scalar, M, N> CMatrixType;
  typedef Matrix<Scalar, N, M> DMatrixType;
  
  MatrixType hessian=MatrixType::Random(n,n);
  hessian*=hessian.transpose();
  CMatrixType A=CMatrixType::Random(m,n);
  
  DMatrixType V1(n,m);
  Eigen::LU<MatrixType>  hessianLU(hessian);
  hessianLU.solve(A.transpose(),&V1);
  std::cout<<"\nLU-------\n"<<hessian*V1-A.transpose()<<"\n-------\n";
  
  DMatrixType V2(n,m);
  Eigen::CholeskyWithoutSquareRoot<MatrixType>  hessianChol(hessian);
  V2=hessianChol.solve(A.transpose());
  std::cout<<"\nChol-------\n"<<hessian*V2-A.transpose()<<"\n-------\n";
  
  DMatrixType V3=hessian.inverse()*A.transpose();
  std::cout<<"\nINV-------\n"<<hessian*V3-A.transpose()<<"\n-------\n";
  
  DMatrixType V4(n,m);
  Eigen::SVD<MatrixType>  hessianSVD(hessian);
  hessianSVD.solve(A.transpose(),&V4);
  std::cout<<"\nSVD-------\n"<<hessian*V4-A.transpose()<<"\n-------\n";

  std::cout<<"\n-------\n"<<V3-V2<<"\n-------\n";
}

template<typename Scalar, int N, int M>
void testChol1(int n, int m)
{
  typedef Matrix<Scalar, N, 1> VectorType;
  typedef Matrix<Scalar, N, N> MatrixType;
  typedef Matrix<Scalar, M, 1> CVectorType;
  typedef Matrix<Scalar, M, N> CMatrixType;
  typedef Matrix<Scalar, N, M> DMatrixType;
  
  MatrixType hessian=MatrixType::Random(n,n);
  hessian*=hessian.transpose();
  CMatrixType A=CMatrixType::Random(m,n);
  
  DMatrixType V2(n,m);
  Eigen::CholeskyWithoutSquareRoot<MatrixType>  hessianChol(hessian);
  
  std::cout<<"\nChol-------\n"<<hessianChol.matrixL()<<"\n-------\n";
  std::cout<<"\nChol-------\n"<<hessianChol.vectorD()<<"\n-------\n";
  MatrixType chol2=hessianChol.matrixL()*(hessianChol.vectorD().asDiagonal())*hessianChol.matrixL().transpose();
  std::cout<<"\nChol-------\n"<<hessian-chol2<<"\n-------\n";
  V2=hessianChol.solve(A.transpose());
  std::cout<<"\nChol-------\n"<<hessian*V2-A.transpose()<<"\n-------\n";
}

int main(int argc, char **argv)
{
  std::cout<<"##Static (10,2)##\n";
  testChol1<double, 10, 2>(10,2);

  std::cout<<"##Dynamic (10,2)##\n";
  testChol1<double, Eigen::Dynamic, Eigen::Dynamic>(10,2);

  std::cout<<"##Static (10,2)##\n";
  testLU1<double, 10, 2>(10,2);

  std::cout<<"##Dynamic (10,2)##\n";
  testLU1<double, Eigen::Dynamic, Eigen::Dynamic>(10,2);

/*  std::cout<<"##Static (5,2)##\n";
  testLU1<double, 5, 2>(5,2);

  std::cout<<"##Dynamic (5,2)##\n";
  testLU1<double, Eigen::Dynamic, Eigen::Dynamic>(5,2);*/
  return 0;
}