Block conjugate gradient in matlab

4
votes

I want to solve a system of linear equations, AX = B, where A is sparse and positive definite. B is a matrix rather than a column vector. So I have to solve multiple system of linear equations (with multiple right hand sides). How can I use conjugate gradient for this in Matlab?

I can use the one that works for a column vector B.

2
matlaboptimizationlinear-algebra
As mentioned in your related post, you should probably not do this when A is small enough to be handled by MLDIVIDE and when B has many columns. MLDIVIDE will be more efficient. - Matt J

2 Answers

2
votes

Supply B is a column vector B(:) instead. Also, supply an efficient implementation of A in functional form,

 [ma,na]=size(A);
 [mb,nb]=size(B);
 afun=@(x)  reshape(A*reshape(x,na,[]),[],1);

 X=pcg(afun,B(:));

  X=reshape(X,na,nb);
2
votes

Solving the linear system of equations AX=B where B is a matrix will result in X also being a matrix. However, the columns of X will each be a solution to the linear system where the right hand side is the corresponding column of B.

So, if you already have a conjugate gradient function that works on a column vector B (which in Matlab is x = pcg(A,b);), then you can find the solution in the case where B is a matrix by looping over the columns:

X = zeros(size(A,2), size(B,2));
for i=1:size(B,2)
    X(:,i) = pcg(A,B(:,i));
end