Question

以下示例来自 Trefethen 和 Bau 的第 9 讲，尽管从 MATLAB 翻译成 Python。

示例：经典与改进的 Gram-Schmidt

这个例子是 Trefethen 第 9 节的实验 2。 我们想要构造一个方阵 A ，它具有随机奇异向量和广泛变化的奇异值，间隔为 和 之间的 2 的倍数。

import matplotlib.pyplot as plt
from matplotlib import rcParams
%matplotlib inline

n = 100
U, X = np.linalg.qr(np.random.randn(n,n))   # 将 U 设为随机正交矩阵
V, X = np.linalg.qr(np.random.randn(n,n))   # 将 V 设为随机正交矩阵
S = np.diag(np.power(2,np.arange(-1,-(n+1),-1), dtype=float))  # 将 S 设为对角矩阵 w/ exp 
                                                              # 值在 2^-1 和 2^-(n+1) 之间

A = np.matmul(U,np.matmul(S,V))

QC, RC = cgs(A)
QM, RM = mgs(A)

plt.figure(figsize=(10,10))
plt.semilogy(np.diag(S), 'r.', basey=2, label="True Singular Values")
plt.semilogy(np.diag(RM), 'go', basey=2, label="Modified Gram-Shmidt")
plt.semilogy(np.diag(RC), 'bx', basey=2, label="Classic Gram-Shmidt")
plt.legend()
rcParams.update({'font.size': 18})

type(A[0,0]), type(RC[0,0]), type(S[0,0])

# (numpy.float64, numpy.float64, numpy.float64)

eps = np.finfo(np.float64).eps; eps

# 2.2204460492503131e-16

np.log2(eps), np.log2(np.sqrt(eps))

# (-52.0, -26.0)

示例：正交性的数值损失

这个例子是 Trefethen 第 9 节的实验 3。

A = np.array([[0.70000, 0.70711], [0.70001, 0.70711]])

A

'''
array([[ 0.7    ,  0.70711],
       [ 0.70001,  0.70711]])
'''

Gram-Schmidt：

Q1, R1 = mgs(A)

Householder：

R2, V, F = householder_lots(A)
Q2T = np.matmul(block_diag(np.eye(1), F[1]), F[0])

NumPy 的 Householder：

Q3, R3 = np.linalg.qr(A)

检查 QR 分解是否能用：

np.matmul(Q1, R1)

'''
array([[ 0.7   ,  0.7071],
       [ 0.7   ,  0.7071]])
'''

np.matmul(Q2T.T, R2)

'''
array([[ 0.7   ,  0.7071],
       [ 0.7   ,  0.7071]])
'''

np.matmul(Q3, R3)

'''
array([[ 0.7   ,  0.7071],
       [ 0.7   ,  0.7071]])
'''

检查 Q 多么接近完美正交。

np.linalg.norm(np.matmul(Q1.T, Q1) - np.eye(2))  # 改进的 Gram-Schmidt

# 3.2547268868202263e-11

np.linalg.norm(np.matmul(Q2T.T, Q2T) - np.eye(2))  # 我们的 Householder 实现

# 1.1110522984689321e-16

np.linalg.norm(np.matmul(Q3.T, Q3) - np.eye(2))  # Numpy（它使用 Householder）

# 2.5020189909116529e-16

GS（Q1）不如 Householder（Q2T，Q3）稳定。