CODEWITHSUNDEEP
matlaboptimization
user34790
user34790

Reputation: 2036

Issues with using pcg

I am trying to use preconditioned conjugate gradient in matlab to speed things up. I am using this iterative method because backslash operator was too time consuming. However, I am having some issues. For e.g. I want to solve this system of linear equations given by

Ax=B

where A is a sparse positive definite matrix and B is a matrix as well. In matlab I can do that simply by

x= A\B

However, if I use pcg function, then I would have to loop over all the columns of B and solve individual 

x(:,i)=pcg(A,B(:,i))

This loop will take more time than x=A\B. If I consider just a single column as b instead of matrix B, then pcg works faster than the backslash operator. However, if I consider the whole matrix B, then pcg is slower than backslash operator. So there is no point of using pcg.

Any suggestions guys?

When using the method as suggested by Mattj, it shows the following error

Error using iterapp (line 60)
user supplied function ==>
@(x)reshape(repmat(A*x,1,nb),[],1)
 failed with the following error:

 Inner matrix dimensions must agree.

Error in pcg (line 161)
r = b -
iterapp('mtimes',afun,atype,afcnstr,x,varargin{:});

Upvotes: 2

Views: 957

Answers (2)

Matt J
Matt J

Reputation: 1137

However, if I consider the whole matrix B, then pcg is slower than backslash operator. So there is no point of using pcg.

That does make a certain amount of sense. When backslash solves the first set of equations A*x=B(:,1), it can recycle pieces of its analysis to the later columns B(:,i), e.g., if it performs an LU decomposition of A.

Conversely, all the work that PCG applies to the different B(:,i) are independent. So, it very well might not make sense to use PCG. The one exception to this is if each B(:,i+1) is similar to B(:,i). In other words, if the columns of B change in a gradual continuous manner. If so, then you should run PCG in a loop as you've been doing, but use the i-th solution x(:,i) to initialize PCG in the next iteration of the loop. That will cut down on the total amount of work PCG must perform.

Upvotes: 0

Matt J
Matt J

Reputation: 1137

I think we need to see more data on your timing tests and the dimensions/sparsities of A and B, and better understand why pcg is faster than mldivide. However, you can implement what you're after this way,

 [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);

Upvotes: 0

Related Questions