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 5package mat 6 7import ( 8 "gonum.org/v1/gonum/blas/blas64" 9 "gonum.org/v1/gonum/floats" 10 "gonum.org/v1/gonum/lapack" 11 "gonum.org/v1/gonum/lapack/lapack64" 12) 13 14// GSVDKind specifies the treatment of singular vectors during a GSVD 15// factorization. 16type GSVDKind int 17 18const ( 19 // GSVDNone specifies that no singular vectors should be computed during 20 // the decomposition. 21 GSVDNone GSVDKind = 0 22 23 // GSVDU specifies that the U singular vectors should be computed during 24 // the decomposition. 25 GSVDU GSVDKind = 1 << iota 26 // GSVDV specifies that the V singular vectors should be computed during 27 // the decomposition. 28 GSVDV 29 // GSVDQ specifies that the Q singular vectors should be computed during 30 // the decomposition. 31 GSVDQ 32 33 // GSVDAll is a convenience value for computing all of the singular vectors. 34 GSVDAll = GSVDU | GSVDV | GSVDQ 35) 36 37// GSVD is a type for creating and using the Generalized Singular Value Decomposition 38// (GSVD) of a matrix. 39// 40// The factorization is a linear transformation of the data sets from the given 41// variable×sample spaces to reduced and diagonalized "eigenvariable"×"eigensample" 42// spaces. 43type GSVD struct { 44 kind GSVDKind 45 46 r, p, c, k, l int 47 s1, s2 []float64 48 a, b, u, v, q blas64.General 49 50 work []float64 51 iwork []int 52} 53 54// succFact returns whether the receiver contains a successful factorization. 55func (gsvd *GSVD) succFact() bool { 56 return gsvd.r != 0 57} 58 59// Factorize computes the generalized singular value decomposition (GSVD) of the input 60// the r×c matrix A and the p×c matrix B. The singular values of A and B are computed 61// in all cases, while the singular vectors are optionally computed depending on the 62// input kind. 63// 64// The full singular value decomposition (kind == GSVDAll) deconstructs A and B as 65// A = U * Σ₁ * [ 0 R ] * Q^T 66// 67// B = V * Σ₂ * [ 0 R ] * Q^T 68// where Σ₁ and Σ₂ are r×(k+l) and p×(k+l) diagonal matrices of singular values, and 69// U, V and Q are r×r, p×p and c×c orthogonal matrices of singular vectors. k+l is the 70// effective numerical rank of the matrix [ A^T B^T ]^T. 71// 72// It is frequently not necessary to compute the full GSVD. Computation time and 73// storage costs can be reduced using the appropriate kind. Either only the singular 74// values can be computed (kind == SVDNone), or in conjunction with specific singular 75// vectors (kind bit set according to matrix.GSVDU, matrix.GSVDV and matrix.GSVDQ). 76// 77// Factorize returns whether the decomposition succeeded. If the decomposition 78// failed, routines that require a successful factorization will panic. 79func (gsvd *GSVD) Factorize(a, b Matrix, kind GSVDKind) (ok bool) { 80 // kill the previous decomposition 81 gsvd.r = 0 82 gsvd.kind = 0 83 84 r, c := a.Dims() 85 gsvd.r, gsvd.c = r, c 86 p, c := b.Dims() 87 gsvd.p = p 88 if gsvd.c != c { 89 panic(ErrShape) 90 } 91 var jobU, jobV, jobQ lapack.GSVDJob 92 switch { 93 default: 94 panic("gsvd: bad input kind") 95 case kind == GSVDNone: 96 jobU = lapack.GSVDNone 97 jobV = lapack.GSVDNone 98 jobQ = lapack.GSVDNone 99 case GSVDAll&kind != 0: 100 if GSVDU&kind != 0 { 101 jobU = lapack.GSVDU 102 gsvd.u = blas64.General{ 103 Rows: r, 104 Cols: r, 105 Stride: r, 106 Data: use(gsvd.u.Data, r*r), 107 } 108 } 109 if GSVDV&kind != 0 { 110 jobV = lapack.GSVDV 111 gsvd.v = blas64.General{ 112 Rows: p, 113 Cols: p, 114 Stride: p, 115 Data: use(gsvd.v.Data, p*p), 116 } 117 } 118 if GSVDQ&kind != 0 { 119 jobQ = lapack.GSVDQ 120 gsvd.q = blas64.General{ 121 Rows: c, 122 Cols: c, 123 Stride: c, 124 Data: use(gsvd.q.Data, c*c), 125 } 126 } 127 } 128 129 // A and B are destroyed on call, so copy the matrices. 130 aCopy := DenseCopyOf(a) 131 bCopy := DenseCopyOf(b) 132 133 gsvd.s1 = use(gsvd.s1, c) 134 gsvd.s2 = use(gsvd.s2, c) 135 136 gsvd.iwork = useInt(gsvd.iwork, c) 137 138 gsvd.work = use(gsvd.work, 1) 139 lapack64.Ggsvd3(jobU, jobV, jobQ, aCopy.mat, bCopy.mat, gsvd.s1, gsvd.s2, gsvd.u, gsvd.v, gsvd.q, gsvd.work, -1, gsvd.iwork) 140 gsvd.work = use(gsvd.work, int(gsvd.work[0])) 141 gsvd.k, gsvd.l, ok = lapack64.Ggsvd3(jobU, jobV, jobQ, aCopy.mat, bCopy.mat, gsvd.s1, gsvd.s2, gsvd.u, gsvd.v, gsvd.q, gsvd.work, len(gsvd.work), gsvd.iwork) 142 if ok { 143 gsvd.a = aCopy.mat 144 gsvd.b = bCopy.mat 145 gsvd.kind = kind 146 } 147 return ok 148} 149 150// Kind returns the GSVDKind of the decomposition. If no decomposition has been 151// computed, Kind returns -1. 152func (gsvd *GSVD) Kind() GSVDKind { 153 if !gsvd.succFact() { 154 return -1 155 } 156 return gsvd.kind 157} 158 159// Rank returns the k and l terms of the rank of [ A^T B^T ]^T. 160func (gsvd *GSVD) Rank() (k, l int) { 161 return gsvd.k, gsvd.l 162} 163 164// GeneralizedValues returns the generalized singular values of the factorized matrices. 165// If the input slice is non-nil, the values will be stored in-place into the slice. 166// In this case, the slice must have length min(r,c)-k, and GeneralizedValues will 167// panic with matrix.ErrSliceLengthMismatch otherwise. If the input slice is nil, 168// a new slice of the appropriate length will be allocated and returned. 169// 170// GeneralizedValues will panic if the receiver does not contain a successful factorization. 171func (gsvd *GSVD) GeneralizedValues(v []float64) []float64 { 172 if !gsvd.succFact() { 173 panic(badFact) 174 } 175 r := gsvd.r 176 c := gsvd.c 177 k := gsvd.k 178 d := min(r, c) 179 if v == nil { 180 v = make([]float64, d-k) 181 } 182 if len(v) != d-k { 183 panic(ErrSliceLengthMismatch) 184 } 185 floats.DivTo(v, gsvd.s1[k:d], gsvd.s2[k:d]) 186 return v 187} 188 189// ValuesA returns the singular values of the factorized A matrix. 190// If the input slice is non-nil, the values will be stored in-place into the slice. 191// In this case, the slice must have length min(r,c)-k, and ValuesA will panic with 192// matrix.ErrSliceLengthMismatch otherwise. If the input slice is nil, 193// a new slice of the appropriate length will be allocated and returned. 194// 195// ValuesA will panic if the receiver does not contain a successful factorization. 196func (gsvd *GSVD) ValuesA(s []float64) []float64 { 197 if !gsvd.succFact() { 198 panic(badFact) 199 } 200 r := gsvd.r 201 c := gsvd.c 202 k := gsvd.k 203 d := min(r, c) 204 if s == nil { 205 s = make([]float64, d-k) 206 } 207 if len(s) != d-k { 208 panic(ErrSliceLengthMismatch) 209 } 210 copy(s, gsvd.s1[k:min(r, c)]) 211 return s 212} 213 214// ValuesB returns the singular values of the factorized B matrix. 215// If the input slice is non-nil, the values will be stored in-place into the slice. 216// In this case, the slice must have length min(r,c)-k, and ValuesB will panic with 217// matrix.ErrSliceLengthMismatch otherwise. If the input slice is nil, 218// a new slice of the appropriate length will be allocated and returned. 219// 220// ValuesB will panic if the receiver does not contain a successful factorization. 221func (gsvd *GSVD) ValuesB(s []float64) []float64 { 222 if !gsvd.succFact() { 223 panic(badFact) 224 } 225 r := gsvd.r 226 c := gsvd.c 227 k := gsvd.k 228 d := min(r, c) 229 if s == nil { 230 s = make([]float64, d-k) 231 } 232 if len(s) != d-k { 233 panic(ErrSliceLengthMismatch) 234 } 235 copy(s, gsvd.s2[k:d]) 236 return s 237} 238 239// ZeroRTo extracts the matrix [ 0 R ] from the singular value decomposition, storing 240// the result in-place into dst. [ 0 R ] is size (k+l)×c. 241// If dst is nil, a new matrix is allocated. The resulting ZeroR matrix is returned. 242// 243// ZeroRTo will panic if the receiver does not contain a successful factorization. 244func (gsvd *GSVD) ZeroRTo(dst *Dense) *Dense { 245 if !gsvd.succFact() { 246 panic(badFact) 247 } 248 r := gsvd.r 249 c := gsvd.c 250 k := gsvd.k 251 l := gsvd.l 252 h := min(k+l, r) 253 if dst == nil { 254 dst = NewDense(k+l, c, nil) 255 } else { 256 dst.reuseAsZeroed(k+l, c) 257 } 258 a := Dense{ 259 mat: gsvd.a, 260 capRows: r, 261 capCols: c, 262 } 263 dst.Slice(0, h, c-k-l, c).(*Dense). 264 Copy(a.Slice(0, h, c-k-l, c)) 265 if r < k+l { 266 b := Dense{ 267 mat: gsvd.b, 268 capRows: gsvd.p, 269 capCols: c, 270 } 271 dst.Slice(r, k+l, c+r-k-l, c).(*Dense). 272 Copy(b.Slice(r-k, l, c+r-k-l, c)) 273 } 274 return dst 275} 276 277// SigmaATo extracts the matrix Σ₁ from the singular value decomposition, storing 278// the result in-place into dst. Σ₁ is size r×(k+l). 279// If dst is nil, a new matrix is allocated. The resulting SigmaA matrix is returned. 280// 281// SigmaATo will panic if the receiver does not contain a successful factorization. 282func (gsvd *GSVD) SigmaATo(dst *Dense) *Dense { 283 if !gsvd.succFact() { 284 panic(badFact) 285 } 286 r := gsvd.r 287 k := gsvd.k 288 l := gsvd.l 289 if dst == nil { 290 dst = NewDense(r, k+l, nil) 291 } else { 292 dst.reuseAsZeroed(r, k+l) 293 } 294 for i := 0; i < k; i++ { 295 dst.set(i, i, 1) 296 } 297 for i := k; i < min(r, k+l); i++ { 298 dst.set(i, i, gsvd.s1[i]) 299 } 300 return dst 301} 302 303// SigmaBTo extracts the matrix Σ₂ from the singular value decomposition, storing 304// the result in-place into dst. Σ₂ is size p×(k+l). 305// If dst is nil, a new matrix is allocated. The resulting SigmaB matrix is returned. 306// 307// SigmaBTo will panic if the receiver does not contain a successful factorization. 308func (gsvd *GSVD) SigmaBTo(dst *Dense) *Dense { 309 if !gsvd.succFact() { 310 panic(badFact) 311 } 312 r := gsvd.r 313 p := gsvd.p 314 k := gsvd.k 315 l := gsvd.l 316 if dst == nil { 317 dst = NewDense(p, k+l, nil) 318 } else { 319 dst.reuseAsZeroed(p, k+l) 320 } 321 for i := 0; i < min(l, r-k); i++ { 322 dst.set(i, i+k, gsvd.s2[k+i]) 323 } 324 for i := r - k; i < l; i++ { 325 dst.set(i, i+k, 1) 326 } 327 return dst 328} 329 330// UTo extracts the matrix U from the singular value decomposition, storing 331// the result in-place into dst. U is size r×r. 332// If dst is nil, a new matrix is allocated. The resulting U matrix is returned. 333// 334// UTo will panic if the receiver does not contain a successful factorization. 335func (gsvd *GSVD) UTo(dst *Dense) *Dense { 336 if !gsvd.succFact() { 337 panic(badFact) 338 } 339 if gsvd.kind&GSVDU == 0 { 340 panic("mat: improper GSVD kind") 341 } 342 r := gsvd.u.Rows 343 c := gsvd.u.Cols 344 if dst == nil { 345 dst = NewDense(r, c, nil) 346 } else { 347 dst.reuseAs(r, c) 348 } 349 350 tmp := &Dense{ 351 mat: gsvd.u, 352 capRows: r, 353 capCols: c, 354 } 355 dst.Copy(tmp) 356 return dst 357} 358 359// VTo extracts the matrix V from the singular value decomposition, storing 360// the result in-place into dst. V is size p×p. 361// If dst is nil, a new matrix is allocated. The resulting V matrix is returned. 362// 363// VTo will panic if the receiver does not contain a successful factorization. 364func (gsvd *GSVD) VTo(dst *Dense) *Dense { 365 if !gsvd.succFact() { 366 panic(badFact) 367 } 368 if gsvd.kind&GSVDV == 0 { 369 panic("mat: improper GSVD kind") 370 } 371 r := gsvd.v.Rows 372 c := gsvd.v.Cols 373 if dst == nil { 374 dst = NewDense(r, c, nil) 375 } else { 376 dst.reuseAs(r, c) 377 } 378 379 tmp := &Dense{ 380 mat: gsvd.v, 381 capRows: r, 382 capCols: c, 383 } 384 dst.Copy(tmp) 385 return dst 386} 387 388// QTo extracts the matrix Q from the singular value decomposition, storing 389// the result in-place into dst. Q is size c×c. 390// If dst is nil, a new matrix is allocated. The resulting Q matrix is returned. 391// 392// QTo will panic if the receiver does not contain a successful factorization. 393func (gsvd *GSVD) QTo(dst *Dense) *Dense { 394 if !gsvd.succFact() { 395 panic(badFact) 396 } 397 if gsvd.kind&GSVDQ == 0 { 398 panic("mat: improper GSVD kind") 399 } 400 r := gsvd.q.Rows 401 c := gsvd.q.Cols 402 if dst == nil { 403 dst = NewDense(r, c, nil) 404 } else { 405 dst.reuseAs(r, c) 406 } 407 408 tmp := &Dense{ 409 mat: gsvd.q, 410 capRows: r, 411 capCols: c, 412 } 413 dst.Copy(tmp) 414 return dst 415} 416