1 // Copyright ©2013 The Gonum Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
11 "golang.org/x/exp/rand"
13 "gonum.org/v1/gonum/blas/testblas"
14 "gonum.org/v1/gonum/floats"
17 func TestCholesky(t *testing.T) {
18 for _, test := range []struct {
26 a: NewSymDense(3, []float64{
32 want: NewTriDense(3, true, []float64{
34 0, 1.3228756555322954, 2.0788046015507495,
35 0, 0, 1.195228609334394,
41 for _, chol := range []*Cholesky{
43 {chol: NewTriDense(n-1, true, nil)},
44 {chol: NewTriDense(n, true, nil)},
45 {chol: NewTriDense(n+1, true, nil)},
47 ok := chol.Factorize(test.a)
48 if ok != test.posdef {
49 t.Errorf("unexpected return from Cholesky factorization: got: ok=%t want: ok=%t", ok, test.posdef)
51 fc := DenseCopyOf(chol.chol)
52 if !Equal(fc, test.want) {
53 t.Error("incorrect Cholesky factorization")
55 if math.Abs(test.cond-chol.cond) > 1e-13 {
56 t.Errorf("Condition number mismatch: Want %v, got %v", test.cond, chol.cond)
59 aCopy := DenseCopyOf(test.a)
62 if !EqualApprox(&a, aCopy, 1e-14) {
63 t.Error("unexpected Cholesky factor product")
68 if !EqualApprox(&a, aCopy, 1e-14) {
69 t.Error("unexpected Cholesky factor product")
75 func TestCholeskySolve(t *testing.T) {
76 for _, test := range []struct {
82 a: NewSymDense(2, []float64{
86 b: NewDense(2, 1, []float64{5, 6}),
87 ans: NewDense(2, 1, []float64{5, 6}),
90 a: NewSymDense(3, []float64{
95 b: NewDense(3, 1, []float64{5, 6, 7}),
96 ans: NewDense(3, 1, []float64{0.20745069393718094, -0.17421475529583694, 0.11577794010226464}),
100 ok := chol.Factorize(test.a)
102 t.Fatal("unexpected Cholesky factorization failure: not positive definite")
106 chol.Solve(&x, test.b)
107 if !EqualApprox(&x, test.ans, 1e-12) {
108 t.Error("incorrect Cholesky solve solution")
113 if !EqualApprox(&ans, test.b, 1e-12) {
114 t.Error("incorrect Cholesky solve solution product")
119 func TestCholeskySolveChol(t *testing.T) {
120 for _, test := range []struct {
124 a: NewSymDense(2, []float64{
128 b: NewSymDense(2, []float64{
134 a: NewSymDense(2, []float64{
138 b: NewSymDense(2, []float64{
144 a: NewSymDense(3, []float64{
149 b: NewSymDense(3, []float64{
156 var chola, cholb Cholesky
157 ok := chola.Factorize(test.a)
159 t.Fatal("unexpected Cholesky factorization failure for a: not positive definite")
161 ok = cholb.Factorize(test.b)
163 t.Fatal("unexpected Cholesky factorization failure for b: not positive definite")
167 chola.SolveChol(&x, &cholb)
171 if !EqualApprox(&ans, test.b, 1e-12) {
173 y.Solve(test.a, test.b)
174 t.Errorf("incorrect Cholesky solve solution product\ngot solution:\n%.4v\nwant solution\n%.4v",
175 Formatted(&x), Formatted(&y))
180 func TestCholeskySolveVec(t *testing.T) {
181 for _, test := range []struct {
187 a: NewSymDense(2, []float64{
191 b: NewVecDense(2, []float64{5, 6}),
192 ans: NewVecDense(2, []float64{5, 6}),
195 a: NewSymDense(3, []float64{
200 b: NewVecDense(3, []float64{5, 6, 7}),
201 ans: NewVecDense(3, []float64{0.20745069393718094, -0.17421475529583694, 0.11577794010226464}),
205 ok := chol.Factorize(test.a)
207 t.Fatal("unexpected Cholesky factorization failure: not positive definite")
211 chol.SolveVec(&x, test.b)
212 if !EqualApprox(&x, test.ans, 1e-12) {
213 t.Error("incorrect Cholesky solve solution")
217 ans.MulVec(test.a, &x)
218 if !EqualApprox(&ans, test.b, 1e-12) {
219 t.Error("incorrect Cholesky solve solution product")
224 func TestCholeskyToSym(t *testing.T) {
225 for _, test := range []*SymDense{
226 NewSymDense(3, []float64{
233 ok := chol.Factorize(test)
235 t.Fatal("unexpected Cholesky factorization failure: not positive definite")
239 if !EqualApprox(s, test, 1e-12) {
240 t.Errorf("Cholesky reconstruction not equal to original matrix.\nWant:\n% v\nGot:\n% v\n", Formatted(test), Formatted(s))
245 func TestCloneCholesky(t *testing.T) {
246 for _, test := range []*SymDense{
247 NewSymDense(3, []float64{
254 ok := chol.Factorize(test)
261 if chol.cond != chol2.cond {
262 t.Errorf("condition number mismatch from zero")
264 if !Equal(chol.chol, chol2.chol) {
265 t.Errorf("chol mismatch from zero")
268 // Corrupt chol2 and try again
269 chol2.cond = math.NaN()
270 chol2.chol = NewTriDense(2, Upper, nil)
272 if chol.cond != chol2.cond {
273 t.Errorf("condition number mismatch from non-zero")
275 if !Equal(chol.chol, chol2.chol) {
276 t.Errorf("chol mismatch from non-zero")
281 func TestCholeskyInverseTo(t *testing.T) {
282 for _, n := range []int{1, 3, 5, 9} {
283 data := make([]float64, n*n)
284 for i := range data {
285 data[i] = rand.NormFloat64()
288 s.SymOuterK(1, NewDense(n, n, data))
291 ok := chol.Factorize(&s)
293 t.Errorf("Bad test, cholesky decomposition failed")
297 chol.InverseTo(&sInv)
301 if !equalApprox(eye(n), &ans, 1e-8, false) {
303 diff.Sub(eye(n), &ans)
304 t.Errorf("SymDense times Cholesky inverse not identity. Norm diff = %v", Norm(&diff, 2))
309 func TestCholeskySymRankOne(t *testing.T) {
311 for _, n := range []int{1, 2, 3, 4, 5, 7, 10, 20, 50, 100} {
312 for k := 0; k < 10; k++ {
313 data := make([]float64, n*n)
314 for i := range data {
315 data[i] = rand.NormFloat64()
319 a.SymOuterK(1, NewDense(n, n, data))
321 xdata := make([]float64, n)
322 for i := range xdata {
323 xdata[i] = rand.NormFloat64()
325 x := NewVecDense(n, xdata)
328 ok := chol.Factorize(&a)
330 t.Errorf("Bad random test, Cholesky factorization failed")
334 alpha := rand.Float64()
335 ok = chol.SymRankOne(&chol, alpha, x)
337 t.Errorf("n=%v, alpha=%v: unexpected failure", n, alpha)
340 a.SymRankOne(&a, alpha, x)
344 if !EqualApprox(&achol, &a, 1e-13) {
345 t.Errorf("n=%v, alpha=%v: mismatch between updated matrix and from Cholesky:\nupdated:\n%v\nfrom Cholesky:\n%v",
346 n, alpha, Formatted(&a), Formatted(&achol))
351 for i, test := range []struct {
359 // Update (to positive definite matrix).
360 a: NewSymDense(4, []float64{
367 x: []float64{0, 0, 0, 1},
371 // Downdate to singular matrix.
372 a: NewSymDense(4, []float64{
379 x: []float64{0, 0, 0, 1},
383 // Downdate to positive definite matrix.
384 a: NewSymDense(4, []float64{
391 x: []float64{0, 0, 0, 1},
396 ok := chol.Factorize(test.a)
398 t.Errorf("Case %v: bad test, Cholesky factorization failed", i)
402 x := NewVecDense(len(test.x), test.x)
403 ok = chol.SymRankOne(&chol, test.alpha, x)
406 t.Errorf("Case %v: unexpected failure from SymRankOne", i)
410 if ok && !test.wantOk {
411 t.Errorf("Case %v: expected a failure from SymRankOne", i)
415 a.SymRankOne(a, test.alpha, x)
419 if !EqualApprox(&achol, a, 1e-13) {
420 t.Errorf("Case %v: mismatch between updated matrix and from Cholesky:\nupdated:\n%v\nfrom Cholesky:\n%v",
421 i, Formatted(a), Formatted(&achol))
426 func TestCholeskyExtendVecSym(t *testing.T) {
427 for cas, test := range []struct {
431 a: NewSymDense(3, []float64{
438 n := test.a.Symmetric()
439 as := test.a.SliceSquare(0, n-1).(*SymDense)
441 // Compute the full factorization to use later (do the full factorization
442 // first to ensure the matrix is positive definite).
443 var cholFull Cholesky
444 ok := cholFull.Factorize(test.a)
446 panic("mat: bad test, matrix not positive definite")
450 ok = chol.Factorize(as)
452 panic("mat: bad test, subset is not positive definite")
454 row := NewVecDense(n, nil)
455 for i := 0; i < n; i++ {
456 row.SetVec(i, test.a.At(n-1, i))
460 ok = cholNew.ExtendVecSym(&chol, row)
462 t.Errorf("cas %v: update not positive definite", cas)
464 a := cholNew.ToSym(nil)
465 if !EqualApprox(a, test.a, 1e-12) {
466 t.Errorf("cas %v: mismatch", cas)
470 ok = chol.ExtendVecSym(&chol, row)
472 t.Errorf("cas %v: in-place update not positive definite", cas)
474 if !equalChol(&chol, &cholNew) {
475 t.Errorf("cas %v: Cholesky different in-place vs. new", cas)
478 // Test that the factorization is about right compared with the direct
479 // full factorization. Use a high tolerance on the condition number
480 // since the condition number with the updated rule is approximate.
481 if !equalApproxChol(&chol, &cholFull, 1e-12, 0.3) {
482 t.Errorf("cas %v: updated Cholesky does not match full", cas)
487 func TestCholeskyScale(t *testing.T) {
488 for cas, test := range []struct {
493 a: NewSymDense(3, []float64{
502 ok := chol.Factorize(test.a)
504 t.Errorf("Case %v: bad test, Cholesky factorization failed", cas)
508 // Compare the update to a new Cholesky to an update in-place.
509 var cholUpdate Cholesky
510 cholUpdate.Scale(test.f, &chol)
511 chol.Scale(test.f, &chol)
512 if !equalChol(&chol, &cholUpdate) {
513 t.Errorf("Case %d: cholesky mismatch new receiver", cas)
518 comp.ScaleSym(test.f, test.a)
519 if !EqualApprox(&comp, &sym, 1e-14) {
520 t.Errorf("Case %d: cholesky reconstruction doesn't match scaled matrix", cas)
523 var cholTest Cholesky
524 cholTest.Factorize(&comp)
525 if !equalApproxChol(&cholTest, &chol, 1e-12, 1e-12) {
526 t.Errorf("Case %d: cholesky mismatch with scaled matrix. %v, %v", cas, cholTest.cond, chol.cond)
531 // equalApproxChol checks that the two Cholesky decompositions are equal.
532 func equalChol(a, b *Cholesky) bool {
533 return Equal(a.chol, b.chol) && a.cond == b.cond
536 // equalApproxChol checks that the two Cholesky decompositions are approximately
537 // the same with the given tolerance on equality for the Triangular component and
539 func equalApproxChol(a, b *Cholesky, matTol, condTol float64) bool {
540 if !EqualApprox(a.chol, b.chol, matTol) {
543 return floats.EqualWithinAbsOrRel(a.cond, b.cond, condTol, condTol)
546 func BenchmarkCholeskySmall(b *testing.B) {
547 benchmarkCholesky(b, 2)
550 func BenchmarkCholeskyMedium(b *testing.B) {
551 benchmarkCholesky(b, testblas.MediumMat)
554 func BenchmarkCholeskyLarge(b *testing.B) {
555 benchmarkCholesky(b, testblas.LargeMat)
558 func benchmarkCholesky(b *testing.B, n int) {
559 base := make([]float64, n*n)
560 for i := range base {
561 base[i] = rand.Float64()
563 bm := NewDense(n, n, base)
565 am := NewSymDense(n, bm.mat.Data)
569 for i := 0; i < b.N; i++ {
570 ok := chol.Factorize(am)