Re: [eigen] TriangularView::solve() interface

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

2009/8/25 Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>:
> On Mon, 24 Aug 2009, Benoit Jacob wrote:
>
>>> Something to perhaps add to the to-do list: a routine to estimate the
>>> condition number of a matrix based on the LU decomposition.
>>
>> Is that possible at all? I assume that by condition number we mean the
>> ratio between the biggest and smallest eigenvalue.
>
> Yes, after replacing "eigenvalue" with "singular value".
>
>> I don't think that the LU decomposition sees that! The SVD, Schur, and
>> diagonalizations are the ones that can be used here.
>
> I hadn't look into it, but matlab prints a warning when you solve a
> ill-conditioned linear system, and the warning contains an estimate for the
> condition number.
>
> ---------------------- [ matlab ] -------------------
>>>
>>> A = hilb(12);
>>> b = randn(12,1);
>>> A \ b
>
> Warning: Matrix is close to singular or badly scaled.
>         Results may be inaccurate. RCOND = 2.458252e-17.
> ans =
>   1.0e+15 *
>    0.0000
>   -0.0000
>    0.0015
>      etc
> ---------------------- [ matlab ] -------------------
>
> I found this quite useful on occasions since it often points to errors in my
> code.
>
> Reading a bit more, RCOND stands for the reciprocal of the condition number
> based on the 1-norm: RCOND = 1 / ( \| A^{-1} \|_1  \| A \|_1 ). Matlab uses
> an algorithm by Nick Higham [1] which is implemented in SGECON/DGECON in
> Lapack. This algorithm computes an estimate of the condition number, and not
> the exact value. Contrary to what I thought, it does not use the LU
> decomposition but solely matrix-vector products (usually four or five).
>
>  [1] http://doi.acm.org/10.1145/50063.214386
>
> An alternative possibility, which does use the LU decomposition, is
> presented in Golub & Van Loan, Section 3.5.4. The references list even more
> methods.

I had a look at this section, but again their condition number isn't
based on eigenvalues/singular values, it is based on the infinity-norm
(biggest coefficient). OK, i didn't remember but it's true that there
are many inequivalent things called "condition number". If there are
use cases where this is one what wants, I can believe that, then it's
OK to include in Eigen if you write the patch.

The documentation just needs a good explanation that there are many
inequivalent things called "condition number", so it must explain what
is called condition number here.

Benoit