[eigen] Contributing condition number estimator |

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

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: [eigen] Contributing condition number estimator*From*: Rasmus Larsen <rmlarsen@xxxxxxxxx>*Date*: Mon, 23 Nov 2015 16:22:57 -0800*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=20120113; h=mime-version:date:message-id:subject:from:to:content-type; bh=QQqwtq8w9tl9qZglkrAuyhY5yloIm1wBzPp2pyCwYpM=; b=VGDGSHTwgWDagIbuniqnsgnSnosszewXqGVf29dwySjYdsiMojlQVMWyObsLztADQA osQuU2IWYYrbQ+iYSc8TT583C1XzqgxdFfIyXizvUW8FVCwqtEqAzHTP2zqSVidwLKOW KX+3P2S4OScRlpMNZAKNA9L92mix9z4lI4MVL3WbYGmzVwCIxAYPbNNp+O7BjcLPpisO 3VOCWxbs5/8EnpHDYmK0me8PVPKYG4KEgSnlaZ66M/XqyC0Th/g062r9jvb3XcJjQvn6 Brtup8m+gQsbZDnM38PHgMrtF3UcPzBFKXf+acQid+EXiOU2cU2q/at2mnxVlilVKk1T EAbg==

Hi Gael & Eigeners,

I would like to contribute a feature that I find lacking in Eigen's linear algebra code: Efficient condition number estimation. I have implemented Higham's improved version of Hager's L1 norm estimation algorithm described in:

N. J. Higham, "FORTRAN Codes for Estimating the One-Norm of a Real or

Complex Matrix with Applications to Condition Estimation", ACM

Transactions on Mathematical Software, Vol. 14, No. 4, December 1988,

Pages 381-396.

The algorithm performs at most 10 matrix-vector multiplies (often much less). It's main application is for condition number estimation and it is used by the xyyCON condition estimation routines in LAPACK. The motivation is that the best trade-off between speed and reliability when solving general linear systems is often to use LU with partial pivoting (which has a fast blocked implementation in Eigen) and then do this fast O(n^2) condition estimation to verify that it is "safe" to use. The L1-estimator can also be used to compute accurate forward and backward error estimates for the solutions.

My implementation uses Eigen but resides in a separate code base. I am starting to move it into Eigen itself, but rather than just submitting a patch, I would like to discuss what would be the cleanest way.

The first step that I would like you opinion on, is to how to add functionality to the non-symmetric decompositions (specifically FullPivLU and PartialPivLU) for solving the system adjoint system A^* x = rhs. Here I see two approaches:

a) Add new methods "solveAdjoint" in the LU classes that return a new type solve_adjoint_retval internal type in Solve.h. (I have this implemented already. Benoit Jacob reviewed it.)

b) Add a template parameter "bool Adjoint" (with default value false) to existing solve methods in the LU classes (and pass it to solve_retval etc.).

The second step would be to add a method "Scalar rcond()" to all the decomposition classes returning the estimate of the reciprocal condition number.

Let me know what you think/prefer.

Thanks,

Rasmus

**Follow-Ups**:**Re: [eigen] Contributing condition number estimator***From:*Gael Guennebaud

**Messages sorted by:**[ date | thread ]- Prev by Date:
**[eigen] Google TensorFlow uses Eigen** - Next by Date:
**Re: [eigen] Contributing condition number estimator** - Previous by thread:
**[eigen] Google TensorFlow uses Eigen** - Next by thread:
**Re: [eigen] Contributing condition number estimator**

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