Created
January 10, 2015 09:38
-
-
Save sasekazu/a406c51adca34a37f2cc to your computer and use it in GitHub Desktop.
jacobi法による対角化 jacobi diagonalization
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// 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
例えばランダムな対称行列を対角化するには以下のようにする。