Re: [eigen] Contributing condition number estimator

Hi,

On Tue, Nov 24, 2015 at 1:22 AM, Rasmus Larsen <rmlarsen@xxxxxxxxx> wrote:

I would like to contribute a feature that I find lacking in Eigen's linear algebra code: Efficient condition number estimation.

yes, this would be a nice addition.

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.)

there is not solve_retval pseudo expressions anymore, but see below...

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.).

I guess that the cleanest way would be to make the public API looks like:

1) b = lu.solve(x);

2) b = lu.transpose().solve(x);

3) b = lu.adjoint().solve(x);

where PartialPivLU::adjoint()/transpose() would be implemented as in DenseBase<>, i.e., transpose() would return a Transpose<PartialPivLU>.

The general idea is that a decomposition object should resemble as close as possible to a matrix but with a special storage scheme (product of factors)...

then via a few specializations of TransposeImpl<> and Solve<> we can easily forward the above calls to:

1) lu._solve_impl(x, b); // already done through the class Solve<> that generalize solve_retval

2) lu._solve_impl_transposed<false>(x, b);

3) lu._solve_impl_transposed<true>(x, b);

I can help with that internal magic, and in the meantime you can directly implement and call _solve_impl_transposed.

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

sounds good.

Gael