Re: rcond, condest, and a block 1-norm estimator

Jason Riedy wrote:
> And David Bateman writes:
>   
>> I believe the rcond function should take rcond from the LAPACK
>> factorization, at least I believe is what matlab does. [...]
>>     
>
> For various reasons, I'm not referring to MATLAB(tm) itself or
> any docs other than what is on-line and trivially findable.  I
> have no idea if MATLAB's rcond applies to sparse matrices.  I
> know LAPACK's does not.
>   
>> help rcond
 RCOND  LAPACK reciprocal condition estimator.
    RCOND(X) is an estimate for the reciprocal of the
    condition of X in the 1-norm obtained by the LAPACK
    condition estimator. If X is well conditioned, RCOND(X)
    is near 1.0. If X is badly conditioned, RCOND(X) is
    near EPS.
 
    See also COND, NORM, CONDEST, NORMEST.

>> rcond(speye(1024))
??? Error using ==> rcond
Use rcond(full(S)) or condest(S).

So the matlab rcond function works only for full matrices and uses the
lapack rcond estimation from the factorization.. There needs to be a
little reorganization to get this to work optimally in Octave. Till then
a minimal rcond implementation for Octave is

function  rc = rcond (x)
  [dummy, rc] = inv (x);
endfunction

As Octave also calculates the condition number for inverse of sparse
matrices if requested, this works also for sparse matrices, though
forming the inverse of a sparse matrix is always a bad idea..

> That said, I have been losing the argument to use the block
> estimator in LAPACK...  We may weaken the routine slightly
> instead.  For "small" problems, computing RCOND is very, very
> expensive right now.  The current norm estimator uses a very safe
> but very slow triangular solve.  We likely will replace that with
> the fast solve, followed by a check for overflow and re-running
> with the slow solve if necessary.  That will catch and re-scale
> overflows, but it won't catch possibly avoidable underflows
> during estimation.
>
> So it's up to y'all what rcond() should be.  The routine I
> included computes the inverse of the 1-norm condition number of A
> (iirc), and LAPACK's RCOND returned by xyySVX is the inverse of
> the inf-norm condition number of the matrix that actually is
> factored.  That matrix may have been equilibrated, so determining
> how LAPACK's RCOND relates to your original matrix requires
> examining a bunch of parameters MATLAB's rcond() doesn't return.
> Without experimenting, I can't tell what MATLAB's rcond() really
> does.
>   
I believe its simpler than that.. I think matlab uses the xGETRF
functions of lapack that return an estimate of rcond, and matlab just
returns that estimate..
> The only use for LAPACK's RCOND directly is to see if the
> factorization algorithm came close to solving with a nearly a
> singular matrix when forming a Schur complement.  It tells you
> very little about the actual problem you are trying to solve or
> even about any errors when using the factorization; RCOND may be
> measured with an entirely different norm (scaled by
> equilibration) than your problem.
>   
xGETRF scales the inverse norm estimation by the 1-norm and so I don't
think this is the case. You're certainly right for the xGESVX LAPACK
routines though. The result is that if the matrix is near singular then
the rcond extimate will be incorrect. However, as rcond is only supposed
to be used to estimate if a matrix is singular, then this error probably
isn't an issue.

> And if you pick up our extra-precise refinement drivers (once we
> actually ship the damned things, see LAWN 165), you will have
> more and more relevant condition numbers from which to chose.
>   
Choice is good for a developer, though how to hide it from the user
behind simple driver routines is always an issue.

>   
>> Matlab calls the normest1 function, that is basically your
>> block_onenorm_est, and then multiples the estimated inverse
>> norm by norm(A,1) to get the condition number estimate.
>>     
>
> >From MathWorks' pages, it seems normest has a very different
> interface.  Hence the cumbersome block_onenorm_est name.
>   
normest is a 2-norm estimate using the power method to find the largest
singular value. Octave has that routine. No your block_onenorm_est is
like the normest1 routine from Matlab. The interface to normest1 has a
function handle to the solver passed as the first argument, whereas your
code has a solver and a solver for the conjugate transpose. Frankly I
don't think normest1 in itself is a useful function and so if its not
implemented who cares. We just need a matlab compatible interface to
condest..

> But your description is not the RCOND returned from LAPACK's
> linear system drivers.  ;)
>   
No, I was trying in my non-mathematicians manner to describe your code :-)

>   
>> condest is one of the annoying missing sparse functions in
>> Octave...
>>     
>
> Amen.  Limits Octave's applicability for many colleagues.
>   

Well, I definitely want condest and your code gives it. I see if I can
convert it to Octave style and sent it back to you for your comments..

>   
>> I'd hesitate to include it though with the rcond function with
>> its current name...
>>     
>
> I have little ego about these things.  Use, change, and rename as
> you see fit.  Just thought I should send along my good-enough
> versions as a starting point.
>
> I'm not in a position to get into the copyright assignment
> argument here right now, so feel free to make sufficient,
> copyrightable changes that are assigned and licensed
> appropriately.  ;) These are just the routines I'm using right
> now.
>   
If I used the BSD license you've chosen it allows relicensing under a
GPL, so I would suggest doing that and keeping you as the author..

Cheers, and thanks for the code.
D.




> Jason
>
>   


-- 
David Bateman                                [email protected]
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary