I need to solve many systems of linear equations with the same matrix. So, I done LU-factorisation of the matrix, and calling .solve for each right-hand side.

It is very slow.

But for each right-hand side I exactly know which variables will be zero (I have about 90% of zeros in each solution). How can I use this information to speed up solution?

What's the average size of your matrices?

Note that you can solve for multiple right hand side at once using a matrix:

Code: Select all

MatrixXd A(N,N), B(N,M), X(N,M);
PartialPivLU<MatrixXd> lu(A);
X = lu.solve(B);

that is significantly faster than solving for each column:

Code: Select all

for(i)
  X.col(i) = lu.solve(B.col(i));

Now, regarding your special use case with 90% of the result being 0, I don't see an obvious solution with a precomputed LU. You could restart from scratch by forming a small matrices like this:

Code: Select all

M = number of non zeros in X
VectorXi indices(M);
int k=0;
for(i=0..N-1)
  if x_i != 0 then indices(k) = i;
A1.resize(M,M);
B1.resize(M);
for(i=0..M-1)
  B1(i) = B(indices(i));
  for(j=0..M-1)
    A1(i,j) = A(indices(i),indicies(j));
   
X1 = A1.lu().solve(B1);
for(i=0..M-1)
  X(indices(i)) = X1(i);

in the near future these for loops could be avoided by simply writing:

Code: Select all

B1 = B(indices);
A1 = A(indices,indices);
X(indices) = X1;

Another solution would be to explicitly form the inverse of A and then do:

Code: Select all

X(indices) = Ainv(indices,All) * B;

Currently you have to manually form Ainv(indices,All) using a for loop. However, explicitly computing the inverse is not recommended for numerical precision and other reasons.

Thank you for you answer!

I trying to implement Kriging interpolation. I have about N=2000 data points, so I have 2000×2000 matrix A. I need to solve Ax=f for every point where interpolation should be calculated (about 1 million of points).

The program works slowly (tens of minutes).

I cannot see the reason for speed gain if I will solve multiple right-hand sides at ones. Can you explain this?

I searched the Internet and found that usually all data points are not considered for each interpolation point. Instead small amount (for example, M=20) of nearest data points is taken into account, and small system of equations is constructed and solved (as you recommended) for each interpolation point.

The problem is that these small systems of equations are solved in cubic time O(M^3). So, this popular method quickly becomes bad as M grows. If M>150 it will become slower than my current O(N^2) for each interpolation point.

My idea is to construct LU decomposition for full matrix A, and then use it to solve subsytems of A in O(M^2). If formulas of LU decomposition will not allow me to do such trick, then I will try to find some other decomposition which can do this.

SAnton wrote:I trying to implement Kriging interpolation. I have about N=2000 data points, so I have 2000×2000 matrix A. I need to solve Ax=f for every point where interpolation should be calculated (about 1 million of points).

The program works slowly (tens of minutes).

I cannot see the reason for speed gain if I will solve multiple right-hand sides at ones. Can you explain this?

Please, try to solve them per packet of about 2000 vectors stored in a matrix and you will see! For vectors of size 2000 I expect a huge speedup since according to our benchmark matrix-vector operations performs very badly for such size because it is not possible to take advantage of the caches.

The reason is that matrix-matrix operations run at a much higher rate because of better cache usage. See e.g. our benchmark: http://eigen.tuxfamily.org/index.php?title=Benchmark