Hello,

I am putting together a short note to answer your problem because it's not the

fist time that we have users asking the question. So just to answer your questions

quickly before the note, here are some facts.

In the overdetermined nonsingular case (or numerically nonsingular if you prefer),

* the solution returns by LAPACK routines DGELSY, DGELSS, DGELSD,

* the solution returned by A\b, lsocv(A , b) (Matlab)

* the solution retruned by pinv(A)*b

are all the same.

In the overdetermined singular case (or numerically rank deficient), this is your case,

the solution returned by LAPACK routine DGELSY, DGELSS, and DGELSD is the same

as the solution return by Matlab pinv(A)*b but is not the same as the solution returned

by Matlab A\b or lscov(A,b). This is normal. In the rank deficient case, there is an inifinite

number of solutions and so Matlab and LAPACK can both give you different answers

while being both correct.

The strategy/philosophy of Matlab is to represent the solution with as few variables as

possible. This is why you have an exact 0 in your x. Only two variables are needed in

your case.

The strategy/philosophy or LAPACK (or x =pinv(A)*b in Matlab) is to minimize the

two-norm of the solution x. Then x is fully populated but has minimum norm.

In any case the residual r = b-Ax is the same (even though the x's are different) and is

guaranteed to have minimal norm. So both LAPACK and Matlab solve the Least Squares problem.

Note that Matlab's approach is slightly faster than LAPACK.

OK so that's the first point. The second point is that I think you are not computing the

estimated standard error of the solution (stdx) in a stable manner for the solution of

LAPACK while Matlab function lscov computes stdx stably.

My suspicion is that you are doing, something like

- Code: Select all
` mse = norm(b-A*x)^2/(m-n);`

S = inv(A'*A)*mse;

stdx = sqrt(diag(S));

This previous code works fine for well-conditionned matrix.

If A is ill-conditionned (a fortiori singular), the inv(A'*A) is going to blow you in the face

and I think this is what you are observing.

Try the following code, you should feel better with what LAPACK is doing after that.

Cheers,

Julien

- Code: Select all
` m = 4;`

n = 3;

A = [1 1 17; 1 2 18; 1 3 19; 1 4 20] ;

b = [18.75715 18.83939 18.91538 18.98599]' ;

x = [ +0.140107 -1.082733 +1.158984 ]'; % This is the LAPACK solution

rankA = 2;

mse = norm(b-A*x)^2/(m-rankA);

stdx = zeros(n,1);

[u,s,v] = svd(A);

s = s(1:rankA,1:rankA);

v = v(:,1:rankA);

for i=1:n,

stdx(i) = norm( s \ (v(i,:)') ) * sqrt(mse);

end

fprintf('rankA = %d\n',rankA);

fprintf('x = [ %+f %+f %+f ]''\n',x(1),x(2),x(3));

fprintf('mse = %f\n',mse);

fprintf('stdx = [ %+f %+f %+f ]''\n',stdx(1),stdx(2),stdx(3));

fprintf('\n');