MathNet.Numerics.MKLをLinux上のF#で使う

# 更新履歴

- 2019-12-13: @htsignさんにコードブロックの色付けをしていただきました。QiitaのMarkdownよくわかっていないので助かりました。ありがとうございました。

# はじめに

私はF#を使って数値計算のプログラムを書いています。数値計算が専門な訳ではなく、ロボットアームやその制御システムのモデリングやシミュレーションに数値計算を使っています。この用途なら線形代数が使えるならばどの言語でも達成することができますが、エコシステムの強さやコードの書きやすさからMATLABやPythonを使ってる人が多いと思います。

MATLABやPythonがあるにも関わらず、なぜ私がF#を使っているのかというと、

- 実機を動かすことがあるので、.NET APIのシリアル通信を使えるF#は魅力的。
- 同じ.NET基盤のC#に比べて書き方がスクリプト的なので、手軽に書くことができる。
- スクリプト的に書けるにも関わらず、それなりに速度が出る。
- 関数型プログラミング言語に興味があった。

などの理由からです。さらにF#を使う追い風として、

- .NetCoreの登場で、WindowsとLinuxで同じコードを動かすことが簡単になった。
- .NetCore 3.0からLinuxでも.NET APIのシリアル通信が部分的に使えるようになった。

といった要素があります。

F#で線形代数を扱うライブラリのデファクトスタンダードはMathNet.Numericsだと思います。このMathNet.NumericsとF#を使うと、次のNumpyとPythonによる線形代数方程式`A x = y`の`x`を求めるスクリプト、

```python:solve_ale.py
import numpy.random as random
import numpy.linalg as linalg

N = 5000

A = random.normal(0, 1, (N, N))
y = random.normal(0, 1, (N,))
x = linalg.solve(A, y)
print(x)
```

と同様のスクリプトをほぼ同等の記述量で次のように書くことができます。

```fsharp:solve_ale.fsx
#r @"/path/to/MathNet.Numerics.dll"

open MathNet.Numerics.LinearAlgebra

[<Literal>]
let N = 5000

let A = Matrix<float>.Build.Random(N, N)
let y = Vector<float>.Build.Random(N)
let x = A.Solve(y)

printfn "%A" x
```

このようにF#とMathNet.Numericsの組み合わせで線形代数を使った数値計算を行うことができますが、この方法はIntel CPU上での実行速度に問題があります（ほかのCPUについては調べていません）。上のF#のコードはPythonのコードに比べて、100倍以上の実行時間がかかってしまいます。これほど大きな速度差が生じるのは、Numpyの`linalg.solve`はIntel MKLで計算されるのに対して、MathNet.Numericsの`A.Solve`は.NET基盤上で実装されたコードで計算されるからです。Intel MKLはIntel CPUに対して最適化された数値計算ライブラリで非常に高速であることが知られています。NumpyとMathNet.Numericsの速度はIntel MKLを使っているかどうかで大きく差がついていると考えられます。

この速度の問題はMathNet.Numerics.MKLを使うことで解決することができます。MathNet.Numerics.MKLはMathNet.Numericsの計算部分をIntel MKLを使ったものに置き換えます。これによりNumpyとMathNet.Numericsの速度差が縮まり、F#とMathNet.Numericsでも線形代数を高速に計算することができます。WindowsではFFTなど線形代数以外も高速化されるのですが、Linuxでは後述するライブラリのバージョン問題で線形代数のみの高速化となることに注意してください。

この記事ではLinux上のF#でMathNet.Numerics.MKLを線形代数に対して使う手順と、その際に生じる問題点の暫定的な解決策を紹介します。上の`solve_ale.fsx`と同等のF#コンソールアプリケーションをMathNet.NumericsとMathNet.Numerics.MKLで動作させることを目的とします。Ubuntu 18.04 64bitと.NET Core SDK 3.1を対象に説明しますので、異なるディストリビューションやSDKでは適宜読み替えてください。

# プロジェクトの作成と計算コードの実装

まずMathNet.Numerics.MKLを使わずに線形代数方程式を解くプログラムを作ります。Fsiスクリプトである`solve_ale.fsx`との違いは、ビルドされるコンソールアプリケーションとして実装される点です。

`dotnet`コマンドでF#コンソールアプリケーションを作成します。`dotnet`便利ですよね。

```console
$ dotnet new console -lang F# -o solveALE
```

プロジェクトディレクトリに移動して、MathNet.Numericsパッケージを追加してからソースコードを編集します。

```console
$ cd solveALE
$ pwd
/path/to/solveALE
$ dotnet add package MathNet.Numerics
$ <好きなエディタ> Program.fs
```

編集したコードを次に示します。

```fsharp:Program.fs
open System

open MathNet.Numerics.LinearAlgebra

[<Literal>]
let N = 5000

[<EntryPoint>]
let main argv =
    let A = Matrix<float>.Build.Random(N, N)
    let y = Vector<float>.Build.Random(N)
    let x = A.Solve(y)

    printfn "%A" x
    0 // return an integer exit code
```

`solve_ale.fsx`と比べると、dllを読み込む行がなくなり、計算部分がエントリーポイントである`main`関数に含まれるようになりました。

# MathNet.Numerics.MKLパッケージの追加

Linuxに対するMathNet.Numerics.MKLのパッケージ名は`MathNet.Numerics.MKL.Linux`ですので、これを`dotnet`で追加します。

```console
$ pwd
/path/to/solveALE
$ dotnet add package MathNet.Numerics.MKL.Linux
```

MathNet.Numerics.MKL.Linuxがある状態で`dotnet build`すると、MathNet.Numerics.MKLに関連するライブラリが`/path/to/solveALE/bin/Debug/netcoreapp3.1/x64`以下に配置されます。Windows版（MathNet.Numerics.MKL.Win）ではこの状態でMKLが自動的に利用されますが、Linux版ではうまくいきません。

Linux版でMKLが利用されない原因は次の2つです。

- ライブラリの探索に失敗して、MathNet.Numerics.MKLが読み込まれない。
- MathNet.Numerics.MKL.Linuxのバージョンが古くて、FFTの最適化に対応していない。

それぞれの原因を次の方法で解決します。

- 原因: ライブラリの探索に失敗して、MathNet.Numerics.MKLが読み込まれない。
    - 解決策: 探索パスに共有ライブラリの場所を追加する。
- 原因: MathNet.Numerics.MKL.Linuxのバージョンが古くて、FFTの最適化に対応していない。
    - 解決策: FFT対応を諦めて線形代数のみMKLを適用する。

## 探索パスへのライブラリの追加

MathNet.NumericsはMathNet.Numerics.MKLのDLLである`/path/to/solveALE/bin/Debug/netcoreapp3.1/x64/MathNet.Numerics.MKL.dll`を自動的に読み込もうとしますが、そのDLLが依存する共有ライブラリ`/path/to/solveALE/bin/Debug/netcoreapp3.1/x64/libiomp5.so`の探索がうまくいっていないようです。このことは`ldd`コマンド、

```console
$ ldd /path/to/solveALE/bin/Debug/netcoreapp3.1/x64/MathNet.Numerics.MKL.dll
```

の`libiomp5.so`についての結果が`not found`となっていることからわかります。したがって`libiomp5.so`をlinuxから見えるようにしなければなりません。

暫定的に共有ライブラリを探索パスへ追加するには`LD_LIBRARY_PATH`環境変数を編集します。`/path/to/solveALE/bin/Debug/netcoreapp3.1/x64`を環境変数に追加します。

```console
$ export LD_LIBRARY_PATH=/path/to/solveALE/bin/Debug/netcoreapp3.1/x64:${LD_LIBRARY_PATH}
```

ライブラリが有効になったことを`ldd`コマンドで確認します。

```console
$ ldd /path/to/solveALE/bin/Debug/netcoreapp3.1/x64/MathNet.Numerics.MKL.dll
```

コマンド結果を確認して、`not found`になっている箇所がなければ設定完了です。

## 線形代数のみへのMKLの適用

ライブラリを読み込める状態にしても、まだMKLが利用されない場合があります（されることもある。差が現れる理由は不明）。このときMathNet.Numericsは処理をMKLに切り替えようと試行してはいるのですが、失敗して標準の計算実装にフォールバックしています。MathNet.Numericsにおいて、MKLを適用できる機能は線形代数とFFTがあります。両方の機能に対して同時に切り替えを試みて、FFT機能について失敗した結果、うまくいっている線形代数機能までも標準実装にフォールバックしてしまいます。

FFTに対してMKLによる計算を有効化できない理由は、`MathNet.Numerics.MKL`のバージョンが古いからです。2.2.0よりも新しければ問題ありませんが、nugetで入手できる`MathNet.Numerics.MKL.Linux`はバージョン2.0.0で更新が止まっています。それに対して`MathNet.Numerics.MKL.Win`は2.3.0まで公開されているので、FFTの問題が発生しなかったわけです。根本的な解決は`MathNet.Numerics.MKL.Linux`の更新を待つことですが、すぐに更新される保証はありません。

そこでFFTについては諦めて、線形代数の計算のみでMKLを使うように手動で切り替えます。計算前に`MathNet.Numerics.Providers.LinearAlgebra.LinearAlgebraControl.UseNativeMKL()`を呼ぶようにコードを変更します。

```fsharp:Program.fs
open System

open MathNet.Numerics.LinearAlgebra

[<Literal>]
let N = 5000

[<EntryPoint>]
let main argv =
    MathNet.Numerics.Providers.LinearAlgebra.LinearAlgebraControl.UseNativeMKL()
    let A = Matrix<float>.Build.Random(N, N)
    let y = Vector<float>.Build.Random(N)
    let x = A.Solve(y)

    printfn "%A" x
    0 // return an integer exit code
```

この関数は線形代数の処理に対してMKLを使用するように設定します。

# MathNet.Numerics.MKL導入前後の速度比較

`solveALE`プロジェクトの実行速度を`MathNet.Numerics.MKL`の導入前後で比較しました。

導入前:

    $ time /path/to/solveALE/bin/Debug/netcoreapp3.1/solveALE
    real    8m8.321s
    user    8m7.923s
    sys     0m0.429s

導入後:

    $ time /path/to/solveALE/bin/Debug/netcoreapp3.1/solveALE
    real    0m4.350s
    user    0m6.742s
    sys     0m0.272s

Macbook air mid 2011のi7モデル上のUbuntuで実行しました。導入後の速度は導入前の約120倍となり、Numpyに近い速度が出るようになったので満足です。