I have to calculate the full inverse of a sparse matrix, and ggael told me on IRC to take a look at test/sparse_solvers.cpp.

I did, and also read the API docs of SparseLU - which gives me .solve, but I found no .inverse method.

So what's the easiest way to get the inverse of a sparse matrix?

(FYI my matrices are about 7k X 7k, with perhaps 50k non-zero complex elements (neither hermitian nor tri-diagonal). The inverse is a full matrix. Currently boost/ublas takes about 11 hours for that task. Anything that's at last 1.5 times faster is worth switching to Eigen for me)

well, unfortunately, none of the supported libraries (SuperLU and UmfPack) provides a solve function (AX=B) with X and B sparse. Therefore there is no efficient way to extract the full inverse matrix. What you can still do is to use dense matrices for the solve step:

Code: Select all

// build your sparse matrix A, then:
SparseLU >,SuperLU> slu(A);
MatrixXcd I = MatrixXcd::Identity(size,size), inv(size,size);
slu.solve(I, &inv);

of course since your matrix is large, you should extract the inverse per column:

Code: Select all

// build your sparse matrix A, then:
SparseMatrix > matInv(size,size);
SparseLU >,SuperLU> slu(A);
VectorXcd base(size), invCol(size);
matInv.startFill();
for (int i=0; i<size; ++i)
{
  base = VectorXcd::UnitX(i,size), invCol(size);
  slu.solve(I, &invCol);
  for (int j=0; j<size; ++j)
    if (ei_abs(invCol[i]) < epsilon)
       matInv.fill(j,i) = invCol[i];
}
matInv.endFill();

this should still be quite slow but much much faster than your current solution !

ggael wrote:

Code: Select all

// build your sparse matrix A, then:
SparseLU >,SuperLU> slu(A);
MatrixXcd I = MatrixXcd::Identity(size,size), inv(size,size);
slu.solve(I, &inv);

There's no need to actually construct the identity matrix, you can pass the Identity expression directly to solve():

Code: Select all

// build your sparse matrix A, then:
SparseLU >,SuperLU> slu(A);
MatrixXcd inv(size,size);
slu.solve(MatrixXcd::Identity(size,size), &inv);

[quote='ggael' pid='44563' dateline='1234372241']
of course since your matrix is large, you should extract the inverse per column:

Code: Select all

// build your sparse matrix A, then:
SparseMatrix > matInv(size,size);
SparseLU >,SuperLU> slu(A);
VectorXcd base(size), invCol(size);
matInv.startFill();
for (int i=0; i > matInv(size,size);
SparseLU >,SuperLU> slu(A);
VectorXcd base(size), invCol(size);
matInv.startFill();
for (int i=0; i<size; ++i)
{
  base = VectorXcd::UnitX(i,size);
  slu.solve(base, &invCol);
  for (int j=0; j<size; ++j)
    if (ei_abs(invCol[i]) < epsilon)
       matInv.fill(j,i) = invCol[i];
}
matInv.endFill();

bjacob wrote:There's no need to actually construct the identity matrix, you can pass the Identity expression directly to solve():

well, superLU needs a plainmatrix, so yes, currently you need to pass an expression with DirectAccessBit flags (e.g., Matrix or Block)

This doesn't compile in my brain I think you meant this:

Code: Select all

// build your sparse matrix A, then:
SparseMatrix > matInv(size,size);
SparseLU >,SuperLU> slu(A);
VectorXcd base(size), invCol(size);
matInv.startFill();
for (int i=0; i<size; ++i)
{
  base = VectorXcd::UnitX(i,size);
  slu.solve(base, &invCol);
  for (int j=0; j<size; ++j)
    if (ei_abs(invCol[i]) < epsilon)
       matInv.fill(j,i) = invCol[i];
}
matInv.endFill();

thanks !

ggael wrote:well, superLU needs a plainmatrix, so yes, currently you need to pass an expression with DirectAccessBit flags (e.g., Matrix or Block)

Oops, I didn't see that you were working with the SparseLU (now i see), as you said we could use a dense matrix for the solve step, i jumped to the "conclusion" that you were using a dense LU.

[quote='bjacob' pid='44591' dateline='1234379370']
This doesn't compile in my brain ;) I think you meant this:

Code: Select all

// build your sparse matrix A, then:
SparseMatrix > matInv(size,size);
SparseLU >,SuperLU> slu(A);
VectorXcd base(size), invCol(size);
matInv.startFill();
for (int i=0; i epsilon)"?

I've tried the code, and as a check multiplied A with the inverse, and the result isn't the unit matrix. This is what I've tried:

[code]
#include
#include
#include
#include

using namespace std;
using namespace Eigen;

typedef double num;
typedef complex cnum;

int Nx = 4;
int size = 2 * Nx * Nx;
double epsilon = 1e-8;

int main(int argc, char** argv) {
    SparseMatrix A(size, size);
    {
        RandomSetter > setter(A);
        for (int i = 0; i  matInv(size, size);
    SparseLU,SuperLU> slu(A);
    VectorXcd base(size), invCol(size);
    matInv.startFill();
    for (int i=0; i epsilon)
                matInv.fill(j,i) = invCol[i];
    }
    matInv.endFill();
    cout  result(size, size);
    result = A * matInv;

    cout << result << endl;
    return 0;
}


Any idea what might be wrong?

(In my real application the matrix has more items set than the upper triangle, this was just to get a matrix that's surely invertible)

moritz wrote:I got it to run by substituting UnitX(i, size) with Unit(size, i)

Right

moritz wrote:Shouldn't it be "if (ei_abs(invCol[i]) > epsilon)"?

No, it's really < because you want to zero the smaller entries, not the bigger ones.

moritz wrote:I've tried the code, and as a check multiplied A with the inverse, and the result isn't the unit matrix.

Did you use a fuzzy comparison? How far away from the identity are you?

[quote='bjacob' pid='44801' dateline='1234444490']

moritz wrote:Shouldn't it be "if (ei_abs(invCol[i]) > epsilon)"?

No, it's really " most elements are either (1, 0) or (2, 0) even though they should be zero, so it's very far away.

moritz wrote:Now I'm really confused - in the following no elements are zeroed, but filled.

Woooops ! Sure you're right, sorry! indeed it's >.

moritz wrote:I've just printed it out, and with "" most elements are either (1, 0) or (2, 0) even though they should be zero, so it's very far away.

Try this actually (j instead of i) :

Code: Select all

    if (ei_abs(invCol[j]) > epsilon)


well, I had better to say that the code I wrote was just to give you an idea about how to do it, and indeed the comparison operator was wrong and the index as well, I did not expected you to simply copy paste it as it!

I think this one should be correct (with i for the row and j for the col => less trouble in mine mind

Code: Select all

for (int j=0; j epsilon)
                matInv.fill(i,j) = invCol[i];
    }

I'd be interested to know how long this takes compared to boost/ublas. If I remember correctly you could also compute inverses blockwise or by cholesky / lu, right? I don't remember in which situations you should do what, though.

To me, some benchmark on factorization and creating inverse matrices of sparse matrices with different sparsity levels would be helpful as well. I guess, as time is approaching...


This is a little offtopic, but when using UMFPACK to create the inverse matrix I get a compilation error.

In function `Eigen::SparseLU, 4>::~SparseLU()':
(.text._ZN5Eigen8SparseLUINS_12SparseMatrixIdLi0EEELi4EED1Ev[Eigen::SparseLU, 4>::~SparseLU()]+0x15): undefined reference to `umfpack_di_free_numeric'

Changing

Code: Select all

extern "C" {
  #include "umfpack.h"
}


in Eigen/Sparse does not seem to help...

The same error goes to Cholmod. Using gcc 4.3.1 Any ideas?

Seb wrote:To me, some benchmark on factorization and creating inverse matrices of sparse matrices with different sparsity levels would be helpful as well. I guess, as time is approaching...

that's very difficult because the number of options is really HUGE.

This is a little offtopic, but when using UMFPACK to create the inverse matrix I get a compilation error.

In function `Eigen::SparseLU, 4>::~SparseLU()':
(.text._ZN5Eigen8SparseLUINS_12SparseMatrixIdLi0EEELi4EED1Ev[Eigen::SparseLU, 4>::~SparseLU()]+0x15): undefined reference to `umfpack_di_free_numeric'

Changing

Code: Select all

extern "C" {
  #include "umfpack.h"
}


in Eigen/Sparse does not seem to help...

The same error goes to Cholmod. Using gcc 4.3.1 Any ideas?

what is your umfpack version ? (here 5.2) perhaps you can check that the library really contains the symbol:

Code: Select all

$ strings /usr/lib64/libumfpack.so.5.2.0 | grep umfpack_di_free_numeric

should return:

Code: Select all

umfpack_di_free_numeric

if not, then you have a problem with your umfpack library.

ggael wrote:

Seb wrote:To me, some benchmark on factorization and creating inverse matrices of sparse matrices with different sparsity levels would be helpful as well. I guess, as time is approaching...

that's very difficult because the number of options is really HUGE.

True. But it would be interesting to understand the overhead by constructing the actual inverse...

if not, then you have a problem with your umfpack library.

True. Was compiled on another system. Thanks for letting me bugging you again!

Thank you very much everyone, it now works as I want it. The execution time dropped from 11 hours to about half an hour, and I'm rather happy with that improvement