Re: [eigen] Contributing condition number estimator

[ Thread Index | Date Index | More Archives ]


Thanks for the quick response.

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


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


Mail converted by MHonArc 2.6.19+