vendor/gonum.org/v1/gonum/mat/qr_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/blas64"
  14 )
  15
  16 func TestQR(t *testing.T) {
  17         for _, test := range []struct {
  18                 m, n int
  19         }{
  20                 {5, 5},
  21                 {10, 5},
  22         } {
  23                 m := test.m
  24                 n := test.n
  25                 a := NewDense(m, n, nil)
  26                 for i := 0; i < m; i++ {
  27                         for j := 0; j < n; j++ {
  28                                 a.Set(i, j, rand.NormFloat64())
  29                         }
  30                 }
  31                 var want Dense
  32                 want.Clone(a)
  33
  34                 var qr QR
  35                 qr.Factorize(a)
  36                 q := qr.QTo(nil)
  37
  38                 if !isOrthonormal(q, 1e-10) {
  39                         t.Errorf("Q is not orthonormal: m = %v, n = %v", m, n)
  40                 }
  41
  42                 r := qr.RTo(nil)
  43
  44                 var got Dense
  45                 got.Mul(q, r)
  46                 if !EqualApprox(&got, &want, 1e-12) {
  47                         t.Errorf("QR does not equal original matrix. \nWant: %v\nGot: %v", want, got)
  48                 }
  49         }
  50 }
  51
  52 func isOrthonormal(q *Dense, tol float64) bool {
  53         m, n := q.Dims()
  54         if m != n {
  55                 return false
  56         }
  57         for i := 0; i < m; i++ {
  58                 for j := i; j < m; j++ {
  59                         dot := blas64.Dot(m,
  60                                 blas64.Vector{Inc: 1, Data: q.mat.Data[i*q.mat.Stride:]},
  61                                 blas64.Vector{Inc: 1, Data: q.mat.Data[j*q.mat.Stride:]},
  62                         )
  63                         // Dot product should be 1 if i == j and 0 otherwise.
  64                         if i == j && math.Abs(dot-1) > tol {
  65                                 return false
  66                         }
  67                         if i != j && math.Abs(dot) > tol {
  68                                 return false
  69                         }
  70                 }
  71         }
  72         return true
  73 }
  74
  75 func TestSolveQR(t *testing.T) {
  76         for _, trans := range []bool{false, true} {
  77                 for _, test := range []struct {
  78                         m, n, bc int
  79                 }{
  80                         {5, 5, 1},
  81                         {10, 5, 1},
  82                         {5, 5, 3},
  83                         {10, 5, 3},
  84                 } {
  85                         m := test.m
  86                         n := test.n
  87                         bc := test.bc
  88                         a := NewDense(m, n, nil)
  89                         for i := 0; i < m; i++ {
  90                                 for j := 0; j < n; j++ {
  91                                         a.Set(i, j, rand.Float64())
  92                                 }
  93                         }
  94                         br := m
  95                         if trans {
  96                                 br = n
  97                         }
  98                         b := NewDense(br, bc, nil)
  99                         for i := 0; i < br; i++ {
 100                                 for j := 0; j < bc; j++ {
 101                                         b.Set(i, j, rand.Float64())
 102                                 }
 103                         }
 104                         var x Dense
 105                         var qr QR
 106                         qr.Factorize(a)
 107                         qr.Solve(&x, trans, b)
 108
 109                         // Test that the normal equations hold.
 110                         // A^T * A * x = A^T * b if !trans
 111                         // A * A^T * x = A * b if trans
 112                         var lhs Dense
 113                         var rhs Dense
 114                         if trans {
 115                                 var tmp Dense
 116                                 tmp.Mul(a, a.T())
 117                                 lhs.Mul(&tmp, &x)
 118                                 rhs.Mul(a, b)
 119                         } else {
 120                                 var tmp Dense
 121                                 tmp.Mul(a.T(), a)
 122                                 lhs.Mul(&tmp, &x)
 123                                 rhs.Mul(a.T(), b)
 124                         }
 125                         if !EqualApprox(&lhs, &rhs, 1e-10) {
 126                                 t.Errorf("Normal equations do not hold.\nLHS: %v\n, RHS: %v\n", lhs, rhs)
 127                         }
 128                 }
 129         }
 130         // TODO(btracey): Add in testOneInput when it exists.
 131 }
 132
 133 func TestSolveQRVec(t *testing.T) {
 134         for _, trans := range []bool{false, true} {
 135                 for _, test := range []struct {
 136                         m, n int
 137                 }{
 138                         {5, 5},
 139                         {10, 5},
 140                 } {
 141                         m := test.m
 142                         n := test.n
 143                         a := NewDense(m, n, nil)
 144                         for i := 0; i < m; i++ {
 145                                 for j := 0; j < n; j++ {
 146                                         a.Set(i, j, rand.Float64())
 147                                 }
 148                         }
 149                         br := m
 150                         if trans {
 151                                 br = n
 152                         }
 153                         b := NewVecDense(br, nil)
 154                         for i := 0; i < br; i++ {
 155                                 b.SetVec(i, rand.Float64())
 156                         }
 157                         var x VecDense
 158                         var qr QR
 159                         qr.Factorize(a)
 160                         qr.SolveVec(&x, trans, b)
 161
 162                         // Test that the normal equations hold.
 163                         // A^T * A * x = A^T * b if !trans
 164                         // A * A^T * x = A * b if trans
 165                         var lhs Dense
 166                         var rhs Dense
 167                         if trans {
 168                                 var tmp Dense
 169                                 tmp.Mul(a, a.T())
 170                                 lhs.Mul(&tmp, &x)
 171                                 rhs.Mul(a, b)
 172                         } else {
 173                                 var tmp Dense
 174                                 tmp.Mul(a.T(), a)
 175                                 lhs.Mul(&tmp, &x)
 176                                 rhs.Mul(a.T(), b)
 177                         }
 178                         if !EqualApprox(&lhs, &rhs, 1e-10) {
 179                                 t.Errorf("Normal equations do not hold.\nLHS: %v\n, RHS: %v\n", lhs, rhs)
 180                         }
 181                 }
 182         }
 183         // TODO(btracey): Add in testOneInput when it exists.
 184 }
 185
 186 func TestSolveQRCond(t *testing.T) {
 187         for _, test := range []*Dense{
 188                 NewDense(2, 2, []float64{1, 0, 0, 1e-20}),
 189                 NewDense(3, 2, []float64{1, 0, 0, 1e-20, 0, 0}),
 190         } {
 191                 m, _ := test.Dims()
 192                 var qr QR
 193                 qr.Factorize(test)
 194                 b := NewDense(m, 2, nil)
 195                 var x Dense
 196                 if err := qr.Solve(&x, false, b); err == nil {
 197                         t.Error("No error for near-singular matrix in matrix solve.")
 198                 }
 199
 200                 bvec := NewVecDense(m, nil)
 201                 var xvec VecDense
 202                 if err := qr.SolveVec(&xvec, false, bvec); err == nil {
 203                         t.Error("No error for near-singular matrix in matrix solve.")
 204                 }
 205         }
 206 }