stdlib/matrix/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