L <- list(matrix(seq(4),4,4),matrix(seq(9),3,3),matrix(seq(25),5,5))
library(Matrix)
bdiag(L)
M<- as.matrix(bdiag(L))
N<-length(M[,1])
M2<-M[sample(1:N),sample(1:N)]
M3<-sign(M2)
par(ask=TRUE)
heatmap(M)
heatmap(M2)
heatmap(M3)
- PDF
- ブロック対角化、無理やり
- ブロック対角化できる行・列の組み合わせを数が小さい方からすべて試して行く
- ブロック対角化は、2ブロックに分けられる行・列の分け方を探索し
- 分割された対角ブロックについて、同じ処理を行い、ブロック対角化できなくなるまで再帰的に繰り返す
- 小さい行列なら動く
- 出力は、再帰処理にあたって、「n行・n列」が取り出されるたびに "number of selected cols/rows"として取り出された数が表示され、
- 最後に並べ替えた行、並べ替えた列の順列が示される
library(gtools)
DivideInto2BlockDiagonals<-function(M){
n<-length(M[,1])
found<-FALSE
seq1<-1:n
seq2<-1:n
ns<-c(n)
if(n==1)return(list(seq1=seq1,seq2=seq2,ns=ns))
for(i in 1:(n-1)){
if(found)break
chosen<-combinations(n,i)
for(j in 1:length(chosen[,1])){
if(found)break
X1<-chosen[j,]
X2<-(1:n)[-chosen[j,]]
for(k in 1:length(chosen[,1])){
if(found)break
Y1<-chosen[k,]
Y2<-(1:n)[-chosen[k,]]
M1<-M[X1,Y2]
M2<-M[X2,Y1]
if(sum(abs(M1))==0 && sum(abs(M2))==0){
seq1<-c(X1,X2)
seq2<-c(Y1,Y2)
print("number of selected cols/rows")
print(i)
ns<-c(i,n-i)
N1<-as.matrix(M[X1,Y1])
N2<-as.matrix(M[X2,Y2])
tmpseqs1<-DivideInto2BlockDiagonals(N1)
tmpseqs2<-DivideInto2BlockDiagonals(N2)
ns<-c(tmpseqs1$ns,tmpseqs2$ns)
tmpX1<-X1[tmpseqs1$seq1]
tmpX2<-X2[tmpseqs2$seq1]
tmpY1<-Y1[tmpseqs1$seq2]
tmpY2<-Y2[tmpseqs2$seq2]
seq1<-c(tmpX1,tmpX2)
seq2<-c(tmpY1,tmpY2)
found<-TRUE
}
}
}
}
list(seq1=seq1,seq2=seq2,ns=ns)
}
N<-10
M<-matrix(sample(c(0,1),N^2,prob=c(0.9,0.1),replace=TRUE),N,N)
diag(M)<-1
M<-M[sample(1:N),sample(1:N)]
outM<-DivideInto2BlockDiagonals(M)
outM
M[outM$seq1,outM$seq2]
nD<-sample(2:4,1)
nD<-3
aa<-c()
for(i in 1:nD){
aa<-c(aa,sample(1:4,1))
}
L<-list()
for(i in 1:nD){
L[[i]]<-matrix(1,aa[i],aa[i])
}
library(Matrix)
bdiag(L)
M4<- as.matrix(bdiag(L))
N<-length(M4[,1])
M5<-M4[sample(1:N),sample(1:N)]
M<-sign(M5)
M
outM<-DivideInto2BlockDiagonals(M)
outM
M[outM$seq1,outM$seq2]