On Mon, Jul 15, 2013 at 5:11 PM, Ted Dunning <[email protected]> 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 <[email protected]>
> 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 <[email protected]>
> > 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.
> > >
> > >
> >
>

Reply via email to