floating point error in Ruby matrix calculation

Question

I'm writing some code that involves finding the eigenvectors of a given matrix, and was surprised that Ruby produces some unreasonable results in simple cases.

For example, the following matrix has an eigenvector associated with eigenvalue 1:

> m = Matrix[[0r, 1/2r, 1/2r, 1/3r],
             [0r,  0r,  1/4r, 1/3r],
             [0r, 1/4r,  0r,  1/3r],
             [1r, 1/4r, 1/4r,  0r]]

Ruby finds the eigenvalues well enough, but the eigenvector explodes:

> m.eigen.eigenvalues[2]
=> 1.0000000000000009

m.eigen.eigenvectors[2]
=> Vector[5.957702309312754e+15, 5.957702309312748e+15, 5.957702309312743e+15, 5.957702309312753e+15]

The actual eigenvector should be (7, 4, 4, 9).

Isn't this troubling? If Ruby can't handle tiny matrices, then how can we trust it at all? Or am I doing something wrong?

Consecrate

The result looks wrong to me too (as we can also see by doing a m*m.eigen.eigenvectors[2]). It would be helpful if you could submit a bug report using your example. Maybe you can find even a small matrix which yields incorrect results? – Overtire

Could this have anything to do with the eigenvalues being complex, rather than real? Your matrix is not symmetric, so is certainly not guaranteed to have real eigenvalues. Doing the eigen-decomposition with numpy.linalg.eig() in Python, one finds that two of the eigenvalues are complex. If R's eigen method is geared towards symmetric matrices, then it would be no surprise if it struggled to correctly find any of the eigenvectors correctly. – Trumpet

It could be the eigenvectors are not normalized, as done in Numpy. Just for the record, I tried your example with Ruby stdlib and Python's Numpy. Eigenvalues match (different precision though), but the vectors don't, even normalized (I just checked the 3rd you mention). – Xeres

@Trumpet I should have mentioned that I was constructing the matrices especially to have an eigenvalue of $1$, corresponding to an eigenvector whose entries are real and positive. This is guaranteed by the Perron-Frobenius theorem. True, the other eigenvalues are complex, but the eigenvalue of $1$ seems more important, as it has the biggest absolute value. – Polluted

@EricPlaton But the vector above is a multiple of $[1,1,1,1]$, not $[7,4,4,9]$, so it isn't just a normalization problem. – Polluted

Yes. My comment is convoluted. Normalization addresses the “explosion”. The problem is the mismatch with Numpy. – Xeres

@EricPlaton Ah, I see. I recall trying this on a different computer several months ago, and it gave the correct result. I thought at the time that the problem was fixed (e.g., from using a more recent version of Ruby), but I tried installing 2.5.1 just now, and I get the same problem. Hmm. – Polluted

My tests were on Ruby 2.3.1. – Xeres

Answer 1

No, this is not troubling. That matrix likely just doesn't work well with that particular eigenvector algorithm implementation. Efficient and stable general eigenvector computation is nontrivial, after all.

The Matrix library is adapted from JAMA, a Java matrix package, which says it does a numerical computation and not a symbolic computation:

Not Covered. JAMA is by no means a complete linear algebra environment ... it focuses on the principle mathematical functionality required to do numerical linear algebra

The QR Algorithm: Numerical Computation

Looking at the source code for Matrix::EigenvalueDecomposition, I've found that it names the usage of the QR algorithm. I don't fully understand the intricacies of the mathematics, but I think I might understand why this computation fails. The mechanism of computation works as stated:

At the k-th step (starting with k = 0), we compute the QR decomposition A_k=Q_kR_k ... Under certain conditions,^[4] the matrices A_k converge to a triangular matrix, the Schur form of A. The eigenvalues of a triangular matrix are listed on the diagonal, and the eigenvalue problem is solved.

In "pseudo" Ruby, this conceptually means:

working_matrix = orig_matrix.dup
all_q_matrices = []
loop do
  q, r = working_matrix.qr_decomposition
  all_q_matrices << q
  next_matrix = r * q
  break if difference_between(working_matrix, next_matrix) < accuracy_threshold
end
eigenvalues = working_matrix.diagonal_values

For eigenvectors, it continues:

upon convergence, AQ = QΛ, where Λ is the diagonal matrix of eigenvalues to which A converged, and where Q is a composite of all the orthogonal similarity transforms required to get there. Thus the columns of Q are the eigenvectors.

In "pseudo" Ruby, continued:

eigenvectors = all_q_matrices.inject(:*).columns

Floating Point Error in Numerical Computations

We can see that an iteration of numerical computations are made to compute the approximate eigenvalues, and as a side-effect, a bunch of approximate Q matrices are collected. Then, these approximated Q matrices are composed together to form the eigenvectors.

The compounding of approximations is what probably caused the wildly inaccurate results. An example of catastrophic cancellation on Math StackExchange shows a simple quadratic computation with 400% relative error. You might be able to imagine how an iterative matrix algorithm with repeated arithmetic operations could do much worse.

A Grain of Salt

Again, I don't have a deep understanding of the mathematics of the algorithm nor the implementation, so I don't know precisely what parts of the computation caused your specific 85110032990182200% error, but I hope you now can understand how it may have happened.

Crept

I'm still troubled, but thank you for the detailed explanation. :) – Polluted