From 1d8caee34e2bf68129643aa1cac98e2ace88b29a Mon Sep 17 00:00:00 2001 From: Brendan Tracey Date: Sat, 23 Mar 2019 17:20:14 +0000 Subject: [PATCH] mat: change factorization inputs to use bit types (#872) * mat: change factorization inputs to use bit types Fixes #756 and #748. --- mat/consts.go | 39 ----------- mat/eigen.go | 48 +++++++++++--- mat/eigen_example_test.go | 2 +- mat/eigen_test.go | 8 +-- mat/gsvd.go | 76 ++++++++++++++++----- mat/hogsvd.go | 19 ++++-- mat/svd.go | 136 +++++++++++++++++++++++++------------- 7 files changed, 205 insertions(+), 123 deletions(-) diff --git a/mat/consts.go b/mat/consts.go index a7b6370d..3de3f5bf 100644 --- a/mat/consts.go +++ b/mat/consts.go @@ -13,42 +13,3 @@ const ( // Lower specifies a lower triangular matrix. Lower TriKind = false ) - -// 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 = iota + 1 - // SVDThin computes the thin singular vectors, that is, it computes - // 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). - SVDThin - // SVDFull computes the full singular value decomposition, - // A = U * Σ * V^T - // where U is of size m×m, Σ is an m×n diagonal matrix, and V is an n×n matrix. - SVDFull -) - -// GSVDKind specifies the treatment of singular vectors during a GSVD -// factorization. -type GSVDKind int - -const ( - // GSVDU specifies that the U singular vectors should be computed during - // the decomposition. - GSVDU GSVDKind = 1 << iota - // GSVDV specifies that the V singular vectors should be computed during - // the decomposition. - GSVDV - // GSVDQ specifies that the Q singular vectors should be computed during - // the decomposition. - GSVDQ - - // GSVDNone specifies that no singular vector should be computed during - // the decomposition. - GSVDNone -) diff --git a/mat/eigen.go b/mat/eigen.go index 90a94ae9..825f459c 100644 --- a/mat/eigen.go +++ b/mat/eigen.go @@ -104,12 +104,25 @@ func (m *Dense) EigenvectorsSym(e *EigenSym) { m.Copy(e.vectors) } +// EigenKind specifies the computation of eigenvectors during factorization. +type EigenKind int + +const ( + // EigenNone specifies to not compute any eigenvectors. + EigenNone EigenKind = 0 + // EigenLeft specifies to compute the left eigenvectors. + EigenLeft EigenKind = 1 << iota + // EigenRight specifies to compute the right eigenvectors. + EigenRight + // EigenBoth is a convenience value for computing both eigenvectors. + EigenBoth EigenKind = EigenLeft | EigenRight +) + // Eigen is a type for creating and using the eigenvalue decomposition of a dense matrix. type Eigen struct { n int // The size of the factorized matrix. - right bool // have the right eigenvectors been computed - left bool // have the left eigenvectors been computed + kind EigenKind values []complex128 rVectors *CDense @@ -118,7 +131,7 @@ type Eigen struct { // succFact returns whether the receiver contains a successful factorization. func (e *Eigen) succFact() bool { - return len(e.values) != 0 + return e.n != 0 } // Factorize computes the eigenvalues of the square matrix a, and optionally @@ -135,13 +148,17 @@ func (e *Eigen) succFact() bool { // // Typically eigenvectors refer to right eigenvectors. // -// In all cases, Factorize computes the eigenvalues of the matrix. If right and left -// are true, then the right and left eigenvectors will be computed, respectively. +// In all cases, Factorize computes the eigenvalues of the matrix. kind +// specifies which of the eigenvectors, if any, to compute. See the EigenKind +// documentation for more information. // Eigen panics if the input matrix is not square. // // Factorize returns whether the decomposition succeeded. If the decomposition // failed, methods that require a successful factorization will panic. -func (e *Eigen) Factorize(a Matrix, left, right bool) (ok bool) { +func (e *Eigen) Factorize(a Matrix, kind EigenKind) (ok bool) { + // kill previous factorization. + e.n = 0 + e.kind = 0 // Copy a because it is modified during the Lapack call. r, c := a.Dims() if r != c { @@ -150,6 +167,9 @@ func (e *Eigen) Factorize(a Matrix, left, right bool) (ok bool) { var sd Dense sd.Clone(a) + left := kind&EigenLeft != 0 + right := kind&EigenRight != 0 + var vl, vr Dense jobvl := lapack.LeftEVNone jobvr := lapack.RightEVNone @@ -183,8 +203,7 @@ func (e *Eigen) Factorize(a Matrix, left, right bool) (ok bool) { return false } e.n = r - e.right = right - e.left = left + e.kind = kind // Construct complex eigenvalues from float64 data. values := make([]complex128, r) @@ -212,6 +231,15 @@ func (e *Eigen) Factorize(a Matrix, left, right bool) (ok bool) { return true } +// Kind returns the EigenKind of the decomposition. If no decomposition has been +// computed, Kind returns -1. +func (e *Eigen) Kind() EigenKind { + if !e.succFact() { + return -1 + } + return e.kind +} + // Values extracts the eigenvalues of the factorized matrix. If dst is // non-nil, the values are stored in-place into dst. In this case // dst must have length n, otherwise Values will panic. If dst is @@ -281,7 +309,7 @@ func (e *Eigen) VectorsTo(dst *CDense) *CDense { if !e.succFact() { panic(badFact) } - if !e.right { + if e.kind&EigenRight == 0 { panic(badNoVect) } if dst == nil { @@ -303,7 +331,7 @@ func (e *Eigen) LeftVectorsTo(dst *CDense) *CDense { if !e.succFact() { panic(badFact) } - if !e.left { + if e.kind&EigenLeft == 0 { panic(badNoVect) } if dst == nil { diff --git a/mat/eigen_example_test.go b/mat/eigen_example_test.go index 6f95bee6..a0175be2 100644 --- a/mat/eigen_example_test.go +++ b/mat/eigen_example_test.go @@ -50,7 +50,7 @@ func ExampleEigen() { fmt.Printf("A = %v\n\n", mat.Formatted(a, mat.Prefix(" "))) var eig mat.Eigen - ok := eig.Factorize(a, true, false) + ok := eig.Factorize(a, mat.EigenLeft) if !ok { log.Fatal("Eigendecomposition failed") } diff --git a/mat/eigen_test.go b/mat/eigen_test.go index e221875b..5dbc61ba 100644 --- a/mat/eigen_test.go +++ b/mat/eigen_test.go @@ -63,13 +63,13 @@ func TestEigen(t *testing.T) { }, } { var e1, e2, e3, e4 Eigen - ok := e1.Factorize(test.a, true, true) + ok := e1.Factorize(test.a, EigenBoth) if !ok { panic("bad factorization") } - e2.Factorize(test.a, false, true) - e3.Factorize(test.a, true, false) - e4.Factorize(test.a, false, false) + e2.Factorize(test.a, EigenRight) + e3.Factorize(test.a, EigenLeft) + e4.Factorize(test.a, EigenNone) v1 := e1.Values(nil) if !cmplxEqualTol(v1, test.values, 1e-14) { diff --git a/mat/gsvd.go b/mat/gsvd.go index f2b82cb5..2de511a9 100644 --- a/mat/gsvd.go +++ b/mat/gsvd.go @@ -11,6 +11,29 @@ import ( "gonum.org/v1/gonum/lapack/lapack64" ) +// GSVDKind specifies the treatment of singular vectors during a GSVD +// factorization. +type GSVDKind int + +const ( + // GSVDNone specifies that no singular vectors should be computed during + // the decomposition. + GSVDNone GSVDKind = 0 + + // GSVDU specifies that the U singular vectors should be computed during + // the decomposition. + GSVDU GSVDKind = 1 << iota + // GSVDV specifies that the V singular vectors should be computed during + // the decomposition. + GSVDV + // GSVDQ specifies that the Q singular vectors should be computed during + // the decomposition. + GSVDQ + + // GSVDAll is a convenience value for computing all of the singular vectors. + GSVDAll = GSVDU | GSVDV | GSVDQ +) + // GSVD is a type for creating and using the Generalized Singular Value Decomposition // (GSVD) of a matrix. // @@ -28,12 +51,17 @@ type GSVD struct { iwork []int } +// succFact returns whether the receiver contains a successful factorization. +func (gsvd *GSVD) succFact() bool { + return gsvd.r != 0 +} + // 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 +// The full singular value decomposition (kind == GSVDAll) deconstructs A and B as // A = U * Σ₁ * [ 0 R ] * Q^T // // B = V * Σ₂ * [ 0 R ] * Q^T @@ -49,6 +77,10 @@ type GSVD struct { // 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) { + // kill the previous decomposition + gsvd.r = 0 + gsvd.kind = 0 + r, c := a.Dims() gsvd.r, gsvd.c = r, c p, c := b.Dims() @@ -64,7 +96,7 @@ func (gsvd *GSVD) Factorize(a, b Matrix, kind GSVDKind) (ok bool) { jobU = lapack.GSVDNone jobV = lapack.GSVDNone jobQ = lapack.GSVDNone - case (GSVDU|GSVDV|GSVDQ)&kind != 0: + case GSVDAll&kind != 0: if GSVDU&kind != 0 { jobU = lapack.GSVDU gsvd.u = blas64.General{ @@ -115,9 +147,12 @@ func (gsvd *GSVD) Factorize(a, b Matrix, kind GSVDKind) (ok bool) { return ok } -// Kind returns the matrix.GSVDKind of the decomposition. If no decomposition has been -// computed, Kind returns 0. +// Kind returns the GSVDKind of the decomposition. If no decomposition has been +// computed, Kind returns -1. func (gsvd *GSVD) Kind() GSVDKind { + if !gsvd.succFact() { + return -1 + } return gsvd.kind } @@ -134,8 +169,8 @@ func (gsvd *GSVD) Rank() (k, l int) { // // 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") + if !gsvd.succFact() { + panic(badFact) } r := gsvd.r c := gsvd.c @@ -159,8 +194,8 @@ func (gsvd *GSVD) GeneralizedValues(v []float64) []float64 { // // 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") + if !gsvd.succFact() { + panic(badFact) } r := gsvd.r c := gsvd.c @@ -184,8 +219,8 @@ func (gsvd *GSVD) ValuesA(s []float64) []float64 { // // 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") + if !gsvd.succFact() { + panic(badFact) } r := gsvd.r c := gsvd.c @@ -207,8 +242,8 @@ func (gsvd *GSVD) ValuesB(s []float64) []float64 { // // 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") + if !gsvd.succFact() { + panic(badFact) } r := gsvd.r c := gsvd.c @@ -245,8 +280,8 @@ func (gsvd *GSVD) ZeroRTo(dst *Dense) *Dense { // // 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") + if !gsvd.succFact() { + panic(badFact) } r := gsvd.r k := gsvd.k @@ -271,8 +306,8 @@ func (gsvd *GSVD) SigmaATo(dst *Dense) *Dense { // // 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") + if !gsvd.succFact() { + panic(badFact) } r := gsvd.r p := gsvd.p @@ -298,6 +333,9 @@ func (gsvd *GSVD) SigmaBTo(dst *Dense) *Dense { // // UTo will panic if the receiver does not contain a successful factorization. func (gsvd *GSVD) UTo(dst *Dense) *Dense { + if !gsvd.succFact() { + panic(badFact) + } if gsvd.kind&GSVDU == 0 { panic("mat: improper GSVD kind") } @@ -324,6 +362,9 @@ func (gsvd *GSVD) UTo(dst *Dense) *Dense { // // VTo will panic if the receiver does not contain a successful factorization. func (gsvd *GSVD) VTo(dst *Dense) *Dense { + if !gsvd.succFact() { + panic(badFact) + } if gsvd.kind&GSVDV == 0 { panic("mat: improper GSVD kind") } @@ -350,6 +391,9 @@ func (gsvd *GSVD) VTo(dst *Dense) *Dense { // // QTo will panic if the receiver does not contain a successful factorization. func (gsvd *GSVD) QTo(dst *Dense) *Dense { + if !gsvd.succFact() { + panic(badFact) + } if gsvd.kind&GSVDQ == 0 { panic("mat: improper GSVD kind") } diff --git a/mat/hogsvd.go b/mat/hogsvd.go index 899ecff5..2f710f40 100644 --- a/mat/hogsvd.go +++ b/mat/hogsvd.go @@ -24,6 +24,11 @@ type HOGSVD struct { err error } +// succFact returns whether the receiver contains a successful factorization. +func (gsvd *HOGSVD) succFact() bool { + return gsvd.n != 0 +} + // Factorize computes the higher order generalized singular value decomposition (HOGSVD) // of the n input r_i×c column tall matrices in m. HOGSV extends the GSVD case from 2 to n // input matrices. @@ -96,7 +101,7 @@ func (gsvd *HOGSVD) Factorize(m ...Matrix) (ok bool) { s.Scale(1/float64(len(m)*(len(m)-1)), s) var eig Eigen - ok = eig.Factorize(s.T(), false, true) + ok = eig.Factorize(s.T(), EigenRight) if !ok { gsvd.err = errors.New("hogsvd: eigen decomposition failed") return false @@ -157,8 +162,8 @@ func (gsvd *HOGSVD) Len() int { // // UTo will panic if the receiver does not contain a successful factorization. func (gsvd *HOGSVD) UTo(dst *Dense, n int) *Dense { - if gsvd.n == 0 { - panic("hogsvd: unsuccessful factorization") + if !gsvd.succFact() { + panic(badFact) } if n < 0 || gsvd.n <= n { panic("hogsvd: invalid index") @@ -187,8 +192,8 @@ func (gsvd *HOGSVD) UTo(dst *Dense, n int) *Dense { // // Values will panic if the receiver does not contain a successful factorization. func (gsvd *HOGSVD) Values(s []float64, n int) []float64 { - if gsvd.n == 0 { - panic("hogsvd: unsuccessful factorization") + if !gsvd.succFact() { + panic(badFact) } if n < 0 || gsvd.n <= n { panic("hogsvd: invalid index") @@ -214,8 +219,8 @@ func (gsvd *HOGSVD) Values(s []float64, n int) []float64 { // // VTo will panic if the receiver does not contain a successful factorization. func (gsvd *HOGSVD) VTo(dst *Dense) *Dense { - if gsvd.n == 0 { - panic("hogsvd: unsuccessful factorization") + if !gsvd.succFact() { + panic(badFact) } if dst == nil { r, c := gsvd.v.Dims() diff --git a/mat/svd.go b/mat/svd.go index 7b4de1dd..17c96588 100644 --- a/mat/svd.go +++ b/mat/svd.go @@ -20,6 +20,35 @@ type SVD struct { 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. @@ -32,61 +61,67 @@ type SVD struct { // The first min(m,n) columns of U and V are, respectively, the left and right // singular vectors of A. // -// It is frequently not necessary to compute the full SVD. Computation time and -// storage costs can be reduced using the appropriate kind. Only the singular -// values can be computed (kind == SVDNone), or a "thin" representation of the -// orthogonal matrices U and V (kind = SVDThin). The thin representation can -// save a significant amount of memory if m >> n or m << n. +// 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 - switch kind { - default: - panic("svd: bad input kind") - case SVDNone: - svd.u.Stride = 1 - svd.vt.Stride = 1 - jobU = lapack.SVDNone - jobVT = lapack.SVDNone - case SVDFull: - // 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. + + // 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), } - svd.vt = blas64.General{ - Rows: n, - Cols: n, - Stride: n, - Data: use(svd.vt.Data, n*n), - } - jobU = lapack.SVDAll - jobVT = lapack.SVDAll - case SVDThin: - // TODO(btracey): This code should be modified to have the larger - // matrix written in-place into aCopy when the lapack/native/dgesvd - // implementation is complete. + 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: + svd.u.Stride = 1 + 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), } - jobU = lapack.SVDStore jobVT = lapack.SVDStore + default: + svd.vt.Stride = 1 + jobVT = lapack.SVDNone } // A is destroyed on call, so copy the matrix. @@ -105,17 +140,20 @@ func (svd *SVD) Factorize(a Matrix, kind SVDKind) (ok bool) { return ok } -// Kind returns the matrix.SVDKind of the decomposition. If no decomposition has been -// computed, Kind returns 0. +// 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.kind == 0 { - panic("svd: no decomposition computed") + if !svd.succFact() { + panic(badFact) } return svd.s[0] / svd.s[len(svd.s)-1] } @@ -129,8 +167,8 @@ func (svd *SVD) Cond() float64 { // // Values will panic if the receiver does not contain a successful factorization. func (svd *SVD) Values(s []float64) []float64 { - if svd.kind == 0 { - panic("svd: no decomposition computed") + if !svd.succFact() { + panic(badFact) } if s == nil { s = make([]float64, len(svd.s)) @@ -147,13 +185,16 @@ func (svd *SVD) Values(s []float64) []float64 { // 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 svd.Kind() == SVDFull, size m×min(m,n) if svd.Kind() == SVDThin, and -// UTo panics otherwise. If dst is nil, a new matrix of the appropriate size is -// allocated and returned. +// 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 != SVDFull && kind != SVDThin { - panic("mat: improper SVD kind") + if kind&SVDThinU == 0 && kind&SVDFullU == 0 { + panic("svd: u not computed during factorization") } r := svd.u.Rows c := svd.u.Cols @@ -178,13 +219,16 @@ func (svd *SVD) UTo(dst *Dense) *Dense { // 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 svd.Kind() == SVDFull, size n×min(m,n) if svd.Kind() == SVDThin, and -// VTo panics otherwise. If dst is nil, a new matrix of the appropriate size is -// allocated and returned. +// 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 != SVDFull && kind != SVDThin { - panic("mat: improper SVD kind") + if kind&SVDThinU == 0 && kind&SVDFullV == 0 { + panic("svd: v not computed during factorization") } r := svd.vt.Rows c := svd.vt.Cols