Solving Ax=b but only for the values in some rows of x

I have a linear system where:

- is square, large (m and n on the order of ), asymmetric, very sparse (around non-zeros)

- is very sparse (also - non-zeros)

Due to the properties of the problem where the matrix and vector come from, I can guarantee the system has one unique solution .

Currently, I'm solving this system using sparse LU decomposition or BiCGSTAB and this is already quite fast. But in the search for more speed, which is my primary goal, I realised that, given how the vector is used once calculated, I only actually look at a very small subset of its values. In other words, if , then I actually only need to know the values of a few (few is about ). I know which values I care about before solving . After solving the system, all other values in can be completely incorrect so long as those important are right.

Is there any faster algorithm that can take advantage of this to reduce the amount of work done when solving ?

Note: An obvious general way to calculate values directly is using Cramer's rule but I don't know how it could be applied to such large sparse matrices quickly.

[John D'Errico]

NO. You do NOT want to even consider using Cramer's rule.

[Tarek Hajj Shehadi]

You are right, I knew it would be a really bad idea

[John D'Errico]

Oh. I would not even call that matrix really very sparse. 0.05% non-zeros is roughly 50 non-zeros per row. The sparsity of b is pretty much irrelevant.

[John D'Errico]

Is A fixed, and you are solving this problem repeatedly for multiple right hand sides?

[Tarek Hajj Shehadi]

Indeed, I am solving for one matrix A and one vector b

[Image Analyst]

1. How long does it actually take? How many minutes or hours? (You can use tic and toc to time code)
2. Can you attach your variables in a .mat file so people can try things with it?

[John D'Errico]

Oh. Multiple b, but fixed A? Then the answer should be easy. Let me write it up as an answer.

[Tarek Hajj Shehadi]

Edit: I am actually taking A as a whole and b as a whole

Answer 1

[John D'Errico]

We are given a problem where A is a large sparse matrix. We wish to compute the solution for the problem

A*X == b

where b will vary, but A is fixed. The thing is, we want to solve for only a few elements of X, a small subset of them.

The trick would seem to be to effectively compute the corresponding rows of the inverse of A. For example, I'll create a very small problem, with A as a 6x6 matrix, so we can see what is happening. Any larger, and it will overflow the window, and what matters is the idea, not the implementation on your real problem.

A = randn(6)
A =
      0.91539       1.2244      0.23128       2.5526    -0.044323      0.07913
      -1.7792       1.1284     -0.54674       -1.255     -0.53844      -2.2881
     -0.71475      0.44224       0.4884       1.0128     -0.53548     -0.43906
    -0.092051        1.516       -1.918      -1.2029      -1.1266     -0.73841
      0.72988    -0.048539     -0.53597      -1.0549     -0.53726     -0.99856
    -0.059816      0.53947       1.7059     -0.21087     -0.86017      0.75975

Now suppose we wish to solve the problem A*X == b, but we care only about elements 1 and 3 of X?

The idea is to compute rows 1 and 3 of inv(A). The following very simple code does exactly that.

function Ainv = Ainvrows(A,rind)
  % generate the rows of the inverse of A, for a small subset of rows
  %   A - a square matrix, presumed to be non-singular, presumed sparse
  %   rind - integer vector listing the rows of the inverse of A that are of interest
  
  Asize = size(A);
  n = Asize(1);
  if n ~= Asize(2)
    error('A must be square')
  end
  
  nr = numel(rind);
  delta = sparse(rind,1:nr,1,n,nr);
  
  % the trick now is to solve the problem A'\delta. That effectively
  % extracts only the rows of the inverse of A that we wish to see.
  Ainv = (A'\delta)';

If you prefer to use some other solver than a sparse backslash, perhaps bicgstab, just modify that last line. I think bicgstab is not set up to handle multiple right hand sides. So you would need to have a loop in there.

Did it work on my example matrix?

inv(A)
ans =
       0.2899    -0.056298     -0.40007    -0.095469      0.47353     0.098636
      0.45274      0.39894     -0.82201      0.11297     -0.36017      0.31568
      0.13989      0.23266     -0.25796     -0.33976      0.12684      0.37352
     0.069612     -0.17474      0.54229    -0.004086   -0.0012222     -0.22568
      0.38001      0.51357      -1.1788     -0.33684     -0.42659    -0.062176
     -0.16319     -0.27714    -0.052687      0.29264     -0.47509      0.12814
Ainv13 = Ainvrows(A,[1 3])
Ainv13 =
       0.2899    -0.056298     -0.40007    -0.095469      0.47353     0.098636
      0.13989      0.23266     -0.25796     -0.33976      0.12684      0.37352

Now, once we have the necessary rows of the inverse of A, all it takes to compute the value of elements 1 and 3 of the solution is a matrix multiply. For example:

b = rand(6,1);

A\b
ans =
      0.10635
      0.10011
      0.19457
    -0.073614
     -0.14575
     -0.17135
     
Ainv13*b
ans =
      0.10635
      0.19457

As you can see, it returns the 1 and 3 elements of the solution.

This code would be efficient if you have MANY vectors b, more than the number of elements of X that you care about the values. It will need to solve for as many right hand sides in delta as there are elements you will want to see the value for.

[Matt J]

This code would be efficient if you have MANY vectors b

Unfortuantely, the OP has backpedalled on that. There is only one b.

Answer 2

[Matt J]

Indeed, I am solving it for multiplie right hand sides with A fixed

If so, then you should be solving with multiple concatenated b,

x=A\[b1,b2,b3,...,bn];

This will make it so the decomposition of A is reused for multiple b_i.

Also, if you sort the columns of A so that the x(i) of interest are the final m variables, then an LU decomposion could be truncated to,

[L,U]=lu(A);
U=U(end+1-m:end, end+1-m);

Q=L\[b1,b2,b3,...,bn];
Q(1:end-m,:)=[];

x=U\Q; %only the x(i) of interest

[Tarek Hajj Shehadi]

I truly apologize for editing my previous comment, I am solving for one matrix A and one matrix b

[Matt J]

Still, the technique of truncating U to consider only the x(i) of interest still applies.

[Tarek Hajj Shehadi]

Alright, I shall take a look at it thank you very much!!