- メモリンク
- "Matrix.Rnw"と"R_Matrix.Rnw"
\documentclass{jsbook}
\usepackage{amsmath,amssymb}
\usepackage{bm}
\usepackage{ascmac}
\usepackage[dvips]{graphicx}
\usepackage{graphicx}
\usepackage{bigdelim,multirow}
\usepackage{amsmath,amsthm,amssymb,cases}
\usepackage{ascmac}
\usepackage{eclbkbox}
\usepackage{wrapfig}
\usepackage{listings, jlisting}
\usepackage[dvips,usenames]{color}
\usepackage{makeidx}
\usepackage{url}
\usepackage[dvipdfmx,setpagesize=false]{hyperref}
\SweaveOpts{eps=TRUE}
\lstdefinestyle{MyFrame}{backgroundcolor=\color{green},frame=shadowbox}
\lstdefinestyle{MyC++Style} {language=C++,style=MyFrame,frame=none}
\lstset{
language=R,
basicstyle=\ttfamily,
backgroundcolor={\color[gray]{.85}},
keywordstyle=\color{blue}\bfseries,
}
\begin{document}
\chapter{行列}\label{Matrix}
\section{行列とは}
行列は数値を長方形に並べたもの。
$n$行$m$列の行列は$n\times m$行列とも言い、次のように書く。
\begin{equation*}
A_{nm} = \begin{pmatrix}
a_{11} & \ldots & a_{1m} \\
\vdots & \ddots & \vdots \\
a_{n1} & \ldots & a_{nm}
\end{pmatrix}
\end{equation*}
$n,m=1,2,...$であるから
$1\times 1$行列は次の通りで、これは行列とみなしてもよいが、ただの数値とみなしてもよい。
\begin{equation*}
A_{11} = \begin{pmatrix}
a_{11}\\
\end{pmatrix}
\end{equation*}
ただの数値はスカラー値と呼ぶこともある。
行列が行列らしさを発揮するのは行も列も複数のときなので、「行列らしさ」の感じられるもっとも小さい行列は、というと、$2 \times 2$行列となる。
\begin{equation*}
A_{22} = \begin{pmatrix}
a_{11} & a_{11} \\
a_{21} & a_{22}
\end{pmatrix}
\end{equation*}
\section{行列の使い方}
\subsection{影響しあう関係〜推移行列}
２種類の生物X,Yが居て、Xは自然増するがYに捕食され、
YはXが居れば、それを捕食して増加するが、放っておけば減少するとする。
そのような関係が１年ごとに繰り返されると、X,Yの頭数はどう変わるかを、
行列を使ってシミュレーションしてみる。
第i年のX,Yの頭数を$x(i),y(i)$と書くことにする。
ここで$k_x,k_{xy},k_{yx},k_y$は係数で、Xの自然増を表す$k_x$は１より大、
YによるXの捕食を表す$k_{xy}$は負、Xに依存するYの増加を表す$k_{yx}$は正、
Yの自然減を表す$k_{y}$は１より小さい正の数であって、関係式は以下の通り。
\begin{eqnarray*}
x(i+1) &=& k_x \times x(i) + k_{xy} \times y(i)\\
y(i+1) &=& k_{yx} \times x(i) + k_y \times y(i)
\end{eqnarray*}
このとき４個の$k_*$の値を$2 \times 2$行列で表せば
\begin{equation*}
K = \begin{pmatrix}
k_x & x_{xy} \\
k_{yx} & k_y
\end{pmatrix}
\end{equation*}
となり、
\begin{equation*}
\begin{pmatrix}
x(i+1)\\
y(i+1)
\end{pmatrix}
= K
\begin{pmatrix}
x(i)\\
y(i)
\end{pmatrix}
\end{equation*}
Rを使ってシミュレーションしてみる。ある状態に収束する様子が見て取れる。
\begin{lstlisting}
n.year <- 50
xy <- matrix(0,n.year,2)
K <- matrix(c(1.1,-0.2,0.1,0.8),byrow=TRUE,2,2)
xy[1,] <- c(100,5)
for(i in 2:n.year){
xy[i,] <- K
}
matplot(1:n.year,xy,type="l")
\end{lstlisting}
\begin{breakbox}
<<fig=TRUE>>=
n.year <- 50
xy <- matrix(0,n.year,2)
K <- matrix(c(1.1,-0.2,0.1,0.8),byrow=TRUE,2,2)
xy[1,] <- c(100,5)
for(i in 2:n.year){
xy[i,] <- K
}
matplot(1:n.year,xy,type="l")
@
\end{breakbox}
少し変えてみよう。
\subsection{回転・拡縮〜３次元画像処理}
\subsubsection{２次元の回転と拡縮}
適当な形を作ってみる。
\begin{lstlisting}
n.pt <- 10000
t <- seq(from = 0, to = 1, length = n.pt) * 2 * pi
t <- t[-1]
R <- rep(0,length(t))
ss <- rep(1,length(t))
s <- c(0,2*pi,runif(sample(3:8,1)) * 10-5)
for(j in 1:length(s)){
ss <- ss * (t-s[j])^2
}
R <- R + runif(1) * ss
t.2 <- c(t,t[1])
R.2 <- c(R,R[1])
R.2 <- (R.2 - min(R.2))/(max(R.2)-min(R.2))+1
xlim <- ylim <- range(X)
plot(R.2*cos(t.2),R.2*sin(t.2),cex=0.1)
\end{lstlising}
\begin{breakbox}
<<fig=TRUE>>=
n.pt <- 10000
t <- seq(from = 0, to = 1, length = n.pt) * 2 * pi
t <- t[-1]
R <- rep(0,length(t))
ss <- rep(1,length(t))
s <- c(0,2*pi,runif(sample(3:8,1)) * 10-5)
for(j in 1:length(s)){
ss <- ss * (t-s[j])^2
}
R <- R + runif(1) * ss
t.2 <- c(t,t[1])
R.2 <- c(R,R[1])
R.2 <- (R.2 - min(R.2))/(max(R.2)-min(R.2))+1
X <- cbind(R.2*cos(t.2),R.2*sin(t.2))
xlim <- ylim <- range(X)
plot(X,cex = 0.1,xlim=xlim,ylim=ylim)
@
\end{breakbox}
これを回転させるには、
\begin{equation*}
M = \begin{pmatrix}
\cos{\theta} & -\sin{\theta} \\
\sin{\theta} & \cos{\theta}
\end{pmatrix}
\end{equation*}
を各点の座標にかけてやればよい。
\begin{equation*}
M v
\end{equation*}
\begin{lstlisting}
theta <- pi/3
M <- matrix(c(cos(theta),-sin(theta),sin(theta),cos(theta)),byrow=TRUE,2,2)
X.rot <- t(M
xlim <- ylim <- range(c(X,X.rot))
plot(X,type="l",cex = 0.1,xlim=xlim,ylim=ylim)
points(X.rot,cex=0.1,col=2)
\end{lstlisting}
\begin{breakbox}
<<fig=TRUE>>=
theta <- pi/3
M <- matrix(c(cos(theta),-sin(theta),sin(theta),cos(theta)),byrow=TRUE,2,2)
X.rot <- t(M
xlim <- ylim <- range(c(X,X.rot))
plot(X,cex = 0.1,xlim=xlim,ylim=ylim)
points(X.rot,cex=0.1,col=2)
@
\end{breakbox}
これをx,y軸方向に$Lx,Ly$倍するには
\begin{equation*}
L = \begin{pmatrix}
Lx & 0 \\
0 & Ly
\end{pmatrix}
\end{equation*}
を各点の座標にかけてやればよい。
\begin{equation*}
L v
\end{equation*}
\begin{lstlisting}
Lx <- 2
Ly <- 0.5
L <- matrix(c(Lx,0,0,Ly),byrow=TRUE,2,2)
X.mult <- t(L
xlim <- ylim <- range(c(X,X.rot,X.mult))
plot(X,cex = 0.1,xlim=xlim,ylim=ylim)
points(X.rot,cex=0.1,col=2)
points(X.mult,cex=0.1,col=3)
\end{lstlisting}
\begin{breakbox}
<<fig=TRUE>>=
Lx <- 2
Ly <- 0.5
L <- matrix(c(Lx,0,0,Ly),byrow=TRUE,2,2)
X.mult <- t(L
xlim <- ylim <- range(c(X,X.rot,X.mult))
plot(X,type="l",cex = 0.1,xlim=xlim,ylim=ylim)
points(X.rot,cex=0.1,col=2)
points(X.mult,cex=0.1,col=3)
@
\end{breakbox}
\subsubsection{行列をかける順番について}
では、次の２つを試してみる。
\begin{itemize}
\item 回転してから、x,y軸方向に伸縮する
\item x,y軸方向に伸縮してから回転させる
\end{itemize}
それぞれ
\begin{eqnarray*}
LM v
ML v
\end{eqnarray*}
となる。
\begin{lstlisting}
X.rot.mult <- t((L
X.mult.rot <- t((M
xlim <- ylim <- range(c(X.rot.mult,X.mult.rot))
plot(X.rot.mult,type="l",cex = 0.1,xlim=xlim,ylim=ylim)
points(X.mult.rot,cex=0.1,col=2)
\end{lstlisting}
\begin{breakbox}
<<fig=TRUE>>=
X.rot.mult <- t((L
X.mult.rot <- t((M
xlim <- ylim <- range(c(X.rot.mult,X.mult.rot))
plot(X.rot.mult,type="l",cex = 0.1,xlim=xlim,ylim=ylim)
points(X.mult.rot,cex=0.1,col=2)
@
\end{breakbox}
２つの動きが異なることがわかる。
行列の掛け算は、$LM$と$ML$とで違うことも意味している。
\subsection{３次元物体の回転}
３次元物体をいろいろな方向から眺めたいことはよくある。
人体画像がそのよい例である。
人体のデータが３次元座標で得られているので、正面から見てみよう。
横からも見て、上からも見てみる。
\begin{lstlisting}
perfumev <- as.matrix(read.table("perfume_v.txt"))
perfumev.2 <- rbind(perfumev,rep(max(perfumev),3))
perfumev.2 <- rbind(perfumev.2,rep(min(perfumev),3))
par(mfcol=c(1,3))
plot(perfumev.2[,c(1,2)],cex=0.01)
plot(perfumev.2[,c(3,2)],cex=0.01)
plot(perfumev.2[,c(1,3)],cex=0.01)
par(mfcol=c(1,1))
\end{lstlisting}
\begin{breakbox}
<<fig=TRUE>>=
perfumev <- as.matrix(read.table("perfume_v.txt"))
perfumev.2 <- rbind(perfumev,rep(max(perfumev),3))
perfumev.2 <- rbind(perfumev.2,rep(min(perfumev),3))
par(mfcol=c(1,3))
plot(perfumev.2[,c(1,2)],cex=0.01)
plot(perfumev.2[,c(3,2)],cex=0.01)
plot(perfumev.2[,c(1,3)],cex=0.01)
par(mfcol=c(1,1))
@
\end{breakbox}
このように３方向から眺められることは素晴らしいが、任意の角度から眺められたらさらによい。
３次元空間での回転を定義して、それを実行しよう。
v1-v3平面で回る(右端の表示での回転)は、頭を上に足を下にしたまま、背中を見せる回転である。
そのような回転をしたうえで、上から覗き込むように回転を加えよう。v2-v3平面で回る(左端の表示での回転)ことになる。
それぞれの回転は、３次元空間の３軸のうち２軸だけを取り出して、それについて、２次元回転をすることだから、次のように書ける。
\begin{equation*}
M1 = \begin{pmatrix}
\cos{\theta} & 0 & -\sin{\theta} \\
0 & 1 & 0 \\
\sin{\theta} & 0 & \cos{\theta}
\end{pmatrix}
\end{equation*}
\begin{equation*}
M2 = \begin{pmatrix}
1 & 0 & 0 \\
0 & \cos{\theta2} & -\sin{\theta2} \\
0 & \sin{\theta2} & \cos{\theta2}
\end{pmatrix}
\end{equation*}
その上で２つの回転をすればよい。
\begin{equation*}
M2 M1
\end{equation*}
\begin{lstlisting}
theta1 <- pi/3
theta2 <- pi/4
M1 <- matrix(c(cos(theta1),0,-sin(theta1),0,1,0,sin(theta1),0,cos(theta1)),byrow=TRUE,3,3)
M2 <- matrix(c(1,0,0,0,cos(theta1),-sin(theta1),0,sin(theta1),cos(theta1)),byrow=TRUE,3,3)
M12 <- M2
perfumev.2.rot <- t(M12
par(mfcol=c(1,3))
plot(perfumev.2.rot[,c(1,2)],cex=0.01)
plot(perfumev.2.rot[,c(3,2)],cex=0.01)
plot(perfumev.2.rot[,c(1,3)],cex=0.01)
par(mfcol=c(1,1))
\end{lstlisting}
\begin{breakbox}
<<fig=TRUE>>=
theta1 <- pi/3
theta2 <- pi/4
M1 <- matrix(c(cos(theta1),0,-sin(theta1),0,1,0,sin(theta1),0,cos(theta1)),byrow=TRUE,3,3)
M2 <- matrix(c(1,0,0,0,cos(theta1),-sin(theta1),0,sin(theta1),cos(theta1)),byrow=TRUE,3,3)
M12 <- M2
perfumev.2.rot <- t(M12
par(mfcol=c(1,3))
plot(perfumev.2.rot[,c(1,2)],cex=0.01)
plot(perfumev.2.rot[,c(3,2)],cex=0.01)
plot(perfumev.2.rot[,c(1,3)],cex=0.01)
par(mfcol=c(1,1))
@
\end{breakbox}
\subsubsection{回転とは、そして、回転行列の特徴}
２次元、３次元の回転行列を扱ってきた。
任意次元の回転行列も２次元の回転の繰り返しで作ることができる。
\begin{lstlisting}
n <- 5
R <- diag(rep(1,n))
for(i in 1:n){
a <- sample(1:n,2)
theta <- runif(1)*2*pi
tmp <- diag(rep(1,n))
tmp[a,a] <- matrix(c(cos(theta),-sin(theta),sin(theta),cos(theta)),byrow=TRUE,2,2)
R <- tmp
}
\end{lstlisting}
２次元の回転や３次元の回転は、像が「きちんと回転している」ということから、使用した
回転行列が「きちんと回転を表している」ことを実感できた。
ではここで作成した例のように高次元になったとき、行列が回転を表しているかどうかはどうやって確認すればよいだろうか。
回転では基底も回転する。
基底のうち、一般によく使われるのは正規直交基底。これは、２次元空間なら、(1,0)と(0,1)の２つ、３次元空間なら、(1,0,0),(0,1,0)と(0,0,1)の３つ、というように、各軸方向の単位ベクトルを合わせたもの。
これらはすべてが単位ベクトルで、どの２つも相互に直交する。
第i要素が1である単位ベクトルに行列を作用させると、行列の第i列が取り出される。
したがって、n個の単位行列からなる正規直交基底に行列を作用させると、正規直交基底が、その行列に移る。
今、回転では、ベクトルの長さも２つのベクトルが作る角度も変わらないから、移った先の基底のベクトルはすべて単位ベクトルで、それらは相互に直交しているはずである。
このような性質を行列の特徴で言い直すと、
\begin{itemize}
\item 直交行列である
\item 行列式が1である
\end{itemize}
直交行列とは、行列自身とその転置行列($a_{ij}$の値と$a_{ji}$の値を入れ替えた行列)との積が単位行列になるような(正方)行列のことであり、行列式とは$2 \times 2$行列で言えば
\begin{equation*}
a_{11}a_{22}-a_{12}a_{21}
\end{equation*}
と計算される行列を特徴づける量である。
上で作った回転行列と目される行列がこの条件を満足しているかどうかを確認してみます。
\begin{lstlisting}
R_t <- t(R)
R
R_t
det(R)
det(R_t)
\end{lstlisting}
\begin{breakbox}
<<fig=FALSE>>=
n <- 5
R <- diag(rep(1,n))
for(i in 1:n){
a <- sample(1:n,2)
theta <- runif(1)*2*pi
tmp <- diag(rep(1,n))
tmp[a,a] <- matrix(c(cos(theta),-sin(theta),sin(theta),cos(theta)),byrow=TRUE,2,2)
R <- tmp
}
R_t <- t(R)
R
R_t
det(R)
det(R_t)
@
\end{breakbox}
\subsubsection{データの特徴を抽出する}
次のような３次元データを考えよう。
３人が立っている。
どの方向から見ると、一番「よく見える」だろうか。
３次元の対象を見るとき、３次元のうち２次元の絵を作って、残りの
１次元はつぶしてしまう。
したがって、２次元を選んだときに一番絵がばらついて見えるときが
「もっともよく見える」状態だと考えるのは、よいやり方の一つだろうと思える。
図のように点がばらついているとき、３つの互いに垂直な軸をとるとする。
このようなとき、第１の軸として、とりえるあらゆる軸のうち、分散が最大になるものを採用し、
ついで、第２の軸をとることにする。そのとき、第１の軸に垂直であることを条件にして、分散が最大になるようにする。多次元の場合には、順次、それまでに採用された軸に垂直であることを条件として順番に選んでいくことにする。
このようにすると、全部の軸について見てみなくても、上位の軸だけで結構、全体像が見えてくる可能性がある。
では、アイディアとして、この方法がよいと思われるが、それを実際にはどうしたらよいか、ということが問題になる。
そのときにデータ点の行列$D$をあるルール(特異値分解という方法)で分解して
\begin{equation*}
D = V \Sigma U^T
\end{equation*}
としてやる。
そうすると、回転行列$U$を使って$D$を回転させた
\begin{equation*}
DU = V \Sigma
\end{equation*}
は、各軸の分散が上で述べたようなルールに従うことが知られている。
\begin{breakbox}
<<fig=TRUE>>=
# ３人の座標を取り込む
#perfumev <- as.matrix(read.table("perfume_v.txt"))
#perfume2v <- as.matrix(read.table("perfume_v2.txt"))
# ３人の立ち位置を少し変える
#perfume2v[,1] <- perfume2v[,1] +100
#perfume2v[,3] <- perfume2v[,3] + 200
#perfume3v <- as.matrix(read.table("perfume_v3.txt"))
#perfume3v[,1] <- perfume3v[,1] + 200
#perfume3v[,3] <- perfume3v[,3] - 200
# ３人の座標を一つの行列にする
#perfume.trio <- rbind(perfumev,perfume2v,perfume3v)
# 回転行列を適当に作って、３人の座標を回す
library(GPArotation)
R <- Random.Start(3)
perfume.trio <- t(R
# 点の数を少なめに取り直す
D <- perfume.trio[sample(1:length(perfume.trio[,1]),5000),]
# ３軸を平均が0になるように調整する
D <- D-mean(D)
mu <- apply(D,2,mean)
D <- t(t(D)-mu)
# 特異値分解する
svdout <-svd(D)
# US = DV のどちらでも回転後の座標が得られるので、２方法で求めておく
D2 <- svdout$u
D3 <- D
par(mfcol=c(3,3))
plot(D[,c(1,2)],cex=0.01)
plot(D[,c(3,2)],cex=0.01)
plot(D[,c(1,3)],cex=0.01)
plot(D2[,c(1,2)],cex=0.01)
plot(D2[,c(3,2)],cex=0.01)
plot(D2[,c(1,3)],cex=0.01)
plot(D3[,c(1,2)],cex=0.01)
plot(D3[,c(3,2)],cex=0.01)
plot(D3[,c(1,3)],cex=0.01)
par(mfcol=c(1,1))
apply(D2,2,var)
apply(D3,2,var)
@
\end{breakbox}
\subsubsection{状態の変化を表す推移行列}
２種類のものA1,A2があって、A1はA2の量に比例して増え、逆にA2はA1の量に比例して減るとする。
\begin{eqnarray*}
\delta a_1 = k a_2\\
\delta a_2 = - k a_1
\end{eqnarray*}
これは時刻$t$の量$a_1(t),a_2(t)$として表すと
\begin{eqnarray*}
a(t+\delta t) = a_1 + k a_2\\
a(t+\delta t) = - k a_1 + a_2
\end{eqnarray*}
これを$(a_1(t),a_2(t))^T$というベクトルで考えれば
\begin{equation*}
\begin{pmatrix}
a_1(t+\delta t)\\
a_2(t+\delta t)
\end{pmatrix}
=
\begin{pmatrix}
1 & k \\
-k & 1
\end{pmatrix}
\begin{pmatrix}
a_1(t)\\
a_2(t)
\end{pmatrix}
\end{equation*}
このように状態をベクトルで表すとき、行列がベクトルの時間変化を表すことができる場合も多い。
Rでやってみよう。
\begin{lstlisting}
k <- 0.001
A <- matrix(c(1,k,-k,1),byrow=TRUE,2,2)
n.t <- 10000
a <- matrix(0,n.t,2)
a[1,] <- c(1,0)
for(i in 2:n.t){
a[i,] <- A
}
matplot(a,type="l")
\end{lstlisting}
\begin{breakbox}
<<fig=TRUE>>=
k <- 0.001
A <- matrix(c(1,k,-k,1),byrow=TRUE,2,2)
n.t <- 10000
a <- matrix(0,n.t,2)
a[1,] <- c(1,0)
for(i in 2:n.t){
a[i,] <- A
}
matplot(a,type="l")
@
\end{breakbox}
要素数は２個に限る必要もなく、細胞内の化学反応系や生態系の捕食関係など、
この枠組みでとらえられる現象は非常に多い。
また、今回のシミュレーションは、短い時間を離散的に計算しているが、連続的な計算をすることもできる(指数行列)。
また、増減が複数の要素の量の積の影響を受けるなどの場合も推移行列で表現できるが、この章では扱わないこととする。
\subsection{行列は便利}
行列は基本的には以下の便利さがあるので、数値をたくさん使う場面(数値を扱う学術分野のすべて)で使われる。
便利さをまとめると以下の通り。
\begin{itemize}
\item $n\times m$個の数字($a_{ij}$)が必要なときに、１個の変数($A_{nm}$)を書くだけで済ませられる
\item たくさんの数字を使った計算を一気に実行できる。特に行列計算機能を持った計算機を持っていれば、楽。
\item 『行列がこれこれの特徴(P)を持っていたら、それそれの特徴(Q)を持つ、と言える』という便利な知識がたくさんある。したがって、Pを確かめる計算が簡単で、Qを確かめる計算が面倒臭いとき、Pをすればよい
\end{itemize}
\section{キーポイント}
\begin{itemize}
\item 行列の和と積。積は順序依存
\item 回転行列
\item 単位行列・転置行列・逆行列
\item 直交
\item 正規直交基底
\item 行列式 determinant
\item 特異値分解
\item 楕円と特異値分解
\item 推移行列
\item 指数行列
\end{itemize}
\end{document}
\documentclass{jsbook}
\usepackage{amsmath,amssymb}
\usepackage{bm}
%\usepackage{graphicx}
\usepackage{ascmac}
\usepackage[dvips]{graphicx}
\usepackage{graphicx}
\usepackage{bigdelim,multirow}
\usepackage{amsmath,amsthm,amssymb,cases}
\usepackage{ascmac}
\usepackage{url}
\usepackage{eclbkbox}
\usepackage{wrapfig}
\usepackage{listings, jlisting}
\usepackage[dvips,usenames]{color}
\usepackage{makeidx}
\SweaveOpts{eps=TRUE}
\lstdefinestyle{MyFrame}{backgroundcolor=\color{green},frame=shadowbox}
\lstdefinestyle{MyC++Style} {language=C++,style=MyFrame,frame=none}
\lstset{
language=R,
basicstyle=\ttfamily,
backgroundcolor={\color[gray]{.85}},
keywordstyle=\color{blue}\bfseries,
}
\begin{document}
\chapter{Rで行列を勉強する}
\section{Matrixオブジェクトの作成}
\subsection{行列を作る}
\begin{lstlisting}
n <- 2
m <- 3
v <- 1:(n*m)
M <- matrix(v,n,m)
M
M2 <- matrix(v,n,m,byrow=TRUE)
M2
\end{lstlisting}
\begin{breakbox}
<<fig=FALSE>>=
n <- 2
m <- 3
v <- 1:(n*m)
M <- matrix(v,n,m)
M
M2 <- matrix(v,n,m,byrow=TRUE)
M2
@
\end{breakbox}
\subsection{行列の四則演算〜数値を納めた箱として扱う}
同じ形の行列は数値を納めた箱として扱える。
\begin{lstlisting}
n <- 2
m <- 3
M1 <- matrix(1:(n*m),n,m)
M2 <- matrix(rep(1:2,(n*m)/2),n,m)
M1
M2
M1 + M2
M1 - M2
M1 * M2
M1 / M2
M1 ^ M2
\end{lstlisting}
\begin{breakbox}
<<fig=FALSE>>=
n <- 2
m <- 3
M1 <- matrix(1:(n*m),n,m)
M2 <- matrix(rep(1:2,(n*m)/2),n,m)
M1
M2
M1 + M2
M1 - M2
M1 * M2
M1 / M2
M1 ^ M2
@
\end{breakbox}
\section{行列の『積』}
\subsection{行列の『いわゆる積』}
$n\times m$ 行列$M_1$を $m \times n$ の行列$M_2$に左からかけることができる。
行列の積$M_1 M_2$はRでは
\begin{lstlisting}
M1 %*% M2
\end{lstlisting}
と書かれこれは、$M1$の第$i$行ベクトルと$M2$の第$j$列ベクトルの内積が
$M_1 M_2$の$(i,j)$成分となるような計算のこと。
\begin{lstlisting}
n <- 2
m <- 3
M1 <- matrix(1:(n*m),n,m)
M2 <- matrix(rep(1:2,(n*m)/2),m,n)
M1
M2
M1 %*% M2
\end{lstlisting}
\begin{breakbox}
<<fig=FALSE>>=
n <- 2
m <- 3
M1 <- matrix(1:(n*m),n,m)
M2 <- matrix(rep(1:2,(n*m)/2),m,n)
M1
M2
M1 %*% M2
@
\end{breakbox}
\subsection{行列の『いろいろな積』}
行列の『積』を理解するために２つの道を取る。
\begin{itemize}
\item 行列をベクトルの束と見て、ベクトルの『積』を考える
\item 行列の『積』を『行列の要素の総当たりの積』から考える
\end{itemize}
\subsubsection{ベクトルとベクトルの『積』}
長さ$n_v$のベクトル$V$と長さ$n_u$のベクトル$U$とを考える。
これらの要素の総当たりペアをとると、$n_v \times n_u$個の値が得られる。
これを$n_v\times n_u$行列として表現することはごく自然である。
($V$を列ベクトル、$U$を行ベクトルとみなしたときの行列の『いわゆる積』
ともみなせる)。
\begin{equation*}
V U^T = \begin{pmatrix}
v_1 u_1 & v_1 u_2 & ... & v_1 u_{n_u}\\
v_2 u_1 & v_2 u_2 & ... & v_2 u_{n_u}\\
... & ... & ... & ... \\
v_{n_v} u_1 & v_{n_v} u_2 & ... & v_{n_v} u_{n_u}
\end{pmatrix}
\end{equation*}
\begin{lstlisting}
nv <- 3
nu <- 4
V <- 1:nv
U <- 1:nu
V %*% t(U)
\end{lstlisting}
\begin{breakbox}
<<fig=FALSE>>=
nv <- 3
nu <- 4
V <- 1:nv
U <- 1:nu
V %*% t(U)
@
\end{breakbox}
これを一番「基本的なベクトルとベクトルの積」としよう。\\
長さが等しい２つのベクトルの場合には、「基本的なベクトルとベクトルの積」は正方行列になる。
その対角成分だけを取り出して和を取ったものが、ベクトルの内積
\begin{breakbox}
<<fig=FALSE>>=
nv <- 3
nu <- nv
V <- 1:nv
U <- 1:nu
W <- V %*% t(U)
diag(W)
sum(diag(W))
sum(V*U)
@
\end{breakbox}
\subsubsection{３次元ベクトルの外積}
３次元ベクトルの間には「ベクトル積」という『積』があって、これは
\begin{equation*}
V U^T = \begin{pmatrix}
v_1 u_1 & v_1 u_2 & v_1 u_3\\
v_2 u_1 & v_2 u_2 & v_2 u_3\\
v_3 u_1 & v_3 u_2 & v_3 u_3
\end{pmatrix}
\end{equation*}
というベクトルの積の９成分のうち、
内積に使われた対角の３成分を使わず、
残りの６成分を、対称な２成分のトリオに分け、結果として、
３つの数字のベクトルにしたものである。
元のベクトルが３次元ベクトルであったが、２つの３次元ベクトルの『積』が
３次元ベクトルになっているので、「うまく閉じた関係」になっていることがわかる。
この『積』である３次元ベクトルには、有用な性質があることから、
名前を「ベクトル積・クロス積」と呼んでいる。
\begin{lstlisting}
my.cross.prod.3d <- function(v,u){
c(v[2]*u[3]-v[3]*u[2], v[3]*u[1]-v[1]*u[3], v[1]*u[2]-v[2]*u[1])
}
\end{lstlisting}
\subsubsection{行列の要素の総当たりの積}
$n_1 \times m_1$ 行列 $M_1$と$n_2 \times m_2$ 行列 $M_2$ との積の計算の仕方は一つではない。
$M_1$行列には$n_1 \times m_1$個の要素があり、$M_2$行列には$n_2 \times m_2$個の要素があるので、
それらの総当たりの積が$(n_1 \times m_1) \times (n_2 \times m_2)$個ある。
これを使って、複数の『行列の積』が定義できる。
この$(n_1 \times m_1) \times (n_2 \times m_2)$個の要素を
$(n_1 \times m_1)$行$\times \times (n_2 \times m_2)$列の行列にするのはクロネッカー積。
\begin{lstlisting}
n1 <- 2
m1 <- 2
n2 <- 3
m2 <-2
M1 <- matrix(1:(n1*n2),n1,m1)
M2 <- matrix(1:(n2*m2),n2,m2)
M1
M2
Kr<- M1 %x% M2
Kr
\end{lstlisting}
\begin{breakbox}
<<fig=FALSE>>=
n1 <- 2
m1 <- 2
n2 <- 3
m2 <-2
M1 <- matrix(1:(n1*n2),n1,m1)
M2 <- matrix(1:(n2*m2),n2,m2)
M1
M2
Kr<- M1 %x% M2
Kr
@
\end{breakbox}
$M_1[i_1,j_2] \times M_2[i_2,j_2]$がクロネッカー積の何行目・何列目になっているかを目視で確認してみるとよい。\\
$(n_1 \times m_1) \times (n_2 \times m_2)$個の要素を
$n_1 \times m_1 \times n_2 \times m_2$と４次元の整理箱にしまう(Rならarrayオブジェクト)のが行列の外積。
これは、$M_1$の行ベクトルと$M_2$の列ベクトルの間のベクトル積を$n_1 \times m_2$個求め、それを$n_1 \times m_2$の行列状に並べたもの。ただし、$M_1$の行ベクトルと$M_2$の列ベクトルの間のベクトル積自体が$m_1 \times n_2$行列なので、出来上がったものは$n_1 \times m_1 \times n_2 \times m_2$と４次元になっている。
Rであれば以下のようにする。
その要素がクロネッカー積の要素と一致することも示しておく(証明ではなく、適当に作ってもいつも要素の値の集合が一致することを示す)。
\begin{lstlisting}
Vpr<- M1 %o% M2
Vpr
range(sort(Kr) - sort(Vpr))
\end{lstlisting}
\begin{lstlisting}
k <- 100
diffs <- rep(0,k)
for(i in 1:k){
n1 <- sample(2:6,1)
m1 <- sample(2:6,1)
n2 <- sample(2:6,1)
m2 <- sample(2:6,1)
M1 <- matrix(runif(n1*m1),n1,m1)
M2 <- matrix(runif(n2*m2),n2,m2)
tmp1 <- M1 %o% M2
tmp2 <- M1 %x% M2
diffs[i] <- max(abs(sort(c(tmp1))-sort(c(tmp2))))
}
range(diffs)
\end{lstlisting}
\begin{breakbox}
<<fig=FALSE>>=
Vpr<- M1 %o% M2
Vpr
range(sort(Kr) - sort(Vpr))
@
\end{breakbox}
\begin{breakbox}
<<fig=FALSE>>=
k <- 100
diffs <- rep(0,k)
for(i in 1:k){
n1 <- sample(2:6,1)
m1 <- sample(2:6,1)
n2 <- sample(2:6,1)
m2 <- sample(2:6,1)
M1 <- matrix(runif(n1*m1),n1,m1)
M2 <- matrix(runif(n2*m2),n2,m2)
tmp1 <- M1 %o% M2
tmp2 <- M1 %x% M2
diffs[i] <- max(abs(sort(c(tmp1))-sort(c(tmp2))))
}
range(diffs)
@
\end{breakbox}
\begin{itemize}
\item $m_1 = n_2$ のとき
\end{itemize}
\section{特徴的な行列}
\section{行列の色々}
\begin{lstlisting}
A <- matrix(1:6,2,3)
t(A)
B <- matrix(1:16,4,4)
B
diag(B)
D <- diag(1:4)
D
D. <- diag(4)
D.
D.. <- diag(rep(1,4))
D..
\end{lstlisting}
連立方程式を解く。
\begin{equation*}
A x = b
\end{equation*}
\begin{lstlisting}
A <- matrix(sample(1:9),3,3)
b <- c(2,3,5)
result <- solve(A,b)
result
A %*% result
F <- matrix(rnorm(9),3,3)
F.inv <- solve(F)
F %*% F.inv
F.inv %*% F
eigen.out <- eigen(F)
V <- eigen.out$vectors
S <- diag(eigen.out$values)
V %*% S %*% solve(V)
V%*%t(V)
det(Re(V))
F. <- F + t(F)
F.
eigen.out.2 <- eigen(F.)
V <- eigen.out.2$vectors
S <- diag(eigen.out.2$values)
V %*% S %*% t(V)
V%*%t(V)
det(Re(V))
G <- matrix(rnorm(12),3,4)
G
svd.out <- svd(G)
svd.out$v %*% diag(svd.out$d) %*% svd.out$u
svd.out$u %*% t(svd.out$u)
\end{lstlisting}
\begin{breakbox}
<<fig=FALSE>>=
A <- matrix(1:6,2,3)
t(A)
B <- matrix(1:16,4,4)
B
diag(B)
D <- diag(1:4)
D
D. <- diag(4)
D.
D.. <- diag(rep(1,4))
D..
A <- matrix(sample(1:9),3,3)
b <- c(2,3,5)
result <- solve(A,b)
result
A %*% result
F <- matrix(rnorm(9),3,3)
F.inv <- solve(F)
F %*% F.inv
F.inv %*% F
eigen.out <- eigen(F)
V <- eigen.out$vectors
S <- diag(eigen.out$values)
V %*% S %*% solve(V)
V%*%t(V)
det(Re(V))
F. <- F + t(F)
F.
eigen.out.2 <- eigen(F.)
V <- eigen.out.2$vectors
S <- diag(eigen.out.2$values)
V %*% S %*% t(V)
V%*%t(V)
det(Re(V))
G <- matrix(rnorm(12),3,4)
G
svd.out <- svd(G)
svd.out$v %*% diag(svd.out$d) %*% svd.out$u
svd.out$u %*% t(svd.out$u)
@
\end{breakbox}
\end{document}
