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

   1 // Copyright ©2017 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         "errors"
   9
  10         "gonum.org/v1/gonum/blas/blas64"
  11 )
  12
  13 // HOGSVD is a type for creating and using the Higher Order Generalized Singular Value
  14 // Decomposition (HOGSVD) of a set of matrices.
  15 //
  16 // The factorization is a linear transformation of the data sets from the given
  17 // variable×sample spaces to reduced and diagonalized "eigenvariable"×"eigensample"
  18 // spaces.
  19 type HOGSVD struct {
  20         n int
  21         v *Dense
  22         b []Dense
  23
  24         err error
  25 }
  26
  27 // Factorize computes the higher order generalized singular value decomposition (HOGSVD)
  28 // of the n input r_i×c column tall matrices in m. HOGSV extends the GSVD case from 2 to n
  29 // input matrices.
  30 //
  31 //  M_0 = U_0 * Σ_0 * V^T
  32 //  M_1 = U_1 * Σ_1 * V^T
  33 //  .
  34 //  .
  35 //  .
  36 //  M_{n-1} = U_{n-1} * Σ_{n-1} * V^T
  37 //
  38 // where U_i are r_i×c matrices of singular vectors, Σ are c×c matrices singular values, and V
  39 // is a c×c matrix of singular vectors.
  40 //
  41 // Factorize returns whether the decomposition succeeded. If the decomposition
  42 // failed, routines that require a successful factorization will panic.
  43 func (gsvd *HOGSVD) Factorize(m ...Matrix) (ok bool) {
  44         // Factorize performs the HOGSVD factorisation
  45         // essentially as described by Ponnapalli et al.
  46         // https://doi.org/10.1371/journal.pone.0028072
  47
  48         if len(m) < 2 {
  49                 panic("hogsvd: too few matrices")
  50         }
  51         gsvd.n = 0
  52
  53         r, c := m[0].Dims()
  54         a := make([]Cholesky, len(m))
  55         var ts SymDense
  56         for i, d := range m {
  57                 rd, cd := d.Dims()
  58                 if rd < cd {
  59                         gsvd.err = ErrShape
  60                         return false
  61                 }
  62                 if rd > r {
  63                         r = rd
  64                 }
  65                 if cd != c {
  66                         panic(ErrShape)
  67                 }
  68                 ts.Reset()
  69                 ts.SymOuterK(1, d.T())
  70                 ok = a[i].Factorize(&ts)
  71                 if !ok {
  72                         gsvd.err = errors.New("hogsvd: cholesky decomposition failed")
  73                         return false
  74                 }
  75         }
  76
  77         s := getWorkspace(c, c, true)
  78         defer putWorkspace(s)
  79         sij := getWorkspace(c, c, false)
  80         defer putWorkspace(sij)
  81         for i, ai := range a {
  82                 for _, aj := range a[i+1:] {
  83                         gsvd.err = ai.SolveChol(sij, &aj)
  84                         if gsvd.err != nil {
  85                                 return false
  86                         }
  87                         s.Add(s, sij)
  88
  89                         gsvd.err = aj.SolveChol(sij, &ai)
  90                         if gsvd.err != nil {
  91                                 return false
  92                         }
  93                         s.Add(s, sij)
  94                 }
  95         }
  96         s.Scale(1/float64(len(m)*(len(m)-1)), s)
  97
  98         var eig Eigen
  99         ok = eig.Factorize(s.T(), false, true)
 100         if !ok {
 101                 gsvd.err = errors.New("hogsvd: eigen decomposition failed")
 102                 return false
 103         }
 104         v := eig.Vectors()
 105         var cv VecDense
 106         for j := 0; j < c; j++ {
 107                 cv.ColViewOf(v, j)
 108                 cv.ScaleVec(1/blas64.Nrm2(c, cv.mat), &cv)
 109         }
 110
 111         b := make([]Dense, len(m))
 112         biT := getWorkspace(c, r, false)
 113         defer putWorkspace(biT)
 114         for i, d := range m {
 115                 // All calls to reset will leave a zeroed
 116                 // matrix with capacity to store the result
 117                 // without additional allocation.
 118                 biT.Reset()
 119                 gsvd.err = biT.Solve(v, d.T())
 120                 if gsvd.err != nil {
 121                         return false
 122                 }
 123                 b[i].Clone(biT.T())
 124         }
 125
 126         gsvd.n = len(m)
 127         gsvd.v = v
 128         gsvd.b = b
 129         return true
 130 }
 131
 132 // Err returns the reason for a factorization failure.
 133 func (gsvd *HOGSVD) Err() error {
 134         return gsvd.err
 135 }
 136
 137 // Len returns the number of matrices that have been factorized. If Len returns
 138 // zero, the factorization was not successful.
 139 func (gsvd *HOGSVD) Len() int {
 140         return gsvd.n
 141 }
 142
 143 // UTo extracts the matrix U_n from the singular value decomposition, storing
 144 // the result in-place into dst. U_n is size r×c.
 145 // If dst is nil, a new matrix is allocated. The resulting U matrix is returned.
 146 //
 147 // UTo will panic if the receiver does not contain a successful factorization.
 148 func (gsvd *HOGSVD) UTo(dst *Dense, n int) *Dense {
 149         if gsvd.n == 0 {
 150                 panic("hogsvd: unsuccessful factorization")
 151         }
 152         if n < 0 || gsvd.n <= n {
 153                 panic("hogsvd: invalid index")
 154         }
 155
 156         if dst == nil {
 157                 r, c := gsvd.b[n].Dims()
 158                 dst = NewDense(r, c, nil)
 159         } else {
 160                 dst.reuseAs(gsvd.b[n].Dims())
 161         }
 162         dst.Copy(&gsvd.b[n])
 163         var v VecDense
 164         for j, f := range gsvd.Values(nil, n) {
 165                 v.ColViewOf(dst, j)
 166                 v.ScaleVec(1/f, &v)
 167         }
 168         return dst
 169 }
 170
 171 // Values returns the nth set of singular values of the factorized system.
 172 // If the input slice is non-nil, the values will be stored in-place into the slice.
 173 // In this case, the slice must have length c, and Values will panic with
 174 // matrix.ErrSliceLengthMismatch otherwise. If the input slice is nil,
 175 // a new slice of the appropriate length will be allocated and returned.
 176 //
 177 // Values will panic if the receiver does not contain a successful factorization.
 178 func (gsvd *HOGSVD) Values(s []float64, n int) []float64 {
 179         if gsvd.n == 0 {
 180                 panic("hogsvd: unsuccessful factorization")
 181         }
 182         if n < 0 || gsvd.n <= n {
 183                 panic("hogsvd: invalid index")
 184         }
 185
 186         r, c := gsvd.b[n].Dims()
 187         if s == nil {
 188                 s = make([]float64, c)
 189         } else if len(s) != c {
 190                 panic(ErrSliceLengthMismatch)
 191         }
 192         var v VecDense
 193         for j := 0; j < c; j++ {
 194                 v.ColViewOf(&gsvd.b[n], j)
 195                 s[j] = blas64.Nrm2(r, v.mat)
 196         }
 197         return s
 198 }
 199
 200 // VTo extracts the matrix V from the singular value decomposition, storing
 201 // the result in-place into dst. V is size c×c.
 202 // If dst is nil, a new matrix is allocated. The resulting V matrix is returned.
 203 //
 204 // VTo will panic if the receiver does not contain a successful factorization.
 205 func (gsvd *HOGSVD) VTo(dst *Dense) *Dense {
 206         if gsvd.n == 0 {
 207                 panic("hogsvd: unsuccessful factorization")
 208         }
 209         if dst == nil {
 210                 r, c := gsvd.v.Dims()
 211                 dst = NewDense(r, c, nil)
 212         } else {
 213                 dst.reuseAs(gsvd.v.Dims())
 214         }
 215         dst.Copy(gsvd.v)
 216         return dst
 217 }