vendor/gonum.org/v1/gonum/mat/cholesky_test.go

   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.
   4
   5 package mat
   6
   7 import (
   8         "math"
   9         "testing"
  10
  11         "golang.org/x/exp/rand"
  12
  13         "gonum.org/v1/gonum/blas/testblas"
  14         "gonum.org/v1/gonum/floats"
  15 )
  16
  17 func TestCholesky(t *testing.T) {
  18         for _, test := range []struct {
  19                 a *SymDense
  20
  21                 cond   float64
  22                 want   *TriDense
  23                 posdef bool
  24         }{
  25                 {
  26                         a: NewSymDense(3, []float64{
  27                                 4, 1, 1,
  28                                 0, 2, 3,
  29                                 0, 0, 6,
  30                         }),
  31                         cond: 37,
  32                         want: NewTriDense(3, true, []float64{
  33                                 2, 0.5, 0.5,
  34                                 0, 1.3228756555322954, 2.0788046015507495,
  35                                 0, 0, 1.195228609334394,
  36                         }),
  37                         posdef: true,
  38                 },
  39         } {
  40                 _, n := test.a.Dims()
  41                 for _, chol := range []*Cholesky{
  42                         {},
  43                         {chol: NewTriDense(n-1, true, nil)},
  44                         {chol: NewTriDense(n, true, nil)},
  45                         {chol: NewTriDense(n+1, true, nil)},
  46                 } {
  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)
  50                         }
  51                         fc := DenseCopyOf(chol.chol)
  52                         if !Equal(fc, test.want) {
  53                                 t.Error("incorrect Cholesky factorization")
  54                         }
  55                         if math.Abs(test.cond-chol.cond) > 1e-13 {
  56                                 t.Errorf("Condition number mismatch: Want %v, got %v", test.cond, chol.cond)
  57                         }
  58                         U := chol.UTo(nil)
  59                         aCopy := DenseCopyOf(test.a)
  60                         var a Dense
  61                         a.Mul(U.TTri(), U)
  62                         if !EqualApprox(&a, aCopy, 1e-14) {
  63                                 t.Error("unexpected Cholesky factor product")
  64                         }
  65
  66                         L := chol.LTo(nil)
  67                         a.Mul(L, L.TTri())
  68                         if !EqualApprox(&a, aCopy, 1e-14) {
  69                                 t.Error("unexpected Cholesky factor product")
  70                         }
  71                 }
  72         }
  73 }
  74
  75 func TestCholeskySolve(t *testing.T) {
  76         for _, test := range []struct {
  77                 a   *SymDense
  78                 b   *Dense
  79                 ans *Dense
  80         }{
  81                 {
  82                         a: NewSymDense(2, []float64{
  83                                 1, 0,
  84                                 0, 1,
  85                         }),
  86                         b:   NewDense(2, 1, []float64{5, 6}),
  87                         ans: NewDense(2, 1, []float64{5, 6}),
  88                 },
  89                 {
  90                         a: NewSymDense(3, []float64{
  91                                 53, 59, 37,
  92                                 0, 83, 71,
  93                                 37, 71, 101,
  94                         }),
  95                         b:   NewDense(3, 1, []float64{5, 6, 7}),
  96                         ans: NewDense(3, 1, []float64{0.20745069393718094, -0.17421475529583694, 0.11577794010226464}),
  97                 },
  98         } {
  99                 var chol Cholesky
 100                 ok := chol.Factorize(test.a)
 101                 if !ok {
 102                         t.Fatal("unexpected Cholesky factorization failure: not positive definite")
 103                 }
 104
 105                 var x Dense
 106                 chol.Solve(&x, test.b)
 107                 if !EqualApprox(&x, test.ans, 1e-12) {
 108                         t.Error("incorrect Cholesky solve solution")
 109                 }
 110
 111                 var ans Dense
 112                 ans.Mul(test.a, &x)
 113                 if !EqualApprox(&ans, test.b, 1e-12) {
 114                         t.Error("incorrect Cholesky solve solution product")
 115                 }
 116         }
 117 }
 118
 119 func TestCholeskySolveChol(t *testing.T) {
 120         for _, test := range []struct {
 121                 a, b *SymDense
 122         }{
 123                 {
 124                         a: NewSymDense(2, []float64{
 125                                 1, 0,
 126                                 0, 1,
 127                         }),
 128                         b: NewSymDense(2, []float64{
 129                                 1, 0,
 130                                 0, 1,
 131                         }),
 132                 },
 133                 {
 134                         a: NewSymDense(2, []float64{
 135                                 1, 0,
 136                                 0, 1,
 137                         }),
 138                         b: NewSymDense(2, []float64{
 139                                 2, 0,
 140                                 0, 2,
 141                         }),
 142                 },
 143                 {
 144                         a: NewSymDense(3, []float64{
 145                                 53, 59, 37,
 146                                 59, 83, 71,
 147                                 37, 71, 101,
 148                         }),
 149                         b: NewSymDense(3, []float64{
 150                                 2, -1, 0,
 151                                 -1, 2, -1,
 152                                 0, -1, 2,
 153                         }),
 154                 },
 155         } {
 156                 var chola, cholb Cholesky
 157                 ok := chola.Factorize(test.a)
 158                 if !ok {
 159                         t.Fatal("unexpected Cholesky factorization failure for a: not positive definite")
 160                 }
 161                 ok = cholb.Factorize(test.b)
 162                 if !ok {
 163                         t.Fatal("unexpected Cholesky factorization failure for b: not positive definite")
 164                 }
 165
 166                 var x Dense
 167                 chola.SolveChol(&x, &cholb)
 168
 169                 var ans Dense
 170                 ans.Mul(test.a, &x)
 171                 if !EqualApprox(&ans, test.b, 1e-12) {
 172                         var y Dense
 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))
 176                 }
 177         }
 178 }
 179
 180 func TestCholeskySolveVec(t *testing.T) {
 181         for _, test := range []struct {
 182                 a   *SymDense
 183                 b   *VecDense
 184                 ans *VecDense
 185         }{
 186                 {
 187                         a: NewSymDense(2, []float64{
 188                                 1, 0,
 189                                 0, 1,
 190                         }),
 191                         b:   NewVecDense(2, []float64{5, 6}),
 192                         ans: NewVecDense(2, []float64{5, 6}),
 193                 },
 194                 {
 195                         a: NewSymDense(3, []float64{
 196                                 53, 59, 37,
 197                                 0, 83, 71,
 198                                 0, 0, 101,
 199                         }),
 200                         b:   NewVecDense(3, []float64{5, 6, 7}),
 201                         ans: NewVecDense(3, []float64{0.20745069393718094, -0.17421475529583694, 0.11577794010226464}),
 202                 },
 203         } {
 204                 var chol Cholesky
 205                 ok := chol.Factorize(test.a)
 206                 if !ok {
 207                         t.Fatal("unexpected Cholesky factorization failure: not positive definite")
 208                 }
 209
 210                 var x VecDense
 211                 chol.SolveVec(&x, test.b)
 212                 if !EqualApprox(&x, test.ans, 1e-12) {
 213                         t.Error("incorrect Cholesky solve solution")
 214                 }
 215
 216                 var ans VecDense
 217                 ans.MulVec(test.a, &x)
 218                 if !EqualApprox(&ans, test.b, 1e-12) {
 219                         t.Error("incorrect Cholesky solve solution product")
 220                 }
 221         }
 222 }
 223
 224 func TestCholeskyToSym(t *testing.T) {
 225         for _, test := range []*SymDense{
 226                 NewSymDense(3, []float64{
 227                         53, 59, 37,
 228                         0, 83, 71,
 229                         0, 0, 101,
 230                 }),
 231         } {
 232                 var chol Cholesky
 233                 ok := chol.Factorize(test)
 234                 if !ok {
 235                         t.Fatal("unexpected Cholesky factorization failure: not positive definite")
 236                 }
 237                 s := chol.ToSym(nil)
 238
 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))
 241                 }
 242         }
 243 }
 244
 245 func TestCloneCholesky(t *testing.T) {
 246         for _, test := range []*SymDense{
 247                 NewSymDense(3, []float64{
 248                         53, 59, 37,
 249                         0, 83, 71,
 250                         0, 0, 101,
 251                 }),
 252         } {
 253                 var chol Cholesky
 254                 ok := chol.Factorize(test)
 255                 if !ok {
 256                         panic("bad test")
 257                 }
 258                 var chol2 Cholesky
 259                 chol2.Clone(&chol)
 260
 261                 if chol.cond != chol2.cond {
 262                         t.Errorf("condition number mismatch from zero")
 263                 }
 264                 if !Equal(chol.chol, chol2.chol) {
 265                         t.Errorf("chol mismatch from zero")
 266                 }
 267
 268                 // Corrupt chol2 and try again
 269                 chol2.cond = math.NaN()
 270                 chol2.chol = NewTriDense(2, Upper, nil)
 271                 chol2.Clone(&chol)
 272                 if chol.cond != chol2.cond {
 273                         t.Errorf("condition number mismatch from non-zero")
 274                 }
 275                 if !Equal(chol.chol, chol2.chol) {
 276                         t.Errorf("chol mismatch from non-zero")
 277                 }
 278         }
 279 }
 280
 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()
 286                 }
 287                 var s SymDense
 288                 s.SymOuterK(1, NewDense(n, n, data))
 289
 290                 var chol Cholesky
 291                 ok := chol.Factorize(&s)
 292                 if !ok {
 293                         t.Errorf("Bad test, cholesky decomposition failed")
 294                 }
 295
 296                 var sInv SymDense
 297                 chol.InverseTo(&sInv)
 298
 299                 var ans Dense
 300                 ans.Mul(&sInv, &s)
 301                 if !equalApprox(eye(n), &ans, 1e-8, false) {
 302                         var diff Dense
 303                         diff.Sub(eye(n), &ans)
 304                         t.Errorf("SymDense times Cholesky inverse not identity. Norm diff = %v", Norm(&diff, 2))
 305                 }
 306         }
 307 }
 308
 309 func TestCholeskySymRankOne(t *testing.T) {
 310         rand.Seed(1)
 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()
 316                         }
 317
 318                         var a SymDense
 319                         a.SymOuterK(1, NewDense(n, n, data))
 320
 321                         xdata := make([]float64, n)
 322                         for i := range xdata {
 323                                 xdata[i] = rand.NormFloat64()
 324                         }
 325                         x := NewVecDense(n, xdata)
 326
 327                         var chol Cholesky
 328                         ok := chol.Factorize(&a)
 329                         if !ok {
 330                                 t.Errorf("Bad random test, Cholesky factorization failed")
 331                                 continue
 332                         }
 333
 334                         alpha := rand.Float64()
 335                         ok = chol.SymRankOne(&chol, alpha, x)
 336                         if !ok {
 337                                 t.Errorf("n=%v, alpha=%v: unexpected failure", n, alpha)
 338                                 continue
 339                         }
 340                         a.SymRankOne(&a, alpha, x)
 341
 342                         var achol SymDense
 343                         chol.ToSym(&achol)
 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))
 347                         }
 348                 }
 349         }
 350
 351         for i, test := range []struct {
 352                 a     *SymDense
 353                 alpha float64
 354                 x     []float64
 355
 356                 wantOk bool
 357         }{
 358                 {
 359                         // Update (to positive definite matrix).
 360                         a: NewSymDense(4, []float64{
 361                                 1, 1, 1, 1,
 362                                 0, 2, 3, 4,
 363                                 0, 0, 6, 10,
 364                                 0, 0, 0, 20,
 365                         }),
 366                         alpha:  1,
 367                         x:      []float64{0, 0, 0, 1},
 368                         wantOk: true,
 369                 },
 370                 {
 371                         // Downdate to singular matrix.
 372                         a: NewSymDense(4, []float64{
 373                                 1, 1, 1, 1,
 374                                 0, 2, 3, 4,
 375                                 0, 0, 6, 10,
 376                                 0, 0, 0, 20,
 377                         }),
 378                         alpha:  -1,
 379                         x:      []float64{0, 0, 0, 1},
 380                         wantOk: false,
 381                 },
 382                 {
 383                         // Downdate to positive definite matrix.
 384                         a: NewSymDense(4, []float64{
 385                                 1, 1, 1, 1,
 386                                 0, 2, 3, 4,
 387                                 0, 0, 6, 10,
 388                                 0, 0, 0, 20,
 389                         }),
 390                         alpha:  -1 / 2,
 391                         x:      []float64{0, 0, 0, 1},
 392                         wantOk: true,
 393                 },
 394         } {
 395                 var chol Cholesky
 396                 ok := chol.Factorize(test.a)
 397                 if !ok {
 398                         t.Errorf("Case %v: bad test, Cholesky factorization failed", i)
 399                         continue
 400                 }
 401
 402                 x := NewVecDense(len(test.x), test.x)
 403                 ok = chol.SymRankOne(&chol, test.alpha, x)
 404                 if !ok {
 405                         if test.wantOk {
 406                                 t.Errorf("Case %v: unexpected failure from SymRankOne", i)
 407                         }
 408                         continue
 409                 }
 410                 if ok && !test.wantOk {
 411                         t.Errorf("Case %v: expected a failure from SymRankOne", i)
 412                 }
 413
 414                 a := test.a
 415                 a.SymRankOne(a, test.alpha, x)
 416
 417                 var achol SymDense
 418                 chol.ToSym(&achol)
 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))
 422                 }
 423         }
 424 }
 425
 426 func TestCholeskyExtendVecSym(t *testing.T) {
 427         for cas, test := range []struct {
 428                 a *SymDense
 429         }{
 430                 {
 431                         a: NewSymDense(3, []float64{
 432                                 4, 1, 1,
 433                                 0, 2, 3,
 434                                 0, 0, 6,
 435                         }),
 436                 },
 437         } {
 438                 n := test.a.Symmetric()
 439                 as := test.a.SliceSquare(0, n-1).(*SymDense)
 440
 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)
 445                 if !ok {
 446                         panic("mat: bad test, matrix not positive definite")
 447                 }
 448
 449                 var chol Cholesky
 450                 ok = chol.Factorize(as)
 451                 if !ok {
 452                         panic("mat: bad test, subset is not positive definite")
 453                 }
 454                 row := NewVecDense(n, nil)
 455                 for i := 0; i < n; i++ {
 456                         row.SetVec(i, test.a.At(n-1, i))
 457                 }
 458
 459                 var cholNew Cholesky
 460                 ok = cholNew.ExtendVecSym(&chol, row)
 461                 if !ok {
 462                         t.Errorf("cas %v: update not positive definite", cas)
 463                 }
 464                 a := cholNew.ToSym(nil)
 465                 if !EqualApprox(a, test.a, 1e-12) {
 466                         t.Errorf("cas %v: mismatch", cas)
 467                 }
 468
 469                 // test in-place
 470                 ok = chol.ExtendVecSym(&chol, row)
 471                 if !ok {
 472                         t.Errorf("cas %v: in-place update not positive definite", cas)
 473                 }
 474                 if !equalChol(&chol, &cholNew) {
 475                         t.Errorf("cas %v: Cholesky different in-place vs. new", cas)
 476                 }
 477
 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)
 483                 }
 484         }
 485 }
 486
 487 func TestCholeskyScale(t *testing.T) {
 488         for cas, test := range []struct {
 489                 a *SymDense
 490                 f float64
 491         }{
 492                 {
 493                         a: NewSymDense(3, []float64{
 494                                 4, 1, 1,
 495                                 0, 2, 3,
 496                                 0, 0, 6,
 497                         }),
 498                         f: 0.5,
 499                 },
 500         } {
 501                 var chol Cholesky
 502                 ok := chol.Factorize(test.a)
 503                 if !ok {
 504                         t.Errorf("Case %v: bad test, Cholesky factorization failed", cas)
 505                         continue
 506                 }
 507
 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)
 514                 }
 515                 var sym SymDense
 516                 chol.ToSym(&sym)
 517                 var comp SymDense
 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)
 521                 }
 522
 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)
 527                 }
 528         }
 529 }
 530
 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
 534 }
 535
 536 // equalApproxChol checks that the two Cholesky decompositions are approximately
 537 // the same with the given tolerance on equality for the Triangular component and
 538 // condition.
 539 func equalApproxChol(a, b *Cholesky, matTol, condTol float64) bool {
 540         if !EqualApprox(a.chol, b.chol, matTol) {
 541                 return false
 542         }
 543         return floats.EqualWithinAbsOrRel(a.cond, b.cond, condTol, condTol)
 544 }
 545
 546 func BenchmarkCholeskySmall(b *testing.B) {
 547         benchmarkCholesky(b, 2)
 548 }
 549
 550 func BenchmarkCholeskyMedium(b *testing.B) {
 551         benchmarkCholesky(b, testblas.MediumMat)
 552 }
 553
 554 func BenchmarkCholeskyLarge(b *testing.B) {
 555         benchmarkCholesky(b, testblas.LargeMat)
 556 }
 557
 558 func benchmarkCholesky(b *testing.B, n int) {
 559         base := make([]float64, n*n)
 560         for i := range base {
 561                 base[i] = rand.Float64()
 562         }
 563         bm := NewDense(n, n, base)
 564         bm.Mul(bm.T(), bm)
 565         am := NewSymDense(n, bm.mat.Data)
 566
 567         var chol Cholesky
 568         b.ResetTimer()
 569         for i := 0; i < b.N; i++ {
 570                 ok := chol.Factorize(am)
 571                 if !ok {
 572                         panic("not pos def")
 573                 }
 574         }
 575 }