You are viewing a plain text version of this content. The canonical link for it is here.
Posted to dev@mahout.apache.org by Dmitriy Lyubimov <dl...@gmail.com> on 2013/07/16 02:00:50 UTC
Cholesky test failures...
Hi Ted,
I am getting Cholesky test failures when trying to solve of Ax=B
The L matrix and solveLeft() output do not make much sense to me. For once,
L doesn't even have the expected L-shape:
Do you have an idea where i go wrong? (the test is wrapped into scala DSL
but it is Mahout's cholesky underwraps) .
test code:
test("chol") {
// try to solve Ax=b with cholesky:
// this requires
// (LL')x = B
// L'x= (L^-1)B
// x=(L'^-1)(L^-1)B
val a = dense((1, 2, 3), (2, 3, 4), (3, 4, 5))
// make sure it is symmetric for a valid solution
a := a.t %*% a
printf("A= \n%s\n", a)
val b = dense((9, 8, 7)).t
printf ("b = \n%s\n", b)
val ch = chol(a)
printf ("L = \n%s\n", ch.getL)
printf ( "(L^-1)b =\n%s\n", ch.solveLeft(b))
val x = ch.solveRight(diag(1,3)) %*% ch.solveLeft(b)
printf("x = \n%s\n", x.toString)
val axmb = (a %*% b) - b
printf("AX - B = \n%s\n", axmb.toString)
assert(axmb.norm < 1e-10)
}
Output:
A=
{
0 => {0:14.0,1:20.0,2:26.0}
1 => {0:20.0,1:29.0,2:38.0}
2 => {0:26.0,1:38.0,2:50.0}
}
b =
{
0 => {0:9.0}
1 => {0:8.0}
2 => {0:7.0}
}
L =
{
0 => {0:0.6928203230275511,2:3.676955262170047}
1 => {0:0.3464101615137781,2:5.374011537017761}
2 => {2:7.0710678118654755}
}
(L^-1)b =
{
0 => {0:1.2727922061357855}
1 => {0:11.547005383792511}
2 => {}
}
X =
{
0 => {0:0.18}
1 => {0:5.119661282874144}
2 => {}
}
AX - B =
{
0 => {0:459.0}
1 => {0:670.0}
2 => {0:881.0}
}
org.scalatest.exceptions.TestFailedException was thrown.
Re: Cholesky test failures...
Posted by Ted Dunning <te...@gmail.com>.
On Thu, Jul 18, 2013 at 12:23 AM, Dmitriy Lyubimov <dl...@gmail.com>wrote:
> You are at HEAD of the trunk?
>
> Ok i will tinker a little more. But i am pretty sure different results are
> produced for the same input if i just flip pivoted flag. I can add a java
> version of the failing test as well if it is easier.
>
I have a java version. Still walking through your steps to see what is
happening.
Having another java version is not a bad thing.
> I also noticed that R version is more accurate in detecting non positive
> definite situations while Mahout version sometimes would accept it and
> produce what seems like a nonsense.
>
Hmm... were these marginal situations? Or is the test in Mahout's version
just stupid?
> oh well. it seems to be working for me in my Bagel-based ALS-WR stuff for
> now if i accurately assert valid input and don't use pivoted=true in
> constructor.
>
OK.
> Also, there seems to be no benefit in building classic thin QR
> implementation then? I suspect A'A step could make this route more
> expensive than a dedicated QR method, but i am not sure.
>
Correct.
But I am not entirely clear how good the numerical properties are. It
would be good to have somebody at least measure this.
Not constructing Q except on the fly as it is needed is very nice.
>
>
> On Wed, Jul 17, 2013 at 3:59 PM, Ted Dunning <te...@gmail.com>
> wrote:
>
> > My time budget for this problem for today is now exhausted.
> >
> > What I have determined so far:
> >
> > 1) CholeskyDecomposition produces correctly triangular matrices in your
> > example (with or without pivoting)
> >
> > 2) multiplication of these matrices reproduces the input as expected
> (with
> > or without pivoting)
> >
> > 3) I am not looking clearly enough at the test of your example to
> > understand the correct behavior.
> >
> > My route from this point is to start with Cholesky as a shortcut for QR
> > decomposition via the identity:
> >
> > A = QR
> > (QR)' (QR) = A' A = R' (Q' Q) R = R' R = L L'
> >
> > Thus A (L')^-1 = Q
> >
> > I will be checking that Q is produced and is orthonormal as expected.
> >
> > I begin to think that the overall problem is a misunderstanding of the
> API
> > rather than a functional error. That should be cured with better
> > documentation, but it isn't clear that this will be easier than a code
> bug.
> >
> >
> >
> >
> > On Wed, Jul 17, 2013 at 3:13 PM, Dmitriy Lyubimov <dl...@gmail.com>
> > wrote:
> >
> > > sorry, i mean ch.solveRight(eye(4)) fails with AOOB
> > >
> > >
> > > On Wed, Jul 17, 2013 at 3:11 PM, Dmitriy Lyubimov <dl...@gmail.com>
> > > wrote:
> > >
> > > > if you have scala plugin installed in idea, i have this scala DSL
> > module
> > > > added to mahout on this branch
> > > > https://github.com/dlyubimov/mahout-commits/tree/dev-0.9.x-scala
> > > >
> > > >
> > > > file
> > >
> >
> /home/dmitriy/projects/asf/mahout-commits/math-scala/src/test/scala/mahout/math/MatrixOpsTest.scala
> > > >
> > > > this test fails if used with chol(a, true) or if
> rightMultiply(eye(4)):
> > > >
> > > > test("chol") {
> > > >
> > > > // try to solve Ax=b with cholesky:
> > > > // this requires
> > > > // (LL')x = B
> > > > // L'x= (L^-1)B
> > > > // x=(L'^-1)(L^-1)B
> > > >
> > > > val a = dense((1, 2, 3), (2, 3, 4), (3, 4, 5.5))
> > > >
> > > > // make sure it is symmetric for a valid solution
> > > > a := a.t %*% a
> > > >
> > > > printf("A= \n%s\n", a)
> > > >
> > > > val b = dense((9, 8, 7)).t
> > > >
> > > > printf("b = \n%s\n", b)
> > > >
> > > > // fails if chol(a,true)
> > > > val ch = chol(a)
> > > >
> > > > printf("L = \n%s\n", ch.getL)
> > > >
> > > > printf("(L^-1)b =\n%s\n", ch.solveLeft(b))
> > > >
> > > > val x = ch.solveRight(eye(3)) %*% ch.solveLeft(b)
> > > >
> > > > printf("x = \n%s\n", x.toString)
> > > >
> > > > val axmb = (a %*% x) - b
> > > >
> > > > printf("AX - B = \n%s\n", axmb.toString)
> > > >
> > > > assert(axmb.norm < 1e-10)
> > > >
> > > > }
> > > >
> > > >
> > > >
> > > >
> > > > On Wed, Jul 17, 2013 at 2:28 PM, Ted Dunning <ted.dunning@gmail.com
> > > >wrote:
> > > >
> > > >> These problems are very strange.
> > > >>
> > > >> I am now looking at the tests for Cholesky and it seems that they
> > cover
> > > >> all
> > > >> of the sorts of things that you are talking about.
> > > >>
> > > >> I can fix the size compatibility test and will add a test that
> > > implements
> > > >> your other issue to see if that helps me understand what is
> happening.
> > > >>
> > > >>
> > > >> On Wed, Jul 17, 2013 at 1:02 PM, Dmitriy Lyubimov <
> dlieu.7@gmail.com>
> > > >> wrote:
> > > >>
> > > >> > btw, another nitpicking, solveRight(eye(n) ) and solveLeft() do
> not
> > > >> check
> > > >> > for cardinality of argument, throwing ArrayOutOfBounds. yes the
> > burden
> > > >> of
> > > >> > dumbness is on the user but the burden of explanation is on
> > > >> implementation
> > > >> > :)
> > > >> >
> > > >>
> > > >
> > > >
> > >
> >
>
Re: Cholesky test failures...
Posted by Dmitriy Lyubimov <dl...@gmail.com>.
You are at HEAD of the trunk?
Ok i will tinker a little more. But i am pretty sure different results are
produced for the same input if i just flip pivoted flag. I can add a java
version of the failing test as well if it is easier.
I also noticed that R version is more accurate in detecting non positive
definite situations while Mahout version sometimes would accept it and
produce what seems like a nonsense.
oh well. it seems to be working for me in my Bagel-based ALS-WR stuff for
now if i accurately assert valid input and don't use pivoted=true in
constructor.
Also, there seems to be no benefit in building classic thin QR
implementation then? I suspect A'A step could make this route more
expensive than a dedicated QR method, but i am not sure.
On Wed, Jul 17, 2013 at 3:59 PM, Ted Dunning <te...@gmail.com> wrote:
> My time budget for this problem for today is now exhausted.
>
> What I have determined so far:
>
> 1) CholeskyDecomposition produces correctly triangular matrices in your
> example (with or without pivoting)
>
> 2) multiplication of these matrices reproduces the input as expected (with
> or without pivoting)
>
> 3) I am not looking clearly enough at the test of your example to
> understand the correct behavior.
>
> My route from this point is to start with Cholesky as a shortcut for QR
> decomposition via the identity:
>
> A = QR
> (QR)' (QR) = A' A = R' (Q' Q) R = R' R = L L'
>
> Thus A (L')^-1 = Q
>
> I will be checking that Q is produced and is orthonormal as expected.
>
> I begin to think that the overall problem is a misunderstanding of the API
> rather than a functional error. That should be cured with better
> documentation, but it isn't clear that this will be easier than a code bug.
>
>
>
>
> On Wed, Jul 17, 2013 at 3:13 PM, Dmitriy Lyubimov <dl...@gmail.com>
> wrote:
>
> > sorry, i mean ch.solveRight(eye(4)) fails with AOOB
> >
> >
> > On Wed, Jul 17, 2013 at 3:11 PM, Dmitriy Lyubimov <dl...@gmail.com>
> > wrote:
> >
> > > if you have scala plugin installed in idea, i have this scala DSL
> module
> > > added to mahout on this branch
> > > https://github.com/dlyubimov/mahout-commits/tree/dev-0.9.x-scala
> > >
> > >
> > > file
> >
> /home/dmitriy/projects/asf/mahout-commits/math-scala/src/test/scala/mahout/math/MatrixOpsTest.scala
> > >
> > > this test fails if used with chol(a, true) or if rightMultiply(eye(4)):
> > >
> > > test("chol") {
> > >
> > > // try to solve Ax=b with cholesky:
> > > // this requires
> > > // (LL')x = B
> > > // L'x= (L^-1)B
> > > // x=(L'^-1)(L^-1)B
> > >
> > > val a = dense((1, 2, 3), (2, 3, 4), (3, 4, 5.5))
> > >
> > > // make sure it is symmetric for a valid solution
> > > a := a.t %*% a
> > >
> > > printf("A= \n%s\n", a)
> > >
> > > val b = dense((9, 8, 7)).t
> > >
> > > printf("b = \n%s\n", b)
> > >
> > > // fails if chol(a,true)
> > > val ch = chol(a)
> > >
> > > printf("L = \n%s\n", ch.getL)
> > >
> > > printf("(L^-1)b =\n%s\n", ch.solveLeft(b))
> > >
> > > val x = ch.solveRight(eye(3)) %*% ch.solveLeft(b)
> > >
> > > printf("x = \n%s\n", x.toString)
> > >
> > > val axmb = (a %*% x) - b
> > >
> > > printf("AX - B = \n%s\n", axmb.toString)
> > >
> > > assert(axmb.norm < 1e-10)
> > >
> > > }
> > >
> > >
> > >
> > >
> > > On Wed, Jul 17, 2013 at 2:28 PM, Ted Dunning <ted.dunning@gmail.com
> > >wrote:
> > >
> > >> These problems are very strange.
> > >>
> > >> I am now looking at the tests for Cholesky and it seems that they
> cover
> > >> all
> > >> of the sorts of things that you are talking about.
> > >>
> > >> I can fix the size compatibility test and will add a test that
> > implements
> > >> your other issue to see if that helps me understand what is happening.
> > >>
> > >>
> > >> On Wed, Jul 17, 2013 at 1:02 PM, Dmitriy Lyubimov <dl...@gmail.com>
> > >> wrote:
> > >>
> > >> > btw, another nitpicking, solveRight(eye(n) ) and solveLeft() do not
> > >> check
> > >> > for cardinality of argument, throwing ArrayOutOfBounds. yes the
> burden
> > >> of
> > >> > dumbness is on the user but the burden of explanation is on
> > >> implementation
> > >> > :)
> > >> >
> > >>
> > >
> > >
> >
>
Re: Cholesky test failures...
Posted by Ted Dunning <te...@gmail.com>.
My time budget for this problem for today is now exhausted.
What I have determined so far:
1) CholeskyDecomposition produces correctly triangular matrices in your
example (with or without pivoting)
2) multiplication of these matrices reproduces the input as expected (with
or without pivoting)
3) I am not looking clearly enough at the test of your example to
understand the correct behavior.
My route from this point is to start with Cholesky as a shortcut for QR
decomposition via the identity:
A = QR
(QR)' (QR) = A' A = R' (Q' Q) R = R' R = L L'
Thus A (L')^-1 = Q
I will be checking that Q is produced and is orthonormal as expected.
I begin to think that the overall problem is a misunderstanding of the API
rather than a functional error. That should be cured with better
documentation, but it isn't clear that this will be easier than a code bug.
On Wed, Jul 17, 2013 at 3:13 PM, Dmitriy Lyubimov <dl...@gmail.com> wrote:
> sorry, i mean ch.solveRight(eye(4)) fails with AOOB
>
>
> On Wed, Jul 17, 2013 at 3:11 PM, Dmitriy Lyubimov <dl...@gmail.com>
> wrote:
>
> > if you have scala plugin installed in idea, i have this scala DSL module
> > added to mahout on this branch
> > https://github.com/dlyubimov/mahout-commits/tree/dev-0.9.x-scala
> >
> >
> > file
> /home/dmitriy/projects/asf/mahout-commits/math-scala/src/test/scala/mahout/math/MatrixOpsTest.scala
> >
> > this test fails if used with chol(a, true) or if rightMultiply(eye(4)):
> >
> > test("chol") {
> >
> > // try to solve Ax=b with cholesky:
> > // this requires
> > // (LL')x = B
> > // L'x= (L^-1)B
> > // x=(L'^-1)(L^-1)B
> >
> > val a = dense((1, 2, 3), (2, 3, 4), (3, 4, 5.5))
> >
> > // make sure it is symmetric for a valid solution
> > a := a.t %*% a
> >
> > printf("A= \n%s\n", a)
> >
> > val b = dense((9, 8, 7)).t
> >
> > printf("b = \n%s\n", b)
> >
> > // fails if chol(a,true)
> > val ch = chol(a)
> >
> > printf("L = \n%s\n", ch.getL)
> >
> > printf("(L^-1)b =\n%s\n", ch.solveLeft(b))
> >
> > val x = ch.solveRight(eye(3)) %*% ch.solveLeft(b)
> >
> > printf("x = \n%s\n", x.toString)
> >
> > val axmb = (a %*% x) - b
> >
> > printf("AX - B = \n%s\n", axmb.toString)
> >
> > assert(axmb.norm < 1e-10)
> >
> > }
> >
> >
> >
> >
> > On Wed, Jul 17, 2013 at 2:28 PM, Ted Dunning <ted.dunning@gmail.com
> >wrote:
> >
> >> These problems are very strange.
> >>
> >> I am now looking at the tests for Cholesky and it seems that they cover
> >> all
> >> of the sorts of things that you are talking about.
> >>
> >> I can fix the size compatibility test and will add a test that
> implements
> >> your other issue to see if that helps me understand what is happening.
> >>
> >>
> >> On Wed, Jul 17, 2013 at 1:02 PM, Dmitriy Lyubimov <dl...@gmail.com>
> >> wrote:
> >>
> >> > btw, another nitpicking, solveRight(eye(n) ) and solveLeft() do not
> >> check
> >> > for cardinality of argument, throwing ArrayOutOfBounds. yes the burden
> >> of
> >> > dumbness is on the user but the burden of explanation is on
> >> implementation
> >> > :)
> >> >
> >>
> >
> >
>
Re: Cholesky test failures...
Posted by Dmitriy Lyubimov <dl...@gmail.com>.
sorry, i mean ch.solveRight(eye(4)) fails with AOOB
On Wed, Jul 17, 2013 at 3:11 PM, Dmitriy Lyubimov <dl...@gmail.com> wrote:
> if you have scala plugin installed in idea, i have this scala DSL module
> added to mahout on this branch
> https://github.com/dlyubimov/mahout-commits/tree/dev-0.9.x-scala
>
>
> file /home/dmitriy/projects/asf/mahout-commits/math-scala/src/test/scala/mahout/math/MatrixOpsTest.scala
>
> this test fails if used with chol(a, true) or if rightMultiply(eye(4)):
>
> test("chol") {
>
> // try to solve Ax=b with cholesky:
> // this requires
> // (LL')x = B
> // L'x= (L^-1)B
> // x=(L'^-1)(L^-1)B
>
> val a = dense((1, 2, 3), (2, 3, 4), (3, 4, 5.5))
>
> // make sure it is symmetric for a valid solution
> a := a.t %*% a
>
> printf("A= \n%s\n", a)
>
> val b = dense((9, 8, 7)).t
>
> printf("b = \n%s\n", b)
>
> // fails if chol(a,true)
> val ch = chol(a)
>
> printf("L = \n%s\n", ch.getL)
>
> printf("(L^-1)b =\n%s\n", ch.solveLeft(b))
>
> val x = ch.solveRight(eye(3)) %*% ch.solveLeft(b)
>
> printf("x = \n%s\n", x.toString)
>
> val axmb = (a %*% x) - b
>
> printf("AX - B = \n%s\n", axmb.toString)
>
> assert(axmb.norm < 1e-10)
>
> }
>
>
>
>
> On Wed, Jul 17, 2013 at 2:28 PM, Ted Dunning <te...@gmail.com>wrote:
>
>> These problems are very strange.
>>
>> I am now looking at the tests for Cholesky and it seems that they cover
>> all
>> of the sorts of things that you are talking about.
>>
>> I can fix the size compatibility test and will add a test that implements
>> your other issue to see if that helps me understand what is happening.
>>
>>
>> On Wed, Jul 17, 2013 at 1:02 PM, Dmitriy Lyubimov <dl...@gmail.com>
>> wrote:
>>
>> > btw, another nitpicking, solveRight(eye(n) ) and solveLeft() do not
>> check
>> > for cardinality of argument, throwing ArrayOutOfBounds. yes the burden
>> of
>> > dumbness is on the user but the burden of explanation is on
>> implementation
>> > :)
>> >
>>
>
>
Re: Cholesky test failures...
Posted by Dmitriy Lyubimov <dl...@gmail.com>.
if you have scala plugin installed in idea, i have this scala DSL module
added to mahout on this branch
https://github.com/dlyubimov/mahout-commits/tree/dev-0.9.x-scala
file /home/dmitriy/projects/asf/mahout-commits/math-scala/src/test/scala/mahout/math/MatrixOpsTest.scala
this test fails if used with chol(a, true) or if rightMultiply(eye(4)):
test("chol") {
// try to solve Ax=b with cholesky:
// this requires
// (LL')x = B
// L'x= (L^-1)B
// x=(L'^-1)(L^-1)B
val a = dense((1, 2, 3), (2, 3, 4), (3, 4, 5.5))
// make sure it is symmetric for a valid solution
a := a.t %*% a
printf("A= \n%s\n", a)
val b = dense((9, 8, 7)).t
printf("b = \n%s\n", b)
// fails if chol(a,true)
val ch = chol(a)
printf("L = \n%s\n", ch.getL)
printf("(L^-1)b =\n%s\n", ch.solveLeft(b))
val x = ch.solveRight(eye(3)) %*% ch.solveLeft(b)
printf("x = \n%s\n", x.toString)
val axmb = (a %*% x) - b
printf("AX - B = \n%s\n", axmb.toString)
assert(axmb.norm < 1e-10)
}
On Wed, Jul 17, 2013 at 2:28 PM, Ted Dunning <te...@gmail.com> wrote:
> These problems are very strange.
>
> I am now looking at the tests for Cholesky and it seems that they cover all
> of the sorts of things that you are talking about.
>
> I can fix the size compatibility test and will add a test that implements
> your other issue to see if that helps me understand what is happening.
>
>
> On Wed, Jul 17, 2013 at 1:02 PM, Dmitriy Lyubimov <dl...@gmail.com>
> wrote:
>
> > btw, another nitpicking, solveRight(eye(n) ) and solveLeft() do not check
> > for cardinality of argument, throwing ArrayOutOfBounds. yes the burden of
> > dumbness is on the user but the burden of explanation is on
> implementation
> > :)
> >
>
Re: Cholesky test failures...
Posted by Ted Dunning <te...@gmail.com>.
These problems are very strange.
I am now looking at the tests for Cholesky and it seems that they cover all
of the sorts of things that you are talking about.
I can fix the size compatibility test and will add a test that implements
your other issue to see if that helps me understand what is happening.
On Wed, Jul 17, 2013 at 1:02 PM, Dmitriy Lyubimov <dl...@gmail.com> wrote:
> btw, another nitpicking, solveRight(eye(n) ) and solveLeft() do not check
> for cardinality of argument, throwing ArrayOutOfBounds. yes the burden of
> dumbness is on the user but the burden of explanation is on implementation
> :)
>
Re: Cholesky test failures...
Posted by Dmitriy Lyubimov <dl...@gmail.com>.
btw, another nitpicking, solveRight(eye(n) ) and solveLeft() do not check
for cardinality of argument, throwing ArrayOutOfBounds. yes the burden of
dumbness is on the user but the burden of explanation is on implementation
:)
Re: Cholesky test failures...
Posted by Ted Dunning <te...@gmail.com>.
On Mon, Jul 15, 2013 at 5:16 PM, Dmitriy Lyubimov <dl...@gmail.com> wrote:
> On Mon, Jul 15, 2013 at 5:13 PM, Ted Dunning <te...@gmail.com>
> wrote:
>
> > It sounds like you missed the part of the API contract that says that you
> > have to fix the bugs before using the code.
> >
> > :-)
> >
> > More seriously, with pivot=true, the result can be hard to understand
> > because they are permuted. Could that be the problem?
> >
>
> ok, but solveLeft() in this case should still work as expected? or not ? I
> don't understand the "permuted" part -- i guess it is some sort of indexing
> overlay over L matrix -- but whatever solution is, shouldn't the api be
> agnostic of the method used?
>
Yes. It should. And there should be a test case that proves it.
I am surprised that there isn't.
The idea behind the pivoting is that you get AP = QR so
(AP)'(AP) = R'R = P' A'A P
thus
A' A = P R' R P' = (R P')' (R P')
What you get in place of R is RP' (R with columns rearranged).
Re: Cholesky test failures...
Posted by Dmitriy Lyubimov <dl...@gmail.com>.
On Mon, Jul 15, 2013 at 5:13 PM, Ted Dunning <te...@gmail.com> wrote:
> It sounds like you missed the part of the API contract that says that you
> have to fix the bugs before using the code.
>
> :-)
>
> More seriously, with pivot=true, the result can be hard to understand
> because they are permuted. Could that be the problem?
>
ok, but solveLeft() in this case should still work as expected? or not ? I
don't understand the "permuted" part -- i guess it is some sort of indexing
overlay over L matrix -- but whatever solution is, shouldn't the api be
agnostic of the method used?
>
>
> On Mon, Jul 15, 2013 at 5:11 PM, Dmitriy Lyubimov <dl...@gmail.com>
> wrote:
>
> > Hm.
> >
> > if i specify pivoted=false, everything works.
> > In addition it seems i have created singular input by chance, but making
> it
> > non-singular
> >
> >
> > val a = dense((1, 2, 3), (2, 3, 4), (3, 4, 5.5))
> >
> > still doesn't help with pivoted=true (which is default) , my test fails.
> >
> > do pivoted results require some special handling?
> >
> >
> >
> >
> > On Mon, Jul 15, 2013 at 5:03 PM, Dmitriy Lyubimov <dl...@gmail.com>
> > wrote:
> >
> > > should read
> > >
> > > val axmb = (a %*% x) - b
> > >
> > > of course but it doesn't make difference, the Cholesky output doesn't
> > make
> > > sense to me even before that
> > >
> > >
> > > On Mon, Jul 15, 2013 at 5:00 PM, Dmitriy Lyubimov <dlieu.7@gmail.com
> > >wrote:
> > >
> > >> Hi Ted,
> > >>
> > >> I am getting Cholesky test failures when trying to solve of Ax=B
> > >>
> > >> The L matrix and solveLeft() output do not make much sense to me. For
> > >> once, L doesn't even have the expected L-shape:
> > >>
> > >> Do you have an idea where i go wrong? (the test is wrapped into scala
> > DSL
> > >> but it is Mahout's cholesky underwraps) .
> > >>
> > >> test code:
> > >>
> > >> test("chol") {
> > >>
> > >> // try to solve Ax=b with cholesky:
> > >> // this requires
> > >> // (LL')x = B
> > >> // L'x= (L^-1)B
> > >> // x=(L'^-1)(L^-1)B
> > >>
> > >> val a = dense((1, 2, 3), (2, 3, 4), (3, 4, 5))
> > >>
> > >> // make sure it is symmetric for a valid solution
> > >> a := a.t %*% a
> > >>
> > >> printf("A= \n%s\n", a)
> > >>
> > >> val b = dense((9, 8, 7)).t
> > >>
> > >> printf ("b = \n%s\n", b)
> > >>
> > >> val ch = chol(a)
> > >>
> > >> printf ("L = \n%s\n", ch.getL)
> > >>
> > >> printf ( "(L^-1)b =\n%s\n", ch.solveLeft(b))
> > >>
> > >> val x = ch.solveRight(diag(1,3)) %*% ch.solveLeft(b)
> > >>
> > >> printf("x = \n%s\n", x.toString)
> > >>
> > >> val axmb = (a %*% b) - b
> > >>
> > >> printf("AX - B = \n%s\n", axmb.toString)
> > >>
> > >> assert(axmb.norm < 1e-10)
> > >>
> > >> }
> > >>
> > >>
> > >>
> > >> Output:
> > >>
> > >> A=
> > >> {
> > >> 0 => {0:14.0,1:20.0,2:26.0}
> > >> 1 => {0:20.0,1:29.0,2:38.0}
> > >> 2 => {0:26.0,1:38.0,2:50.0}
> > >> }
> > >> b =
> > >> {
> > >> 0 => {0:9.0}
> > >> 1 => {0:8.0}
> > >> 2 => {0:7.0}
> > >> }
> > >>
> > >> L =
> > >> {
> > >> 0 => {0:0.6928203230275511,2:3.676955262170047}
> > >> 1 => {0:0.3464101615137781,2:5.374011537017761}
> > >> 2 => {2:7.0710678118654755}
> > >> }
> > >> (L^-1)b =
> > >> {
> > >> 0 => {0:1.2727922061357855}
> > >> 1 => {0:11.547005383792511}
> > >> 2 => {}
> > >> }
> > >> X =
> > >> {
> > >> 0 => {0:0.18}
> > >> 1 => {0:5.119661282874144}
> > >> 2 => {}
> > >> }
> > >> AX - B =
> > >> {
> > >> 0 => {0:459.0}
> > >> 1 => {0:670.0}
> > >> 2 => {0:881.0}
> > >> }
> > >>
> > >> org.scalatest.exceptions.TestFailedException was thrown.
> > >>
> > >>
> > >
> >
>
Re: Cholesky test failures...
Posted by Ted Dunning <te...@gmail.com>.
It sounds like you missed the part of the API contract that says that you
have to fix the bugs before using the code.
:-)
More seriously, with pivot=true, the result can be hard to understand
because they are permuted. Could that be the problem?
On Mon, Jul 15, 2013 at 5:11 PM, Dmitriy Lyubimov <dl...@gmail.com> wrote:
> Hm.
>
> if i specify pivoted=false, everything works.
> In addition it seems i have created singular input by chance, but making it
> non-singular
>
>
> val a = dense((1, 2, 3), (2, 3, 4), (3, 4, 5.5))
>
> still doesn't help with pivoted=true (which is default) , my test fails.
>
> do pivoted results require some special handling?
>
>
>
>
> On Mon, Jul 15, 2013 at 5:03 PM, Dmitriy Lyubimov <dl...@gmail.com>
> wrote:
>
> > should read
> >
> > val axmb = (a %*% x) - b
> >
> > of course but it doesn't make difference, the Cholesky output doesn't
> make
> > sense to me even before that
> >
> >
> > On Mon, Jul 15, 2013 at 5:00 PM, Dmitriy Lyubimov <dlieu.7@gmail.com
> >wrote:
> >
> >> Hi Ted,
> >>
> >> I am getting Cholesky test failures when trying to solve of Ax=B
> >>
> >> The L matrix and solveLeft() output do not make much sense to me. For
> >> once, L doesn't even have the expected L-shape:
> >>
> >> Do you have an idea where i go wrong? (the test is wrapped into scala
> DSL
> >> but it is Mahout's cholesky underwraps) .
> >>
> >> test code:
> >>
> >> test("chol") {
> >>
> >> // try to solve Ax=b with cholesky:
> >> // this requires
> >> // (LL')x = B
> >> // L'x= (L^-1)B
> >> // x=(L'^-1)(L^-1)B
> >>
> >> val a = dense((1, 2, 3), (2, 3, 4), (3, 4, 5))
> >>
> >> // make sure it is symmetric for a valid solution
> >> a := a.t %*% a
> >>
> >> printf("A= \n%s\n", a)
> >>
> >> val b = dense((9, 8, 7)).t
> >>
> >> printf ("b = \n%s\n", b)
> >>
> >> val ch = chol(a)
> >>
> >> printf ("L = \n%s\n", ch.getL)
> >>
> >> printf ( "(L^-1)b =\n%s\n", ch.solveLeft(b))
> >>
> >> val x = ch.solveRight(diag(1,3)) %*% ch.solveLeft(b)
> >>
> >> printf("x = \n%s\n", x.toString)
> >>
> >> val axmb = (a %*% b) - b
> >>
> >> printf("AX - B = \n%s\n", axmb.toString)
> >>
> >> assert(axmb.norm < 1e-10)
> >>
> >> }
> >>
> >>
> >>
> >> Output:
> >>
> >> A=
> >> {
> >> 0 => {0:14.0,1:20.0,2:26.0}
> >> 1 => {0:20.0,1:29.0,2:38.0}
> >> 2 => {0:26.0,1:38.0,2:50.0}
> >> }
> >> b =
> >> {
> >> 0 => {0:9.0}
> >> 1 => {0:8.0}
> >> 2 => {0:7.0}
> >> }
> >>
> >> L =
> >> {
> >> 0 => {0:0.6928203230275511,2:3.676955262170047}
> >> 1 => {0:0.3464101615137781,2:5.374011537017761}
> >> 2 => {2:7.0710678118654755}
> >> }
> >> (L^-1)b =
> >> {
> >> 0 => {0:1.2727922061357855}
> >> 1 => {0:11.547005383792511}
> >> 2 => {}
> >> }
> >> X =
> >> {
> >> 0 => {0:0.18}
> >> 1 => {0:5.119661282874144}
> >> 2 => {}
> >> }
> >> AX - B =
> >> {
> >> 0 => {0:459.0}
> >> 1 => {0:670.0}
> >> 2 => {0:881.0}
> >> }
> >>
> >> org.scalatest.exceptions.TestFailedException was thrown.
> >>
> >>
> >
>
Re: Cholesky test failures...
Posted by Dmitriy Lyubimov <dl...@gmail.com>.
Hm.
if i specify pivoted=false, everything works.
In addition it seems i have created singular input by chance, but making it
non-singular
val a = dense((1, 2, 3), (2, 3, 4), (3, 4, 5.5))
still doesn't help with pivoted=true (which is default) , my test fails.
do pivoted results require some special handling?
On Mon, Jul 15, 2013 at 5:03 PM, Dmitriy Lyubimov <dl...@gmail.com> wrote:
> should read
>
> val axmb = (a %*% x) - b
>
> of course but it doesn't make difference, the Cholesky output doesn't make
> sense to me even before that
>
>
> On Mon, Jul 15, 2013 at 5:00 PM, Dmitriy Lyubimov <dl...@gmail.com>wrote:
>
>> Hi Ted,
>>
>> I am getting Cholesky test failures when trying to solve of Ax=B
>>
>> The L matrix and solveLeft() output do not make much sense to me. For
>> once, L doesn't even have the expected L-shape:
>>
>> Do you have an idea where i go wrong? (the test is wrapped into scala DSL
>> but it is Mahout's cholesky underwraps) .
>>
>> test code:
>>
>> test("chol") {
>>
>> // try to solve Ax=b with cholesky:
>> // this requires
>> // (LL')x = B
>> // L'x= (L^-1)B
>> // x=(L'^-1)(L^-1)B
>>
>> val a = dense((1, 2, 3), (2, 3, 4), (3, 4, 5))
>>
>> // make sure it is symmetric for a valid solution
>> a := a.t %*% a
>>
>> printf("A= \n%s\n", a)
>>
>> val b = dense((9, 8, 7)).t
>>
>> printf ("b = \n%s\n", b)
>>
>> val ch = chol(a)
>>
>> printf ("L = \n%s\n", ch.getL)
>>
>> printf ( "(L^-1)b =\n%s\n", ch.solveLeft(b))
>>
>> val x = ch.solveRight(diag(1,3)) %*% ch.solveLeft(b)
>>
>> printf("x = \n%s\n", x.toString)
>>
>> val axmb = (a %*% b) - b
>>
>> printf("AX - B = \n%s\n", axmb.toString)
>>
>> assert(axmb.norm < 1e-10)
>>
>> }
>>
>>
>>
>> Output:
>>
>> A=
>> {
>> 0 => {0:14.0,1:20.0,2:26.0}
>> 1 => {0:20.0,1:29.0,2:38.0}
>> 2 => {0:26.0,1:38.0,2:50.0}
>> }
>> b =
>> {
>> 0 => {0:9.0}
>> 1 => {0:8.0}
>> 2 => {0:7.0}
>> }
>>
>> L =
>> {
>> 0 => {0:0.6928203230275511,2:3.676955262170047}
>> 1 => {0:0.3464101615137781,2:5.374011537017761}
>> 2 => {2:7.0710678118654755}
>> }
>> (L^-1)b =
>> {
>> 0 => {0:1.2727922061357855}
>> 1 => {0:11.547005383792511}
>> 2 => {}
>> }
>> X =
>> {
>> 0 => {0:0.18}
>> 1 => {0:5.119661282874144}
>> 2 => {}
>> }
>> AX - B =
>> {
>> 0 => {0:459.0}
>> 1 => {0:670.0}
>> 2 => {0:881.0}
>> }
>>
>> org.scalatest.exceptions.TestFailedException was thrown.
>>
>>
>
Re: Cholesky test failures...
Posted by Dmitriy Lyubimov <dl...@gmail.com>.
On Mon, Jul 15, 2013 at 5:11 PM, Ted Dunning <te...@gmail.com> wrote:
> Does it help to know that the matrix A has rank 2 instead of 3?
>
yes, already figured that in my previous email, but full-rank doesn't work
with pivoted=true
>
> > chol(a)
> > Error in chol.default(a) :
> > the leading minor of order 3 is not positive definite
> > > qr.R(qr(a))
> > [,1] [,2] [,3]
> > [1,] -35.66511 -51.8153470 -6.796559e+01
> > [2,] 0.00000 -0.4120817 -8.241634e-01
> > [3,] 0.00000 0.0000000 4.218847e-15
> > > svd(a)$d
> > [1] 9.261128e+01 3.887216e-01 5.102882e-15
>
>
> I don't remember if I implemented a pivoting Cholesky or not.
>
ok so pivoted=true is the default and if it doesn't work, it probably
shouldn't be the default..
>
>
> On Mon, Jul 15, 2013 at 5:03 PM, Dmitriy Lyubimov <dl...@gmail.com>
> wrote:
>
> > should read
> >
> > val axmb = (a %*% x) - b
> >
> > of course but it doesn't make difference, the Cholesky output doesn't
> make
> > sense to me even before that
> >
> >
> > On Mon, Jul 15, 2013 at 5:00 PM, Dmitriy Lyubimov <dl...@gmail.com>
> > wrote:
> >
> > > Hi Ted,
> > >
> > > I am getting Cholesky test failures when trying to solve of Ax=B
> > >
> > > The L matrix and solveLeft() output do not make much sense to me. For
> > > once, L doesn't even have the expected L-shape:
> > >
> > > Do you have an idea where i go wrong? (the test is wrapped into scala
> DSL
> > > but it is Mahout's cholesky underwraps) .
> > >
> > > test code:
> > >
> > > test("chol") {
> > >
> > > // try to solve Ax=b with cholesky:
> > > // this requires
> > > // (LL')x = B
> > > // L'x= (L^-1)B
> > > // x=(L'^-1)(L^-1)B
> > >
> > > val a = dense((1, 2, 3), (2, 3, 4), (3, 4, 5))
> > >
> > > // make sure it is symmetric for a valid solution
> > > a := a.t %*% a
> > >
> > > printf("A= \n%s\n", a)
> > >
> > > val b = dense((9, 8, 7)).t
> > >
> > > printf ("b = \n%s\n", b)
> > >
> > > val ch = chol(a)
> > >
> > > printf ("L = \n%s\n", ch.getL)
> > >
> > > printf ( "(L^-1)b =\n%s\n", ch.solveLeft(b))
> > >
> > > val x = ch.solveRight(diag(1,3)) %*% ch.solveLeft(b)
> > >
> > > printf("x = \n%s\n", x.toString)
> > >
> > > val axmb = (a %*% b) - b
> > >
> > > printf("AX - B = \n%s\n", axmb.toString)
> > >
> > > assert(axmb.norm < 1e-10)
> > >
> > > }
> > >
> > >
> > >
> > > Output:
> > >
> > > A=
> > > {
> > > 0 => {0:14.0,1:20.0,2:26.0}
> > > 1 => {0:20.0,1:29.0,2:38.0}
> > > 2 => {0:26.0,1:38.0,2:50.0}
> > > }
> > > b =
> > > {
> > > 0 => {0:9.0}
> > > 1 => {0:8.0}
> > > 2 => {0:7.0}
> > > }
> > >
> > > L =
> > > {
> > > 0 => {0:0.6928203230275511,2:3.676955262170047}
> > > 1 => {0:0.3464101615137781,2:5.374011537017761}
> > > 2 => {2:7.0710678118654755}
> > > }
> > > (L^-1)b =
> > > {
> > > 0 => {0:1.2727922061357855}
> > > 1 => {0:11.547005383792511}
> > > 2 => {}
> > > }
> > > X =
> > > {
> > > 0 => {0:0.18}
> > > 1 => {0:5.119661282874144}
> > > 2 => {}
> > > }
> > > AX - B =
> > > {
> > > 0 => {0:459.0}
> > > 1 => {0:670.0}
> > > 2 => {0:881.0}
> > > }
> > >
> > > org.scalatest.exceptions.TestFailedException was thrown.
> > >
> > >
> >
>
Re: Cholesky test failures...
Posted by Ted Dunning <te...@gmail.com>.
Does it help to know that the matrix A has rank 2 instead of 3?
> chol(a)
> Error in chol.default(a) :
> the leading minor of order 3 is not positive definite
> > qr.R(qr(a))
> [,1] [,2] [,3]
> [1,] -35.66511 -51.8153470 -6.796559e+01
> [2,] 0.00000 -0.4120817 -8.241634e-01
> [3,] 0.00000 0.0000000 4.218847e-15
> > svd(a)$d
> [1] 9.261128e+01 3.887216e-01 5.102882e-15
I don't remember if I implemented a pivoting Cholesky or not.
On Mon, Jul 15, 2013 at 5:03 PM, Dmitriy Lyubimov <dl...@gmail.com> wrote:
> should read
>
> val axmb = (a %*% x) - b
>
> of course but it doesn't make difference, the Cholesky output doesn't make
> sense to me even before that
>
>
> On Mon, Jul 15, 2013 at 5:00 PM, Dmitriy Lyubimov <dl...@gmail.com>
> wrote:
>
> > Hi Ted,
> >
> > I am getting Cholesky test failures when trying to solve of Ax=B
> >
> > The L matrix and solveLeft() output do not make much sense to me. For
> > once, L doesn't even have the expected L-shape:
> >
> > Do you have an idea where i go wrong? (the test is wrapped into scala DSL
> > but it is Mahout's cholesky underwraps) .
> >
> > test code:
> >
> > test("chol") {
> >
> > // try to solve Ax=b with cholesky:
> > // this requires
> > // (LL')x = B
> > // L'x= (L^-1)B
> > // x=(L'^-1)(L^-1)B
> >
> > val a = dense((1, 2, 3), (2, 3, 4), (3, 4, 5))
> >
> > // make sure it is symmetric for a valid solution
> > a := a.t %*% a
> >
> > printf("A= \n%s\n", a)
> >
> > val b = dense((9, 8, 7)).t
> >
> > printf ("b = \n%s\n", b)
> >
> > val ch = chol(a)
> >
> > printf ("L = \n%s\n", ch.getL)
> >
> > printf ( "(L^-1)b =\n%s\n", ch.solveLeft(b))
> >
> > val x = ch.solveRight(diag(1,3)) %*% ch.solveLeft(b)
> >
> > printf("x = \n%s\n", x.toString)
> >
> > val axmb = (a %*% b) - b
> >
> > printf("AX - B = \n%s\n", axmb.toString)
> >
> > assert(axmb.norm < 1e-10)
> >
> > }
> >
> >
> >
> > Output:
> >
> > A=
> > {
> > 0 => {0:14.0,1:20.0,2:26.0}
> > 1 => {0:20.0,1:29.0,2:38.0}
> > 2 => {0:26.0,1:38.0,2:50.0}
> > }
> > b =
> > {
> > 0 => {0:9.0}
> > 1 => {0:8.0}
> > 2 => {0:7.0}
> > }
> >
> > L =
> > {
> > 0 => {0:0.6928203230275511,2:3.676955262170047}
> > 1 => {0:0.3464101615137781,2:5.374011537017761}
> > 2 => {2:7.0710678118654755}
> > }
> > (L^-1)b =
> > {
> > 0 => {0:1.2727922061357855}
> > 1 => {0:11.547005383792511}
> > 2 => {}
> > }
> > X =
> > {
> > 0 => {0:0.18}
> > 1 => {0:5.119661282874144}
> > 2 => {}
> > }
> > AX - B =
> > {
> > 0 => {0:459.0}
> > 1 => {0:670.0}
> > 2 => {0:881.0}
> > }
> >
> > org.scalatest.exceptions.TestFailedException was thrown.
> >
> >
>
Re: Cholesky test failures...
Posted by Dmitriy Lyubimov <dl...@gmail.com>.
should read
val axmb = (a %*% x) - b
of course but it doesn't make difference, the Cholesky output doesn't make
sense to me even before that
On Mon, Jul 15, 2013 at 5:00 PM, Dmitriy Lyubimov <dl...@gmail.com> wrote:
> Hi Ted,
>
> I am getting Cholesky test failures when trying to solve of Ax=B
>
> The L matrix and solveLeft() output do not make much sense to me. For
> once, L doesn't even have the expected L-shape:
>
> Do you have an idea where i go wrong? (the test is wrapped into scala DSL
> but it is Mahout's cholesky underwraps) .
>
> test code:
>
> test("chol") {
>
> // try to solve Ax=b with cholesky:
> // this requires
> // (LL')x = B
> // L'x= (L^-1)B
> // x=(L'^-1)(L^-1)B
>
> val a = dense((1, 2, 3), (2, 3, 4), (3, 4, 5))
>
> // make sure it is symmetric for a valid solution
> a := a.t %*% a
>
> printf("A= \n%s\n", a)
>
> val b = dense((9, 8, 7)).t
>
> printf ("b = \n%s\n", b)
>
> val ch = chol(a)
>
> printf ("L = \n%s\n", ch.getL)
>
> printf ( "(L^-1)b =\n%s\n", ch.solveLeft(b))
>
> val x = ch.solveRight(diag(1,3)) %*% ch.solveLeft(b)
>
> printf("x = \n%s\n", x.toString)
>
> val axmb = (a %*% b) - b
>
> printf("AX - B = \n%s\n", axmb.toString)
>
> assert(axmb.norm < 1e-10)
>
> }
>
>
>
> Output:
>
> A=
> {
> 0 => {0:14.0,1:20.0,2:26.0}
> 1 => {0:20.0,1:29.0,2:38.0}
> 2 => {0:26.0,1:38.0,2:50.0}
> }
> b =
> {
> 0 => {0:9.0}
> 1 => {0:8.0}
> 2 => {0:7.0}
> }
>
> L =
> {
> 0 => {0:0.6928203230275511,2:3.676955262170047}
> 1 => {0:0.3464101615137781,2:5.374011537017761}
> 2 => {2:7.0710678118654755}
> }
> (L^-1)b =
> {
> 0 => {0:1.2727922061357855}
> 1 => {0:11.547005383792511}
> 2 => {}
> }
> X =
> {
> 0 => {0:0.18}
> 1 => {0:5.119661282874144}
> 2 => {}
> }
> AX - B =
> {
> 0 => {0:459.0}
> 1 => {0:670.0}
> 2 => {0:881.0}
> }
>
> org.scalatest.exceptions.TestFailedException was thrown.
>
>