パッケージの宣言
# We need the following packages:
# install.packages(
# "http://mayekawa.in.coocan.jp/Rpackages/lazy.tools_0.1.6.zip", repos=NULL)
# install.packages(
# "http://mayekawa.in.coocan.jp/Rpackages/lazy.mat_0.1.4.zip", repos=NULL)
# install.packages(
# "http://mayekawa.in.coocan.jp/Rpackages/lazy.symbolic_0.1.4.zip", repos=NULL)
# install.packages(
# "http://mayekawa.in.coocan.jp/Rpackages/lazy.symbolic2_0.1.4.zip", repos=NULL)
#
library(lazy.symbolic2)\(n \times n\) の行列 の行列式を \(n-1 \times n-1\) の行列の行列式を用いて 計算する方法を行列式の余因子展開、または、ラプラス展開と呼ぶ。
行列からその \(i\) 行 \(j\) 列を除いた \(n-1 \times n-1\) の小行列を \(\mathbf A_{ij}\) とし、この小行列の行列式 \(\left|\mathbf A_{ij}\right|\) を の \((ij)\) 小行列式 (minor) と呼ぶ。また、\((ij)\) 小行列式に \(ij\) の組み合わせで定まる符号をつけたものを の \((ij)\) 余因子 (cofactor) と呼び、 と書く。すなわち、 \[\begin{align} \mbox{$\widetilde{\mathbf A_{ij}}$} = (-1)^{i+j} \; \left|\mathbf A_{ij}\right| \end{align}\] である。
行列の行列式は、この余因子を用いて、任意の列 \(j\) を用いて \[\begin{align} \text{第 $j$ 列に関する余因子展開} \\ \left|\mathbf A\right| &= \sum_{i=1}^n a_{ij} \; \mbox{$\widetilde{\mathbf A_{ij}}$} \;=\; \sum_{i=1}^n a_{ij} \; (-1)^{i+j} \; \left|\mathbf A_{ij}\right| \end{align}\] または、任意の行 \(i\) を用いて \[\begin{align} \text{第 $i$ 行に関する余因子展開} \\ \left|\mathbf A\right| &= \sum_{j=1}^n a_{ij} \; \mbox{$\widetilde{\mathbf A_{ij}}$} \;=\; \sum_{j=1}^n a_{ij} \; (-1)^{i+j} \; \left|\mathbf A_{ij}\right| \end{align}\] と書くことが出来る。
以下の出力では の \((ij)\) 小行列式 \(\left|\mathbf A_{ij}\right|\) を DetL(A[-i,-j]) と 示している。
# 2 x 2 の行列
Aa2=demomat(2,2,root="a")
det=DetL(Aa2)
printm(Aa2,det)## Aa2 det
## c1 c2 [,1]
## r1 a11 a12 [1,] a11*a22-a12*a21
## r2 a21 a22DetL 関数では、一般に \(n \times n\) の行列の行列式は、 デフォールトでは第 1 行に関する 余因子展開を に \(n-1 \times n-1\) の小行列の余因子を用いて表現した形で 表示される。 余因子展開をどのレベルまで行うかのコントロールを maxlevel で 指定することが出来る。デフォールトでは maxlevel=1 である。
Aa3 <-demomat(3,3,root="a")
# 3 x 3 の行列を 2 x 2 の小行列の余因子を用いて表現
DetL(Aa3)## [,1]
## [1,] "a11*DetL(A[-1,-1])-a12*DetL(A[-1,-2])+a13*DetL(A[-1,-3])"# 最後まで展開した形
d1=DetL(Aa3, maxlevel=2)
printm(Aa3,d1)## Aa3
## c1 c2 c3
## r1 a11 a12 a13
## r2 a21 a22 a23
## r3 a31 a32 a33
##
## d1
## [,1]
## [1,] a11*(a22*a33-a23*a32)-a12*(a21*a33-a23*a31)+a13*(a21*a32-a22*a31)第 2 行に関する余因子展開
# 3 x 3 の行列を 2 x 2 の小行列の余因子を用いて表現
DetL(Aa3, row=2)## [,1]
## [1,] "-a21*DetL(A[-2,-1])+a22*DetL(A[-2,-2])-a23*DetL(A[-2,-3])"# 最後まで展開した形
d2=DetL(Aa3, row=2, maxlevel=2)
printm(d2)## d2
## [,1]
## [1,] -a21*(a12*a33-a13*a32)+a22*(a11*a33-a13*a31)-a23*(a11*a32-a12*a31)第 1 列に関する余因子展開
# 3 x 3 の行列を 2 x 2 の小行列の余因子を用いて表現
DetL(Aa3, row=NULL, col=1)## [,1]
## [1,] "a11*DetL(A[-1,-1])-a21*DetL(A[-2,-1])+a31*DetL(A[-3,-1])"# 最後まで展開した形
d3=DetL(Aa3, row=NULL, col=1, maxlevel=2)
printm(d3)## d3
## [,1]
## [1,] a11*(a22*a33-a23*a32)-a21*(a12*a33-a13*a32)+a31*(a12*a23-a13*a22)数値的検証
数値的行列 An を作成し、Aa3 行列の各要素を対応する An の要素に置き換えて 上記で計算された行列式の数値を計算する。
# Numerical verification
set.seed(1701)
An <- matrix( runif(9), 3, 3 )
printm("R の関数で計算した値", det(An))## R の関数で計算した値 det(An)
## 0.111813110076844printm("Sweep 演算して計算した値", Det(An))## Sweep 演算して計算した値
##
## Det(An)
## (0.00231331307440996-0.0544216749258339*0.665117462631315/0.127393427770585-(0.296445079380646-0.0544216749258339*0.561384594766423/0.127393427770585)*(0.652215607231483-0.766723201842979*0.665117462631315/0.127393427770585)/(0.93759793555364-0.766723201842979*0.561384594766423/0.127393427770585))*(0.93759793555364-0.766723201842979*0.561384594766423/0.127393427770585)*0.127393427770585printm( Eval(d1, Aa3=An), Eval(d2, Aa3=An), Eval(d3, Aa3=An) )## Eval(d1, Aa3 = An) Eval(d2, Aa3 = An) Eval(d3, Aa3 = An)
## [,1] [,1] [,1]
## [1,] 0.111813110076844 [1,] 0.111813110076844 [1,] 0.111813110076844三角行列の行列式は対角要素の積である。
# 2 x 2 の下三角行列
Al=demomat(2,2,root="a", shape="lowert")
det=DetL(Al, maxlevel=3)
printm(Al,det)## Al det
## c1 c2 [,1]
## r1 a11 0 [1,] a11*a22
## r2 a21 a22# 任意の三角行列の行列式
L=demomat(4,4,root="L", shape="lowert")
DetL(L, maxlevel=3)## [,1]
## [1,] "L11*L22*L33*L44"U=demomat(4,4,root="U", shape="uppert")
DetL(U, maxlevel=3)## [,1]
## [1,] "U11*U22*U33*U44"逆行列を求める方法には様々なのもがあるが、ここではガウスの消去法と呼ばれる方法 に関して説明する。 ガウスの消去法は \(\mathbf T_1, \mathbf T_2, \ldots, \mathbf T_r\) を行に関する基本演算を行う行列とした時に、逆行列を \[\begin{eqnarray} \mathbf A^{-1} = \mathbf T_r \times \mathbf T_{r-1} \times \cdots \times \mathbf T_2 \times \mathbf T_1 \tag{GE} \end{eqnarray}\] として求める方法であるが、 \(\mathbf T_k, \; k=1,2,\ldots,r\) は、 \(\mathbf B_0=\mathbf A\) として \[\begin{eqnarray} \mathbf B_{k} = \mathbf T_k \times \mathbf B_{k-1} \end{eqnarray}\] の第 \(k\) 列が第 番目の要素が 1 でそれ以外が 0 となる様に定められる。 すなわち、\(\mathbf T_k\) 行列は単位行列の第 \(k\) 列を以下に示すように書き換えたものとなっている。
このようにして \(\mathbf B_k\) 行列を順次求めていくと、最終的に得られる \(\mathbf B_r\) 行列は単位行列となっており、 \[\begin{eqnarray} \mathbf B_r = \mathbf T_r \times \mathbf T_{r-1} \times \cdots \times \mathbf T_2 \times \mathbf T_1 \times \mathbf A\;=\; \mathbf I \end{eqnarray}\] であるから (GE) 式を得る。
計算機上では 行列の右に同じサイズの単位行列を並べたものを \(\mathbf B_0\) とし、 以下の作業を行えば最終的に \(\mathbf B_n\) が、単位行列の右に を並べたもの となっている。 \[\begin{align} \mathbf B_0 &\leftarrow \text{cbind}(\mathbf A,\mathbf I) \\ \mathbf B_1 &\leftarrow \mathbf T_1 \; \mathbf B_0 \\ \mathbf B_2 &\leftarrow \mathbf T_2 \; \mathbf B_1 \\ : & : \\ \mathbf B_n &\leftarrow \mathbf T_n \; \mathbf B_{n-1}, \hspace{2ex}\hspace{2ex}\mathbf B_n=\text{cbind}(\mathbf I,\mathbf A^{-1}) \end{align}\]
2 x 2 の行列の例
invGE 関数では、上記の \(\mathbf B\) 行列を AI という行列で示している。 また、\(\mathbf B\) 行列の最初の \(n\) 列を AI[,1:n] で、最後の \(n\) 列を AI[,-(1:n)] で 示している。
options(width=150)
# 2 x 2 の行列
n=2
Aa=demomat(n,n,root="a")
res=invGE(Aa, print=1)##
## Initial augumented matrix
## AI
## c1 c2 I1 I2
## r1 a11 a12 1 0
## r2 a21 a22 0 1
##
##
## after pivoting at i = 1 ( 1 )
## T AI[, 1:n] AI[, -(1:n)]
## [,1] [,2] c1 c2 I1 I2
## [1,] 1/a11 0 [1,] 1 1/a11*a12 [1,] 1/a11 0
## [2,] -a21/a11 1 [2,] 0 -a21/a11*a12+a22 [2,] -a21/a11 1
##
##
## after pivoting at i = 2 ( 2 )
## T AI[, 1:n] AI[, -(1:n)]
## [,1] [,2] c1 c2 I1 I2
## [1,] 1 -1/a11*a12/(-a21/a11*a12+a22) [1,] 1 0 [1,] 1/a11+1/a11*a12/(-a21/a11*a12+a22)*a21/a11 -1/a11*a12/(-a21/a11*a12+a22)
## [2,] 0 1/(-a21/a11*a12+a22) [2,] 0 1 [2,] -1/(-a21/a11*a12+a22)*a21/a11 1/(-a21/a11*a12+a22)invAa=res[,-c(1:n)]数値的検証
set.seed(1701)
An <- matrix(rnorm(n*n),n,n)
Eval(Aa%*%invAa, Aa=An) |> round( 4 )## [,1] [,2]
## [1,] 1 0
## [2,] 0 13 x 3 の行列でも大きすぎるので、以下に \(T_1\) 行列のみを示す。
options(width=150)
# 3 x 3 の行列
n=3
Aa=demomat(n,n,root="a")
res=invGE(Aa, print=1, loc=1)##
## Initial augumented matrix
## AI
## c1 c2 c3 I1 I2 I3
## r1 a11 a12 a13 1 0 0
## r2 a21 a22 a23 0 1 0
## r3 a31 a32 a33 0 0 1
##
##
## after pivoting at i = 1 ( 1 )
## T AI[, 1:n] AI[, -(1:n)]
## [,1] [,2] [,3] c1 c2 c3 I1 I2 I3
## [1,] 1/a11 0 0 [1,] 1 1/a11*a12 1/a11*a13 [1,] 1/a11 0 0
## [2,] -a21/a11 1 0 [2,] 0 -a21/a11*a12+a22 -a21/a11*a13+a23 [2,] -a21/a11 1 0
## [3,] -a31/a11 0 1 [3,] 0 -a31/a11*a12+a32 -a31/a11*a13+a33 [3,] -a31/a11 0 1数値的検証
set.seed(1701)
An <- matrix(rnorm(n*n),n,n)
res=invGE(Aa)
invAa=res[,-c(1:n)]
printm( solve(An)-Eval(invAa, Aa=An), fmt=".5")## solve(An) - Eval(invAa, Aa = An)
## [,1] [,2] [,3]
## [1,] -0.00000 0.00000 0.00000
## [2,] 0.00000 0.00000 -0.00000
## [3,] 0.00000 0.00000 0.00000数値的行列を用いた方が分かりやすいかも。
set.seed(1702)
n <- 3
A <- matrix( sample(0:9,n*n,replace=1), n,n )
# result
res <- invGE(A, print=1)##
## Initial augumented matrix
## AI
## c1 c2 c3 I1 I2 I3
## r1 3.00 8.00 3.00 1.00 0.00 0.00
## r2 6.00 8.00 5.00 0.00 1.00 0.00
## r3 5.00 5.00 1.00 0.00 0.00 1.00
##
##
## after pivoting at i = 1 ( 1 )
## T AI[, 1:n] AI[, -(1:n)]
## [,1] [,2] [,3] c1 c2 c3 I1 I2 I3
## [1,] 0.33 0.00 0.00 [1,] 1.00 2.67 1.00 [1,] 0.33 0.00 0.00
## [2,] -2.00 1.00 0.00 [2,] 0.00 -8.00 -1.00 [2,] -2.00 1.00 0.00
## [3,] -1.67 0.00 1.00 [3,] 0.00 -8.33 -4.00 [3,] -1.67 0.00 1.00
##
##
## after pivoting at i = 2 ( 2 )
## T AI[, 1:n] AI[, -(1:n)]
## [,1] [,2] [,3] c1 c2 c3 I1 I2 I3
## [1,] 1.00 0.33 0.00 [1,] 1.00 0.00 0.67 [1,] -0.33 0.33 0.00
## [2,] 0.00 -0.12 0.00 [2,] 0.00 1.00 0.12 [2,] 0.25 -0.12 0.00
## [3,] 0.00 -1.04 1.00 [3,] 0.00 0.00 -2.96 [3,] 0.42 -1.04 1.00
##
##
## after pivoting at i = 3 ( 3 )
## T AI[, 1:n] AI[, -(1:n)]
## [,1] [,2] [,3] c1 c2 c3 I1 I2 I3
## [1,] 1.00 0.00 0.23 [1,] 1.00 0.00 0.00 [1,] -0.24 0.10 0.23
## [2,] 0.00 1.00 0.04 [2,] 0.00 1.00 0.00 [2,] 0.27 -0.17 0.04
## [3,] 0.00 0.00 -0.34 [3,] 0.00 0.00 1.00 [3,] -0.14 0.35 -0.34invA <- res[,-(1:n)]
printm(A, invA, A%*%invA, fmt=".2")## A invA A %*% invA
## [,1] [,2] [,3] I1 I2 I3 I1 I2 I3
## [1,] 3.00 8.00 3.00 [1,] -0.24 0.10 0.23 [1,] 1.00 0.00 0.00
## [2,] 6.00 8.00 5.00 [2,] 0.27 -0.17 0.04 [2,] -0.00 1.00 0.00
## [3,] 5.00 5.00 1.00 [3,] -0.14 0.35 -0.34 [3,] 0.00 -0.00 1.00第 (1,1) 要素から始める必要は無い。以下に \(B_1\) 行列の第 2 列を単位ベクトル とするような変換を示す。
options(width=150)
res=invGE(Aa, print=1, loc=2)##
## Initial augumented matrix
## AI
## c1 c2 c3 I1 I2 I3
## r1 a11 a12 a13 1 0 0
## r2 a21 a22 a23 0 1 0
## r3 a31 a32 a33 0 0 1
##
##
## after pivoting at i = 1 ( 2 )
## T AI[, 1:n] AI[, -(1:n)]
## [,1] [,2] [,3] c1 c2 c3 I1 I2 I3
## [1,] 1 -a12/a22 0 [1,] a11-a12/a22*a21 0 a13-a12/a22*a23 [1,] 1 -a12/a22 0
## [2,] 0 1/a22 0 [2,] 1/a22*a21 1 1/a22*a23 [2,] 0 1/a22 0
## [3,] 0 -a32/a22 1 [3,] -a32/a22*a21+a31 0 -a32/a22*a23+a33 [3,] 0 -a32/a22 1