- 正方行列は一次独立な線形変換することで、単純な形にすることができる
- ここで言う単純な形とは、対角成分付近にのみ値があって、それ以外の成分が0であるような行列のことである
- 実際、線形変換によって、対角成分勝ちにしたときに、その対角成分勝ちな成分が同一にできるような行列は、相互に「本質的に同じ・座標変換すれば同じ」とみなされることから、「相似な行列」と言われる
- では、その「対角成分勝ち」に線形変換することが大事なわけだが、それに関する話である
- 「対角成分勝ち」には大きく二つの「対角成分勝ちな状態」がある
- １つは対角成分のみがあるような場合で、対角行列と呼ばれる
- もう１つは、対角成分勝ちではあるが、非対角成分もあるような行列である。非対角成分もある、とは言うが、たとえ、完全な対角化ができないような行列であっても、１段だけオフダイアゴナルなところに１が立つような行列にまでは「対角成分勝ち」にできることが知られていて、その形をジョルダン標準形と言う
- 言い方を変えると
- すべての正方行列は、「ジョルダン標準形」という相当しっかりと「対角成分勝ち」にすることができるが、ある特定の正方行列の中には、オフダイアゴナル成分に１が立たないようなジョルダン標準形になるものがある
- 対角化できる場合というのは、ジョルダン標準形の特殊形なので、ジョルダン標準形がどういうものか、から考える
- $n\times n$ 正方行列を固有値分解すると、解の重複を考慮すればいつも$n$個の固有値がある
- ジョルダン標準形では、重複を考慮した$n$個の固有値が対角成分になっている
- オフダイアゴナルな１がどこに立つか、と言えば、重複している固有値の場合は、その重複のされ方に応じて１が立つ
- 今、$n$個の固有値がすべて異なれば、すべての固有値は重複していない(重複度が1なので、そのような場合には、オフダイアゴナルに１が立つことはない
- このことから、きれいに対角化できる、というジョルダン標準形の特殊な場合は、重複のない固有値を持つ場合に相当することもある、また、重複のある固有値のジョルダン細胞が、1つずつに分割される、という特別な場合でも、「きれいに対角化」される
- ある固有値が重複度kだと、$k\times k$行列であって、その対角成分は、すべてその固有値で、オフダイアゴナルに1がきれいに並ぶような、そんな$k\times k$行列が、ジョルダン標準形に現れる(こともある)
- 「現れる『こともある』」というのは、重複度kのときに$k\times k$のジョルダン標準形の部品(ジョルダン細胞)が立つ場合もあれば、重複度kがいくつかの小さ目のジョルダン細胞に分かれることもあるために、そのように記載した。ある固有値に対して、サイズがいくつのジョルダン細胞が何個立つかについては、こちらのサイトにあるような関係式があるらしい
- そんなこんなをRの関数にしてみる
- 例題的な行列ではうまく行くが、任意の行列となると、なにかおかしい。ジョルダンブロックの数の特定に失敗している模様
- だが、まあ、目的はジョルダン標準形ってなんだっけ、がわかることだったので、不備のある関数でも、まあ、良しとする
- 参考にしたサイト
my.jordan <- function(A){
n <- dim(A)[1]
eigen.out <- eigen(A)
lambdas <- eigen.out[[1]]
lambdas.re <- Re(lambdas)
lambdas.im <- round(Im(lambdas),6)
lambdas2 <- lambdas.re + 1i * lambdas.im
unique.lambdas <- unique(lambdas2)
num.j.block <- matrix(0,length(unique.lambdas),n)
for(i in 1:length(unique.lambdas)){
Ai <- A - unique.lambdas[i]*diag(rep(1,n))
ranks <- rep(0,n+2)
tmpAi <- diag(rep(1,n))
for(j in 1:(n+2)){
ranks[j] <- rankMatrix(tmpAi)
tmpAi <- Ai %*% tmpAi
}
for(j in 1:n){
num.j.block[i,j] <- ranks[j+2] - 2 * ranks[j+1] + ranks[j]
}
}
J <- matrix(0,n,n)
st <- 1
for(i in 1:length(unique.lambdas)){
ones <- which(num.j.block[i,]>0)
for(j in 1:length(ones)){
jordan.cell <- matrix(0,ones[j],ones[j])
diag(jordan.cell) <- unique.lambdas[i]
if(ones[j]>1){
for(k in 1:(ones[j]-1)){
jordan.cell[k,k+1] <- 1
}
}
for(k in 1:num.j.block[i,ones[j]]){
J[st:(st+ones[j]-1),st:(st+ones[j]-1)] <- jordan.cell
st <- st+ones[j]
}
}
}
P <- matrix(0,n,n)
I <- diag(rep(1,n))
cnt <- 1
for(i in 1:n){
b <- rep(0,n)
pure.eigen <- TRUE
if(i>1){
if(J[i-1,i]==1){
b <- P[i-1,]
pure.eigen <- FALSE
}else{
cnt <- 1
}
}
if(pure.eigen){
tmp <- which(lambdas2==J[i,i])[1]
P[,i] <- eigen.out[[2]][,tmp]
}else{
cnt <- cnt+1
tmpA <- (A-J[i,i]*I) %^% cnt
tmpNull <- Null(t(tmpA))
dim <- length(tmpNull[1,])
coef <- solve(((A-J[i,i]*I) %*% tmpNull)[1:dim,],P[1:dim,i-1])
P[,i] <- tmpNull %*% coef
}
}
return(list(J=J,P=P,num.j.block=num.j.block))
}
A <- matrix(sample(1:50,16),4,4)
A <- matrix(c(1,3,-4,1,0,3,1,-3,6),3,3)
A <- matrix(c(5,0,-1,1,4,1,-1,1,2,-1,3,-1,1,-1,0,2),4,4)
library(Matrix)
library(expm)
library(MASS)
out <- my.jordan(A)
P <- out$P
J <- out$J
A %*% P
P %*% J
A <- matrix(c(1,-1,2,-3,1,4,-2,-1,5),3,3)
my.jordan(A)
A <- matrix(c(1,-1,2,-3,1,4,-2,-1,4),3,3)
my.jordan(A)
A <- matrix(sample(1:40,4^2,replace=TRUE),4,4)
my.jordan(A)
