文章来源于网络收集而来，版权归原创者所有，如有侵权请及时联系！

Householder

发布于 2025-01-01 12:38:42 字数 5648 浏览 0 评论 0 收藏 0

引言

$\begin{array}{ l | l | c } \hline Gram-Schmidt & Triangular\, Orthogonalization & A R_1 R_2 \cdots R_n = Q \\ Householder & Orthogonal\, Triangularization & Q_n \cdots Q_2 Q_1 A = R \\ \hline \end{array}$

Householder 反射产生更接近正交的矩阵 Q ，具有舍入误差

Gram-Schmidt 可以部分停止，留下 A 的前 n 列的简化 QR。

初始化

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

Q = np.zeros([n,n], dtype=np.float64)
R = np.zeros([n,n], dtype=np.float64)

A

'''
array([[ 0.5435,  0.6379,  0.4011,  0.5773],
       [ 0.0054,  0.8049,  0.6804,  0.0821],
       [ 0.2832,  0.2416,  0.8656,  0.8099],
       [ 0.1139,  0.9621,  0.7623,  0.5648]])
'''

from scipy.linalg import block_diag

np.set_printoptions(5)

算法

我添加了更多的计算和更多的信息，因为它说明了算法的工作原理。此版本也返回 Householder 反射。

def householder_lots(A):
    m, n = A.shape
    R = np.copy(A)
    V = []
    Fs = []
    for k in range(n):
        v = np.copy(R[k:,k])
        v = np.reshape(v, (n-k, 1))
        v[0] += np.sign(v[0]) * np.linalg.norm(v)
        v /= np.linalg.norm(v)
        R[k:,k:] = R[k:,k:] - 2*np.matmul(v, np.matmul(v.T, R[k:,k:]))
        V.append(v)
        F = np.eye(n-k) - 2 * np.matmul(v, v.T)/np.matmul(v.T, v)
        Fs.append(F)
    return R, V, Fs

检查 R 是上三角。

R

'''
array([[-0.62337, -0.84873, -0.88817, -0.97516],
       [ 0.     , -1.14818, -0.86417, -0.30109],
       [ 0.     ,  0.     , -0.64691, -0.45234],
       [-0.     ,  0.     ,  0.     , -0.26191]])
'''

作为检查，我们将使用分块矩阵 F 计算 $Q^T$ 和 R 。矩阵 F 是 householder 反射。

请注意，这不是一种处理 Q 的有效计算方式。在大多数情况下，你实际上并不需要 Q 。例如，如果你使用 QR 来求解最小二乘，则只需要 Q * b 。

对于隐式计算乘积 Q * b 或 Qx 的技巧，请参阅 Trefethen 第 74 页。
请参阅这些讲义，了解 Householder 的不同实现，它同时计算 Q ，作为 R 的一部分。

QT = np.matmul(block_diag(np.eye(3), F[3]), 
               np.matmul(block_diag(np.eye(2), F[2]), 
                         np.matmul(block_diag(np.eye(1), F[1]), F[0])))

F[1]

'''
array([[-0.69502,  0.10379, -0.71146],
       [ 0.10379,  0.99364,  0.04356],
       [-0.71146,  0.04356,  0.70138]])
'''

block_diag(np.eye(1), F[1])

'''
array([[ 1.     ,  0.     ,  0.     ,  0.     ],
       [ 0.     , -0.69502,  0.10379, -0.71146],
       [ 0.     ,  0.10379,  0.99364,  0.04356],
       [ 0.     , -0.71146,  0.04356,  0.70138]])
'''

block_diag(np.eye(2), F[2])

'''
array([[ 1.     ,  0.     ,  0.     ,  0.     ],
       [ 0.     ,  1.     ,  0.     ,  0.     ],
       [ 0.     ,  0.     , -0.99989,  0.01452],
       [ 0.     ,  0.     ,  0.01452,  0.99989]])
'''

block_diag(np.eye(3), F[3])

'''
array([[ 1.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.],
       [ 0.,  0.,  1.,  0.],
       [ 0.,  0.,  0., -1.]])
'''

np.matmul(block_diag(np.eye(1), F[1]), F[0])

'''
array([[-0.87185, -0.00861, -0.45431, -0.18279],
       [ 0.08888, -0.69462,  0.12536, -0.70278],
       [-0.46028,  0.10167,  0.88193, -0.00138],
       [-0.14187, -0.71211,  0.00913,  0.68753]])
'''

QT

'''
array([[-0.87185, -0.00861, -0.45431, -0.18279],
       [ 0.08888, -0.69462,  0.12536, -0.70278],
       [ 0.45817, -0.112  , -0.88171,  0.01136],
       [ 0.14854,  0.71056, -0.02193, -0.68743]])
'''

R2 = np.matmul(block_diag(np.eye(3), F[3]), 
               np.matmul(block_diag(np.eye(2), F[2]),
                         np.matmul(block_diag(np.eye(1), F[1]), 
                                   np.matmul(F[0], A))))

np.allclose(A, np.matmul(np.transpose(QT), R2))

# True

np.allclose(R, R2)

# True

这是 Householder 的简洁版本（尽管我创建了一个新的 R ，而不是覆盖 A ，并原地计算它）。

def householder(A):
    m, n = A.shape
    R = np.copy(A)
    Q = np.eye(m)
    V = []
    for k in range(n):
        v = np.copy(R[k:,k])
        v = np.reshape(v, (n-k, 1))
        v[0] += np.sign(v[0]) * np.linalg.norm(v)
        v /= np.linalg.norm(v)
        R[k:,k:] = R[k:,k:] - 2 * v @ v.T @ R[k:,k:]
        V.append(v)
    return R, V

RH, VH = householder(A)

检查 R 是对角的。

RH

'''
array([[-0.62337, -0.84873, -0.88817, -0.97516],
       [-0.     , -1.14818, -0.86417, -0.30109],
       [-0.     , -0.     , -0.64691, -0.45234],
       [-0.     ,  0.     ,  0.     , -0.26191]])
'''

VH

'''
[array([[ 0.96743],
        [ 0.00445],
        [ 0.2348 ],
        [ 0.09447]]), array([[ 0.9206 ],
        [-0.05637],
        [ 0.38641]]), array([[ 0.99997],
        [-0.00726]]), array([[ 1.]])]
'''

np.allclose(R, RH)

# True

def implicit_Qx(V,x):
    n = len(x)
    for k in range(n-1,-1,-1):
        x[k:n] -= 2*np.matmul(v[-k], np.matmul(v[-k], x[k:n]))

A

'''
array([[ 0.54348,  0.63791,  0.40114,  0.57728],
       [ 0.00537,  0.80485,  0.68037,  0.0821 ],
       [ 0.2832 ,  0.24164,  0.86556,  0.80986],
       [ 0.11395,  0.96205,  0.76232,  0.56475]])
'''

经典和改良的 Gram-Schmidt 都需要 2mn^2 个浮点运算。