vendor/gonum.org/v1/gonum/mat/svd.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         "gonum.org/v1/gonum/blas/blas64"
   9         "gonum.org/v1/gonum/lapack"
  10         "gonum.org/v1/gonum/lapack/lapack64"
  11 )
  12
  13 // SVD is a type for creating and using the Singular Value Decomposition (SVD)
  14 // of a matrix.
  15 type SVD struct {
  16         kind SVDKind
  17
  18         s  []float64
  19         u  blas64.General
  20         vt blas64.General
  21 }
  22
  23 // Factorize computes the singular value decomposition (SVD) of the input matrix
  24 // A. The singular values of A are computed in all cases, while the singular
  25 // vectors are optionally computed depending on the input kind.
  26 //
  27 // The full singular value decomposition (kind == SVDFull) deconstructs A as
  28 //  A = U * Σ * V^T
  29 // where Σ is an m×n diagonal matrix of singular vectors, U is an m×m unitary
  30 // matrix of left singular vectors, and V is an n×n matrix of right singular vectors.
  31 //
  32 // It is frequently not necessary to compute the full SVD. Computation time and
  33 // storage costs can be reduced using the appropriate kind. Only the singular
  34 // values can be computed (kind == SVDNone), or a "thin" representation of the
  35 // singular vectors (kind = SVDThin). The thin representation can save a significant
  36 // amount of memory if m >> n. See the documentation for the lapack.SVDKind values
  37 // for more information.
  38 //
  39 // Factorize returns whether the decomposition succeeded. If the decomposition
  40 // failed, routines that require a successful factorization will panic.
  41 func (svd *SVD) Factorize(a Matrix, kind SVDKind) (ok bool) {
  42         m, n := a.Dims()
  43         var jobU, jobVT lapack.SVDJob
  44         switch kind {
  45         default:
  46                 panic("svd: bad input kind")
  47         case SVDNone:
  48                 jobU = lapack.SVDNone
  49                 jobVT = lapack.SVDNone
  50         case SVDFull:
  51                 // TODO(btracey): This code should be modified to have the smaller
  52                 // matrix written in-place into aCopy when the lapack/native/dgesvd
  53                 // implementation is complete.
  54                 svd.u = blas64.General{
  55                         Rows:   m,
  56                         Cols:   m,
  57                         Stride: m,
  58                         Data:   use(svd.u.Data, m*m),
  59                 }
  60                 svd.vt = blas64.General{
  61                         Rows:   n,
  62                         Cols:   n,
  63                         Stride: n,
  64                         Data:   use(svd.vt.Data, n*n),
  65                 }
  66                 jobU = lapack.SVDAll
  67                 jobVT = lapack.SVDAll
  68         case SVDThin:
  69                 // TODO(btracey): This code should be modified to have the larger
  70                 // matrix written in-place into aCopy when the lapack/native/dgesvd
  71                 // implementation is complete.
  72                 svd.u = blas64.General{
  73                         Rows:   m,
  74                         Cols:   min(m, n),
  75                         Stride: min(m, n),
  76                         Data:   use(svd.u.Data, m*min(m, n)),
  77                 }
  78                 svd.vt = blas64.General{
  79                         Rows:   min(m, n),
  80                         Cols:   n,
  81                         Stride: n,
  82                         Data:   use(svd.vt.Data, min(m, n)*n),
  83                 }
  84                 jobU = lapack.SVDInPlace
  85                 jobVT = lapack.SVDInPlace
  86         }
  87
  88         // A is destroyed on call, so copy the matrix.
  89         aCopy := DenseCopyOf(a)
  90         svd.kind = kind
  91         svd.s = use(svd.s, min(m, n))
  92
  93         work := []float64{0}
  94         lapack64.Gesvd(jobU, jobVT, aCopy.mat, svd.u, svd.vt, svd.s, work, -1)
  95         work = getFloats(int(work[0]), false)
  96         ok = lapack64.Gesvd(jobU, jobVT, aCopy.mat, svd.u, svd.vt, svd.s, work, len(work))
  97         putFloats(work)
  98         if !ok {
  99                 svd.kind = 0
 100         }
 101         return ok
 102 }
 103
 104 // Kind returns the matrix.SVDKind of the decomposition. If no decomposition has been
 105 // computed, Kind returns 0.
 106 func (svd *SVD) Kind() SVDKind {
 107         return svd.kind
 108 }
 109
 110 // Cond returns the 2-norm condition number for the factorized matrix. Cond will
 111 // panic if the receiver does not contain a successful factorization.
 112 func (svd *SVD) Cond() float64 {
 113         if svd.kind == 0 {
 114                 panic("svd: no decomposition computed")
 115         }
 116         return svd.s[0] / svd.s[len(svd.s)-1]
 117 }
 118
 119 // Values returns the singular values of the factorized matrix in decreasing order.
 120 // If the input slice is non-nil, the values will be stored in-place into the slice.
 121 // In this case, the slice must have length min(m,n), and Values will panic with
 122 // matrix.ErrSliceLengthMismatch otherwise. If the input slice is nil,
 123 // a new slice of the appropriate length will be allocated and returned.
 124 //
 125 // Values will panic if the receiver does not contain a successful factorization.
 126 func (svd *SVD) Values(s []float64) []float64 {
 127         if svd.kind == 0 {
 128                 panic("svd: no decomposition computed")
 129         }
 130         if s == nil {
 131                 s = make([]float64, len(svd.s))
 132         }
 133         if len(s) != len(svd.s) {
 134                 panic(ErrSliceLengthMismatch)
 135         }
 136         copy(s, svd.s)
 137         return s
 138 }
 139
 140 // UTo extracts the matrix U from the singular value decomposition, storing
 141 // the result in-place into dst. U is size m×m if svd.Kind() == SVDFull,
 142 // of size m×min(m,n) if svd.Kind() == SVDThin, and UTo panics otherwise.
 143 func (svd *SVD) UTo(dst *Dense) *Dense {
 144         kind := svd.kind
 145         if kind != SVDFull && kind != SVDThin {
 146                 panic("mat: improper SVD kind")
 147         }
 148         r := svd.u.Rows
 149         c := svd.u.Cols
 150         if dst == nil {
 151                 dst = NewDense(r, c, nil)
 152         } else {
 153                 dst.reuseAs(r, c)
 154         }
 155
 156         tmp := &Dense{
 157                 mat:     svd.u,
 158                 capRows: r,
 159                 capCols: c,
 160         }
 161         dst.Copy(tmp)
 162
 163         return dst
 164 }
 165
 166 // VTo extracts the matrix V from the singular value decomposition, storing
 167 // the result in-place into dst. V is size n×n if svd.Kind() == SVDFull,
 168 // of size n×min(m,n) if svd.Kind() == SVDThin, and VTo panics otherwise.
 169 func (svd *SVD) VTo(dst *Dense) *Dense {
 170         kind := svd.kind
 171         if kind != SVDFull && kind != SVDThin {
 172                 panic("mat: improper SVD kind")
 173         }
 174         r := svd.vt.Rows
 175         c := svd.vt.Cols
 176         if dst == nil {
 177                 dst = NewDense(c, r, nil)
 178         } else {
 179                 dst.reuseAs(c, r)
 180         }
 181
 182         tmp := &Dense{
 183                 mat:     svd.vt,
 184                 capRows: r,
 185                 capCols: c,
 186         }
 187         dst.Copy(tmp.T())
 188
 189         return dst
 190 }