+++ /dev/null
-// 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"
- "testing"
-
- "golang.org/x/exp/rand"
-
- "gonum.org/v1/gonum/blas/testblas"
- "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)
- }
- U := chol.UTo(nil)
- aCopy := DenseCopyOf(test.a)
- var a Dense
- a.Mul(U.TTri(), U)
- if !EqualApprox(&a, aCopy, 1e-14) {
- t.Error("unexpected Cholesky factor product")
- }
-
- L := chol.LTo(nil)
- a.Mul(L, L.TTri())
- if !EqualApprox(&a, aCopy, 1e-14) {
- t.Error("unexpected Cholesky factor product")
- }
- }
- }
-}
-
-func TestCholeskySolve(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
- chol.Solve(&x, test.b)
- 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 TestCholeskySolveChol(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
- chola.SolveChol(&x, &cholb)
-
- var ans Dense
- ans.Mul(test.a, &x)
- if !EqualApprox(&ans, test.b, 1e-12) {
- var y Dense
- y.Solve(test.a, test.b)
- t.Errorf("incorrect Cholesky solve solution product\ngot solution:\n%.4v\nwant solution\n%.4v",
- Formatted(&x), Formatted(&y))
- }
- }
-}
-
-func TestCholeskySolveVec(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
- chol.SolveVec(&x, test.b)
- 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")
- }
- s := chol.ToSym(nil)
-
- 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 zero")
- }
- if !Equal(chol.chol, chol2.chol) {
- t.Errorf("chol mismatch from zero")
- }
-
- // 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-zero")
- }
- if !Equal(chol.chol, chol2.chol) {
- t.Errorf("chol mismatch from non-zero")
- }
- }
-}
-
-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
- chol.InverseTo(&sInv)
-
- 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 < 10; k++ {
- data := make([]float64, n*n)
- for i := range data {
- data[i] = rand.NormFloat64()
- }
-
- var a SymDense
- a.SymOuterK(1, NewDense(n, n, data))
-
- xdata := make([]float64, n)
- for i := range xdata {
- xdata[i] = rand.NormFloat64()
- }
- x := NewVecDense(n, xdata)
-
- var chol Cholesky
- ok := chol.Factorize(&a)
- if !ok {
- t.Errorf("Bad random test, Cholesky factorization failed")
- continue
- }
-
- alpha := rand.Float64()
- ok = chol.SymRankOne(&chol, alpha, x)
- if !ok {
- t.Errorf("n=%v, alpha=%v: unexpected failure", n, alpha)
- continue
- }
- a.SymRankOne(&a, alpha, x)
-
- var achol SymDense
- chol.ToSym(&achol)
- if !EqualApprox(&achol, &a, 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(&a), Formatted(&achol))
- }
- }
- }
-
- 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,
- },
- } {
- 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.SliceSquare(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)
- }
- a := cholNew.ToSym(nil)
- 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 BenchmarkCholeskySmall(b *testing.B) {
- benchmarkCholesky(b, 2)
-}
-
-func BenchmarkCholeskyMedium(b *testing.B) {
- benchmarkCholesky(b, testblas.MediumMat)
-}
-
-func BenchmarkCholeskyLarge(b *testing.B) {
- benchmarkCholesky(b, testblas.LargeMat)
-}
-
-func benchmarkCholesky(b *testing.B, n int) {
- base := make([]float64, n*n)
- for i := range base {
- base[i] = rand.Float64()
- }
- bm := NewDense(n, n, base)
- bm.Mul(bm.T(), bm)
- am := NewSymDense(n, bm.mat.Data)
-
- var chol Cholesky
- b.ResetTimer()
- for i := 0; i < b.N; i++ {
- ok := chol.Factorize(am)
- if !ok {
- panic("not pos def")
- }
- }
-}
OSDN Git Service
