You are viewing a plain text version of this content. The canonical link for it is here.
Posted to dev@commons.apache.org by Gilles Sadowski <gi...@harfang.homelinux.org> on 2011/09/07 12:45:43 UTC
[Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Hello.
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).
Regards,
Gilles
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Greg Sterijevski <gs...@gmail.com>.
If speed is not a factor, why not do the best thing and calculate the
eigenvalue decomposition and then form a moore-penrose inverse. The
inversion will never fail (unless all eigenvalues are zero). You will also
have an immediate diagnostic of the stability of the system (through the
condition or inverse condition number).
On Thu, Sep 8, 2011 at 2:41 AM, Phil Steitz <ph...@gmail.com> wrote:
> On 9/7/11 9:26 PM, Greg Sterijevski wrote:
> > I thought that QR is of O(n^3) complexity, while LU is probably in the
> > vicinity of O(n^2.5). While this is not a barn burning improvement, it is
> an
> > improvement. Consider that if you are using the covariance of the
> parameters
> > in deciding how to make the next step (some kind of inverse hessian
> rule),
> > the repeated calls to inverse add up?
>
> Sounds a little contrived. I can't see the speed difference as
> really material given that we are talking about parameter covariance
> matrices. I would be fine with any of the three - LU, Cholesky
> (which may also be fastest, btw) or QR.
>
> Phil
> >
> > Would it be a terrific eyesore to include a constructor with a boolean
> > argument that says bUseQR? The call to LU is localized, so there is
> minimal
> > refactoring. One could also pass the boolean into the getCovariances
> > method.
> >
> >
> >
> >
> > On Wed, Sep 7, 2011 at 4:48 PM, Ted Dunning <te...@gmail.com>
> wrote:
> >
> >> It shouldn't be all that different from QR (at most about 2x different).
> >>
> >> On Wed, Sep 7, 2011 at 1:16 PM, Greg Sterijevski <
> gsterijevski@gmail.com
> >>> wrote:
> >>> It is also my recollection that LU is very quick to calculate. Would it
> >> be
> >>> possible to allow users to choose?
> >>>
> >>> On Wed, Sep 7, 2011 at 3:07 PM, Ted Dunning <te...@gmail.com>
> >> wrote:
> >>>> Does the LUDecomposition not use pivots? LU should always do so since
> >> it
> >>>> is
> >>>> numerically unstable otherwise. I would be surprised if it doesn't
> >> given
> >>>> the normal level of quality in commons math.
> >>>>
> >>>> QR is, of course, almost always preferable to LU as you note. But I
> >>> would
> >>>> be surprised at radically different answers.
> >>>>
> >>>> Perhaps the only real difference between the two methods in this one
> >> case
> >>>> is
> >>>> a difference in threshold.
> >>>>
> >>>> What is the condition number of your matrix?
> >>>>
> >>>> On Wed, Sep 7, 2011 at 6:34 AM, Gilles Sadowski <
> >>>> gilles@harfang.homelinux.org> 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)?
> >>>>> 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
>
>
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Phil Steitz <ph...@gmail.com>.
On 9/7/11 9:26 PM, Greg Sterijevski wrote:
> I thought that QR is of O(n^3) complexity, while LU is probably in the
> vicinity of O(n^2.5). While this is not a barn burning improvement, it is an
> improvement. Consider that if you are using the covariance of the parameters
> in deciding how to make the next step (some kind of inverse hessian rule),
> the repeated calls to inverse add up?
Sounds a little contrived. I can't see the speed difference as
really material given that we are talking about parameter covariance
matrices. I would be fine with any of the three - LU, Cholesky
(which may also be fastest, btw) or QR.
Phil
>
> Would it be a terrific eyesore to include a constructor with a boolean
> argument that says bUseQR? The call to LU is localized, so there is minimal
> refactoring. One could also pass the boolean into the getCovariances
> method.
>
>
>
>
> On Wed, Sep 7, 2011 at 4:48 PM, Ted Dunning <te...@gmail.com> wrote:
>
>> It shouldn't be all that different from QR (at most about 2x different).
>>
>> On Wed, Sep 7, 2011 at 1:16 PM, Greg Sterijevski <gsterijevski@gmail.com
>>> wrote:
>>> It is also my recollection that LU is very quick to calculate. Would it
>> be
>>> possible to allow users to choose?
>>>
>>> On Wed, Sep 7, 2011 at 3:07 PM, Ted Dunning <te...@gmail.com>
>> wrote:
>>>> Does the LUDecomposition not use pivots? LU should always do so since
>> it
>>>> is
>>>> numerically unstable otherwise. I would be surprised if it doesn't
>> given
>>>> the normal level of quality in commons math.
>>>>
>>>> QR is, of course, almost always preferable to LU as you note. But I
>>> would
>>>> be surprised at radically different answers.
>>>>
>>>> Perhaps the only real difference between the two methods in this one
>> case
>>>> is
>>>> a difference in threshold.
>>>>
>>>> What is the condition number of your matrix?
>>>>
>>>> On Wed, Sep 7, 2011 at 6:34 AM, Gilles Sadowski <
>>>> gilles@harfang.homelinux.org> 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)?
>>>>> 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
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Ted Dunning <te...@gmail.com>.
LU is theoretically lower if you have a faster matrix multiply.
Otherwise it is about 2n^3/3 + n^2 + n/3, I think.
On Wed, Sep 7, 2011 at 9:26 PM, Greg Sterijevski <gs...@gmail.com>wrote:
> I thought that QR is of O(n^3) complexity, while LU is probably in the
> vicinity of O(n^2.5). While this is not a barn burning improvement, it is
> an
> improvement. Consider that if you are using the covariance of the
> parameters
> in deciding how to make the next step (some kind of inverse hessian rule),
> the repeated calls to inverse add up?
>
> Would it be a terrific eyesore to include a constructor with a boolean
> argument that says bUseQR? The call to LU is localized, so there is minimal
> refactoring. One could also pass the boolean into the getCovariances
> method.
>
>
>
>
> On Wed, Sep 7, 2011 at 4:48 PM, Ted Dunning <te...@gmail.com> wrote:
>
> > It shouldn't be all that different from QR (at most about 2x different).
> >
> > On Wed, Sep 7, 2011 at 1:16 PM, Greg Sterijevski <gsterijevski@gmail.com
> > >wrote:
> >
> > > It is also my recollection that LU is very quick to calculate. Would it
> > be
> > > possible to allow users to choose?
> > >
> > > On Wed, Sep 7, 2011 at 3:07 PM, Ted Dunning <te...@gmail.com>
> > wrote:
> > >
> > > > Does the LUDecomposition not use pivots? LU should always do so
> since
> > it
> > > > is
> > > > numerically unstable otherwise. I would be surprised if it doesn't
> > given
> > > > the normal level of quality in commons math.
> > > >
> > > > QR is, of course, almost always preferable to LU as you note. But I
> > > would
> > > > be surprised at radically different answers.
> > > >
> > > > Perhaps the only real difference between the two methods in this one
> > case
> > > > is
> > > > a difference in threshold.
> > > >
> > > > What is the condition number of your matrix?
> > > >
> > > > On Wed, Sep 7, 2011 at 6:34 AM, Gilles Sadowski <
> > > > gilles@harfang.homelinux.org> 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)?
> > > > > 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
> > > > >
> > > > >
> > > >
> > >
> >
>
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Greg Sterijevski <gs...@gmail.com>.
I thought that QR is of O(n^3) complexity, while LU is probably in the
vicinity of O(n^2.5). While this is not a barn burning improvement, it is an
improvement. Consider that if you are using the covariance of the parameters
in deciding how to make the next step (some kind of inverse hessian rule),
the repeated calls to inverse add up?
Would it be a terrific eyesore to include a constructor with a boolean
argument that says bUseQR? The call to LU is localized, so there is minimal
refactoring. One could also pass the boolean into the getCovariances
method.
On Wed, Sep 7, 2011 at 4:48 PM, Ted Dunning <te...@gmail.com> wrote:
> It shouldn't be all that different from QR (at most about 2x different).
>
> On Wed, Sep 7, 2011 at 1:16 PM, Greg Sterijevski <gsterijevski@gmail.com
> >wrote:
>
> > It is also my recollection that LU is very quick to calculate. Would it
> be
> > possible to allow users to choose?
> >
> > On Wed, Sep 7, 2011 at 3:07 PM, Ted Dunning <te...@gmail.com>
> wrote:
> >
> > > Does the LUDecomposition not use pivots? LU should always do so since
> it
> > > is
> > > numerically unstable otherwise. I would be surprised if it doesn't
> given
> > > the normal level of quality in commons math.
> > >
> > > QR is, of course, almost always preferable to LU as you note. But I
> > would
> > > be surprised at radically different answers.
> > >
> > > Perhaps the only real difference between the two methods in this one
> case
> > > is
> > > a difference in threshold.
> > >
> > > What is the condition number of your matrix?
> > >
> > > On Wed, Sep 7, 2011 at 6:34 AM, Gilles Sadowski <
> > > gilles@harfang.homelinux.org> 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)?
> > > > 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
> > > >
> > > >
> > >
> >
>
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Ted Dunning <te...@gmail.com>.
It shouldn't be all that different from QR (at most about 2x different).
On Wed, Sep 7, 2011 at 1:16 PM, Greg Sterijevski <gs...@gmail.com>wrote:
> It is also my recollection that LU is very quick to calculate. Would it be
> possible to allow users to choose?
>
> On Wed, Sep 7, 2011 at 3:07 PM, Ted Dunning <te...@gmail.com> wrote:
>
> > Does the LUDecomposition not use pivots? LU should always do so since it
> > is
> > numerically unstable otherwise. I would be surprised if it doesn't given
> > the normal level of quality in commons math.
> >
> > QR is, of course, almost always preferable to LU as you note. But I
> would
> > be surprised at radically different answers.
> >
> > Perhaps the only real difference between the two methods in this one case
> > is
> > a difference in threshold.
> >
> > What is the condition number of your matrix?
> >
> > On Wed, Sep 7, 2011 at 6:34 AM, Gilles Sadowski <
> > gilles@harfang.homelinux.org> 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)?
> > > 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
> > >
> > >
> >
>
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Greg Sterijevski <gs...@gmail.com>.
It is also my recollection that LU is very quick to calculate. Would it be
possible to allow users to choose?
On Wed, Sep 7, 2011 at 3:07 PM, Ted Dunning <te...@gmail.com> wrote:
> Does the LUDecomposition not use pivots? LU should always do so since it
> is
> numerically unstable otherwise. I would be surprised if it doesn't given
> the normal level of quality in commons math.
>
> QR is, of course, almost always preferable to LU as you note. But I would
> be surprised at radically different answers.
>
> Perhaps the only real difference between the two methods in this one case
> is
> a difference in threshold.
>
> What is the condition number of your matrix?
>
> On Wed, Sep 7, 2011 at 6:34 AM, Gilles Sadowski <
> gilles@harfang.homelinux.org> 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)?
> > 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
> >
> >
>
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Phil Steitz <ph...@gmail.com>.
On 9/7/11 1:07 PM, Ted Dunning wrote:
> Does the LUDecomposition not use pivots?
Yes
> LU should always do so since it is
> numerically unstable otherwise. I would be surprised if it doesn't given
> the normal level of quality in commons math.
>
> QR is, of course, almost always preferable to LU as you note. But I would
> be surprised at radically different answers.
Me to.
>
> Perhaps the only real difference between the two methods in this one case is
> a difference in threshold.
Yup. I think QR is effectively bugged because there is no threshold.
>
> What is the condition number of your matrix?
>
> On Wed, Sep 7, 2011 at 6:34 AM, Gilles Sadowski <
> gilles@harfang.homelinux.org> 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)?
>> 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
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Ted Dunning <te...@gmail.com>.
Does the LUDecomposition not use pivots? LU should always do so since it is
numerically unstable otherwise. I would be surprised if it doesn't given
the normal level of quality in commons math.
QR is, of course, almost always preferable to LU as you note. But I would
be surprised at radically different answers.
Perhaps the only real difference between the two methods in this one case is
a difference in threshold.
What is the condition number of your matrix?
On Wed, Sep 7, 2011 at 6:34 AM, Gilles Sadowski <
gilles@harfang.homelinux.org> 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)?
> 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
>
>
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Ted Dunning <te...@gmail.com>.
Rank revealing Cholesky should be effectively about as good as an eigen
solution. The pivoting results in non-increasing diagonals and a sharp rank
indicator. You can do pseudo inverse things with the result as well.
Getting the actual pseudo inverse isn't a bad thing. I just don't see it as
entirely necessary.
On Thu, Sep 8, 2011 at 7:49 AM, Greg Sterijevski <gs...@gmail.com>wrote:
> Cholesky, in my opinion, is not robust as you have discovered. When it
> encounters a non-psd matrix it gives up. Maybe that is the correct course
> of
> action, but I still think that when you are using the getCovariance to
> estimate the curvature in the neighborhood of a point it would be okay to
> take the generalized inverse and not worry too much about why your matrix
> is
> bordering on non-PSDness... If your optimization stops on that point, that
> is another story and should be flagged.
>
>
>
> > 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.
> >
> > I've tried it with my problem, and it also throws an exception.
> > However, I would like to obtain the covariance matrix anyway, because
> I've
> > no other clue as to what might be wrong.
> > So I think that, at least, users should be able to set the positive
> > definiteness threshold in order to avoid raising an exception.
> >
> >
> >
>
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Ted Dunning <te...@gmail.com>.
I offer the Mahout implementation from MAHOUT-792 (not yet committed) as a
reference that does pivoting and has a settable threshold. It also uses
vectorized view operations to avoid touching the internal representation if
that helps.
On Thu, Sep 8, 2011 at 8:22 AM, Gilles Sadowski <
gilles@harfang.homelinux.org> wrote:
> > > I've tried it with my problem, and it also throws an exception.
> > > However, I would like to obtain the covariance matrix anyway, because
> I've
> > > no other clue as to what might be wrong.
> > > So I think that, at least, users should be able to set the positive
> > > definiteness threshold in order to avoid raising an exception.
>
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Phil Steitz <ph...@gmail.com>.
On 9/9/11 7:31 AM, Gilles Sadowski wrote:
> Hi.
>
>>> On Thu, Sep 08, 2011 at 10:27:00AM -0500, Greg Sterijevski wrote:
>>>> Could you be contaminating the cross product matrix when you do your
>>>> accumulation? I noticed one of the utility classes has methods for doing
>>>> safe accumulations. I know that this was part of the problem when I was
>>>> trying to estimate the Wampler data with the MillerRegression. Just a
>>>> thought...
>>> To make sure, I've tried replacing the "sum" calculation in "getCovariances"
>>> with one using "linearCombination for "MathUtils" but it still fails, so it
>>> seems that the (default) threshold is too small, and I need to be able to
>>> change it.
>> Could be the best thing to do here is to make the
>> DecompositionSolver pluggable
> That's what I tried to do a couple of hours ago ;-)
> I.e. I thought that we could have
>
> public double[][] getCovariances(DecompositionSolver solver) {
> // ...
> }
>
> But that won't work with the current design because the solver is an inner
> class of the various "...Decomposition" classes. And not a static one
> because it needs to access the results of the decomposition stored in the
> enclosing class. [Or did I miss something?]
Yeah, sorry, you are right this does not quite work. Another
logical possibility is to define an interface to be implemented by
the decomposition impls to support making the decomposition
algorithm pluggable in situations like this where you really just
want a solver factory. This interface could consist of the static
fectory method:
static DecompositionSolver getSolver(RealMatrix, singularityThreshold)
Note that for this to work, we would have to agree to add a
singularityThreshold constructor parameter for each of the
decomposition impls. The problem with that is that, as discussed on
other threads, the thresholds mean different things for different
decomposition algorithms.
I guess the simplest, for the present problem, is just to leave it
hard-coded to LU, but add a version of getCovariances that takes a
threshold argument.
>
>> (maybe just for the LMOptimizer, but
>> with hooks in the base class)?
> A more flexible "getCovariances" should not be restricted to a single
> optimizer. In fact, currently, I extend "AbstractLeastSquaresOptimizer"
> directly, just to retrieve the covariance matrix; the optimum is not
> computed by a subclass of it. [In the unit test which started this thread,
> the optimum is known (thus: no optimization at all); the next step will be
> to compute an optimum with an optimizer that does not use derivatives (e.g.
> "SimplexOptimizer") but the covariance matrix must still be computed at the
> solution found by that optimizer.]
Yes, that would be the most flexible.
Phil
>
>> That way you could provide an LU
>> solver with whatever threshold you want or choose among the others.
>> Do you know if the jTj matrix you end up trying to invert is
>> analytically singular (like the one in the testNonInvertible test
>> case for LMOptimizer)?
> No, I don't know. But if it is, then there must be a bug somewhere, I think
>
>
> 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
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Gilles Sadowski <gi...@harfang.homelinux.org>.
Hi.
> >
> > On Thu, Sep 08, 2011 at 10:27:00AM -0500, Greg Sterijevski wrote:
> >> Could you be contaminating the cross product matrix when you do your
> >> accumulation? I noticed one of the utility classes has methods for doing
> >> safe accumulations. I know that this was part of the problem when I was
> >> trying to estimate the Wampler data with the MillerRegression. Just a
> >> thought...
> > To make sure, I've tried replacing the "sum" calculation in "getCovariances"
> > with one using "linearCombination for "MathUtils" but it still fails, so it
> > seems that the (default) threshold is too small, and I need to be able to
> > change it.
>
> Could be the best thing to do here is to make the
> DecompositionSolver pluggable
That's what I tried to do a couple of hours ago ;-)
I.e. I thought that we could have
public double[][] getCovariances(DecompositionSolver solver) {
// ...
}
But that won't work with the current design because the solver is an inner
class of the various "...Decomposition" classes. And not a static one
because it needs to access the results of the decomposition stored in the
enclosing class. [Or did I miss something?]
> (maybe just for the LMOptimizer, but
> with hooks in the base class)?
A more flexible "getCovariances" should not be restricted to a single
optimizer. In fact, currently, I extend "AbstractLeastSquaresOptimizer"
directly, just to retrieve the covariance matrix; the optimum is not
computed by a subclass of it. [In the unit test which started this thread,
the optimum is known (thus: no optimization at all); the next step will be
to compute an optimum with an optimizer that does not use derivatives (e.g.
"SimplexOptimizer") but the covariance matrix must still be computed at the
solution found by that optimizer.]
> That way you could provide an LU
> solver with whatever threshold you want or choose among the others.
> Do you know if the jTj matrix you end up trying to invert is
> analytically singular (like the one in the testNonInvertible test
> case for LMOptimizer)?
No, I don't know. But if it is, then there must be a bug somewhere, I think.
Gilles
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Phil Steitz <ph...@gmail.com>.
On 9/9/11 6:00 AM, Gilles Sadowski wrote:
> Hi.
>
> On Thu, Sep 08, 2011 at 10:27:00AM -0500, Greg Sterijevski wrote:
>> Could you be contaminating the cross product matrix when you do your
>> accumulation? I noticed one of the utility classes has methods for doing
>> safe accumulations. I know that this was part of the problem when I was
>> trying to estimate the Wampler data with the MillerRegression. Just a
>> thought...
> To make sure, I've tried replacing the "sum" calculation in "getCovariances"
> with one using "linearCombination for "MathUtils" but it still fails, so it
> seems that the (default) threshold is too small, and I need to be able to
> change it.
Could be the best thing to do here is to make the
DecompositionSolver pluggable (maybe just for the LMOptimizer, but
with hooks in the base class)? That way you could provide an LU
solver with whatever threshold you want or choose among the others.
Do you know if the jTj matrix you end up trying to invert is
analytically singular (like the one in the testNonInvertible test
case for LMOptimizer)?
Phil
>
>
> 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
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Gilles Sadowski <gi...@harfang.homelinux.org>.
Hi.
On Thu, Sep 08, 2011 at 10:27:00AM -0500, Greg Sterijevski wrote:
> Could you be contaminating the cross product matrix when you do your
> accumulation? I noticed one of the utility classes has methods for doing
> safe accumulations. I know that this was part of the problem when I was
> trying to estimate the Wampler data with the MillerRegression. Just a
> thought...
To make sure, I've tried replacing the "sum" calculation in "getCovariances"
with one using "linearCombination for "MathUtils" but it still fails, so it
seems that the (default) threshold is too small, and I need to be able to
change it.
Regards,
Gilles
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Greg Sterijevski <gs...@gmail.com>.
Could you be contaminating the cross product matrix when you do your
accumulation? I noticed one of the utility classes has methods for doing
safe accumulations. I know that this was part of the problem when I was
trying to estimate the Wampler data with the MillerRegression. Just a
thought...
On Thu, Sep 8, 2011 at 10:22 AM, Gilles Sadowski <
gilles@harfang.homelinux.org> wrote:
> On Thu, Sep 08, 2011 at 09:49:12AM -0500, Greg Sterijevski wrote:
> > Cholesky, in my opinion, is not robust as you have discovered. When it
> > encounters a non-psd matrix it gives up. Maybe that is the correct course
> of
> > action, but I still think that when you are using the getCovariance to
> > estimate the curvature in the neighborhood of a point it would be okay to
> > take the generalized inverse and not worry too much about why your matrix
> is
> > bordering on non-PSDness... If your optimization stops on that point,
> that
> > is another story and should be flagged.
> >
>
> In fact, in my current case (which is a unit test), the point is supposed
> to
> be the minimum, by construction. I'm trying to figure out where the
> problem comes from (namely, whether the Jacobian matrix is correct)...
>
> Gilles
>
> >
> > > 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.
> > >
> > > I've tried it with my problem, and it also throws an exception.
> > > However, I would like to obtain the covariance matrix anyway, because
> I've
> > > no other clue as to what might be wrong.
> > > So I think that, at least, users should be able to set the positive
> > > definiteness threshold in order to avoid raising an exception.
> > >
> > >
> > >
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
>
>
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Gilles Sadowski <gi...@harfang.homelinux.org>.
On Thu, Sep 08, 2011 at 09:49:12AM -0500, Greg Sterijevski wrote:
> Cholesky, in my opinion, is not robust as you have discovered. When it
> encounters a non-psd matrix it gives up. Maybe that is the correct course of
> action, but I still think that when you are using the getCovariance to
> estimate the curvature in the neighborhood of a point it would be okay to
> take the generalized inverse and not worry too much about why your matrix is
> bordering on non-PSDness... If your optimization stops on that point, that
> is another story and should be flagged.
>
In fact, in my current case (which is a unit test), the point is supposed to
be the minimum, by construction. I'm trying to figure out where the
problem comes from (namely, whether the Jacobian matrix is correct)...
Gilles
>
> > 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.
> >
> > I've tried it with my problem, and it also throws an exception.
> > However, I would like to obtain the covariance matrix anyway, because I've
> > no other clue as to what might be wrong.
> > So I think that, at least, users should be able to set the positive
> > definiteness threshold in order to avoid raising an exception.
> >
> >
> >
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Greg Sterijevski <gs...@gmail.com>.
Cholesky, in my opinion, is not robust as you have discovered. When it
encounters a non-psd matrix it gives up. Maybe that is the correct course of
action, but I still think that when you are using the getCovariance to
estimate the curvature in the neighborhood of a point it would be okay to
take the generalized inverse and not worry too much about why your matrix is
bordering on non-PSDness... If your optimization stops on that point, that
is another story and should be flagged.
> 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.
>
> I've tried it with my problem, and it also throws an exception.
> However, I would like to obtain the covariance matrix anyway, because I've
> no other clue as to what might be wrong.
> So I think that, at least, users should be able to set the positive
> definiteness threshold in order to avoid raising an exception.
>
>
>
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Gilles Sadowski <gi...@harfang.homelinux.org>.
On Thu, Sep 08, 2011 at 12:15:20AM -0700, Phil Steitz wrote:
> 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.
Still, I don't understand: The "problem" is passed to the optimizer which
returns an answer. Is the answer correct? If not, what if the user does not
call "getCovariances()"?
The failure should be detected earlier (possibly by the optimization
algorithm). And if it is not possible during the search, shouldn't there be
a utility method explicitly aimed at checking the quality of the solution?
> 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.
I've tried it with my problem, and it also throws an exception.
However, I would like to obtain the covariance matrix anyway, because I've
no other clue as to what might be wrong.
So I think that, at least, users should be able to set the positive
definiteness threshold in order to avoid raising an exception.
>
> 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.
What do you mean by "with the defaults"?
Regards,
Gilles
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Gilles Sadowski <gi...@harfang.homelinux.org>.
Hi.
> Is the test checked in? I am interested in looking at it.
Thanks for the interest, but it is a unit test for the software which I
work on. I'm afraid that a discussion about this will be fairly OT here.
Basically, I need to compute the covariance matrix, but the model function
is the result of a simulation: thus, currently I compute the Jacobian
numerically (using a trivial approximation of the derivatives).
One of the issues is how to make sure that the resulting matrix is
meaningful.
Another issue, further up, is how to make sure that the Jacobian is
meaningful. We also want to test optimizers that use derivatives (i.e.
"LevenbergMarquardtOptimizer") but, since the parameters are subject to
(simple) boundary conditions (e.g. orbit eccentricity, temperature, ...),
I use a "logit" transform when passing the values to the optimizer. This
leads to a (transformed) Jacobian whose columns have very different orders
of magnitudes. So I wonder whether the "transform" approach is usable at
all...
Best,
Gilles
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Greg Sterijevski <gs...@gmail.com>.
Gilles,
Is the test checked in? I am interested in looking at it.
Thank you,
-Greg
On Thu, Sep 8, 2011 at 2:51 PM, Ted Dunning <te...@gmail.com> wrote:
> On Thu, Sep 8, 2011 at 12:42 PM, Luc Maisonobe <Luc.Maisonobe@free.fr
> >wrote:
>
> > ... 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).
>
>
> This is my impression as well.
>
> In addition, Cholesky very nearly is a QR decomposition through the back
> door. That is, if
>
> Q R = A
>
> Then
>
> R' R = A' A
>
> is the Cholesky decomposition. The algorithms are quite similar when
> viewed
> this way. Cholesky does not produce as accurate a result for R if you are
> given A, but if you are given A'A as in the discussion here, it is pretty
> much just as good, I think.
>
> I use this property for very large SVD's via stochastic projection because
> Q
> is much larger than R in that context so computing just R can be done
> in-core while the full QR would require out-of-core operations.
>
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Ted Dunning <te...@gmail.com>.
On Thu, Sep 8, 2011 at 12:42 PM, Luc Maisonobe <Lu...@free.fr>wrote:
> ... 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).
This is my impression as well.
In addition, Cholesky very nearly is a QR decomposition through the back
door. That is, if
Q R = A
Then
R' R = A' A
is the Cholesky decomposition. The algorithms are quite similar when viewed
this way. Cholesky does not produce as accurate a result for R if you are
given A, but if you are given A'A as in the discussion here, it is pretty
much just as good, I think.
I use this property for very large SVD's via stochastic projection because Q
is much larger than R in that context so computing just R can be done
in-core while the full QR would require out-of-core operations.
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Luc Maisonobe <Lu...@free.fr>.
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
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Phil Steitz <ph...@gmail.com>.
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.
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
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Gilles Sadowski <gi...@harfang.homelinux.org>.
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)?
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
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Luc Maisonobe <Lu...@free.fr>.
Hi Gilles,
Le 07/09/2011 15:05, Gilles Sadowski a écrit :
> On Wed, Sep 07, 2011 at 02:46:59PM +0200, Luc Maisonobe wrote:
>> Le 07/09/2011 12:45, Gilles Sadowski a écrit :
>>> Hello.
>>>
>>> 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.
Luc
>
> Thanks,
> 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
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Gilles Sadowski <gi...@harfang.homelinux.org>.
On Wed, Sep 07, 2011 at 02:46:59PM +0200, Luc Maisonobe wrote:
> Le 07/09/2011 12:45, Gilles Sadowski a écrit :
> >Hello.
> >
> >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?
Thanks,
Gilles
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [Math] "LUDecomposition" in "AbstractLeastSquaresOptimizer"
Posted by Luc Maisonobe <Lu...@free.fr>.
Le 07/09/2011 12:45, Gilles Sadowski a écrit :
> Hello.
>
> 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.
Luc
>
>
> 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