summaryrefslogtreecommitdiff log msg author committer range
diff options
 context: 12345678910152025303540 space: includeignore mode: unifiedssdiffstat only
Diffstat (limited to 'lib/matrix/eigenvalue_decomposition.rb')
-rw-r--r--lib/matrix/eigenvalue_decomposition.rb886
1 files changed, 886 insertions, 0 deletions
 diff --git a/lib/matrix/eigenvalue_decomposition.rb b/lib/matrix/eigenvalue_decomposition.rbnew file mode 100644index 0000000..55e41aa--- /dev/null+++ b/lib/matrix/eigenvalue_decomposition.rb@@ -0,0 +1,886 @@+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_size+ @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://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], -@v[j][i-1])}+ 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+ if (x != 0.0)+ p /= x+ q /= x+ r /= x+ end+ end+ if (x == 0.0)+ break+ 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 |n|+ p = @d[n]+ q = @e[n]++ # Real vector++ if (q == 0)+ l = n+ @h[n][n] = 1.0+ (n-1).downto(0) do |i|+ w = @h[i][i] - p+ r = 0.0+ l.upto(n) do |j|+ r += @h[i][j] * @h[j][n]+ end+ if (@e[i] < 0.0)+ z = w+ s = r+ else+ l = i+ if (@e[i] == 0.0)+ if (w != 0.0)+ @h[i][n] = -r / w+ else+ @h[i][n] = -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][n] = t+ if (x.abs > z.abs)+ @h[i+1][n] = (-r - w * t) / x+ else+ @h[i+1][n] = (-s - y * t) / z+ end+ end++ # Overflow control++ t = @h[i][n].abs+ if ((eps * t) * t > 1)+ i.upto(n) do |j|+ @h[j][n] = @h[j][n] / 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