CGNE and underdetermined least squares

From: spasmous (spasmous_at_yahoo.com)
Date: 03/16/05 

Date: 16 Mar 2005 13:36:46 -0800

I'm solving an underdetermined system of equations Ax = b. The matrix
is around 1e4 x 1e5 or so and sparse, ~20 non-zeros per row. I've had
success with lsqr for solving this problem, however lsqr solves the
normal equations

A'*A*x = A'*b

which is a 1e5 x 1e5 system. It seems more natural in this case to
solve the other form of the normal equations

A*A'*y = b followed by x = A'*y

since this is a 1e4 x 1e4 system.

I coded in CGNE to do this but it seems to suffer from roundoff - or at
least tend to infinity given the slightest encouragement ;) while lsqr
still manages to return meaningful results.

My code is below, can anyone suggest an improvement or point me towards
a "cgne" variant of lsqr?

Thanks.

function x = cgne(A,b,tol,kmax)

% Conjugate gradient on the normal equations.
% "Naive" implementation more susceptible to
% roundoff than lsqr.
%
% A is an m x n matrix
% b is an m vector
% tol is the tolerance
% kmax is the maximum iterations
% damp is a regularization term

% size
[m n] = size(A);

% initialise
y = zeros(m,1);
r = b;
d = r;
dnew = r'*r;
d0 = dnew;

% begin iterations
for k = 1:kmax

        q = A*(A'*d);
        alpha = dnew./real(d'*q);
        y = y + alpha.*d;
        r = r - alpha.*q;
        dold = dnew;
        dnew = r'*r;
        if dnew<tol.^2*dold
                break
        end
        beta = dnew./dold;
        d = r + beta.*d;

end

% recover x
x = A'*y;

Loading