PuLP で変数の和や内積を計算する際の注意点

## TL; DR
PuLP で大きなモデル作るなら、numpy や pandas の sum や dot の使用は避ける。最低でも `pulp.lpSum` と `pulp.lpDot` を使い、場合によっては `LpAffineExpression` を自前で定義する。

## はじめに

数理最適化、特に MILP のモデリングツールとして知られている PuLP だが、Python 標準の sum や numpy.sum を使うと、モデルの構築が非常に遅くなるケースがある。今回、次の計算の速度を測定した。

- ベクトルの要素の総和
- ベクトル同士の内積

実験環境は、Google Colaboratory の無償版。Jupyter notebook 上のセルで、10 回走らせたうちの最良の計算時間を採用した（`%%timeit -r 10 -n 1`）。Python のバージョンは 3.7 で、各ライブラリのバージョンは下記の通り。

- numpy==1.19.5
- pandas==1.1.5
- PuLP==2.5.0


## ベクトルの要素の総和

PuLP の変数が入ったベクトル $\mathbf{y}$ があるときに、その要素の総和 $\sum_i y_i$ は、以下の方法で計算できる。

```python
z1 = sum(y)
z2 = np.sum(y)  # pandas.Series.sum と同じ
z3 = pl.lpSum(y)
z4 = pl.LpAffineExpression((yi, 1) for yi in y)  # y が独立変数の場合
z5 = pl.LpAffineExpression((k, v) for yi in y for var, coef in yi.items())  # y が従属変数の場合
```

z4 と z5 は少しトリッキーだが、PuLP の総和計算を高速化する方法として、以下の stack overflow の質問で紹介されている。

https://stackoverflow.com/questions/26643421/pulp-very-slow-when-adding-many-constraints

※独立変数と従属変数が混在している場合や、重複した変数が y に含まれる場合は、上記の方法では計算できない。その場合、要素の型を判定して場合分けする必要がある。

※この方法では、定数部分は無視されるので注意。微分すれば $0$ になるので最適化結果には影響しないが、目的関数の値に意味がある場合は修正の必要がある。定数も扱いたい場合は、`LpAffineExpression.constant` を取り出して、その総和を新しい `LpAffineExpression` に入れればよい。


## 実験結果

次のようにデータを作成し、前節の方法で計算した。

```python
import pandas as pd
import pulp as pl
import numpy as np

# Constants
a = np.random.rand(100)

# Variables
x = pd.Series(pl.LpVariable.dicts("x", range(1000)), index=range(1000))
X = pd.DataFrame(pl.LpVariable.dicts("X", (range(100), range(1000))), index=range(1000), columns=range(100))

# 速度比較を行う対象
y1 = x
y2 = 2 * x
y3 = X @ a
```

結果は次の通り。

| 対象 | x | 2 * x | X @ a |
|-|-|-|-|-|-|
| sum  | 239 [ms] | 271 [ms] | 92 [s] |
| numpy.sum  | 238 [ms] | 273 [ms] | 90 [s] |
| pulp.lpSum  | 982 [μs] | 1.84 [ms] | 150 [ms] |
| LpAffineExpression 再定義 | 405 [μs] | 909 [μs] | 128 [ms] |

標準の sum 関数や numpy.sum は、下 2 つに比べて数百倍の時間がかかってしまうことが分かる。また、`pulp.lpSum` よりも、`LpAffineExpression` を再定義する方が数十パーセント早い。


## ベクトル同士の内積

PuLP の変数が入ったベクトル $\mathbf{y}$ と、定数ベクトル $\mathbf{b}$ があるときに、その内積 $\mathbf{y}^\top \mathbf{b}$ は、以下の方法で計算できる。

```python
z1 = y @ b
z2 = pl.lpDot(y, b)
z3 = pl.LpAffineExpression(zip(y, b))  # y が全て独立変数の場合
z4 = pl.LpAffineExpression((var, bi * coef / 2) for bi, yi in zip(b, y) for var, coef in yi.items())  # y が全て従属変数の場合
```

※総和と同様に、独立変数と従属変数が混在している場合や、重複した変数が x に含まれる場合は、上記の方法では計算できない。

※z4 で 2 で割っているのは暫定的な処理。これをやらないとなぜか値が二倍されてしまう。PuLP 公式に [issue](https://github.com/coin-or/pulp/issues/483) は出しているが、特に反応はないため、これが仕様なのかバグなのかは分からない。


## 実験結果

次のようにデータを作成し、`y` と `b` の内積を前節の方法で計算した。

```python
import pandas as pd
import pulp as pl
import numpy as np

# Constants
a = np.random.rand(100)
b = np.random.rand(1000)

# Variables
x = pd.Series(pl.LpVariable.dicts("x", range(1000)), index=range(1000))
X = pd.DataFrame(pl.LpVariable.dicts("X", (range(100), range(1000))), index=range(1000), columns=range(100))

# 速度比較を行う対象
y1 = x
y2 = 2 * x
y3 = X @ a
```

結果は次の通り。

| 対象 | x | 2 * x | X @ a |
|-|-|-|-|
| np.dot | 241 [ms] | 234 [ms] | 58.7 [s] |
| lpDot | 10.4 [ms] | 19.7 [ms] | 1.18 [s] |
| LpAffineExpression 再定義 | 443 [μs] | 1.75 [ms] | 172 [ms] |

やはり numpy の内積だとかなり遅い。そして、`LpAffineExpression` を自前で定義するのが最速。


## 結論
和や内積の計算は、PuLP の関数を使うか、できれば `LpAffineExpression` を再定義しよう！

`lpSum` より `LpAffineExpression` を再定義する方が早い理由は分からないが、`lpSum` の方は汎用性を上げるために、型の判定などの処理が入っているのが原因かもしれない。