さんすう教室w

.NET 用の数学ライブラリ Math.Net Numerics を使った行列計算から最小二乗近似まで

「プログラムをしていると行列の計算をしたいことがありませんか?」

え、ないって?w
まあまあ、あるとして聴いて下さいませw

「必要なら数値計算ライブラリ買うかフリーの探せばいいじゃん?」

と思うかもしれませんが C# であんまりいいのないんですよね・・・ もちろんフリーでも C/C++ ならば結構あるんですが、 C# 用でないといまいち使いにくかったりして結局、自作ライブラリを C# で作って使っていました。 あと、仕事の場合には有料のライブラリでも構わないのですが、 その場合にしても「ちょっと逆行列が計算したい」程度の用途に お金を払うのは勿体ないというのもあります。

というわけで、今までは自作ライブラリで用は足りていたのですが、 所詮自作は自作、精度面や速度面で良い選択といえないと思っていました。 なんせプログラムや数学の専門家でもないので、 参考書片手に古典的かつ教科書的なアルゴリズムを採用しているわけですから どうしても見劣りしてしまいます。

でも、見つけました! そんで、見つけた数値演算ライブラリ Math.Net Numerics が 思いのほか使いやすいライブラリでしたのでご紹介したいと思います。 このライブラリ、個人で使用する分には十分な機能がそろっていて機能的に申し分ありません。 これを機会に自作の行列演算クラスを Math.Net Numerics へ置き換えようかなと思っています。

この小ネタでは、よく使う以下の機能、主に線形代数関連についての簡単な説明と、 プログラム例を私の忘備録として書いておきます。

補足 : 数値を扱うプログラムで行列は非常に強力な道具となるので、 行列をあまり知らない方もこれを機会に勉強するのも良いと思います。
私が思うにドラクエで例えると、知ってるのと知らないでは、”銅の剣”と”破壊ハヤブサ”くらいに攻撃力の差があると思いますw (わかりにくい?w)

(1) ダウンロードと使い方

まあ、ここで躓く人はいないと思いますが、とりあえず・・・

http://mathnetnumerics.codeplex.com/

へ行って「DOWNLOADS」のタブをクリックして、「Binaries」というところをクリックしますとダウンロードできます。


2014/06/06 現在のバージョンです
補足 : NuGetといわれるもので簡単にインストールできるそうですが、そんなの知らないんで手動っすw

で、「MathNet.Numerics-2.6.1.30.zip」というファイルがダウンロードできたと思います。 その中の「Net40」フォルダが目的のファイルとなります。 ちなみに、名前から想像できるかと思いますが .NET のバージョンは 4.0 以降となります。

補足 : 現在進行中の ver3.0 では .NET3.5 版が実装されるようです。ベータ版には確かに .NET3.5 のフォルダがありましたが全ての機能が実装されていないようです。

あとは解凍してどっかに置いとけばいいです。私の場合は、テスト用に作成したプロジェクト 「TestMath.NetNum」のプロジェクト直下に展開しました。


フォルダの中身

で、あとは参照設定右クリック「参照の追加」で Net40 の中の「MathNet.Numerics.dll」を指定するだけでOKです!

補足 : 今回 Math.Net numerics を一押しするのはこの手軽さです。
    ・でっかいパッケージをダウンロードしてインストールする必要がない。
    ・path の設定などのめんどくさい設定がない。
    ・さらに、C/C++ 用ライブラリの C# 向けのラッパーでありがちな、 元になったライブラリの dll が必要だったりせずに、 純粋に .NET のクラスライブラリのみで動作する。
    ・もちろんフリーで、商用もOK!!(MITライセンスは GNU より緩いので非常に使い良いライセンスです。)
3Dゲームを作りたい時には DirectX、OpenGL、画像処理がメインなら OpenCV など、 同様なことができそうなフリーのライブラリは多々ありますが、 例えば「ちょっとこのデータの近似が必要なんだけど・・・」という様な用途には重すぎます。
逆行列の計算をするだけなのに DirectX のインストールを要求するのはダサすぎますよね?w

これらが私の評価ポイントです。

(2) まずは逆行列を求めてみよう

まあ、足し算、かけ算なんかは普通にできますので、スルーの方向で・・・ まずは逆行列から行ってみましょう! 逆行列さんには連立方程式を解く度にお世話になっておりますw
さて、例として以下の計算をやってみます。

\[ M_{1} = \left( \begin{array}{ccc} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 9 & 8 \end{array} \right) \] \[ M_{2} = M^{-1} \] \[ M_{2}M_{1} = \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array} \right) ← 単位行列 I になることを確認してみる! \]
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

using MathNet.Numerics.LinearAlgebra.Double;

namespace TestMath.NetNum
{
    class Program
    {
        static void Main(string[] args)
        {
            // こんな風に行列を初期化できる
            var M1 = DenseMatrix.OfArray(new double[,] { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 9, 8 } });
            // 逆行列
            var M2 = M1.Inverse();

            Console.WriteLine("元の行列 M1");
            Console.WriteLine(M1);
            Console.WriteLine("M1 の逆行列");
            Console.WriteLine(M2);
            Console.WriteLine("M2 * M1 は単位行列 I になっているはず!");
            Console.WriteLine(M2 * M1);
        }
    }
}
		

実行結果です

元の行列 M1
DenseMatrix 3x3-Double
           1            2            3
           4            5            6
           7            9            8
M1 の逆行列
DenseMatrix 3x3-Double
    -1.55556      1.22222    -0.333333
     1.11111     -1.44444     0.666667
    0.111111     0.555556    -0.333333
M2 * M1 は単位行列 I になっているはず!
DenseMatrix 3x3-Double
           1  2.22045E-15  2.22045E-15
-1.77636E-15            1 -8.88178E-16
           0 -4.44089E-16            1
続行するには何かキーを押してください . . .
		

まあ、このくらいは当然ですね。
次行ってみましょう!

(3) 疑似逆行列を求めてみる

手っ取り早く近似したい時にお世話になる疑似逆行列ですが、ずばりそのものは実装されていないようです(されてたらごめんw)。 ですが、特異値分解が実装されているので簡単に計算できます。

疑似逆行列というのは逆行列を一般化したものです。 皆さんもご存じの通り逆行列は n×n みたいな正方行列しか求められませんが、 疑似逆行列は n×m のような正方ではない行列についての逆行列を与えます。 ちなみに、記号は逆行列は \(M^{-1}\) でしたが、疑似逆行列は \(M^{+}\) と書きます。

\[ \left( \begin{array}{ccc} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 9 & 8 \end{array} \right)^{-1}~~こういうのは逆行列があるが、 \] \[ \left( \begin{array}{ccc} 1 & 2 \\ 4 & 5 \\ 7 & 8 \end{array} \right)^{-1} ~~こういうのはダメ!絶対ダメ!! ~~でも、~~ \left( \begin{array}{ccc} 1 & 2 \\ 4 & 5 \\ 7 & 8 \end{array} \right)^{+}~~これはOKなんだ! \]
さて、\(M\) の疑似逆行列 \(M^{+}\) は以下の様に計算できます。
\[ M^{+} = (M^{T}M)^{-1}M^{T} \] (\(M^{T}\) は \(M\) の転置行列です)

この疑似逆行列は、以下の様に計算すると 単位行列 \(I\) になるので逆行列のように機能することがわかります。 \[ M^{+}M = I \]

まあ、さっきの様に計算してもいいのですが、一般的には特異値分解というものを使用して計算します。
特異値分解とは行列 \(M\) を以下の様に分解することです。詳しくは教科書やwikiでも見てもらいたいです。

\[ M = U \Sigma V^{T} \]

この様に分解した行列で以下の様に計算すると \(M\) の疑似逆行列 \( M^{+} \) を計算できます。

\[ M^{+} = V \Sigma^{+} U^{T} \]

ちなみに、Math.Net Numerics で \( \Sigma^{+} \) を計算する際には、 \( \Sigma \) が S というベクトルで出力されますので、ベクトルの段階で逆数を計算して 対角行列に変換し転置させると良いと思います。

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

using MathNet.Numerics.LinearAlgebra.Double;

namespace TestMath.NetNum
{
    class Program
    {
        static void Main(string[] args)
        {
            // こんな風にも行列を初期化できる
            var M1 = new DenseMatrix(3, 2);

            M1[0, 0] = 1;
            M1[0, 1] = 2;
            
            M1[1, 0] = 4;
            M1[1, 1] = 5;

            M1[2, 0] = 7;
            M1[2, 1] = 8;

            // 特異値分解(true にしてね)
            var svd = M1.Svd(true);

            // ベクトル S の逆数から対角行列を作る(もっといい方法がありそう?)
            var S = new DiagonalMatrix(M1.RowCount, M1.ColumnCount, (1/svd.S()).ToArray());
            
            // 疑似逆行列を計算する
            var M2 = svd.VT().Transpose() * S.Transpose() * svd.U().Transpose();

            Console.WriteLine("元の行列 M1");
            Console.WriteLine(M1);
            Console.WriteLine("M1 の疑似逆行列");
            Console.WriteLine(M2);
            Console.WriteLine("M2 * M1 は単位行列になっているはず!");
            Console.WriteLine(M2 * M1);
        }
    }
}
		
// ベクトル S の逆数から対角行列を作る(もっといい方法がありそう?)
var S = new DiagonalMatrix(M1.RowCount, M1.ColumnCount, (1/svd.S()).ToArray());
		

DiagonalMatrix は対角行列クラスで、まず \(M_{1}\) と同じサイズで初期化しています。 次にベクトル S の逆数を対角成分に取りたいので、(1/svd.S()).ToArray()の部分で、 ベクトル S の逆数をとってさらに double の配列化しています。
(実は W も定義されているのですが、効率の良い逆数の取り方がわからなかった・・・w)

実行結果です

元の行列 M1
DenseMatrix 3x2-Double
           1            2
           4            5
           7            8
M1 の疑似逆行列
DenseMatrix 2x3-Double
    -1.16667    -0.333333          0.5
           1     0.333333    -0.333333
M2 * M1 は単位行列になっているはず!
DenseMatrix 2x2-Double
           1            0
-1.33227E-15            1
続行するには何かキーを押してください . . .
		

うーん、いい感じですね。
まあ、3行で書けるのでそう大変でもないですし、必要ならクラス化すれば良いかと思います。

class MyMath
{
    public static Matrix Pinv(Matrix M)
    {
        var svd = M.Svd(true);
        var S = new DiagonalMatrix(M.RowCount, M.ColumnCount, (1 / svd.S()).ToArray());
        return (Matrix)(svd.VT().Transpose() * S.Transpose() * svd.U().Transpose());
    }
}
		

こんな感じでよく使う機能をまとめたヘルパークラスを用意しておいても良いと思います。 疑似逆行列は便利なものなので、このページの最後の最小二乗近似での使用例も見て下さい。

補足 : ちなみに、特異値分解や次に紹介する固有値、あとLU分解などは、”using MathNet.Numerics.LinearAlgebra.Double;” で宣言しないと、例のように”M.Svd(true);”というように使用することができません。 まあ、素直に using すればよいのでしょうが、using すると都合が悪いようでしたら以下の様にすればよいです。

var svd = M.Svd(true);
			

の変わりに、

var svd = MathNet.Numerics.LinearAlgebra.Double.ExtensionMethods.Svd(M, true);
			

この様にすればよい。

では、次行ってみましょう!

(4) 固有値、固有ベクトルを求めてみる

「なんのこっちゃ?」な人は教科書やwikiを見てもらうとして、 「何に使うの?」ですが、私の最近の例では三次元点群の姿勢を求めるために使用しました。

例えば、

\[ M_{1} = \left( \begin{array}{ccc} 2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{array} \right) \]

のような対称行列の固有ベクトルは必ず実数の直交行列となります。この固有ベクトルを \(R\) とすると、

\[ R^{T} R = I \]

この様に単位行列 \(I\) となり、直交行列であることを確認できます。
また、固有値を \(A\) とすると・・・

\[ R^{T} M_{1} R = diag(A) \] \[ ちなみに~~diag(A) = \left( \begin{array}{ccc} A_{0} & 0 & 0 \\ 0 & A_{1} & 0 \\ 0 & 0 & A_{2} \end{array} \right)~~ということです。 \]

つまり、\(M_{1}\) を \(R\) で \(A\) を対角成分とする対角行列へ対角化できるはずなので、これを確認してみます。

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

using MathNet.Numerics.LinearAlgebra.Double;

namespace TestMath.NetNum
{
    class Program
    {
        static void Main(string[] args)
        {
            // とりあえず対称行列で初期化
            var M1 = DenseMatrix.OfArray(new double[,] { { 2, 1, 1 }, { 1, 2, 1 }, { 1, 1, 2 } });

            // 固有値を計算
            var evd = M1.Evd();

            // 固有ベクトル R
            var R = evd.EigenVectors();
            // 固有値 A
            var A = evd.EigenValues();

            Console.WriteLine("元の行列 M1");
            Console.WriteLine(M1);
            Console.WriteLine("固有値 A");
            Console.WriteLine(A);
            Console.WriteLine("固有ベクトル R");
            Console.WriteLine(R);
            Console.WriteLine("");
            Console.WriteLine("RT * R = I ← 直行行列って事だね");
            Console.WriteLine(R.Transpose() * R);
            Console.WriteLine("RT * M1 * R = A");
            Console.WriteLine(R.Transpose() * M1 * R);
        }
    }
}
		

実行結果です

元の行列 M1
DenseMatrix 3x3-Double
           2            1            1
           1            2            1
           1            1            2
固有値 A
DenseVector 3-Complex
      (1, 0)       (1, 0)       (4, 0)
固有ベクトル R
DenseMatrix 3x3-Double
   0.0244713     -0.81613      0.57735
   -0.719025     0.386872      0.57735
    0.694553     0.429258      0.57735

RT * R = I ← 直行行列って事だね
DenseMatrix 3x3-Double
           1  2.22045E-16  2.77556E-16
 2.22045E-16            1  2.77556E-17
 2.77556E-16  2.77556E-17            1
RT * M1 * R = A
DenseMatrix 3x3-Double
           1  5.55112E-17  -9.4369E-16
 2.22045E-16            1  2.77556E-17
-8.88178E-16            0            4
続行するには何かキーを押してください . . .
		

物体の姿勢を求めるには、慣性テンソルの固有ベクトルを求める必要がありますので、興味がある方はググってみるといいかも? です。

(5) おまけで、最小二乗近似ってやつをやってみる

最後に少し実用的な例をということで、疑似逆行列で最小二乗近似をやってみましょう。

例題は簡単に以下の1次関数へ、

\[ y=ax+b \]
以下の x, y それぞれ11組のデータをフィッティングさせてみます。
x = 0,1,2,3,4,5,6,7,8,9,10
y = -0.577, 2.512, 3.000, 5.660, 8.330, 10.256, 12.699, 14.371, 16.756, 17.136, 20.121

計算の詳細は教科書でも読んでもらうとして、求めたいパラメータ a, b は以下の様に計算すれば近似できます。

x, y を以下の様に行列をおきます。 \[ X = \left( \begin{array}{ccc} 0 & 1 \\ 1 & 1 \\ \vdots & \vdots \\ 9 & 1 \\ 10 & 1 \\ \end{array} \right) ~~,~~ Y = \left( \begin{array}{ccc} -0.577 \\ 2.512 \\ \vdots \\ 17.136 \\ 20.121 \\ \end{array} \right) \] これで以下の様な式を考えて、a, b を求めるわけです。 \[ X\left( \begin{array}{ccc} a \\ b \\ \end{array} \right)=Y\\ \] ちなみに、行列を計算すると x, y 11個分についての11組の1次関数になっていることを確認してください。 \[ ax_{0}+b=y_{0}\\ ax_{1}+b=y_{1}\\ \vdots\\ ax_{9}+b=y_{9}\\ ax_{10}+b=y_{10}\\ \] あとは両辺に \( X \) の疑似逆行列 \( X^{+} \) を掛ければ、 \[ \left( \begin{array}{ccc} a \\ b \\ \end{array} \right)=X^{+}Y\\ \] この様に、パラメータ a, b が求まるわけです。
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

using MathNet.Numerics.LinearAlgebra.Double;

namespace TestMath.NetNum
{
    class Program
    {
        static void Main(string[] args)
        {
            // X, Y を用意する
            var X = DenseMatrix.Create(11, 2, (r, c) => 1);
            X.SetColumn(0, new double[] { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 });

            var Y = DenseMatrix.OfColumnVectors(new DenseVector(new double[] { -0.577, 2.512, 3.000, 5.660, 8.330, 10.256, 12.699, 14.371, 16.756, 17.136, 20.121 }));

            // 疑似逆行列を計算する
            var svd = X.Svd(true);
            var S = new DiagonalMatrix(X.RowCount, X.ColumnCount, (1 / svd.S()).ToArray());
            var pinv = svd.VT().Transpose() * S.Transpose() * svd.U().Transpose();

            Console.WriteLine("X");
            Console.WriteLine(X);
            Console.WriteLine("Y");
            Console.WriteLine(Y);
            Console.WriteLine("pinv X の逆行列");
            Console.WriteLine(pinv);
            Console.WriteLine("pinv * Y で a, b を計算する");
            Console.WriteLine(pinv * Y);
        }
    }
}
		
var X = DenseMatrix.Create(11, 2, (r, c) => 1);
X.SetColumn(0, new double[] { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 });
		

この部分は、11行2列の行列を 1 で初期化(ラムダ式 (r, c) => 1 で全要素へ 1 を設定しています)した行列を いったん用意して、改めて1列目へ X を設定しています。

var XA = new double[] { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };

// X, Y を用意する
var X = DenseMatrix.Create(11, 2, (r, c) => (c == 0) ? XA[r] : 1);
		

ちなみにラムダ式でいろいろできます。事前に X の配列が準備されていれば、こういうのもありです。

var Y = DenseMatrix.OfColumnVectors(new DenseVector(new double[] { -0.577, 2.512, 3.000, 5.660, 8.330, 10.256, 12.699, 14.371, 16.756, 17.136, 20.121 }));
		

この部分は、double 配列でベクトルを生成、それから行列 Y を生成しています。

実行結果です

X
DenseMatrix 11x2-Double
           0            1
           1            1
           2            1
           3            1
           4            1
           5            1
           6            1
         ...          ...
          10            1
Y
DenseMatrix 11x1-Double
      -0.577
       2.512
           3
        5.66
        8.33
      10.256
      12.699
         ...
      20.121
pinv X の逆行列
DenseMatrix 2x11-Double
  -0.0454545   -0.0363636   -0.0272727   -
5
    0.318182     0.272727     0.227273
4
pinv * Y で a, b を計算する
DenseMatrix 2x1-Double
     2.04586
   -0.205318
続行するには何かキーを押してください . . .
		

というわけで、求めたい直線の式は、

\[ y = 2.04593x - 0.205172 \]

ですね。


近似結果
ちなみに 2次式 \( y = ax^{2}+bx+c \) で近似したかったら以下の様になります。 \[ X = \left( \begin{array}{ccc} 0 & 0 & 1 \\ 1 & 1 & 1 \\ \vdots & \vdots & \vdots \\ 81 & 9 & 1 \\ 100 & 10 & 1 \\ \end{array} \right) ~~,~~ Y = \left( \begin{array}{ccc} -0.577 \\ 2.512 \\ \vdots \\ 17.136 \\ 20.121 \\ \end{array} \right) \] あとは同様に考えて、a, b, c を求めます。 \[ \left( \begin{array}{ccc} a \\ b \\ c \\ \end{array} \right)=X^{+}Y\\ \]

どうでしょうか? 普通の多項式近似ならこれだけ覚えておけば困ることはないでしょう。

Math.Net Numerics にはこれ以外にも一杯機能があるみたいですが、私的にこの程度で満足レベルです。 問題点は日本語の説明が少ないのと、私自身の数学のレベルが低いので使いこなせないと言うことでしょうか?w でも、ちょっと行列の計算がしたい時には結構使えそうですよねこれ?

前回へ戻る

次回へすすむ