vendor/gonum.org/v1/gonum/mat/gsvd.go

   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.
   4
   5 package mat
   6
   7 import (
   8         "gonum.org/v1/gonum/blas/blas64"
   9         "gonum.org/v1/gonum/floats"
  10         "gonum.org/v1/gonum/lapack"
  11         "gonum.org/v1/gonum/lapack/lapack64"
  12 )
  13
  14 // GSVD is a type for creating and using the Generalized Singular Value Decomposition
  15 // (GSVD) of a matrix.
  16 //
  17 // The factorization is a linear transformation of the data sets from the given
  18 // variable×sample spaces to reduced and diagonalized "eigenvariable"×"eigensample"
  19 // spaces.
  20 type GSVD struct {
  21         kind GSVDKind
  22
  23         r, p, c, k, l int
  24         s1, s2        []float64
  25         a, b, u, v, q blas64.General
  26
  27         work  []float64
  28         iwork []int
  29 }
  30
  31 // Factorize computes the generalized singular value decomposition (GSVD) of the input
  32 // the r×c matrix A and the p×c matrix B. The singular values of A and B are computed
  33 // in all cases, while the singular vectors are optionally computed depending on the
  34 // input kind.
  35 //
  36 // The full singular value decomposition (kind == GSVDU|GSVDV|GSVDQ) deconstructs A and B as
  37 //  A = U * Σ₁ * [ 0 R ] * Q^T
  38 //
  39 //  B = V * Σ₂ * [ 0 R ] * Q^T
  40 // where Σ₁ and Σ₂ are r×(k+l) and p×(k+l) diagonal matrices of singular values, and
  41 // U, V and Q are r×r, p×p and c×c orthogonal matrices of singular vectors. k+l is the
  42 // effective numerical rank of the matrix [ A^T B^T ]^T.
  43 //
  44 // It is frequently not necessary to compute the full GSVD. Computation time and
  45 // storage costs can be reduced using the appropriate kind. Either only the singular
  46 // values can be computed (kind == SVDNone), or in conjunction with specific singular
  47 // vectors (kind bit set according to matrix.GSVDU, matrix.GSVDV and matrix.GSVDQ).
  48 //
  49 // Factorize returns whether the decomposition succeeded. If the decomposition
  50 // failed, routines that require a successful factorization will panic.
  51 func (gsvd *GSVD) Factorize(a, b Matrix, kind GSVDKind) (ok bool) {
  52         r, c := a.Dims()
  53         gsvd.r, gsvd.c = r, c
  54         p, c := b.Dims()
  55         gsvd.p = p
  56         if gsvd.c != c {
  57                 panic(ErrShape)
  58         }
  59         var jobU, jobV, jobQ lapack.GSVDJob
  60         switch {
  61         default:
  62                 panic("gsvd: bad input kind")
  63         case kind == GSVDNone:
  64                 jobU = lapack.GSVDNone
  65                 jobV = lapack.GSVDNone
  66                 jobQ = lapack.GSVDNone
  67         case (GSVDU|GSVDV|GSVDQ)&kind != 0:
  68                 if GSVDU&kind != 0 {
  69                         jobU = lapack.GSVDU
  70                         gsvd.u = blas64.General{
  71                                 Rows:   r,
  72                                 Cols:   r,
  73                                 Stride: r,
  74                                 Data:   use(gsvd.u.Data, r*r),
  75                         }
  76                 }
  77                 if GSVDV&kind != 0 {
  78                         jobV = lapack.GSVDV
  79                         gsvd.v = blas64.General{
  80                                 Rows:   p,
  81                                 Cols:   p,
  82                                 Stride: p,
  83                                 Data:   use(gsvd.v.Data, p*p),
  84                         }
  85                 }
  86                 if GSVDQ&kind != 0 {
  87                         jobQ = lapack.GSVDQ
  88                         gsvd.q = blas64.General{
  89                                 Rows:   c,
  90                                 Cols:   c,
  91                                 Stride: c,
  92                                 Data:   use(gsvd.q.Data, c*c),
  93                         }
  94                 }
  95         }
  96
  97         // A and B are destroyed on call, so copy the matrices.
  98         aCopy := DenseCopyOf(a)
  99         bCopy := DenseCopyOf(b)
 100
 101         gsvd.s1 = use(gsvd.s1, c)
 102         gsvd.s2 = use(gsvd.s2, c)
 103
 104         gsvd.iwork = useInt(gsvd.iwork, c)
 105
 106         gsvd.work = use(gsvd.work, 1)
 107         lapack64.Ggsvd3(jobU, jobV, jobQ, aCopy.mat, bCopy.mat, gsvd.s1, gsvd.s2, gsvd.u, gsvd.v, gsvd.q, gsvd.work, -1, gsvd.iwork)
 108         gsvd.work = use(gsvd.work, int(gsvd.work[0]))
 109         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)
 110         if ok {
 111                 gsvd.a = aCopy.mat
 112                 gsvd.b = bCopy.mat
 113                 gsvd.kind = kind
 114         }
 115         return ok
 116 }
 117
 118 // Kind returns the matrix.GSVDKind of the decomposition. If no decomposition has been
 119 // computed, Kind returns 0.
 120 func (gsvd *GSVD) Kind() GSVDKind {
 121         return gsvd.kind
 122 }
 123
 124 // Rank returns the k and l terms of the rank of [ A^T B^T ]^T.
 125 func (gsvd *GSVD) Rank() (k, l int) {
 126         return gsvd.k, gsvd.l
 127 }
 128
 129 // GeneralizedValues returns the generalized singular values of the factorized matrices.
 130 // If the input slice is non-nil, the values will be stored in-place into the slice.
 131 // In this case, the slice must have length min(r,c)-k, and GeneralizedValues will
 132 // panic with matrix.ErrSliceLengthMismatch otherwise. If the input slice is nil,
 133 // a new slice of the appropriate length will be allocated and returned.
 134 //
 135 // GeneralizedValues will panic if the receiver does not contain a successful factorization.
 136 func (gsvd *GSVD) GeneralizedValues(v []float64) []float64 {
 137         if gsvd.kind == 0 {
 138                 panic("gsvd: no decomposition computed")
 139         }
 140         r := gsvd.r
 141         c := gsvd.c
 142         k := gsvd.k
 143         d := min(r, c)
 144         if v == nil {
 145                 v = make([]float64, d-k)
 146         }
 147         if len(v) != d-k {
 148                 panic(ErrSliceLengthMismatch)
 149         }
 150         floats.DivTo(v, gsvd.s1[k:d], gsvd.s2[k:d])
 151         return v
 152 }
 153
 154 // ValuesA returns the singular values of the factorized A matrix.
 155 // If the input slice is non-nil, the values will be stored in-place into the slice.
 156 // In this case, the slice must have length min(r,c)-k, and ValuesA will panic with
 157 // matrix.ErrSliceLengthMismatch otherwise. If the input slice is nil,
 158 // a new slice of the appropriate length will be allocated and returned.
 159 //
 160 // ValuesA will panic if the receiver does not contain a successful factorization.
 161 func (gsvd *GSVD) ValuesA(s []float64) []float64 {
 162         if gsvd.kind == 0 {
 163                 panic("gsvd: no decomposition computed")
 164         }
 165         r := gsvd.r
 166         c := gsvd.c
 167         k := gsvd.k
 168         d := min(r, c)
 169         if s == nil {
 170                 s = make([]float64, d-k)
 171         }
 172         if len(s) != d-k {
 173                 panic(ErrSliceLengthMismatch)
 174         }
 175         copy(s, gsvd.s1[k:min(r, c)])
 176         return s
 177 }
 178
 179 // ValuesB returns the singular values of the factorized B matrix.
 180 // If the input slice is non-nil, the values will be stored in-place into the slice.
 181 // In this case, the slice must have length min(r,c)-k, and ValuesB will panic with
 182 // matrix.ErrSliceLengthMismatch otherwise. If the input slice is nil,
 183 // a new slice of the appropriate length will be allocated and returned.
 184 //
 185 // ValuesB will panic if the receiver does not contain a successful factorization.
 186 func (gsvd *GSVD) ValuesB(s []float64) []float64 {
 187         if gsvd.kind == 0 {
 188                 panic("gsvd: no decomposition computed")
 189         }
 190         r := gsvd.r
 191         c := gsvd.c
 192         k := gsvd.k
 193         d := min(r, c)
 194         if s == nil {
 195                 s = make([]float64, d-k)
 196         }
 197         if len(s) != d-k {
 198                 panic(ErrSliceLengthMismatch)
 199         }
 200         copy(s, gsvd.s2[k:d])
 201         return s
 202 }
 203
 204 // ZeroRTo extracts the matrix [ 0 R ] from the singular value decomposition, storing
 205 // the result in-place into dst. [ 0 R ] is size (k+l)×c.
 206 // If dst is nil, a new matrix is allocated. The resulting ZeroR matrix is returned.
 207 //
 208 // ZeroRTo will panic if the receiver does not contain a successful factorization.
 209 func (gsvd *GSVD) ZeroRTo(dst *Dense) *Dense {
 210         if gsvd.kind == 0 {
 211                 panic("gsvd: no decomposition computed")
 212         }
 213         r := gsvd.r
 214         c := gsvd.c
 215         k := gsvd.k
 216         l := gsvd.l
 217         h := min(k+l, r)
 218         if dst == nil {
 219                 dst = NewDense(k+l, c, nil)
 220         } else {
 221                 dst.reuseAsZeroed(k+l, c)
 222         }
 223         a := Dense{
 224                 mat:     gsvd.a,
 225                 capRows: r,
 226                 capCols: c,
 227         }
 228         dst.Slice(0, h, c-k-l, c).(*Dense).
 229                 Copy(a.Slice(0, h, c-k-l, c))
 230         if r < k+l {
 231                 b := Dense{
 232                         mat:     gsvd.b,
 233                         capRows: gsvd.p,
 234                         capCols: c,
 235                 }
 236                 dst.Slice(r, k+l, c+r-k-l, c).(*Dense).
 237                         Copy(b.Slice(r-k, l, c+r-k-l, c))
 238         }
 239         return dst
 240 }
 241
 242 // SigmaATo extracts the matrix Σ₁ from the singular value decomposition, storing
 243 // the result in-place into dst. Σ₁ is size r×(k+l).
 244 // If dst is nil, a new matrix is allocated. The resulting SigmaA matrix is returned.
 245 //
 246 // SigmaATo will panic if the receiver does not contain a successful factorization.
 247 func (gsvd *GSVD) SigmaATo(dst *Dense) *Dense {
 248         if gsvd.kind == 0 {
 249                 panic("gsvd: no decomposition computed")
 250         }
 251         r := gsvd.r
 252         k := gsvd.k
 253         l := gsvd.l
 254         if dst == nil {
 255                 dst = NewDense(r, k+l, nil)
 256         } else {
 257                 dst.reuseAsZeroed(r, k+l)
 258         }
 259         for i := 0; i < k; i++ {
 260                 dst.set(i, i, 1)
 261         }
 262         for i := k; i < min(r, k+l); i++ {
 263                 dst.set(i, i, gsvd.s1[i])
 264         }
 265         return dst
 266 }
 267
 268 // SigmaBTo extracts the matrix Σ₂ from the singular value decomposition, storing
 269 // the result in-place into dst. Σ₂ is size p×(k+l).
 270 // If dst is nil, a new matrix is allocated. The resulting SigmaB matrix is returned.
 271 //
 272 // SigmaBTo will panic if the receiver does not contain a successful factorization.
 273 func (gsvd *GSVD) SigmaBTo(dst *Dense) *Dense {
 274         if gsvd.kind == 0 {
 275                 panic("gsvd: no decomposition computed")
 276         }
 277         r := gsvd.r
 278         p := gsvd.p
 279         k := gsvd.k
 280         l := gsvd.l
 281         if dst == nil {
 282                 dst = NewDense(p, k+l, nil)
 283         } else {
 284                 dst.reuseAsZeroed(p, k+l)
 285         }
 286         for i := 0; i < min(l, r-k); i++ {
 287                 dst.set(i, i+k, gsvd.s2[k+i])
 288         }
 289         for i := r - k; i < l; i++ {
 290                 dst.set(i, i+k, 1)
 291         }
 292         return dst
 293 }
 294
 295 // UTo extracts the matrix U from the singular value decomposition, storing
 296 // the result in-place into dst. U is size r×r.
 297 // If dst is nil, a new matrix is allocated. The resulting U matrix is returned.
 298 //
 299 // UTo will panic if the receiver does not contain a successful factorization.
 300 func (gsvd *GSVD) UTo(dst *Dense) *Dense {
 301         if gsvd.kind&GSVDU == 0 {
 302                 panic("mat: improper GSVD kind")
 303         }
 304         r := gsvd.u.Rows
 305         c := gsvd.u.Cols
 306         if dst == nil {
 307                 dst = NewDense(r, c, nil)
 308         } else {
 309                 dst.reuseAs(r, c)
 310         }
 311
 312         tmp := &Dense{
 313                 mat:     gsvd.u,
 314                 capRows: r,
 315                 capCols: c,
 316         }
 317         dst.Copy(tmp)
 318         return dst
 319 }
 320
 321 // VTo extracts the matrix V from the singular value decomposition, storing
 322 // the result in-place into dst. V is size p×p.
 323 // If dst is nil, a new matrix is allocated. The resulting V matrix is returned.
 324 //
 325 // VTo will panic if the receiver does not contain a successful factorization.
 326 func (gsvd *GSVD) VTo(dst *Dense) *Dense {
 327         if gsvd.kind&GSVDV == 0 {
 328                 panic("mat: improper GSVD kind")
 329         }
 330         r := gsvd.v.Rows
 331         c := gsvd.v.Cols
 332         if dst == nil {
 333                 dst = NewDense(r, c, nil)
 334         } else {
 335                 dst.reuseAs(r, c)
 336         }
 337
 338         tmp := &Dense{
 339                 mat:     gsvd.v,
 340                 capRows: r,
 341                 capCols: c,
 342         }
 343         dst.Copy(tmp)
 344         return dst
 345 }
 346
 347 // QTo extracts the matrix Q from the singular value decomposition, storing
 348 // the result in-place into dst. Q is size c×c.
 349 // If dst is nil, a new matrix is allocated. The resulting Q matrix is returned.
 350 //
 351 // QTo will panic if the receiver does not contain a successful factorization.
 352 func (gsvd *GSVD) QTo(dst *Dense) *Dense {
 353         if gsvd.kind&GSVDQ == 0 {
 354                 panic("mat: improper GSVD kind")
 355         }
 356         r := gsvd.q.Rows
 357         c := gsvd.q.Cols
 358         if dst == nil {
 359                 dst = NewDense(r, c, nil)
 360         } else {
 361                 dst.reuseAs(r, c)
 362         }
 363
 364         tmp := &Dense{
 365                 mat:     gsvd.q,
 366                 capRows: r,
 367                 capCols: c,
 368         }
 369         dst.Copy(tmp)
 370         return dst
 371 }