jacobi法による対角化 jacobi diagonalization

gistfile1.js

// 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};
}

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

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);