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.
10 "gonum.org/v1/gonum/blas/blas64"
13 // HOGSVD is a type for creating and using the Higher Order Generalized Singular Value
14 // Decomposition (HOGSVD) of a set of matrices.
16 // The factorization is a linear transformation of the data sets from the given
17 // variable×sample spaces to reduced and diagonalized "eigenvariable"×"eigensample"
27 // Factorize computes the higher order generalized singular value decomposition (HOGSVD)
28 // of the n input r_i×c column tall matrices in m. HOGSV extends the GSVD case from 2 to n
31 // M_0 = U_0 * Σ_0 * V^T
32 // M_1 = U_1 * Σ_1 * V^T
36 // M_{n-1} = U_{n-1} * Σ_{n-1} * V^T
38 // where U_i are r_i×c matrices of singular vectors, Σ are c×c matrices singular values, and V
39 // is a c×c matrix of singular vectors.
41 // Factorize returns whether the decomposition succeeded. If the decomposition
42 // failed, routines that require a successful factorization will panic.
43 func (gsvd *HOGSVD) Factorize(m ...Matrix) (ok bool) {
44 // Factorize performs the HOGSVD factorisation
45 // essentially as described by Ponnapalli et al.
46 // https://doi.org/10.1371/journal.pone.0028072
49 panic("hogsvd: too few matrices")
54 a := make([]Cholesky, len(m))
69 ts.SymOuterK(1, d.T())
70 ok = a[i].Factorize(&ts)
72 gsvd.err = errors.New("hogsvd: cholesky decomposition failed")
77 s := getWorkspace(c, c, true)
79 sij := getWorkspace(c, c, false)
80 defer putWorkspace(sij)
81 for i, ai := range a {
82 for _, aj := range a[i+1:] {
83 gsvd.err = ai.SolveChol(sij, &aj)
89 gsvd.err = aj.SolveChol(sij, &ai)
96 s.Scale(1/float64(len(m)*(len(m)-1)), s)
99 ok = eig.Factorize(s.T(), false, true)
101 gsvd.err = errors.New("hogsvd: eigen decomposition failed")
106 for j := 0; j < c; j++ {
108 cv.ScaleVec(1/blas64.Nrm2(c, cv.mat), &cv)
111 b := make([]Dense, len(m))
112 biT := getWorkspace(c, r, false)
113 defer putWorkspace(biT)
114 for i, d := range m {
115 // All calls to reset will leave a zeroed
116 // matrix with capacity to store the result
117 // without additional allocation.
119 gsvd.err = biT.Solve(v, d.T())
132 // Err returns the reason for a factorization failure.
133 func (gsvd *HOGSVD) Err() error {
137 // Len returns the number of matrices that have been factorized. If Len returns
138 // zero, the factorization was not successful.
139 func (gsvd *HOGSVD) Len() int {
143 // UTo extracts the matrix U_n from the singular value decomposition, storing
144 // the result in-place into dst. U_n is size r×c.
145 // If dst is nil, a new matrix is allocated. The resulting U matrix is returned.
147 // UTo will panic if the receiver does not contain a successful factorization.
148 func (gsvd *HOGSVD) UTo(dst *Dense, n int) *Dense {
150 panic("hogsvd: unsuccessful factorization")
152 if n < 0 || gsvd.n <= n {
153 panic("hogsvd: invalid index")
157 r, c := gsvd.b[n].Dims()
158 dst = NewDense(r, c, nil)
160 dst.reuseAs(gsvd.b[n].Dims())
164 for j, f := range gsvd.Values(nil, n) {
171 // Values returns the nth set of singular values of the factorized system.
172 // If the input slice is non-nil, the values will be stored in-place into the slice.
173 // In this case, the slice must have length c, and Values will panic with
174 // matrix.ErrSliceLengthMismatch otherwise. If the input slice is nil,
175 // a new slice of the appropriate length will be allocated and returned.
177 // Values will panic if the receiver does not contain a successful factorization.
178 func (gsvd *HOGSVD) Values(s []float64, n int) []float64 {
180 panic("hogsvd: unsuccessful factorization")
182 if n < 0 || gsvd.n <= n {
183 panic("hogsvd: invalid index")
186 r, c := gsvd.b[n].Dims()
188 s = make([]float64, c)
189 } else if len(s) != c {
190 panic(ErrSliceLengthMismatch)
193 for j := 0; j < c; j++ {
194 v.ColViewOf(&gsvd.b[n], j)
195 s[j] = blas64.Nrm2(r, v.mat)
200 // VTo extracts the matrix V from the singular value decomposition, storing
201 // the result in-place into dst. V is size c×c.
202 // If dst is nil, a new matrix is allocated. The resulting V matrix is returned.
204 // VTo will panic if the receiver does not contain a successful factorization.
205 func (gsvd *HOGSVD) VTo(dst *Dense) *Dense {
207 panic("hogsvd: unsuccessful factorization")
210 r, c := gsvd.v.Dims()
211 dst = NewDense(r, c, nil)
213 dst.reuseAs(gsvd.v.Dims())