解析エンジニアの自動化 blog

コツコツと自動化した方法を残す blog

C# で 逆行列を求める



こんにちは。
仕事の自動化にやりがいと達成感を感じるガッくんです。


この記事の目次



背景・目的


C# ってベクトル計算、行列計算のライブラリが豊富ではありません。
私が知らないだけかもしれませんが。

実に簡単なアルゴリズムですが逆行列を解くプログラムを書いたので、まとめておきます。



動作環境


Windows 7
Visual Studio 2017



プログラム

ソースコード


using System;
using System.Collections.Generic;
 
namespace inverse_matrix
{
    class Program
    {
        static void Main(string[] args)
        {
            List> matrix = new List>();
 
            matrix.Add(new List());
            matrix[matrix.Count - 1].Add(1);
            matrix[matrix.Count - 1].Add(1);
            matrix[matrix.Count - 1].Add(-1);
 
            matrix.Add(new List());
            matrix[matrix.Count - 1].Add(-2);
            matrix[matrix.Count - 1].Add(-1);
            matrix[matrix.Count - 1].Add(1);
 
            matrix.Add(new List());
            matrix[matrix.Count - 1].Add(-1);
            matrix[matrix.Count - 1].Add(-2);
            matrix[matrix.Count - 1].Add(1);
 
            List> inv = matrixinv(matrix);
 
            for (int i = 0; i < inv.Count; i++)
            {
                for (int j = 0; j < inv[i].Count; j++)
                {
                    if (j == 0) Console.Write(string.Format("{0,10}", inv[i][j]));
                    else Console.Write(" , " + string.Format("{0,10}", inv[i][j]));
                }
                Console.WriteLine("");
            }
 
            Console.ReadKey();
        }
 
        private static List> matrixinv(List> matrix)
        {
            // 戻り値 List
            List> matrixinv = new List>();
 
            // matrix の行数取得
            int n = matrix.Count;
 
            // matrix の列数が行数と同じであることの確認
            int m = -1;
            for (int i = 0; i < n; i++)
            {
                m = matrix[i].Count;
                if(m != n)
                {
                    return matrixinv;
                }
            }
 
            int max;
 
            // matrixinv に単位行列を格納する
            for (int i = 0; i < n; i++)
            {
                // matrixinv に行を追加する
                matrixinv.Add(new List());
 
                // 1 行ずつ値の代入
                for (int j = 0; j < n; j++)
                {
                    if (i == j) matrixinv[matrixinv.Count - 1].Add(1.0); // 1 を代入
                    else matrixinv[matrixinv.Count - 1].Add(0);          // 0 を代入
                }
            }
 
            // 逆行列の計算
            for(int k = 0; k < n; k++)
            {
                max = k;
                double tmp;
 
                for(int j = k + 1; j < n; j++)
                {
                    if (Math.Abs(matrix[j][k]) > Math.Abs(matrix[max][k]))
                    {
                        max = j;
                    }
                }
 
                if (max != k)
                {
                    for (int i = 0; i < n; i++)
                    {
                        // 入力行列側
                        tmp = matrix[max][i];
                        matrix[max][i] = matrix[k][i];
                        matrix[k][i] = tmp;
                       
                        // 単位行列側
                        tmp = matrixinv[max][i];
                        matrixinv[max][i] = matrixinv[k][i];
                        matrixinv[k][i] = tmp;
                    }
                }
 
                tmp = matrix[k][k];
 
                for (int i = 0; i < n; i++)
                {
                    matrix[k][i] /= tmp;
                    matrixinv[k][i] /= tmp;
                }
 
                for (int j = 0; j < n; j++)
                {
                    if (j != k)
                    {
                        tmp = matrix[j][k] / matrix[k][k];
                        for (int i = 0; i < n; i++)
                        {
                            matrix[j][i] = matrix[j][i] - matrix[k][i] * tmp;
                            matrixinv[j][i] = matrixinv[j][i] - matrixinv[k][i] * tmp;
                        }
                    }
                }
            }
 
            //逆行列が計算できなかった時
            for (int j = 0; j < n; j++)
            {
                for (int i = 0; i < n; i++)
                {
                    if (double.IsNaN(matrixinv[j][i]))
                    {
                        matrixinv[j][i] = 0; // NaN を ゼロ に置き換える
                    }
                }
            }
 
            // 戻り値
            return matrixinv;
        }
    }
} 



結果


図1 結果



コメント

他にも固有値解析や特異値分解などもまとめておきたいですね。



以上