Let's 'C#で数値計算' with Mol & Dsl

行列・ベクトルの高速演算

行列には「密・疎」、「一般・対称・三角・帯」、さらに double・Complex 等が組み合わさって 様々な種類があります。個々のデータ構造等の違いは行列・ベクトルのデータ構造と配列演算を参照してください。

行列のデータ構造が異なっていても、 Mol では宣言時の宣言文が異なるだけで、使用方法に違いはありません。 例えば、行列とベクトルの(A が行列、e、x、それに b はベクトル）演算 e = A*x - b は行列 A が何であれ、このように 書くことができるし、配列要素 A[i,j] にアクセスするには v = A[i,j] と書くだけです。

Mol ではこのようなデータ構造の違いをプログラマーから隠蔽、全ての行列やベクトルの演算（四則演算 + - * / ）に対応しています が以下の、_Blas や _VML 等の MKL に用意された専用演算メソッドに比べるとかなり低速になります。

行列サイズが非常に大きい場合は C = A*B-C 等と単純に記述せず、以下のメソッドの使用を検討してください。
C = A*B-C等のように四則演算子を用いればプログラムが分かりやすくなりますが、途中結果を保持する作業用行列等が 見えないところで作成される結果、メモリー不足に陥る場合等もあり得ます。

ベクトル・行列のノルム

ベクトルや行列の計算には結果等を確かめるためにノルムを計算することがよくあります。
以下、A は任意のベクトル、または行列。
呼出し形式名称ベクトル行列備考
v = A.Norm1()１－ノルム‖A‖1=Σi｜A[i]｜||A||1 = Maxj(Σi(|A[i,j]|)
v = A.Norm2()２－ノルム‖A‖2=(Σi(｜A[i]｜）2)1/2　--　（定義： A の最大特異値）スペクトルノルム、ベクトルはＥ－ノルムと同じ
v = A.NormE()Ｅ－ノルム‖A‖E=(Σi(｜A[i]｜）2)1/2‖A‖E = (ΣiΣj(|A[i,j]|)2)1/2ユークリッドノルム／フロベニウスノルム、ベクトルは２－ノルムと同じ
v = A.NormI()∞－ノルム‖A‖∞=Maxi(｜a[i]｜)‖A‖∞ = Maxi(Σj(|A[i,j]|))無限大ノルム

_Blas

_Blas クラスには行列同士、行列とベクトルの演算が（static メソッドとして）いくつか用意されています。

以下、

A は以下で表記する行列タイプです xxxxx の部分は Double か Complex で置き換えてください。
表記行列タイプ宣言
G密な一般行列MatrixDenseGeneralxxxxx
U密な上三角行列MatrixDenseUpperTrianglexxxxx
L密な下三角行列MatrixDenseUpperTrianglexxxxx
B密な帯行列MatrixDenseBandxxxxx
S密な対称行列MatrixDenseSymmetricxxxxx
H密なエルミート行列MatrixDenseHermite
SG疎な一般行列MatrixSparseGeneralxxxxx
SU疎な上三角行列MatrixSparseUpperTrianglexxxxx
SL疎な下三角行列MatrixSparseLowerTrianglexxxxx
SS疎な対称行列MatrixSparseSymmetricxxxxx
SH疎なエルミート行列MatrixSparseHermite

計算式：　C = α・opA(A)・opB(B) + β・C

呼出形式： _Blas.LeMm(C,A,B,opA,opB,α,β)
行列　　A:　G

計算式：　C = α・op(A)・B + β・C

呼出形式： _Blas.LeMm(C,A,B,op,α,β)
行列　　A:　SG,SS,SU,SL,SH

計算式：　y = α・op(A)・x + β・y

呼出形式： _Blas.LeMv(y,A,x,op,α,β)
行列　　A:　G,B,S(doubleのみ),H,SG,SS,SH,SU,SL

計算式：　A = A + α・x・y* + α*・y・x*

呼出形式： _Blas.LeMvv(A,α,x,y)
行列　　A:　S,H
y* はベクトル y の複素共役転置（実数対称行列の場合は転置）を意味し、結果は行ベクトルになります。従って、x・y* は行列になります。 α* は α の複素共役を意味します。

計算式：　A = A + α・x・y'

呼出形式： _Blas.LeMvv(A,α,x,y[,fConjY])
行列　　A:　G
y' はベクトル y の転置を意味し、結果は行ベクトルになります。従って、x・y' は行列になります。 A,x,y が複素数の時、fConjY = true とすれば y の複素共役転置を指定できます。

計算式：　A = A + α・x・x*

呼出形式： _Blas.LeMvv(A,α,x)
行列　　A:　S,H
x* はベクトル x の複素共役転置（実数なら単に転置）を意味し、結果は行ベクトルになります。従って、x* は行列になります。

計算式：　y = op(A)・x

呼出形式： _Blas.Mv(y,A,x,op)
行列　　A:　G

計算式：　y = op(A)・y

呼出形式： _Blas.Mv(y,A,op)
行列　　A:　U,L

計算式：　y = A・x

呼出形式： _Blas.Mv(y,A,x)
行列　　A:　SS

_VML

_VML はベクトルのマルチスレッド処理を強く意識したクラスです。

例えば y、x1、x2は全てサイズ(n)が等しい密なベクトル（double）として、 y = _VML.Add(y,x1,x2); は

 if (y == null || y.Count < n) y = new VectorDenseDouble(n);
 for(int i=1;i<=n;++i) y[i] = x1[i] + x2[i];
 return y;
と同等の計算を実行します。違いは y[i] = x1[i] + x2[i] を同時並行的に実行可能ということです。
例えば n=4 で使用できるスレッド数が 4 以上なら、これら要素ごとの 足し算は全て異なるスレッドに振り分けられて一度に計算されることになります。

以下、_VML クラスで用意されている代表的な(static)メソッドの一覧です。

また（下記テーブルにはない）y を省略した x1 = _VML.Add(x1,x2); という形式もあります。 これは

 for(int i=1;i<=n;++i) x1[i] = x1[i] + x2[i];
と同じで、第一引数の内容が上書きされます。

以下、各メソッド呼出しは特に備考に記述が無い限り double と Complex のベクトルが使用可能です。
呼出形式説明備考
y = _VML.Add(y,x1,x2) y[i] = x1[i] + x2[i]
y = _VML.Sub(y,x1,x2) y[i] = x1[i] - x2[i]
y = _VML.Mul(y,x1,x2) y[i] = x1[i] * x2[i]
y = _VML.Div(y,x1,x2) y[i] = x1[i] / x2[i]
y = _VML.Sqr(y,x) y[i] = x[i] * x[i]
y = _VML.MulByConj(y,x1,x2) y[i] = x1[i] * Complex.Conjugate(x2[i]) Complex のみ
y = _VML.Conjugate(y,x) y[i] = Complex.Conjugate(x[i]) Complex のみ
y = _VML.Abs(y,x) y[i] = |x[i]| double のみ
y = _VML.Inv(y,x) y[i] = 1.0/x[i] double のみ
y = _VML.LinearFrac(y,x1,a,b,x2,c,d) y[i] = (a*x1[i]+b)/(c*x2[i]+d) double のみ、a,b,c,dはスカラー
y = _VML.Sqrt(y,x) y[i] = (x[i])1/2
y = _VML.InvSqrt(y,x) y[i] = (x[i])-1/2 double のみ
y = _VML.Cbrt(y,x) y[i] = (x[i])1/3 double のみ
y = _VML.InvCbrt(y,x) y[i] = (x[i])-1/3 double のみ
y = _VML.Pow2o3(y,x) y[i] = (x[i])2/3 double のみ
y = _VML.Pow3o2(y,x) y[i] = (x[i])3/2 double のみ
y = _VML.Pow(y,x1,x2) y[i] = (x1[i])x2[i]
y = _VML.Powx(y,x,b) y[i] = (x1[i])b b はスカラー
y = _VML.Hypot(y,x1,x2) y[i] = (x1[i]2+x2[i]2)1/2 double のみ
y = _VML.Exp(y,x) y[i] = ex[i]
y = _VML.Log(y,x) y[i] = Loge(x[i])
y = _VML.Log10(y,x) y[i] = Log10(x[i])
y = _VML.Cos(y,x) y[i] = Cos(x[i])
y = _VML.Sin(y,x) y[i] = Sin(x[i])
y = _VML.Tan(y,x) y[i] = Tan(x[i])
y = _VML.CiS(y,x) y[i] = Cos(x[i])+i*Sin(X[i])iは虚数単位。y はComplex、x は double
y = _VML.Acos(y,x) y[i] = Cos-1(x[i])
y = _VML.Asin(y,x) y[i] = Sin-1(x[i])
y = _VML.Atan(y,x) y[i] = Tan-1(x[i])
y = _VML.Atan2(y,x1,x2) y[i] = Tan-1(x1[i]/x2[i]) double のみ
y = _VML.Cosh(y,x) y[i] = (ex[i]+e-x[i])/2双曲線余弦関数
y = _VML.Sinh(y,x) y[i] = (ex[i]+e-x[i])/2双曲線正弦関数
y = _VML.Tanh(y,x) y[i] = Sinh(x[i])/Cosh(x[i])双曲線正接関数
y = _VML.Acosh(y,x) y[i] = Acosh-1(x[i])
y = _VML.Asinh(y,x) y[i] = Asinh-1(x[i])
y = _VML.Atanh(y,x) y[i] = Atanh-1(x[i])
y = _VML.Erf(y,x) y[i] = (2/π1/2))∫0x[i]e-t*tdt誤差関数、double のみ
y = _VML.ErfInv(y,x) y[i] = Erf-1(x[i])逆誤差関数、double のみ
y = _VML.Erfc(y,x) y[i] = 1-Erf(x[i]) 相補誤差関数、double のみ
y = _VML.ErfcInv(y,x) y[i] = Erfc-1(x[i]) 逆相補誤差関数、double のみ
y = _VML.CdfNorm(y,x) y[i] = 1/(2π)1/2∫-∞x[i]e-t*t/2dt累積正規分布関数、double のみ
y = _VML.CdfNormInv(y,x) y[i] = CdfNorm-1(x[i]) 逆累積正規分布関数、double のみ
y = _VML.Gamma(y,x) y[i] = ガンマ関数値(x[i]) double のみ
y = _VML.Floor(y,x) y[i] = 負の無限方向に切り捨て(x[i]) 複素数は実数部と虚数部に適用されます。
y = _VML.Ceil(y,x) y[i] = 正の無限方向に切り上げ(x[i]) 複素数は実数部と虚数部に適用されます。
y = _VML.Trunc(y,x) y[i] = ゼロ方向の切り上げ（切り捨て）(x[i]) 複素数は実数部と虚数部に適用されます。
y = _VML.Round(y,x) y[i] = 最も近い整数に切り上げ（切り捨て）(x[i]) 複素数は実数部と虚数部に適用されます。
y = _VML.Frac(y,x) y[i] = 整数部分の削除(x[i])、（y[i]=x[i]-⎿x[i]⏌ x[i]≧0の場合、y[i]=x[i]-⎾x[i]⏋ x[i]<0の場合） 複素数は実数部と虚数部に適用されます。
y = _VML.Convolut(y,x1,x2) x1 と x2 の線形畳み込み計算
y = _VML.Correlat(y,x1,x2) x1 と x2 の線形相関計算