commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Luc Maisonobe <Luc.Maison...@free.fr>
Subject Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Date Thu, 08 Sep 2011 19:42:33 GMT
Le 08/09/2011 09:15, Phil Steitz a écrit :
> On 9/7/11 6:34 AM, Gilles Sadowski wrote:
>> Hi.
>>
>>>>>> In class "AbstractLeastSquaresOptimizer" (in "o.a.c.m.optimization.general"),
>>>>>> the method "getCovariances()" uses "LUDecompositionImpl" to compute
the
>>>>>> inverse of a matrix.
>>>>>> In my application, this leads to a "SingularMatrixException". If
I change
>>>>>> "LUDecompositionImpl" to "QRDecompositionImpl", no exception is raised.
>>>>>> Also, keeping "LUDecompositionImpl" but passing a much lower singularity
>>>>>> threshold, does not raise the exception either.
>>>>>>
>>>>>> Thus, I wonder whether there was a reason for using "LU", and if
not,
>>>>>> whether I could change the decomposition solver to "QR" (as this
is a
>>>>>> cleaner solution than guessing a good value for the threshold).
>>>>> There are no reason for LU decomposition, and QR decomposition is
>>>>> known to be more stable. So I would also consider switching to this
>>>>> algorithm is a cleaner solution.
>>>> Fine. I'll open a JIRA issue.
>>>>
>>>> A unit test "testNonInvertible" in "LevenbergMarquardtOptimizerTest" fails
>>>> with the change to "QRDecomposition" because no "SingularMatrixException"
>>>> is raised anymore.
>>>> What was the purpose of that test?
>>> The purpose was to check that impossible problems were detected properly.
>> My question should have been clearer: Was the previous behaviour correct
>> (i.e. an exception *must* be raised somehow)?
>
> Yes.  The getCovariances method is trying to invert a singular
> matrix.  SingularMatrixException is appropriate there.  The problem
> is that the QR decomp method has no threshold defined for
> singularity, doing an exact check against 0 on the elements of
> rDiag.  Many (most?) singular matrices or
> near-singular-enough-to-lead-to-meaningless-results matrices will
> not lead to exact 0's there due to rounding in the computations.
> That's why the LU impl defines a threshold.   For the (exactly)
> singular matrix in the test case above, LU does not get exact 0's
> either on pivots, but the default singularity threshold (correctly)
> detects the singularity.  IMO, QR should do the same thing;
> otherwise we will happily return meaningless results, as in this
> case.  In this case, there is no question, since the matrix is
> exactly singular.  In cases where matrices are near-singular, we are
> doing users a favor by throwing, since inversion results are going
> to be unstable.
>
> A possibly more robust option here is to use Cholesky decomposition,
> which is known to be stable for symmetric positive definite
> matrices, which the covariance matrix being inverted here should
> be.  The exceptions thrown will be different; but they will give
> more specific information about what is wrong with the covariance
> matrix.
>
> Luc - are there other reasons that QR would be better for cov
> matrices?   I would have to play with a bunch of examples, but I
> suspect with the defaults, Cholesky may do the best job detecting
> singular problems.

I'm not sure about Cholesky, but I have always thought that at least QR 
was better than LU for near singular matrices, with only a factor 2 
overhead in number of operations (but number of operations is not the 
main bottleneck in modern computers, cache behavior is more important).

Luc

>
> Phil
>> The switch to "QR" seems to imply that a previously impossible problem is
>> now quite possible.  So, is the problem really impossible or was the test
>> actually testing a fragile implementation of "getCovariances()"?
>>
>>
>> Regards,
>> Gilles
>>
>> ---------------------------------------------------------------------
>> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
>> For additional commands, e-mail: dev-help@commons.apache.org
>>
>>
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
>
>


---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org


Mime
View raw message