github.com/gonum/matrix@v0.0.0-20181209220409-c518dec07be9/mat64/cholesky_example_test.go (about)

     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 mat64_test
     6  
     7  import (
     8  	"fmt"
     9  
    10  	"github.com/gonum/matrix/mat64"
    11  )
    12  
    13  func ExampleCholesky() {
    14  	// Construct a symmetric positive definite matrix.
    15  	tmp := mat64.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 mat64.SymDense
    22  	a.SymOuterK(1, tmp)
    23  
    24  	fmt.Printf("a = %0.4v\n", mat64.Formatted(&a, mat64.Prefix("    ")))
    25  
    26  	// Compute the cholesky factorization.
    27  	var chol mat64.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 := mat64.NewVector(4, []float64{1, 2, 3, 4})
    37  	var x mat64.Vector
    38  	if err := x.SolveCholeskyVec(&chol, 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", mat64.Formatted(&x, mat64.Prefix("    ")))
    43  
    44  	// Extract the factorization and check that it equals the original matrix.
    45  	var t mat64.TriDense
    46  	t.LFromCholesky(&chol)
    47  	var test mat64.Dense
    48  	test.Mul(&t, t.T())
    49  	fmt.Println()
    50  	fmt.Printf("L * L^T = %0.4v\n", mat64.Formatted(&a, mat64.Prefix("          ")))
    51  
    52  	// Output:
    53  	// a = ⎡120  114   -4  -16⎤
    54  	//     ⎢114  118   11  -24⎥
    55  	//     ⎢ -4   11   58   17⎥
    56  	//     ⎣-16  -24   17   73⎦
    57  	//
    58  	// The determinant of a is 1.543e+06
    59  	//
    60  	// Solve a * x = b
    61  	// x = ⎡  -0.239⎤
    62  	//     ⎢  0.2732⎥
    63  	//     ⎢-0.04681⎥
    64  	//     ⎣  0.1031⎦
    65  	//
    66  	// L * L^T = ⎡120  114   -4  -16⎤
    67  	//           ⎢114  118   11  -24⎥
    68  	//           ⎢ -4   11   58   17⎥
    69  	//           ⎣-16  -24   17   73⎦
    70  }
    71  
    72  func ExampleCholeskySymRankOne() {
    73  	a := mat64.NewSymDense(4, []float64{
    74  		1, 1, 1, 1,
    75  		0, 2, 3, 4,
    76  		0, 0, 6, 10,
    77  		0, 0, 0, 20,
    78  	})
    79  	fmt.Printf("A = %0.4v\n", mat64.Formatted(a, mat64.Prefix("    ")))
    80  
    81  	// Compute the Cholesky factorization.
    82  	var chol mat64.Cholesky
    83  	if ok := chol.Factorize(a); !ok {
    84  		fmt.Println("matrix a is not positive definite.")
    85  	}
    86  
    87  	x := mat64.NewVector(4, []float64{0, 0, 0, 1})
    88  	fmt.Printf("\nx = %0.4v\n", mat64.Formatted(x, mat64.Prefix("    ")))
    89  
    90  	// Rank-1 update the factorization.
    91  	chol.SymRankOne(&chol, 1, x)
    92  	// Rank-1 update the matrix a.
    93  	a.SymRankOne(a, 1, x)
    94  
    95  	var au mat64.SymDense
    96  	au.FromCholesky(&chol)
    97  
    98  	// Print the matrix that was updated directly.
    99  	fmt.Printf("\nA' =        %0.4v\n", mat64.Formatted(a, mat64.Prefix("            ")))
   100  	// Print the matrix recovered from the factorization.
   101  	fmt.Printf("\nU'^T * U' = %0.4v\n", mat64.Formatted(&au, mat64.Prefix("            ")))
   102  
   103  	// Output:
   104  	// A = ⎡ 1   1   1   1⎤
   105  	//     ⎢ 1   2   3   4⎥
   106  	//     ⎢ 1   3   6  10⎥
   107  	//     ⎣ 1   4  10  20⎦
   108  	//
   109  	// x = ⎡0⎤
   110  	//     ⎢0⎥
   111  	//     ⎢0⎥
   112  	//     ⎣1⎦
   113  	//
   114  	// A' =        ⎡ 1   1   1   1⎤
   115  	//             ⎢ 1   2   3   4⎥
   116  	//             ⎢ 1   3   6  10⎥
   117  	//             ⎣ 1   4  10  21⎦
   118  	//
   119  	// U'^T * U' = ⎡ 1   1   1   1⎤
   120  	//             ⎢ 1   2   3   4⎥
   121  	//             ⎢ 1   3   6  10⎥
   122  	//             ⎣ 1   4  10  21⎦
   123  }