JezK
Edit File: eigenvalue_decomposition.rb
# frozen_string_literal: false class Matrix # Adapted from JAMA: http://math.nist.gov/javanumerics/jama/ # Eigenvalues and eigenvectors of a real matrix. # # Computes the eigenvalues and eigenvectors of a matrix A. # # If A is diagonalizable, this provides matrices V and D # such that A = V*D*V.inv, where D is the diagonal matrix with entries # equal to the eigenvalues and V is formed by the eigenvectors. # # If A is symmetric, then V is orthogonal and thus A = V*D*V.t class EigenvalueDecomposition # Constructs the eigenvalue decomposition for a square matrix +A+ # def initialize(a) # @d, @e: Arrays for internal storage of eigenvalues. # @v: Array for internal storage of eigenvectors. # @h: Array for internal storage of nonsymmetric Hessenberg form. raise TypeError, "Expected Matrix but got #{a.class}" unless a.is_a?(Matrix) @size = a.row_count @d = Array.new(@size, 0) @e = Array.new(@size, 0) if (@symmetric = a.symmetric?) @v = a.to_a tridiagonalize diagonalize else @v = Array.new(@size) { Array.new(@size, 0) } @h = a.to_a @ort = Array.new(@size, 0) reduce_to_hessenberg hessenberg_to_real_schur end end # Returns the eigenvector matrix +V+ # def eigenvector_matrix Matrix.send(:new, build_eigenvectors.transpose) end alias v eigenvector_matrix # Returns the inverse of the eigenvector matrix +V+ # def eigenvector_matrix_inv r = Matrix.send(:new, build_eigenvectors) r = r.transpose.inverse unless @symmetric r end alias v_inv eigenvector_matrix_inv # Returns the eigenvalues in an array # def eigenvalues values = @d.dup @e.each_with_index{|imag, i| values[i] = Complex(values[i], imag) unless imag == 0} values end # Returns an array of the eigenvectors # def eigenvectors build_eigenvectors.map{|ev| Vector.send(:new, ev)} end # Returns the block diagonal eigenvalue matrix +D+ # def eigenvalue_matrix Matrix.diagonal(*eigenvalues) end alias d eigenvalue_matrix # Returns [eigenvector_matrix, eigenvalue_matrix, eigenvector_matrix_inv] # def to_ary [v, d, v_inv] end alias_method :to_a, :to_ary private def build_eigenvectors # JAMA stores complex eigenvectors in a strange way # See http://web.archive.org/web/20111016032731/http://cio.nist.gov/esd/emaildir/lists/jama/msg01021.html @e.each_with_index.map do |imag, i| if imag == 0 Array.new(@size){|j| @v[j][i]} elsif imag > 0 Array.new(@size){|j| Complex(@v[j][i], @v[j][i+1])} else Array.new(@size){|j| Complex(@v[j][i-1], -@v[j][i])} end end end # Complex scalar division. def cdiv(xr, xi, yr, yi) if (yr.abs > yi.abs) r = yi/yr d = yr + r*yi [(xr + r*xi)/d, (xi - r*xr)/d] else r = yr/yi d = yi + r*yr [(r*xr + xi)/d, (r*xi - xr)/d] end end # Symmetric Householder reduction to tridiagonal form. def tridiagonalize # This is derived from the Algol procedures tred2 by # Bowdler, Martin, Reinsch, and Wilkinson, Handbook for # Auto. Comp., Vol.ii-Linear Algebra, and the corresponding # Fortran subroutine in EISPACK. @size.times do |j| @d[j] = @v[@size-1][j] end # Householder reduction to tridiagonal form. (@size-1).downto(0+1) do |i| # Scale to avoid under/overflow. scale = 0.0 h = 0.0 i.times do |k| scale = scale + @d[k].abs end if (scale == 0.0) @e[i] = @d[i-1] i.times do |j| @d[j] = @v[i-1][j] @v[i][j] = 0.0 @v[j][i] = 0.0 end else # Generate Householder vector. i.times do |k| @d[k] /= scale h += @d[k] * @d[k] end f = @d[i-1] g = Math.sqrt(h) if (f > 0) g = -g end @e[i] = scale * g h -= f * g @d[i-1] = f - g i.times do |j| @e[j] = 0.0 end # Apply similarity transformation to remaining columns. i.times do |j| f = @d[j] @v[j][i] = f g = @e[j] + @v[j][j] * f (j+1).upto(i-1) do |k| g += @v[k][j] * @d[k] @e[k] += @v[k][j] * f end @e[j] = g end f = 0.0 i.times do |j| @e[j] /= h f += @e[j] * @d[j] end hh = f / (h + h) i.times do |j| @e[j] -= hh * @d[j] end i.times do |j| f = @d[j] g = @e[j] j.upto(i-1) do |k| @v[k][j] -= (f * @e[k] + g * @d[k]) end @d[j] = @v[i-1][j] @v[i][j] = 0.0 end end @d[i] = h end # Accumulate transformations. 0.upto(@size-1-1) do |i| @v[@size-1][i] = @v[i][i] @v[i][i] = 1.0 h = @d[i+1] if (h != 0.0) 0.upto(i) do |k| @d[k] = @v[k][i+1] / h end 0.upto(i) do |j| g = 0.0 0.upto(i) do |k| g += @v[k][i+1] * @v[k][j] end 0.upto(i) do |k| @v[k][j] -= g * @d[k] end end end 0.upto(i) do |k| @v[k][i+1] = 0.0 end end @size.times do |j| @d[j] = @v[@size-1][j] @v[@size-1][j] = 0.0 end @v[@size-1][@size-1] = 1.0 @e[0] = 0.0 end # Symmetric tridiagonal QL algorithm. def diagonalize # This is derived from the Algol procedures tql2, by # Bowdler, Martin, Reinsch, and Wilkinson, Handbook for # Auto. Comp., Vol.ii-Linear Algebra, and the corresponding # Fortran subroutine in EISPACK. 1.upto(@size-1) do |i| @e[i-1] = @e[i] end @e[@size-1] = 0.0 f = 0.0 tst1 = 0.0 eps = Float::EPSILON @size.times do |l| # Find small subdiagonal element tst1 = [tst1, @d[l].abs + @e[l].abs].max m = l while (m < @size) do if (@e[m].abs <= eps*tst1) break end m+=1 end # If m == l, @d[l] is an eigenvalue, # otherwise, iterate. if (m > l) iter = 0 begin iter = iter + 1 # (Could check iteration count here.) # Compute implicit shift g = @d[l] p = (@d[l+1] - g) / (2.0 * @e[l]) r = Math.hypot(p, 1.0) if (p < 0) r = -r end @d[l] = @e[l] / (p + r) @d[l+1] = @e[l] * (p + r) dl1 = @d[l+1] h = g - @d[l] (l+2).upto(@size-1) do |i| @d[i] -= h end f += h # Implicit QL transformation. p = @d[m] c = 1.0 c2 = c c3 = c el1 = @e[l+1] s = 0.0 s2 = 0.0 (m-1).downto(l) do |i| c3 = c2 c2 = c s2 = s g = c * @e[i] h = c * p r = Math.hypot(p, @e[i]) @e[i+1] = s * r s = @e[i] / r c = p / r p = c * @d[i] - s * g @d[i+1] = h + s * (c * g + s * @d[i]) # Accumulate transformation. @size.times do |k| h = @v[k][i+1] @v[k][i+1] = s * @v[k][i] + c * h @v[k][i] = c * @v[k][i] - s * h end end p = -s * s2 * c3 * el1 * @e[l] / dl1 @e[l] = s * p @d[l] = c * p # Check for convergence. end while (@e[l].abs > eps*tst1) end @d[l] = @d[l] + f @e[l] = 0.0 end # Sort eigenvalues and corresponding vectors. 0.upto(@size-2) do |i| k = i p = @d[i] (i+1).upto(@size-1) do |j| if (@d[j] < p) k = j p = @d[j] end end if (k != i) @d[k] = @d[i] @d[i] = p @size.times do |j| p = @v[j][i] @v[j][i] = @v[j][k] @v[j][k] = p end end end end # Nonsymmetric reduction to Hessenberg form. def reduce_to_hessenberg # This is derived from the Algol procedures orthes and ortran, # by Martin and Wilkinson, Handbook for Auto. Comp., # Vol.ii-Linear Algebra, and the corresponding # Fortran subroutines in EISPACK. low = 0 high = @size-1 (low+1).upto(high-1) do |m| # Scale column. scale = 0.0 m.upto(high) do |i| scale = scale + @h[i][m-1].abs end if (scale != 0.0) # Compute Householder transformation. h = 0.0 high.downto(m) do |i| @ort[i] = @h[i][m-1]/scale h += @ort[i] * @ort[i] end g = Math.sqrt(h) if (@ort[m] > 0) g = -g end h -= @ort[m] * g @ort[m] = @ort[m] - g # Apply Householder similarity transformation # @h = (I-u*u'/h)*@h*(I-u*u')/h) m.upto(@size-1) do |j| f = 0.0 high.downto(m) do |i| f += @ort[i]*@h[i][j] end f = f/h m.upto(high) do |i| @h[i][j] -= f*@ort[i] end end 0.upto(high) do |i| f = 0.0 high.downto(m) do |j| f += @ort[j]*@h[i][j] end f = f/h m.upto(high) do |j| @h[i][j] -= f*@ort[j] end end @ort[m] = scale*@ort[m] @h[m][m-1] = scale*g end end # Accumulate transformations (Algol's ortran). @size.times do |i| @size.times do |j| @v[i][j] = (i == j ? 1.0 : 0.0) end end (high-1).downto(low+1) do |m| if (@h[m][m-1] != 0.0) (m+1).upto(high) do |i| @ort[i] = @h[i][m-1] end m.upto(high) do |j| g = 0.0 m.upto(high) do |i| g += @ort[i] * @v[i][j] end # Double division avoids possible underflow g = (g / @ort[m]) / @h[m][m-1] m.upto(high) do |i| @v[i][j] += g * @ort[i] end end end end end # Nonsymmetric reduction from Hessenberg to real Schur form. def hessenberg_to_real_schur # This is derived from the Algol procedure hqr2, # by Martin and Wilkinson, Handbook for Auto. Comp., # Vol.ii-Linear Algebra, and the corresponding # Fortran subroutine in EISPACK. # Initialize nn = @size n = nn-1 low = 0 high = nn-1 eps = Float::EPSILON exshift = 0.0 p = q = r = s = z = 0 # Store roots isolated by balanc and compute matrix norm norm = 0.0 nn.times do |i| if (i < low || i > high) @d[i] = @h[i][i] @e[i] = 0.0 end ([i-1, 0].max).upto(nn-1) do |j| norm = norm + @h[i][j].abs end end # Outer loop over eigenvalue index iter = 0 while (n >= low) do # Look for single small sub-diagonal element l = n while (l > low) do s = @h[l-1][l-1].abs + @h[l][l].abs if (s == 0.0) s = norm end if (@h[l][l-1].abs < eps * s) break end l-=1 end # Check for convergence # One root found if (l == n) @h[n][n] = @h[n][n] + exshift @d[n] = @h[n][n] @e[n] = 0.0 n-=1 iter = 0 # Two roots found elsif (l == n-1) w = @h[n][n-1] * @h[n-1][n] p = (@h[n-1][n-1] - @h[n][n]) / 2.0 q = p * p + w z = Math.sqrt(q.abs) @h[n][n] = @h[n][n] + exshift @h[n-1][n-1] = @h[n-1][n-1] + exshift x = @h[n][n] # Real pair if (q >= 0) if (p >= 0) z = p + z else z = p - z end @d[n-1] = x + z @d[n] = @d[n-1] if (z != 0.0) @d[n] = x - w / z end @e[n-1] = 0.0 @e[n] = 0.0 x = @h[n][n-1] s = x.abs + z.abs p = x / s q = z / s r = Math.sqrt(p * p+q * q) p /= r q /= r # Row modification (n-1).upto(nn-1) do |j| z = @h[n-1][j] @h[n-1][j] = q * z + p * @h[n][j] @h[n][j] = q * @h[n][j] - p * z end # Column modification 0.upto(n) do |i| z = @h[i][n-1] @h[i][n-1] = q * z + p * @h[i][n] @h[i][n] = q * @h[i][n] - p * z end # Accumulate transformations low.upto(high) do |i| z = @v[i][n-1] @v[i][n-1] = q * z + p * @v[i][n] @v[i][n] = q * @v[i][n] - p * z end # Complex pair else @d[n-1] = x + p @d[n] = x + p @e[n-1] = z @e[n] = -z end n -= 2 iter = 0 # No convergence yet else # Form shift x = @h[n][n] y = 0.0 w = 0.0 if (l < n) y = @h[n-1][n-1] w = @h[n][n-1] * @h[n-1][n] end # Wilkinson's original ad hoc shift if (iter == 10) exshift += x low.upto(n) do |i| @h[i][i] -= x end s = @h[n][n-1].abs + @h[n-1][n-2].abs x = y = 0.75 * s w = -0.4375 * s * s end # MATLAB's new ad hoc shift if (iter == 30) s = (y - x) / 2.0 s *= s + w if (s > 0) s = Math.sqrt(s) if (y < x) s = -s end s = x - w / ((y - x) / 2.0 + s) low.upto(n) do |i| @h[i][i] -= s end exshift += s x = y = w = 0.964 end end iter = iter + 1 # (Could check iteration count here.) # Look for two consecutive small sub-diagonal elements m = n-2 while (m >= l) do z = @h[m][m] r = x - z s = y - z p = (r * s - w) / @h[m+1][m] + @h[m][m+1] q = @h[m+1][m+1] - z - r - s r = @h[m+2][m+1] s = p.abs + q.abs + r.abs p /= s q /= s r /= s if (m == l) break end if (@h[m][m-1].abs * (q.abs + r.abs) < eps * (p.abs * (@h[m-1][m-1].abs + z.abs + @h[m+1][m+1].abs))) break end m-=1 end (m+2).upto(n) do |i| @h[i][i-2] = 0.0 if (i > m+2) @h[i][i-3] = 0.0 end end # Double QR step involving rows l:n and columns m:n m.upto(n-1) do |k| notlast = (k != n-1) if (k != m) p = @h[k][k-1] q = @h[k+1][k-1] r = (notlast ? @h[k+2][k-1] : 0.0) x = p.abs + q.abs + r.abs next if x == 0 p /= x q /= x r /= x end s = Math.sqrt(p * p + q * q + r * r) if (p < 0) s = -s end if (s != 0) if (k != m) @h[k][k-1] = -s * x elsif (l != m) @h[k][k-1] = -@h[k][k-1] end p += s x = p / s y = q / s z = r / s q /= p r /= p # Row modification k.upto(nn-1) do |j| p = @h[k][j] + q * @h[k+1][j] if (notlast) p += r * @h[k+2][j] @h[k+2][j] = @h[k+2][j] - p * z end @h[k][j] = @h[k][j] - p * x @h[k+1][j] = @h[k+1][j] - p * y end # Column modification 0.upto([n, k+3].min) do |i| p = x * @h[i][k] + y * @h[i][k+1] if (notlast) p += z * @h[i][k+2] @h[i][k+2] = @h[i][k+2] - p * r end @h[i][k] = @h[i][k] - p @h[i][k+1] = @h[i][k+1] - p * q end # Accumulate transformations low.upto(high) do |i| p = x * @v[i][k] + y * @v[i][k+1] if (notlast) p += z * @v[i][k+2] @v[i][k+2] = @v[i][k+2] - p * r end @v[i][k] = @v[i][k] - p @v[i][k+1] = @v[i][k+1] - p * q end end # (s != 0) end # k loop end # check convergence end # while (n >= low) # Backsubstitute to find vectors of upper triangular form if (norm == 0.0) return end (nn-1).downto(0) do |k| p = @d[k] q = @e[k] # Real vector if (q == 0) l = k @h[k][k] = 1.0 (k-1).downto(0) do |i| w = @h[i][i] - p r = 0.0 l.upto(k) do |j| r += @h[i][j] * @h[j][k] end if (@e[i] < 0.0) z = w s = r else l = i if (@e[i] == 0.0) if (w != 0.0) @h[i][k] = -r / w else @h[i][k] = -r / (eps * norm) end # Solve real equations else x = @h[i][i+1] y = @h[i+1][i] q = (@d[i] - p) * (@d[i] - p) + @e[i] * @e[i] t = (x * s - z * r) / q @h[i][k] = t if (x.abs > z.abs) @h[i+1][k] = (-r - w * t) / x else @h[i+1][k] = (-s - y * t) / z end end # Overflow control t = @h[i][k].abs if ((eps * t) * t > 1) i.upto(k) do |j| @h[j][k] = @h[j][k] / t end end end end # Complex vector elsif (q < 0) l = n-1 # Last vector component imaginary so matrix is triangular if (@h[n][n-1].abs > @h[n-1][n].abs) @h[n-1][n-1] = q / @h[n][n-1] @h[n-1][n] = -(@h[n][n] - p) / @h[n][n-1] else cdivr, cdivi = cdiv(0.0, -@h[n-1][n], @h[n-1][n-1]-p, q) @h[n-1][n-1] = cdivr @h[n-1][n] = cdivi end @h[n][n-1] = 0.0 @h[n][n] = 1.0 (n-2).downto(0) do |i| ra = 0.0 sa = 0.0 l.upto(n) do |j| ra = ra + @h[i][j] * @h[j][n-1] sa = sa + @h[i][j] * @h[j][n] end w = @h[i][i] - p if (@e[i] < 0.0) z = w r = ra s = sa else l = i if (@e[i] == 0) cdivr, cdivi = cdiv(-ra, -sa, w, q) @h[i][n-1] = cdivr @h[i][n] = cdivi else # Solve complex equations x = @h[i][i+1] y = @h[i+1][i] vr = (@d[i] - p) * (@d[i] - p) + @e[i] * @e[i] - q * q vi = (@d[i] - p) * 2.0 * q if (vr == 0.0 && vi == 0.0) vr = eps * norm * (w.abs + q.abs + x.abs + y.abs + z.abs) end cdivr, cdivi = cdiv(x*r-z*ra+q*sa, x*s-z*sa-q*ra, vr, vi) @h[i][n-1] = cdivr @h[i][n] = cdivi if (x.abs > (z.abs + q.abs)) @h[i+1][n-1] = (-ra - w * @h[i][n-1] + q * @h[i][n]) / x @h[i+1][n] = (-sa - w * @h[i][n] - q * @h[i][n-1]) / x else cdivr, cdivi = cdiv(-r-y*@h[i][n-1], -s-y*@h[i][n], z, q) @h[i+1][n-1] = cdivr @h[i+1][n] = cdivi end end # Overflow control t = [@h[i][n-1].abs, @h[i][n].abs].max if ((eps * t) * t > 1) i.upto(n) do |j| @h[j][n-1] = @h[j][n-1] / t @h[j][n] = @h[j][n] / t end end end end end end # Vectors of isolated roots nn.times do |i| if (i < low || i > high) i.upto(nn-1) do |j| @v[i][j] = @h[i][j] end end end # Back transformation to get eigenvectors of original matrix (nn-1).downto(low) do |j| low.upto(high) do |i| z = 0.0 low.upto([j, high].min) do |k| z += @v[i][k] * @h[k][j] end @v[i][j] = z end end end end end