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 5package mat 6 7import ( 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. 15type SVD struct { 16 kind SVDKind 17 18 s []float64 19 u blas64.General 20 vt blas64.General 21} 22 23// SVDKind specifies the treatment of singular vectors during an SVD 24// factorization. 25type SVDKind int 26 27const ( 28 // SVDNone specifies that no singular vectors should be computed during 29 // the decomposition. 30 SVDNone SVDKind = 0 31 32 // SVDThinU specifies the thin decomposition for U should be computed. 33 SVDThinU SVDKind = 1 << (iota - 1) 34 // SVDFullU specifies the full decomposition for U should be computed. 35 SVDFullU 36 // SVDThinV specifies the thin decomposition for V should be computed. 37 SVDThinV 38 // SVDFullV specifies the full decomposition for V should be computed. 39 SVDFullV 40 41 // SVDThin is a convenience value for computing both thin vectors. 42 SVDThin SVDKind = SVDThinU | SVDThinV 43 // SVDThin is a convenience value for computing both full vectors. 44 SVDFull SVDKind = SVDFullU | SVDFullV 45) 46 47// succFact returns whether the receiver contains a successful factorization. 48func (svd *SVD) succFact() bool { 49 return len(svd.s) != 0 50} 51 52// Factorize computes the singular value decomposition (SVD) of the input matrix A. 53// The singular values of A are computed in all cases, while the singular 54// vectors are optionally computed depending on the input kind. 55// 56// The full singular value decomposition (kind == SVDFull) is a factorization 57// of an m×n matrix A of the form 58// A = U * Σ * V^T 59// where Σ is an m×n diagonal matrix, U is an m×m orthogonal matrix, and V is an 60// n×n orthogonal matrix. The diagonal elements of Σ are the singular values of A. 61// The first min(m,n) columns of U and V are, respectively, the left and right 62// singular vectors of A. 63// 64// Significant storage space can be saved by using the thin representation of 65// the SVD (kind == SVDThin) instead of the full SVD, especially if 66// m >> n or m << n. The thin SVD finds 67// A = U~ * Σ * V~^T 68// where U~ is of size m×min(m,n), Σ is a diagonal matrix of size min(m,n)×min(m,n) 69// and V~ is of size n×min(m,n). 70// 71// Factorize returns whether the decomposition succeeded. If the decomposition 72// failed, routines that require a successful factorization will panic. 73func (svd *SVD) Factorize(a Matrix, kind SVDKind) (ok bool) { 74 // kill previous factorization 75 svd.s = svd.s[:0] 76 svd.kind = kind 77 78 m, n := a.Dims() 79 var jobU, jobVT lapack.SVDJob 80 81 // TODO(btracey): This code should be modified to have the smaller 82 // matrix written in-place into aCopy when the lapack/native/dgesvd 83 // implementation is complete. 84 switch { 85 case kind&SVDFullU != 0: 86 jobU = lapack.SVDAll 87 svd.u = blas64.General{ 88 Rows: m, 89 Cols: m, 90 Stride: m, 91 Data: use(svd.u.Data, m*m), 92 } 93 case kind&SVDThinU != 0: 94 jobU = lapack.SVDStore 95 svd.u = blas64.General{ 96 Rows: m, 97 Cols: min(m, n), 98 Stride: min(m, n), 99 Data: use(svd.u.Data, m*min(m, n)), 100 } 101 default: 102 jobU = lapack.SVDNone 103 } 104 switch { 105 case kind&SVDFullV != 0: 106 svd.vt = blas64.General{ 107 Rows: n, 108 Cols: n, 109 Stride: n, 110 Data: use(svd.vt.Data, n*n), 111 } 112 jobVT = lapack.SVDAll 113 case kind&SVDThinV != 0: 114 svd.vt = blas64.General{ 115 Rows: min(m, n), 116 Cols: n, 117 Stride: n, 118 Data: use(svd.vt.Data, min(m, n)*n), 119 } 120 jobVT = lapack.SVDStore 121 default: 122 jobVT = lapack.SVDNone 123 } 124 125 // A is destroyed on call, so copy the matrix. 126 aCopy := DenseCopyOf(a) 127 svd.kind = kind 128 svd.s = use(svd.s, min(m, n)) 129 130 work := []float64{0} 131 lapack64.Gesvd(jobU, jobVT, aCopy.mat, svd.u, svd.vt, svd.s, work, -1) 132 work = getFloats(int(work[0]), false) 133 ok = lapack64.Gesvd(jobU, jobVT, aCopy.mat, svd.u, svd.vt, svd.s, work, len(work)) 134 putFloats(work) 135 if !ok { 136 svd.kind = 0 137 } 138 return ok 139} 140 141// Kind returns the SVDKind of the decomposition. If no decomposition has been 142// computed, Kind returns -1. 143func (svd *SVD) Kind() SVDKind { 144 if !svd.succFact() { 145 return -1 146 } 147 return svd.kind 148} 149 150// Cond returns the 2-norm condition number for the factorized matrix. Cond will 151// panic if the receiver does not contain a successful factorization. 152func (svd *SVD) Cond() float64 { 153 if !svd.succFact() { 154 panic(badFact) 155 } 156 return svd.s[0] / svd.s[len(svd.s)-1] 157} 158 159// Values returns the singular values of the factorized matrix in descending order. 160// 161// If the input slice is non-nil, the values will be stored in-place into 162// the slice. In this case, the slice must have length min(m,n), and Values will 163// panic with ErrSliceLengthMismatch otherwise. If the input slice is nil, a new 164// slice of the appropriate length will be allocated and returned. 165// 166// Values will panic if the receiver does not contain a successful factorization. 167func (svd *SVD) Values(s []float64) []float64 { 168 if !svd.succFact() { 169 panic(badFact) 170 } 171 if s == nil { 172 s = make([]float64, len(svd.s)) 173 } 174 if len(s) != len(svd.s) { 175 panic(ErrSliceLengthMismatch) 176 } 177 copy(s, svd.s) 178 return s 179} 180 181// UTo extracts the matrix U from the singular value decomposition. The first 182// min(m,n) columns are the left singular vectors and correspond to the singular 183// values as returned from SVD.Values. 184// 185// If dst is not nil, U is stored in-place into dst, and dst must have size 186// m×m if the full U was computed, size m×min(m,n) if the thin U was computed, 187// and UTo panics otherwise. If dst is nil, a new matrix of the appropriate size 188// is allocated and returned. 189func (svd *SVD) UTo(dst *Dense) *Dense { 190 if !svd.succFact() { 191 panic(badFact) 192 } 193 kind := svd.kind 194 if kind&SVDThinU == 0 && kind&SVDFullU == 0 { 195 panic("svd: u not computed during factorization") 196 } 197 r := svd.u.Rows 198 c := svd.u.Cols 199 if dst == nil { 200 dst = NewDense(r, c, nil) 201 } else { 202 dst.reuseAs(r, c) 203 } 204 205 tmp := &Dense{ 206 mat: svd.u, 207 capRows: r, 208 capCols: c, 209 } 210 dst.Copy(tmp) 211 212 return dst 213} 214 215// VTo extracts the matrix V from the singular value decomposition. The first 216// min(m,n) columns are the right singular vectors and correspond to the singular 217// values as returned from SVD.Values. 218// 219// If dst is not nil, V is stored in-place into dst, and dst must have size 220// n×n if the full V was computed, size n×min(m,n) if the thin V was computed, 221// and VTo panics otherwise. If dst is nil, a new matrix of the appropriate size 222// is allocated and returned. 223func (svd *SVD) VTo(dst *Dense) *Dense { 224 if !svd.succFact() { 225 panic(badFact) 226 } 227 kind := svd.kind 228 if kind&SVDThinU == 0 && kind&SVDFullV == 0 { 229 panic("svd: v not computed during factorization") 230 } 231 r := svd.vt.Rows 232 c := svd.vt.Cols 233 if dst == nil { 234 dst = NewDense(c, r, nil) 235 } else { 236 dst.reuseAs(c, r) 237 } 238 239 tmp := &Dense{ 240 mat: svd.vt, 241 capRows: r, 242 capCols: c, 243 } 244 dst.Copy(tmp.T()) 245 246 return dst 247} 248