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

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

All in all, it's more complicated than I thought ...

Cheers,
Jitse