lazy.symbolic
パッケージを利用しての
記号的行列演算
この文書は R の Markdown で作成しています。灰色に囲まれた部分が R のコードで、 白枠の部分がその出力です。コードをコピーして R に貼り付けると同じ結果が 出力されるはずです。ただし、文書の最初から順に実行しないと、未定義のオブジェクト やパッケージが出てくる可能性があります。
コピーをするためには灰色部分の左上のボタンをクリックして下さい。
主なパッケージは
lazy.symbolic
ですが、以下のパッケージ
lazy.mat, and lazy.tools
も必要となります。これらのパッケージはここに置いてあります。
http://mayekawa.in.coocan.jp/Rpackages/Rpackages.html
また、CRAN の MASS パッケージも必要です。
パッケージの宣言
library(lazy.symbolic)（ここでいくつかのメッセージが出ますが省略してあります。）
このパッケージを利用すると、通常の R とは異なり、 記号を要素として持つ行列やベクトル間の四則演算をふくむいろいろな演算が可能です。
まず、以下の様に modify_math_operators 関数を実行して、R の数値演算子
+, -, *, /, %*%, t が以下に示す記号行列を取り扱えるようにします。
# set output line width
options(width=500)
# 数値演算子の拡張:
modify_math_operators()##
##
## **
## The following mathematical operators have been modified
## and made available in the following environment.
## <environment: R_GlobalEnv>これらの演算子を使うと、たとえば、
( A=demomat(2,3,root="a") )## c1 c2 c3
## r1 "a11" "a12" "a13"
## r2 "a21" "a22" "a23"( b=demomat(3,1,root="b") )## c1
## r1 "b11"
## r2 "b21"
## r3 "b31"( A %*% b )## c1
## r1 "a11*b11+a12*b21+a13*b31"
## r2 "a21*b11+a22*b21+a23*b31"の様なことが可能となります。
以下では、この記号的な行列を使って行列演算を学習しますが、 文字列（記号）を要素として持つ行列を記号行列 (semi-symbolic matrix)、 要素は定義されていないが、行列と見なすべきものを完全記号行列 (fully symbolic matrix) と呼びます。
ベクトルとは数値を縦一列もしくは横一行に並べたもので、 前者を縦ベクトル（もしくは列ベクトル）、後者は横ベクトル（もしくは行ベクトル） と呼びます。
行列はサイズが同じベクトルを並べたものです。 ただし、サイズが同じという場合には、縦ベクトルと横ベクトルを区別します。 縦ベクトルの左右には縦ベクトルが、横ベクトルの上下には横ベクトルが並びます。
パイプ演算子は優先順位が高いので括弧が必要です。
# A=3+4 を印刷したい場合。
# wrong!
A=3+4 |> print()## [1] 4# correct
(A=3+4) |> print()## [1] 7demomat 関数を用いて記号行列を定義することが出来る。 demomat 関数を使わず、中身を直接指定して記号的行列を作成することも可能。
# set output line width
options(width=500)
# 行列とベクトルの定義
#
# ベクトル
v <- demomat(3,1,root="v",vec=1) # 3 x 1 の縦ベクトル
u <- demomat(3,1,root="u",vec=1)
u4 <- demomat(4,1,root="u",vec=1) # 4 x 1 の縦ベクトル
v4 <- demomat(4,1,root="v",vec=1)
# 定数
one3 <- matrix("1",3,1) # col vector of length 3 consisting of 1s.
one4 <- matrix("1",4,1) # col vector of length 4 consisting of 1s.
zero3 <- matrix("0",3,1) # col vector of length 3 consisting of 0s.
zero4 <- matrix("0",4,1) # col vector of length 4 consisting of 0s.
# 行列
A <- demomat(4,3,root="a") # 4 x 3 の行列
B <- demomat(4,3,root="b")
X <- demomat(3,2,root="x") # 3 x 2 の行列
Y <- demomat(3,2,root="y")
Z <- demomat(3,2,root="z")
Y2 <- demomat(2,3,root="y") # 2 x 3 の行列
Z2 <- demomat(3,3,root="y") # 3 x 3 の行列
# 直接に要素を定義して 2 x 3 の行列を作る
X2 <- matrix( paste("x",outer(1:2,1:3, paste, sep=""), sep=""), 2,3 )
# 対角行列
D3 <- demomat(3,3,root="d", shape="diag") # 3 x 3 の対角行列
D4 <- demomat(4,4,root="d", shape="diag") # 4 x 4 の対角行列
# 対称行列
S <- demomat(3,3,root="s", shape="sym") # 3 x 3 の対称行列
Su <- demomat(3,3,root="s", shape="symU") # 3 x 3 の対称行列
# 三角行列
L <- demomat(3,3,root="s", shape="lowert") # 3 x 3 の下三角行列
U <- demomat(3,3,root="s", shape="uppert", nodiag=1) # 3 x 3 の上三角行列
# 単位行列
I3 <- diag(3)
I4 <- diag(4)
# 零行列
Zero3 <- matrix("0",3,3)これらの行列は、通常の数値行列と同様に、コンソールにその名前をタイプすれば 表示される。また、以下の様に printm または Print 関数等を用いることで、 複数の行列を名前付きで出力することが出来る。
# set output line width
options(width=500)
X; X2## c1 c2
## r1 "x11" "x12"
## r2 "x21" "x22"
## r3 "x31" "x32"## [,1] [,2] [,3]
## [1,] "x11" "x12" "x13"
## [2,] "x21" "x22" "x23"printm(A, D3, I4, Zero3)## A D3 I4 Zero3
## c1 c2 c3 [,1] [,2] [,3] [,1] [,2] [,3] [,4] [,1] [,2] [,3]
## r1 a11 a12 a13 [1,] d11 0 0 [1,] 1 0 0 0 [1,] 0 0 0
## r2 a21 a22 a23 [2,] 0 d22 0 [2,] 0 1 0 0 [2,] 0 0 0
## r3 a31 a32 a33 [3,] 0 0 d33 [3,] 0 0 1 0 [3,] 0 0 0
## r4 a41 a42 a43 [4,] 0 0 0 1printm( S, Su, L, U)## S Su L U
## c1 c2 c3 c1 c2 c3 c1 c2 c3 c1 c2 c3
## r1 s11 s21 s31 r1 s11 s12 s13 r1 s11 0 0 r1 0 s12 s13
## r2 s21 s22 s32 r2 s12 s22 s23 r2 s21 s22 0 r2 0 0 s23
## r3 s31 s32 s33 r3 s13 s23 s33 r3 s31 s32 s33 r3 0 0 0identical(X,X2)## [1] FALSE数値行列と同様に記号行列の要素を抽出することが出来る。
# set output line width
options(width=500)
# これは R native のベクトル
(xv_1=X[,1])## r1 r2 r3
## "x11" "x21" "x31"# これらは 3 x 1 の縦ベクトル
xv1=X[,1,drop=0]
xv2=X[,2,drop=0]
printm(xv_1,xv1,xv2) # 並べて印刷するとほとんど同じ。## xv_1 xv1 xv2
## c1 c2
## r1 x11 r1 x11 r1 x12
## r2 x21 r2 x21 r2 x22
## r3 x31 r3 x31 r3 x32# これらから元のデータ行列を再現
X2=cbind(xv1,xv2)
# チェック
identical(X,X2)## [1] TRUE対角化演算子 Diag は行列の非対角要素を 0 とする。
転置演算子 t は、正方行列をその主対角線 （左上から右下の (1,1), (2,2), (3,3) … 要素を繋ぐ線）の周りに回転する。 すなわち、縦長行列は横長に、横長行列は縦長に、縦ベクトルは横ベクトルに、 横ベクトルは縦ベクトルになる。
# set output line width
options(width=500)
# 対角化演算子（対角要素からなる対角行列を作成する。diag(diag(S)) は使えない。）
Diag(S)## [,1] [,2] [,3]
## [1,] "s11" "0" "0"
## [2,] "0" "s22" "0"
## [3,] "0" "0" "s33"# 対角要素を R のベクトルとして抽出
diag(S)## [1] "s11" "s22" "s33"# 転置演算
printm(v,t(v), t(t(v)))## v t(v) t(t(v))
## [,1] [,1] [,2] [,3] [,1]
## [1,] v1 [1,] v1 v2 v3 [1,] v1
## [2,] v2 [2,] v2
## [3,] v3 [3,] v3printm(A,t(A), t(t(A)))## A t(A) t(t(A))
## c1 c2 c3 r1 r2 r3 r4 c1 c2 c3
## r1 a11 a12 a13 c1 a11 a21 a31 a41 r1 a11 a12 a13
## r2 a21 a22 a23 c2 a12 a22 a32 a42 r2 a21 a22 a23
## r3 a31 a32 a33 c3 a13 a23 a33 a43 r3 a31 a32 a33
## r4 a41 a42 a43 r4 a41 a42 a43printm(S,t(S))## S t(S)
## c1 c2 c3 r1 r2 r3
## r1 s11 s21 s31 c1 s11 s21 s31
## r2 s21 s22 s32 c2 s21 s22 s32
## r3 s31 s32 s33 c3 s31 s32 s33行列は、ベクトルを並べたものと見なすことが出来る。
# set output line width
options(width=500)
#
# 行列を列ベクトルを cbind 関数で横に並べたものとみる場合。
# cbind は（縦、列）ベクトルを横に並べて行列を作る R の関数
E=cbind(u,v)
printm(u,v,E)## u v E
## [,1] [,1] [,1] [,2]
## [1,] u1 [1,] v1 [1,] u1 v1
## [2,] u2 [2,] v2 [2,] u2 v2
## [3,] u3 [3,] v3 [3,] u3 v3# その転置は
t(E)## [,1] [,2] [,3]
## [1,] "u1" "u2" "u3"
## [2,] "v1" "v2" "v3"# それぞれの縦ベクトルの転置を rbind 関数で縦に並べたものとなる。
rbind( t(u), t(v) )## [,1] [,2] [,3]
## [1,] "u1" "u2" "u3"
## [2,] "v1" "v2" "v3"# 行列を行ベクトルを rbind 関数で縦に並べたものとみる場合。
# rbind は（横、行）ベクトルを縦に並べて行列を作る R の関数
F=rbind(t(u4), t(v4))
printm(u4,v4,F)## u4 v4 F
## [,1] [,1] [,1] [,2] [,3] [,4]
## [1,] u1 [1,] v1 [1,] u1 u2 u3 u4
## [2,] u2 [2,] v2 [2,] v1 v2 v3 v4
## [3,] u3 [3,] v3
## [4,] u4 [4,] v4# その転置は
t(F)## [,1] [,2]
## [1,] "u1" "v1"
## [2,] "u2" "v2"
## [3,] "u3" "v3"
## [4,] "u4" "v4"# それぞれのベクトルの転置（縦ベクトル）を cbind 関数で横に並べたものとなる。
cbind( u4, v4 )## [,1] [,2]
## [1,] "u1" "v1"
## [2,] "u2" "v2"
## [3,] "u3" "v3"
## [4,] "u4" "v4"行列の行に個体を列に変数を配置したものをデータ行列と呼ぶ。 行列の第 i 行には個体 i の全変数による測定値が含まれる。 行列の第 i 行は横ベクトルであるが、 データサイエンスでは、それを転置した縦ベクトルを個体 i の観測値ベクトルとする 場合が多い。 逆に行列の第 j 列には変数 j による全個体の測定値が入るがそれもやはり 縦ベクトルとして取り扱う。
# set output line width
options(width=500)
# データ行列（個体 x 変数）
# 各個体の観測値(縦ベクトル）を転置して得られる行ベクトルをを縦に並べたもの
# X=demomat(3,2,root="x")
# 個体のデータの抽出
# これは R native のベクトル
(x1_0=X[1,])## c1 c2
## "x11" "x12"# これらは 2 x 1 の縦ベクトル
x1=t( X[1,,drop=0] )
x2=t( X[2,,drop=0] )
x3=t( X[3,,drop=0] )
printm(x1,x2,x3)## x1 x2 x3
## r1 r2 r3
## c1 x11 c1 x21 c1 x31
## c2 x12 c2 x22 c2 x32# これらから元のデータ行列を再現
X1=rbind(t(x1), t(x2), t(x3))
# チェック
X-X1## c1 c2
## r1 "0" "0"
## r2 "0" "0"
## r3 "0" "0"identical(X,X1)## [1] TRUE# 個体 i のデータ
for( i in 1:nrow(X) ){
xi=t(X[i,,drop=0])
printm(i, xi)
}## i xi
## r1
## 1 c1 x11
## c2 x12
##
## i xi
## r2
## 2 c1 x21
## c2 x22
##
## i xi
## r3
## 3 c1 x31
## c2 x32# 変数の抽出
# これは R native のベクトル
(xv_1=X[,1])## r1 r2 r3
## "x11" "x21" "x31"# これらは 3 x 1 の縦ベクトル
xv1=X[,1,drop=0]
xv2=X[,2,drop=0]
printm(xv1,xv2)## xv1 xv2
## c1 c2
## r1 x11 r1 x12
## r2 x21 r2 x22
## r3 x31 r3 x32# これらから元のデータ行列を再現
X2=cbind(xv1,xv2)
# チェック
X-X2## c1 c2
## r1 "0" "0"
## r2 "0" "0"
## r3 "0" "0"identical(X,X2)## [1] TRUE# 変数 j
for( j in 1:ncol(X) ){
xvj=X[,j,drop=0]
printm(j, xvj)
}## j xvj
## c1
## 1 r1 x11
## r2 x21
## r3 x31
##
## j xvj
## c2
## 2 r1 x12
## r2 x22
## r3 x32以下にいくつかの行列演算を定義する。
ベクトルや行列にスカラーの定数 c を掛けたものは、 元のベクトルや行列の要素をその定数倍したものとなる。
表記としては定数が左に来る場合が多いが、交換の法則が成り立つため、右からスカラーを掛けても良い。要素の掛け算の記号は * を用いる。
# set output line width
options(width=500)
# 定数倍 c をスカラーの定数とする。
"c" * v## [,1]
## [1,] "c*v1"
## [2,] "c*v2"
## [3,] "c*v3"3 * v## [,1]
## [1,] "3*v1"
## [2,] "3*v2"
## [3,] "3*v3"-1 * v## [,1]
## [1,] "-v1"
## [2,] "-v2"
## [3,] "-v3"v * 3## [,1]
## [1,] "3*v1"
## [2,] "3*v2"
## [3,] "3*v3""c" * A## c1 c2 c3
## r1 "c*a11" "c*a12" "c*a13"
## r2 "c*a21" "c*a22" "c*a23"
## r3 "c*a31" "c*a32" "c*a33"
## r4 "c*a41" "c*a42" "c*a43"A * 4## c1 c2 c3
## r1 "4*a11" "4*a12" "4*a13"
## r2 "4*a21" "4*a22" "4*a23"
## r3 "4*a31" "4*a32" "4*a33"
## r4 "4*a41" "4*a42" "4*a43"サイズが同じベクトルや行列の加算と減算は、要素同士の加算や減算となる。
# set output line width
options(width=500)
# サイズが同じベクトル同士
v + u## [,1]
## [1,] "v1+u1"
## [2,] "v2+u2"
## [3,] "v3+u3"v - u## [,1]
## [1,] "v1-u1"
## [2,] "v2-u2"
## [3,] "v3-u3"# サイズが同じ行列同士
A + B## c1 c2 c3
## r1 "a11+b11" "a12+b12" "a13+b13"
## r2 "a21+b21" "a22+b22" "a23+b23"
## r3 "a31+b31" "a32+b32" "a33+b33"
## r4 "a41+b41" "a42+b42" "a43+b43"A - B## c1 c2 c3
## r1 "a11-b11" "a12-b12" "a13-b13"
## r2 "a21-b21" "a22-b22" "a23-b23"
## r3 "a31-b31" "a32-b32" "a33-b33"
## r4 "a41-b41" "a42-b42" "a43-b43"# 全ての要素が 0 である行列をゼロ行列と呼ぶ。
A - A## c1 c2 c3
## r1 "0" "0" "0"
## r2 "0" "0" "0"
## r3 "0" "0" "0"
## r4 "0" "0" "0"# 対称行列
S - t(S)## c1 c2 c3
## r1 "0" "0" "0"
## r2 "0" "0" "0"
## r3 "0" "0" "0"サイズが同じ二つのベクトルの対応する要素を掛け合わせて加えたたものを ベクトルの内積と呼ぶ。記法としては、縦ベクトルの左から横ベクトルが掛かった 形、すなわち、横ベクトルの右から縦ベクトルが掛かった形となるが、 掛け算の記号としては %*% を用いる。
# set output line width
options(width=500)
# ベクトルの内積もしくはスカラー積
t(u) %*% v## [,1]
## [1,] "u1*v1+u2*v2+u3*v3"t(v) %*% u## [,1]
## [1,] "v1*u1+v2*u2+v3*u3"# これは、関数型言語(Higher-Order Functions in Functional Programming Languages)
# を用いると以下の様に書くことが出来る。
Reduce ("+", u*v )## [,1]
## [1,] "u1*v1+u2*v2+u3*v3"Sum(u*v)## [,1]
## [1,] "u1*v1+u2*v2+u3*v3"# set output line width
options(width=500)
# ベクトルの要素の和
t(one3) %*% v## [,1]
## [1,] "v1+v2+v3"t(v) %*% one3## [,1]
## [1,] "v1+v2+v3"# ベクトルの要素の２乗和: 自分自身との内積
t(v) %*% v## [,1]
## [1,] "v1*v1+v2*v2+v3*v3"t(one3) %*% one3 # ちゃんと計算すれば 3 になるがそこまでは出来ない。## [,1]
## [1,] "2*1+1"縦ベクトルの右から横ベクトルを掛けたものは サイズが（縦ベクトルの行数 x 横ベクトルの列数）の行列となる。
なお、outer product は外積と訳せそうだが、 日本語の外積は英語では cross product に相当し、直積とは異なることに注意。
# set output line width
options(width=500)
u %*% t(v)## [,1] [,2] [,3]
## [1,] "u1*v1" "u1*v2" "u1*v3"
## [2,] "u2*v1" "u2*v2" "u2*v3"
## [3,] "u3*v1" "u3*v2" "u3*v3"v %*% t(u4)## [,1] [,2] [,3] [,4]
## [1,] "v1*u1" "v1*u2" "v1*u3" "v1*u4"
## [2,] "v2*u1" "v2*u2" "v2*u3" "v2*u4"
## [3,] "v3*u1" "v3*u2" "v3*u3" "v3*u4"one3 %*% t(one3)## [,1] [,2] [,3]
## [1,] "1" "1" "1"
## [2,] "1" "1" "1"
## [3,] "1" "1" "1"# クロネッカー積を使うと以下の通り
u %@% t(v)## c1-1 c2-2 c3-3
## r1-1 "u1*v1" "u1*v2" "u1*v3"
## r1-2 "u2*v1" "u2*v2" "u2*v3"
## r1-3 "u3*v1" "u3*v2" "u3*v3"
## attr(,"SMdim")
## attr(,"SMdim")$nrows
## [1] 1 1 1
##
## attr(,"SMdim")$ncols
## [1] 3この行列をデータ行列の左から掛けたものは、その列平均が全て 0 となる。
# set output line width
options(width=500)
# 中心化行列: J = I_n - (1/n) * one_n %*% t(one_n)
( J3=I3 - "(1/3)" * one3 %*% t(one3) )## [,1] [,2] [,3]
## [1,] "1-1/3" "-1/3" "-1/3"
## [2,] "-1/3" "1-1/3" "-1/3"
## [3,] "-1/3" "-1/3" "1-1/3"mat2sum 関数を使うと、ベクトルと行列の演算の結果として得られるベクトルの要素を 総和記号で表現することが出来る。 たとえば、 \boldsymbol A + \boldsymbol B の i,j 要素は a_{i,j} + b_{i,j} であり \boldsymbol v' \boldsymbol u は \sum_{s=1}^{r} v_s u_s であるが これらは mat2sum 関数では以下の様に表現される。 なお、simple オプションを省略した場合には添字の範囲は示されない。
# set output line width
options(width=500)
# A + B の総和記号による表現:
mat2sum( "A + B", simple=0 )## [1] "A[i,j] + B[i,j]"# 内積の総和記号による表現
mat2sum( "t(v) %*% u", simple=0 )## [1] "s(v[s1]*u[s1],{s1,1,nrow(v)})"# 外積
mat2sum( "v %*% t(u)" ) ## [1] "v[i]*u[j]"# 複雑な式
mat2sum("A+B-diag(D)+u%*%t(v)")## [1] "A[i,j] + B[i,j] - diag(D)[i,j] + u[i]*v[j]"r x c の行列の右から c x 1 の縦ベクトルを掛けることが出来る。 このとき、A 行列を r 個の 1 x c の横ベクトルを縦に並べたものと考えれば、 この積は r 個の内積を縦に並べた r x 1 の縦ベクトルとなる。
# set output line width
options(width=500)
# 行列とベクトルの積
( Av=A %*% v )## [,1]
## r1 "a11*v1+a12*v2+a13*v3"
## r2 "a21*v1+a22*v2+a23*v3"
## r3 "a31*v1+a32*v2+a33*v3"
## r4 "a41*v1+a42*v2+a43*v3"# A %*% v は A 行列の各行と v との内積を縦に並べたもの
Av1=rbind(A[1,,drop=0]%*%v, A[2,,drop=0]%*%v, A[3,,drop=0]%*%v,A[4,,drop=0]%*%v)
Av-Av1## [,1]
## r1 "0"
## r2 "0"
## r3 "0"
## r4 "0"同様に、1 x r の横ベクトルの右から r x c の行列を掛けることが出来る。 このとき、A 行列を c 個の r x 1 の縦ベクトルを横に並べたものと考えれば、 この積は c 個の内積を横に並べた c x 1 の横ベクトルとなる。
# set output line width
options(width=500)
# 行列とベクトルの積
( tvA=t(v4) %*% A )## c1 c2 c3
## [1,] "v1*a11+v2*a21+v3*a31+v4*a41" "v1*a12+v2*a22+v3*a32+v4*a42" "v1*a13+v2*a23+v3*a33+v4*a43"# A %*% v は A 行列の各行と v との内積を縦に並べたもの
tvA1=cbind(t(v4)%*%A[,1], t(v4)%*%A[,2], t(v4)%*%A[,3])
tvA-tvA1## c1 c2 c3
## [1,] "0" "0" "0"mat2sum 関数を使うと、ベクトルと行列の積の結果として得られるベクトルの要素を 総和記号で表現することが出来る。 たとえば、 \boldsymbol A \; \boldsymbol v の第 i 要素は \sum_{s=1}^{c} a_{i,s} v_{s} であり \boldsymbol v' \boldsymbol A の第 j 要素は \sum_{s=1}^{r} v_s a_{s,j} であるが これらは mat2sum 関数では以下の様に表現される。 なお、simple オプションを省略した場合には添字の範囲は示されない。
# set output line width
options(width=500)
# A %*% v の総和記号による表現: \sum_{s=1}^{r} a_{i,s} v_{s}
mat2sum( "A %*% v", simple=0 )## [1] "s(A[i,s1]*v[s1],{s1,1,ncol(A)})"# t(v) %*% A の総和記号による表現: \sum_{s=1}^{r} v[s] a_{s,j}
mat2sum( "t(v) %*% A", simple=0 )## [1] "s(v[s1]*A[s1,j],{s1,1,nrow(v)})"mat2sum( "t(v) %*% A" ) # simple version## [1] "s(v[s1]*A[s1,j],{s1})"演算結果が複雑な場合には、 得られた結果に Expand 関数や Simplify() 関数を作用させると、 結果が簡略化される場合がある。
また、結果が数値であることが明らかな場合、 Eval 関数を作用させると簡略化される。
# set output line width
options(width=500)
# 単位行列との積
I3 %*% v## [,1]
## [1,] "v1"
## [2,] "v2"
## [3,] "v3"t(v) %*% I3## [,1] [,2] [,3]
## [1,] "v1" "v2" "v3"# 行の和 row sum
A %*% one3## [,1]
## r1 "a11+a12+a13"
## r2 "a21+a22+a23"
## r3 "a31+a32+a33"
## r4 "a41+a42+a43"# R style
apply( A, 1, Sum )## r1 r2 r3 r4
## "a11+a12+a13" "a21+a22+a23" "a31+a32+a33" "a41+a42+a43"#
# 列の和 row sum
t(one4) %*% A## c1 c2 c3
## [1,] "a11+a21+a31+a41" "a12+a22+a32+a42" "a13+a23+a33+a43"# 列の平均
(1/4) * t(one4) %*% A## c1 c2 c3
## [1,] "0.25*(a11+a21+a31+a41)" "0.25*(a12+a22+a32+a42)" "0.25*(a13+a23+a33+a43)"# R style
apply( A, 2, function(x) Sum(x)*(1/4) )## c1 c2 c3
## "(a11+a21+a31+a41)*0.25" "(a12+a22+a32+a42)*0.25" "(a13+a23+a33+a43)*0.25"# 中心化行列の行和と列和
t(one3)%*%J3 |> print() |> Expand() |> print() |> Eval() |> round(2)## [,1] [,2] [,3]
## [1,] "(1-1/3)-2*1/3" "-2*1/3+(1-1/3)" "-2*1/3+(1-1/3)"
## [,1] [,2] [,3]
## [1,] "1-1/3-2/3" "-2/3+1-1/3" "-2/3+1-1/3"## [,1] [,2] [,3]
## [1,] 0 0 0apply(J3,1,Sum) |> print() |> Eval() |> round(2)## [1] "1-2/3-1/3" "-2/3+1-1/3" "-2/3+1-1/3"## [1] 0 0 0J3 %*% one3 |> print() |> Eval() |> round(2)## [,1]
## [1,] "(1-1/3)-2*1/3"
## [2,] "-2*1/3+(1-1/3)"
## [3,] "-2*1/3+(1-1/3)"## [,1]
## [1,] 0
## [2,] 0
## [3,] 0apply(J3,2,Sum) |> print() |> Eval() |> round(2)## [1] "1-2/3-1/3" "-2/3+1-1/3" "-2/3+1-1/3"## [1] 0 0 0# 対角行列とベクトルの積
D3 %*% v## [,1]
## [1,] "d11*v1"
## [2,] "d22*v2"
## [3,] "d33*v3"t(v) %*% D3## [,1] [,2] [,3]
## [1,] "v1*d11" "v2*d22" "v3*d33"# D3 %*% v の第 i 要素
mat2sum("D3%*%v")## [1] "s(D3[i,s1]*v[s1],{s1})"# t(v) %*% D3 の第 j 要素
mat2sum("t(v)%*%D3")## [1] "s(v[s1]*D3[s1,j],{s1})"# 対角行列
mat2sum("diag(D)%*%v")## [1] "D[i,i]*v[i]"mat2sum("t(v)%*%diag(D)")## [1] "v[j]*D[j,j]"正方行列 \boldsymbol S とベクトル \boldsymbol v が与えられた時に \boldsymbol v' \boldsymbol S \boldsymbol v を \boldsymbol S を用いた \boldsymbol v の２次形式と呼ぶが その値はスカラーで \sum_{k1=1}^r \sum_{k2=1}^c \; v_{k1} s_{k1,k2} v_{k2} となる。
# set output line width
options(width=500)
t(v) %*% S %*% v |> Expand()## [,1]
## [1,] "v1*s11*v1+v2*s21*v1+v3*s31*v1+v1*s21*v2+v2*s22*v2+v3*s32*v2+v1*s31*v3+v2*s32*v3+v3*s33*v3"mat2sum("t(v) %*% S %*% v", root="k",simple=0)## [1] "s(v[k1]*s(S[k1,k2]*v[k2],{k2,1,ncol(S)}),{k1,1,nrow(v)})"a x b の行列 A と c x d の行列 B が与えられた場合、 もし b=c であれば A の右から B を掛けることが出来る。 また、もし d = a であれば A の左から B を掛けることが出来る。
# set output line width
options(width=500)
# 行列と行列の積 1: C = A %*% E
# 右から掛かる行列を列ベクトルの集まりと見做す。
E=cbind(u,v)
Print(u,v,E)## c1 c1 c1 c2
## u = r1 u1 , v = r1 v1 , E = r1 u1 v1
## r2 u2 r2 v2 r2 u2 v2
## r3 u3 r3 v3 r3 u3 v3# 行列とそれぞれのベクトルの積を横に並べる。
options(width=500)
c1=A%*%u; c2=A%*%v
C0=cbind(c1,c2)
Print(c1,c2,C0)## c1 c1 c1 c2
## c1 = r1 a11*u1+a12*u2+a13*u3 , c2 = r1 a11*v1+a12*v2+a13*v3 , C0 = r1 a11*u1+a12*u2+a13*u3 a11*v1+a12*v2+a13*v3
## r2 a21*u1+a22*u2+a23*u3 r2 a21*v1+a22*v2+a23*v3 r2 a21*u1+a22*u2+a23*u3 a21*v1+a22*v2+a23*v3
## r3 a31*u1+a32*u2+a33*u3 r3 a31*v1+a32*v2+a33*v3 r3 a31*u1+a32*u2+a33*u3 a31*v1+a32*v2+a33*v3
## r4 a41*u1+a42*u2+a43*u3 r4 a41*v1+a42*v2+a43*v3 r4 a41*u1+a42*u2+a43*u3 a41*v1+a42*v2+a43*v3# 上記の演算を行列の積で書くと
C=A %*% E
Print(C,"\n", C0)## C =
## c1 c2
## r1 a11*u1+a12*u2+a13*u3 a11*v1+a12*v2+a13*v3
## r2 a21*u1+a22*u2+a23*u3 a21*v1+a22*v2+a23*v3
## r3 a31*u1+a32*u2+a33*u3 a31*v1+a32*v2+a33*v3
## r4 a41*u1+a42*u2+a43*u3 a41*v1+a42*v2+a43*v3
##
## C0 =
## c1 c2
## r1 a11*u1+a12*u2+a13*u3 a11*v1+a12*v2+a13*v3
## r2 a21*u1+a22*u2+a23*u3 a21*v1+a22*v2+a23*v3
## r3 a31*u1+a32*u2+a33*u3 a31*v1+a32*v2+a33*v3
## r4 a41*u1+a42*u2+a43*u3 a41*v1+a42*v2+a43*v3C-C0## [,1] [,2]
## r1 "0" "0"
## r2 "0" "0"
## r3 "0" "0"
## r4 "0" "0"# set output line width
options(width=500)
# 行列と行列の積 2: R = F %*% A
# 左から掛かる行列を行ベクトルの集まりと見做す。
F=rbind(t(u4), t(v4))
Print(t(u4),t(v4),F)## c1 c2 c3 c4 c1 c2 c3 c4 c1 c2 c3 c4
## t(u4) = r1 u1 u2 u3 u4 , t(v4) = r1 v1 v2 v3 v4 , F = r1 u1 u2 u3 u4
## r2 v1 v2 v3 v4# 行列とそれぞれのベクトルの転置の積（行ベクトル）を縦に並べる。
r1=t(u4)%*%A; r2=t(v4)%*%A
R0=rbind(r1,r2)
Print(r1, "\n", r2,"\n",R0)## r1 =
## c1 c2 c3
## r1 u1*a11+u2*a21+u3*a31+u4*a41 u1*a12+u2*a22+u3*a32+u4*a42 u1*a13+u2*a23+u3*a33+u4*a43
##
## r2 =
## c1 c2 c3
## r1 v1*a11+v2*a21+v3*a31+v4*a41 v1*a12+v2*a22+v3*a32+v4*a42 v1*a13+v2*a23+v3*a33+v4*a43
##
## R0 =
## c1 c2 c3
## r1 u1*a11+u2*a21+u3*a31+u4*a41 u1*a12+u2*a22+u3*a32+u4*a42 u1*a13+u2*a23+u3*a33+u4*a43
## r2 v1*a11+v2*a21+v3*a31+v4*a41 v1*a12+v2*a22+v3*a32+v4*a42 v1*a13+v2*a23+v3*a33+v4*a43# 上記の演算を行列の積で書くと
R=F %*% A
Print(R, "\n", R0)## R =
## c1 c2 c3
## r1 u1*a11+u2*a21+u3*a31+u4*a41 u1*a12+u2*a22+u3*a32+u4*a42 u1*a13+u2*a23+u3*a33+u4*a43
## r2 v1*a11+v2*a21+v3*a31+v4*a41 v1*a12+v2*a22+v3*a32+v4*a42 v1*a13+v2*a23+v3*a33+v4*a43
##
## R0 =
## c1 c2 c3
## r1 u1*a11+u2*a21+u3*a31+u4*a41 u1*a12+u2*a22+u3*a32+u4*a42 u1*a13+u2*a23+u3*a33+u4*a43
## r2 v1*a11+v2*a21+v3*a31+v4*a41 v1*a12+v2*a22+v3*a32+v4*a42 v1*a13+v2*a23+v3*a33+v4*a43R-R0## c1 c2 c3
## [1,] "0" "0" "0"
## [2,] "0" "0" "0"# set output line width
options(width=500)
# A %*% B の i-j 要素は
mat2sum("A%*%B", simple=0)## [1] "s(A[i,s1]*B[s1,j],{s1,1,ncol(A)})"# set output line width
options(width=500)
# 単位行列との積
A %*% I3## [,1] [,2] [,3]
## r1 "a11" "a12" "a13"
## r2 "a21" "a22" "a23"
## r3 "a31" "a32" "a33"
## r4 "a41" "a42" "a43"I4 %*% A## c1 c2 c3
## [1,] "a11" "a12" "a13"
## [2,] "a21" "a22" "a23"
## [3,] "a31" "a32" "a33"
## [4,] "a41" "a42" "a43"# set output line width
options(width=500)
# 対角行列との積：列や行の定数倍
A %*% D3 # A の各列を対応する D3 の対角要素倍したもの## [,1] [,2] [,3]
## r1 "a11*d11" "a12*d22" "a13*d33"
## r2 "a21*d11" "a22*d22" "a23*d33"
## r3 "a31*d11" "a32*d22" "a33*d33"
## r4 "a41*d11" "a42*d22" "a43*d33"D4 %*% A # A の各行を対応する D4 の対角要素倍したもの## c1 c2 c3
## [1,] "d11*a11" "d11*a12" "d11*a13"
## [2,] "d22*a21" "d22*a22" "d22*a23"
## [3,] "d33*a31" "d33*a32" "d33*a33"
## [4,] "d44*a41" "d44*a42" "d44*a43"# 総和記号による表現
mat2sum("A%*%diag(D)")## [1] "A[i,j]*D[j,j]"mat2sum("diag(D)%*%A")## [1] "D[i,i]*A[i,j]"mat2sum("diag(D1)%*%A%*%diag(D2)")## [1] "D1[i,i]*A[i,j]*D2[j,j]"以下にいくつかの行列演算の規則を列挙する。 結果が複雑な場合、 目視で二つの結果が等しいことを確認することが難しい場合があるが、 そのような場合には、記号的行列と同じサイズの数値的行列を定義し、 得られた結果をそれらの数値的行列を用いて評価することにより確認することが出来る。 そのためには Eval 関数を用いる。
Eval 関数の利用法は
Eval( expr, pats1=patn1, pats2=patn2, … )
であり、expr は評価したい記号演算式、 pats1 は expr に含まれる記号的行列やベクトルの名称、 patn1 はそれらと同じサイズの数値的行列やベクトルである。 記号的行列やベクトルが複数ある場合には pats2, patn2 等を指定する。
# set output line width
options(width=500)
# 結合法則 associative law
# (A+B) + C = A + (B+C) = A+B+C
X+Y+Z## c1 c2
## r1 "x11+y11+z11" "x12+y12+z12"
## r2 "x21+y21+z21" "x22+y22+z22"
## r3 "x31+y31+z31" "x32+y32+z32"XY=(X+Y)
XY+Z## c1 c2
## r1 "x11+y11+z11" "x12+y12+z12"
## r2 "x21+y21+z21" "x22+y22+z22"
## r3 "x31+y31+z31" "x32+y32+z32"YZ=Y+Z
X+YZ## c1 c2
## r1 "x11+y11+z11" "x12+y12+z12"
## r2 "x21+y21+z21" "x22+y22+z22"
## r3 "x31+y31+z31" "x32+y32+z32"# (A %*% B) %*% C = A %*% (B %*% C) = A %*% B %*% C
res1=t(X) %*% Y %*% t(Z)
res2=t(X) %*% Y
res21=res2 %*% t(Z)
res1-res21## r1 r2 r3
## c1 "0" "0" "0"
## c2 "0" "0" "0"res3=Y %*% t(Z)
res31=t(X) %*% res3
diff=res1-res31
diff |> Expand() |> Expand()## r1 r2 r3
## c1 "0" "0" "0"
## c2 "0" "0" "0"# set output line width
options(width=500)
# 分配法則 distributional law Expansion or Expanding
# A %*% (B + C) = A %*% B + A %*% C
# (A + B) %*% C = A %*% C + B %*% C
res1=t(X) %*% (Y+Z)
res2=t(X) %*% Y + t(X) %*% Z
(res1-res2) |> Expand()## c1 c2
## c1 "0" "0"
## c2 "0" "0"# 数値を代入
Eval(res1-res2, X=demomat(3,2), Y=demomat(3,2), Z=demomat(3,2) )## [,1] [,2]
## [1,] 0 0
## [2,] 0 0# set output line width
options(width=500)
# #交換法則 commutative law
# A + B = B + A
X+Y## c1 c2
## r1 "x11+y11" "x12+y12"
## r2 "x21+y21" "x22+y22"
## r3 "x31+y31" "x32+y32"Y+X## c1 c2
## r1 "y11+x11" "y12+x12"
## r2 "y21+x21" "y22+x22"
## r3 "y31+x31" "y32+x32"diff=X+Y-(Y+X)
Expand(diff)## c1 c2
## r1 "0" "0"
## r2 "0" "0"
## r3 "0" "0"# set output line width
options(width=500)
# 行列とベクトルの積の転置
( tAv = t(A%*%v) )## r1 r2 r3 r4
## [1,] "a11*v1+a12*v2+a13*v3" "a21*v1+a22*v2+a23*v3" "a31*v1+a32*v2+a33*v3" "a41*v1+a42*v2+a43*v3"( tvtA = t(v)%*%t(A) )## r1 r2 r3 r4
## [1,] "v1*a11+v2*a12+v3*a13" "v1*a21+v2*a22+v3*a23" "v1*a31+v2*a32+v3*a33" "v1*a41+v2*a42+v3*a43"# set output line width
options(width=500)
# 行列の積の転置
# 行列を列ベクトルの集合と見た場合、
# 転置結果は行ベクトルを縦に並べたものとなる。
options(width=500)
E = cbind(u,v)
# C = A %*% E = A %*% cbind(u,v) = cbind( A%*%u, A%*%v )
c1=A%*%u; c2=A%*%v
C=cbind(c1,c2)
tC=t(C)
tC1=t( cbind( c1, c2 ) )
tC2=rbind( t(c1), t(c2) )
Print(tC, "\n", tC1, "\n", tC2)## tC =
## r1 r2 r3 r4
## r1 a11*u1+a12*u2+a13*u3 a21*u1+a22*u2+a23*u3 a31*u1+a32*u2+a33*u3 a41*u1+a42*u2+a43*u3
## r2 a11*v1+a12*v2+a13*v3 a21*v1+a22*v2+a23*v3 a31*v1+a32*v2+a33*v3 a41*v1+a42*v2+a43*v3
##
## tC1 =
## r1 r2 r3 r4
## r1 a11*u1+a12*u2+a13*u3 a21*u1+a22*u2+a23*u3 a31*u1+a32*u2+a33*u3 a41*u1+a42*u2+a43*u3
## r2 a11*v1+a12*v2+a13*v3 a21*v1+a22*v2+a23*v3 a31*v1+a32*v2+a33*v3 a41*v1+a42*v2+a43*v3
##
## tC2 =
## r1 r2 r3 r4
## r1 a11*u1+a12*u2+a13*u3 a21*u1+a22*u2+a23*u3 a31*u1+a32*u2+a33*u3 a41*u1+a42*u2+a43*u3
## r2 a11*v1+a12*v2+a13*v3 a21*v1+a22*v2+a23*v3 a31*v1+a32*v2+a33*v3 a41*v1+a42*v2+a43*v3# 行列を行ベクトルの集合と見た場合、転置結果は列ベクトルを横に並べたものとなる。
# R0=rbind( r1=t(u4)%*%A, r2=t(v4)%*%A )
# 省略
# まとめると: それぞれを転置したものを逆順に掛ける。
C=A %*% E
tC=t( A %*% E )
tC1=t(E) %*% t(A)
Print(tC, "\n", tC1)## tC =
## r1 r2 r3 r4
## r1 a11*u1+a12*u2+a13*u3 a21*u1+a22*u2+a23*u3 a31*u1+a32*u2+a33*u3 a41*u1+a42*u2+a43*u3
## r2 a11*v1+a12*v2+a13*v3 a21*v1+a22*v2+a23*v3 a31*v1+a32*v2+a33*v3 a41*v1+a42*v2+a43*v3
##
## tC1 =
## r1 r2 r3 r4
## r1 u1*a11+u2*a12+u3*a13 u1*a21+u2*a22+u3*a23 u1*a31+u2*a32+u3*a33 u1*a41+u2*a42+u3*a43
## r2 v1*a11+v2*a12+v3*a13 v1*a21+v2*a22+v3*a23 v1*a31+v2*a32+v3*a33 v1*a41+v2*a42+v3*a43# 数値的に検算
tCn=Eval( tC, A=demomat(4,3), u=demomat(3,1), v=demomat(3,1) )
tC1n=Eval( tC1, A=demomat(4,3), u=demomat(3,1), v=demomat(3,1) )
Print(tCn,tC1n)## c1 c2 c3 c4 c1 c2 c3 c4
## tCn = r1 776 1406 2036 2666 , tC1n = r1 776 1406 2036 2666
## r2 776 1406 2036 2666 r2 776 1406 2036 2666Eval(tC-tC1, A=demomat(4,3), u=demomat(3,1), v=demomat(3,1) )## [,1] [,2] [,3] [,4]
## [1,] 0 0 0 0
## [2,] 0 0 0 0# 自分自身の転置との積
t(A) %*% A## c1 c2 c3
## c1 "a11*a11+a21*a21+a31*a31+a41*a41" "a11*a12+a21*a22+a31*a32+a41*a42" "a11*a13+a21*a23+a31*a33+a41*a43"
## c2 "a12*a11+a22*a21+a32*a31+a42*a41" "a12*a12+a22*a22+a32*a32+a42*a42" "a12*a13+a22*a23+a32*a33+a42*a43"
## c3 "a13*a11+a23*a21+a33*a31+a43*a41" "a13*a12+a23*a22+a33*a32+a43*a42" "a13*a13+a23*a23+a33*a33+a43*a43"A %*% t(A)## r1 r2 r3 r4
## r1 "a11*a11+a12*a12+a13*a13" "a11*a21+a12*a22+a13*a23" "a11*a31+a12*a32+a13*a33" "a11*a41+a12*a42+a13*a43"
## r2 "a21*a11+a22*a12+a23*a13" "a21*a21+a22*a22+a23*a23" "a21*a31+a22*a32+a23*a33" "a21*a41+a22*a42+a23*a43"
## r3 "a31*a11+a32*a12+a33*a13" "a31*a21+a32*a22+a33*a23" "a31*a31+a32*a32+a33*a33" "a31*a41+a32*a42+a33*a43"
## r4 "a41*a11+a42*a12+a43*a13" "a41*a21+a42*a22+a43*a23" "a41*a31+a42*a32+a43*a33" "a41*a41+a42*a42+a43*a43"正方行列の対角要素の和を行列のトレースと呼ぶが、それはスカラー値である。
# set output line width
options(width=500)
# トレース： 正方行列の対角要素の和
tr(S)## [1] "s11+s22+s33"Sum(Diag(S))## [,1]
## [1,] "s11+s22+s33"# 行列の対応する要素の積の和
tr( t(A) %*% B )## [1] "a11*b11+a21*b21+a31*b31+a41*b41+a12*b12+a22*b22+a32*b32+a42*b42+a13*b13+a23*b23+a33*b33+a43*b43"# 行列の要素の２乗和
ssq1=tr( t(A) %*% A )
ssq2=tr( A %*% t(A) )
Print(ssq1,ssq2)##
## ssq1 = a11*a11+a21*a21+a31*a31+a41*a41+a12*a12+a22*a22+a32*a32+a42*a42+a13*a13+a23*a23+a33*a33+a43*a43 , ssq2 = a11*a11+a12*a12+a13*a13+a21*a21+a22*a22+a23*a23+a31*a31+a32*a32+a33*a33+a41*a41+a42*a42+a43*a43# 数値的に検算
Eval(ssq1-ssq2, A=demomat(4,3))## [,1]
## [1,] 0# set output line width
options(width=500)
# 線形性
# tr(A+B)=tr(A)+tr(B), d tr(A)=tr(d*A)
A2=demomat(2,2,root="a"); B2=demomat(2,2,root="b"); C2=demomat(2,2,root="c")
tr(A2+B2)## [1] "a11+b11+a22+b22"tr(A2)+tr(B2)## [,1]
## [1,] "a11+a22+b11+b22"tr("d"*A2)## [1] "d*a11+d*a22""d"*tr(A2)## [1] "d*(a11+a22)"# 積のトレース:循環性
# tr(AB) = tr(BA)
tr(A2%*%B2)## [1] "a11*b11+a12*b21+a21*b12+a22*b22"tr(B2%*%A2)## [1] "b11*a11+b12*a21+b21*a12+b22*a22"# tr(ABC) = tr(CAB) = tr(CAB)
tr1=tr(A2%*%B2%*%C2) |> Expand() |> print()## [1] "a11*b11*c11+a12*b21*c11+a11*b12*c21+a12*b22*c21+a21*b11*c12+a22*b21*c12+a21*b12*c22+a22*b22*c22"tr2=tr(B2%*%C2%*%A2) |> Expand() |> print()## [1] "b11*c11*a11+b12*c21*a11+b11*c12*a21+b12*c22*a21+b21*c11*a12+b22*c21*a12+b21*c12*a22+b22*c22*a22"tr3=tr(C2%*%A2%*%B2) |> Expand() |> print()## [1] "c11*a11*b11+c12*a21*b11+c11*a12*b21+c12*a22*b21+c21*a11*b12+c22*a21*b12+c21*a12*b22+c22*a22*b22"# 数値的検証
seed=1701
A2n=matrix(sample(0:9,2*2, replace=1),2,2)
B2n=matrix(sample(0:9,2*2, replace=1),2,2)
C2n=matrix(sample(0:9,2*2, replace=1),2,2)
Eval(tr1-tr2, A2=A2n, B2=B2n, C2=C2n)## [,1]
## [1,] 0Eval(tr1-tr3, A2=A2n, B2=B2n, C2=C2n)## [,1]
## [1,] 0# 転置のトレース
tr(t(A2))## [1] "a11+a22"tr(A2)## [1] "a11+a22"# 転置と循環性、
tr1=tr(t(A2)%*%B2) |> Expand() |> print()## [1] "a11*b11+a21*b21+a12*b12+a22*b22"tr2=tr(t(B2)%*%A2) |> Expand() |> print()## [1] "b11*a11+b21*a21+b12*a12+b22*a22"tr3=tr(B2%*%t(A2)) |> Expand() |> print()## [1] "b11*a11+b12*a12+b21*a21+b22*a22"tr4=tr(A2%*%t(B2)) |> Expand() |> print()## [1] "a11*b11+a12*b12+a21*b21+a22*b22"Eval(tr1-tr2, A2=A2n, B2=B2n, C2=C2n)## [,1]
## [1,] 0Eval(tr1-tr3, A2=A2n, B2=B2n, C2=C2n)## [,1]
## [1,] 0# set output line width
options(width=500)
mat2sum("A%*%B%*%C%*%D%*%E")## [1] "s(A[i,s1]*s(B[s1,s2]*s(C[s2,s3]*s(D[s3,s4]*E[s4,j],{s4}),{s3}),{s2}),{s1})"mat2sum("A%*%B%*%C%*%Diag(D)%*%E%*%f")## [1] "s(A[i,s1]*s(B[s1,s2]*s(C[s2,s3]*D[s3,s3]*s(E[s3,s4]*f[s4],{s4}),{s3}),{s2}),{s1})"mat2sum("t(f)%*%A%*%t(B)%*%C%*%Diag(D)%*%t(E)")## [1] "s(f[s1]*s(A[s1,s2]*s(B[s3,s2]*s(C[s3,s4]*D[s4,s4]*E[j,s4],{s4}),{s3}),{s2}),{s1})"# set output line width
options(width=500)
options(width=500)
# 三角行列: 掛けても足しても三角のまま。
L1=demomat(3,3,root="L1_",shape="lowert")
L2=demomat(3,3,root="L2_",shape="lowert")
Print(L1,L2)## c1 c2 c3 c1 c2 c3
## L1 = r1 L1_11 0 0 , L2 = r1 L2_11 0 0
## r2 L1_21 L1_22 0 r2 L2_21 L2_22 0
## r3 L1_31 L1_32 L1_33 r3 L2_31 L2_32 L2_33L1 %*% L2## c1 c2 c3
## r1 "L1_11*L2_11" "0" "0"
## r2 "L1_21*L2_11+L1_22*L2_21" "L1_22*L2_22" "0"
## r3 "L1_31*L2_11+L1_32*L2_21+L1_33*L2_31" "L1_32*L2_22+L1_33*L2_32" "L1_33*L2_33"自分自身の転置を左から掛けたときに単位行列となる正方行列を直交行列と呼ぶ。
# set output line width
options(width=500)
#
# ハウスホルダー変換
# どのような 3 x 1 の縦ベクトル v を用いても以下の行列は直交行列となる。
den=paste("2/(",Sum(v*v),")", sep="")
O3=I3-v%*%t(v)*den |> print()## [1] "2/(v1*v1+v2*v2+v3*v3)"# 数値的チェック
O3n=Eval(O3, v1=1, v2=2, v3=3)
Print(t(O3n)%*%O3n, O3n%*%t(O3n), fmt="6.2")## c1 c2 c3 c1 c2 c3
## t(O3n) %*% O3n = r1 1.00 -0.00 0.00 , O3n %*% t(O3n) = r1 1.00 -0.00 0.00
## r2 -0.00 1.00 0.00 r2 -0.00 1.00 0.00
## r3 0.00 0.00 1.00 r3 0.00 0.00 1.00O2=matrix(c("cos(theta)", "-sin(theta)"
, "sin(theta)", "cos(theta)"), 2,2, byrow=1) |> print()## [,1] [,2]
## [1,] "cos(theta)" "-sin(theta)"
## [2,] "sin(theta)" "cos(theta)"tO2O2=(t(O2)%*%O2)
O2tO2=O2%*%t(O2)
tO2O2n=Eval(tO2O2, theta=5, check=0)
O2tO2n=Eval(O2tO2, theta=5, check=0)
Print(tO2O2n, O2tO2n, fmt="6.2")## c1 c2 c1 c2
## tO2O2n = r1 1.00 0.00 , O2tO2n = r1 1.00 0.00
## r2 0.00 1.00 r2 0.00 1.00# 列が直交する行列: semi-orthogonal matrix
Z=matrix(c(1,1,0,0,2,-1),3,2)
Z=matrix(c("z1","-z1","0", "z1","z1", "z2"), 3,2)
dZ=matrix("0",2,2)
dZ[1,1]="1/sqrt(2*z1^2)"
dZ[2,2]="1/sqrt(2*z1^2+z2^2)"
Z=Z%*%dZ
Print(Z)## Z =
## c1 c2
## r1 z1/sqrt(2*z1^2) z1/sqrt(2*z1^2+z2^2)
## r2 -z1/sqrt(2*z1^2) z1/sqrt(2*z1^2+z2^2)
## r3 0 z2/sqrt(2*z1^2+z2^2)# 列は直交する
tZZ=t(Z)%*%Z
Eval(tZZ, z1=1, z2=2, check=0)## [,1] [,2]
## [1,] 1 0
## [2,] 0 1# 行は直交しない
ZtZ=Z%*%t(Z)
Eval(ZtZ, z1=1, z2=2, check=0)## [,1] [,2] [,3]
## [1,] 0.6666667 -0.3333333 0.3333333
## [2,] -0.3333333 0.6666667 0.3333333
## [3,] 0.3333333 0.3333333 0.6666667# set output line width
options(width=500)
# 行と列の入れ替え：右から掛けると列の、左から掛けると行の入れ替え
P3=pmat(3, 1,2)
P4=pmat(4, 2,4)
Print(P3,P4)## c1 c2 c3 c1 c2 c3 c4
## P3 = r1 0 1 0 , P4 = r1 1 0 0 0
## r2 1 0 0 r2 0 0 0 1
## r3 0 0 1 r3 0 0 1 0
## r4 0 1 0 0P3 %*% v## [,1]
## [1,] "v2"
## [2,] "v1"
## [3,] "v3"A %*% P3## [,1] [,2] [,3]
## r1 "a12" "a11" "a13"
## r2 "a22" "a21" "a23"
## r3 "a32" "a31" "a33"
## r4 "a42" "a41" "a43"P4 %*% v4## [,1]
## [1,] "v1"
## [2,] "v4"
## [3,] "v3"
## [4,] "v2"P4 %*% A## c1 c2 c3
## [1,] "a11" "a12" "a13"
## [2,] "a41" "a42" "a43"
## [3,] "a31" "a32" "a33"
## [4,] "a21" "a22" "a23"P4 %*% A %*% P3## [,1] [,2] [,3]
## [1,] "a12" "a11" "a13"
## [2,] "a42" "a41" "a43"
## [3,] "a32" "a31" "a33"
## [4,] "a22" "a21" "a23"# 入れ替え行列は直交行列である。
Print(t(P3)%*%P3, P3%*%t(P3))## c1 c2 c3 c1 c2 c3
## t(P3) %*% P3 = r1 1 0 0 , P3 %*% t(P3) = r1 1 0 0
## r2 0 1 0 r2 0 1 0
## r3 0 0 1 r3 0 0 1# set output line width
options(width=500)
# 列方向の中心化 column centering
# 中心化行列: J = I_n - (1/n) * one_n %*% t(one_n)
JX=J3 %*% X |> print() |> Expand() |> print()## c1 c2
## [1,] "(1-1/3)*x11-1/3*x21-1/3*x31" "(1-1/3)*x12-1/3*x22-1/3*x32"
## [2,] "-1/3*x11+(1-1/3)*x21-1/3*x31" "-1/3*x12+(1-1/3)*x22-1/3*x32"
## [3,] "-1/3*x11-1/3*x21+(1-1/3)*x31" "-1/3*x12-1/3*x22+(1-1/3)*x32"
## c1 c2
## [1,] "x11-1/3*x11-1/3*x21-1/3*x31" "x12-1/3*x12-1/3*x22-1/3*x32"
## [2,] "-1/3*x11+x21-1/3*x21-1/3*x31" "-1/3*x12+x22-1/3*x22-1/3*x32"
## [3,] "-1/3*x11-1/3*x21+x31-1/3*x31" "-1/3*x12-1/3*x22+x32-1/3*x32"# 数値的に確かめる。
Eval(JX, X=demomat(3,2)) ## [,1] [,2]
## [1,] -10 -1.000000e+01
## [2,] 0 1.776357e-15
## [3,] 10 1.000000e+01seed0=1701
Eval(JX, X=matrix(rnorm(3*2),3,2)) |> round(2)## [,1] [,2]
## [1,] 1.04 1.23
## [2,] -0.02 -0.74
## [3,] -1.02 -0.49# 列平均を求める。
t(one3) %*% JX |> print() |> Expand()## c1 c2
## [1,] "(x11-1/3*x11-1/3*x21-1/3*x31)+(-1/3*x11+x21-1/3*x21-1/3*x31)+(-1/3*x11-1/3*x21+x31-1/3*x31)" "(x12-1/3*x12-1/3*x22-1/3*x32)+(-1/3*x12+x22-1/3*x22-1/3*x32)+(-1/3*x12-1/3*x22+x32-1/3*x32)"## c1 c2
## [1,] "x11-2*1/3*x11-2*1/3*x21-2*1/3*x31+x21-1/3*x11-1/3*x21+x31-1/3*x31" "x12-2*1/3*x12-2*1/3*x22-2*1/3*x32+x22-1/3*x12-1/3*x22+x32-1/3*x32"apply(JX,2,Sum)## c1 c2
## "x11-2/3*x11-2/3*x21-2/3*x31+x21-1/3*x11-1/3*x21+x31-1/3*x31" "x12-2/3*x12-2/3*x22-2/3*x32+x22-1/3*x12-1/3*x22+x32-1/3*x32"# 数値的に確かめる。
Eval(t(one3) %*% JX, X=demomat(3,2)) |> round(2)## [,1] [,2]
## [1,] 0 0Eval(t(one3) %*% JX, X=matrix(rnorm(3*2),3,2)) |> round(2)## [,1] [,2]
## [1,] 0 0Eval(apply(JX,2,Sum), X=matrix(rnorm(3*2),3,2)) |> round(2)## [1] 0 0# set output line width
options(width=500)
# 分散共分散行列
# var_j = (1/n) \sum_{i=1}^n (a_{ij}-m_j)^2
# = (1/n) \sum_{i=1}^n a_{ij}^2 - m_j^2
# cov_{jk} = (1/n) \sum_{i=1}^n (a_{ij}-m_j)*(a_{ik}-m_k)
#
DispX=(1/nrow(X))*t(JX) %*% JX
# 分散と共分散
DispX[1,1]## [1] "0.333333333333333*((x11-1/3*x11-1/3*x21-1/3*x31)*(x11-1/3*x11-1/3*x21-1/3*x31)+(-1/3*x11+x21-1/3*x21-1/3*x31)*(-1/3*x11+x21-1/3*x21-1/3*x31)+(-1/3*x11-1/3*x21+x31-1/3*x31)*(-1/3*x11-1/3*x21+x31-1/3*x31))"DispX[2,2]## [1] "0.333333333333333*((x12-1/3*x12-1/3*x22-1/3*x32)*(x12-1/3*x12-1/3*x22-1/3*x32)+(-1/3*x12+x22-1/3*x22-1/3*x32)*(-1/3*x12+x22-1/3*x22-1/3*x32)+(-1/3*x12-1/3*x22+x32-1/3*x32)*(-1/3*x12-1/3*x22+x32-1/3*x32))"DispX[1,2]## [1] "0.333333333333333*((x11-1/3*x11-1/3*x21-1/3*x31)*(x12-1/3*x12-1/3*x22-1/3*x32)+(-1/3*x11+x21-1/3*x21-1/3*x31)*(-1/3*x12+x22-1/3*x22-1/3*x32)+(-1/3*x11-1/3*x21+x31-1/3*x31)*(-1/3*x12-1/3*x22+x32-1/3*x32))"# 数値的に確かめる。R の cov 関数は n-1 で割るため修正が必要。
seed0=1701
Xn=matrix(rnorm(nrow(X)*ncol(X)), nrow(X), ncol(X))
Eval(DispX, X=Xn) |> round(2)## [,1] [,2]
## [1,] 2.54 -0.84
## [2,] -0.84 0.28( (nrow(X)-1)/nrow(X)*cov(Xn) ) |> round(2)## [,1] [,2]
## [1,] 2.54 -0.84
## [2,] -0.84 0.28# 列方向の標準偏差： standardization
( V=Diag(DispX) )## [,1] [,2]
## [1,] "0.333333333333333*((x11-1/3*x11-1/3*x21-1/3*x31)*(x11-1/3*x11-1/3*x21-1/3*x31)+(-1/3*x11+x21-1/3*x21-1/3*x31)*(-1/3*x11+x21-1/3*x21-1/3*x31)+(-1/3*x11-1/3*x21+x31-1/3*x31)*(-1/3*x11-1/3*x21+x31-1/3*x31))" "0"
## [2,] "0" "0.333333333333333*((x12-1/3*x12-1/3*x22-1/3*x32)*(x12-1/3*x12-1/3*x22-1/3*x32)+(-1/3*x12+x22-1/3*x22-1/3*x32)*(-1/3*x12+x22-1/3*x22-1/3*x32)+(-1/3*x12-1/3*x22+x32-1/3*x32)*(-1/3*x12-1/3*x22+x32-1/3*x32))"# 標準偏差を対角要素に持つ行列： Sd=sqrt(V)=sqrt( (1/n)*Diag(t(JX) %*% JX))
# symbolic な sqrt 関数は未定義なので s_{ii} = sqrt(v_{ii}) とし、
# それを対角要素に持つ対角行列を Sd とする。
Sd=demomat(2,2,root="s",shape="diag")
# 標準偏差の逆数を対角要素として持つ行列を invSd とする。
# まだ Inv 関数が定義されていないのでとりあえず以下の様にする。
invSd=gsub("(s[0-9]+)","1/\\1",Sd) # = Inv(Sd)
# 標準化：中心化した後に各列を対応する標準偏差で割る。
ZX=JX %*% invSd
print(ZX)## [,1] [,2]
## [1,] "(x11-1/3*x11-1/3*x21-1/3*x31)/s11" "(x12-1/3*x12-1/3*x22-1/3*x32)/s22"
## [2,] "(-1/3*x11+x21-1/3*x21-1/3*x31)/s11" "(-1/3*x12+x22-1/3*x22-1/3*x32)/s22"
## [3,] "(-1/3*x11-1/3*x21+x31-1/3*x31)/s11" "(-1/3*x12-1/3*x22+x32-1/3*x32)/s22"# 相関係数を計算して数値的に確かめる。
Eval(1/nrow(ZX)*t(ZX)%*%ZX, X=Xn, s11=sd(Xn[,1]), s22=sd(Xn[,2]) ) |> round(2)## [,1] [,2]
## [1,] 0.67 -0.67
## [2,] -0.67 0.67( (nrow(X)-1)/nrow(X)*cor(Xn) ) |> round(2)## [,1] [,2]
## [1,] 0.67 -0.67
## [2,] -0.67 0.67推移核行列 \boldsymbol K の i-j 要素、k_{ij}、 は、状態 j から 状態 i への 推移確率を示す。各状態は互いに背反で全てを尽くすため \sum_{i=1}^{nrow(K)} k_{ij} = 1 である。 初期状態のよらず、推移を繰り返すと最終的に定常状態に達するが、 そこでの状態の割合は、 \boldsymbol K 行列の固有値 1 に対応する固有ベクトルを その和が 1 になるように基準化したものとなる。 なお、推移各行列を行から列への推移を表すように定義する場合があるが、 その時は推移各行列の行和が 1 となり、状態ベクトルの転置を左から掛けることになる。
# set output line width
options(width=500)
# K 推移核行列
# K[row,col] = Prob( from col to row )
K=t( matrix(c(
0.5, 0.3, 0.2, 0.1, 0.6, 0.3, 0.1, 0.1, 0.8
),3,3, byrow=1) )
colSums(K)## [1] 1 1 1 # 状態ベクトルの初期値
p=matrix(c(1,0,0),3,1)
# next
for( rep in 1:10 ){
Print( rep, p )
p=K%*%p
}## c1
## rep = 1 , p = r1 1
## r2 0
## r3 0
## c1
## rep = 2 , p = r1 0.5
## r2 0.3
## r3 0.2
## c1
## rep = 3 , p = r1 0.30
## r2 0.35
## r3 0.35
## c1
## rep = 4 , p = r1 0.220
## r2 0.335
## r3 0.445
## c1
## rep = 5 , p = r1 0.1880
## r2 0.3115
## r3 0.5005
## c1
## rep = 6 , p = r1 0.17520
## r2 0.29335
## r3 0.53145
## c1
## rep = 7 , p = r1 0.170080
## r2 0.281715
## r3 0.548205
## c1
## rep = 8 , p = r1 0.1680320
## r2 0.2748735
## r3 0.5570945
## c1
## rep = 9 , p = r1 0.1672128
## r2 0.2710432
## r3 0.5617441
## c1
## rep = 10 , p = r1 0.1668851
## r2 0.2689641
## r3 0.5641507 # 収束先は K の固有値 1 に対応する固有ベクトルを
# 足して 1 になるように正規化したもの。
ev=eigen(K)
print(ev)## eigen() decomposition
## $values
## [1] 1.0 0.5 0.4
##
## $vectors
## [,1] [,2] [,3]
## [1,] 0.2571722 6.480563e-16 0.4082483
## [2,] 0.4114756 -7.071068e-01 -0.8164966
## [3,] 0.8743856 7.071068e-01 0.4082483 ev1=ev$vectors[,1]
ev1=ev1/sum(ev1)
Print( p,ev1)## c1
## p = r1 0.1667540 , ev1 = 0.1666667
## r2 0.2678591 0.2666667
## r3 0.5653869 0.5666667横軸を x1 軸、縦軸を x2 軸とし、 原点から点 \boldsymbol x=(x_1, x_2)' まで届くベクトルを、 反時計回りに \theta だけ回転する。 そのときの回転後の点の新しい座標 \boldsymbol y は \boldsymbol y = \boldsymbol T \boldsymbol x で与えられるが、この \boldsymbol T を回転行列と呼ぶ。 なお、このことは、 時計回りに座標を \theta だけ回転することと同じであるため、 \boldsymbol x の新しい座標での座標は \boldsymbol T' \boldsymbol x となる。
# set output line width
options(width=500)
# 回転行列： planar rotation matrix
options(width=500)
T=matrix(c("cos(theta)", "-sin(theta)"
, "sin(theta)", "cos(theta)"), 2,2, byrow=1)
T## [,1] [,2]
## [1,] "cos(theta)" "-sin(theta)"
## [2,] "sin(theta)" "cos(theta)"Eval(T, theta=pi/4, check=0)## [,1] [,2]
## [1,] 0.7071068 -0.7071068
## [2,] 0.7071068 0.7071068# 回転行列は直交行列である。
tTT = t(T) %*% T |> print()## [,1] [,2]
## [1,] "cos(theta)*cos(theta)+sin(theta)*sin(theta)" "-cos(theta)*sin(theta)+sin(theta)*cos(theta)"
## [2,] "-sin(theta)*cos(theta)+cos(theta)*sin(theta)" "sin(theta)*sin(theta)+cos(theta)*cos(theta)"Eval(tTT, theta=pi/5, check=0)## [,1] [,2]
## [1,] 1 0
## [2,] 0 1cos(\theta_1+\theta_2) = cos(\theta_1)cos(\theta_2)-sin(\theta_1)sin(\theta_2) sin(\theta_1+\theta_2) = cos(\theta_1)sin(\theta_2)+sin(\theta_1)cos(\theta_2)
# set output line width
options(width=500)
# 反時計回りに theta1 だけ回転する回転行列
options(width=500)
T1=matrix(c("cos(theta1)", "-sin(theta1)"
, "sin(theta1)", "cos(theta1)"), 2,2, byrow=1)
# 反時計回りに theta2 だけ回転する回転行列
T2=matrix(c("cos(theta2)", "-sin(theta2)"
, "sin(theta2)", "cos(theta2)"), 2,2, byrow=1)
# まず theta1 回転し、それを次いで theta2 回転する時の回転行列はこの二つの積。
T2T1=T2 %*% T1
# 直接に theta1+theta2 回転する際の回転行列。
T12=matrix(c("cos(theta1+theta2)", "-sin(theta1+theta2)"
, "sin(theta1+theta2)", "cos(theta1+theta2)"), 2,2, byrow=1)
# この二つの行列の第１列を比べてみると加法定理を得る。
Print(T2T1, "/", T12)## T2T1 =
## c1 c2
## r1 cos(theta2)*cos(theta1)-sin(theta2)*sin(theta1) -cos(theta2)*sin(theta1)-sin(theta2)*cos(theta1)
## r2 sin(theta2)*cos(theta1)+cos(theta2)*sin(theta1) -sin(theta2)*sin(theta1)+cos(theta2)*cos(theta1)
##
## T12 =
## c1 c2
## r1 cos(theta1+theta2) -sin(theta1+theta2)
## r2 sin(theta1+theta2) cos(theta1+theta2)# この二つが等しいかの数値的検証
Eval(T12-T2T1, theta1=pi/2, theta2=pi/3, check=0 )## [,1] [,2]
## [1,] 0.000000e+00 5.551115e-17
## [2,] -5.551115e-17 0.000000e+00# ２次元平面での回転は交換可能 commutative !!
# まず theta2 回転し、その後 theta1 回転しても同じ。
T1T2=T1 %*% T2
# 数値的検証
dif=Simplify(Expand(T2T1-T1T2)) |> print()## [,1] [,2]
## [1,] "cos(theta2)*cos(theta1)-sin(theta2)*sin(theta1)-cos(theta1)*cos(theta2)+sin(theta1)*sin(theta2)" "-cos(theta2)*sin(theta1)-sin(theta2)*cos(theta1)+cos(theta1)*sin(theta2)+sin(theta1)*cos(theta2)"
## [2,] "sin(theta2)*cos(theta1)+cos(theta2)*sin(theta1)-sin(theta1)*cos(theta2)-cos(theta1)*sin(theta2)" "-sin(theta2)*sin(theta1)+cos(theta2)*cos(theta1)+sin(theta1)*sin(theta2)-cos(theta1)*cos(theta2)"Eval(dif, theta1=pi/2, theta2=pi/3, check=0 ) |> round(3)## [,1] [,2]
## [1,] 0 0
## [2,] 0 0# set output line width
options(width=500)
# x 軸は奥から手前に、y 軸は左から右へ、z 軸は下から上へ配置。
# x 軸の周りの回転： x 座標は不変 roll
options(width=9000)
temp <- c(
"1", "0", "0",
"0", "cos(theta32)", "-sin(theta32)",
"0", "sin(theta32)", "cos(theta32)"
)
Tx <- matrix(temp,3,3, byrow=1)
# y 軸の周りの回転： y 座標は不変 pitch
temp <- c(
"cos(theta31)", "0", "-sin(theta31)",
"0", "1", "0",
"sin(theta31)", "0", "cos(theta31)"
)
Ty <- matrix(temp,3,3, byrow=1)
# z 軸の周りの回転： z 座標は不変 yaw
temp <- c(
"cos(theta21)", "-sin(theta21)", "0",
"sin(theta21)", "cos(theta21)", "0",
"0", "0", "1"
)
Tz <- matrix(temp,3,3, byrow=1)
# これらを掛け合わせたものが３次元空間での回転行列となる。
T <- Tz %*% Ty %*% Tx |> Simplify()
# 幅広なので列ごとに印刷する
Print(T[,1, drop=0], "/", T[,2, drop=0], "\n", T[,3, drop=0])## T[, 1, drop = 0] =
## c1
## r1 cos(theta21)*cos(theta31)
## r2 sin(theta21)*cos(theta31)
## r3 sin(theta31)
##
## T[, 2, drop = 0] =
## c1
## r1 -sin(theta21)*cos(theta32)-cos(theta21)*sin(theta31)*sin(theta32)
## r2 cos(theta21)*cos(theta32)-sin(theta21)*sin(theta31)*sin(theta32)
## r3 cos(theta31)*sin(theta32)
##
## T[, 3, drop = 0] =
## c1
## r1 sin(theta21)*sin(theta32)-cos(theta21)*sin(theta31)*cos(theta32)
## r2 -cos(theta21)*sin(theta32)-sin(theta21)*sin(theta31)*cos(theta32)
## r3 cos(theta31)*cos(theta32)# たとえば以下の数値を入れたものは
tval="theta32=30/180*pi; theta31=45/180*pi; theta21=30/180*pi"
Tn <- Eval(T, values=tval, check=0)
Print(Tn)## Tn =
## c1 c2 c3
## r1 0.6123724 -0.7391989 -0.2803301
## r2 0.3535534 0.5732233 -0.7391989
## r3 0.7071068 0.3535534 0.6123724# 回転行列は直交行列である。
t(Tn) %*% Tn |> round(3)## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1Tn%*%t(Tn) |> round(3)## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1options(width=9999)
# 回転の順序
Txyz=Tx%*%Ty%*%Tz
Tyzx=Ty%*%Tz%*%Tx
Tzxy=Tz%*%Tx%*%Ty
Txzy=Tx%*%Tz%*%Ty
Tyxz=Ty%*%Tx%*%Tz
Tzyx=Tz%*%Ty%*%Tx
Txy=Tx%*%Ty
Tyx=Ty%*%Tx
Txz=Tx%*%Tz
Tzx=Tz%*%Tx
# 数値的検証
tval="theta32=30/180*pi; theta31=45/180*pi; theta21=60/180*pi"
Txyzn=Eval(Txyz,values=tval, check=0)
Tyzxn=Eval(Tyzx,values=tval, check=0)
Tzxyn=Eval(Tzxy,values=tval, check=0)
Txzyn=Eval(Txzy,values=tval, check=0)
Tyxzn=Eval(Tyxz,values=tval, check=0)
Tzyxn=Eval(Tzyx,values=tval, check=0)
Print(Txyzn,Tyzxn,Tzxyn, fmt="5.2")## c1 c2 c3 c1 c2 c3 c1 c2 c3
## Txyzn = r1 0.35 -0.61 -0.71 , Tyzxn = r1 0.35 -0.88 -0.31 , Tzxyn = r1 0.66 -0.75 -0.05
## r2 0.57 0.74 -0.35 r2 0.87 0.43 -0.25 r2 0.44 0.43 -0.79
## r3 0.74 -0.28 0.61 r3 0.35 -0.18 0.92 r3 0.61 0.50 0.61Print(Txzyn,Tyxzn,Tzyxn, fmt="5.2")## c1 c2 c3 c1 c2 c3 c1 c2 c3
## Txzyn = r1 0.35 -0.87 -0.35 , Tyxzn = r1 0.05 -0.79 -0.61 , Tzyxn = r1 0.35 -0.93 0.13
## r2 0.18 0.43 -0.88 r2 0.75 0.43 -0.50 r2 0.61 0.13 -0.78
## r3 0.92 0.25 0.31 r3 0.66 -0.44 0.61 r3 0.71 0.35 0.61Txyn=Eval(Txy,values=tval, check=0)
Tyxn=Eval(Tyx,values=tval, check=0)
Txzn=Eval(Txz,values=tval, check=0)
Tzxn=Eval(Tzx,values=tval, check=0)
Print(Txyn,Tyxn, fmt="5.2")## c1 c2 c3 c1 c2 c3
## Txyn = r1 0.71 0.00 -0.71 , Tyxn = r1 0.71 -0.35 -0.61
## r2 -0.35 0.87 -0.35 r2 0.00 0.87 -0.50
## r3 0.61 0.50 0.61 r3 0.71 0.35 0.61Print(Txzn, Tzxn, fmt="5.2")## c1 c2 c3 c1 c2 c3
## Txzn = r1 0.50 -0.87 0.00 , Tzxn = r1 0.50 -0.75 0.43
## r2 0.75 0.43 -0.50 r2 0.87 0.43 -0.25
## r3 0.43 0.25 0.87 r3 0.00 0.50 0.87# しかし直交行列ではある。
t(Txyzn)%*%Txyzn |> round(2)## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1t(Tyzxn)%*%Tyzxn |> round(2)## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1正方行列 A に関して、もしも、 A %% B と B %% A が単位行列になるなら B は A の逆行列と呼ぶ。
# set output line width
options(width=500)
# 対角行列の逆行列は対角要素をその逆数に置き換えたもの。
D3=diag(3); diag(D3)=paste("d",1:3,sep="")
invD3=diag(3); diag(invD3)=paste("1/d",1:3,sep="")
Print(D3, invD3)## c1 c2 c3 c1 c2 c3
## D3 = r1 d1 0 0 , invD3 = r1 1/d1 0 0
## r2 0 d2 0 r2 0 1/d2 0
## r3 0 0 d3 r3 0 0 1/d3# 検証
D3 %*% invD3## [,1] [,2] [,3]
## [1,] "1" "0" "0"
## [2,] "0" "1" "0"
## [3,] "0" "0" "1"invD3 %*% D3 |> Expand()## [,1] [,2] [,3]
## [1,] "1" "0" "0"
## [2,] "0" "1" "0"
## [3,] "0" "0" "1"# set output line width
options(width=500)
# 2 x 2 の行列の逆行列の行列式による計算
# 2 x 2 の行列
A2=demomat(2,2,root="a")
detA2="a11*a22-a21*a12"
invA2d=(paste("1/(",detA2,")",sep=""))*matrix(c("a22","-a21","-a12","a11"),2,2)
Print(invA2d)## invA2d =
## c1 c2
## r1 1/(a11*a22-a21*a12)*a22 -1/(a11*a22-a21*a12)*a12
## r2 -1/(a11*a22-a21*a12)*a21 1/(a11*a22-a21*a12)*a11# 数値的検証: A %*% invA ならびに invA%*%A が単位行列になるか否か
A2invA2=A2%*%invA2d
invA2A2=invA2d%*%A2
Print(A2invA2, "\n", invA2A2)## A2invA2 =
## c1 c2
## r1 a11/(a11*a22-a21*a12)*a22-a12/(a11*a22-a21*a12)*a21 -a11/(a11*a22-a21*a12)*a12+a12/(a11*a22-a21*a12)*a11
## r2 a21/(a11*a22-a21*a12)*a22-a22/(a11*a22-a21*a12)*a21 -a21/(a11*a22-a21*a12)*a12+a22/(a11*a22-a21*a12)*a11
##
## invA2A2 =
## c1 c2
## r1 1/(a11*a22-a21*a12)*a22*a11-1/(a11*a22-a21*a12)*a12*a21 1/(a11*a22-a21*a12)*a22*a12-1/(a11*a22-a21*a12)*a12*a22
## r2 -1/(a11*a22-a21*a12)*a21*a11+1/(a11*a22-a21*a12)*a11*a21 -1/(a11*a22-a21*a12)*a21*a12+1/(a11*a22-a21*a12)*a11*a22# 数値的に検証
set.seed(1701)
A2n=matrix(rnorm(nrow(A2)*ncol(A2)),nrow(A2))
Eval( A2invA2, A2=A2n ) |> round(4)## [,1] [,2]
## [1,] 1 0
## [2,] 0 1Eval( invA2A2, A2=A2n ) |> round(4)## [,1] [,2]
## [1,] 1 0
## [2,] 0 1スィープ演算子を用いれば、どのような大きさの記号的行列の逆行列でも 計算することが可能であるが、n >= 5 の場合には時間が掛かり実用的ではない。
# set output line width
options(width=500)
# 逆行列
invA2=matSwp(A2) |> print()## c1 c2
## r1 "1/a11-(a12/a11)*(-a21/a11)/(a22-a21*a12/a11)" "-(a12/a11)/(a22-a21*a12/a11)"
## r2 "(-a21/a11)/(a22-a21*a12/a11)" "1/(a22-a21*a12/a11)"上式は一見 invA2d と異なるように見えるが、各要素を整理すると同じになる｡ 以下にはそれを数値的に検証する。
Eval(invA2d, A2=A2n)## [,1] [,2]
## [1,] 0.2158860 -0.7769990
## [2,] 0.8116729 -0.5764779Eval(invA2, A2=A2n)## [,1] [,2]
## [1,] 0.2158860 -0.7769990
## [2,] 0.8116729 -0.5764779また、それが本当に A 行列の逆行列になっているかを数値的に検証する。
# 数値的検証: A %*% invA ならびに invA%*%A が単位行列になるか否か
A2invA2=A2%*%invA2
Print(A2invA2)## A2invA2 =
## c1 c2
## r1 a11*(1/a11-(a12/a11)*(-a21/a11)/(a22-a21*a12/a11))+a12*(-a21/a11)/(a22-a21*a12/a11) -a11*(a12/a11)/(a22-a21*a12/a11)+a12/(a22-a21*a12/a11)
## r2 a21*(1/a11-(a12/a11)*(-a21/a11)/(a22-a21*a12/a11))+a22*(-a21/a11)/(a22-a21*a12/a11) -a21*(a12/a11)/(a22-a21*a12/a11)+a22/(a22-a21*a12/a11)# 数値的に検証
Eval(A2invA2, A2=A2n) |> round(4)## [,1] [,2]
## [1,] 1 0
## [2,] 0 1invA2A2=invA2%*%A2
Print(invA2A2)## invA2A2 =
## c1 c2
## r1 (1/a11-(a12/a11)*(-a21/a11)/(a22-a21*a12/a11))*a11-(a12/a11)/(a22-a21*a12/a11)*a21 (1/a11-(a12/a11)*(-a21/a11)/(a22-a21*a12/a11))*a12-(a12/a11)/(a22-a21*a12/a11)*a22
## r2 (-a21/a11)/(a22-a21*a12/a11)*a11+1/(a22-a21*a12/a11)*a21 (-a21/a11)/(a22-a21*a12/a11)*a12+1/(a22-a21*a12/a11)*a22# 数値的に検証
Eval(invA2A2, A2=A2n) |> round(4)## [,1] [,2]
## [1,] 1 0
## [2,] 0 1以下に 3 x 3 の A 行列の例を示す。
# 3 x 3 の行列
set.seed(1701)
A3=demomat(3,3,root="a")
A3n=matrix(sample(0:9,9,replace=1),3,3)
invA3=matSwp(A3)
A3invA3=A3%*%invA3
# 数値的に検証
Eval(A3invA3, A3=A3n) |> round(4)## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1# 数値的に検証
invA3A3=invA3%*%A3
Eval(invA3A3, A3=A3n) |> round(4)## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1また、 n x n の A 行列の例を示すが n >= 5 の場合は時間が掛かる。
# set output line width
options(width=500)
# n x n の行列
n=4
set.seed(1701)
An=demomat(n,n,root="a")
Ann=matrix(sample(0:n^2,n^2,replace=1),n,n)
invAn=matSwp(An)
AninvAn=An%*%invAn
# 数値的に検証
Eval(AninvAn, An=Ann) |> round(4)## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 1 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 1# 数値的に検証
invAnAn=invAn%*%An
Eval(invAnAn, An=Ann) |> round(4)## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 1 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 1# set output line width
options(width=500)
# 三角行列
n=3
L=demomat(n,n,root="L",shape="lowert")
U=demomat(n,n,root="U",shape="uppert")
Print(L,U)## c1 c2 c3 c1 c2 c3
## L = r1 L11 0 0 , U = r1 U11 U12 U13
## r2 L21 L22 0 r2 0 U22 U23
## r3 L31 L32 L33 r3 0 0 U33# 三角行列の逆行列は三角行列となる。
invL=matSwp(L) |> Expand()
invU=matSwp(U) |> Expand()
Print(invL, "\n", invU)## invL =
## c1 c2 c3
## r1 1/L11 0 0
## r2 -L21/L11/L22 1/L22 0
## r3 -L31/L11/L33+L32*L21/L11/L22/L33 -L32/L22/L33 1/L33
##
## invU =
## c1 c2 c3
## r1 1/U11 -U12/U11/U22 -U13/U11/U33+U12/U11*U23/U22/U33
## r2 0 1/U22 -U23/U22/U33
## r3 0 0 1/U33LinvL=L%*%invL |> Expand() |> Simplify() |> print()## c1 c2 c3
## r1 "1" "0" "0"
## r2 "0" "1" "0"
## r3 "0" "0" "1"invLL=invL%*%L |> Expand() |> Simplify() |> print()## c1 c2 c3
## r1 "1" "0" "0"
## r2 "-L21/L22+1/L22*L21" "1" "0"
## r3 "-L31/L33+L32*L21/L22/L33-L32/L22/L33*L21+1/L33*L31" "-L32/L33+1/L33*L32" "1"# 数値的検証
set.seed(1701)
Ln=matrix(sample(1:n^2, replace=1),n,n)
Ln[upper.tri(L)]=0
Eval(invLL, L=Ln) |> round(4)## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1UinvU=U%*%invU |> Expand() |> Simplify() |> print()## c1 c2 c3
## r1 "1" "0" "0"
## r2 "0" "1" "0"
## r3 "0" "0" "1"invUU=invU%*%U |> Expand() |> Simplify() |> print()## c1 c2 c3
## r1 "1" "1/U11*U12-U12/U11" "1/U11*U13-U12/U11/U22*U23-U13/U11+U12/U11*U23/U22"
## r2 "0" "1" "1/U22*U23-U23/U22"
## r3 "0" "0" "1"# 数値的検証
set.seed(1701)
Un=matrix(sample(1:n^2, replace=1),n,n)
Un[lower.tri(U)]=0
Eval(invUU, U=Un) |> round(4)## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1以下に係数行列が正則の場合の例を示す。 詳細は LinearEq2.Rmd を参照のこと。
options(width=999)
# 連立方程式
n=2
An=demomat(n,n,root="a")
bn=demomat(n,1,vec=1,root="b")
x=demomat(n,1,vec=1,root="x")
# 方程式 Ax=b
Ax=An %*% x
Print(Ax,bn)## c1 c1
## Ax = r1 a11*x1+a12*x2 , bn = r1 b1
## r2 a21*x1+a22*x2 r2 b2# 解
# 方程式の両辺の左から Inv(A) を掛ける。
# Inv(A) A x = Inv(A) b
# ~~~ ~~
# x = Inv(A) b
xx=matSwp(An) %*% bn
Print(xx)## xx =
## c1
## r1 (1/a11-(a12/a11)*(-a21/a11)/(a22-a21*a12/a11))*b1-(a12/a11)/(a22-a21*a12/a11)*b2
## r2 (-a21/a11)/(a22-a21*a12/a11)*b1+1/(a22-a21*a12/a11)*b2# 検算
# これが b に等しくなるはず。
Axx=An %*% xx
Print(Axx)## Axx =
## c1
## r1 a11*((1/a11-(a12/a11)*(-a21/a11)/(a22-a21*a12/a11))*b1-(a12/a11)/(a22-a21*a12/a11)*b2)+a12*((-a21/a11)/(a22-a21*a12/a11)*b1+1/(a22-a21*a12/a11)*b2)
## r2 a21*((1/a11-(a12/a11)*(-a21/a11)/(a22-a21*a12/a11))*b1-(a12/a11)/(a22-a21*a12/a11)*b2)+a22*((-a21/a11)/(a22-a21*a12/a11)*b1+1/(a22-a21*a12/a11)*b2)# 数値的検証
set.seed(1701)
Ann=matrix(sample(0:9,n^2,replace=1),n,n)
bnn=matrix(sample(0:9,n,replace=1),n)
Print(Ann,bnn,Eval( Axx, An=Ann, bn=bnn ))## c1 c2 c1 c1
## Ann = r1 7 6 , bnn = r1 5 , Eval(Axx, An = Ann, bn = bnn) = r1 5
## r2 6 3 r2 7 r2 7階数が q の r x c の行列 \boldsymbol A に関して、c x r の行列 \boldsymbol A^{-} が \boldsymbol A=\boldsymbol A\boldsymbol A^{-}\boldsymbol A を満たすならば、 \boldsymbol A^{-} を \boldsymbol A の一般逆行列 (g-inverse, g1-inverse, conditional inverse, generalized inverse) という。
一般逆行列の求め方は様々なものがあるが、その行と列を並べ替行列 \boldsymbol P と \boldsymbol Q を用いて \boldsymbol B=\boldsymbol P\boldsymbol A\boldsymbol Q と並べ替えた行列 \boldsymbol B を作り、 その最初の q 行 q 列の部分に階数 q の行列 \boldsymbol B_{11} が入るようにしたものを \boldsymbol B=\begin{bmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{bmatrix} とする。すると、\boldsymbol B の一般逆行列 \boldsymbol B^{-} は \boldsymbol B^{-}=\begin{bmatrix} B_{11}^{-1} & \mathfrak{\Large 0}_{21} \\ \mathfrak{\Large 0}_{12} & \mathfrak{\Large 0}_{22} \end{bmatrix} で与えられる。したがって、\boldsymbol A\ 行列の一般逆行列は\boldsymbol A^{-}=\boldsymbol Q\boldsymbol B^{-} \boldsymbol P となる。
以下に、\boldsymbol A が 2 x 2 の行列でその 1,1 要素以外が 0 の場合の例を示す。 この場合は、行と列を並べ替える必要はなく \boldsymbol B=\boldsymbol A である。
# set output line width
options(width=500)
# A 行列とその一般逆行列
A=demomat(2,2,root="a")
A[2,]=0
A[,2]=0
ginvA=A
ginvA[1,1]="1/a11"
Print(A,ginvA)## c1 c2 c1 c2
## A = r1 a11 0 , ginvA = r1 1/a11 0
## r2 0 0 r2 0 0# ginvA が A の一般逆行列であることの検証
A%*%ginvA%*%A |> print()## c1 c2
## r1 "a11" "0"
## r2 "0" "0"ginvA%*%A%*%ginvA |> print()## c1 c2
## r1 "1/a11" "0"
## r2 "0" "0"次に、より一般的に、\boldsymbol B^{-} が \boldsymbol B の一般逆行列であることを示す。 この場合、\boldsymbol B 行列の第 2 列は第 1 列に線形従属であるため、 \boldsymbol C 行列を係数行列として \boldsymbol B_{12}=\boldsymbol B_{11} \boldsymbol C と \boldsymbol B_{22}=\boldsymbol B_{21} \boldsymbol C と書くことが出来る。 なお、\boldsymbol B 行列を表す B 行列は完全記号的行列としているが、 その詳細は、以下のブロック行列の節を参照されたい。
B=demomat( 2,2, root="B", fullsymb=1 )
B[1,2]="B11%*%C"
B[2,2]="B21%*%C"
ginvB=matSwp(B,1)
ginvB[,2]=ginvB[2,]=0
Print(B,ginvB)## c1 c2 c1 c2
## B = r1 B11 B11%*%C , ginvB = r1 inv(B11) 0
## r2 B21 B21%*%C r2 0 0BginvBB=B%*%ginvB%*%B |> print()## c1 c2
## r1 "B11" "B11%*%C"
## r2 "B21" "B21%*%C"
## attr(,"fullsymb")
## [1] 1ginvBBginvB=ginvB%*%B%*%ginvB |> print()## c1 c2
## r1 "inv(B11)" "0"
## r2 "0" "0"
## attr(,"fullsymb")
## [1] 1\boldsymbol P(\boldsymbol A) = \boldsymbol A(\boldsymbol A'\boldsymbol A)^{-}\boldsymbol A' を \boldsymbol A の列空間への直交射影演算子 （直交射影子）と呼ぶが、この行列は対称かつ冪等である。 また、 \boldsymbol P(\boldsymbol A) \; \boldsymbol A= \boldsymbol A である。 なお、\boldsymbol A が列フルランクの場合には一般逆行列でなく逆行列を用いて \boldsymbol P(\boldsymbol A) = \boldsymbol A(\boldsymbol A'\boldsymbol A)^{-1}\boldsymbol A' となる。以下に列フルランクの場合の証明を記すが、一般逆行列を用いる場合の証明は Harville (1997) を参照のこと。
# set output line width
options(width=500)
# 列フルランクの場合
PA="A%*%inv(t(A)%*%A)%*%t(A)"
attributes(PA)$fullsymb=1
PAA=PA%*%A |> print()## c1 c2
## [1,] "A%*%inv(t(A)%*%A)%*%t(A)%*%a11" "0"
## [2,] "0" "0"
## attr(,"fullsymb")
## [1] 1# 逆行列の性質を使うと
gsub("inv(t(A)%*%A)%*%t(A)%*%A", "I", PAA, fixed=1) ## c1 c2
## [1,] "A%*%inv(t(A)%*%A)%*%t(A)%*%a11" "0"
## [2,] "0" "0"
## attr(,"fullsymb")
## [1] 1# 列フルランクでは無い場合
PA="A%*%ginv(t(A)%*%A)%*%t(A)"
attributes(PA)$fullsymb=1
# 冪等性の証明： PA%*%A=A を利用する。
PA2=PA%*%PA |> print()## [,1]
## [1,] "A%*%ginv(t(A)%*%A)%*%t(A)%*%A%*%ginv(t(A)%*%A)%*%t(A)"
## attr(,"fullsymb")
## [1] 1# PA%*%A=A を利用する。
PA2=gsub("A%*%ginv(t(A)%*%A)%*%t(A)%*%A", "A", PA2, fixed=1) |> print()## [,1]
## [1,] "A%*%ginv(t(A)%*%A)%*%t(A)"
## attr(,"fullsymb")
## [1] 1\boldsymbol Q(\boldsymbol A)=I-\boldsymbol P(\boldsymbol A) を \boldsymbol A の列空間の直交補空間への射影演算子と呼ぶが、 \boldsymbol P(\boldsymbol A) \boldsymbol Q(\boldsymbol A)=\boldsymbol Q(\boldsymbol A) \boldsymbol P(\boldsymbol A)=\mathfrak{\Large 0} である。
# set output line width
options(width=500)
QA="I"-PA |> print()## [1] "A%*%ginv(t(A)%*%A)%*%t(A)"
## attr(,"fullsymb")
## [1] 1PQ=PA%*%QA |> Expand() |> print()## [,1]
## [1,] "A%*%ginv(t(A)%*%A)%*%t(A)-A%*%ginv(t(A)%*%A)%*%t(A)%*%A%*%ginv(t(A)%*%A)%*%t(A)"
## attr(,"fullsymb")
## [1] 1# PA%*%A=A を利用する。
gsub("A%*%ginv(t(A)%*%A)%*%t(A)%*%A", "A", PQ, fixed=1) |> Simplify()## [,1]
## [1,] "0"
## attr(,"fullsymb")
## [1] 1\boldsymbol 1 の直交補空間への射影子
# set output line width
options(width=500)
# 中心化行列: J = I_n - (1/n) * one_n %*% t(one_n)
# t(one_n)%*%one_n = n であるから、
( J3=I3 - "(1/3)" * one3 %*% t(one3) )## [,1] [,2] [,3]
## [1,] "1-1/3" "-1/3" "-1/3"
## [2,] "-1/3" "1-1/3" "-1/3"
## [3,] "-1/3" "-1/3" "1-1/3"# set output line width
options(width=500)
A="A"
attributes(A)$fullsymb=1
ginvA1="ginv(t(A)%*%A)%*%t(A)" |> print()## [1] "ginv(t(A)%*%A)%*%t(A)"A%*%ginvA1%*%A## [,1]
## [1,] "A%*%ginv(t(A)%*%A)%*%t(A)%*%A"
## attr(,"fullsymb")
## [1] 1二つのベクトルや行列の間のクロネッカー積 A %@% B は、 左側に現れる A の各要素に右側に現れる B を掛けたものを 左側の形に並べたものとなる。 クロネッカー積の結果を表示するためには printb 関数を使うとわかりやすい。
# set output line width
options(width=500)
A <- demomat(4,3,root="a") # 4 x 3 の行列
B <- demomat(4,3,root="b")
# 縦ベクトルと縦ベクトル
printb(v %@% u)## %@% v u
## c1-1
## r1-1 v1*u1
## r2-2 v1*u2
## r3-3 v1*u3
##
## r1-4 v2*u1
## r2-5 v2*u2
## r3-6 v2*u3
##
## r1-7 v3*u1
## r2-8 v3*u2
## r3-9 v3*u3rbind( v[1]*u, v[2]*u, v[3]*u )## [,1]
## [1,] "v1*u1"
## [2,] "v1*u2"
## [3,] "v1*u3"
## [4,] "v2*u1"
## [5,] "v2*u2"
## [6,] "v2*u3"
## [7,] "v3*u1"
## [8,] "v3*u2"
## [9,] "v3*u3"# 縦ベクトルと横ベクトル
printb(v %@% t(u))## %@% v t(u)
## c1-1 c2-2 c3-3
## r1-1 v1*u1 v1*u2 v1*u3
##
## r1-2 v2*u1 v2*u2 v2*u3
##
## r1-3 v3*u1 v3*u2 v3*u3rbind( v[1]*t(u), v[2]*t(u), v[3]*t(u) )## [,1] [,2] [,3]
## [1,] "v1*u1" "v1*u2" "v1*u3"
## [2,] "v2*u1" "v2*u2" "v2*u3"
## [3,] "v3*u1" "v3*u2" "v3*u3"# 横ベクトルと縦ベクトル
printb(t(v) %@% u)## %@% t(v) u
## c1-1 c1-2 c1-3
## r1-1 v1*u1 v2*u1 v3*u1
## r2-2 v1*u2 v2*u2 v3*u2
## r3-3 v1*u3 v2*u3 v3*u3cbind( v[1]*u, v[2]*u, v[3]*u )## [,1] [,2] [,3]
## [1,] "v1*u1" "v2*u1" "v3*u1"
## [2,] "v1*u2" "v2*u2" "v3*u2"
## [3,] "v1*u3" "v2*u3" "v3*u3"# 横ベクトルと横ベクトル
printb(t(v) %@% t(u))## %@% t(v) t(u)
## c1-1 c2-2 c3-3 c1-4 c2-5 c3-6 c1-7 c2-8 c3-9
## r1-1 v1*u1 v1*u2 v1*u3 v2*u1 v2*u2 v2*u3 v3*u1 v3*u2 v3*u3matrix(c( v[1]*u, v[2]*u, v[3]*u ),1)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] "v1*u1" "v1*u2" "v1*u3" "v2*u1" "v2*u2" "v2*u3" "v3*u1" "v3*u2" "v3*u3"# ベクトルと行列
v %@% A |> printb()## %@% v A
## c1-1 c2-2 c3-3
## r1-1 v1*a11 v1*a12 v1*a13
## r2-2 v1*a21 v1*a22 v1*a23
## r3-3 v1*a31 v1*a32 v1*a33
## r4-4 v1*a41 v1*a42 v1*a43
##
## r1-5 v2*a11 v2*a12 v2*a13
## r2-6 v2*a21 v2*a22 v2*a23
## r3-7 v2*a31 v2*a32 v2*a33
## r4-8 v2*a41 v2*a42 v2*a43
##
## r1-9 v3*a11 v3*a12 v3*a13
## r2-10 v3*a21 v3*a22 v3*a23
## r3-11 v3*a31 v3*a32 v3*a33
## r4-12 v3*a41 v3*a42 v3*a43rbind( v[1]*A, v[2]*A, v[3]*A )## c1 c2 c3
## r1 "v1*a11" "v1*a12" "v1*a13"
## r2 "v1*a21" "v1*a22" "v1*a23"
## r3 "v1*a31" "v1*a32" "v1*a33"
## r4 "v1*a41" "v1*a42" "v1*a43"
## r1 "v2*a11" "v2*a12" "v2*a13"
## r2 "v2*a21" "v2*a22" "v2*a23"
## r3 "v2*a31" "v2*a32" "v2*a33"
## r4 "v2*a41" "v2*a42" "v2*a43"
## r1 "v3*a11" "v3*a12" "v3*a13"
## r2 "v3*a21" "v3*a22" "v3*a23"
## r3 "v3*a31" "v3*a32" "v3*a33"
## r4 "v3*a41" "v3*a42" "v3*a43"# 行列とベクトル
A %@% v |> printb()## %@% A v
## c1-1 c1-2 c1-3
## r1-1 a11*v1 a12*v1 a13*v1
## r2-2 a11*v2 a12*v2 a13*v2
## r3-3 a11*v3 a12*v3 a13*v3
##
## r1-4 a21*v1 a22*v1 a23*v1
## r2-5 a21*v2 a22*v2 a23*v2
## r3-6 a21*v3 a22*v3 a23*v3
##
## r1-7 a31*v1 a32*v1 a33*v1
## r2-8 a31*v2 a32*v2 a33*v2
## r3-9 a31*v3 a32*v3 a33*v3
##
## r1-10 a41*v1 a42*v1 a43*v1
## r2-11 a41*v2 a42*v2 a43*v2
## r3-12 a41*v3 a42*v3 a43*v3cbind(
rbind( A[1,1]*v, A[2,1]*v, A[3,1]*v, A[4,1]*v ),
rbind( A[1,2]*v, A[2,2]*v, A[3,2]*v, A[4,2]*v ),
rbind( A[1,3]*v, A[2,3]*v, A[3,3]*v, A[4,3]*v )
)## [,1] [,2] [,3]
## [1,] "a11*v1" "a12*v1" "a13*v1"
## [2,] "a11*v2" "a12*v2" "a13*v2"
## [3,] "a11*v3" "a12*v3" "a13*v3"
## [4,] "a21*v1" "a22*v1" "a23*v1"
## [5,] "a21*v2" "a22*v2" "a23*v2"
## [6,] "a21*v3" "a22*v3" "a23*v3"
## [7,] "a31*v1" "a32*v1" "a33*v1"
## [8,] "a31*v2" "a32*v2" "a33*v2"
## [9,] "a31*v3" "a32*v3" "a33*v3"
## [10,] "a41*v1" "a42*v1" "a43*v1"
## [11,] "a41*v2" "a42*v2" "a43*v2"
## [12,] "a41*v3" "a42*v3" "a43*v3"# 行列と行列
( A %@% X ) |> Printb()## ( A %@% X
## c1-1 c2-2 c1-3 c2-4 c1-5 c2-6
## r1-1 a11*x11 a11*x12 a12*x11 a12*x12 a13*x11 a13*x12
## r2-2 a11*x21 a11*x22 a12*x21 a12*x22 a13*x21 a13*x22
## r3-3 a11*x31 a11*x32 a12*x31 a12*x32 a13*x31 a13*x32
##
## r1-4 a21*x11 a21*x12 a22*x11 a22*x12 a23*x11 a23*x12
## r2-5 a21*x21 a21*x22 a22*x21 a22*x22 a23*x21 a23*x22
## r3-6 a21*x31 a21*x32 a22*x31 a22*x32 a23*x31 a23*x32
##
## r1-7 a31*x11 a31*x12 a32*x11 a32*x12 a33*x11 a33*x12
## r2-8 a31*x21 a31*x22 a32*x21 a32*x22 a33*x21 a33*x22
## r3-9 a31*x31 a31*x32 a32*x31 a32*x32 a33*x31 a33*x32
##
## r1-10 a41*x11 a41*x12 a42*x11 a42*x12 a43*x11 a43*x12
## r2-11 a41*x21 a41*x22 a42*x21 a42*x22 a43*x21 a43*x22
## r3-12 a41*x31 a41*x32 a42*x31 a42*x32 a43*x31 a43*x32cbind(
rbind( A[1,1]*X, A[2,1]*X, A[3,1]*X, A[4,1]*X ),
rbind( A[1,2]*X, A[2,2]*X, A[3,2]*X, A[4,2]*X ),
rbind( A[1,3]*X, A[2,3]*X, A[3,3]*X, A[4,3]*X )
)## c1 c2 c1 c2 c1 c2
## r1 "a11*x11" "a11*x12" "a12*x11" "a12*x12" "a13*x11" "a13*x12"
## r2 "a11*x21" "a11*x22" "a12*x21" "a12*x22" "a13*x21" "a13*x22"
## r3 "a11*x31" "a11*x32" "a12*x31" "a12*x32" "a13*x31" "a13*x32"
## r1 "a21*x11" "a21*x12" "a22*x11" "a22*x12" "a23*x11" "a23*x12"
## r2 "a21*x21" "a21*x22" "a22*x21" "a22*x22" "a23*x21" "a23*x22"
## r3 "a21*x31" "a21*x32" "a22*x31" "a22*x32" "a23*x31" "a23*x32"
## r1 "a31*x11" "a31*x12" "a32*x11" "a32*x12" "a33*x11" "a33*x12"
## r2 "a31*x21" "a31*x22" "a32*x21" "a32*x22" "a33*x21" "a33*x22"
## r3 "a31*x31" "a31*x32" "a32*x31" "a32*x32" "a33*x31" "a33*x32"
## r1 "a41*x11" "a41*x12" "a42*x11" "a42*x12" "a43*x11" "a43*x12"
## r2 "a41*x21" "a41*x22" "a42*x21" "a42*x22" "a43*x21" "a43*x22"
## r3 "a41*x31" "a41*x32" "a42*x31" "a42*x32" "a43*x31" "a43*x32"vec 演算子は行列の要素を縦に引き延ばしてベクトルを作成する。 通常は行列の各列を縦に並べた形でベクトルが作成されるが、これを column major な方法と呼ぶ。したがって、転置行列を与えると row major になる。 正方行列を与えた場合に、その下半分のみを引き延ばしてベクトルを作る演算子は vech と呼ばれる。
# set output line width
options(width=500)
# vec 演算子: column-wise roll-out: column major
vec(A)## [,1]
## r1-c1 "a11"
## r2-c1 "a21"
## r3-c1 "a31"
## r4-c1 "a41"
## r1-c2 "a12"
## r2-c2 "a22"
## r3-c2 "a32"
## r4-c2 "a42"
## r1-c3 "a13"
## r2-c3 "a23"
## r3-c3 "a33"
## r4-c3 "a43"vec(t(A))## [,1]
## c1-r1 "a11"
## c2-r1 "a12"
## c3-r1 "a13"
## c1-r2 "a21"
## c2-r2 "a22"
## c3-r2 "a23"
## c1-r3 "a31"
## c2-r3 "a32"
## c3-r3 "a33"
## c1-r4 "a41"
## c2-r4 "a42"
## c3-r4 "a43"vec(A,byrow=1)## [,1]
## r1-c1 "a11"
## r1-c2 "a12"
## r1-c3 "a13"
## r2-c1 "a21"
## r2-c2 "a22"
## r2-c3 "a23"
## r3-c1 "a31"
## r3-c2 "a32"
## r3-c3 "a33"
## r4-c1 "a41"
## r4-c2 "a42"
## r4-c3 "a43"vec(S)## [,1]
## r1-c1 "s11"
## r2-c1 "s21"
## r3-c1 "s31"
## r1-c2 "s21"
## r2-c2 "s22"
## r3-c2 "s32"
## r1-c3 "s31"
## r2-c3 "s32"
## r3-c3 "s33"vech(S)## [,1]
## r1-c1 "s11"
## r2-c1 "s21"
## r2-c2 "s22"
## r3-c1 "s31"
## r3-c2 "s32"
## r3-c3 "s33"vech(S, nodiag=1)## [,1]
## r2-c1 "s21"
## r3-c1 "s31"
## r3-c2 "s32"vec(Y2)## [,1]
## r1-c1 "y11"
## r2-c1 "y21"
## r1-c2 "y12"
## r2-c2 "y22"
## r1-c3 "y13"
## r2-c3 "y23"vech(Y2)## error(vech): input matrix is not square.
## c1 c2 c3
## r1 "y11" "y12" "y13"
## r2 "y21" "y22" "y23"## NULL# set output line width
options(width=500)
# vec(a %*% t(b)) = b %@% a
u%*%t(v)## [,1] [,2] [,3]
## [1,] "u1*v1" "u1*v2" "u1*v3"
## [2,] "u2*v1" "u2*v2" "u2*v3"
## [3,] "u3*v1" "u3*v2" "u3*v3"vec(u%*%t(v))## [,1]
## [1,] "u1*v1"
## [2,] "u2*v1"
## [3,] "u3*v1"
## [4,] "u1*v2"
## [5,] "u2*v2"
## [6,] "u3*v2"
## [7,] "u1*v3"
## [8,] "u2*v3"
## [9,] "u3*v3"v%@%u |> printb()## %@% v u
## c1-1
## r1-1 v1*u1
## r2-2 v1*u2
## r3-3 v1*u3
##
## r1-4 v2*u1
## r2-5 v2*u2
## r3-6 v2*u3
##
## r1-7 v3*u1
## r2-8 v3*u2
## r3-9 v3*u3# 行列の対応する要素の積の和
tr( t(X) %*% Y )## [1] "x11*y11+x21*y21+x31*y31+x12*y12+x22*y22+x32*y32"t(vec(X)) %*% vec(Y)## [,1]
## [1,] "x11*y11+x21*y21+x31*y31+x12*y12+x22*y22+x32*y32"# vec( X %*% Y %*% Z) == ( t(Z) %@% X ) %*% vec(Y)
vecXYZ <- vec(X %*% Y2 %*% Z) |> Expand()
tZKXvY <- (t(Z)%@%X) %*% vec(Y2)
Print(vecXYZ,"\n",tZKXvY)## vecXYZ =
## c1
## r1 x11*y11*z1/sqrt(2*z1^2)+x12*y21*z1/sqrt(2*z1^2)-x11*y12*z1/sqrt(2*z1^2)-x12*y22*z1/sqrt(2*z1^2)
## r2 x21*y11*z1/sqrt(2*z1^2)+x22*y21*z1/sqrt(2*z1^2)-x21*y12*z1/sqrt(2*z1^2)-x22*y22*z1/sqrt(2*z1^2)
## r3 x31*y11*z1/sqrt(2*z1^2)+x32*y21*z1/sqrt(2*z1^2)-x31*y12*z1/sqrt(2*z1^2)-x32*y22*z1/sqrt(2*z1^2)
## r4 x11*y11*z1/sqrt(2*z1^2+z2^2)+x12*y21*z1/sqrt(2*z1^2+z2^2)+x11*y12*z1/sqrt(2*z1^2+z2^2)+x12*y22*z1/sqrt(2*z1^2+z2^2)+x11*y13*z2/sqrt(2*z1^2+z2^2)+x12*y23*z2/sqrt(2*z1^2+z2^2)
## r5 x21*y11*z1/sqrt(2*z1^2+z2^2)+x22*y21*z1/sqrt(2*z1^2+z2^2)+x21*y12*z1/sqrt(2*z1^2+z2^2)+x22*y22*z1/sqrt(2*z1^2+z2^2)+x21*y13*z2/sqrt(2*z1^2+z2^2)+x22*y23*z2/sqrt(2*z1^2+z2^2)
## r6 x31*y11*z1/sqrt(2*z1^2+z2^2)+x32*y21*z1/sqrt(2*z1^2+z2^2)+x31*y12*z1/sqrt(2*z1^2+z2^2)+x32*y22*z1/sqrt(2*z1^2+z2^2)+x31*y13*z2/sqrt(2*z1^2+z2^2)+x32*y23*z2/sqrt(2*z1^2+z2^2)
##
## tZKXvY =
## c1
## r1-1 z1/sqrt(2*z1^2)*x11*y11+z1/sqrt(2*z1^2)*x12*y21-z1/sqrt(2*z1^2)*x11*y12-z1/sqrt(2*z1^2)*x12*y22
## r2-2 z1/sqrt(2*z1^2)*x21*y11+z1/sqrt(2*z1^2)*x22*y21-z1/sqrt(2*z1^2)*x21*y12-z1/sqrt(2*z1^2)*x22*y22
## r3-3 z1/sqrt(2*z1^2)*x31*y11+z1/sqrt(2*z1^2)*x32*y21-z1/sqrt(2*z1^2)*x31*y12-z1/sqrt(2*z1^2)*x32*y22
## r1-4 z1/sqrt(2*z1^2+z2^2)*x11*y11+z1/sqrt(2*z1^2+z2^2)*x12*y21+z1/sqrt(2*z1^2+z2^2)*x11*y12+z1/sqrt(2*z1^2+z2^2)*x12*y22+z2/sqrt(2*z1^2+z2^2)*x11*y13+z2/sqrt(2*z1^2+z2^2)*x12*y23
## r2-5 z1/sqrt(2*z1^2+z2^2)*x21*y11+z1/sqrt(2*z1^2+z2^2)*x22*y21+z1/sqrt(2*z1^2+z2^2)*x21*y12+z1/sqrt(2*z1^2+z2^2)*x22*y22+z2/sqrt(2*z1^2+z2^2)*x21*y13+z2/sqrt(2*z1^2+z2^2)*x22*y23
## r3-6 z1/sqrt(2*z1^2+z2^2)*x31*y11+z1/sqrt(2*z1^2+z2^2)*x32*y21+z1/sqrt(2*z1^2+z2^2)*x31*y12+z1/sqrt(2*z1^2+z2^2)*x32*y22+z2/sqrt(2*z1^2+z2^2)*x31*y13+z2/sqrt(2*z1^2+z2^2)*x32*y23
## attr(,"SMdim")
## attr(,"SMdim")$nrows
## [1] 3 3
##
## attr(,"SMdim")$ncols
## NULL# 数値的に検証する。
Eval(vecXYZ - tZKXvY, X=demomat(3,2), Y2=demomat(2,3)*10, Z=demomat(3,2)+5)## [,1]
## [1,] 0
## [2,] 0
## [3,] 0
## [4,] 0
## [5,] 0
## [6,] 0# vec( t(X) %*% Y %*% t(Z) == ( Z %@% t(X) ) %*% vec(Y)
vecXYZ <- vec(t(X)%*%Y%*%t(Z))
tZKXvY <- (Z%@%t(X)) %*% vec(Y)
Print(vecXYZ,"\n",tZKXvY)## vecXYZ =
## c1
## r1 (x11*y11+x21*y21+x31*y31)*z1/sqrt(2*z1^2)+(x11*y12+x21*y22+x31*y32)*z1/sqrt(2*z1^2+z2^2)
## r2 (x12*y11+x22*y21+x32*y31)*z1/sqrt(2*z1^2)+(x12*y12+x22*y22+x32*y32)*z1/sqrt(2*z1^2+z2^2)
## r3 -(x11*y11+x21*y21+x31*y31)*z1/sqrt(2*z1^2)+(x11*y12+x21*y22+x31*y32)*z1/sqrt(2*z1^2+z2^2)
## r4 -(x12*y11+x22*y21+x32*y31)*z1/sqrt(2*z1^2)+(x12*y12+x22*y22+x32*y32)*z1/sqrt(2*z1^2+z2^2)
## r5 (x11*y12+x21*y22+x31*y32)*z2/sqrt(2*z1^2+z2^2)
## r6 (x12*y12+x22*y22+x32*y32)*z2/sqrt(2*z1^2+z2^2)
##
## tZKXvY =
## c1
## r1-1 z1/sqrt(2*z1^2)*x11*y11+z1/sqrt(2*z1^2)*x21*y21+z1/sqrt(2*z1^2)*x31*y31+z1/sqrt(2*z1^2+z2^2)*x11*y12+z1/sqrt(2*z1^2+z2^2)*x21*y22+z1/sqrt(2*z1^2+z2^2)*x31*y32
## r2-2 z1/sqrt(2*z1^2)*x12*y11+z1/sqrt(2*z1^2)*x22*y21+z1/sqrt(2*z1^2)*x32*y31+z1/sqrt(2*z1^2+z2^2)*x12*y12+z1/sqrt(2*z1^2+z2^2)*x22*y22+z1/sqrt(2*z1^2+z2^2)*x32*y32
## r1-3 -z1/sqrt(2*z1^2)*x11*y11-z1/sqrt(2*z1^2)*x21*y21-z1/sqrt(2*z1^2)*x31*y31+z1/sqrt(2*z1^2+z2^2)*x11*y12+z1/sqrt(2*z1^2+z2^2)*x21*y22+z1/sqrt(2*z1^2+z2^2)*x31*y32
## r2-4 -z1/sqrt(2*z1^2)*x12*y11-z1/sqrt(2*z1^2)*x22*y21-z1/sqrt(2*z1^2)*x32*y31+z1/sqrt(2*z1^2+z2^2)*x12*y12+z1/sqrt(2*z1^2+z2^2)*x22*y22+z1/sqrt(2*z1^2+z2^2)*x32*y32
## r1-5 z2/sqrt(2*z1^2+z2^2)*x11*y12+z2/sqrt(2*z1^2+z2^2)*x21*y22+z2/sqrt(2*z1^2+z2^2)*x31*y32
## r2-6 z2/sqrt(2*z1^2+z2^2)*x12*y12+z2/sqrt(2*z1^2+z2^2)*x22*y22+z2/sqrt(2*z1^2+z2^2)*x32*y32
## attr(,"SMdim")
## attr(,"SMdim")$nrows
## [1] 2 2 2
##
## attr(,"SMdim")$ncols
## NULL# 数値的に検証する。
Eval(vecXYZ - tZKXvY, X=demomat(3,2), Y=demomat(3,2)*3, Z=demomat(3,2))## [,1]
## [1,] 0
## [2,] 0
## [3,] 0
## [4,] 0
## [5,] 0
## [6,] 0# トレース、vec 演算子とクロネッカー積:
# tr(XYZ) == t(vec(t(X)) I%@%Y vec(Z)
trXYZ <- tr( X%*%Y2%*%Z2 )
vXetc <- t(vec(t(X))) %*% (diag(3)%@%Y2) %*% vec(Z2)
# 数値的に検証する。
Eval(trXYZ - vXetc, X=demomat(3,2), Y2=demomat(2,3)*10, Z2=demomat(3,3)+5)## [,1]
## [1,] 0# set output line width
options(width=500)
one3 %@% v |> printb()## %@% one3 v
## c1-1
## r1-1 v1
## r2-2 v2
## r3-3 v3
##
## r1-4 v1
## r2-5 v2
## r3-6 v3
##
## r1-7 v1
## r2-8 v2
## r3-9 v3v %@% one3 |> printb()## %@% v one3
## c1-1
## r1-1 v1
## r2-2 v1
## r3-3 v1
##
## r1-4 v2
## r2-5 v2
## r3-6 v2
##
## r1-7 v3
## r2-8 v3
## r3-9 v3t(one3) %@% v |> printb()## %@% t(one3) v
## c1-1 c1-2 c1-3
## r1-1 v1 v1 v1
## r2-2 v2 v2 v2
## r3-3 v3 v3 v3v %@% t(one3) |> printb()## %@% v t(one3)
## c1-1 c2-2 c3-3
## r1-1 v1 v1 v1
##
## r1-2 v2 v2 v2
##
## r1-3 v3 v3 v3one3 %@% A |> printb()## %@% one3 A
## c1-1 c2-2 c3-3
## r1-1 a11 a12 a13
## r2-2 a21 a22 a23
## r3-3 a31 a32 a33
## r4-4 a41 a42 a43
##
## r1-5 a11 a12 a13
## r2-6 a21 a22 a23
## r3-7 a31 a32 a33
## r4-8 a41 a42 a43
##
## r1-9 a11 a12 a13
## r2-10 a21 a22 a23
## r3-11 a31 a32 a33
## r4-12 a41 a42 a43A %@% one3 |> printb()## %@% A one3
## c1-1 c1-2 c1-3
## r1-1 a11 a12 a13
## r2-2 a11 a12 a13
## r3-3 a11 a12 a13
##
## r1-4 a21 a22 a23
## r2-5 a21 a22 a23
## r3-6 a21 a22 a23
##
## r1-7 a31 a32 a33
## r2-8 a31 a32 a33
## r3-9 a31 a32 a33
##
## r1-10 a41 a42 a43
## r2-11 a41 a42 a43
## r3-12 a41 a42 a43# 対角化？
I3 %@% A |> printb()## %@% I3 A
## c1-1 c2-2 c3-3 c1-4 c2-5 c3-6 c1-7 c2-8 c3-9
## r1-1 a11 a12 a13 0 0 0 0 0 0
## r2-2 a21 a22 a23 0 0 0 0 0 0
## r3-3 a31 a32 a33 0 0 0 0 0 0
## r4-4 a41 a42 a43 0 0 0 0 0 0
##
## r1-5 0 0 0 a11 a12 a13 0 0 0
## r2-6 0 0 0 a21 a22 a23 0 0 0
## r3-7 0 0 0 a31 a32 a33 0 0 0
## r4-8 0 0 0 a41 a42 a43 0 0 0
##
## r1-9 0 0 0 0 0 0 a11 a12 a13
## r2-10 0 0 0 0 0 0 a21 a22 a23
## r3-11 0 0 0 0 0 0 a31 a32 a33
## r4-12 0 0 0 0 0 0 a41 a42 a43A %@% I3 |> printb()## %@% A I3
## c1-1 c2-2 c3-3 c1-4 c2-5 c3-6 c1-7 c2-8 c3-9
## r1-1 a11 0 0 a12 0 0 a13 0 0
## r2-2 0 a11 0 0 a12 0 0 a13 0
## r3-3 0 0 a11 0 0 a12 0 0 a13
##
## r1-4 a21 0 0 a22 0 0 a23 0 0
## r2-5 0 a21 0 0 a22 0 0 a23 0
## r3-6 0 0 a21 0 0 a22 0 0 a23
##
## r1-7 a31 0 0 a32 0 0 a33 0 0
## r2-8 0 a31 0 0 a32 0 0 a33 0
## r3-9 0 0 a31 0 0 a32 0 0 a33
##
## r1-10 a41 0 0 a42 0 0 a43 0 0
## r2-11 0 a41 0 0 a42 0 0 a43 0
## r3-12 0 0 a41 0 0 a42 0 0 a43行列を要素として持つ行列をブロック行列 (block matrix) または 分割行列 (partitioned matrix) と呼ぶ。
以下の例では bigA や bigB 等の行列の要素である A や B は行列であると 見なしているが、そのサイズや要素は未定義である。このような記号的行列を 完全記号行列と呼んでいるが、それらの行列には fullsymb という attribute が 付随し、その値が 1 となっている。
# set output line width
options(width=500)
# ブロック行列、区分行列、分割行列
# block matrix, partitioned matrix
bigA=demomat(2,2,root="A", fullsymb=1) |> print()## c1 c2
## r1 "A11" "A12"
## r2 "A21" "A22"
## attr(,"fullsymb")
## [1] 1bigB=demomat(2,2,root="B", fullsymb=1) |> print()## c1 c2
## r1 "B11" "B12"
## r2 "B21" "B22"
## attr(,"fullsymb")
## [1] 1attributes(bigA)## $dim
## [1] 2 2
##
## $dimnames
## $dimnames[[1]]
## [1] "r1" "r2"
##
## $dimnames[[2]]
## [1] "c1" "c2"
##
##
## $fullsymb
## [1] 1# 転置行列：あたかも各要素をスカラーだと思って転置した後で、
# それぞれの行列も転置する。
t(bigA)## r1 r2
## c1 "t(A11)" "t(A21)"
## c2 "t(A12)" "t(A22)"
## attr(,"fullsymb")
## [1] 1# トレース：各要素のトレースの和
tr(bigA)## [,1]
## [1,] "tr(A11)+tr(A22)"# 行列の和と差：要素の和と差
bigA+bigB## c1 c2
## r1 "A11+B11" "A12+B12"
## r2 "A21+B21" "A22+B22"
## attr(,"fullsymb")
## [1] 1bigA-bigB## c1 c2
## r1 "A11-B11" "A12-B12"
## r2 "A21-B21" "A22-B22"
## attr(,"fullsymb")
## [1] 1# 行列の積：各要素を行列の掛け算を使って掛ける。
bigA%*%bigB## c1 c2
## r1 "A11%*%B11+A12%*%B21" "A11%*%B12+A12%*%B22"
## r2 "A21%*%B11+A22%*%B21" "A21%*%B12+A22%*%B22"
## attr(,"fullsymb")
## [1] 1t(bigA)%*%bigB## c1 c2
## c1 "t(A11)%*%B11+t(A21)%*%B21" "t(A11)%*%B12+t(A21)%*%B22"
## c2 "t(A12)%*%B11+t(A22)%*%B21" "t(A12)%*%B12+t(A22)%*%B22"
## attr(,"fullsymb")
## [1] 1# set output line width
options(width=500)
# 逆行列
invbigA=matSwp(bigA) |> print()## c1 c2
## r1 "inv(A11)-inv(A11)%*%A12%*%inv(A22-A21%*%inv(A11)%*%A12)%*%(-A21%*%inv(A11))" "-inv(A11)%*%A12%*%inv(A22-A21%*%inv(A11)%*%A12)"
## r2 "-inv(A22-A21%*%inv(A11)%*%A12)%*%A21%*%inv(A11)" "inv(A22-A21%*%inv(A11)%*%A12)"
## attr(,"fullsymb")
## [1] 1# 検証
bigAinvbigA=bigA%*%invbigA
bigAinvbigA=bigAinvbigA |> Expand() |> Expand() |> Simplify() |> Simplify()
Print(bigAinvbigA)## bigAinvbigA =
## c1 c2
## r1 I 0
## r2 A21%*%inv(A11)+A21%*%inv(A11)%*%A12%*%inv(A22-A21%*%inv(A11)%*%A12)%*%A21%*%inv(A11)-A22%*%inv(A22-A21%*%inv(A11)%*%A12)%*%A21%*%inv(A11) -A21%*%inv(A11)%*%A12%*%inv(A22-A21%*%inv(A11)%*%A12)+A22%*%inv(A22-A21%*%inv(A11)%*%A12)
## attr(,"fullsymb")
## [1] 1この行列の [1,1] 要素は I、 [1,2] 要素は 0 である。 また、[2,1] 要素は、スペースの省略のために行列の掛け算記号を # で書くと
A21#inv(A11)+A21#inv(A11)#A12#inv(A22-A21#inv(A11)#A12)#A21#inv(A11)-A22#inv(A22-A21#inv(A11)#A12)#A21#inv(A11)
であるが、これを右端に現れる inv(A22-A21#inv(A11)#A12)#A21#inv(A11) で括ると
A21#inv(A11)+A21#inv(A11)#A12#inv(A22-A21#inv(A11)#A12)#A21#inv(A11) - A22#inv(A22-A21#inv(A11)#A12)#A21#inv(A11)
= A21#inv(A11) + (A21#inv(A11)#A12 - A22) #inv(A22-A21#inv(A11)#A12)#A21#inv(A11)
= A21#inv(A11) - (A22-A21#inv(A11)#A12) #inv(A22-A21#inv(A11)#A12)#A21#inv(A11)
= A21#inv(A11) - I #A21#inv(A11) = 0
となる。 また、 [2,2] 要素は
-A21#inv(A11)#A12#inv(A22-A21#inv(A11)#A12)+A22#inv(A22-A21#inv(A11)#A12)
であるが、それを inv(A22-A21#inv(A11)#A12) で括ると
-A21#inv(A11)#A12#inv(A22-A21#inv(A11)#A12) + A22#inv(A22-A21#inv(A11)#A12)
= (-A21#inv(A11)#A12 + A22) #inv(A22-A21#inv(A11)#A12)
= -(A22-A21#inv(A11)#A12) #inv(A22-A21#inv(A11)#A12) = I
となり、結局 bigAinvbigA は単位行列であることが解る。
# set output line width
options(width=500)
# 検証その 2
invbigAbigA=invbigA%*%bigA
invbigAbigA=invbigAbigA |> Expand() |> Expand() |> Simplify() |> Simplify()
Print(invbigAbigA)## invbigAbigA =
## c1 c2
## r1 I inv(A11)%*%A12+inv(A11)%*%A12%*%inv(A22-A21%*%inv(A11)%*%A12)%*%A21%*%inv(A11)%*%A12-inv(A11)%*%A12%*%inv(A22-A21%*%inv(A11)%*%A12)%*%A22
## r2 0 -inv(A22-A21%*%inv(A11)%*%A12)%*%A21%*%inv(A11)%*%A12+inv(A22-A21%*%inv(A11)%*%A12)%*%A22
## attr(,"fullsymb")
## [1] 1この行列の [1,1] 要素は I、 [2,1] 要素は 0 である。 また、[1,2] 要素は、スペースの省略のために行列の掛け算記号を # で書くと
inv(A11)#A12+inv(A11)#A12#inv(A22-A21#inv(A11)#A12)#A21#inv(A11)#A12-inv(A11)#A12#inv(A22-A21#inv(A11)#A12)#A22
であるが、これを左端に現れる inv(A11)#A12#inv(A22-A21#inv(A11)#A12) で括ると
inv(A11)#A12+inv(A11)#A12#inv(A22-A21#inv(A11)#A12)#A21#inv(A11)#A12 - inv(A11)#A12#inv(A22-A21#inv(A11)#A12)#A22
= inv(A11)#A12 + inv(A11)#A12#inv(A22-A21#inv(A11)#A12) # ( A21#inv(A11)#A12 - A22 )
= inv(A11)#A12 - inv(A11)#A12#inv(A22-A21#inv(A11)#A12) # ( A22-A21#inv(A11)#A12 )
= inv(A11)#A12 - inv(A11)#A12#I = 0 を得る。
また、[2,2] 要素は、スペースの省略のために行列の掛け算記号を # で書くと
-inv(A22-A21#inv(A11)#A12)#A21#inv(A11)#A12+inv(A22-A21#inv(A11)#A12)#A22
であるが、これを左端に現れる inv(A22-A21#inv(A11)#A12) で括ると
-inv(A22-A21#inv(A11)#A12)#A21#inv(A11)#A12 + inv(A22-A21#inv(A11)#A12)#A22
= inv(A22-A21#inv(A11)#A12)# (-A21#inv(A11)#A12 + A22) = I
となり、こちらも単位行列となることが分かる。
ちなみに以下のものも bigA 行列の逆行列となっている。
# set output line width
options(width=500)
# 逆行列の別の表現
invbigA2=matSwp(bigA, 2:1) |> print()## c1 c2
## r1 "inv(A11-A12%*%inv(A22)%*%A21)" "-inv(A11-A12%*%inv(A22)%*%A21)%*%A12%*%inv(A22)"
## r2 "-inv(A22)%*%A21%*%inv(A11-A12%*%inv(A22)%*%A21)" "inv(A22)-inv(A22)%*%A21%*%inv(A11-A12%*%inv(A22)%*%A21)%*%(-A12%*%inv(A22))"
## attr(,"fullsymb")
## [1] 1したがって、 inv(A11)-inv(A11)%%A12%%inv(A22-A21%%inv(A11)%%A12)%%(-A21%%inv(A11)) = inv(A11-A12%%inv(A22)%%A21) という関係を得る。
# set output line width
options(width=500)
# ブロック対角行列
bigD=demomat(2,2,root="B", fullsymb=1, shape="diag") |> print()## [,1] [,2]
## [1,] "B11" "0"
## [2,] "0" "B22"
## attr(,"fullsymb")
## [1] 1# ブロック対角行列の逆行列は対角要素をその逆行列に置き換えたもの
invbigD=matSwp(bigD) |> print()## [,1] [,2]
## [1,] "inv(B11)" "0"
## [2,] "0" "inv(B22)"
## attr(,"fullsymb")
## [1] 1# 検証
bigD%*%invbigD## [,1] [,2]
## [1,] "I" "0"
## [2,] "0" "I"
## attr(,"fullsymb")
## [1] 1# set output line width
options(width=500)
# 2 x 2 の分割行列：行列を要素として持つ行列
A11=demomat(2,2,root="a11_")
A12=demomat(2,3,root="a12_")
A21=demomat(3,2,root="a21_")
A22=demomat(3,3,root="a22_")
bigA=rbind( cbind(A11,A12), cbind(A21,A22) )
( bigA0=demomat(c(2,3),c(2,3), root="a") ) |> printb()## ( bigA0 = demomat(c(2, 3), c(2, 3), root = "a")
## c1-1 c1-2 c2-1 c2-2 c2-3
## r1-1 a11_11 a11_12 a12_11 a12_12 a12_13
## r1-2 a11_21 a11_22 a12_21 a12_22 a12_23
##
## r2-1 a21_11 a21_12 a22_11 a22_12 a22_13
## r2-2 a21_21 a21_22 a22_21 a22_22 a22_23
## r2-3 a21_31 a21_32 a22_31 a22_32 a22_33# 2 x 2 の分割行列
B11=demomat(2,2,root="b11_")
B12=demomat(2,3,root="b12_")
B21=demomat(3,2,root="b21_")
B22=demomat(3,3,root="b22_")
bigB=rbind( cbind(B11,B12), cbind(B21,B22) )
bigB0=demomat(c(2,3),c(2,3), root="b") |> printb()## demomat c(2, 3) c(2, 3) b
## c1-1 c1-2 c2-1 c2-2 c2-3
## r1-1 b11_11 b11_12 b12_11 b12_12 b12_13
## r1-2 b11_21 b11_22 b12_21 b12_22 b12_23
##
## r2-1 b21_11 b21_12 b22_11 b22_12 b22_13
## r2-2 b21_21 b21_22 b22_21 b22_22 b22_23
## r2-3 b21_31 b21_32 b22_31 b22_32 b22_33# 転置は要素をスカラーだと思って転置した後に要素を転置
# 要素の行列の転置をスカラーだと思って転置
t(bigA) |> Printb()## t bigA
## r1 r2 r1 r2 r3
## c1 "a11_11" "a11_21" "a21_11" "a21_21" "a21_31"
## c2 "a11_12" "a11_22" "a21_12" "a21_22" "a21_32"
## c1 "a12_11" "a12_21" "a22_11" "a22_21" "a22_31"
## c2 "a12_12" "a12_22" "a22_12" "a22_22" "a22_32"
## c3 "a12_13" "a12_23" "a22_13" "a22_23" "a22_33"tbigA=rbind( cbind(t(A11),t(A21)), cbind(t(A12),t(A22)) )
t(bigA)-tbigA## r1 r2 r1 r2 r3
## c1 "0" "0" "0" "0" "0"
## c2 "0" "0" "0" "0" "0"
## c1 "0" "0" "0" "0" "0"
## c2 "0" "0" "0" "0" "0"
## c3 "0" "0" "0" "0" "0"# 積
bigAB=bigA%*%bigB
# 要素の行列をあたかもスカラーだと思って積を計算
bigAB2=rbind(
cbind( A11%*%B11+A12%*%B21, A11%*%B12+A12%*%B22 )
, cbind( A21%*%B11+A22%*%B21, A21%*%B12+A22%*%B22 )
)
bigAB-bigAB2## c1 c2 c1 c2 c3
## r1 "0" "0" "0" "0" "0"
## r2 "0" "0" "0" "0" "0"
## r1 "0" "0" "0" "0" "0"
## r2 "0" "0" "0" "0" "0"
## r3 "0" "0" "0" "0" "0"# 数値的検算
Eval( bigAB-bigAB2
, A11=demomat(2,2),A12=demomat(2,3), A21=demomat(3,2), A22=demomat(3,3)
, B11=demomat(2,2),B12=demomat(2,3), B21=demomat(3,2), B22=demomat(3,3)
)## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 0 0 0 0
## [2,] 0 0 0 0 0
## [3,] 0 0 0 0 0
## [4,] 0 0 0 0 0
## [5,] 0 0 0 0 0# 積の転置
tbigAB=t(bigAB)
tbigAB2=t(bigB)%*%t(bigA)
tbigAB-tbigAB2## r1 r2 r1
## c1 "a11_11*b11_11+a11_12*b11_21+a12_11*b21_11+a12_12*b21_21+a12_13*b21_31-(b11_11*a11_11+b11_21*a11_12+b21_11*a12_11+b21_21*a12_12+b21_31*a12_13)" "a11_21*b11_11+a11_22*b11_21+a12_21*b21_11+a12_22*b21_21+a12_23*b21_31-(b11_11*a11_21+b11_21*a11_22+b21_11*a12_21+b21_21*a12_22+b21_31*a12_23)" "a21_11*b11_11+a21_12*b11_21+a22_11*b21_11+a22_12*b21_21+a22_13*b21_31-(b11_11*a21_11+b11_21*a21_12+b21_11*a22_11+b21_21*a22_12+b21_31*a22_13)"
## c2 "a11_11*b11_12+a11_12*b11_22+a12_11*b21_12+a12_12*b21_22+a12_13*b21_32-(b11_12*a11_11+b11_22*a11_12+b21_12*a12_11+b21_22*a12_12+b21_32*a12_13)" "a11_21*b11_12+a11_22*b11_22+a12_21*b21_12+a12_22*b21_22+a12_23*b21_32-(b11_12*a11_21+b11_22*a11_22+b21_12*a12_21+b21_22*a12_22+b21_32*a12_23)" "a21_11*b11_12+a21_12*b11_22+a22_11*b21_12+a22_12*b21_22+a22_13*b21_32-(b11_12*a21_11+b11_22*a21_12+b21_12*a22_11+b21_22*a22_12+b21_32*a22_13)"
## c1 "a11_11*b12_11+a11_12*b12_21+a12_11*b22_11+a12_12*b22_21+a12_13*b22_31-(b12_11*a11_11+b12_21*a11_12+b22_11*a12_11+b22_21*a12_12+b22_31*a12_13)" "a11_21*b12_11+a11_22*b12_21+a12_21*b22_11+a12_22*b22_21+a12_23*b22_31-(b12_11*a11_21+b12_21*a11_22+b22_11*a12_21+b22_21*a12_22+b22_31*a12_23)" "a21_11*b12_11+a21_12*b12_21+a22_11*b22_11+a22_12*b22_21+a22_13*b22_31-(b12_11*a21_11+b12_21*a21_12+b22_11*a22_11+b22_21*a22_12+b22_31*a22_13)"
## c2 "a11_11*b12_12+a11_12*b12_22+a12_11*b22_12+a12_12*b22_22+a12_13*b22_32-(b12_12*a11_11+b12_22*a11_12+b22_12*a12_11+b22_22*a12_12+b22_32*a12_13)" "a11_21*b12_12+a11_22*b12_22+a12_21*b22_12+a12_22*b22_22+a12_23*b22_32-(b12_12*a11_21+b12_22*a11_22+b22_12*a12_21+b22_22*a12_22+b22_32*a12_23)" "a21_11*b12_12+a21_12*b12_22+a22_11*b22_12+a22_12*b22_22+a22_13*b22_32-(b12_12*a21_11+b12_22*a21_12+b22_12*a22_11+b22_22*a22_12+b22_32*a22_13)"
## c3 "a11_11*b12_13+a11_12*b12_23+a12_11*b22_13+a12_12*b22_23+a12_13*b22_33-(b12_13*a11_11+b12_23*a11_12+b22_13*a12_11+b22_23*a12_12+b22_33*a12_13)" "a11_21*b12_13+a11_22*b12_23+a12_21*b22_13+a12_22*b22_23+a12_23*b22_33-(b12_13*a11_21+b12_23*a11_22+b22_13*a12_21+b22_23*a12_22+b22_33*a12_23)" "a21_11*b12_13+a21_12*b12_23+a22_11*b22_13+a22_12*b22_23+a22_13*b22_33-(b12_13*a21_11+b12_23*a21_12+b22_13*a22_11+b22_23*a22_12+b22_33*a22_13)"
## r2 r3
## c1 "a21_21*b11_11+a21_22*b11_21+a22_21*b21_11+a22_22*b21_21+a22_23*b21_31-(b11_11*a21_21+b11_21*a21_22+b21_11*a22_21+b21_21*a22_22+b21_31*a22_23)" "a21_31*b11_11+a21_32*b11_21+a22_31*b21_11+a22_32*b21_21+a22_33*b21_31-(b11_11*a21_31+b11_21*a21_32+b21_11*a22_31+b21_21*a22_32+b21_31*a22_33)"
## c2 "a21_21*b11_12+a21_22*b11_22+a22_21*b21_12+a22_22*b21_22+a22_23*b21_32-(b11_12*a21_21+b11_22*a21_22+b21_12*a22_21+b21_22*a22_22+b21_32*a22_23)" "a21_31*b11_12+a21_32*b11_22+a22_31*b21_12+a22_32*b21_22+a22_33*b21_32-(b11_12*a21_31+b11_22*a21_32+b21_12*a22_31+b21_22*a22_32+b21_32*a22_33)"
## c1 "a21_21*b12_11+a21_22*b12_21+a22_21*b22_11+a22_22*b22_21+a22_23*b22_31-(b12_11*a21_21+b12_21*a21_22+b22_11*a22_21+b22_21*a22_22+b22_31*a22_23)" "a21_31*b12_11+a21_32*b12_21+a22_31*b22_11+a22_32*b22_21+a22_33*b22_31-(b12_11*a21_31+b12_21*a21_32+b22_11*a22_31+b22_21*a22_32+b22_31*a22_33)"
## c2 "a21_21*b12_12+a21_22*b12_22+a22_21*b22_12+a22_22*b22_22+a22_23*b22_32-(b12_12*a21_21+b12_22*a21_22+b22_12*a22_21+b22_22*a22_22+b22_32*a22_23)" "a21_31*b12_12+a21_32*b12_22+a22_31*b22_12+a22_32*b22_22+a22_33*b22_32-(b12_12*a21_31+b12_22*a21_32+b22_12*a22_31+b22_22*a22_32+b22_32*a22_33)"
## c3 "a21_21*b12_13+a21_22*b12_23+a22_21*b22_13+a22_22*b22_23+a22_23*b22_33-(b12_13*a21_21+b12_23*a21_22+b22_13*a22_21+b22_23*a22_22+b22_33*a22_23)" "a21_31*b12_13+a21_32*b12_23+a22_31*b22_13+a22_32*b22_23+a22_33*b22_33-(b12_13*a21_31+b12_23*a21_32+b22_13*a22_31+b22_23*a22_32+b22_33*a22_33)"Eval(tbigAB-tbigAB2
, A11=demomat(2,2),A12=demomat(2,3), A21=demomat(3,2), A22=demomat(3,3)
, B11=demomat(2,2),B12=demomat(2,3), B21=demomat(3,2), B22=demomat(3,3)
)## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 0 0 0 0
## [2,] 0 0 0 0 0
## [3,] 0 0 0 0 0
## [4,] 0 0 0 0 0
## [5,] 0 0 0 0 0# ブロック対角行列
bigA=b_diag(A,B)
bigX=b_diag(X,Y)
# 転置
tbigA=b_diag(t(A),t(B))
t(bigA)-tbigA## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] "0" "0" "0" "0" "0" "0" "0" "0"
## [2,] "0" "0" "0" "0" "0" "0" "0" "0"
## [3,] "0" "0" "0" "0" "0" "0" "0" "0"
## [4,] "0" "0" "0" "0" "0" "0" "0" "0"
## [5,] "0" "0" "0" "0" "0" "0" "0" "0"
## [6,] "0" "0" "0" "0" "0" "0" "0" "0"# かけ算
bigAX=b_diag(A %*% X, B %*% Y)
bigA%*%bigX-bigAX## [,1] [,2] [,3] [,4]
## [1,] "0" "0" "0" "0"
## [2,] "0" "0" "0" "0"
## [3,] "0" "0" "0" "0"
## [4,] "0" "0" "0" "0"
## [5,] "0" "0" "0" "0"
## [6,] "0" "0" "0" "0"
## [7,] "0" "0" "0" "0"
## [8,] "0" "0" "0" "0"# set output line width
options(width=500)
# ブロック対角行列の逆行列はブロックごとに逆行列を計算したもの
B2=demomat(2,2,root="b")
bigA=b_diag(A2,B2) |> print()## [,1] [,2] [,3] [,4]
## [1,] "a11" "a12" "0" "0"
## [2,] "a21" "a22" "0" "0"
## [3,] "0" "0" "b11" "b12"
## [4,] "0" "0" "b21" "b22"invBigA=b_diag(matSwp(A2),matSwp(B2)) |> print()## [,1] [,2] [,3] [,4]
## [1,] "1/a11-(a12/a11)*(-a21/a11)/(a22-a21*a12/a11)" "-(a12/a11)/(a22-a21*a12/a11)" "0" "0"
## [2,] "(-a21/a11)/(a22-a21*a12/a11)" "1/(a22-a21*a12/a11)" "0" "0"
## [3,] "0" "0" "1/b11-(b12/b11)*(-b21/b11)/(b22-b21*b12/b11)" "-(b12/b11)/(b22-b21*b12/b11)"
## [4,] "0" "0" "(-b21/b11)/(b22-b21*b12/b11)" "1/(b22-b21*b12/b11)"AinvA=bigA%*%invBigA |> print()## [,1] [,2] [,3] [,4]
## [1,] "a11*(1/a11-(a12/a11)*(-a21/a11)/(a22-a21*a12/a11))+a12*(-a21/a11)/(a22-a21*a12/a11)" "-a11*(a12/a11)/(a22-a21*a12/a11)+a12/(a22-a21*a12/a11)" "0" "0"
## [2,] "a21*(1/a11-(a12/a11)*(-a21/a11)/(a22-a21*a12/a11))+a22*(-a21/a11)/(a22-a21*a12/a11)" "-a21*(a12/a11)/(a22-a21*a12/a11)+a22/(a22-a21*a12/a11)" "0" "0"
## [3,] "0" "0" "b11*(1/b11-(b12/b11)*(-b21/b11)/(b22-b21*b12/b11))+b12*(-b21/b11)/(b22-b21*b12/b11)" "-b11*(b12/b11)/(b22-b21*b12/b11)+b12/(b22-b21*b12/b11)"
## [4,] "0" "0" "b21*(1/b11-(b12/b11)*(-b21/b11)/(b22-b21*b12/b11))+b22*(-b21/b11)/(b22-b21*b12/b11)" "-b21*(b12/b11)/(b22-b21*b12/b11)+b22/(b22-b21*b12/b11)"Eval(AinvA, A2=demomat(2,2), B2=demomat(2,2)+10) |> round(4)## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 1 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 1invbigA=matSwp(bigA) |> print()## [,1] [,2] [,3] [,4]
## [1,] "1/a11-(a12/a11)*(-a21/a11)/(a22-a21*a12/a11)" "-(a12/a11)/(a22-a21*a12/a11)" "0" "0"
## [2,] "(-a21/a11)/(a22-a21*a12/a11)" "1/(a22-a21*a12/a11)" "0" "0"
## [3,] "0" "0" "1/b11-(b12/b11)*(-b21/b11)/(b22-b21*b12/b11)" "-(b12/b11)/(b22-b21*b12/b11)"
## [4,] "0" "0" "(-b21/b11)/(b22-b21*b12/b11)" "1/(b22-b21*b12/b11)"invbigA-invBigA## [,1] [,2] [,3] [,4]
## [1,] "0" "0" "0" "0"
## [2,] "0" "0" "0" "0"
## [3,] "0" "0" "0" "0"
## [4,] "0" "0" "0" "0"# set output line width
options(width=500)
# u=c("f1","z","e1","e2") の共分散構造
Lambda=demomat(2,4,root="A", fullsymb=1)
Lambda[1,]=c("A1","0","I", "0")
Lambda[2,]=c("A2%*%G","A2","0", "I")
Du=demomat(4,4,root="d",fullsymb=1)
Du[]=0
diag(Du)=c("Phi","Sz","Psi1","Psi2")
Sigma=Lambda%*%Du%*%t(Lambda)
Print(Lambda, Du, "\n", Sigma)## c1 c2 c3 c4 c1 c2 c3 c4
## Lambda = r1 A1 0 I 0 , Du = r1 Phi 0 0 0
## r2 A2%*%G A2 0 I r2 0 Sz 0 0
## r3 0 0 Psi1 0
## r4 0 0 0 Psi2
## Sigma =
## r1 r2
## r1 A1%*%Phi%*%t(A1)+Psi1 A1%*%Phi%*%t(A2%*%G)
## r2 A2%*%G%*%Phi%*%t(A1) A2%*%G%*%Phi%*%t(A2%*%G)+A2%*%Sz%*%t(A2)+Psi2
## attr(,"fullsymb")
## [1] 1# set output line width
options(width=500)
# v=c("xi","eta","y")
# u=c("xi","z","eps")
Lambda=demomat(4,2,root="a",fullsymb=1)
Lambda[]=0
Lambda[1,1]="L1"; Lambda[2,2]="L2"
Kappa=demomat(4,2,root="a",fullsymb=1)
Kappa[]=0
Kappa[3,1]="K1"; Kappa[4,2]="K2"
Gamma=demomat(2,2,root="a",fullsymb=1)
Gamma[]=0
Gamma[1,1]="G1"; Gamma[2,1]="G2"; Gamma[2,2]="G3"
Print(Lambda,Kappa,Gamma)## c1 c2 c1 c2 c1 c2
## Lambda = r1 L1 0 , Kappa = r1 0 0 , Gamma = r1 G1 0
## r2 0 L2 r2 0 0 r2 G2 G3
## r3 0 0 r3 K1 0
## r4 0 0 r4 0 K2Bss=demomat(2,2,fullsymb=1)
Bss[]=0
Bss[2,1]="Bs"
Bs=demomat(3,3,fullsymb=1)
Bs[]=0
Bs[2,1]="G"
Bs[2,2]="Bss"
Bs[3,1]="L"
Bs[3,2]="K"
Print(Bss,Bs)## c1 c2 c1 c2 c3
## Bss = r1 0 0 , Bs = r1 0 0 0
## r2 Bs 0 r2 G Bss 0
## r3 L K 0I3=demomat(3,3,shape="diag",fullsymb=1)
diag(I3)="I"
I3mBs=I3-Bs
Print(I3mBs)## I3mBs =
## c1 c2 c3
## r1 I 0 0
## r2 -G I-Bss 0
## r3 -L -K I
## attr(,"fullsymb")
## [1] 1invI3mBs=matSweep(I3mBs)
Print(invI3mBs)## invI3mBs =
## c1 c2 c3
## r1 I 0 0
## r2 inv(I-Bss)%*%G inv(I-Bss) 0
## r3 L+K%*%inv(I-Bss)%*%G K%*%inv(I-Bss) I
## attr(,"fullsymb")
## [1] 1# dispersin of u
Du=demomat(3,3,fullsymb=1)
Du[]=0
Du[1,1]="Phi"
Du[2,2]="Del"
Du[3,3]="Psi"
# dispersion of v
Sigma=invI3mBs%*%Du%*%t(invI3mBs)
S33=Sigma[3,3]
Print(S33)## S33 = (L+K%*%inv(I-Bss)%*%G)%*%Phi%*%t(L+K%*%inv(I-Bss)%*%G)+K%*%inv(I-Bss)%*%Del%*%t(K%*%inv(I-Bss))+Psifind_terms(S33)$tab## oper term
## term_1 "" "(L+K%*%inv(I-Bss)%*%G)%*%Phi%*%t(L+K%*%inv(I-Bss)%*%G)"
## term_2 "+" "K%*%inv(I-Bss)%*%Del%*%t(K%*%inv(I-Bss))"
## term_3 "+" "Psi"Harville, D. A. (1997) Matrix Algebra From a Statistician’s Perspective. Springer.