Re: [eigen] Documentation: Partial vs Full pivot in LU decomposition |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen <eigen@xxxxxxxxxxxxxxxxxxx>*Subject*: Re: [eigen] Documentation: Partial vs Full pivot in LU decomposition*From*: Gael Guennebaud <gael.guennebaud@xxxxxxxxx>*Date*: Mon, 21 Jan 2019 15:28:44 +0100*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=20161025; h=mime-version:references:in-reply-to:from:date:message-id:subject:to; bh=fBgYU4Ah60ymMWW9swb9WaqjhNnOdAYdnisFUJvu1sk=; b=XXi4wvf9OmEHX1YBIF6uaxy0bxQVKpzYcs7cN1hqh9wIXk/VdNHvDweUYG1r0xxdIT qxS92dmRm4B3bt52bda1jccwwMbmaBhF/AlGGupyEP3v2UxhqBjBsVG9ROfC0bQThRAL I91BFiDMB4vabHv2hxULeIuupU6TH2PQA/mOs+Qhrovfe51wQqGg+bE3kN5/QhaTmVoO ioxnNEA/35TKswkRHeHPAQJAJREbv10Lp22ydX8uJkfnH42e1yoEdo9+eyLo7lToQiII pSJIA4Xa0kZ+Rk3KXJHKWY1DC+NBbpEmF4mBVqQmPpv/tu/W0ukFCe4f+5XMQVzjBI3A 49vg==

Hi,

I too find this column "Algorithm reliability and accuracy" somewhat subjective and neither "proven" nor "depends on conditional number" are really satisfactory. I would not put them on equal foots though. Depending on your context PartialPivLU might either turn out to be 100% successful or always failing. This is why people investigated efficient alternative to full-piv-LU like Panel-Rank-Revealing and rook pivoting. This distinction is, however, already clear from the "matrix requirements" and "Rank-revealing" columns. I wonder whether it would not be better to simply drop it in favor of more useful information, like least-squares or minimal-norm ability, and/or recommended usages.

Since we're at it, I would also drop the "Maturity of Eigen's implementation" column which is very subjective too, and kind of meaningless. On the other hand, we could add a column to point to the LAPACK equivalent.

This table is also missing CompleteOrthogonalDecomposition.

gael

On Fri, Jan 18, 2019 at 1:57 PM Francois Fayard <fayard@xxxxxxxxxxxxx> wrote:

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 9On the webGrowth factor for full pivoting strategy: http://math.mit.edu/~edelman/publications/complete_pivoting.pdfCleve Moler blog on the subject: https://blogs.mathworks.com/cleve/2014/08/04/gaussian-elimination-with-partial-pivoting/Average case stability of Gaussian elimination: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.157.8369&rep=rep1&type=pdf

**References**:**[eigen] Documentation: Partial vs Full pivot in LU decomposition***From:*Francois Fayard

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] C++11 Style Initialisation** - Next by Date:
**[eigen] Why Eigen::PartialPivLU::setThreshold() and Eigen::HouseholderQR::setThreshold() does not exist ?** - Previous by thread:
**[eigen] Documentation: Partial vs Full pivot in LU decomposition** - Next by thread:
**[eigen] Why Eigen::PartialPivLU::setThreshold() and Eigen::HouseholderQR::setThreshold() does not exist ?**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |