Thread Subject: pinv(A)*b

A is a (m x n) sparse matrix with:
  m > n (m and n are about few thousands), and
  rank A < n (rank A ~ n-10).

b is (m x 1) vector;

What is a good way to compute pinv(A)*B? Currently I use
lsqr(A,b). I would like to have your inputs on the subject.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in
message <g9pau7$nma$1@fred.mathworks.com>...
> A is a (m x n) sparse matrix with:
> m > n (m and n are about few thousands), and
> rank A < n (rank A ~ n-10).
>
> b is (m x 1) vector;
>
> What is a good way to compute pinv(A)*B? Currently I use
> lsqr(A,b). I would like to have your inputs on the
subject.
>
> Bruno

p = colamd(A');
[Q,R]=qr(A(p,:)');
x = Q*(R'\b(p)) ;

I'm not sure how well that will handle rank deficiency,
however. lsqr(A,b) handles rank deficiency pretty well.

oh ... wait ... the code above is for m < n. Should work
for m > n too, but I'll have to double check.

I'm working on a sparse QR that I will release shortly. Can
you send me a samples of your matrices? (Both A and b?).
You can upload them to

http://www.cise.ufl.edu/~web-gfs

for user "davis". I'll add them to my collection.

"Tim Davis" <davis@cise.ufl.edu> wrote in message
<g9q3uq$ibh$1@fred.mathworks.com>...

> p = colamd(A');
> [Q,R]=qr(A(p,:)');
> x = Q*(R'\b(p)) ;
>

Hello Tim,

I have uploaded the matrix A, and b (small problem from
photogrammetry, m = 1388, n = 390, rank = 293).

In your method, how can I fix a tolerance for null space?
with pinv and lsqr, I can do somthing like this:

% load A and b
load BrunoPhotogrammetry.mat;

tic;
x1=lsqr(A,b,1e-6,1000);
toc; % Elapsed time is 0.116326 seconds

tic;
p = colamd(A');
[Q,R]=qr(A(p,:)');
x2 = Q*(R'\b(p)) ;
toc; % Elapsed time is 0.252176 seconds.

tic;
x3=pinv(full(A),1e-6)*b;
toc; % Elapsed time is 0.560355 seconds.

subplot(3,1,1); plot(x1);
subplot(3,1,2); plot(x2);
subplot(3,1,3); plot(x3);

Is there an equivalent threshold with qr? You will see that
without setting the threshold, the solution is not as I
expected.

Thanks,

Bruno

Here is another timing for A of dimension (4472 x 936):

lsqr: Elapsed time is 0.237809 seconds.
qr (Tim): Elapsed time is 5.333321 seconds.
Brute force pinv: Elapsed time is 17.796502 seconds.

qr does not work for my application because the cutting rank threshold cannot be not adjusted.

Bruno