Question

我们可以将 Arnoldi 迭代用于阶段 1（以及 QR 算法用于阶段 2）。

初始化

import numpy as np
n = 5
A0 = np.random.rand(n,n)  #.astype(np.float64)
A = A0 @ A0.T

np.set_printoptions(precision=5, suppress=True)

线性代数回顾：投影

当向量 b 投影到直线 a 上时，其投影 p 是 b 沿着直线 a 的一部分。

让我们看看 沉浸式线性代数在线版 的 第 3.2.2 节：投影 的交互图。

来源： 沉浸式数学

以下是将向量投影到平面上的样子：

来源： 最小二乘回归的线性代数视角

当向量 b 投影到直线 a 上时，其投影 p 是 b 沿着直线 a 的一部分。 所以 p 是 a 的一些倍数。 设 其中 是标量。

正交性

投影的关键是正交性：从 b 到 p 的直线（可以写成 ）垂直于 a 。

这意味着：

所以：

算法

动机

我们想要 Q 中的正交列和海森堡矩阵 H ，使得 AQ = QH 。

迭代式思考它：

其中 为 ， 为 。 这会创建一个可解决的递归关系。

来源：Trefethen 第 30 讲

Arnoldi 算法的伪代码

对于 Q 的第一列，以任意向量开始（标准化以便范数为 1）
for n=1,2,3...
    v = A @ nth col of Q
    for j=1,...n
        将 v 投影到 q_j，从 v 减去投影
        打算捕获 v 的一部分，它没有由 Q 的前一列跨越
        将系数储存在 H 中
    标准化 v, 之后使它为 Q 的第 (n+1) 列

请注意，我们将 A 乘以 Q 中的前一个向量，并删除与 Q 的现有列不正交的分量。

问题： A 的重复乘法会让你想起什么吗？

答案：幂方法也包括 A 的迭代乘法！

Arnoldi 迭代如何工作

通过 Arnoldi 迭代，我们找到了 Krylov 子空间的标准正交基。 Krylov 矩阵：

具有 QR 分解：

这与在 Arnoldi 迭代中发现的 Q 相同。 请注意，Arnoldi 迭代不会明确计算 K 或 R 。

直觉： K 包含关于 A 的最大特征值的良好信息，并且 QR 分解通过一次剥离一个近似特征向量，来揭示该信息。

Arnoldi 迭代是两件事：

• 许多数值线性代数迭代算法的基础
• 一种寻找非厄米矩阵特征值的技术（Trefethen，第 257 页）

Arnoldi 如何定位特征值

• 执行 Arnoldi 迭代
• 使用 QR 算法定期计算海森堡矩阵 H 的特征值（称为 Arnoldi 估计值或 Ritz 值）
• 检查这些值是否收敛。 如果是，它们可能是 A 的特征值。

实现

# 分解方阵 A @ Q ~= Q @ H
def arnoldi(A):
    m, n = A.shape
    assert(n <= m)

    # 海森堡矩阵
    H = np.zeros([n+1,n]) #, dtype=np.float64)
    # 正交列
    Q = np.zeros([m,n+1]) #, dtype=np.float64)
    # Q 的第 1 列是具有单位范数的随机列
    b = np.random.rand(m)
    Q[:,0] = b / np.linalg.norm(b)
    for j in range(n):
        v = A @ Q[:,j]
        for i in range(j+1):
            # 这来自 v 投影到 q 的公式。
            # 由于列 q 是正交的，q dot q = 1
            H[i,j] = np.dot(Q[:,i], v)
            v = v - (H[i,j] * Q[:,i])
        H[j+1,j] = np.linalg.norm(v)
        Q[:,j+1] = v / H[j+1,j]

        # 打印这个看到收敛，在实践中使用会很慢
        print(np.linalg.norm(A @ Q[:,:-1] - Q @ H))
    return Q[:,:-1], H[:-1,:]

Q, H = arnoldi(A)

'''
8.59112969133
4.45398729097
0.935693639985
3.36613943339
0.817740180293
'''

检查 H 是三对角的。

H

'''
array([[ 5.62746,  4.05085, -0.     ,  0.     , -0.     ],
       [ 4.05085,  3.07109,  0.33036,  0.     , -0.     ],
       [ 0.     ,  0.33036,  0.98297,  0.11172, -0.     ],
       [ 0.     ,  0.     ,  0.11172,  0.29777,  0.07923],
       [ 0.     ,  0.     ,  0.     ,  0.07923,  0.06034]])
'''

练习

编写代码来验证：

• AQ = QH
• Q 是正交的

答案

# 练习：

# True

# 练习：

# True

一般案例：

一般矩阵：现在我们可以在我们的一般矩阵 A （非对称）上执行此操作。 在这种情况下，我们得到的是海森堡而不三对角。

Q0, H0 = arnoldi(A0)

'''
1.44287067354
1.06234006889
0.689291414367
0.918098818651
4.7124490411e-16
'''

检查 H 是海森堡的。

H0

'''
array([[ 1.67416,  0.83233, -0.39284,  0.10833,  0.63853],
       [ 1.64571,  1.16678, -0.54779,  0.50529,  0.28515],
       [ 0.     ,  0.16654, -0.22314,  0.08577, -0.02334],
       [ 0.     ,  0.     ,  0.79017,  0.11732,  0.58978],
       [ 0.     ,  0.     ,  0.     ,  0.43238, -0.07413]])
'''

np.allclose(A0 @ Q0, Q0 @ H0)

# True

np.allclose(np.eye(len(Q0)), Q0.T @ Q0), np.allclose(np.eye(len(Q0)), Q0 @ Q0.T)

# (True, True)