Ruby 矩阵计算中的浮点错误

floating point error in Ruby matrix calculation

我正在编写一些涉及查找给定矩阵的特征向量的代码,令我惊讶的是 Ruby 在简单情况下会产生一些不合理的结果。

例如,以下矩阵具有与特征值 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 很好地找到了特征值,但是特征向量爆炸了:

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

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

实际的特征向量应该是(7, 4, 4, 9).

这不是很麻烦吗?如果 Ruby 不能处理微小的矩阵,那我们怎么能相信它呢?还是我做错了什么?

不,这并不麻烦。该矩阵可能不适用于特定的特征向量算法实现。 Efficient and stable general eigenvector computation is nontrivial,毕竟

Matrix 库改编自 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

二维码算法:数值计算

正在查看 Matrix::EigenvalueDecomposition, I've found that it names the usage of the QR algorithm 的源代码。我不完全理解数学的复杂性,但我想我可能理解为什么 这个计算失败了。计算机制如下所述:

At the k-th step (starting with k = 0), we compute the QR decomposition Ak=QkRk ... Under certain conditions,[4] the matrices Ak 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.

在 "pseudo" Ruby 中,这在概念上意味着:

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

对于特征向量,它继续:

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.

在"pseudo"Ruby中,续:

eigenvectors = all_q_matrices.inject(:*).columns

数值计算中的浮点数错误

我们可以看到,进行了一次数值计算迭代来计算近似特征值,作为副作用,收集了一堆近似 Q 矩阵。然后,将这些近似的 Q 矩阵组合在一起形成特征向量。

近似值的复合可能是造成极其不准确结果的原因。 Math StackExchange 上的灾难性取消示例显示 simple quadratic computation with 400% relative error。您或许可以想象具有重复算术运算的迭代矩阵算法如何做得更糟。

一粒盐

同样,我对算法的数学原理和实现都没有深入的了解,所以我不知道确切地计算的哪些部分导致了您的特定问题85110032990182200% 错误,但我希望您现在可以了解它是如何发生的。