どうも、木村（@kimu3_slime）です。
今回は、Julia（SymPy）で行列の標準形、ランク、LU分解、線形独立性を計算する方法を紹介します。
準備
SymPyを使うので、持っていなければインストールしておきましょう。
1 2 | using Pkg Pkg.add("SymPy") |
準備として、以下のコードを実行しておきます。
1 | using SymPy |
行列の標準形、ランク、LU分解
「行列.echelon_form()」で行階段標準形（row echelon form）、「行列.rank()」でそのランク（階数）が求められます。
1 2 3 | B = Sym[1 2; 3 4] B.echelon_form() B.rank() |
\[ \begin{aligned}\left[ \begin{array}{rr}1&2\\3&4\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rr}1&2\\0&-2\end{array}\right]\end{aligned} \]
\[ \begin{aligned}2\end{aligned} \]
「行列.rref()」で、簡約な行階段標準形（reduced row-echelon form）です。
1 | B.rref() |
1 | (Sym[1 0; 0 1], (0, 1)) |
戻り値の第2成分は、掃き出しに使った列（要 pivot）の添字です。Python流の書き方になっており、1列、2列目の掃き出しとなっています。Julia流でない（(1,2)でない）のはわかりにくいですね。
「行列.LUdecomposition()」でLU分解が求められます。
1 | B.LUdecomposition() |
戻り値の第3成分は、置換行列\(P\)を表しています。一般に\(PA=LU\)となる分解は、PLU分解と呼ばれるものです。今回の結果は空なので、気にしなくて良いでしょう。
\(L,U\)はそれぞれ下三角、上三角行列であり、その積はきちんと元の行列になっています。
1 2 3 | B.LUdecomposition()[1] B.LUdecomposition()[2] B.LUdecomposition()[1] *B.LUdecomposition()[2] |
\[ \begin{aligned}\left[ \begin{array}{rr}1&0\\3&1\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rr}1&2\\0&-2\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rr}1&2\\3&4\end{array}\right]\end{aligned} \]
記号を含む行列に対して、階段標準形、ランク、LU分解を求めてみます。
1 2 3 4 5 | @vars a11 a12 a13 a21 a22 a23 a31 a32 a33 x y z n α β R_α= Sym[cos(α) -sin(α); sin(α) cos(α)] R_α.echelon_form() R_α.rank() R_α.LUdecomposition() |
\[ \begin{aligned}\left[ \begin{array}{rr}\cos{\left(α \right)}&- \sin{\left(α \right)}\\\sin{\left(α \right)}&\cos{\left(α \right)}\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rr}\cos{\left(α \right)}&- \sin{\left(α \right)}\\0&1\end{array}\right]\end{aligned} \]
\[ \begin{aligned}2\end{aligned} \]
回転行列のランクは2であり、可逆であることがわかります。
上の表示式では、\(\cos \alpha =0\)となるケースでランクが落ちるように見えますが、それが起こらない（別のランク2の階段形となる）ことを代入によってチェックしましょう。
1 2 | R_α.subs(α,PI/2).echelon_form() R_α.subs(α,3/2*PI).echelon_form() |
\[ \begin{aligned}\left[ \begin{array}{rr}1&0\\0&-1\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rr}-1&0\\0&1\end{array}\right]\end{aligned} \]
1 2 3 | G = Sym[x y z; y z x; z x y] G.echelon_form() G.rank() |
\[ \begin{aligned}\left[ \begin{array}{rrr}x&y&z\\y&z&x\\z&x&y\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rrr}x&y&z\\0&x z – y^{2}&x^{2} – y z\\0&0&- x^{4} + 3 x^{2} y z – x y^{3} – x z^{3}\end{array}\right]\end{aligned} \]
\[ \begin{aligned}3\end{aligned} \]
階段形にしたときに、文字式を含む成分は0でないものとしてランクが計算されるようです。この結果は、変数の値によっては正しくないケースもあることに注意しましょう。
例えば、\(x=y=z\)のときを考えると、
1 2 3 | G.subs([(y,x),(z,x)]) G.subs([(y,x),(z,x)]).echelon_form() G.subs([(y,x),(z,x)]).rank() |
\[ \begin{aligned}\left[ \begin{array}{rrr}x&x&x\\x&x&x\\x&x&x\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rrr}x&x&x\\0&0&0\\0&0&0\end{array}\right]\end{aligned} \]
\[ \begin{aligned}1\end{aligned} \]
ランクは1に落ちています。
線形方程式\(Ax=b\)が解を持つことは、\(A\)とその拡大係数行列のランクが等しいことが必要十分であることが知られています。
このような問題は、さきほどまでの方法で簡単に解けますね。
1 2 3 4 5 6 | C = Sym[1 2 1;2 1 3;1 3 2] C.echelon_form() C.rank() D = Sym[1 2 1 3;2 1 3 1;1 3 2 2] D.echelon_form() D.rank() |
\[ \begin{aligned}\left[ \begin{array}{rrr}1&2&1\\2&1&3\\1&3&2\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rrr}1&2&1\\0&-3&1\\0&0&-4\end{array}\right]\end{aligned} \]
\[ \begin{aligned}3\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rrrr}1&2&1&3\\2&1&3&1\\1&3&2&2\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rrrr}1&2&1&3\\0&-3&1&-5\\0&0&-4&8\end{array}\right]\end{aligned} \]
\[ \begin{aligned}3\end{aligned} \]
ランクがともに3なので、これらによって定まる線形方程式は解を持つことがわかりました。
線形独立性の判定
与えられたベクトルの組が線形独立か線形従属であるかは、それを並べてできる行列のランクがベクトルの個数に等しいかで判定できることが知られています。
手計算では少し手間がかかりますが、簡単に判定できます。やってみましょう。
\[ \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} \]
1 2 3 | E = Sym[1 2 2;2 1 1;1 2 1] E.echelon_form() E.rank() |
\[ \begin{aligned}\left[ \begin{array}{rrr}1&2&2\\2&1&1\\1&2&1\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rrr}1&2&2\\0&-3&-3\\0&0&-1\end{array}\right]\end{aligned} \]
\[ \begin{aligned}3\end{aligned} \]
ベクトルの個数は3、ランクは3なので、これらは線形独立です。
\[ \begin{aligned}b_1,b_2,b_3,b_4 =\begin{pmatrix} 2\\ 0 \\-3 \\-4\end{pmatrix},\begin{pmatrix} 0\\ 3 \\1 \\ 2\end{pmatrix}, \begin{pmatrix} 1\\ 2 \\1 \\ 0 \end{pmatrix} , \begin{pmatrix} 0\\ -1 \\-4 \\ -2 \end{pmatrix}\end{aligned} \]
1 2 3 | F = Sym[2 0 1 0; 0 3 2 -1;-3 1 1 -4; -4 2 0 -2] F.echelon_form() F.rank() |
\[ \begin{aligned}\left[ \begin{array}{rrrr}2&0&1&0\\0&3&2&-1\\-3&1&1&-4\\-4&2&0&-2\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rrrr}2&0&1&0\\0&3&2&-1\\0&0&11&-22\\0&0&0&0\end{array}\right]\end{aligned} \]
\[ \begin{aligned}3\end{aligned} \]
ベクトルの個数は4、ランクは3と等しくないので、これらは線形従属です。
次元定理によると、\((m,n)\)行列に対し\(\operatorname{rank} (A) + \dim( \ker (A))=n\)が成り立ちます。そして、核の次元\(\dim( \ker (A))\)を退化次数（nullity）と呼びます。それを求めてみましょう。
1 2 3 4 5 | function nullity(A::Matrix{Sym}) return A.shape[2]-A.rank() end nullity(F) |
さきほどの4つのベクトルが定める行列の退化次数は1であり、4次元の空間を3次元へ写して1次元分落としていることがわかります。
核や像、その次元の求め方については、別記事でも紹介予定です。
以上、Julia（SymPy）で行列の標準形、ランク、LU分解、線形独立性を計算する方法を紹介してきました。
階段標準形やランクを求める方法は、線形方程式の可解性や線形独立性の判定など応用が幅広いです。手計算だと計算ミスしやすいので、コンピュータで計算できると嬉しいですね。
木村すらいむ（@kimu3_slime）でした。ではでは。
コロナ社 (2020-03-26T00:00:01Z)
￥7,353 (コレクター商品)
こちらもおすすめ
線形方程式の解き方：ガウスの消去法と基本変形・ランク、LU分解
