123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248 |
- // Copyright ©2013 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/lapack"
- "gonum.org/v1/gonum/lapack/lapack64"
- )
- // SVD is a type for creating and using the Singular Value Decomposition (SVD)
- // of a matrix.
- type SVD struct {
- kind SVDKind
- s []float64
- u blas64.General
- vt blas64.General
- }
- // SVDKind specifies the treatment of singular vectors during an SVD
- // factorization.
- type SVDKind int
- const (
- // SVDNone specifies that no singular vectors should be computed during
- // the decomposition.
- SVDNone SVDKind = 0
- // SVDThinU specifies the thin decomposition for U should be computed.
- SVDThinU SVDKind = 1 << (iota - 1)
- // SVDFullU specifies the full decomposition for U should be computed.
- SVDFullU
- // SVDThinV specifies the thin decomposition for V should be computed.
- SVDThinV
- // SVDFullV specifies the full decomposition for V should be computed.
- SVDFullV
- // SVDThin is a convenience value for computing both thin vectors.
- SVDThin SVDKind = SVDThinU | SVDThinV
- // SVDThin is a convenience value for computing both full vectors.
- SVDFull SVDKind = SVDFullU | SVDFullV
- )
- // succFact returns whether the receiver contains a successful factorization.
- func (svd *SVD) succFact() bool {
- return len(svd.s) != 0
- }
- // Factorize computes the singular value decomposition (SVD) of the input matrix A.
- // The singular values of A are computed in all cases, while the singular
- // vectors are optionally computed depending on the input kind.
- //
- // The full singular value decomposition (kind == SVDFull) is a factorization
- // of an m×n matrix A of the form
- // A = U * Σ * V^T
- // where Σ is an m×n diagonal matrix, U is an m×m orthogonal matrix, and V is an
- // n×n orthogonal matrix. The diagonal elements of Σ are the singular values of A.
- // The first min(m,n) columns of U and V are, respectively, the left and right
- // singular vectors of A.
- //
- // Significant storage space can be saved by using the thin representation of
- // the SVD (kind == SVDThin) instead of the full SVD, especially if
- // m >> n or m << n. The thin SVD finds
- // A = U~ * Σ * V~^T
- // where U~ is of size m×min(m,n), Σ is a diagonal matrix of size min(m,n)×min(m,n)
- // and V~ is of size n×min(m,n).
- //
- // Factorize returns whether the decomposition succeeded. If the decomposition
- // failed, routines that require a successful factorization will panic.
- func (svd *SVD) Factorize(a Matrix, kind SVDKind) (ok bool) {
- // kill previous factorization
- svd.s = svd.s[:0]
- svd.kind = kind
- m, n := a.Dims()
- var jobU, jobVT lapack.SVDJob
- // TODO(btracey): This code should be modified to have the smaller
- // matrix written in-place into aCopy when the lapack/native/dgesvd
- // implementation is complete.
- switch {
- case kind&SVDFullU != 0:
- jobU = lapack.SVDAll
- svd.u = blas64.General{
- Rows: m,
- Cols: m,
- Stride: m,
- Data: use(svd.u.Data, m*m),
- }
- case kind&SVDThinU != 0:
- jobU = lapack.SVDStore
- svd.u = blas64.General{
- Rows: m,
- Cols: min(m, n),
- Stride: min(m, n),
- Data: use(svd.u.Data, m*min(m, n)),
- }
- default:
- jobU = lapack.SVDNone
- }
- switch {
- case kind&SVDFullV != 0:
- svd.vt = blas64.General{
- Rows: n,
- Cols: n,
- Stride: n,
- Data: use(svd.vt.Data, n*n),
- }
- jobVT = lapack.SVDAll
- case kind&SVDThinV != 0:
- svd.vt = blas64.General{
- Rows: min(m, n),
- Cols: n,
- Stride: n,
- Data: use(svd.vt.Data, min(m, n)*n),
- }
- jobVT = lapack.SVDStore
- default:
- jobVT = lapack.SVDNone
- }
- // A is destroyed on call, so copy the matrix.
- aCopy := DenseCopyOf(a)
- svd.kind = kind
- svd.s = use(svd.s, min(m, n))
- work := []float64{0}
- lapack64.Gesvd(jobU, jobVT, aCopy.mat, svd.u, svd.vt, svd.s, work, -1)
- work = getFloats(int(work[0]), false)
- ok = lapack64.Gesvd(jobU, jobVT, aCopy.mat, svd.u, svd.vt, svd.s, work, len(work))
- putFloats(work)
- if !ok {
- svd.kind = 0
- }
- return ok
- }
- // Kind returns the SVDKind of the decomposition. If no decomposition has been
- // computed, Kind returns -1.
- func (svd *SVD) Kind() SVDKind {
- if !svd.succFact() {
- return -1
- }
- return svd.kind
- }
- // Cond returns the 2-norm condition number for the factorized matrix. Cond will
- // panic if the receiver does not contain a successful factorization.
- func (svd *SVD) Cond() float64 {
- if !svd.succFact() {
- panic(badFact)
- }
- return svd.s[0] / svd.s[len(svd.s)-1]
- }
- // Values returns the singular values of the factorized matrix in descending order.
- //
- // 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(m,n), and Values will
- // panic with ErrSliceLengthMismatch otherwise. If the input slice is nil, a new
- // slice of the appropriate length will be allocated and returned.
- //
- // Values will panic if the receiver does not contain a successful factorization.
- func (svd *SVD) Values(s []float64) []float64 {
- if !svd.succFact() {
- panic(badFact)
- }
- if s == nil {
- s = make([]float64, len(svd.s))
- }
- if len(s) != len(svd.s) {
- panic(ErrSliceLengthMismatch)
- }
- copy(s, svd.s)
- return s
- }
- // UTo extracts the matrix U from the singular value decomposition. The first
- // min(m,n) columns are the left singular vectors and correspond to the singular
- // values as returned from SVD.Values.
- //
- // If dst is not nil, U is stored in-place into dst, and dst must have size
- // m×m if the full U was computed, size m×min(m,n) if the thin U was computed,
- // and UTo panics otherwise. If dst is nil, a new matrix of the appropriate size
- // is allocated and returned.
- func (svd *SVD) UTo(dst *Dense) *Dense {
- if !svd.succFact() {
- panic(badFact)
- }
- kind := svd.kind
- if kind&SVDThinU == 0 && kind&SVDFullU == 0 {
- panic("svd: u not computed during factorization")
- }
- r := svd.u.Rows
- c := svd.u.Cols
- if dst == nil {
- dst = NewDense(r, c, nil)
- } else {
- dst.reuseAs(r, c)
- }
- tmp := &Dense{
- mat: svd.u,
- capRows: r,
- capCols: c,
- }
- dst.Copy(tmp)
- return dst
- }
- // VTo extracts the matrix V from the singular value decomposition. The first
- // min(m,n) columns are the right singular vectors and correspond to the singular
- // values as returned from SVD.Values.
- //
- // If dst is not nil, V is stored in-place into dst, and dst must have size
- // n×n if the full V was computed, size n×min(m,n) if the thin V was computed,
- // and VTo panics otherwise. If dst is nil, a new matrix of the appropriate size
- // is allocated and returned.
- func (svd *SVD) VTo(dst *Dense) *Dense {
- if !svd.succFact() {
- panic(badFact)
- }
- kind := svd.kind
- if kind&SVDThinU == 0 && kind&SVDFullV == 0 {
- panic("svd: v not computed during factorization")
- }
- r := svd.vt.Rows
- c := svd.vt.Cols
- if dst == nil {
- dst = NewDense(c, r, nil)
- } else {
- dst.reuseAs(c, r)
- }
- tmp := &Dense{
- mat: svd.vt,
- capRows: r,
- capCols: c,
- }
- dst.Copy(tmp.T())
- return dst
- }
|