最適輸送「モンジュの問題」と「カントロビッチの問題」

最適輸送問題とは?

ある砂山の形  \mu を別の砂山の形  \nu に移すことを考える.
砂山の移動に掛かる労力が「質量」と「始点と終点の位置」で決まっているとき, トータルの労力が最小になるような移し方はどんなものだろうか?

数学的に定式化すると

砂がありえる場所の全体を領域  \Omega \subset \mathbb{R}^d として表し, 砂山の形をそれぞれ確率測度  \mu, \nu \in P(\Omega) で表し, 関数  c:\Omega\times\Omega \rightarrow [0,+\infty] として単位体積を始点から終点へ輸送したときの労力を表す.

モンジュの問題

次の最小化問題を満たす  T を見つけ出せ;
\begin{equation}
\inf\left\{\int c(x,T(x))\,d\mu(x)\,|\, T:\Omega\rightarrow\Omega :\text{可測}, T_{\#}\mu =\nu\right\}.
\end{equation}
ここで,  T_{\#}\mu
\begin{equation}
T_{\#}\mu(A) = \mu(T^{-1}(A)),\quad A\in \mathcal{B}(\Omega)
\end{equation}
で与えられる測度であり,  \mu T による押し出し (push forward) という.

また, この最小化問題の解  T は ( \mu から  \nu への) 最適輸送写像と呼ばれる.

カントロビッチの問題

次の最小化問題を満たす  \gamma を見つけ出せ;
\begin{equation}
\inf\left\{\int\int_{\Omega\times\Omega} c\,d\gamma\,|\, \gamma\in\Pi(\mu,\nu)\right\}.
\end{equation}
ここで,  \Pi(\mu,\nu)
\begin{equation}
\Pi(\mu,\nu) = \left\{ \gamma\in P(\Omega\times\Omega)\,|\, (\pi_1)_{\#}\gamma = \mu, (\pi_2)_{\#}\gamma = \nu\right\}
\end{equation}
により与えられる測度の族であり, この要素を ( \mu から  \nu への) 輸送計画と呼ぶ.
ただしここで,  \pi_1, \pi_2 はそれぞれ  \pi_1(A\times B) = A, \pi_2(A\times B) = B なる射影である.

また, この最小化問題の解  \gamma は ( \mu から  \nu への) 最適輸送計画と呼ばれる.

カントロビッチの問題は割と解ける

関数 c :X\times Y\rightarrow \mathbb{R}\cup \{+\infty\} が下に有界かつ下半連続であるとする.
このとき, カントロビッチの問題は解を持つ.

モンジュの問題が解けるならば, カントロビッチの問題の最適輸送計画はそれ以下のコストを持つ

輸送写像  T がモンジュの問題の解とすると,
\begin{equation}
\gamma_T = (\text{id}, T)_{\#}\mu
\end{equation}
 \gamma_T \in \Pi(\mu,\nu) により,
\begin{align}
\int c(x,T(x))\,d\mu(x)
&= \int\int_{\Omega\times\Omega} c(x,y)\,d\gamma_T(x,y)\\
&\ge \inf\left\{\int\int_{\Omega\times\Omega} c\,d\gamma\,|\, \gamma\in\Pi(\mu,\nu)\right\}.
\end{align}

モンジュの問題は解けないが, カントロビッチの問題は解ける場合がある

これを見るためには,
\begin{align}
\mu &= \delta_0,\\
\nu &= \dfrac{\delta_0 + \delta_1}{2},\\
c(x,y)&= |x-y|
\end{align}
の条件で, それぞれの問題を考察すればよい.
つまり,  0 にあった砂のうちの半分を  1 に輸送することを考える.
このとき, モンジュの問題は  T(0) の値がいずれであっても輸送写像が作れず, 一方でカントロビッチの問題は
\begin{equation}
\gamma(\{0\}\times \{0\}) = \gamma(\{0\}\times\{1\}) = 1/2
\end{equation}
となるようなものを取れば, これが最適輸送計画となる.

無限次元ヒルベルト空間の有界点列は収束部分列を持たないが、弱収束部分列なら持つ

この記事では, 無限次元ヒルベルト空間  H について,

  1. 有界な点列であっても収束するような部分列がない場合があること
  2. 弱収束を考えれば収束するような部分列があること

を述べていきます.

続きを読む

AGC044 の A 問題「Pay to Win」を Python で解いてみる

再帰を書くのが苦手です.
問題はこちら:
atcoder.jp

問題概要

次の選択肢を繰り返すことで初期値  0 N にするための最小コストを求めよ.

  • コスト  A を払うことで, 数値を  2 倍する
  • コスト  B を払うことで, 数値を  3 倍する
  • コスト  C を払うことで, 数値を  5 倍する
  • コスト  D を払うことで, 数値を  1 増やす
  • コスト  D を払うことで, 数値を  1 減らす

ただし, テストケースは  T 個ある

制約


\begin{align*}
&1\le T \le 10 \\
&1\le N\le 10^{18}\\
&1\le A,B,C,D \le 10^9
\end{align*}

続きを読む

L^2 で弱収束するが強収束しない関数列の3パターン

関数列が弱収束するのに強収束しない場合は, この記事で述べるような3パターンの例があります*1.
この記事では, 2乗可積分ルベーグ空間  L^2(\Omega) を対象とし*2,  0 に弱収束する関数の列  \{f_k\}_k \subset L^2(\Omega) の例*3を考えていきます.

例1. めっちゃ振動する

 L^2(0,\pi) における関数列  \{f_k\}_k \subset L^2(0,\pi) を,
\begin{equation}
f_k(x) = \sin{kx},\quad x \in (0, \pi)
\end{equation}
と定めると, 定数関数  0 に弱収束するが, 強収束しない.

例2. どこかに凝集する

 L^2(0,1) における関数列  \{f_k\}_k \subset L^2(0,1) を,
\begin{equation}
f_k(x) =
\begin{cases}
k ,\quad x \in \left(0, \dfrac{1}{k^{2}}\right],\\
0, \quad x \in \left(\dfrac{1}{k^{2}}, 1\right)
\end{cases}
\end{equation}
と定めると, 定数関数  0 に弱収束するが, 強収束しない.

例3. 無限遠へさまよう

 L^2(\mathbb{R}) における関数列  \{f_k\}_k \in L^2(\mathbb{R}) を,
\begin{equation}
f_k(x) =
\begin{cases}
1 ,\quad x \in \left(k, k+1 \right),\\
0, \quad otherwise
\end{cases}
\end{equation}
と定めると, 定数関数  0 に弱収束するが, 強収束しない.

これらの事実について述べます.

*1:もちろん考えている関数空間や領域に依ります

*2:指数  1 < p < \infty に対する  L^p(\Omega) 空間でも同様の例が作れます

*3: 0 でない場合は平行移動すればよい

続きを読む

ABC 162 A - Lucky 7 を Python で解いてみる

問題概要

atcoder.jp

3桁の数字が与えられたとき, 一度でも 7 があったら "Yes" を, そうでないなら "No" を出力せよ.

考えたこと

  • 100 の位の数字が 7 である
  • 10 の位の数字が 7 である
  • 1 の位の数字が 7 である

の3つの条件のいずれかが成り立てば"Yes"を出力すればよいことがわかります.
回答を行うには, 次の 2 つの手順が必要になります:

  • or 演算子を使って 3 つの条件をまとめる
  • if 文を使って条件が成り立つ場合とそうでない場合に分ける

また, 入力を 3 桁の数字としてではなく, 3文字の文字列として扱うことにしましょう.
こうすることで, 3 つの条件は,

  • 1文字目が '7' である
  • 2文字目が '7' である
  • 3文字目が '7' である

のようにそれぞれ書き表すことができます.

入力

入力を文字列として受け取りましょう.

N = input()

これで, 入力を変数  N に 3 文字の文字列として代入することができました.

出力

  • 3つの条件のどれかが成り立つならば, "Yes" を出力する
  • そうでないならば, "No" を出力する

により回答を行います.

if N[0] == '7' or N[1] == '7' or N[2] == '7':
  print("Yes")
else:
  print("No")

Submission #11861851 - AtCoder Beginner Contest 162

一般化パレート分布の平均と分散をレベルセットを使って計算する

積分の計算をレベルセットの積分に変形する公式として, Layer Cake Representation が知られている.

Layer Cake Representation

定理 1. (Layer Cake Representation)

 \nu を非負実数直線  [0,+\infty) 上のBorel 集合族の上の測度とし, 単調増大関
\begin{equation}\label{eq_LCR_1}
\phi(\lambda) = \nu ([0,\lambda))
\end{equation}
が任意の  \lambda>0 に対して有限であるとする. また,  (\Omega,\Sigma, \mu) \sigma -有限な測度空間とし,  f を非負な可測関数とする.
このとき,
\begin{equation}
\int_\Omega \phi(f(x))\,d\mu(x)
= \int_0^\infty \mu\left(x\in \Omega; f(x)>\lambda\right)\,d\nu (\lambda)
\end{equation}
が成り立つ. 特に,  a>0 に対して  \phi(\lambda) = \lambda^a とすると,
\begin{equation}
\int_\Omega f(x)^a\,d\mu(x)
= a\int_0^\infty \lambda^{a-1}\mu\left(x\in \Omega; f(x)>\lambda\right)\,d\lambda
\end{equation}
が成り立つ.

注意 2.

  1.  \mu は符号付き測度でもよい. このとき,  \phi有界変動な関数であり, 第2式は正負のいずれかが有限確定の場合に成り立つ.
  2. 第3式を確率論っぽく書くと,  f の代わりに確率変数 X を用いて,

\begin{equation}
E[X^a]
= \int_\Omega X(\omega)^a\,dp(\omega)
= a\int_0^\infty \lambda^{a-1} P(X>\lambda)\,d\lambda
\end{equation}
と書ける. ただし  P\left(\omega\in \Omega; X(\omega)>\lambda\right) P(X>\lambda) と表している.

続きを読む

「割った余り」での割り算は、既約分数にしなくても答えは変わらない

この記事では,  P素数とし,  \textrm{mod } P の下での割り算について述べます.

整数  a,b に対する  \textrm{mod } P の下での割り算  a\times b^{-1} では,
分数 a/b が既約分数になっていなくても計算結果は変わりません.

そのことを述べていきます.

復習

wakabame.hatenablog.com
でも触れましたが, まずは逆元とは何かを述べます.
 P の倍数ではない  a=0,1, \cdots に対して,  a の法  P の下での逆元  a^{-1} を求めましょう.
まず,

 x =0,1, \cdots , P-1  a の法  P の下での逆元であるとは,
ある  y_a \in \mathbb{N} が存在して,
\begin{equation}
ax+Py_a =1
\end{equation}
を満たすこと

と定式化できたのでした.
そのような  x a P の倍数でない限り一意に存在し, 以後これを  a^{-1} と表します.

このような  a^{-1} がなぜとれるかについては,
拡張ユークリッドの互除法について述べられている @drken さんの記事
qiita.com
を参考にしてみてください.

示すべきこと

 a\times b^{-1} \textrm{mod } P は既約でなくても計算結果が変わらないことを示すために,

 P 未満の自然数  k に対して,  ka\times (kb)^{-1} を計算して, その差を観察しましょう.

まずは, 拡張ユークリッドの互除法から  b^{-1},  (kb)^{-1} はある整数  y_b, y_{kb} を用いて,
\begin{equation}
b^{-1} = \dfrac{1-Py_b}{b},\quad (kb)^{-1} = \dfrac{1-Py_{kb}}{kb}
\end{equation}
と表せたことに注意します.

これにより,
\begin{align*}
ab^{-1} - ka(kb)^{-1}
&= a\dfrac{1-Py_b}{b} -ka \dfrac{1-Py_{kb}}{kb}\\
&= \dfrac{a(y_{kb}-y_{b})}{b}P
\end{align*}

を得ます.

ここで, 左辺は整数だから右辺も整数であり,  b P の倍数ではないので
 \dfrac{a(y_{kb}-y_{b})}{b}
は整数です.

したがって両辺は  P の倍数であることから,
\begin{equation}
ab^{-1} \textrm{mod } P = ka(kb)^{-1} \textrm{mod } P
\end{equation}

が従います.

別の説明

 kb についての逆元は,
\begin{equation}
(kb)^{-1} \equiv b^{-1}k^{-1} \pmod P
\end{equation}
を満たします.

これにより,
\begin{equation}
ka(kb)^{-1} \equiv kab^{-1}k^{-1} \equiv ab kk^{-1} \equiv ab\pmod P
\end{equation}
を示すこともできます(指摘されるまで気が付きませんでいた。。).