+++ /dev/null
-// 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
-}
OSDN Git Service
