blob: 663a07299bcbd8fc4d3947ce9b79b51f8d925b60 [file] [log] [blame]
 // Copyright ©2017 The gonum Authors. All rights reserved. // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. package mat import ( "gonum.org/v1/gonum/blas/blas64" "gonum.org/v1/gonum/floats" "gonum.org/v1/gonum/lapack" "gonum.org/v1/gonum/lapack/lapack64" ) // GSVD is a type for creating and using the Generalized Singular Value Decomposition // (GSVD) of a matrix. // // The factorization is a linear transformation of the data sets from the given // variable×sample spaces to reduced and diagonalized "eigenvariable"×"eigensample" // spaces. type GSVD struct { kind GSVDKind r, p, c, k, l int s1, s2 []float64 a, b, u, v, q blas64.General work []float64 iwork []int } // Factorize computes the generalized singular value decomposition (GSVD) of the input // the r×c matrix A and the p×c matrix B. The singular values of A and B are computed // in all cases, while the singular vectors are optionally computed depending on the // input kind. // // The full singular value decomposition (kind == GSVDU|GSVDV|GSVDQ) deconstructs A and B as // A = U * Σ₁ * [ 0 R ] * Q^T // // B = V * Σ₂ * [ 0 R ] * Q^T // where Σ₁ and Σ₂ are r×(k+l) and p×(k+l) diagonal matrices of singular values, and // U, V and Q are r×r, p×p and c×c orthogonal matrices of singular vectors. k+l is the // effective numerical rank of the matrix [ A^T B^T ]^T. // // It is frequently not necessary to compute the full GSVD. Computation time and // storage costs can be reduced using the appropriate kind. Either only the singular // values can be computed (kind == SVDNone), or in conjunction with specific singular // vectors (kind bit set according to matrix.GSVDU, matrix.GSVDV and matrix.GSVDQ). // // Factorize returns whether the decomposition succeeded. If the decomposition // failed, routines that require a successful factorization will panic. func (gsvd *GSVD) Factorize(a, b Matrix, kind GSVDKind) (ok bool) { r, c := a.Dims() gsvd.r, gsvd.c = r, c p, c := b.Dims() gsvd.p = p if gsvd.c != c { panic(ErrShape) } var jobU, jobV, jobQ lapack.GSVDJob switch { default: panic("gsvd: bad input kind") case kind == GSVDNone: jobU = lapack.GSVDNone jobV = lapack.GSVDNone jobQ = lapack.GSVDNone case (GSVDU|GSVDV|GSVDQ)&kind != 0: if GSVDU&kind != 0 { jobU = lapack.GSVDU gsvd.u = blas64.General{ Rows: r, Cols: r, Stride: r, Data: use(gsvd.u.Data, r*r), } } if GSVDV&kind != 0 { jobV = lapack.GSVDV gsvd.v = blas64.General{ Rows: p, Cols: p, Stride: p, Data: use(gsvd.v.Data, p*p), } } if GSVDQ&kind != 0 { jobQ = lapack.GSVDQ gsvd.q = blas64.General{ Rows: c, Cols: c, Stride: c, Data: use(gsvd.q.Data, c*c), } } } // A and B are destroyed on call, so copy the matrices. aCopy := DenseCopyOf(a) bCopy := DenseCopyOf(b) gsvd.s1 = use(gsvd.s1, c) gsvd.s2 = use(gsvd.s2, c) gsvd.iwork = useInt(gsvd.iwork, c) gsvd.work = use(gsvd.work, 1) lapack64.Ggsvd3(jobU, jobV, jobQ, aCopy.mat, bCopy.mat, gsvd.s1, gsvd.s2, gsvd.u, gsvd.v, gsvd.q, gsvd.work, -1, gsvd.iwork) gsvd.work = use(gsvd.work, int(gsvd.work[0])) 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) if ok { gsvd.a = aCopy.mat gsvd.b = bCopy.mat gsvd.kind = kind } return ok } // Kind returns the matrix.GSVDKind of the decomposition. If no decomposition has been // computed, Kind returns 0. func (gsvd *GSVD) Kind() GSVDKind { return gsvd.kind } // Rank returns the k and l terms of the rank of [ A^T B^T ]^T. func (gsvd *GSVD) Rank() (k, l int) { return gsvd.k, gsvd.l } // GeneralizedValues returns the generalized singular values of the factorized matrices. // If the input slice is non-nil, the values will be stored in-place into the slice. // In this case, the slice must have length min(r,c)-k, and GeneralizedValues will // panic with matrix.ErrSliceLengthMismatch otherwise. If the input slice is nil, // a new slice of the appropriate length will be allocated and returned. // // GeneralizedValues will panic if the receiver does not contain a successful factorization. func (gsvd *GSVD) GeneralizedValues(v []float64) []float64 { if gsvd.kind == 0 { panic("gsvd: no decomposition computed") } r := gsvd.r c := gsvd.c k := gsvd.k d := min(r, c) if v == nil { v = make([]float64, d-k) } if len(v) != d-k { panic(ErrSliceLengthMismatch) } floats.DivTo(v, gsvd.s1[k:d], gsvd.s2[k:d]) return v } // ValuesA returns the singular values of the factorized A matrix. // If the input slice is non-nil, the values will be stored in-place into the slice. // In this case, the slice must have length min(r,c)-k, and ValuesA will panic with // matrix.ErrSliceLengthMismatch otherwise. If the input slice is nil, // a new slice of the appropriate length will be allocated and returned. // // ValuesA will panic if the receiver does not contain a successful factorization. func (gsvd *GSVD) ValuesA(s []float64) []float64 { if gsvd.kind == 0 { panic("gsvd: no decomposition computed") } r := gsvd.r c := gsvd.c k := gsvd.k d := min(r, c) if s == nil { s = make([]float64, d-k) } if len(s) != d-k { panic(ErrSliceLengthMismatch) } copy(s, gsvd.s1[k:min(r, c)]) return s } // ValuesB returns the singular values of the factorized B matrix. // If the input slice is non-nil, the values will be stored in-place into the slice. // In this case, the slice must have length min(r,c)-k, and ValuesB will panic with // matrix.ErrSliceLengthMismatch otherwise. If the input slice is nil, // a new slice of the appropriate length will be allocated and returned. // // ValuesB will panic if the receiver does not contain a successful factorization. func (gsvd *GSVD) ValuesB(s []float64) []float64 { if gsvd.kind == 0 { panic("gsvd: no decomposition computed") } r := gsvd.r c := gsvd.c k := gsvd.k d := min(r, c) if s == nil { s = make([]float64, d-k) } if len(s) != d-k { panic(ErrSliceLengthMismatch) } copy(s, gsvd.s2[k:d]) return s } // ZeroRTo extracts the matrix [ 0 R ] from the singular value decomposition, storing // the result in-place into dst. [ 0 R ] is size (k+l)×c. // If dst is nil, a new matrix is allocated. The resulting ZeroR matrix is returned. // // ZeroRTo will panic if the receiver does not contain a successful factorization. func (gsvd *GSVD) ZeroRTo(dst *Dense) *Dense { if gsvd.kind == 0 { panic("gsvd: no decomposition computed") } r := gsvd.r c := gsvd.c k := gsvd.k l := gsvd.l h := min(k+l, r) if dst == nil { dst = NewDense(k+l, c, nil) } else { dst.reuseAsZeroed(k+l, c) } a := Dense{ mat: gsvd.a, capRows: r, capCols: c, } dst.Slice(0, h, c-k-l, c).(*Dense). Copy(a.Slice(0, h, c-k-l, c)) if r < k+l { b := Dense{ mat: gsvd.b, capRows: gsvd.p, capCols: c, } dst.Slice(r, k+l, c+r-k-l, c).(*Dense). Copy(b.Slice(r-k, l, c+r-k-l, c)) } return dst } // SigmaATo extracts the matrix Σ₁ from the singular value decomposition, storing // the result in-place into dst. Σ₁ is size r×(k+l). // If dst is nil, a new matrix is allocated. The resulting SigmaA matrix is returned. // // SigmaATo will panic if the receiver does not contain a successful factorization. func (gsvd *GSVD) SigmaATo(dst *Dense) *Dense { if gsvd.kind == 0 { panic("gsvd: no decomposition computed") } r := gsvd.r k := gsvd.k l := gsvd.l if dst == nil { dst = NewDense(r, k+l, nil) } else { dst.reuseAsZeroed(r, k+l) } for i := 0; i < k; i++ { dst.set(i, i, 1) } for i := k; i < min(r, k+l); i++ { dst.set(i, i, gsvd.s1[i]) } return dst } // SigmaBTo extracts the matrix Σ₂ from the singular value decomposition, storing // the result in-place into dst. Σ₂ is size p×(k+l). // If dst is nil, a new matrix is allocated. The resulting SigmaB matrix is returned. // // SigmaBTo will panic if the receiver does not contain a successful factorization. func (gsvd *GSVD) SigmaBTo(dst *Dense) *Dense { if gsvd.kind == 0 { panic("gsvd: no decomposition computed") } r := gsvd.r p := gsvd.p k := gsvd.k l := gsvd.l if dst == nil { dst = NewDense(p, k+l, nil) } else { dst.reuseAsZeroed(p, k+l) } for i := 0; i < min(l, r-k); i++ { dst.set(i, i+k, gsvd.s2[k+i]) } for i := r - k; i < l; i++ { dst.set(i, i+k, 1) } return dst } // UTo extracts the matrix U from the singular value decomposition, storing // the result in-place into dst. U is size r×r. // If dst is nil, a new matrix is allocated. The resulting U matrix is returned. // // UTo will panic if the receiver does not contain a successful factorization. func (gsvd *GSVD) UTo(dst *Dense) *Dense { if gsvd.kind&GSVDU == 0 { panic("mat: improper GSVD kind") } r := gsvd.u.Rows c := gsvd.u.Cols if dst == nil { dst = NewDense(r, c, nil) } else { dst.reuseAs(r, c) } tmp := &Dense{ mat: gsvd.u, capRows: r, capCols: c, } dst.Copy(tmp) return dst } // VTo extracts the matrix V from the singular value decomposition, storing // the result in-place into dst. V is size p×p. // If dst is nil, a new matrix is allocated. The resulting V matrix is returned. // // VTo will panic if the receiver does not contain a successful factorization. func (gsvd *GSVD) VTo(dst *Dense) *Dense { if gsvd.kind&GSVDV == 0 { panic("mat: improper GSVD kind") } r := gsvd.v.Rows c := gsvd.v.Cols if dst == nil { dst = NewDense(r, c, nil) } else { dst.reuseAs(r, c) } tmp := &Dense{ mat: gsvd.v, capRows: r, capCols: c, } dst.Copy(tmp) return dst } // QTo extracts the matrix Q from the singular value decomposition, storing // the result in-place into dst. Q is size c×c. // If dst is nil, a new matrix is allocated. The resulting Q matrix is returned. // // QTo will panic if the receiver does not contain a successful factorization. func (gsvd *GSVD) QTo(dst *Dense) *Dense { if gsvd.kind&GSVDQ == 0 { panic("mat: improper GSVD kind") } r := gsvd.q.Rows c := gsvd.q.Cols if dst == nil { dst = NewDense(r, c, nil) } else { dst.reuseAs(r, c) } tmp := &Dense{ mat: gsvd.q, capRows: r, capCols: c, } dst.Copy(tmp) return dst }