Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0

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

2010/7/17 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:
> 2010/7/17 Aron Ahmadia <aja2111@xxxxxxxxxxxx>:
>> Hi List,
>>
>> Just responding to this discussion on condition number functions and
>> their availability in similar packages to Eigen..
>>
>> LAPACK's expert driver routines do calculate the condition number of
>> the matrix (http://www.netlib.org/lapack/single/).
>
> Sorry I shouldn't have written that "LAPACK offers any condition
> number", instead I should have written that LAPACK doesn't standardize
> on a single notion of condition number.
>
> LAPACK does exposes the condition number of certain decompositions in
> certain expert solvers. However, these condition numbers are mutually
> inequivalent in different expert solvers, and are not the notion of
> condition number that you describe below.
>
> For example, the first link on that page mentioning condition numbers:
>
> http://www.netlib.org/lapack/single/sgesvx.f
>
> is actually doing "LU condition numbers". Namely, it is calling SGECON,
>
> http://www.netlib.org/lapack/single/sgecon.f
>
> which is doing this:
>
> *  An estimate is obtained for norm(inv(A)), and the reciprocal of the
> *  condition number is computed as
> *     RCOND = 1 / ( norm(A) * norm(inv(A)) ).
> *
> *  Arguments
> *  =========
> *
> *  NORM    (input) CHARACTER*1
> *          Specifies whether the 1-norm condition number or the
> *          infinity-norm condition number is required:
> *          = '1' or 'O':  1-norm;
> *          = 'I':         Infinity-norm.
>
> This is not the same as the condition number computed from the
> quotient of singular values (that would be the above formula with the
> spectral norm instead of L1 or Linf norm, which can be very different)
>
>>  In Numerical
>> Linear Algebra, I would consider the condition number's definition
>> under the 2-norm as the largest singular value divided by the smallest
>> to be well-accepted (see Trefethen and Bau, "Numerical Linear
>> Algebra", Strang, "Introduction to Linear Algebra")
>
> From the above, we already see that it's not "well accepted" by
> LAPACK! I'm not even sure if any LAPACK function is using this notion
> of condition number.

To be clear: I do agree that the notion of condition number that you
mention here is the only "geometrically" meaningful one in general,
and it is the same as I had mentioned in my previous email; however,
we see that LAPACK uses decomposition-specific notions of condition
numbers instead of this one. Which makes sense because 1) this notion
of condition number requires computing a SVD in general and 2) for
many decompositions, other more decomposition-specific notions are
more relevant.


>
> Here is yet another LAPACK function doing totally different condition
> numbers than what we saw above:
>
> http://www.netlib.org/lapack/single/sgeesx.f
>
> quote:
>
> *  SENSE   (input) CHARACTER*1
> *          Determines which reciprocal condition numbers are computed.
> *          = 'N': None are computed;
> *          = 'E': Computed for average of selected eigenvalues only;
> *          = 'V': Computed for selected right invariant subspace only;
> *          = 'B': Computed for both.
> *          If SENSE = 'E', 'V' or 'B', SORT must equal 'S'.
>
>>
>> http://books.google.com/books?id=Gv4pCVyoUVYC&pg=PA462&lpg=PA462&dq=strang+condition+number&source=bl&ots=SVf4gZC6WK&sig=gXUXnq2IossfbUyIXZ7ogDrDm8I&hl=en&ei=ucdBTJOoOZWC4Qbw-5icDg&sa=X&oi=book_result&ct=result&resnum=5&ved=0CCgQ6AEwBA#v=onepage&q=condition%20number&f=false
>>
>> It is also offered in
>>
>> NumPy
>> http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.cond.html
>>
>> MATLAB
>> http://www.mathworks.com/access/helpdesk/help/techdoc/ref/cond.html
>
> Hm, these aren't very big authorities in numerical linear algebra,
> compared to LAPACK...
>
>> I know that condition number refers to problems, not specific
>> operators, however the common parlance, at least in American textbooks
>> and numerical computing software, is to treat "condition number", as
>> an alias for "condition number of matrix inversion".
>
> The problem is that even "condition number of matrix inversion" is
> very vague, because there are different ways to invert a matrix, and
> depending one the one that you use, different notions of condition
> number are relevant. For example, if you invert a matrix using LU,
> what matters is to avoid having near-zero coefficients on the diagonal
> of U. If you invert a matrix using SVD, what matters is the ratio of
> max/min singular values. If you invert using cofactors, what matters
> is the ratio between the determinnat and the cofactors.
>
> Cheers,
> Benoit
>
>>
>> Hope this helps,
>> Aron
>>
>> On Sat, Jul 17, 2010 at 6:01 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:
>>> 2010/7/17 Robert Lupton the Good <rhl@xxxxxxxxxxxxxxxxxxx>:
>>>> Thanks for the quick response and analysis.
>>>>
>>>>> Your matrix is indeed not very well conditioned for floats. Are you
>>>>> using fixed size matrices (Matrix4f) or dynamically sized ones
>>>>> (MatrixXf) ? I'm asking because there is a special path to invert 4x4
>>>>> fixed size matrices, and as far as I remember there was bug for some
>>>>> special configurations in the early versions of Eigen 2.0, so I would
>>>>> recommend you to upgrade it to the latest.
>>>>
>>>> Fixed size, Matrix4f (I thought I included that in the post; sorry).  Yes, I should upgrade but that's a little inconvenient in quite a large software system.  Doable, and will be done!
>>>>
>>>>>
>>>>> determinant != condition:
>>>>> http://en.wikipedia.org/wiki/Condition_number
>>>>>
>>>>> Matlab says the condition of your matrix is 3.8e8, which is pretty bad. (btw: is there a way to calculate condition with Eigen?)
>>>>
>>>> I know that I wanted a condition number, and quoted the determinant as I couldn't find any eigen way of finding either the condition number, or whether it thinks that the inversion succeeded.
>>>
>>> Eigen doesn't offer any "conditionNumber()" function because there is
>>> no such thing as the "condition number" of a matrix in general. There
>>> are many mutually inequivalent things that are called "the condition
>>> number of a matrix" in different contexts.
>>>
>>> The most general notion of condition number, the only one to make
>>> sense for all matrices, is as the quotient of the biggest singular
>>> value over the smallest singular value.
>>>
>>> However, in the setting of LU/LLT decomposition and inverting
>>> matrices, this is NOT what is usually called the "condition number".
>>> Instead, in the context of LU, one usually calls "condition number"
>>> the quotient of the biggest absolute value of a diagonal coeff of U,
>>> over the smallest. Thus this "condition number" doesn't have a simple
>>> meaning in terms of the original matrix, and is only interesting to
>>> know in the context of LU (it indicates how instable LU decomposition
>>> is, without full pivoting).
>>>
>>> By the way, I don't think that LAPACK offers any "condition number"
>>> computation either.
>>>
>>> In short, it seems to be a case where it's best to let the user
>>> compute his own "condition number" for himself, which is easy from the
>>> decompositions, as that ensures that he understands what he's
>>> computing.
>>>
>>> All that being said, a condition number of 3e+8 for a float matrix
>>> (so, bigger than 1/epsilon) really means that you shouldn't try
>>> inverting it; it's true for both notions of condition number discussed
>>> above.
>>>
>>> Benoit
>>>
>>>>
>>>>                                        R
>>>>
>>>>
>>>>
>>>>
>>>
>>>
>>>
>>
>>
>>
>