Re: [eigen] Contributing condition number estimator

Hi,

Thanks for the quick response.

On Tue, Nov 24, 2015 at 6:04 AM, Gael Guennebaud <gael.guennebaud@xxxxxxxxx> wrote:

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

That would indeed be a more elegant solution.

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.

Sounds good. We are still on 3.2.7 at Google, so I will have to prepare separate patches. But as you say, I can call something like _solve_impl_transposed directly in my 3.2 code.

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