Excelで数値計算～ALGLIBの導入～

スポンサーリンク

Excelで数値計算をしようと思うと、自分で１からそれぞれのプログラムを組む必要があります。

手軽に計算できる方法を調べたところ、ALGLIBというライブラリではExcelでも呼び出しができる数値計算・データ処理ライブラリも提供していました。

導入方法を紹介します。

ALGLIBとは

ALGLIB は、数値計算や最適化アルゴリズムを提供するライブラリで、多言語対応かつ移植性の高い構成が特徴です。特に VBA を含む複数のプログラミング言語で使用可能であり、LAPACK と同様の機能を一部提供しています。

日本語のWikiがありますのでもう少し詳しく知りたい方は、参照してみてください。

ALGLIB の主な特徴

  1. 多言語対応
    • C++、C#、Python、VBA など、さまざまな言語向けに設計されています。VBA 用のインターフェースも用意されており、Excel VBA プロジェクトに組み込みやすいです。
  2. 機能の幅広さ
    ALGLIB は LAPACK に匹敵する数値計算機能を多く備えています。一部の機能例:
    • 線形代数 (LU 分解、QR 分解、特異値分解など)
    • 最適化 (線形最適化、非線形最適化)
    • 方程式ソルバー (線形方程式、非線形方程式)
    • 固有値計算
    • 数値積分
    • 補間 (スプライン、ポリノミアル補間)
  3. DLL の利用で環境依存性が低い
    ALGLIB の DLL をプロジェクトに追加するだけで利用可能。インストールが不要で、VBA から簡単に呼び出せます。
  4. オープンソース & 商用ライセンス
    • 無料版 (GPL ライセンス)
    • 商用版 (独自ライセンス) が選べます。商用利用を予定していない場合、無料版で十分です。

無料版でも機能は変わらないようですので、大企業での利用や大規模の計算に使わなければ無料版を選択してください。

ALGLIBのダウンロード

こちらのページのからダウンロードできます。

ALGLIB Free Edition
www.alglib.net

残念ながら、VBA用は旧バージョンの2.6で止まっているようですが、使う分には問題無さそうです。

alglib-2.6.0.vb6.zip

解凍した\alglib-2.6.0.vb6\vb6\srcにbasファイルを提供してくれています。

必要なbasファイルをExcelにインポートして使います。

主なパッケージ

パッケージ名ユニット名日本語の説明対応する .bas ファイル
DataAnalysisdforest決定木（ランダムフォレスト）分類器dforest.bas
kmeansK-means++ クラスタリングkmeans.bas
lda線形判別分析lda.bas
linreg線形モデルlinreg.bas
logitロジットモデルlogit.bas
mlpbase基本的なニューラルネットワーク操作mlpbase.bas
mlpeニューラルネットワークのアンサンブルモデルmlpe.bas
mlptrainニューラルネットワークの学習mlptrain.bas
pca主成分分析pca.bas
DiffEquationsodesolver常微分方程式ソルバーodesolver.bas
FastTransformsconv高速な実数/複素数の畳み込み演算conv.bas
corr高速な実数/複素数の相互相関演算corr.bas
fft実数/複素数のFFT（高速フーリエ変換）fft.bas
fht実数の高速ハートレー変換fht.bas
Integrationautogk自動適応型1次元積分autogk.bas
gkqガウス-クロード積分公式の生成gkq.bas
gqガウス型積分公式の生成gq.bas
Interpolationidwint逆距離加重法による補間/フィッティングidwint.bas
lsfit線形および非線形の最小二乗法lsfit.bas
polint多項式補間/フィッティングpolint.bas
psplineパラメトリックスプライン補間pspline.bas
ratint有理補間/フィッティングratint.bas
spline1d1次元スプライン補間spline1d.bas
spline2d2次元スプライン補間spline2d.bas
LinAlgablasLevel 2およびLevel 3のBLAS演算ablas.bas
bdsvd双対化SVDbdsvd.bas
evd固有値ソルバーevd.bas
inverseupdate行列逆行列のSherman-Morrison更新inverseupdate.bas
ldltLDLT分解ldlt.bas
matdet行列式の計算matdet.bas
matgenランダム行列の生成matgen.bas
matinv行列の逆行列matinv.bas
ortfac実数/複素数のQR, LQ, 二重三重対角化, Hessenberg分解ortfac.bas
rcond条件数の推定rcond.bas
schurシュール分解schur.bas
sdet対称行列の行列式sdet.bas
sinverse対称行列の逆行列sinverse.bas
spdgevd一般化対称固有値問題のソルバーspdgevd.bas
srcond対称行列の条件数推定srcond.bas
svd特異値分解svd.bas
trfacLUおよびコレスキー分解trfac.bas

例題

実数の非対称行列\(A\)の固有値を求めてみます。

$$
A=
\begin{bmatrix}
1& 2& 3\\
0& 4& 5\\
1& 0& 6
\end{bmatrix}
$$

ファイルのインポート

srcフォルダ内の必要なファイルをVBエディター上でインポートします。

  • Excel を開き、Alt + F11 を押して VBA エディタを起動します。
  • 「ファイル」→「ファイルのインポート」 を選択。
  • 下表のファイルをインポートします。

それぞれの役割は以下の通りです。

ファイル名名称説明関連する関数
ap.basALGLIB Primitives?基本的なデータ型やユーティリティ関数を含むComplex 型定義, C_Add, C_Sub など
blas.basBasic Linear Algebra Subprograms?線形代数の基本操作を含む複素数対応の BLAS 関数
evd.basEigenvalue Decomposition?固有値分解や複素数を扱う高レベルな関数を含むHMatrixEVD など
hsschur.basHessenberg-Schur?Hessenberg化してSchur 分解InternalSchurDecompositionなど
reflections.basReflections?行列計算における反射計算GenerateReflectionなど
rotations.basRotaions?行列の回転ApplyRotationsFromTheLeftなど

手間ですが、インポートしたモジュールは名前をModule*から変更しておいた方が良いでしょう。

evd.basモジュールのRMatrixEVD関数を使って計算しますが、この関数から様々な関数を呼び出すため、これら全てのモジュールが必要になります。

RMatrixEVD関数

実行列に対する固有値を求めるには、こちらの関数を使うことになります。

入力パラメータ

  1. A:
    • 行列。
    • サイズは N×N
    • 配列のインデックスは [0～N-1, 0～N-1] で指定される。
  2. N:
    • 行列 A のサイズ（行列は正方行列）。
  3. VNeeded:
    • 固有ベクトルの計算に関するフラグ。
    • 値によって計算する内容が変わる:
      • 0: 固有値のみを計算。
      • 1: 右固有ベクトルを計算。
      • 2: 左固有ベクトルを計算。
      • 3: 左右両方の固有ベクトルを計算。

出力パラメータ

  1. WR:
    • 固有値の実数部分。
    • 配列 [0..N-1] に格納される。
  2. WI:
    • 固有値の虚数部分。
    • 配列 [0..N-1] に格納される。
  3. VL:
    • 左固有ベクトル（必要な場合）。
    • 配列 [0..N-1, 0..N-1] に格納される。
  4. VR:
    • 右固有ベクトル（必要な場合）。
    • 配列 [0..N-1, 0..N-1] に格納される。

固有値と固有ベクトルの形式

  1. 実数固有値の場合:
    • WI[i] = 0 の場合、固有値は実数。
    • 固有ベクトルは対応する行列の i 列に格納される。
  2. 複素数固有値の場合:
    • WI[i] > 0 の場合、固有値は次のようにペアで格納される:
      • \(WR[i] + \sqrt{-1} \cdot WI[i]\)
      • \(WR[i+1] + \sqrt{-1} \cdot WI[i+1]\)（これは前者の共役）
    • 対応する固有ベクトルは、以下のように格納される:
      • 行列の i 列に固有ベクトルの実部。
      • 行列の i+1 列に固有ベクトルの虚部。

サンプルプログラム

インポートしたモジュール上だと作業しにくいので、Module1（もしくはプロジェクトウインドウの標準モジュールを右クリックし、「挿入」→「標準モジュール」で別のモジュールを追加）にプログラムを記入していきます。

Sub ExampleRMatrixEVD()
    Dim A(0 To 2, 0 To 2) As Double
    Dim WR() As Double
    Dim WI() As Double
    Dim VR() As Double
    Dim VL() As Double '左固有ベクトルが不要な場合は省略可能
    Dim N As Long
    Dim VNeeded As Long
    Dim I As Long, J As Long

    ' 行列 A の定義
    A(0, 0) = 1: A(0, 1) = 2: A(0, 2) = 3
    A(1, 0) = 0: A(1, 1) = 4: A(1, 2) = 5
    A(2, 0) = 1: A(2, 1) = 0: A(2, 2) = 6

    ' パラメータの設定
    N = 3
    VNeeded = 3 ' 左右固有ベクトル両方を計算

    ' 固有値と右固有ベクトルを計算
    Call RMatrixEVD(A, N, VNeeded, WR, WI, VL, VR)

    ' 結果の出力
    Debug.Print "固有値 (実数部, 虚数部):"
    For I = 0 To N - 1
        Debug.Print WR(I) & ", " & WI(I)
    Next I

    Debug.Print "右固有ベクトル:"
    For I = 0 To N - 1
        For J = 0 To N - 1
            Debug.Print "VR(" & I & "," & J & ") = " & VR(I, J)
        Next J
    Next I

    Debug.Print "左固有ベクトル:"
    For I = 0 To N - 1
        For J = 0 To N - 1
            Debug.Print "VL(" & I & "," & J & ") = " & VL(I, J)
        Next J
    Next I
End Sub

計算結果

イミディエイトウインドウに以下が出力されます。

固有値 (実数部, 虚数部):
7.04096459487746, 0
1.08849668592497, 0
2.87053871919756, 0
右固有ベクトル:
VR(0,0) = 0.633107295508666
VR(0,1) = 1
VR(0,2) = -0.706921069287353
VR(1,0) = 1
VR(1,1) = 0.34965381743329
VR(1,2) = -1
VR(2,0) = 0.608192918975493
VR(2,1) = -0.203603649647201
VR(2,2) = 0.225892256160487
左固有ベクトル:
VL(0,0) = 0.165536477543332
VL(0,1) = 1
VL(0,2) = 0.534605346436768
VL(1,0) = 0.108871032449493
VL(1,1) = -0.686930353241035
VL(1,2) = -0.946655464022549
VL(2,0) = 1
VL(2,1) = 8.84966859249752E-02
VL(2,2) = 1

固有値

$$
\lambda =
\begin{bmatrix}
7.04096459487746\\
1.08849668592497\\
2.87053871919756
\end{bmatrix}
$$

右固有ベクトル
$$
VR =
\begin{bmatrix}
0.633107295508666& 1& -0.706921069287353\\
1& 0.34965381743329& -1\\
0.608192918975493& -0.203603649647201& 0.225892256160487
\end{bmatrix}
$$

左固有ベクトル
$$
VL =
\begin{bmatrix}
0.165536477543332& 1& 0.534605346436768\\
0.108871032449493& -0.686930353241035& -0.946655464022549\\
1& 8.84966859249752E-02& 1
\end{bmatrix}
$$

固有ベクトルは最大値１で正規化されずに出力されているようです。

固有ベクトルは、共役ペアを考慮した形で圧縮されて出力されているので、どれが虚数かを判断するのには慣れが必要になりそうです。なお、今回の例では実数のみで構成された固有ベクトルです。

Numpyとの比較

計算結果を比較するために、PythonでNumpyを使って確認してみます。

Pythonのコード

確認に使ったコードは以下。

import numpy as np

# Define the example matrix
A = np.array([[1, 2, 3],
              [0, 4, 5],
              [1, 0, 6]])

# Compute eigenvalues and right eigenvectors for A
eigenvalues_A, right_eigenvectors_A = np.linalg.eig(A)

# Compute eigenvalues and right eigenvectors for A^T (transpose of A)
eigenvalues_AT, right_eigenvectors_AT = np.linalg.eig(A.T)

# Display results
print("Original Matrix A:")
print(A)
print("\nEigenvalues of A:")
print(eigenvalues_A)
print("\nRight Eigenvectors of A (columns):")
print(right_eigenvectors_A)

print("\nTranspose of Matrix A (A^T):")
print(A.T)
print("\nEigenvalues of A^T:")
print(eigenvalues_AT)
print("\nRight Eigenvectors of A^T (columns):")
print(right_eigenvectors_AT)

計算結果

Original Matrix A:
[[1 2 3]
 [0 4 5]
 [1 0 6]]

Eigenvalues of A:
[7.04096459 1.08849669 2.87053872]

Right Eigenvectors of A (columns):
[[ 0.47577536  0.92699459 -0.56767265]
 [ 0.75149246  0.3241272  -0.80302126]
 [ 0.45705239 -0.18873948  0.18139628]]

Transpose of Matrix A (A^T):
[[1 0 1]
 [2 4 0]
 [3 5 6]]

Eigenvalues of A^T:
[7.04096459 1.08849669 2.87053872]

Right Eigenvectors of A^T (columns):
[[ 0.16238002  0.82207688  0.36191781]
 [ 0.10679508 -0.56470956 -0.64086803]
 [ 0.98093197  0.07275108  0.67698129]]

Right Eigenvectors of A^Tが左固有ベクトルに相当します。（行列Aの転置に対する右固有ベクトル）

比較

NumpyとALGLIBでは、出力されている桁数が違うものの答えは同等でした。

固有ベクトルについては、Numpyが大きさ１（ノルム１）で正規化されているのに対し、ALGLIBでは最大値１（もしくは-1）で正規化されているため、見た目の数字としては違いがあります。これは正規化方法が違うだけで、計算結果は同等と言えます。

まとめ

導入が少しだけ面倒なものの、他のライブラリを使うよりは手軽にVBAで計算環境が構築できると思います。

試してみる価値はあると思います。

スポンサーリンク

コメント

タイトルとURLをコピーしました