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

   1 // Copyright ©2015 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_test
   6
   7 import (
   8         "fmt"
   9
  10         "gonum.org/v1/gonum/mat"
  11 )
  12
  13 func ExampleCholesky() {
  14         // Construct a symmetric positive definite matrix.
  15         tmp := mat.NewDense(4, 4, []float64{
  16                 2, 6, 8, -4,
  17                 1, 8, 7, -2,
  18                 2, 2, 1, 7,
  19                 8, -2, -2, 1,
  20         })
  21         var a mat.SymDense
  22         a.SymOuterK(1, tmp)
  23
  24         fmt.Printf("a = %0.4v\n", mat.Formatted(&a, mat.Prefix("    ")))
  25
  26         // Compute the cholesky factorization.
  27         var chol mat.Cholesky
  28         if ok := chol.Factorize(&a); !ok {
  29                 fmt.Println("a matrix is not positive semi-definite.")
  30         }
  31
  32         // Find the determinant.
  33         fmt.Printf("\nThe determinant of a is %0.4g\n\n", chol.Det())
  34
  35         // Use the factorization to solve the system of equations a * x = b.
  36         b := mat.NewVecDense(4, []float64{1, 2, 3, 4})
  37         var x mat.VecDense
  38         if err := chol.SolveVec(&x, b); err != nil {
  39                 fmt.Println("Matrix is near singular: ", err)
  40         }
  41         fmt.Println("Solve a * x = b")
  42         fmt.Printf("x = %0.4v\n", mat.Formatted(&x, mat.Prefix("    ")))
  43
  44         // Extract the factorization and check that it equals the original matrix.
  45         t := chol.LTo(nil)
  46         var test mat.Dense
  47         test.Mul(t, t.T())
  48         fmt.Println()
  49         fmt.Printf("L * L^T = %0.4v\n", mat.Formatted(&a, mat.Prefix("          ")))
  50
  51         // Output:
  52         // a = ⎡120  114   -4  -16⎤
  53         //     ⎢114  118   11  -24⎥
  54         //     ⎢ -4   11   58   17⎥
  55         //     ⎣-16  -24   17   73⎦
  56         //
  57         // The determinant of a is 1.543e+06
  58         //
  59         // Solve a * x = b
  60         // x = ⎡  -0.239⎤
  61         //     ⎢  0.2732⎥
  62         //     ⎢-0.04681⎥
  63         //     ⎣  0.1031⎦
  64         //
  65         // L * L^T = ⎡120  114   -4  -16⎤
  66         //           ⎢114  118   11  -24⎥
  67         //           ⎢ -4   11   58   17⎥
  68         //           ⎣-16  -24   17   73⎦
  69 }
  70
  71 func ExampleCholesky_SymRankOne() {
  72         a := mat.NewSymDense(4, []float64{
  73                 1, 1, 1, 1,
  74                 0, 2, 3, 4,
  75                 0, 0, 6, 10,
  76                 0, 0, 0, 20,
  77         })
  78         fmt.Printf("A = %0.4v\n", mat.Formatted(a, mat.Prefix("    ")))
  79
  80         // Compute the Cholesky factorization.
  81         var chol mat.Cholesky
  82         if ok := chol.Factorize(a); !ok {
  83                 fmt.Println("matrix a is not positive definite.")
  84         }
  85
  86         x := mat.NewVecDense(4, []float64{0, 0, 0, 1})
  87         fmt.Printf("\nx = %0.4v\n", mat.Formatted(x, mat.Prefix("    ")))
  88
  89         // Rank-1 update the factorization.
  90         chol.SymRankOne(&chol, 1, x)
  91         // Rank-1 update the matrix a.
  92         a.SymRankOne(a, 1, x)
  93
  94         au := chol.ToSym(nil)
  95
  96         // Print the matrix that was updated directly.
  97         fmt.Printf("\nA' =        %0.4v\n", mat.Formatted(a, mat.Prefix("            ")))
  98         // Print the matrix recovered from the factorization.
  99         fmt.Printf("\nU'^T * U' = %0.4v\n", mat.Formatted(au, mat.Prefix("            ")))
 100
 101         // Output:
 102         // A = ⎡ 1   1   1   1⎤
 103         //     ⎢ 1   2   3   4⎥
 104         //     ⎢ 1   3   6  10⎥
 105         //     ⎣ 1   4  10  20⎦
 106         //
 107         // x = ⎡0⎤
 108         //     ⎢0⎥
 109         //     ⎢0⎥
 110         //     ⎣1⎦
 111         //
 112         // A' =        ⎡ 1   1   1   1⎤
 113         //             ⎢ 1   2   3   4⎥
 114         //             ⎢ 1   3   6  10⎥
 115         //             ⎣ 1   4  10  21⎦
 116         //
 117         // U'^T * U' = ⎡ 1   1   1   1⎤
 118         //             ⎢ 1   2   3   4⎥
 119         //             ⎢ 1   3   6  10⎥
 120         //             ⎣ 1   4  10  21⎦
 121 }