如何在 Python 中实现基于 Householder 的 QR 分解？
How can you implement Householder based QR decomposition in Python?
我目前正在尝试实现矩形矩阵的基于 Householder 的 QR 分解，如 http://eprints.ma.man.ac.uk/1192/1/qrupdating_12nov08.pdf（第 3、4、5 页）中所述。
显然我弄错了一些伪代码，因为 (1) 我的结果不同于 numpy.qr.linalg() 和 (2) 我的例程生成的矩阵 R 不是上三角矩阵。
我的代码（也可在 https://pyfiddle.io/fiddle/afcc2e0e-0857-4cb2-adb5-06ff9b80c9d3/?i=true 下获得）
import math
import argparse
import numpy as np
from typing import Union
def householder(alpha: float, x: np.ndarray) -> Union[np.ndarray, int]:
"""
Computes Householder vector for alpha and x.
:param alpha:
:param x:
:return:
"""
s = math.pow(np.linalg.norm(x, ord=2), 2)
v = x
if s == 0:
tau = 0
else:
t = math.sqrt(alpha * alpha + s)
v_one = alpha - t if alpha <= 0 else -s / (alpha + t)
tau = 2 * v_one * v_one / (s + v_one * v_one)
v /= v_one
return v, tau
def qr_decomposition(A: np.ndarray, m: int, n: int) -> Union[np.ndarray, np.ndarray]:
"""
Applies Householder-based QR decomposition on specified matrix A.
:param A:
:param m:
:param n:
:return:
"""
H = []
R = A
Q = A
I = np.eye(m, m)
for j in range(0, n):
# Apply Householder transformation.
x = A[j + 1:m, j]
v_householder, tau = householder(np.linalg.norm(x), x)
v = np.zeros((1, m))
v[0, j] = 1
v[0, j + 1:m] = v_householder
res = I - tau * v * np.transpose(v)
R = np.matmul(res, R)
H.append(res)
return Q, R
m = 10
n = 8
A = np.random.rand(m, n)
q, r = np.linalg.qr(A)
Q, R = qr_decomposition(A, m, n)
print("*****")
print(Q)
print(q)
print("-----")
print(R)
print(r)
所以我不清楚如何将零引入我的 R matrix/about 我的代码的哪一部分不正确。我会很高兴任何指示！
非常感谢您的宝贵时间。
您链接到的注释中有一堆 problems/missing 详细信息。在咨询 few other sources (including this very useful textbook) 之后，我能够想出一个类似的工作实现。
工作算法
这是 qr_decomposition 工作版本的代码：
import numpy as np
from typing import Union
def householder(x: np.ndarray) -> Union[np.ndarray, int]:
alpha = x[0]
s = np.power(np.linalg.norm(x[1:]), 2)
v = x.copy()
if s == 0:
tau = 0
else:
t = np.sqrt(alpha**2 + s)
v[0] = alpha - t if alpha <= 0 else -s / (alpha + t)
tau = 2 * v[0]**2 / (s + v[0]**2)
v /= v[0]
return v, tau
def qr_decomposition(A: np.ndarray) -> Union[np.ndarray, np.ndarray]:
m,n = A.shape
R = A.copy()
Q = np.identity(m)
for j in range(0, n):
# Apply Householder transformation.
v, tau = householder(R[j:, j])
H = np.identity(m)
H[j:, j:] -= tau * v.reshape(-1, 1) @ v
R = H @ R
Q = H @ Q
return Q[:n].T, R[:n]
m = 5
n = 4
A = np.random.rand(m, n)
q, r = np.linalg.qr(A)
Q, R = qr_decomposition(A)
with np.printoptions(linewidth=9999, precision=20, suppress=True):
print("**** Q from qr_decomposition")
print(Q)
print("**** Q from np.linalg.qr")
print(q)
print()
print("**** R from qr_decomposition")
print(R)
print("**** R from np.linalg.qr")
print(r)
输出：
**** Q from qr_decomposition
[[ 0.5194188817843675 -0.10699353671401633 0.4322294754656072 -0.7293293270703678 ]
[ 0.5218635773595086 0.11737804362574514 -0.5171653705211056 0.04467925806590414]
[ 0.34858177783013133 0.6023104248793858 -0.33329256746256875 -0.03450824948274838]
[ 0.03371048915852807 0.6655221685383623 0.6127023580593225 0.28795294754791 ]
[ 0.5789790833500734 -0.411189947884951 0.24337120818874305 0.618041080584351 ]]
**** Q from np.linalg.qr
[[-0.5194188817843672 0.10699353671401617 0.4322294754656068 0.7293293270703679 ]
[-0.5218635773595086 -0.11737804362574503 -0.5171653705211053 -0.044679258065904115]
[-0.3485817778301313 -0.6023104248793857 -0.33329256746256863 0.03450824948274819 ]
[-0.03371048915852807 -0.665522168538362 0.6127023580593226 -0.2879529475479097 ]
[-0.5789790833500733 0.41118994788495106 0.24337120818874317 -0.6180410805843508 ]]
**** R from qr_decomposition
[[ 0.6894219296137802 1.042676051151294 1.3418719684631446 1.2498925815126485 ]
[ 0.00000000000000000685 0.7076056836914905 0.29883043386651403 0.41955370595004277 ]
[-0.0000000000000000097 -0.00000000000000007292 0.5304551654027297 0.18966088433421135 ]
[-0.00000000000000000662 0.00000000000000008718 0.00000000000000002322 0.6156558913022807 ]]
**** R from np.linalg.qr
[[-0.6894219296137803 -1.042676051151294 -1.3418719684631442 -1.2498925815126483 ]
[ 0. -0.7076056836914905 -0.29883043386651376 -0.4195537059500425 ]
[ 0. 0. 0.53045516540273 0.18966088433421188]
[ 0. 0. 0. -0.6156558913022805 ]]
此版本的 qr_decomposition 几乎完全再现了 np.linalg.qr 的输出。差异在下面进行了评论。
输出的数值精度
np.linalg.qr 和 qr_decomposition 的输出值与高精度匹配。但是，qr_decomposition 用于生成 R 中的零的计算组合并没有完全抵消，因此零实际上并不完全等于零。
事实证明，np.linalg.qr 并没有使用任何花哨的浮点技巧来确保其输出中的零为 0.0。它 just calls np.triu，强制将这些值设置为 0.0。所以要达到相同的结果，只需将 qr_decomposition 中的 return 行更改为：
return Q[:n].T, np.triu(R[:n])
输出中的符号 (+/-)
Q 和 R 中的一些 +/- 符号在 np.linalg.qr 和 qr_decomposition 的输出中是不同的，但这不是真正的问题，因为有很多有效的选择标志（见此discussion of the uniqueness of Q and R）。您可以通过使用替代算法生成 v 和 tau:
来完全匹配 np.linalg.qr 的符号约定
def householder_vectorized(a):
"""Use this version of householder to reproduce the output of np.linalg.qr
exactly (specifically, to match the sign convention it uses)
based on https://rosettacode.org/wiki/QR_decomposition#Python
"""
v = a / (a[0] + np.copysign(np.linalg.norm(a), a[0]))
v[0] = 1
tau = 2 / (v.T @ v)
return v,tau
与np.linalg.qr
的输出完全匹配
将它们放在一起，这个版本的 qr_decomposition 将与 np.linalg.qr 的输出完全匹配：
import numpy as np
from typing import Union
def qr_decomposition(A: np.ndarray) -> Union[np.ndarray, np.ndarray]:
m,n = A.shape
R = A.copy()
Q = np.identity(m)
for j in range(0, n):
# Apply Householder transformation.
v, tau = householder_vectorized(R[j:, j, np.newaxis])
H = np.identity(m)
H[j:, j:] -= tau * (v @ v.T)
R = H @ R
Q = H @ Q
return Q[:n].T, np.triu(R[:n])
m = 5
n = 4
A = np.random.rand(m, n)
q, r = np.linalg.qr(A)
Q, R = qr_decomposition(A)
with np.printoptions(linewidth=9999, precision=20, suppress=True):
print("**** Q from qr_decomposition")
print(Q)
print("**** Q from np.linalg.qr")
print(q)
print()
print("**** R from qr_decomposition")
print(R)
print("**** R from np.linalg.qr")
print(r)
输出：
**** Q from qr_decomposition
[[-0.10345123000824041 0.6455437884382418 0.44810714367794663 -0.03963544711256745 ]
[-0.55856415402318 -0.3660716543156899 0.5953932791844518 0.43106504879433577 ]
[-0.30655198880585594 0.6606757192118904 -0.21483067305535333 0.3045011114089389 ]
[-0.48053620675695174 -0.11139783377793576 -0.6310958848894725 0.2956864520726446 ]
[-0.5936453158283703 -0.01904935140131578 -0.016510508076204543 -0.79527388379824 ]]
**** Q from np.linalg.qr
[[-0.10345123000824041 0.6455437884382426 0.44810714367794663 -0.039635447112567376]
[-0.5585641540231802 -0.3660716543156898 0.5953932791844523 0.4310650487943359 ]
[-0.30655198880585594 0.6606757192118904 -0.21483067305535375 0.30450111140893893 ]
[-0.48053620675695186 -0.1113978337779356 -0.6310958848894725 0.29568645207264455 ]
[-0.5936453158283704 -0.01904935140131564 -0.0165105080762043 -0.79527388379824 ]]
**** R from qr_decomposition
[[-1.653391466100325 -1.0838054573405895 -1.0632037969249921 -1.1825735233596888 ]
[ 0. 0.7263519982452554 0.7798481878600413 0.5496287509656425 ]
[ 0. 0. -0.26840760341581243 -0.2002757085967938 ]
[ 0. 0. 0. 0.48524469321440966]]
**** R from np.linalg.qr
[[-1.6533914661003253 -1.0838054573405895 -1.0632037969249923 -1.182573523359689 ]
[ 0. 0.7263519982452559 0.7798481878600418 0.5496287509656428]
[ 0. 0. -0.2684076034158126 -0.2002757085967939]
[ 0. 0. 0. 0.4852446932144096]]
除了尾随数字中不可避免的舍入误差外，输出现在匹配。
我目前正在尝试实现矩形矩阵的基于 Householder 的 QR 分解，如 http://eprints.ma.man.ac.uk/1192/1/qrupdating_12nov08.pdf（第 3、4、5 页）中所述。
显然我弄错了一些伪代码，因为 (1) 我的结果不同于 numpy.qr.linalg() 和 (2) 我的例程生成的矩阵 R 不是上三角矩阵。
我的代码（也可在 https://pyfiddle.io/fiddle/afcc2e0e-0857-4cb2-adb5-06ff9b80c9d3/?i=true 下获得）
import math
import argparse
import numpy as np
from typing import Union
def householder(alpha: float, x: np.ndarray) -> Union[np.ndarray, int]:
"""
Computes Householder vector for alpha and x.
:param alpha:
:param x:
:return:
"""
s = math.pow(np.linalg.norm(x, ord=2), 2)
v = x
if s == 0:
tau = 0
else:
t = math.sqrt(alpha * alpha + s)
v_one = alpha - t if alpha <= 0 else -s / (alpha + t)
tau = 2 * v_one * v_one / (s + v_one * v_one)
v /= v_one
return v, tau
def qr_decomposition(A: np.ndarray, m: int, n: int) -> Union[np.ndarray, np.ndarray]:
"""
Applies Householder-based QR decomposition on specified matrix A.
:param A:
:param m:
:param n:
:return:
"""
H = []
R = A
Q = A
I = np.eye(m, m)
for j in range(0, n):
# Apply Householder transformation.
x = A[j + 1:m, j]
v_householder, tau = householder(np.linalg.norm(x), x)
v = np.zeros((1, m))
v[0, j] = 1
v[0, j + 1:m] = v_householder
res = I - tau * v * np.transpose(v)
R = np.matmul(res, R)
H.append(res)
return Q, R
m = 10
n = 8
A = np.random.rand(m, n)
q, r = np.linalg.qr(A)
Q, R = qr_decomposition(A, m, n)
print("*****")
print(Q)
print(q)
print("-----")
print(R)
print(r)
所以我不清楚如何将零引入我的 R matrix/about 我的代码的哪一部分不正确。我会很高兴任何指示！
非常感谢您的宝贵时间。
您链接到的注释中有一堆 problems/missing 详细信息。在咨询 few other sources (including this very useful textbook) 之后，我能够想出一个类似的工作实现。
工作算法
这是 qr_decomposition 工作版本的代码：
import numpy as np
from typing import Union
def householder(x: np.ndarray) -> Union[np.ndarray, int]:
alpha = x[0]
s = np.power(np.linalg.norm(x[1:]), 2)
v = x.copy()
if s == 0:
tau = 0
else:
t = np.sqrt(alpha**2 + s)
v[0] = alpha - t if alpha <= 0 else -s / (alpha + t)
tau = 2 * v[0]**2 / (s + v[0]**2)
v /= v[0]
return v, tau
def qr_decomposition(A: np.ndarray) -> Union[np.ndarray, np.ndarray]:
m,n = A.shape
R = A.copy()
Q = np.identity(m)
for j in range(0, n):
# Apply Householder transformation.
v, tau = householder(R[j:, j])
H = np.identity(m)
H[j:, j:] -= tau * v.reshape(-1, 1) @ v
R = H @ R
Q = H @ Q
return Q[:n].T, R[:n]
m = 5
n = 4
A = np.random.rand(m, n)
q, r = np.linalg.qr(A)
Q, R = qr_decomposition(A)
with np.printoptions(linewidth=9999, precision=20, suppress=True):
print("**** Q from qr_decomposition")
print(Q)
print("**** Q from np.linalg.qr")
print(q)
print()
print("**** R from qr_decomposition")
print(R)
print("**** R from np.linalg.qr")
print(r)
输出：
**** Q from qr_decomposition
[[ 0.5194188817843675 -0.10699353671401633 0.4322294754656072 -0.7293293270703678 ]
[ 0.5218635773595086 0.11737804362574514 -0.5171653705211056 0.04467925806590414]
[ 0.34858177783013133 0.6023104248793858 -0.33329256746256875 -0.03450824948274838]
[ 0.03371048915852807 0.6655221685383623 0.6127023580593225 0.28795294754791 ]
[ 0.5789790833500734 -0.411189947884951 0.24337120818874305 0.618041080584351 ]]
**** Q from np.linalg.qr
[[-0.5194188817843672 0.10699353671401617 0.4322294754656068 0.7293293270703679 ]
[-0.5218635773595086 -0.11737804362574503 -0.5171653705211053 -0.044679258065904115]
[-0.3485817778301313 -0.6023104248793857 -0.33329256746256863 0.03450824948274819 ]
[-0.03371048915852807 -0.665522168538362 0.6127023580593226 -0.2879529475479097 ]
[-0.5789790833500733 0.41118994788495106 0.24337120818874317 -0.6180410805843508 ]]
**** R from qr_decomposition
[[ 0.6894219296137802 1.042676051151294 1.3418719684631446 1.2498925815126485 ]
[ 0.00000000000000000685 0.7076056836914905 0.29883043386651403 0.41955370595004277 ]
[-0.0000000000000000097 -0.00000000000000007292 0.5304551654027297 0.18966088433421135 ]
[-0.00000000000000000662 0.00000000000000008718 0.00000000000000002322 0.6156558913022807 ]]
**** R from np.linalg.qr
[[-0.6894219296137803 -1.042676051151294 -1.3418719684631442 -1.2498925815126483 ]
[ 0. -0.7076056836914905 -0.29883043386651376 -0.4195537059500425 ]
[ 0. 0. 0.53045516540273 0.18966088433421188]
[ 0. 0. 0. -0.6156558913022805 ]]
此版本的 qr_decomposition 几乎完全再现了 np.linalg.qr 的输出。差异在下面进行了评论。
输出的数值精度
np.linalg.qr 和 qr_decomposition 的输出值与高精度匹配。但是，qr_decomposition 用于生成 R 中的零的计算组合并没有完全抵消，因此零实际上并不完全等于零。
事实证明，np.linalg.qr 并没有使用任何花哨的浮点技巧来确保其输出中的零为 0.0。它 just calls np.triu，强制将这些值设置为 0.0。所以要达到相同的结果，只需将 qr_decomposition 中的 return 行更改为：
return Q[:n].T, np.triu(R[:n])
输出中的符号 (+/-)
Q 和 R 中的一些 +/- 符号在 np.linalg.qr 和 qr_decomposition 的输出中是不同的，但这不是真正的问题，因为有很多有效的选择标志（见此discussion of the uniqueness of Q and R）。您可以通过使用替代算法生成 v 和 tau:
np.linalg.qr 的符号约定
def householder_vectorized(a):
"""Use this version of householder to reproduce the output of np.linalg.qr
exactly (specifically, to match the sign convention it uses)
based on https://rosettacode.org/wiki/QR_decomposition#Python
"""
v = a / (a[0] + np.copysign(np.linalg.norm(a), a[0]))
v[0] = 1
tau = 2 / (v.T @ v)
return v,tau
与np.linalg.qr
的输出完全匹配
将它们放在一起，这个版本的 qr_decomposition 将与 np.linalg.qr 的输出完全匹配：
import numpy as np
from typing import Union
def qr_decomposition(A: np.ndarray) -> Union[np.ndarray, np.ndarray]:
m,n = A.shape
R = A.copy()
Q = np.identity(m)
for j in range(0, n):
# Apply Householder transformation.
v, tau = householder_vectorized(R[j:, j, np.newaxis])
H = np.identity(m)
H[j:, j:] -= tau * (v @ v.T)
R = H @ R
Q = H @ Q
return Q[:n].T, np.triu(R[:n])
m = 5
n = 4
A = np.random.rand(m, n)
q, r = np.linalg.qr(A)
Q, R = qr_decomposition(A)
with np.printoptions(linewidth=9999, precision=20, suppress=True):
print("**** Q from qr_decomposition")
print(Q)
print("**** Q from np.linalg.qr")
print(q)
print()
print("**** R from qr_decomposition")
print(R)
print("**** R from np.linalg.qr")
print(r)
输出：
**** Q from qr_decomposition
[[-0.10345123000824041 0.6455437884382418 0.44810714367794663 -0.03963544711256745 ]
[-0.55856415402318 -0.3660716543156899 0.5953932791844518 0.43106504879433577 ]
[-0.30655198880585594 0.6606757192118904 -0.21483067305535333 0.3045011114089389 ]
[-0.48053620675695174 -0.11139783377793576 -0.6310958848894725 0.2956864520726446 ]
[-0.5936453158283703 -0.01904935140131578 -0.016510508076204543 -0.79527388379824 ]]
**** Q from np.linalg.qr
[[-0.10345123000824041 0.6455437884382426 0.44810714367794663 -0.039635447112567376]
[-0.5585641540231802 -0.3660716543156898 0.5953932791844523 0.4310650487943359 ]
[-0.30655198880585594 0.6606757192118904 -0.21483067305535375 0.30450111140893893 ]
[-0.48053620675695186 -0.1113978337779356 -0.6310958848894725 0.29568645207264455 ]
[-0.5936453158283704 -0.01904935140131564 -0.0165105080762043 -0.79527388379824 ]]
**** R from qr_decomposition
[[-1.653391466100325 -1.0838054573405895 -1.0632037969249921 -1.1825735233596888 ]
[ 0. 0.7263519982452554 0.7798481878600413 0.5496287509656425 ]
[ 0. 0. -0.26840760341581243 -0.2002757085967938 ]
[ 0. 0. 0. 0.48524469321440966]]
**** R from np.linalg.qr
[[-1.6533914661003253 -1.0838054573405895 -1.0632037969249923 -1.182573523359689 ]
[ 0. 0.7263519982452559 0.7798481878600418 0.5496287509656428]
[ 0. 0. -0.2684076034158126 -0.2002757085967939]
[ 0. 0. 0. 0.4852446932144096]]
除了尾随数字中不可避免的舍入误差外，输出现在匹配。