Hi,
In the catalogue of decompositions offered by Eigen, two types of LU decomposition are provided : one with partial pivoting and one with full pivoting. The same documentation states that the accuracy and reliability of the partial pivoting strategy depends upon the condition number whereas the full pivoting strategy has a reliability that is “proven”. I believe that those statements are misleading. In “Matrix computation” by Gene Golub, 3.5.1 (page 124), there are two heuristics for Gaussian elimination (with partial pivoting or full pivoting). If we want to solve the problem A.x = b, that x denotes the exact solution and x0 denotes the solution found by Gaussian elimination (with partial or full pivoting), and u denotes the unit roundoff of the floating point type (u = 10^(-16) for double, u = 10^(-7) for float), we have :
Heuristic 1: Gaussian elimination produces a solution with a relatively small residual. In other words || b - A.x0 || <= u || A || || x0 ||
Heuristic 2: || x - x0 || <= u K(A) || x0 || where K(A) = || A || . || A^(-1) || is the condition number for A. Another way to state this is that if K(A) = 10^q and u = 10^(-d), Gaussian elimination produces a solution with only d - q digits of accuracy.
The reason why they are called heuristics and not theorem is because there is some kind of lie in those formula. If you do the math, one should replace u with 8.n^3.rho(A).u where rho(A) is called the growth factor. The n^3 factor is considered pessimistic in practice as floating point roundoff does not always push errors in one direction but should be considered as a “random walk”. As a consequence, it is often replaced by n. The growth factor is the component that differs with pivoting strategies:
- For partial pivoting: One can prove that rho(A) <= 2^n and that this bound is optimal as some pathological matrices can give such a growth factor.
- For full pivoting: The growth factor is lower than for partial pivoting which shows that the algorithm is “more stable”. Some bounds are known but they are not optimal. It was one thought that rho(A) <= n but it has been proven to be wrong.
In practice wether we use partial or full pivoting, the growth factor is almost never larger than 10. As a consequence, the u in the heuristics can be “safely” replaced by c(n).u where c(n) is said to be a “modest linear function” as said in the LAPACK documentation. The difference in between the 2 pivoting strategies is embedded in this “modest linear function”. The heuristic goes further and replace c(n) by 1 ! But it seems to be a good heuristic though. So yes, one can construct some pathological matrices where the growth factor for partial pivoting goes as 2^n which makes the heuristics fail for those matrices. But in practice, partial pivoting is fine. Here are two quotes on the subject :
"Is Gaussian elimination with partial pivoting stable on average? Everything we know on the subject indicates that the answer is emphatically yes, and that one needs no hypotheses beyond statistical properties to account for the success of this algorithm during nearly half a century of digital computation.” Lloyd N. Trefethen (1990)
“Intolerable pivot-growth with partial pivoting is a phenomenon that happens only to numerical analysts who are looking for that phenomenon”, William M. Kahan (1964)
As a consequence, I would suggest that both algorithms should be flagged as “stable in practice" with a focus on the 2 heuristics as they are both backward stable in practice (heuristic 1) and are forward stable (heuristic 2) only for low condition numbers. In practice full pivoting should be recommended only for rank revealing decomposition where partial pivoting is known to be bad.
Just to be clear, I am not a numerical analyst and all those claims come from my reading of the chapters from the books given below that I have read carefully.
Best regards,
François Fayard
=====
Books:
"Matrix computation”, Gene Golub, Chapter 3
“Accuracy and stability of numerical algorithms”, Nicholas Higham, Chapter 9
On the web