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

   1 // Copyright ©2013 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/lapack"
   9         "gonum.org/v1/gonum/lapack/lapack64"
  10 )
  11
  12 const (
  13         badFact   = "mat: use without successful factorization"
  14         badNoVect = "mat: eigenvectors not computed"
  15 )
  16
  17 // EigenSym is a type for creating and manipulating the Eigen decomposition of
  18 // symmetric matrices.
  19 type EigenSym struct {
  20         vectorsComputed bool
  21
  22         values  []float64
  23         vectors *Dense
  24 }
  25
  26 // Factorize computes the eigenvalue decomposition of the symmetric matrix a.
  27 // The Eigen decomposition is defined as
  28 //  A = P * D * P^-1
  29 // where D is a diagonal matrix containing the eigenvalues of the matrix, and
  30 // P is a matrix of the eigenvectors of A. If the vectors input argument is
  31 // false, the eigenvectors are not computed.
  32 //
  33 // Factorize returns whether the decomposition succeeded. If the decomposition
  34 // failed, methods that require a successful factorization will panic.
  35 func (e *EigenSym) Factorize(a Symmetric, vectors bool) (ok bool) {
  36         n := a.Symmetric()
  37         sd := NewSymDense(n, nil)
  38         sd.CopySym(a)
  39
  40         jobz := lapack.EVJob(lapack.None)
  41         if vectors {
  42                 jobz = lapack.ComputeEV
  43         }
  44         w := make([]float64, n)
  45         work := []float64{0}
  46         lapack64.Syev(jobz, sd.mat, w, work, -1)
  47
  48         work = getFloats(int(work[0]), false)
  49         ok = lapack64.Syev(jobz, sd.mat, w, work, len(work))
  50         putFloats(work)
  51         if !ok {
  52                 e.vectorsComputed = false
  53                 e.values = nil
  54                 e.vectors = nil
  55                 return false
  56         }
  57         e.vectorsComputed = vectors
  58         e.values = w
  59         e.vectors = NewDense(n, n, sd.mat.Data)
  60         return true
  61 }
  62
  63 // succFact returns whether the receiver contains a successful factorization.
  64 func (e *EigenSym) succFact() bool {
  65         return len(e.values) != 0
  66 }
  67
  68 // Values extracts the eigenvalues of the factorized matrix. If dst is
  69 // non-nil, the values are stored in-place into dst. In this case
  70 // dst must have length n, otherwise Values will panic. If dst is
  71 // nil, then a new slice will be allocated of the proper length and filled
  72 // with the eigenvalues.
  73 //
  74 // Values panics if the Eigen decomposition was not successful.
  75 func (e *EigenSym) Values(dst []float64) []float64 {
  76         if !e.succFact() {
  77                 panic(badFact)
  78         }
  79         if dst == nil {
  80                 dst = make([]float64, len(e.values))
  81         }
  82         if len(dst) != len(e.values) {
  83                 panic(ErrSliceLengthMismatch)
  84         }
  85         copy(dst, e.values)
  86         return dst
  87 }
  88
  89 // EigenvectorsSym extracts the eigenvectors of the factorized matrix and stores
  90 // them in the receiver. Each eigenvector is a column corresponding to the
  91 // respective eigenvalue returned by e.Values.
  92 //
  93 // EigenvectorsSym panics if the factorization was not successful or if the
  94 // decomposition did not compute the eigenvectors.
  95 func (m *Dense) EigenvectorsSym(e *EigenSym) {
  96         if !e.succFact() {
  97                 panic(badFact)
  98         }
  99         if !e.vectorsComputed {
 100                 panic(badNoVect)
 101         }
 102         m.reuseAs(len(e.values), len(e.values))
 103         m.Copy(e.vectors)
 104 }
 105
 106 // Eigen is a type for creating and using the eigenvalue decomposition of a dense matrix.
 107 type Eigen struct {
 108         n int // The size of the factorized matrix.
 109
 110         right bool // have the right eigenvectors been computed
 111         left  bool // have the left eigenvectors been computed
 112
 113         values   []complex128
 114         rVectors *Dense
 115         lVectors *Dense
 116 }
 117
 118 // succFact returns whether the receiver contains a successful factorization.
 119 func (e *Eigen) succFact() bool {
 120         return len(e.values) != 0
 121 }
 122
 123 // Factorize computes the eigenvalues of the square matrix a, and optionally
 124 // the eigenvectors.
 125 //
 126 // A right eigenvalue/eigenvector combination is defined by
 127 //  A * x_r = λ * x_r
 128 // where x_r is the column vector called an eigenvector, and λ is the corresponding
 129 // eigenvector.
 130 //
 131 // Similarly, a left eigenvalue/eigenvector combination is defined by
 132 //  x_l * A = λ * x_l
 133 // The eigenvalues, but not the eigenvectors, are the same for both decompositions.
 134 //
 135 // Typically eigenvectors refer to right eigenvectors.
 136 //
 137 // In all cases, Eigen computes the eigenvalues of the matrix. If right and left
 138 // are true, then the right and left eigenvectors will be computed, respectively.
 139 // Eigen panics if the input matrix is not square.
 140 //
 141 // Factorize returns whether the decomposition succeeded. If the decomposition
 142 // failed, methods that require a successful factorization will panic.
 143 func (e *Eigen) Factorize(a Matrix, left, right bool) (ok bool) {
 144         // TODO(btracey): Change implementation to store VecDenses as a *CMat when
 145         // #308 is resolved.
 146
 147         // Copy a because it is modified during the Lapack call.
 148         r, c := a.Dims()
 149         if r != c {
 150                 panic(ErrShape)
 151         }
 152         var sd Dense
 153         sd.Clone(a)
 154
 155         var vl, vr Dense
 156         var jobvl lapack.LeftEVJob = lapack.None
 157         var jobvr lapack.RightEVJob = lapack.None
 158         if left {
 159                 vl = *NewDense(r, r, nil)
 160                 jobvl = lapack.ComputeLeftEV
 161         }
 162         if right {
 163                 vr = *NewDense(c, c, nil)
 164                 jobvr = lapack.ComputeRightEV
 165         }
 166
 167         wr := getFloats(c, false)
 168         defer putFloats(wr)
 169         wi := getFloats(c, false)
 170         defer putFloats(wi)
 171
 172         work := []float64{0}
 173         lapack64.Geev(jobvl, jobvr, sd.mat, wr, wi, vl.mat, vr.mat, work, -1)
 174         work = getFloats(int(work[0]), false)
 175         first := lapack64.Geev(jobvl, jobvr, sd.mat, wr, wi, vl.mat, vr.mat, work, len(work))
 176         putFloats(work)
 177
 178         if first != 0 {
 179                 e.values = nil
 180                 return false
 181         }
 182         e.n = r
 183         e.right = right
 184         e.left = left
 185         e.lVectors = &vl
 186         e.rVectors = &vr
 187         values := make([]complex128, r)
 188         for i, v := range wr {
 189                 values[i] = complex(v, wi[i])
 190         }
 191         e.values = values
 192         return true
 193 }
 194
 195 // Values extracts the eigenvalues of the factorized matrix. If dst is
 196 // non-nil, the values are stored in-place into dst. In this case
 197 // dst must have length n, otherwise Values will panic. If dst is
 198 // nil, then a new slice will be allocated of the proper length and
 199 // filed with the eigenvalues.
 200 //
 201 // Values panics if the Eigen decomposition was not successful.
 202 func (e *Eigen) Values(dst []complex128) []complex128 {
 203         if !e.succFact() {
 204                 panic(badFact)
 205         }
 206         if dst == nil {
 207                 dst = make([]complex128, e.n)
 208         }
 209         if len(dst) != e.n {
 210                 panic(ErrSliceLengthMismatch)
 211         }
 212         copy(dst, e.values)
 213         return dst
 214 }
 215
 216 // Vectors returns the right eigenvectors of the decomposition. Vectors
 217 // will panic if the right eigenvectors were not computed during the factorization,
 218 // or if the factorization was not successful.
 219 //
 220 // The returned matrix will contain the right eigenvectors of the decomposition
 221 // in the columns of the n×n matrix in the same order as their eigenvalues.
 222 // If the j-th eigenvalue is real, then
 223 //  u_j = VL[:,j],
 224 //  v_j = VR[:,j],
 225 // and if it is not real, then j and j+1 form a complex conjugate pair and the
 226 // eigenvectors can be recovered as
 227 //  u_j     = VL[:,j] + i*VL[:,j+1],
 228 //  u_{j+1} = VL[:,j] - i*VL[:,j+1],
 229 //  v_j     = VR[:,j] + i*VR[:,j+1],
 230 //  v_{j+1} = VR[:,j] - i*VR[:,j+1],
 231 // where i is the imaginary unit. The computed eigenvectors are normalized to
 232 // have Euclidean norm equal to 1 and largest component real.
 233 //
 234 // BUG: This signature and behavior will change when issue #308 is resolved.
 235 func (e *Eigen) Vectors() *Dense {
 236         if !e.succFact() {
 237                 panic(badFact)
 238         }
 239         if !e.right {
 240                 panic(badNoVect)
 241         }
 242         return DenseCopyOf(e.rVectors)
 243 }
 244
 245 // LeftVectors returns the left eigenvectors of the decomposition. LeftVectors
 246 // will panic if the left eigenvectors were not computed during the factorization.
 247 // or if the factorization was not successful.
 248 //
 249 // See the documentation in lapack64.Geev for the format of the vectors.
 250 //
 251 // BUG: This signature and behavior will change when issue #308 is resolved.
 252 func (e *Eigen) LeftVectors() *Dense {
 253         if !e.succFact() {
 254                 panic(badFact)
 255         }
 256         if !e.left {
 257                 panic(badNoVect)
 258         }
 259         return DenseCopyOf(e.lVectors)
 260 }