how to solve A*x=b, if A is a singular matrix;

hello, everyone! I am using hte following code to build the matrix New_P and New_Pg, in order to solve New_P*h=New_Pg; Because, the matrix New_P has a possibility of singular, so I use the pinv function. However, it turns out this make the computation very solw, escepically when the size of the matrix growing largely. I wonder if there exist any way to replace the pinv function or to optimize the code. Thank you very much!

S=size(P,1);
if n<=10
   New_P=P(sub2ind(size(P),repmat((1:size(P,1))',1, size(P,2)), repmat(1:size(P,2), size(P,1), 1), repmat(U0, 1,size(P,2))));
else
   New_P=zeros(S,S);
    for s=1:S
      New_P(s,:)=P(s,:,U0(s));
    end;
end;
New_Pg=g(:,1);
New_P=cat(2,ones(S,1),eye(S)-New_P);
New_P(S+1,t+1)=1;  % set state 1 as the reference state;
New_Pg=[New_Pg;0];
if rcond(New_P)>=eps('double')
   h0=New_P\New_Pg;
else
   h0=pinv(New_P)*New_Pg;   % pinv is very unefficient, but it is mostly used in sigular matrix. 
end;
h0=h0(2:end);

Answer 1

[John D'Errico]

Sorry. You have a fast solution, or you have one that is robust to singular systems. Sometimes you cannot have everything you want. :)

No. You cannot optimize pinv.

However, in this case, you may be lucky. One alternative (that may have some quirks of its own of course, because it is an iterative solver) is lsqr. lsqr should generate the same solution as I recall as pinv, when applied to a singular system.

For example...

A = ones(5);
b = ones(5,1);

A\b
Warning: Matrix is singular to working precision. 

ans =
   NaN
   NaN
   NaN
   NaN
   NaN

pinv(A)*b
ans =
          0.2
          0.2
          0.2
          0.2
          0.2

lsqr(A,b)
lsqr converged at iteration 1 to a solution with relative residual 4e-16.
ans =
          0.2
          0.2
          0.2
          0.2
          0.2

You can turn off the displayed response if you return a second output argument for lsqr.

A = rand(2000,2000);
b = rand(2000,1);

tic,x = A\b;toc
Elapsed time is 0.227208 seconds.

tic,x = pinv(A)*b;toc
Elapsed time is 6.194447 seconds.

tic,[x,flag] = lsqr(A,b);toc
Elapsed time is 0.122928 seconds.

[Jason]

Thank you very much, it helps for reducing computational time! Besides, do you know if there a method in MATLAB, where we can set the maxime time limit. when the program achieve the time limit, it gives outputs no matter it gets optimal or not.

[John D'Errico]

No. There is no time limit on lsqr. In fact, due to the use of multiple processors on systems today, you could easily find that two different calls to LSQR will take significantly different amounts of time to execute, just because some of those CPUs were off doing something else at the time. Checking your mail on the side? Using your web browser? These things suck CPU cycles away from MATLAB.

There is a limit on the number of iterations taken, but bigger problems MAY take more iterations, and bigger problems will certainly take more time per iteration. I have no idea what is a reasonable number of iterations to use as a max for LSQR either.