ruby - Ruby 矩阵计算中的浮点错误

Question

我正在编写一些涉及查找给定矩阵的特征向量的代码，令我惊讶的是 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 不能处理微小的矩阵，那我们怎么能相信它呢？还是我做错了什么？

Answer 1

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

Matrix 库改编自 JAMA, a Java matrix package ，它说它做了一个 numerical computation而不是 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

QR 算法:数值计算

查看 Matrix::EigenvalueDecomposition 的源代码，我发现它命名了 QR algorithm 的用法.我不完全理解数学的复杂性，但我想我可能理解为什么这个计算失败了。计算机制如下所述:

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.

在“伪”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.

在“伪”Ruby 中，继续:

eigenvectors = all_q_matrices.inject(:*).columns

数值计算中的浮点错误

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

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

一粒盐

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