Re: [eigen] Contributing condition number estimator

Hi Gael,

I realized that I did not export this from Mercurial in the recommended way. Attached is the same patch exported with "hg export top > myfile". Would you have time to add the template machinery around this to make expressions like

1) b = lu.solve(x);

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

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

work? It would be preferable if .transpose() and .adjoint() is also defined for symmetric factorizations (LLT, LDLT) as a no-op / identity operator. If you don't have time, could you give me a few pointers on what to modify?

Thanks,

Rasmus

On Wed, Nov 25, 2015 at 11:52 AM, Rasmus Larsen <rmlarsen@xxxxxxxxxx> wrote:

Hi Gael,

Please find attached a patch to add _solve_impl_transposed to the two LU classes, and associated additions to existing unit tests.
This is with respect to the default branch in the Mercurial repository.

Best,
Rasmus

On Tue, Nov 24, 2015 at 8:37 AM, Rasmus Larsen <rmlarsen@xxxxxxxxxx> wrote:
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