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

[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 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
>> >
>>
>>
>
>