逆変換、特に行列の平方根

C# には、手続きを引数にする「delegate」の機能があります。この英単語は「使節」という意味らしいですが、実体は「candidate（候補）」に近いのではないでしょうか。稀にしか使わない機能だと思いますが、今回、逆関数を求めるために使ってみました。

例えば、次のプログラムでは、標準の関数 Math.Tan を引数とし、その逆関数（Math.Atan に相当）を実行しています。同じように、自作のものなど任意の関数（引数リストと返す変数の形が同じもの）を引数にできます。

x と y の両方がベクトルである場合の y = f(x) という変換で、y = y₀ を満たす x を求める場合、 y をノルムにしてしまって、その極値（普通は極小値ゼロ）を、最急降下法や共役勾配法で求めることがよく行われるように思います。ですが、ノルムしか知られない場合にそういう計算をするのは仕方ないものの、折角目的のベクトルが知られているなら、その豊富な情報を使わないのは勿体ないです。ベクトル値のまま、目的の点近傍につき逆行列を作って解きましょう。

次のプログラムの中の「vtovfuncinv」は、ベクトルからベクトルへの変換に対する逆変換を与えるもので、その第一引数が引き渡される元の変換を行う手続きです。第二引数には、初期条件とする（なるべく）近傍の点、第三引数には目的とする y₀ を置きます。

using System;
// 
public delegate void vfromvfunc(double[] y, double[] x);
// 
public class Test2 {       
    public void test2() {  
        double[][] M = new double[][] { 
            new double[]{ 8, 3, -2, 1 }, 
            new double[]{ 0, 7, -3, 2 }, 
            new double[]{ 1, 1, 6, -1 }, 
            new double[]{ -1, -2, -2, 4 } };
        int isz = M.Length; 
        double[][] N = Util.newmatelement(isz);
        Util.matsqrt(N, M); // === N が平方根行列として得られる 
        double[][] P = Util.newmatclear(isz, isz); 
        Util.matmulti(P, isz, isz, isz, N, N); // === 確かに平方根かチェック
        Util.matround(P, 12); // P == {{8,3,-2,1},{0,7,-3,2},{1,1,6,-1},{-1,-2,-2,4}}
    } // ===================================================================== 
} // ######################################################################### 
//
public class Util {       
    public static void matsqrt(double[][] N, double[][] M) {
        int isz = M.Length;
        matelement(N, isz); 
        double[] x = new double[isz*isz], y = new double[isz*isz];
        matsqrt_tovect(x, N); matsqrt_tovect(y, M);
        vfromvfuncinv(matsqrt_dblmat, x, y);
            // === x から y を得る関数（matsqrt_dblmat）を引数に取り、
            // y から x への計算を実行する。 
        matsqrt_tomat(N, x); 
    } // ===================================================================== 
    static void matsqrt_tovect(double[] vect, double[][] mat) {
        int isz = mat.Length;
        for (int i=0, ct=0; i<isz; i++) { 
            for (int j=0; j<isz; j++, ct++) {
                vect[ct] = mat[i][j]; } }
    } // ===================================================================== 
    static void matsqrt_dblmat(double[] y, double[] x) { // y <--- x calculation 
        int isz = (int)Math.Sqrt(x.Length);
        double[][] mat = newmatclear(isz, isz); 
        matsqrt_tomat(mat, x); 
        matmulti(mat, isz, isz, isz, mat, mat);
        matsqrt_tovect(y, mat);  
    } // ===================================================================== 
    static void matsqrt_tomat(double[][] mat, double[] vect) {
        int isz = mat.Length;
        for (int i=0, ct=0; i<isz; i++) { 
            for (int j=0; j<isz; j++, ct++) {
                mat[i][j] = vect[ct]; } }         
    } // =====================================================================     
    public static void vfromvfuncinv(vfromvfunc onemat, double[] x, double[] y) {
        int isz = x.Length;
        double[] xorig = new double[isz], yex = new double[isz + 1], 
            xtobeex = new double[isz + 1]; 
        vectcopy(xorig, isz, x); 
        vectcopy(yex, isz, y); yex[isz]=1; 
        vectcopy(xtobeex, isz, x); xtobeex[isz]=1;
        double[][] fctsA = newmatclear(isz + 1, isz + 1); 
        double[][] matB = newmatclear(isz + 1, isz + 1); 
        double[][] ivmatB = newmatclear(isz + 1, isz + 1); 
        double[][] ivQmok = newmatclear(isz + 1, isz + 1);
        double dconv = double.MaxValue; 
        //
        for (int ct=0; ct<4000; ct++) {
            for (int i=0; i<isz; i++) { 
                x[i] += 0.01; fctsA[i][i] = 0.01; fctsA[i][isz] = 1; 
                onemat(matB[i], x); matB[i][isz] = 1; 
                x[i] = xorig[i]; }
            onemat(matB[isz], x);  fctsA[isz][isz] = 1; matB[isz][isz] = 1; 
            matinv(ivmatB, isz + 1, matB);
            matmulti(ivQmok, isz + 1, isz + 1, isz+1, ivmatB, fctsA); 
            vectmultimat(xtobeex, isz+1, isz+1, yex, ivQmok);
            vectaddReform(x, isz, xtobeex); 
            vectcopy(xorig, isz, x);
            double d = getlen2(isz, xtobeex); 
            if (d >= dconv) { 
                break; }
            dconv = d; } 
    } // =====================================================================  
    public static void vectcopy(double[] x, int isz, double[] a) {
        for (int i = 0; i < isz; i++) { x[i] = a[i]; }
    } // ===================================================================== 
    public static double[][] newmatclear(int isz, int jsz) {
        double[][] mattemp = new double[isz][];
        for (int i = 0; i < isz; i++) { mattemp[i] = new double[jsz]; }
        return mattemp;
    } // ===================================================================== 
    public static double[][] newmatelement(int isz) {
        double[][] mattemp = new double[isz][];
        for (int i = 0; i < isz; i++) {
            mattemp[i] = new double[isz];
            vectclear(mattemp[i], isz); mattemp[i][i] = 1.0;
        }
        return mattemp;
    } // ===================================================================== 
    public static void vectclear(double[] x, int isz) {
        for (int i = 0; i < isz; i++) { x[i] = 0.0; }
    } // ===================================================================== 
    static public void matelement(double[][] matx, int vdim) {
        for (int i = 0; i < vdim; i++) {
            vectclear(matx[i], vdim); matx[i][i] = 1.0;
        }
    } // ===================================================================== 
    public static void mataddmultiReform(double[][] x,
            int isz, int msz, int jsz, double[][] a, double[][] b) {
        for (int i = 0; i < isz; i++) {
            for (int m = 0; m < msz; m++) {
                vectaddscmReform(x[i], jsz, b[m], a[i][m]); } }
    } // =====================================================================  
    public static void vectaddscmReform(
            double[] x, int isz, double[] a, double f) {
        for (int i = 0; i < isz; i++) { x[i] += a[i] * f; }
    } // =====================================================================  
    public static void matcopy(double[][] matx,
            int isz, int jsz, double[][] mat1) {
        for (int i = 0; i < isz; i++) {
            vectcopy(matx[i], jsz, mat1[i]); }
    } // ===================================================================== 
    public static void matmulti(double[][] x,
                int isz, int msz, int jsz, double[][] a, double[][] b) {
        double[][] mattemp = newmatclear(isz, jsz);
        mataddmultiReform(mattemp, isz, msz, jsz, a, b);
        matcopy(x, isz, jsz, mattemp);
    } // ===================================================================== 
    public static void vectaddmultimatReform( 
            double[] x, int isz, int jsz, double[] v, double[][] m) {
        for (int j = 0; j < jsz; j++) {
            vectaddscmReform(x, isz, m[j], v[j]); }
    } // ===================================================================== 
    public static void vectmultimat( 
            double[] x, int isz, int jsz, double[] v, double[][] m) { 
        // please take special care that m size is (jsz, isz) !!! 
        double[] vecttemp = new double[isz];
        vectaddmultimatReform(vecttemp, isz, jsz, v, m);
        vectcopy(x, isz, vecttemp);
    } // ===================================================================== 
    public static void vectaddReform(
            double[] x, int isz, double[] a) {
        for (int i = 0; i < isz; i++) { x[i] += a[i]; }
    } // ===================================================================== 
    public static void vectaddrangeReform(
            double[] x, int i0, int isz, double[] a) {
        for (int i = i0; i < isz; i++) { x[i] += a[i]; }
    } // ===================================================================== 
    public static void vectsubtrangeReform(
            double[] x, int i0, int isz, double[] a) {
        for (int i = i0; i < isz; i++) { x[i] -= a[i]; }
    } // ===================================================================== 
    public static double getlen2(int isz, double[] a) {
        double sum1 = 0.0;
        for (int i = 0; i < isz; i++) { sum1 += a[i] * a[i]; }
        return sum1;
    } // ===================================================================== 
    public static bool matinv(double[][] x, int isz, double[][] a) { 
	    int[] Ci=new int[isz], jaex=new int[isz]; 
		double[][] inva=newmatclear(isz, isz); 
		matcopy(inva,   isz,isz,   a); 
		for (int ii=0; ii<isz; ii++) { 
			int iconvent=ii;  double dconvent=Math.Abs(inva[ii][ii]); 
			for (int i=ii+1; i<isz; i++) { 
				double dnw=Math.Abs(inva[i][ii]); 
				if (dnw > dconvent) {iconvent=i;   dconvent=dnw; } }  
			if (iconvent == ii) { 
				Ci[ii]=0; 
				if (dconvent == 0) {
					return false; } 
			} else if (inva[ii][ii]*inva[iconvent][ii] >= 0) { 
				vectaddrangeReform(
						inva[ii],   ii, isz,   inva[iconvent]); 
				Ci[ii]=iconvent; 
			} else {
				vectsubtrangeReform(
						inva[ii],   ii, isz,   inva[iconvent]); 
				Ci[ii]=-iconvent; }    
			double invaovii = 1/inva[ii][ii]; 
			for (int j = ii + 1; j<isz; j++) {inva[ii][j] *= invaovii; } 
			inva[ii][ii]=invaovii; // +++ 
			// ### know non zero ### 
			int jasz=0; 
			for (int j=ii+1; j<isz; j++) { 
				if (inva[ii][j] != 0) {jaex[jasz++]=j; } } 
			///////////////////////////////  
			for (int i=ii+1; i<isz; i++) {  double invaivii=inva[i][ii];  
				if (invaivii != 0) { 
					for (int j1=0; j1<jasz; j1++) { 
						int j=jaex[j1]; inva[i][j] -= inva[ii][j]*invaivii; 
					} } } }
		double[] col=new double[isz]; 
		for (int ii=0; ii<isz; ii++) { 
			vectclear(col,  isz);  
			col[ii]=1.0; 
		    for (int i=0; i<isz; i++) {  
			    if (Ci[i] > 0) {    // row addition again 
				    col[i] += col[Ci[i]]; 
			    } else if (Ci[i] < 0) {  
				    col[i] -= col[-Ci[i]]; } 
			    if (col[i] == 0) {continue; } 
			    col[i] *= inva[i][i];  // row multiplication again 
			    for (int k=i+1; k<isz; k++) {   // row subtraction again 
				    double invaivii=inva[k][i]; 
				    if (invaivii != 0) {col[k] -= col[i]*invaivii; } } } 
		    for (int i=isz-1; i >= 0; i--) { 
			    double sum1=0.0; 
			    double[] invaii=inva[i]; 
			    for (int j=i+1; j<isz; j++) {sum1 += invaii[j]*col[j]; } 
			    col[i] -= sum1; } 
			for (int j=0; j<isz; j++) {x[j][ii]=col[j]; } } 
        return true; 
    } // ===================================================================== 
	public static void matround(double[][] mat, int n) {
		for (int i=0; i<mat.Length; i++) {
			for (int j=0; j<mat[0].Length; j++) {
				double d = mat[i][j], rd = Math.Round(d, n);
				mat[i][j] = rd; } } 
	} // ===================================================================== 
} // #########################################################################

上のプログラムの中で、「matsqrt」は行列の平方根を与えます。ここに平方根とは、A、B を（正方）行列で B = A・A とするとき、B から A を求めるものです。ここでは行列の平方というものを、行列の各要素（を列挙したベクトル）から、平方たる行列の各要素（を列挙したベクトル）への変換と捉え、その逆を求める形で実装しています。これらの要素の関係は、行列の次元を N とすると、N×N 個の変数による、N×N 個連立の二次式ということになります。この計算例は四次元ですが、他の次元でも勿論できます。