Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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};
}
@sasekazu

This comment has been minimized.

Copy link
Owner Author

sasekazu commented Jan 10, 2015

例えばランダムな対称行列を対角化するには以下のようにする。

var maxiteration = 1000;
var size = 30;
var A=numeric.random([size, size]);
for(var i=1; i<size; ++i) {
    for(var j=0; j<i; ++j) {
        A[i][j]=A[j][i];
    }
}
var result=jacobiDiagonalization(size, A, maxIteration, threshold);
var D = numeric.clone(result.DiagonalMatrix);
var P = numeric.clone(result.OrthodonalMatrix);
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.