Question

经典的 Gram-Schmidt（不稳定）

对于每列 j ，计算单一投影：

其中 与 的跨度正交的空间。

def cgs(A):
    m, n = A.shape
    Q = np.zeros([m,n], dtype=np.float64)
    R = np.zeros([n,n], dtype=np.float64)
    for j in range(n):
        v = A[:,j]
        for i in range(j):
            R[i,j] = np.dot(Q[:,i], A[:,j])
            v = v - (R[i,j] * Q[:,i])
        R[j,j] = np.linalg.norm(v)
        Q[:, j] = v / R[j,j]
    return Q, R

Q, R = cgs(A)

np.allclose(A, Q @ R)

# True

检查 Q 是酉矩阵。

np.allclose(np.eye(len(Q)), Q.dot(Q.T))

# True

np.allclose(npQ, -Q)

# True

R

'''
array([[ 0.02771,  0.02006, -0.0164 , ...,  0.00351,  0.00198,  0.00639],
       [ 0.     ,  0.10006, -0.00501, ...,  0.07689, -0.0379 , -0.03095],
       [ 0.     ,  0.     ,  0.01229, ...,  0.01635,  0.02988,  0.01442],
       ..., 
       [ 0.     ,  0.     ,  0.     , ...,  0.     , -0.     , -0.     ],
       [ 0.     ,  0.     ,  0.     , ...,  0.     ,  0.     , -0.     ],
       [ 0.     ,  0.     ,  0.     , ...,  0.     ,  0.     ,  0.     ]])
'''

Gram-Schmidt 应该让你想起来一点 Arnoldi 迭代（用于将矩阵转换为海森堡形式），因为它也是一个结构化的正交化。

改进版 Gram-Schmidt

经典（不稳定的）Gram-Schmidt：对于每列 j ，计算单一投影：

其中 与 的跨度正交的空间。

改进版 Gram-Schmidt：对于每列 j ，计算 n - 1 个投影：

import numpy as np
n = 3
A = np.random.rand(n,n).astype(np.float64)

def cgs(A):
    m, n = A.shape
    Q = np.zeros([m,n], dtype=np.float64)
    R = np.zeros([n,n], dtype=np.float64)
    for j in range(n):
        v = A[:,j]
        for i in range(j):
            R[i,j] = np.dot(Q[:,i], A[:,j])
            v = v - (R[i,j] * Q[:,i])
        R[j,j] = np.linalg.norm(v)
        Q[:, j] = v / R[j,j]
    return Q, R

def mgs(A):
    V = A.copy()
    m, n = A.shape
    Q = np.zeros([m,n], dtype=np.float64)
    R = np.zeros([n,n], dtype=np.float64)
    for i in range(n):
        R[i,i] = np.linalg.norm(V[:,i])
        Q[:,i] = V[:,i] / R[i,i]
        for j in range(i, n):
            R[i,j] = np.dot(Q[:,i],V[:,j])
            V[:,j] = V[:,j] - R[i,j]*Q[:,i]
    return Q, R

Q, R = mgs(A)

np.allclose(np.eye(len(Q)), Q.dot(Q.T.conj()))

# True

np.allclose(A, np.matmul(Q,R))

# True