ライブラリを使わずにC#で行列計算（転置，積，逆行列，擬似逆行列）

行列 $A \in \mathbb{R}^{n \times n}$と行列 $B \in \mathbb{R}^{n \times m}$をC#で記述すると次のようになる．

C#

int n = 2; //ここでは取り敢えず2を格納
int m = 3; //ここでは取り敢えず3を格納
double[,] A = new double[n, n]; //左が行，右が列
double[,] B = new double[n, m]; //左が行，右が列

---

行列 $A$ の転置 $A^{\rm T}$ をC#で記述するには，以下の関数をコピペして Transpose(A)．

C#

double[,] Transpose(double[,] A){

            double[,] AT = new double[A.GetLength(1), A.GetLength(0)];

            for (int i = 0; i < A.GetLength(1); i++)
            {
                for (int j = 0; j < A.GetLength(0); j++)
                {
                    AT[i, j] = A[j, i];
                }
            }

            return AT;

        }

Transpose(A)の使用例

C#

double[,] AT = new double[n, n];
AT = Transpose(A);

---

行列同士の積 $AB$ をC#で記述するには，以下の関数をコピペして MatrixTimesMatrix(A, B)．

C#

double[,] MatrixTimesMatrix(double[,] A, double[,] B){

            double[,] product = new double[A.GetLength(0), B.GetLength(1)];

            for (int i = 0; i < A.GetLength(0); i++)
            {
                for (int j = 0; j < B.GetLength(1); j++)
                {
                    for (int k = 0; k < A.GetLength(1); k++)
                    {
                        product[i, j] += A[i, k] * B[k, j];
                    }
                }
            }

            return product;

        }

MatrixTimesMatrix(A, B)の使用例( $C = AB$ )

C#

double[,] C = new double[n, m];
C = MatrixTimesMatrix(A,B);

---

逆行列 $A^{-1}$ をC#で記述するには，以下の関数をコピペしてinverseMatrix(A)．

C#

double[,] inverseMatrix(double[,] A){

            int n = A.GetLength(0);
            int m = A.GetLength(1);

            double[,] invA = new double[n, m];

            if (n == m)
            {

                int max;
                double tmp;

                for (int j = 0; j < n; j++)
                {
                    for (int i = 0; i < n; i++)
                    {
                        invA[j, i] = (i == j) ? 1 : 0;
                    }
                }

                for (int k = 0; k < n; k++)
                {
                    max = k;
                    for (int j = k + 1; j < n; j++)
                    {
                        if (Math.Abs(A[j, k]) > Math.Abs(A[max, k]))
                        {
                            max = j;
                        }
                    }

                    if (max != k)
                    {
                        for (int i = 0; i < n; i++)
                        {
                            // 入力行列側
                            tmp = A[max, i];
                            A[max, i] = A[k, i];
                            A[k, i] = tmp;
                            // 単位行列側
                            tmp = invA[max, i];
                            invA[max, i] = invA[k, i];
                            invA[k, i] = tmp;
                        }
                    }

                    tmp = A[k, k];

                    for (int i = 0; i < n; i++)
                    {
                        A[k, i] /= tmp;
                        invA[k, i] /= tmp;
                    }

                    for (int j = 0; j < n; j++)
                    {
                        if (j != k)
                        {
                            tmp = A[j, k] / A[k, k];
                            for (int i = 0; i < n; i++)
                            {
                                A[j, i] = A[j, i] - A[k, i] * tmp;
                                invA[j, i] = invA[j, i] - invA[k, i] * tmp;
                            }
                        }
                    }

                }

                
                //逆行列が計算できなかった時の措置
                for (int j = 0; j < n; j++)
                {
                    for (int i = 0; i < n; i++)
                    {
                        if (double.IsNaN(invA[j, i]))
                        {
                            Console.WriteLine("Error : Unable to compute inverse matrix");
                            invA[j, i] = 0;//ここでは，とりあえずゼロに置き換えることにする
                        }
                    }
                }
                

                return invA;

            }
            else
            {
                Console.WriteLine("Error : It is not a square matrix");
                return invA;
            }

        }

inverseMatrix(A)の使用例

double[,] invA = new double[n, n];
invA = inverseMatrix(A);

---

前述した関数たちを使って，行列 $B^{\rm T}AB$ を計算してみよう．

C#

double[,] BTAB = new double[m, m];
BTAB = MatrixTimesMatrix(Transpose(B),MatrixTimesMatrix(A,B));

前述した関数たちを使って，行列 $B^{\rm T}(BB^{\rm T})^{-1}$ を計算してみよう．

C#

double[,] pseudoB = new double[m, n];
double[,] BT = new double[m, n];
BT = Transpose(B);
pseudoB = MatrixTimesMatrix(BT,inverseMatrix(MatrixTimesMatrix(B,BT)));

なお，$n<m$ かつ行列 $B$ が行フルランクであれば，行列 $B^{\rm T}(BB^{\rm T})^{-1}$ は $B$ の擬似逆行列である．