Re: [eigen] Re: news from the LU / rank / solving API front

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

2009/10/26 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:
> 2009/10/26 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
>>
>>
>> On Mon, Oct 26, 2009 at 4:29 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx>
>> wrote:
>>>
>>> Hi again,
>>>
>>> i can finally give some more news. I've essentially done Inverse.h now
>>> (in my fork) which was the big scary part (as far as API design is
>>> concerned). i've also fixed the 4x4 inverse (Euler's 2x2 trick) so it
>>> should now be numerically stable. So I expect the rest of this API
>>> reworking to go fast now.
>>>
>>> But while i'm at it i'd like to discuss one more change, so that users
>>> of the devel branch dont have to adapt again soon.
>>>
>>> About the naming of the decomposition class, we need some homogeneity.
>>> currently, "LU" does full pivoting, on the other hand the equivanlent
>>> in QR is called "FullPivotingQR"
>>> I propose:
>>> * harmonize that: always mention the pivoting in the name
>>> * abbreviate Pivoting as Piv
>>> * keep shorthand methods such as lu(), even introduce qr() etc...
>>
>> ok for me.
>>
>>>
>>> So:
>>> LU--->FullPivLU
>>> PartialLU--->PartialPivLU
>>> in QR, replace Pivoting by Piv
>>> add fullPivLu()
>>> keep lu() as shortcut to fullPivLu()
>>> partialLU() ---> partialPivLu()
>>> add qr() as shortcut to householderQr()
>>> add rrqr() as shortcut to colPivHouseholderQr()
>>>  (note: rr for rank-revealing, that's standard)
>>
>> why does lu() would refer to full pivoting, and qr() to no pivoting ?
>> Perhaps, both should refer to their respective partial pivoting variants
>> which seem to be good default choices: good accuracy and good scalability.
>
> The idea was to default to general, full-featured decompositions.
> Partial-pivoting QR is numerically stable

i mean numerically stable regardless of the condition number, e.g.
doesnt have to assume that the matrix is e.g. invertible

> and rank-revealing.
> Partial-pivoting LU is neither. Also partial-pivoting LU currently
> only applies to square matrices, and could only be extended to M>=N or
> M <=N depending on whether  we do row piv or col piv, but cant work
> for all rectangular sizes. As such, it is pretty much an "optimized
> special case", at the level of the algorithm itself. So i didn't want
> to default to it. However, inverse() and determinant() default to it
> because they are in the right special case (square invertible).
>
> ok?
>
> Benoit
>
>>
>> gael.
>>
>>>
>>> Opinions?
>>>
>>> Benoit
>>>
>>> 2009/10/18 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:
>>> > Hi,
>>> >
>>> > This is mostly for Gael and Hauke (if you're busy, skip to 5. below)
>>> > who I bugged recently about that, and a followup to the goog old
>>> > thread "LU precision tuning" on this list.
>>> >
>>> > It's been very very tedious but i think we're now on the good side of
>>> > the API refactoring in LU (and once I've finished LU the rest should
>>> > be the same).
>>> > See my fork http://bitbucket.org/bjacob/refactor-solve-api/
>>> >
>>> > 1) What Eigen 2.0 des.
>>> > In Eigen 2.0, we determine the rank in the LU decomposition itself. As
>>> > soon as some pivot is small enough to be considered zero (we have a
>>> > threshold for that, see "LU precision tuning" on this list), we stop
>>> > the LU decomposition. In LU::solve(), we also use that threshold to
>>> > determine whether a solution exists. If no solution exists according
>>> > to that test, then we don't try to return a solution.
>>> >
>>> > 2) What LAPACK does.
>>> > In dgesvx, there are 2 cases. If a pivot is exactly zero, so the
>>> > matrix is exactly singular, dgesvx bails out and doesn't give a
>>> > solution. This means that dgesvx can't be used on simple 2x2 examples
>>> > like the matrix (1&0\\0&0) and the right-hand-side (1\\0). If no pivot
>>> > is exactly zero, dgesvx just computes a solution without worrying
>>> > about accuracy if a pivot was near zero.
>>> >
>>> > 3) What the new Eigen (in the fork) *cough* 3.0 *cough* does.
>>> > As previously discussed, i realized that we should try as hard as
>>> > possible to get a solution, even approximate; the only case in which
>>> > we need to bail out is if continuing would generate inf/nan values
>>> > (that would be worthless to the user anyway); and the precision check
>>> > is something that the user can always do by himself afterwards. So
>>> > that means that we get much closer to LAPACK. On the other hand, I
>>> > don't see the need to restrict to nonsingular matrices like LAPACK
>>> > does. If some pivot is zero, we have an equation of the form 0*X=Y, we
>>> > can just ignore Y and set X=0, which can give a useful approximate
>>> > solution if Y was itself meant to be zero.
>>> >
>>> > 4) precision of the rank determination (problem solved, was a PEBKAC)
>>> > As you may know, in the past we've found it very difficult to have
>>> > good accuracy for rank() methods in e.g. LU. We always had a small
>>> > amount of failures in the LU test.
>>> > I just found the reason. Our LU with full pivoting algo consists in
>>> > taking the biggest element in the remaining corner and moving it to
>>> > become the next pivot. From there we assumed that the pivots were
>>> > automatically sorted in decreasing order. That was approximately true,
>>> > but not exactly true. In practice it was possible to get cases where
>>> > the (n+1)-th pivot was several times bigger than the n-th pivot. So
>>> > now we no longer make this assumption, we manually sort the pivots as
>>> > needed (not possible to do in the decomposition itself), and bingo! I
>>> > can now run "test_lu r1000" and it succeeds!
>>> >
>>> > 5) The threshold to determine when a pivot is nonzero is, by default,
>>> > the same as before (Higham formula). The user can override that by
>>> > calling usePrecision(value) and I introduced a special type Default_t
>>> > so one can come back to the default by doing useThreshold(Default). Ah
>>> > yes and nevermind what i said about virtual methods, I had been
>>> > smoking some good raccoon fur blend (a Canadian Thanksgiving
>>> > tradition).
>>> >
>>> > Benoit
>>> >
>>>
>>>
>>
>>
>