どうも、木村（@kimu3_slime）です。
今回は、Julia（SymPy）でQR分解、正規直交系を計算する方法を紹介します。
準備
SymPyを使うので、持っていなければインストールしておきましょう。
1 2 | using Pkg Pkg.add("SymPy") |
準備として、以下のコードを実行しておきます。
1 | using SymPy |
QR分解、正規直交系を計算する方法
\[ \begin{aligned}a_1,a_2,a_3 =\begin{pmatrix} 1\\ 2 \\1\end{pmatrix},\begin{pmatrix} 2\\ 1 \\2\end{pmatrix}, \begin{pmatrix} 2\\ 1 \\1\end{pmatrix}\end{aligned} \]
といったベクトルを、グラム・シュミットの方法によって正規直交化したいとしましょう。
それはこれらのベクトルを並べた行列を作り、それをQR分解することで得られます。「行列.QRdecomposition()」でQR分解です。
1 2 3 4 | A = Sym[1 2 2; 2 1 1; 1 2 1] A.QRdecomposition() A.QRdecomposition()[1] A.QRdecomposition()[2] |
\[ \begin{aligned}\left[ \begin{array}{rrr}1&2&2\\2&1&1\\1&2&1\end{array}\right]\end{aligned} \]
1 | (Sym[sqrt(6)/6 sqrt(3)/3 sqrt(2)/2; sqrt(6)/3 -sqrt(3)/3 0; sqrt(6)/6 sqrt(3)/3 -sqrt(2)/2], Sym[sqrt(6) sqrt(6) 5*sqrt(6)/6; 0 sqrt(3) 2*sqrt(3)/3; 0 0 sqrt(2)/2]) |
\[ \begin{aligned}\left[ \begin{array}{rrr}\frac{\sqrt{6}}{6}&\frac{\sqrt{3}}{3}&\frac{\sqrt{2}}{2}\\\frac{\sqrt{6}}{3}&- \frac{\sqrt{3}}{3}&0\\\frac{\sqrt{6}}{6}&\frac{\sqrt{3}}{3}&- \frac{\sqrt{2}}{2}\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rrr}\sqrt{6}&\sqrt{6}&\frac{5 \sqrt{6}}{6}\\0&\sqrt{3}&\frac{2 \sqrt{3}}{3}\\0&0&\frac{\sqrt{2}}{2}\end{array}\right]\end{aligned} \]
戻り値の第1成分が直交行列\(Q\)、第2成分が上三角行列\(R\)で、\(A=QR\)という関係になっています。
\(Q\)の列ベクトルが、ベクトル\(a_1,a_2,a_3\)を正規直交化したものです。
1 2 3 | A.QRdecomposition()[1][1:3,1] A.QRdecomposition()[1][1:3,2] A.QRdecomposition()[1][1:3,3] |
\[ \begin{aligned}\left[ \begin{array}{r}\frac{\sqrt{6}}{6}\\\frac{\sqrt{6}}{3}\\\frac{\sqrt{6}}{6}\end{array} \right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{r}\frac{\sqrt{3}}{3}\\- \frac{\sqrt{3}}{3}\\\frac{\sqrt{3}}{3}\end{array} \right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{r}\frac{\sqrt{2}}{2}\\0\\- \frac{\sqrt{2}}{2}\end{array} \right]\end{aligned} \]
これらが正規直交系であること、\(Q^{\top} Q =I\)となることは計算で確かめられます。
1 | A.QRdecomposition()[1].T * A.QRdecomposition()[1] |
\[ \begin{aligned}\left[ \begin{array}{rrr}1&0&0\\0&1&0\\0&0&1\end{array}\right]\end{aligned} \]
確かに\(QR =A\)となっていることも、計算で検証できます。
1 | A.QRdecomposition()[1] * A.QRdecomposition()[2] |
\[ \begin{aligned}\left[ \begin{array}{rrr}1&2&2\\2&1&1\\1&2&1\end{array}\right]\end{aligned} \]
別の例として、4次のアダマール行列をQR分解してみましょう。
1 2 3 | B = Sym[1 1 1 1; 1 -1 1 -1; 1 1 -1 -1; 1 -1 -1 1] B.QRdecomposition()[1] B.QRdecomposition()[2] |
\[ \begin{aligned}\left[ \begin{array}{rrrr}1&1&1&1\\1&-1&1&-1\\1&1&-1&-1\\1&-1&-1&1\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rrrr}\frac{1}{2}&\frac{1}{2}&\frac{1}{2}&\frac{1}{2}\\\frac{1}{2}&- \frac{1}{2}&\frac{1}{2}&- \frac{1}{2}\\\frac{1}{2}&\frac{1}{2}&- \frac{1}{2}&- \frac{1}{2}\\\frac{1}{2}&- \frac{1}{2}&- \frac{1}{2}&\frac{1}{2}\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rrrr}2&0&0&0\\0&2&0&0\\0&0&2&0\\0&0&0&2\end{array}\right]\end{aligned} \]
検算してみると、
1 2 | B.QRdecomposition()[1].T * B.QRdecomposition()[1] B.QRdecomposition()[1]*B.QRdecomposition()[2] |
\[ \begin{aligned}\left[ \begin{array}{rrrr}1&0&0&0\\0&1&0&0\\0&0&1&0\\0&0&0&1\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rrrr}1&1&1&1\\1&-1&1&-1\\1&1&-1&-1\\1&-1&-1&1\end{array}\right]\end{aligned} \]
QR分解がきちんと得られていることがわかります。
以上、Julia（SymPy）でQR分解、正規直交系を計算する方法を紹介してきました。
シュミットの直交化法は、手計算だとそれなりに計算量があるので、コンピュータのQR分解で求められると嬉しいですね。
木村すらいむ（@kimu3_slime）でした。ではでは。
コロナ社 (2020-03-26T00:00:01Z)
￥7,353 (コレクター商品)
こちらもおすすめ
直交ベクトルの線形独立性、正規直交基底、直交行列について解説
