vendor/gonum.org/v1/gonum/mat/eigen_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         "testing"
   9
  10         "golang.org/x/exp/rand"
  11
  12         "gonum.org/v1/gonum/floats"
  13 )
  14
  15 func TestEigen(t *testing.T) {
  16         for i, test := range []struct {
  17                 a *Dense
  18
  19                 values []complex128
  20                 left   *Dense
  21                 right  *Dense
  22         }{
  23                 {
  24                         a: NewDense(3, 3, []float64{
  25                                 1, 0, 0,
  26                                 0, 1, 0,
  27                                 0, 0, 1,
  28                         }),
  29                         values: []complex128{1, 1, 1},
  30                         left: NewDense(3, 3, []float64{
  31                                 1, 0, 0,
  32                                 0, 1, 0,
  33                                 0, 0, 1,
  34                         }),
  35                         right: NewDense(3, 3, []float64{
  36                                 1, 0, 0,
  37                                 0, 1, 0,
  38                                 0, 0, 1,
  39                         }),
  40                 },
  41         } {
  42                 var e1, e2, e3, e4 Eigen
  43                 ok := e1.Factorize(test.a, true, true)
  44                 if !ok {
  45                         panic("bad factorization")
  46                 }
  47                 e2.Factorize(test.a, false, true)
  48                 e3.Factorize(test.a, true, false)
  49                 e4.Factorize(test.a, false, false)
  50
  51                 v1 := e1.Values(nil)
  52                 if !cmplxEqual(v1, test.values) {
  53                         t.Errorf("eigenvector mismatch. Case %v", i)
  54                 }
  55                 if !Equal(e1.LeftVectors(), test.left) {
  56                         t.Errorf("left eigenvector mismatch. Case %v", i)
  57                 }
  58                 if !Equal(e1.Vectors(), test.right) {
  59                         t.Errorf("right eigenvector mismatch. Case %v", i)
  60                 }
  61
  62                 // Check that the eigenvectors and values are the same in all combinations.
  63                 if !cmplxEqual(v1, e2.Values(nil)) {
  64                         t.Errorf("eigenvector mismatch. Case %v", i)
  65                 }
  66                 if !cmplxEqual(v1, e3.Values(nil)) {
  67                         t.Errorf("eigenvector mismatch. Case %v", i)
  68                 }
  69                 if !cmplxEqual(v1, e4.Values(nil)) {
  70                         t.Errorf("eigenvector mismatch. Case %v", i)
  71                 }
  72                 if !Equal(e1.Vectors(), e2.Vectors()) {
  73                         t.Errorf("right eigenvector mismatch. Case %v", i)
  74                 }
  75                 if !Equal(e1.LeftVectors(), e3.LeftVectors()) {
  76                         t.Errorf("right eigenvector mismatch. Case %v", i)
  77                 }
  78
  79                 // TODO(btracey): Also add in a test for correctness when #308 is
  80                 // resolved and we have a CMat.Mul().
  81         }
  82 }
  83
  84 func cmplxEqual(v1, v2 []complex128) bool {
  85         for i, v := range v1 {
  86                 if v != v2[i] {
  87                         return false
  88                 }
  89         }
  90         return true
  91 }
  92
  93 func TestSymEigen(t *testing.T) {
  94         // Hand coded tests with results from lapack.
  95         for _, test := range []struct {
  96                 mat *SymDense
  97
  98                 values  []float64
  99                 vectors *Dense
 100         }{
 101                 {
 102                         mat:    NewSymDense(3, []float64{8, 2, 4, 2, 6, 10, 4, 10, 5}),
 103                         values: []float64{-4.707679201365891, 6.294580208480216, 17.413098992885672},
 104                         vectors: NewDense(3, 3, []float64{
 105                                 -0.127343483135656, -0.902414161226903, -0.411621572466779,
 106                                 -0.664177720955769, 0.385801900032553, -0.640331827193739,
 107                                 0.736648893495999, 0.191847792659746, -0.648492738712395,
 108                         }),
 109                 },
 110         } {
 111                 var es EigenSym
 112                 ok := es.Factorize(test.mat, true)
 113                 if !ok {
 114                         t.Errorf("bad factorization")
 115                 }
 116                 if !floats.EqualApprox(test.values, es.values, 1e-14) {
 117                         t.Errorf("Eigenvalue mismatch")
 118                 }
 119                 if !EqualApprox(test.vectors, es.vectors, 1e-14) {
 120                         t.Errorf("Eigenvector mismatch")
 121                 }
 122
 123                 var es2 EigenSym
 124                 es2.Factorize(test.mat, false)
 125                 if !floats.EqualApprox(es2.values, es.values, 1e-14) {
 126                         t.Errorf("Eigenvalue mismatch when no vectors computed")
 127                 }
 128         }
 129
 130         // Randomized tests
 131         rnd := rand.New(rand.NewSource(1))
 132         for _, n := range []int{3, 5, 10, 70} {
 133                 for cas := 0; cas < 10; cas++ {
 134                         a := make([]float64, n*n)
 135                         for i := range a {
 136                                 a[i] = rnd.NormFloat64()
 137                         }
 138                         s := NewSymDense(n, a)
 139                         var es EigenSym
 140                         ok := es.Factorize(s, true)
 141                         if !ok {
 142                                 t.Errorf("Bad test")
 143                         }
 144
 145                         // Check that the eigenvectors are orthonormal.
 146                         if !isOrthonormal(es.vectors, 1e-8) {
 147                                 t.Errorf("Eigenvectors not orthonormal")
 148                         }
 149
 150                         // Check that the eigenvalues are actually eigenvalues.
 151                         for i := 0; i < n; i++ {
 152                                 v := NewVecDense(n, Col(nil, i, es.vectors))
 153                                 var m VecDense
 154                                 m.MulVec(s, v)
 155
 156                                 var scal VecDense
 157                                 scal.ScaleVec(es.values[i], v)
 158
 159                                 if !EqualApprox(&m, &scal, 1e-8) {
 160                                         t.Errorf("Eigenvalue does not match")
 161                                 }
 162                         }
 163                 }
 164         }
 165 }