QR 分解
QR 分解とは、行列を直交行列と上三角行列の積に分解することである。具体的に、行列の各列にグラムシュミット直交化法を適用すると直交行列が得られ、その直交行列の転置を元の行列に左から乗じると上三角行列が得られる。よって、グラムシュミット直交化法を一意に定めれば QR 分解も一意に定まる。ただし QR 分解を得る方法は他にも色々ある。正方行列でない行列の QR 分解の定義は筆者が興味がないので考えない。
$$
\begin{align}
&
\begin{pmatrix}
a_{1,1} & a_{1,2} & \cdots & a_{1, n} \\
a_{2,1} & a_{2,2} & \cdots & a_{2, n} \\
\vdots &\vdots & \ddots & \vdots \\
a_{n,1} & a_{n,2} & \cdots & a_{n, n} \\
\end{pmatrix}
\xrightarrow[\text{Gram–Schmidt process}]{}
\begin{pmatrix}
q_{1,1} & q_{1,2} & \cdots & q_{1, n} \\
q_{2,1} & q_{2,2} & \cdots & q_{2, n} \\
\vdots &\vdots & \ddots & \vdots \\
q_{n,1} & q_{n,2} & \cdots & q_{n, n} \\
\end{pmatrix}
\\
&
\begin{pmatrix}
q_{1,1} & q_{2,1} & \cdots & q_{n,1} \\
q_{1,2} & q_{2,2} & \cdots & q_{n,2} \\
\vdots &\vdots & \ddots & \vdots \\
q_{1,n} & q_{2,n} & \cdots & q_{n,n} \\
\end{pmatrix}
\begin{pmatrix}
a_{1,1} & a_{1,2} & \cdots & a_{1, n} \\
a_{2,1} & a_{2,2} & \cdots & a_{2, n} \\
\vdots &\vdots & \ddots & \vdots \\
a_{n,1} & a_{n,2} & \cdots & a_{n, n} \\
\end{pmatrix}
=
\begin{pmatrix}
r_{1,1} & r_{1,2} & \cdots & r_{1, n} \\
0 & r_{2,2} & \cdots & r_{2, n} \\
\vdots &\vdots & \ddots & \vdots \\
0 & 0 & \cdots & r_{n, n} \\
\end{pmatrix}
\\
&
\quad
\\
&
\begin{pmatrix}
a_{1,1} & a_{1,2} & \cdots & a_{1, n} \\
a_{2,1} & a_{2,2} & \cdots & a_{2, n} \\
\vdots &\vdots & \ddots & \vdots \\
a_{n,1} & a_{n,2} & \cdots & a_{n, n} \\
\end{pmatrix}
=
\begin{pmatrix}
q_{1,1} & q_{1,2} & \cdots & q_{1, n} \\
q_{2,1} & q_{2,2} & \cdots & q_{2, n} \\
\vdots &\vdots & \ddots & \vdots \\
q_{n,1} & q_{n,2} & \cdots & q_{n, n} \\
\end{pmatrix}
\begin{pmatrix}
r_{1,1} & r_{1,2} & \cdots & r_{1, n} \\
0 & r_{2,2} & \cdots & r_{2, n} \\
\vdots &\vdots & \ddots & \vdots \\
0 & 0 & \cdots & r_{n, n} \\
\end{pmatrix}
\end{align}
$$
参考文献
- QR分解 - Wikipedia. https://ja.wikipedia.org/wiki/QR%E5%88%86%E8%A7%A3（2022年3月21日参照）.
- 本記事はこの文献にある方法を numpy で確認しているだけである。
QR 分解を得る方法
参考文献 [1] と同じ行列に numpy.linalg.qr を適用すると以下である。素直にグラムシュミット直交化を適用するときと比べ $Q$ の符号が反転しているが理由は調べていないのでわからない。
import numpy as np
A = np.array([
[12, -51, 4],
[ 6, 167, -68],
[-4, 24, -41]
])
Q, R = np.linalg.qr(A)
print(Q)
print(R)
[[-0.85714286 0.39428571 0.33142857] [-0.42857143 -0.90285714 -0.03428571] [ 0.28571429 -0.17142857 0.94285714]] [[ -14. -21. 14.] [ 0. -175. 70.] [ 0. 0. -35.]]
グラムシュミットの直交化法による方法
グラムシュミットの直交化法をする関数がないのだろうかと思ったがそれこそが numpy.linalg.qr とのことである。
def qr(A):
A = A.astype(np.float64)
m, n = A.shape
# グラムシュミットの直交化法
Q = []
for j in range(n):
a = A[:,j]
v = a.copy() # コピーしないと A が破壊される
for q in Q:
v -= np.dot(a, q) * q
Q.append(v / np.linalg.norm(v))
Q = np.array(Q).T
print(Q)
# A = QR ⇒ R = Q^T A
R = np.matmul(Q.T, A)
print(R)
qr(A)
[[ 0.85714286 -0.39428571 -0.33142857] [ 0.42857143 0.90285714 0.03428571] [-0.28571429 0.17142857 -0.94285714]] [[ 1.40000000e+01 2.10000000e+01 -1.40000000e+01] [-1.11022302e-16 1.75000000e+02 -7.00000000e+01] [-1.77635684e-15 -1.59872116e-14 3.50000000e+01]]
ハウスホルダー変換による方法
元の行列の1列目を $(r_{1, 1}, 0, \cdots)$ に鏡映し、1行目と1列目を除いた右下ブロックを取り、その行列の1列目を $(r_{2, 2}, 0, \cdots)$ に鏡映し、……と繰り返すことでも QR 分解できる。
$$
\begin{align}
&
\begin{pmatrix}
a_{1,1} & a_{1,2} & a_{1,3} & \cdots & a_{1, n} \\
a_{2,1} & a_{2,2} & a_{2,3} & \cdots & a_{2, n} \\
a_{3,1} & a_{3,2} & a_{3,3} & \cdots & a_{3, n} \\
\vdots &\vdots &\vdots & \ddots & \vdots \\
a_{n,1} & a_{n,2} & a_{n,3} & \cdots & a_{n, n} \\
\end{pmatrix}
\xrightarrow[\text{Householder}]{}
\begin{pmatrix}
r_{1,1} & r_{1,2} & r_{1,3} & \cdots & r_{1, n} \\
0 & a_{2,2}^{(1)} & a_{2,3}^{(1)} & \cdots & a_{2, n}^{(1)} \\
0 & a_{3,2}^{(1)} & a_{3,3}^{(1)} & \cdots & a_{3, n}^{(1)} \\
\vdots &\vdots &\vdots & \ddots & \vdots \\
0 & a_{n,2}^{(1)} & a_{n,3}^{(1)} & \cdots & a_{n, n}^{(1)} \\
\end{pmatrix}
\\
&
\begin{pmatrix}
a_{2,2}^{(1)} & a_{2,3}^{(1)} & \cdots & a_{2, n}^{(1)} \\
a_{3,2}^{(1)} & a_{3,3}^{(1)} & \cdots & a_{3, n}^{(1)} \\
\vdots &\vdots & \ddots & \vdots \\
a_{n,2}^{(1)} & a_{n,3}^{(1)} & \cdots & a_{n, n}^{(1)} \\
\end{pmatrix}
\xrightarrow[\text{Householder}]{}
\begin{pmatrix}
r_{2,2} & r_{2,3} & \cdots & r_{2, n} \\
0 & a_{3,3}^{(2)} & \cdots & a_{3, n}^{(2)} \\
\vdots &\vdots & \ddots & \vdots \\
0 & a_{n,3}^{(2)} & \cdots & a_{n, n}^{(2)} \\
\end{pmatrix}
\end{align}
$$
def qr(A):
m, n = A.shape # ただし m = n を想定
Q = None
A = A.astype(np.float64)
for i in range(m - 1):
a = A[i:,i]
target = np.zeros(m - i)
target[0] = np.linalg.norm(a)
u = a - target
u = u / np.linalg.norm(u)
if i > 0:
u = np.concatenate((np.zeros(i), u))
H = np.identity(m) - 2 * np.matmul(np.reshape(u, (m, 1)), np.reshape(u, (1, m)))
if i == 0:
Q = H
else:
Q = np.matmul(H, Q)
A = np.matmul(H, A)
Q = Q.T
print(Q)
print(A)
print(np.matmul(Q, A))
qr(A)
[[ 0.85714286 -0.39428571 0.33142857] [ 0.42857143 0.90285714 -0.03428571] [-0.28571429 0.17142857 0.94285714]] [[ 1.40000000e+01 2.10000000e+01 -1.40000000e+01] [ 5.50670620e-16 1.75000000e+02 -7.00000000e+01] [-3.01980663e-16 -8.88178420e-15 -3.50000000e+01]] [[ 12. -51. 4.] [ 6. 167. -68.] [ -4. 24. -41.]]