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

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

*To*: eigen <eigen@xxxxxxxxxxxxxxxxxxx>*Subject*: [eigen] news from the LU / rank / solving API front*From*: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>*Date*: Sun, 18 Oct 2009 01:14:48 -0400*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:date:message-id:subject :from:to:content-type; bh=SM4IT62BNi27VhkfaENFwze1kkGUUYrI6TanmR5yZJo=; b=h4WvPuwb1jk3ReLxJlILS1loBpnoHU76xybQ2/CK26713rGHDUBhTNxyC2SQHrN2ii t9qLEjME6bmLOErCmyfwk4Pin8baedji5IAHptWfm9MZrxfjl9HVYYSHWszPnKVedAh6 agJaEmgwVgAFv8gmt+Tpw+l98IC9bjA8OfOFo=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:date:message-id:subject:from:to:content-type; b=Cw9+SesXHUM3f/CVJH/VNZZe2MpG/Yu8L1/ag2X5IUXvB1wZApm5umvV5A1iL409m5 i7Wu48GxVFU1RS9LmtBnM0sq3jlL4TaM/Se8TqsBwoXUvwF2Ccb1xT1qE+zvvdWvEXtS mnKPdgMSB/zrZNgBLmI2AV9rsedPBdF40KGfM=

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

**Follow-Ups**:**[eigen] Re: news from the LU / rank / solving API front***From:*Benoit Jacob

**Re: [eigen] news from the LU / rank / solving API front***From:*Hauke Heibel

**[eigen] Re: news from the LU / rank / solving API front***From:*Benoit Jacob

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] Re: Debugging in VS.net (autoexp.dat)** - Next by Date:
**[eigen] Re: news from the LU / rank / solving API front** - Previous by thread:
**Re: [eigen] AutoDiffScalar** - Next by thread:
**[eigen] Re: news from the LU / rank / solving API front**

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