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 <[email protected]> wrote:
> sorry, i mean ch.solveRight(eye(4)) fails with AOOB
>
>
> On Wed, Jul 17, 2013 at 3:11 PM, Dmitriy Lyubimov <[email protected]>
> 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 <[email protected]
> >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 <[email protected]>
> >> 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
> >> > :)
> >> >
> >>
> >
> >
>
