[eigen] cache-friendly matrix inverse

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

Hi,

Just for reference, i post that here. There are 2 approaches to
cache-friendly matrix inversion, the first is to make a cache-friendly
partial LU + triangular solver, the second to do directly the
inversion by blocks with Euler 2x2 trick and rely on cache-friendly
matrix product; I tried the 2nd approach, find it attached to this
email; I don't commit to svn because the results are not very good and
show that it's probably better to do the 1st approach instead. The
performance sure is good (it spends most of its time in cache friendly
matrix product), but the precision is terrible, making it usable only
with double's, and even then with 1000x1000 matrices i get relative
errors of 1e-5.

Cheers,
Benoit

#include <Eigen/LU>
#include <Eigen/Array>

using namespace Eigen;
using namespace std;

template<typename MatrixType>
void compute_fast_inverse(const MatrixType& matrix, MatrixType *result)
{
  typedef typename MatrixType::Scalar Scalar;
  ei_assert(matrix.rows() == matrix.cols());
  const int size = matrix.rows();
  result->resize(size,size);
  const int cachefriendlysize = ei_meta_sqrt<EIGEN_TUNE_FOR_CPU_CACHE_SIZE/sizeof(Scalar)>::ret / 2;
  if(size <= cachefriendlysize) {
    matrix.partialLu().computeInverse(result);
  }
  else {
    /* Let's split M into four blocks:
     * (P Q)
     * (R S)
     * If P is invertible, with inverse denoted by P_inverse, and if
     * (S - R*P_inverse*Q) is also invertible, then the inverse of M is
     * (P' Q')
     * (R' S')
     * where
     * S' = (S - R*P_inverse*Q)^(-1)
     * P' = P1 + (P1*Q) * S' *(R*P_inverse)
     * Q' = -(P_inverse*Q) * S'
     * R' = -S' * (R*P_inverse)
     */
    const int half1 = size/2;
    const int half2 = size - half1; // careful, size may be odd
    MatrixType P = matrix.corner(TopLeft, half1, half1);
    MatrixType P_inverse;
    compute_fast_inverse(P, &P_inverse);
    MatrixType Q = matrix.corner(TopRight, half1, half2);
    MatrixType P_inverse_times_Q = P_inverse * Q;
    MatrixType R = matrix.corner(BottomLeft, half2, half1);
    MatrixType R_times_P_inverse = R * P_inverse;
    MatrixType R_times_P_inverse_times_Q = R_times_P_inverse * Q;
    MatrixType S = matrix.corner(BottomRight, half2, half2);
    MatrixType X = S - R_times_P_inverse_times_Q;
    MatrixType Y;
    compute_fast_inverse(X, &Y);
    result->corner(BottomRight, half2, half2) = Y;
    result->corner(BottomLeft, half2, half1) = - Y * R_times_P_inverse;
    MatrixType Z = P_inverse_times_Q * Y;
    result->corner(TopRight, half1, half2) = - Z;
    result->corner(TopLeft, half1, half1) = P_inverse + Z * R_times_P_inverse;    
  }
}

int main()
{
  const int size=1515;
  MatrixXd m = MatrixXd::Random(size,size);
  MatrixXd minv;
  compute_fast_inverse(m, &minv);
  //minv = m.partialLu().inverse();
  //minv = m.lu().inverse();
  cout << "relative error " << (m*minv - MatrixXd::Identity(size,size)).norm();
}