jacobi法による対角化 jacobi diagonalization
// jacobi法による実対称行列の対角化 | |
// numeric.jsが必要 | |
// 他の関数との依存性排除ver | |
// size: 行列サイズ | |
// intputA: 対角化される入力行列 ※対称行列でなければならない | |
// 返り値 {DiagonalMatrix:A, OrthodonalMatrix:P}; | |
// A(入力時) = P A(対角化後) P.transpose が成り立つ | |
function jacobiDiagonalization(size, inputA, maxIteration, threshold) { | |
var A=numeric.clone(inputA); | |
// Pを単位行列で初期化 | |
var P=new Array(size); | |
for(var i=0; i<size; ++i) { | |
P[i]=new Array(size); | |
for(var j=0; j<size; ++j) { | |
if(i==j) { | |
P[i][j]=1; | |
} else { | |
P[i][j]=0; | |
} | |
} | |
} | |
var R=new Array(size); | |
for(var i=0; i<size; ++i) { | |
R[i]=new Array(size); | |
} | |
var count=0; | |
while(1) { | |
// search max |A[p][q]| | |
var p=0; | |
var q=1; | |
var val=Math.abs(A[0][1]); | |
for(var i=0; i<size-1; ++i) { | |
for(var j=i+1; j<size; ++j) { | |
if(val<Math.abs(A[i][j])) { | |
p=i; | |
q=j; | |
val=Math.abs(A[i][j]); | |
} | |
} | |
} | |
if(val<threshold) { | |
break; | |
} | |
// make rotation matrix | |
for(var i=0; i<size; ++i) { | |
for(var j=0; j<size; ++j) { | |
if(i==j) { | |
R[i][j]=1; | |
} else { | |
R[i][j]=0; | |
} | |
} | |
} | |
var th=0.5*Math.atan2(-2*A[p][q], A[p][p]-A[q][q]); | |
R[p][p]=Math.cos(th); | |
R[p][q]=Math.sin(th); | |
R[q][p]=-Math.sin(th); | |
R[q][q]=Math.cos(th); | |
// A <- R.transpose * A * R | |
A=numeric.dot(numeric.transpose(R), A); | |
A=numeric.dot(A, R); | |
// P <- P * R | |
P=numeric.dot(P, R); | |
++count; | |
} | |
return {DiagonalMatrix:A, OrthodonalMatrix:P}; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
This comment has been minimized.
例えばランダムな対称行列を対角化するには以下のようにする。