blob: 5a090b1ddbb29e76a5ac80a0077df918b6c3d781 [file] [log] [blame]
 // Copyright ©2013 The Gonum Authors. All rights reserved. // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. package mat import ( "math" "strconv" "testing" "golang.org/x/exp/rand" "gonum.org/v1/gonum/floats" ) func TestCholesky(t *testing.T) { for _, test := range []struct { a *SymDense cond float64 want *TriDense posdef bool }{ { a: NewSymDense(3, []float64{ 4, 1, 1, 0, 2, 3, 0, 0, 6, }), cond: 37, want: NewTriDense(3, true, []float64{ 2, 0.5, 0.5, 0, 1.3228756555322954, 2.0788046015507495, 0, 0, 1.195228609334394, }), posdef: true, }, } { _, n := test.a.Dims() for _, chol := range []*Cholesky{ {}, {chol: NewTriDense(n-1, true, nil)}, {chol: NewTriDense(n, true, nil)}, {chol: NewTriDense(n+1, true, nil)}, } { ok := chol.Factorize(test.a) if ok != test.posdef { t.Errorf("unexpected return from Cholesky factorization: got: ok=%t want: ok=%t", ok, test.posdef) } fc := DenseCopyOf(chol.chol) if !Equal(fc, test.want) { t.Error("incorrect Cholesky factorization") } if math.Abs(test.cond-chol.cond) > 1e-13 { t.Errorf("Condition number mismatch: Want %v, got %v", test.cond, chol.cond) } var U TriDense chol.UTo(&U) aCopy := DenseCopyOf(test.a) var a Dense a.Mul(U.TTri(), &U) if !EqualApprox(&a, aCopy, 1e-14) { t.Error("unexpected Cholesky factor product") } var L TriDense chol.LTo(&L) a.Mul(&L, L.TTri()) if !EqualApprox(&a, aCopy, 1e-14) { t.Error("unexpected Cholesky factor product") } } } } func TestCholeskyAt(t *testing.T) { for _, test := range []*SymDense{ NewSymDense(3, []float64{ 53, 59, 37, 59, 83, 71, 37, 71, 101, }), } { var chol Cholesky ok := chol.Factorize(test) if !ok { t.Fatalf("Matrix not positive definite") } n := test.Symmetric() cn := chol.Symmetric() if cn != n { t.Errorf("Cholesky size does not match. Got %d, want %d", cn, n) } for i := 0; i < n; i++ { for j := 0; j < n; j++ { got := chol.At(i, j) want := test.At(i, j) if math.Abs(got-want) > 1e-12 { t.Errorf("Cholesky at does not match at %d, %d. Got %v, want %v", i, j, got, want) } } } } } func TestCholeskySolveTo(t *testing.T) { for _, test := range []struct { a *SymDense b *Dense ans *Dense }{ { a: NewSymDense(2, []float64{ 1, 0, 0, 1, }), b: NewDense(2, 1, []float64{5, 6}), ans: NewDense(2, 1, []float64{5, 6}), }, { a: NewSymDense(3, []float64{ 53, 59, 37, 0, 83, 71, 37, 71, 101, }), b: NewDense(3, 1, []float64{5, 6, 7}), ans: NewDense(3, 1, []float64{0.20745069393718094, -0.17421475529583694, 0.11577794010226464}), }, } { var chol Cholesky ok := chol.Factorize(test.a) if !ok { t.Fatal("unexpected Cholesky factorization failure: not positive definite") } var x Dense err := chol.SolveTo(&x, test.b) if err != nil { t.Errorf("unexpected error from Cholesky solve: %v", err) } if !EqualApprox(&x, test.ans, 1e-12) { t.Error("incorrect Cholesky solve solution") } var ans Dense ans.Mul(test.a, &x) if !EqualApprox(&ans, test.b, 1e-12) { t.Error("incorrect Cholesky solve solution product") } } } func TestCholeskySolveCholTo(t *testing.T) { for _, test := range []struct { a, b *SymDense }{ { a: NewSymDense(2, []float64{ 1, 0, 0, 1, }), b: NewSymDense(2, []float64{ 1, 0, 0, 1, }), }, { a: NewSymDense(2, []float64{ 1, 0, 0, 1, }), b: NewSymDense(2, []float64{ 2, 0, 0, 2, }), }, { a: NewSymDense(3, []float64{ 53, 59, 37, 59, 83, 71, 37, 71, 101, }), b: NewSymDense(3, []float64{ 2, -1, 0, -1, 2, -1, 0, -1, 2, }), }, } { var chola, cholb Cholesky ok := chola.Factorize(test.a) if !ok { t.Fatal("unexpected Cholesky factorization failure for a: not positive definite") } ok = cholb.Factorize(test.b) if !ok { t.Fatal("unexpected Cholesky factorization failure for b: not positive definite") } var x Dense err := chola.SolveCholTo(&x, &cholb) if err != nil { t.Errorf("unexpected error from Cholesky solve: %v", err) } var ans Dense ans.Mul(test.a, &x) if !EqualApprox(&ans, test.b, 1e-12) { var y Dense err := y.Solve(test.a, test.b) if err != nil { t.Errorf("unexpected error from dense solve: %v", err) } t.Errorf("incorrect Cholesky solve solution product\ngot solution:\n%.4v\nwant solution\n%.4v", Formatted(&x), Formatted(&y)) } } } func TestCholeskySolveVecTo(t *testing.T) { for _, test := range []struct { a *SymDense b *VecDense ans *VecDense }{ { a: NewSymDense(2, []float64{ 1, 0, 0, 1, }), b: NewVecDense(2, []float64{5, 6}), ans: NewVecDense(2, []float64{5, 6}), }, { a: NewSymDense(3, []float64{ 53, 59, 37, 0, 83, 71, 0, 0, 101, }), b: NewVecDense(3, []float64{5, 6, 7}), ans: NewVecDense(3, []float64{0.20745069393718094, -0.17421475529583694, 0.11577794010226464}), }, } { var chol Cholesky ok := chol.Factorize(test.a) if !ok { t.Fatal("unexpected Cholesky factorization failure: not positive definite") } var x VecDense err := chol.SolveVecTo(&x, test.b) if err != nil { t.Errorf("unexpected error from Cholesky solve: %v", err) } if !EqualApprox(&x, test.ans, 1e-12) { t.Error("incorrect Cholesky solve solution") } var ans VecDense ans.MulVec(test.a, &x) if !EqualApprox(&ans, test.b, 1e-12) { t.Error("incorrect Cholesky solve solution product") } } } func TestCholeskyToSym(t *testing.T) { for _, test := range []*SymDense{ NewSymDense(3, []float64{ 53, 59, 37, 0, 83, 71, 0, 0, 101, }), } { var chol Cholesky ok := chol.Factorize(test) if !ok { t.Fatal("unexpected Cholesky factorization failure: not positive definite") } var s SymDense chol.ToSym(&s) if !EqualApprox(&s, test, 1e-12) { t.Errorf("Cholesky reconstruction not equal to original matrix.\nWant:\n% v\nGot:\n% v\n", Formatted(test), Formatted(&s)) } } } func TestCloneCholesky(t *testing.T) { for _, test := range []*SymDense{ NewSymDense(3, []float64{ 53, 59, 37, 0, 83, 71, 0, 0, 101, }), } { var chol Cholesky ok := chol.Factorize(test) if !ok { panic("bad test") } var chol2 Cholesky chol2.Clone(&chol) if chol.cond != chol2.cond { t.Errorf("condition number mismatch from empty") } if !Equal(chol.chol, chol2.chol) { t.Errorf("chol mismatch from empty") } // Corrupt chol2 and try again chol2.cond = math.NaN() chol2.chol = NewTriDense(2, Upper, nil) chol2.Clone(&chol) if chol.cond != chol2.cond { t.Errorf("condition number mismatch from non-empty") } if !Equal(chol.chol, chol2.chol) { t.Errorf("chol mismatch from non-empty") } } } func TestCholeskyInverseTo(t *testing.T) { for _, n := range []int{1, 3, 5, 9} { data := make([]float64, n*n) for i := range data { data[i] = rand.NormFloat64() } var s SymDense s.SymOuterK(1, NewDense(n, n, data)) var chol Cholesky ok := chol.Factorize(&s) if !ok { t.Errorf("Bad test, cholesky decomposition failed") } var sInv SymDense err := chol.InverseTo(&sInv) if err != nil { t.Errorf("unexpected error from Cholesky inverse: %v", err) } var ans Dense ans.Mul(&sInv, &s) if !equalApprox(eye(n), &ans, 1e-8, false) { var diff Dense diff.Sub(eye(n), &ans) t.Errorf("SymDense times Cholesky inverse not identity. Norm diff = %v", Norm(&diff, 2)) } } } func TestCholeskySymRankOne(t *testing.T) { rand.Seed(1) for _, n := range []int{1, 2, 3, 4, 5, 7, 10, 20, 50, 100} { for k := 0; k < 50; k++ { // Construct a random positive definite matrix. data := make([]float64, n*n) for i := range data { data[i] = rand.NormFloat64() } var a SymDense a.SymOuterK(1, NewDense(n, n, data)) // Construct random data for updating. xdata := make([]float64, n) for i := range xdata { xdata[i] = rand.NormFloat64() } x := NewVecDense(n, xdata) alpha := rand.NormFloat64() // Compute the updated matrix directly. If alpha > 0, there are no // issues. If alpha < 0, it could be that the final matrix is not // positive definite, so instead switch the two matrices. aUpdate := NewSymDense(n, nil) if alpha > 0 { aUpdate.SymRankOne(&a, alpha, x) } else { aUpdate.CopySym(&a) a.Reset() a.SymRankOne(aUpdate, -alpha, x) } // Compare the Cholesky decomposition computed with Cholesky.SymRankOne // with that computed from updating A directly. var chol Cholesky ok := chol.Factorize(&a) if !ok { t.Errorf("Bad random test, Cholesky factorization failed") continue } var cholUpdate Cholesky ok = cholUpdate.SymRankOne(&chol, alpha, x) if !ok { t.Errorf("n=%v, alpha=%v: unexpected failure", n, alpha) continue } var aCompare SymDense cholUpdate.ToSym(&aCompare) if !EqualApprox(&aCompare, aUpdate, 1e-13) { t.Errorf("n=%v, alpha=%v: mismatch between updated matrix and from Cholesky:\nupdated:\n%v\nfrom Cholesky:\n%v", n, alpha, Formatted(aUpdate), Formatted(&aCompare)) } } } for i, test := range []struct { a *SymDense alpha float64 x []float64 wantOk bool }{ { // Update (to positive definite matrix). a: NewSymDense(4, []float64{ 1, 1, 1, 1, 0, 2, 3, 4, 0, 0, 6, 10, 0, 0, 0, 20, }), alpha: 1, x: []float64{0, 0, 0, 1}, wantOk: true, }, { // Downdate to singular matrix. a: NewSymDense(4, []float64{ 1, 1, 1, 1, 0, 2, 3, 4, 0, 0, 6, 10, 0, 0, 0, 20, }), alpha: -1, x: []float64{0, 0, 0, 1}, wantOk: false, }, { // Downdate to positive definite matrix. a: NewSymDense(4, []float64{ 1, 1, 1, 1, 0, 2, 3, 4, 0, 0, 6, 10, 0, 0, 0, 20, }), alpha: -1 / 2, x: []float64{0, 0, 0, 1}, wantOk: true, }, { // Issue #453. a: NewSymDense(1, []float64{1}), alpha: -1, x: []float64{0.25}, wantOk: true, }, } { var chol Cholesky ok := chol.Factorize(test.a) if !ok { t.Errorf("Case %v: bad test, Cholesky factorization failed", i) continue } x := NewVecDense(len(test.x), test.x) ok = chol.SymRankOne(&chol, test.alpha, x) if !ok { if test.wantOk { t.Errorf("Case %v: unexpected failure from SymRankOne", i) } continue } if ok && !test.wantOk { t.Errorf("Case %v: expected a failure from SymRankOne", i) } a := test.a a.SymRankOne(a, test.alpha, x) var achol SymDense chol.ToSym(&achol) if !EqualApprox(&achol, a, 1e-13) { t.Errorf("Case %v: mismatch between updated matrix and from Cholesky:\nupdated:\n%v\nfrom Cholesky:\n%v", i, Formatted(a), Formatted(&achol)) } } } func TestCholeskyExtendVecSym(t *testing.T) { for cas, test := range []struct { a *SymDense }{ { a: NewSymDense(3, []float64{ 4, 1, 1, 0, 2, 3, 0, 0, 6, }), }, } { n := test.a.Symmetric() as := test.a.SliceSym(0, n-1).(*SymDense) // Compute the full factorization to use later (do the full factorization // first to ensure the matrix is positive definite). var cholFull Cholesky ok := cholFull.Factorize(test.a) if !ok { panic("mat: bad test, matrix not positive definite") } var chol Cholesky ok = chol.Factorize(as) if !ok { panic("mat: bad test, subset is not positive definite") } row := NewVecDense(n, nil) for i := 0; i < n; i++ { row.SetVec(i, test.a.At(n-1, i)) } var cholNew Cholesky ok = cholNew.ExtendVecSym(&chol, row) if !ok { t.Errorf("cas %v: update not positive definite", cas) } var a SymDense cholNew.ToSym(&a) if !EqualApprox(&a, test.a, 1e-12) { t.Errorf("cas %v: mismatch", cas) } // test in-place ok = chol.ExtendVecSym(&chol, row) if !ok { t.Errorf("cas %v: in-place update not positive definite", cas) } if !equalChol(&chol, &cholNew) { t.Errorf("cas %v: Cholesky different in-place vs. new", cas) } // Test that the factorization is about right compared with the direct // full factorization. Use a high tolerance on the condition number // since the condition number with the updated rule is approximate. if !equalApproxChol(&chol, &cholFull, 1e-12, 0.3) { t.Errorf("cas %v: updated Cholesky does not match full", cas) } } } func TestCholeskyScale(t *testing.T) { for cas, test := range []struct { a *SymDense f float64 }{ { a: NewSymDense(3, []float64{ 4, 1, 1, 0, 2, 3, 0, 0, 6, }), f: 0.5, }, } { var chol Cholesky ok := chol.Factorize(test.a) if !ok { t.Errorf("Case %v: bad test, Cholesky factorization failed", cas) continue } // Compare the update to a new Cholesky to an update in-place. var cholUpdate Cholesky cholUpdate.Scale(test.f, &chol) chol.Scale(test.f, &chol) if !equalChol(&chol, &cholUpdate) { t.Errorf("Case %d: cholesky mismatch new receiver", cas) } var sym SymDense chol.ToSym(&sym) var comp SymDense comp.ScaleSym(test.f, test.a) if !EqualApprox(&comp, &sym, 1e-14) { t.Errorf("Case %d: cholesky reconstruction doesn't match scaled matrix", cas) } var cholTest Cholesky cholTest.Factorize(&comp) if !equalApproxChol(&cholTest, &chol, 1e-12, 1e-12) { t.Errorf("Case %d: cholesky mismatch with scaled matrix. %v, %v", cas, cholTest.cond, chol.cond) } } } // equalApproxChol checks that the two Cholesky decompositions are equal. func equalChol(a, b *Cholesky) bool { return Equal(a.chol, b.chol) && a.cond == b.cond } // equalApproxChol checks that the two Cholesky decompositions are approximately // the same with the given tolerance on equality for the Triangular component and // condition. func equalApproxChol(a, b *Cholesky, matTol, condTol float64) bool { if !EqualApprox(a.chol, b.chol, matTol) { return false } return floats.EqualWithinAbsOrRel(a.cond, b.cond, condTol, condTol) } func BenchmarkCholeskyFactorize(b *testing.B) { for _, n := range []int{10, 100, 1000} { b.Run("n="+strconv.Itoa(n), func(b *testing.B) { rnd := rand.New(rand.NewSource(1)) data := make([]float64, n*n) for i := range data { data[i] = rnd.NormFloat64() } var a SymDense a.SymOuterK(1, NewDense(n, n, data)) var chol Cholesky b.ResetTimer() for i := 0; i < b.N; i++ { ok := chol.Factorize(&a) if !ok { panic("not positive definite") } } }) } } func BenchmarkCholeskyToSym(b *testing.B) { for _, n := range []int{10, 100, 1000} { b.Run("n="+strconv.Itoa(n), func(b *testing.B) { rnd := rand.New(rand.NewSource(1)) data := make([]float64, n*n) for i := range data { data[i] = rnd.NormFloat64() } var a SymDense a.SymOuterK(1, NewDense(n, n, data)) var chol Cholesky ok := chol.Factorize(&a) if !ok { panic("not positive definite") } dst := NewSymDense(n, nil) b.ResetTimer() for i := 0; i < b.N; i++ { chol.ToSym(dst) } }) } } func BenchmarkCholeskyInverseTo(b *testing.B) { for _, n := range []int{10, 100, 1000} { b.Run("n="+strconv.Itoa(n), func(b *testing.B) { rnd := rand.New(rand.NewSource(1)) data := make([]float64, n*n) for i := range data { data[i] = rnd.NormFloat64() } var a SymDense a.SymOuterK(1, NewDense(n, n, data)) var chol Cholesky ok := chol.Factorize(&a) if !ok { panic("not positive definite") } dst := NewSymDense(n, nil) b.ResetTimer() for i := 0; i < b.N; i++ { err := chol.InverseTo(dst) if err != nil { b.Fatalf("unexpected error from Cholesky inverse: %v", err) } } }) } }