離散カントロビッチ問題のエントロピー正則化とその双対問題

離散カントロビッチ問題とは?

商品を  n 個の工場から,  m 個のお店に輸送することを考えます.
商品の移動に掛かる労力が「質量」と「始点と終点の距離」の積で決まっているとき, トータルの労力が最小になるような移し方はどんなものでしょうか?

数学的に定式化すると

それぞれの工場での生産量を  a \in \mathbb{R}_+^n として表し, それぞれのお店での需要量を  b \in \mathbb{R}_+^m で表します,
行列  C = (c_{ij})_{ij} として単位体積を工場  i からお店  j へ輸送したときの距離(労力)を表します.

さらに商品全体の総重量を  1 で正規化し,

\begin{align*}
\sum_{i}a_i = 1,\quad \sum_{j}b_j = 1.
\end{align*}

と仮定して一般性を失いません*1.

数理計画問題としての定式化

離散カントロビッチ問題は, 以下の最小化問題により定式化されます*2:

\begin{equation}\tag{P}\label{P}
\mathop{\rm minimize}_{P\in H} \langle C,P \rangle
\end{equation}

ここで, 実行可能領域  H

\begin{align*}
H &= H_1 \cap H_2 \cap H_3,\\
H_1 &= \left\{p \in \mathbb{R}^{n\times m} \, |\, p_{ij} \ge 0\right\},\\
H_2 &= \left\{p \in \mathbb{R}^{n\times m} \, |\, P 1_m = a \right\},\\
H_3 &= \left\{p \in \mathbb{R}^{n\times m} \, |\, P^T1_n = b\right\}.
\end{align*}

により記述されます.

  • 1つ目の条件は輸送量はマイナス値を取らないこと,
  • 2つ目の条件はそれぞれの輸送元  i について, 送り出した商品の総量が  a_i であること,
  • 3つ目の条件はそれぞれの輸送先  j について, 受けとった商品の総量が  b_j であること

をそれぞれ表します.

 H_2, H_3 を等式制約と見なすと,
対応する Lagrange 関数  L: H_1 \times \mathbb{R}^n\times \mathbb{R}^m\rightarrow \mathbb{R} を,

\begin{equation}
L(P,f,g)
= \langle{C},{P}\rangle - \langle{f},{P1_m -a}\rangle - \langle{g},{P^T1_n- b}\rangle
\end{equation}

と定めることにより, 問題 (\ref{P}) を

\begin{align}\tag{P}
\mathop{\rm minimize}_{P\in H_1} \max_{(f,g)\in \mathbb{R}^n\times \mathbb{R}^m}L(P,f,g)
\end{align}

と書き直すことができます.
以後は最小化問題と呼んだときには, こちらの表現を念頭に置いて考察しましょう.

したがってその双対問題は, min と max をひっくり返して計算することにより,

\begin{align}\tag{Q}\label{Q}
&\mathop{\rm maximize}_{(f,g)\in \mathbb{R}^n\times \mathbb{R}^m}
\min_{P\in H_1}L(P,f,g)\\
&\quad = \mathop{\rm maximize}_{(f,g)\in \mathbb{R}^n\times \mathbb{R}^m} \left\{
\langle{f},{a}\rangle +\langle{g},{b}\rangle -\min_{P\in H_1} \langle{C-f1_m^T-1_ng^T},{P}\rangle\right\}\\
&\quad = \mathop{\rm maximize}_{(f,g)\in R(C)} \left\{
\langle{f},{a}\rangle +\langle{g},{b}\rangle\right\}, \quad R(C) = \left\{(f,g)\in \mathbb{R}^n\times \mathbb{R}^m| f_i + g_j \le C_{ij} \right\}
\end{align}

となる制約条件付き最大化問題として表すことができることが, 数理最適化の一般論からわかります.

エントロピー正則化

エントロピー関数を

\begin{equation}
H(P) = -\sum_{i,j} p_{ij}\left(\log{(p_{ij})}-1\right)
\end{equation}

により定義します*3.

さらにパラメータ  \varepsilon >0 を固定します.

最小化問題 \ref{P} におけるコスト関数にエントロピー関数  H(P) を加えて,

\begin{equation}\tag{P$_\varepsilon$}\label{P_eps}
\mathop{\rm minimize}_{P\in H'} \left\{\langle{C},{P}\rangle -\varepsilon H(P)\right\}
\end{equation}

により最小化問題を緩和することをエントロピー正則化と呼びます*4.

このとき,  H' =H_2 \cap H_3 で, 不等号制約  P \in H_1 は自動的に外せます.
なぜならば  P\not\in H_1 のとき,  p_{ij}<0 より  -\varepsilon H(P) = +\infty であり, 最小化問題の解の候補になりえないからです.

対応する Lagrange 関数  \mathcal{E}: \mathbb{R}^{n\times m} \times \mathbb{R}^n\times \mathbb{R}^m\rightarrow \mathbb{R} は,

\begin{equation}
\mathcal{E}(P,f,g)
= \langle{C},{P}\rangle - \varepsilon H(P)- \langle{f},{P1_m -a}\rangle - \langle{g},{P^T1_n- b}\rangle
\end{equation}

なので, 双対問題は

\begin{align}\tag{Q$_\varepsilon$}\label{Q_eps}
&\mathop{\rm maximize}_{(f,g)\in \mathbb{R}^n\times \mathbb{R}^m}
\min_{P\in \mathbb{R}^{n\times m}}\mathcal{E}(P,f,g)\\
&\quad = \mathop{\rm maximize}_{(f,g)\in \mathbb{R}^n\times\mathbb{R}^m} \left\{
\langle{f},{a}\rangle +\langle{g},{b}\rangle - \varepsilon\langle{\exp{(f/\varepsilon)}},{K\exp{(g/\varepsilon)}}\rangle\right\}
\end{align}

であり, こちらは制約条件のない最大化問題となっています.

ただしここで,  K_{ij} = \exp{(-C_{ij}/\varepsilon)} であり, 式変形には停留条件:

\begin{align*}
&\dfrac{\partial \mathcal{E}(P,f,g)}{\partial P_{ij}} = C_{ij} +\varepsilon \log(P_{ij}) -f_i -g_j = 0\\
&\Leftrightarrow P_{ij} = \exp{\left(\frac{f_i+g_j - C_{ij}}{\varepsilon}\right)}
\end{align*}

を用いました.
エントロピー関数の形状が凸であり, また境界近傍の傾きが  \log により発散していることから, 微分により停留条件を求めることができています.

以上をまとめると,

  1. 問題 (\ref{P}) ... affine なので, 最適値は頂点に存在し, 一意とは限らない
  2. 問題 (\ref{Q}) ... 線形計画法の双対問題なので双対ギャップは現れないが, 制約条件が残る
  3. 問題 (\ref{P_eps}) ... strictly convex なので最適値は一意に存在し, エントロピー関数のおかげで最適値は内部に存在する
  4. 問題 (\ref{Q_eps}) ... strictly convex の双対問題なので双対ギャップは現れず, 制約条件もない

ことがわかります.

さらに, 問題 (\ref{P}) は最大流問題に帰着させることができ, 問題 (\ref{P_eps}), (\ref{Q_eps}) はSinkhorn のアルゴリズムによって効率よく近似解を求めることができることが知られています.

参考文献

arxiv.org

*1:総需要量と総供給量が一致していることをあらかじめ仮定しています

*2:行列の  C, P内積は, \begin{equation} \langle C,P \rangle = \sum_{i,j} c_{ij} P_{ij}\end{equation} により定めます

*3:ただし,  0\log{0} = 0 とみなし,  r<0 のとき  r\log{r} = -\infty として考える.

*4: \varepsilon- 緩和ともいう

Python で mod. p での二項係数 (nCk mod. p) の求め方

過去記事, Python で mod の下での逆元を計算するテーブルを作成する - わかばめにっき では, modの下での階乗の逆元を求めました.

本記事では, 以下の記事を参考にして二項係数  {}_nC_k を前処理  O(N), 実行  O(1) で行うアルゴリズムを紹介します.
drken1215.hatenablog.com

続きを読む

AHC001(マラソンマッチ)の参加者の使用言語の分布を調べてみた

競技プログラミングコンテストサイトである AtCoder のマラソン部門のコンテストに参加してみました: atcoder.jp 本記事では、コンテスト参加者の使用しているプログラミング言語はどういうものがあるのかを順位から調査してみました。

続きを読む

ABC192 D - Base n を Python で解いてみる

  • 問題概要

atcoder.jp

数字列  X n 進法表記と見なしたとき,  M 以下になるような値が何種類あるか求めよ.

ただし,  n X に含まれる数字より大きいものとする.

考えたこと

入力例 1

22
10
  • 22 を 3 進数と見なすと,  2\times 3 + 2 = 8 \leqq 10
  • 22 を 4 進数と見なすと,  2\times 4 + 2 = 10 \leqq 10
  • 22 を 5 進数と見なすと,  2\times 5 + 2 = 12 > 10
  • 22 を 6 進数と見なすと,  2\times 6 + 2 = 14 > 10

により, 最初の2つだけが条件を満たします.

このことから,  n が小さいときには条件を満たし,  n を大きくするにつれて  M を超えてしまう閾値を見つければよさそうです.

以下,  X の長さが  K であるとし,  k 桁目を  a_k (   a_k\leqq d)と表します. このとき,  X n(>d) 進法表記と見なすと,

[tex: \displaystyle X_n = \sum_{k=0}^{K-1} nk a_k ]

と表すことができ, これは  X が二桁以上の場合には,  n について単調増大であり, したがっていつかは  M よりも大きくなります.

一方で,  X が一桁の場合には, 何進法表記であっても

 \displaystyle
X_n = a_0

です. いくら  n を大きくしても  M を超えないことがあり得ますが, 問題文で問われているのは  M 以下になるような値が何種類あるかであるため,  X そのものが  M 以下であれば  1 が, そうでなければ  0 が解答になります

解いていく

まずは, 入力を受け取ります. 数字列  X は整数として入力されますが,  n進数展開を行うことを見越してリストとして受け取りましょう.

X = [int(k) for k in list(input())]
d = max(X)
M = int(input())

また, 数字列  X n 進数展開したときの値を p_adic(X, n) により計算します.

def p_adic(X, p):
  curr = 0
  for k in X:
    curr += k
    curr *= p
  curr //= p
  return curr
以下のようなコードだと, TLEしてしまいます

for ループが  O(M) の例 1

ans = 0
for j in range(d+1, M+1):
  if p_adic(X,j) <= M:
    ans += 1

print(ans)

for ループが  O(M) の例 2

ans = 0
for j in range(d+1, M+1):
  if p_adic(X,j) <= M:
    ans = j - d
  else:
    break

print(ans)

今回の制約では  M \leqq 10^{18} のため, より高速なアルゴリズムを構築する必要があります.

二分探索

今回の問題では  n が小さいときには条件を満たし,  n を大きくするにつれて, ある閾値を超えると条件を満たさなくなるという単調性があります.

したがって, その閾値を探す問題に帰着することができます. その際には,

  •  n = M/2 は条件を満たすか?
    • 満たさないならば閾値 0から  M/2の間にある
    • 満たすならば閾値 M/2から  Mの間にある

という性質を使って, 条件の範囲を1ステップごとに半分にしていくことができます. 効率よく探すことで,  O(\log M) の計算量で以上の閾値を求めていきましょう.

一般の二分探索については qiita.com

特に Python の二分探索をうまく書く方法については www.forcia.com

がそれぞれ読みやすく, 参考になります. ポイントとしては,

  • 左端が True, 右端が False となるように初期値 ok, ngを設定する
  • 中央値 (ok + ng)//2 での条件を判定する
    • 左端が True, 右端が False となる状態を維持しながら ok, ngを変更する
  • okng の差が1になったとき, ok の値が「条件を満たす  j の上限」である

ことになります.

解いていく

以上を踏まえて, -  X が一桁の場合がコーナーケースになっていることに注意する - 「  X j 進数展開した結果が $M$ 以下である」をis_ok(j)と表記し, 二分探索で最大の $j$ を求める - 答えとして  j - d を出力する

を行うことにより, ACすることができます

X = [int(k) for k in list(input())]
d = max(X)
M = int(input())

def p_adic(X, p):
  curr = 0
  for k in X:
    curr += k
    curr *= p
  curr //= p
  return curr

def is_ok(j):
  curr = p_adic(X, j)
  if curr <= M:
    return True
  else:
    return False

if len(X) == 1:
  if X[0] <= M:
    print(1)
  else:
    print(0)
  import sys
  sys.exit()

ok = d
ng = M + 1
while ng - ok > 1:
  mid = (ok + ng) // 2 # 平均(小数切り捨て)
  if is_ok(mid):
    ok = mid
  else:
    ng = mid

print(ok-d)

atcoder.jp

SIRモデル方程式を数値的に解いて、gifで表示させてみる

伝染病の拡大過程を記述するSIRモデル

\begin{equation}
\left\{
  \begin{split}
    S′(t)&=−βS(t)I(t),\\
    I′(t)&=βS(t)I(t)−γI(t),\\
    R′(t)&=γI(t)
  \end{split}
\right.
\end{equation}
の初期値問題を Python で解き, 得られた解をgifにします.
f:id:wakabame:20210223221730g:plain

PCとインターネットがあればどなたでも, 以下のノートブックを逐次実行することでgifファイルを生成することができます.
colab.research.google.com

続きを読む

水素原子の基底状態解は原点での値を用いずに特徴付けられる

3次元空間における水素原子ハミルトニアンに対する固有値問題*1

\begin{equation}
-\dfrac{1}{2}\Delta \psi - \dfrac{1}{|x|} \psi = E \psi,\quad x \in \mathbb{R}^3
\end{equation}

について, 1s軌道に対応する基底状態解を考察します.

よく知られているように, 基底状態解は球対称であるため,
 u(|x|) = \psi(x) と変換することにより, 未知関数  u: [0,+\infty)\rightarrow \mathbb{C} についての方程式

\begin{equation}
\left(-\dfrac{1}{2}\dfrac{d^2}{dr^2} -\dfrac{1}{r}\dfrac{d}{dr} -\dfrac{1}{r}\right)u = E u,\quad r \in (0,\infty) \tag{1}
\end{equation}

と変形できます.
さらに, (1) に対して, 形式的(部分積分から出る境界項が消えるよう)に  r^2 \bar{\phi} をかけて部分積分することにより弱形式

\begin{equation}
\left\langle \dfrac{r}{2}\dfrac{du}{dr}, \dfrac{r}{2}\dfrac{d\phi}{dr}\right\rangle
- \left\langle \sqrt{r}u, \sqrt{r}\phi \right\rangle
=E \left\langle ru, r\phi\right\rangle \tag{2}
\end{equation}

を得ます.

本記事では, 1次元における固有値問題 (1) に焦点を当てて, 適当な関数空間での変分問題と見なすことで, 特に陽に境界条件を与えることなく基底状態解が得られることを解説していきます.

また, 本記事は有限要素法による数値計算について紹介している, dc1394さんの記事
qiita.com
についての考察を兼ねています.
モデリング数値計算の箇所を読んでいて, 境界条件を設定していないように見えるのに, 正しく基底状態解が得られていることが気になり調べてみました.
なお, 弱形式の導出などについても紹介記事にて丁寧に行われているため, 気になる方はそちらを参考にしてみてください.

*1:単に, シュレディンガー方程式とも呼びます

続きを読む

変分法の直接法で、ラプラシアンの第一固有関数と第一固有値を手に入れる

この記事では, 滑らかな有界領域  \Omega \subset \mathbb{R}^n に対して Dirichlet Laplacian の第一固有関数  u_1 と第一固有値  \lambda_1 を変分的な手法で得ます.
すなわち, 次の Dirichlet 条件付き楕円型偏微分方程式を満たす固有値  \lambda\in \mathbb{R} と固有関数  u=u(x) を求める問題を考えましょう.

\begin{equation}
\begin{cases}
- \Delta u = \lambda u, &\text{in } \Omega,\\
u = 0, &\text{on } \partial \Omega.
\end{cases}\tag{P}
\end{equation}

簡単のため, 固有値の中でも最小の第一固有値と, 第一固有関数に着目して解いていきます.

レイリー・リッツの特徴づけ

第一固有値  \lambda_1 は次のように特徴づけられる.
\begin{equation}
\lambda_1= \inf\left\{ \frac{\|\nabla v\|^2_2}{\| v\|^2_2} ;v\in H_0^1(\Omega)\setminus\{0\}\right\} \tag{M}
\end{equation}

ただし,  H_0^1(\Omega)ソボレフ空間,  \|\cdot\|_2 L^2 ノルムを表します.
以下では, まずはこの最小化問題が解を持つことを解説していき, 次いでその解がもとの固有値問題の解になっていることを述べます.

続きを読む