Linear Least-square solvers

[anslab]

Hi everyone,

I am fairly new here, and this post may or may not be in the right section. I apologize if it's not where it's supposed to be.

In my work I need to solve a linear least square unconstrained problem. I found three LAPACK routines that can be useful :

Code: Select all

SGELS(..)
SGELSD(..)
SGELSY(..)

I compared them in terms of precision and execution time, and it allowed me to choose one, but now i am looking for a more in depth understanding.


Can any one explain me the differences between these methods?

Thanks
A.

[admin]

Did you take a look at the LAPACK user guide?
http://www.netlib.org/lapack/lug/node1.html
Julie

[anslab]

Yes I did, I also browsed the sources of the aforementionned functions, but I must admit I'm still in the dark. I also asked my good friend wikipedia, but he wasn't of much help. Maybe you can point me to a paper or any other good ressource?

[Julien Langou]

* xGELS solves the linear least squares problem using a QR factorization.
* xGELSY solves the linear least squares problem using a QR factorization with column pivoting.
* xGELSD solves the linear least squares problem using the singular value decompostion.

In term of speed, you should get xGELS faster than xGELSY faster than xGELSD.
In term of functionality, xGELS only handles nonsingular matrices, xGELSY and xGELSD handles any matrix. xGELSY might have some problems for some pathological matrices. xGELSD works for all matrices.

J

[anslab]

Okay,

When i did some testing, i noticed that there is a discrepancy in precision, i.e. xGELSD seems to find a more accurate solution than the other two. Is there any theoretical explanation, or is it implementation-related?

[Julien Langou]

1) What do you mean by "xGELSD seems to find a more accurate solution than the other two"? Measuring the "accuracy" of a linear least squares solver is tricky ... Do you check which one has the smallest 2-norm residual ( || b - A x || )? Do you check which one has the most orthogonal residual (i.e. smallest A^T (b- Ax ) ) ? I assume the first. (Which is fine.)

2) In term of accuracy, you should expect xGELSD to be better than the two other two when the problem becomes ill-conditionned. (The order should actually be: xGELS is the worse, xGELSY should come second, then xGELSD.) Of course there is a price (i.e., time) to this improved accuracy. In term of time, you should expect the reverse order.

Cheers,
Julien.

[anslab]

1) Yes, i did the first solution. I generated a bunch of system from a given solution in order to measure the accuracy (in the 2-norm sense) and the execution time.

What I need to understand now is, why xGELSD is more accurate that the other on the same problem than the other solvers?

Thanks for taking the time to answers my questions, I think you just earned a line of thanks in my PhD.