your question amounts to ask:
what is the difference between solving a square nonsingular system of
equations Ax=b with an LU with partial pivoting factorization (say LU-GEPP
for Gaussian Elimination with Partial Pivoting) of the matrix A (ZGESV) or
a QR factorization using Householder reflections (say QR-HH) of the matrix
1- The cost of LU-GEPP is 2/3 n^3 while the cost of the QR-HH
factorization 4/3 n^3 (for square n-by-n matrix). So LU-GEPP (DGESV) is
expected to be twice as fast as QR-HH (DGELS).
2- QR-HH is unconditionnally backward stable. This means that's for any
input matrix A, the computed solution _x_ from QR-HH is the exact solution
of a slightly perturbed problem. In other words, _x_ is such that:
( A + Delta_A ) _x_ = ( b + delta_b )
where || Delta_A || <= small_constant * machine_eps * || A ||
and || delta_b || <= small_constant * machine_eps * || b ||.
Such inequality (backward stability) is the minimum one ask for any
LU-GEPP has such a good behavior for all the matrix that we have
encountered except a few. So LU-GEPP is not backward stable in theory but
is in practice. So in practice both methods will provide you with a
backward stable solution.
3- if the system is close to singular (we say that it is ill-conditionned)
then you can have a lot of backward stable solutions. Basically the
|| _x_ - true_solution || <= condition_number_A * bwd_error
Both LU_GEPP and QR-HH will give you a small bwd_error but if the
condition number of A is large then the computed solution _x_ can be
pretty far from the true solution. Moreover there is no reason for the
computed solution of LU-GEPP to be the same as the computed solution of
QR-HH and in practice you will observe that they differ by a factor
machine_precision * condition_number_A.
However you can not tell one solution is better than the other. They are
both backward stable solutions. So they are both 'acceptable' answer.
4- If the matrix is ill-conditionned, instead of using ZGELS, you can use
other algorithms. They will filter the small singular values out of your
initial matrix A, once those are removed your system becomes
overdetermined and we solve a least squares problem on it. Those routines
are ZGELSY (QR with column pivoting), ZGELSD (SVD with Divide and
Conquer), ZGELSS (SVD with QR).
If you want to know more, you might want to look at a dense linear algebra
book. (For example: Demmel, Trefethen and Bau, or Golub and Van Loan).
On Fri, 24 Aug 2007, Sorkin Dima wrote:
What is the difference between solving a _square_
system of equations with ZGELS vs. ZGESV.
Can they give different results ?
What happens when matrix is very close to singular ?
Thank you, regards,
Lapack mailing list