[計算ノート] 圧縮センシングのレプリカ対称解

圧縮センシング (compressed sensing) をレプリカ法の処方箋に従って解析したときの計算ノートです。

問題設定

まず各値を以下のように定める。

• $N$: 信号の次元
• $M$: 観測データの個数
• $K$: 原信号の非ゼロ要素の数

ただし

1. $M \lt N$
2. $K \ll N$ 、つまり原信号はゼロの要素の占める割合がかなり大きいスパースな信号である
3. $\rho = K / N$ は非零要素の割合を示す
4. $N$ は非常に大きい

と仮定する。このような仮定のもと、以下のような設定を考える。

• $\bm x^{(0)} \in \mathbb R^N$: 真の信号 (未知)
• $\bm A \in \mathbb R^{M \times N}$: 観測行列 (既知)
• $\bm y = \bm A \bm x^{(0)} \in \mathbb R^M$: 観測データ (既知)

信号と観測行列について以下のような性質を定める。

1. 真の信号は以下の分布に独立同分布で従う:

$$\begin{aligned} P(x^{(0)}_i) &= (1 - \rho) \delta(x^{(0)}_i) + \rho \mathcal N(x^{(0)}_i \mid 0, 1) \end{aligned}$$

2. 観測行列は以下の分布に独立同分布で従う:

$$\begin{aligned} P(A_{ij}) &= \mathcal N\left( A_{ij} \middle| 0, \frac{1}{N} \right) \end{aligned}$$

3. 観測ノイズは限りなく小さい。すなわち $\sigma^2 \to 0$ の極限を考える

以上のような設定のもと、

$$\begin{aligned} \bm y = \bm A \bm x \end{aligned}$$

を満たすような $\bm x$ を求めるという問題を考える。ただし劣決定系 $M < N$ を仮定する。線形方程式を解きたいが、方程式の個数 $M$ が自由度 $N$ よりも少ないので、一般には解が複数存在するという状況である。

圧縮センシングでは、真の解がスパースだと信じ、できるだけゼロの多い解を求めようとする。本ノートではL1ノルムを最小化することでこれを達成する場合を考える。具体的には、最小化問題

$$\begin{aligned} & \operatorname*{minimize}_{\bm x} && \|\bm x\|_1 \\ & \operatorname{subject~to} && \bm y = \bm A \bm x \end{aligned}$$

の解 $\bm x^\ast$ が、真の信号 $\bm x^{(0)}$ に十分に一致するのはどのような場合か、ということを考える。

この最小化問題の目的関数を $\|\bm x\|_1$ とする代わりに、もっと性質の良い関数を使うこともできるが、本ノートではL1ノルムを使うバージョンについて議論する。

罰金法による定式化とBoltzmann分布

等式制約 $\bm y = \bm A \bm x$ を罰金法によって目的関数に組み込むことで、以下のような最小化問題を考えることができる。

$$\begin{aligned} & \operatorname*{minimize}_{\bm x} && H(\bm x) = \|\bm x\|_1 + \frac{\lambda}{2} \|\bm y - \bm A \bm x\|_2^2. \end{aligned}$$

$\lambda > 0$ は罰金項の重みを調整するパラメータである。このパラメータを十分に大きくすれば、等式制約が満たされるようになる。のちほど $\lambda \to \infty$ という極限を取る。以上の目的関数 $H(\bm x)$ を「ハミルトニアン」とみなし、以下のようなBoltzmann分布を考える。

$$\begin{aligned} P(\bm x \mid \bm A, \bm y, \lambda) &= \frac{1}{Z(\bm A, \bm y)} \exp\left( -\beta H(\bm x) \right) \\ &= \frac{1}{Z(\bm A, \bm y)} \exp\left( - \beta \|\bm x\|_1 - \frac{\lambda}{2} \|\bm y - \bm A \bm x\|_2^2 \right) \end{aligned}$$

$\beta > 0$ は「逆温度」とよばれるパラメータで、系が「熱平衡状態」に達したときに $\| \bm x \|_1$ がどれほど最小化されるかを調整する効果がある。$\beta \to \infty$ で目的関数が最小の状態に収束していく。

Boltzmann分布にエネルギー関数を埋め込むにあたって、罰金係数 $\lambda$ を $\beta$ で割った。これには $\beta$ の大小 (=最小化問題の目的関数の値の小ささ) にかかわらず、制約条件である $\bm y = \bm A \bm x$ は満たされていなくてはならないという気持ちが込められている。

正規化定数 $Z(\bm A, \bm y)$ は

$$\begin{aligned} Z(\bm A, \bm y) &= \int d \bm x \exp\left( - \beta \|\bm x\|_1 - \frac{\lambda}{2} \|\bm y - \bm A \bm x\|_2^2 \right) \end{aligned}$$

で定義される量で、統計力学の文脈では分配関数と呼ばれる。

Bayes推定的なアプローチ

前節で定式化された分配関数は、以下のようにBayes推定の枠組みを使って尤度関数と事前分布を組み合わせて導出することもできる。

尤度関数

$\bm x \in \mathbb R^N$ に関する尤度関数を以下のように定める。

$$\begin{aligned} P(\bm y \mid \bm A, \bm x, \sigma) &= \mathcal N(\bm y \mid \bm A \bm x, \sigma^2 \bm I_M) \\ &= \frac{1}{\sqrt{(2 \pi \sigma^2)^M}} \exp\left( -\frac{1}{2\sigma^2} \|\bm y - \bm A \bm x\|_2^2 \right). \end{aligned}$$

この式は、信号 $\bm x$ を観測行列 $\bm A$ によって観測する際に、何らかの観測ノイズ $\bm \varepsilon \sim \mathcal N(\bm 0, \sigma^2 \bm I_M)$ が加わったものが $\bm y$ として得られる状況をモデル化している。

$$\begin{aligned} \bm y &= \bm A \bm x + \bm \varepsilon, \\ \bm \varepsilon &\sim \mathcal N(\bm 0, \sigma^2 \bm I_M). \end{aligned}$$

新たに追加されたパラメータ $\sigma$ は観測ノイズの強さを表す。あとで $\sigma \to 0$ とすることで、ノイズがない状況へと近づけていく。

事前分布

Bayes推定の枠組みでは、推定すべきパラメータ $\bm x$ について何らかの特徴を仮定し、事前分布として定式化する。本ノートでは、$\bm x$ がスパースであるという特徴を仮定するためにLaplace分布を事前分布として使うことにする。Laplace事前分布を以下のように定義する。

$$\begin{aligned} P(\bm x) &= \mathrm{Laplace}(\bm x \mid \bm 0, \beta^{-1}) \\ &= \left( \frac{\beta}{2} \right)^N \exp\left( -\beta \|\bm x\|_1 \right). \end{aligned}$$

パラメータ $\beta > 0$ はL1ノルムの圧縮の程度を調整する。$\beta \to \infty$ とすると、所与の条件下で $\| \bm x \|_1$ が最小の解へと収束していく。

事後分布は

$$\begin{aligned} P(\bm x \mid \bm A, \bm y, \sigma) &\propto P(\bm y \mid \bm A, \bm x, \sigma) P(\bm x) \\ &\propto \frac{1}{Z(\bm A, \bm y)} \exp\left( - \beta \|\bm x\|_1 - \frac{1}{2\sigma^2} \|\bm y - \bm A \bm x\|_2^2 \right) \end{aligned}$$

となる。前節で定義したBoltzmann分布を $\lambda \mapsto 1 / \sigma^2$ と書き換えた形をしている。$\sigma \to 0$ としてノイズなしの状況に近づけることは、前節の定式化において $\lambda \to \infty$ として制約条件を厳密に満たすことに対応する。

前節と同様に、分配関数が

$$\begin{aligned} Z(\bm A, \bm y) &= \int d \bm x \exp\left( - \beta \|\bm x\|_1 - \frac{\lambda}{2} \|\bm y - \bm A \bm x\|_2^2 \right) \end{aligned}$$

と定義される。

レプリカ解析

ここまででハミルトニアン、Boltzmann分布、分配関数を以下のように定めた。

$$\begin{aligned} H(\bm x) &= \|\bm x\|_1 + \frac{\lambda}{2 \beta} \|\bm y - \bm A \bm x\|_2^2, \\ P(\bm x \mid \bm A, \bm y, \beta) &= \frac{1}{Z(\bm A, \bm y)} \exp\left( - \beta H(\bm x) \right), \\ Z(\bm A, \bm y) &= \int d \bm x~ \exp\left( - \beta H(\bm x) \right). \end{aligned}$$

これは復元信号 $\bm x$ が逆温度 $\beta$ のBoltzmann分布に従う系を考えていることに対応する。一方で観測行列 $\bm A$ と観測値 $\bm y$ に関しては、別の逆温度 $\tilde \beta \sim \mathcal O(\beta)$ に従うと考え、以下のような系を考える。

$$\begin{aligned} \tilde H(\bm A, \bm y) &= - \frac{1}{\beta} \log Z(\bm A, \bm y) - \frac{1}{\tilde \beta} \log P(\bm A) P(\bm y) , \\ \tilde P(\bm A, \bm y \mid \tilde \beta) &= \frac{1}{\tilde Z} \exp\left( - \tilde \beta \tilde H(\bm A, \bm y) \right) \\ &= \frac{1}{\tilde Z} Z^{\tilde \beta / \beta}(\bm A, \bm y) P(\bm A) P(\bm y) , \\ \tilde Z &= \int d \bm A~ \int d \bm y~ \exp\left( - \tilde \beta \tilde H(\bm A, \bm y) \right) \\ &= \mathbb E_{\bm A, \bm y}~ Z^{\tilde \beta / \beta}(\bm A, \bm y). \end{aligned}$$

ただし期待値は $\bm A \sim P(\bm A)$ と $\bm y \sim P(\bm y)$ によるものである。「配位平均」などとも呼ばれる。この系に対して、自由エネルギーという量が以下のように定義される。

$$\begin{aligned} \mathcal F(\beta, \tilde \beta) &= -\frac{1}{\tilde \beta} \log \tilde Z \\ &= -\frac{1}{\tilde \beta} \mathbb E_{\bm A, \bm y}~ Z^{\tilde \beta / \beta}(\bm A, \bm y). \end{aligned}$$

自由エネルギーは系の状態の安定さを表す量である。Helmholtzの自由エネルギー最小の原理によれば、系が熱平衡状態にあるときには自由エネルギーが最小値を取る。$\beta \to \infty$ の極限において自由エネルギーが最小値を取るときの $\bm x$ が、最小化問題の最適解に対応する。

更に系の自由度 $N$ で割れば、1自由度あたりの自由エネルギーが得られる。

$$\begin{aligned} f(\beta, \tilde \beta) &= -\frac{1}{N \tilde \beta} \mathbb E_{\bm A, \bm y}~ Z^{\tilde \beta / \beta}(\bm A, \bm y). \end{aligned}$$

以降はこの自由エネルギーが解析の対象となる。

クエンチ系のレプリカ解析

今考えている系では、観測行列 $\bm A$ と観測値 $\bm y$ が

$$\begin{aligned} \tilde P(\bm A, \bm y \mid \tilde \beta) &= \frac{1}{Z(\bm A, \bm y)} \exp\left( \frac{\tilde \beta}{\beta} \log Z(\bm A, \bm y) \right) \times P(\bm A) P(\bm y) , \end{aligned}$$

なる分布に従うことになっている。逆温度の比 $\tilde \beta / \beta$ は $\bm A$, $\bm y$ がBoltzmann分布

$$\begin{aligned} \frac{1}{Z(\bm A, \bm y)} \exp\left( \frac{\tilde \beta}{\beta} \log Z(\bm A, \bm y) \right) \end{aligned}$$

にどれほど依存するかを表す特徴的な量であり、

$$\begin{aligned} n \coloneqq \frac{\tilde \beta}{\beta} \end{aligned}$$

と表される。

今回の設定では観測行列 $\bm A$, 観測ベクトル $\bm y$ はランダムに固定されているべきだから、$n \to +0$ 極限を取るのが適当である。このような系は「クエンチ系」と呼ばれ、そのような配位平均 $\mathbb E_{\bm A, \bm y}$ の平均の取り方は「クエンチ平均」と呼ばれる^[1]。

$$\begin{aligned} f &= - \frac{1}{N \beta} \lim_{n \to +0} \frac{1}{n} \mathbb E_{\bm A, \bm y}~ Z^n(\bm A, \bm y). \end{aligned}$$

これを計算するために以下のことをする。

1. $n$ を正の整数だと思い込む。
2. 系の $n$ 個のコピー (レプリカ) $\bm x^{(1)}, \ldots, \bm x^{(n)}$ に対応する分配関数 $Z^{(1)}, \ldots, Z^{(n)}$ の積 $\prod_a Z^{(a)}$ の期待値 $\mathbb E_{\bm A, \bm y}~Z^n(\bm A, \bm y)$ を計算する。
3. $n$ を実数に解析接続し、$n \to +0$ の極限を取る。

この手法をレプリカ解析という。解析接続の正当性は一般に保証されないから、結果の解釈には注意が必要である。

分配関数のn乗の配位平均の計算

分配関数の $n$ 乗は

$$\begin{aligned} Z^n(\bm A, \bm y) &= \left\{ \int d \bm x \exp\left( -\beta \| \bm x \|_1 - \frac{\lambda}{2} \| \bm y - \bm A \bm x \|_2^2 \right) \right\}^n \\ &= \int d \bm x^{(1)} \cdots \int d \bm x^{(n)} \prod_{a=1}^n \exp\left( -\beta \|\bm x^{(a)}\|_1 - \frac{\lambda}{2} \|\bm y - \bm A \bm x^{(a)}\|_2^2 \right) \\ &= \int \prod_{a=1}^n d \bm x^{(a)} \exp\left( -\beta \sum_{a=1}^n \|\bm x^{(a)}\|_1 - \frac{\lambda}{2} \sum_{a=1}^n \|\bm y - \bm A \bm x^{(a)}\|_2^2 \right) \end{aligned}$$

となる。ただし $n$ 重積分を

$$\begin{aligned} \int \prod_{a=1}^n d \bm x^{(a)} &= \int d \bm x^{(1)} \cdots \int d \bm x^{(n)} \end{aligned}$$

と略記する。配位平均は

$$\begin{aligned} & \mathbb E_{\bm A, \bm y}~ Z^n(\bm A, \bm y) \\ &= \mathbb E_{\bm A, \bm y}~ \int \prod_{a=1}^n d \bm x^{(a)}~ \exp\left( -\beta \sum_{a=1}^n \|\bm x^{(a)}\|_1 - \frac{\lambda}{2} \sum_{a=1}^n \|\bm y - \bm A \bm x^{(a)}\|_2^2 \right) \end{aligned}$$

である。以降はこれの計算を進めていく。

配位平均の変形

配位平均はそのままだと扱いづらいので、別の形に変形していく。

原信号に対してy=Axの適用

まず $\bm y = \bm A \bm x^{(0)}$ を適用するために、$\bm y$ の分布を

$$\begin{aligned} P(\bm y) &= \int d \bm x^{(0)} \delta(\bm y - \bm A \bm x^{(0)}) P(\bm x^{(0)}) \end{aligned}$$

と定義する。配位平均はDiracのデルタ関数の性質 $\int d a \delta(a - b) f(a) = \int db f(b)$ を用いて以下のように書ける。

$$\begin{aligned} &\mathbb E_{\bm A, \bm y}~ Z^n(\bm A, \bm y) \\ &= \mathbb E_{\bm A, \bm y}~ \int \prod_{a=1}^n d \bm x^{(a)}~ \exp\left( -\beta \sum_{a=1}^n \|\bm x^{(a)}\|_1 - \frac{\lambda}{2} \sum_{a=1}^n \|\bm y - \bm A \bm x^{(a)}\|_2^2 \right) \\ &= \mathbb E_{\bm A, \bm x^{(0)}}~ \int \prod_{a=1}^n d \bm x^{(a)}~ \exp\left( -\beta \sum_{a=1}^n \|\bm x^{(a)}\|_1 - \frac{\lambda}{2} \sum_{a=1}^n \|\bm A (\bm x^{(0)} - \bm x^{(a)})\|_2^2 \right). \end{aligned}$$

Aに関する平均

$\bm A$ の分布は

$$\begin{aligned} P(\bm A) &= \prod_{i=1}^M \prod_{j=1}^N\mathcal N\left( A_{ij} \middle| 0, \frac{1}{N} \right) \end{aligned}$$

である。分配関数の式中では $\bm A$ は

$$\begin{aligned} \bm t^{(a)} &= \bm A (\bm x^{(0)} - \bm x^{(a)}) \\ &= \left[ \sum_{j=1}^N A_{ij} (x^{(a)}_j - x^{(0)}_j) \right]_{i=1}^M \in \mathbb R^M \end{aligned}$$

の形で現れる。そこで $P(\bm A)$ の代わりに、$t^{(a)}_i$ をあらゆる $(a, i)$ に関して並べた行列

$$\begin{aligned} \bm T \coloneqq \left[ t^{(a)}_i \right]_{a, i} \in \mathbb R^{n \times M} \end{aligned}$$

の確率分布 $P(\bm T \mid \{\bm x^{(a)}\}, \bm x^{(0)})$ を考える。配位平均は次のように変形される。

$$\begin{aligned} & \mathbb E_{\bm A, \bm y}~ Z^n(\bm A, \bm y) \\ &= \mathbb E_{\bm A, \bm x^{(0)}}~ \int \prod_{a=1}^n d \bm x^{(a)}~ \exp\left( -\beta \sum_{a=1}^n \|\bm x^{(a)}\|_1 - \frac{\lambda}{2} \sum_{a=1}^n \|\bm A \bm x^{(0)} - \bm A \bm x^{(a)}\|_2^2 \right) \\ &= \mathbb E_{\bm x^{(0)}}~ \int \prod_{a=1}^n d \bm x^{(a)}~ \mathbb E_{\bm T \mid \{\bm x^{(a)}\}, \bm x^{(0)}}~ \exp\left( -\beta \sum_{a=1}^n \|\bm x^{(a)}\|_1 - \frac{\lambda}{2} \| \bm T \|_F^2 \right). \end{aligned}$$

まだ $\bm T$ は $\bm x^{(a)}$ と $\bm x^{(0)}$ に依存しているので、期待値を取る順序は変更できないことに注意する。

Tの分布の考察

$P(\bm T \mid \{\bm x^{(a)}\}, \bm x^{(0)})$ が何者かを考える。$A_{ij}$ は正規分布に従うので、その線形結合である $t^{(a)}_i$ も何らかの多変量正規分布に従う。$A_{ij}$ は $i$ と $j$ に関して独立であることに注意しながら、$\bm T$ の1次と2次のモーメントを計算する。

$$\begin{aligned} \mathbb E_{\bm A} ~ t^{(a)}_i &= \mathbb E_{\bm A} \left[ \sum_{j=1}^N A_{ij} (x^{(0)}_j - x^{(a)}_j) \right] \\ &= 0, \\ \mathbb E_{\bm A} ~ t^{(a)}_i t^{(b)}_j &= \mathbb E_{\bm A} ~ \sum_{k=1}^N A_{ik} (x^{(0)}_k - x^{(a)}_k) \sum_{l=1}^N A_{jl} (x^{(0)}_l - x^{(b)}_l) \\ &= \sum_{k=1}^N \sum_{l=1}^N \mathbb E_{\bm A} \left[ A_{ik} A_{jl} \right] (x^{(0)}_k - x^{(a)}_k) (x^{(0)}_l - x^{(b)}_l) \\ &= \sum_{k=1}^N \sum_{l=1}^N \frac{ \delta_{ij} \delta_{kl} }{N} (x^{(0)}_k - x^{(a)}_k) (x^{(0)}_l - x^{(b)}_l) \\ &= \frac{\delta_{ij}}{N} \left( \bm x^{(0) \top} \bm x^{(0)} - \bm x^{(0) \top} \bm x^{(a)} - \bm x^{(0) \top} \bm x^{(b)} + \bm x^{(a) \top} \bm x^{(b)} \right). \end{aligned}$$

よって $t^{(a)}_i$ には

• データ次元 $i \in \{1, \dots, M \}$ に関しては独立同分布に従い、
• レプリカ番号 $a \in \{ 1, \dots, n \}$ に関しては相関を持つ

という特徴があるので、

$$\begin{aligned} P(\bm T \mid \{\bm x^{(a)}\}, \bm x^{(0)}) &= \prod_{i=1}^M P( \bm t^{(:)}_i \mid \{ \bm x^{(a)} \}, \bm x^{(0)} ) \end{aligned}$$

と分解できる。ただし

$$\begin{aligned} \bm t^{(:)}_i \coloneqq \left[ t^{(a)}_i \right]_a \in \mathbb R^n \end{aligned}$$

と定める。そしてレプリカ番号 $a$ に関する共分散行列を $\bm \Sigma \in \mathbb R^{n \times n}$ と書けば、その $(a, b)$ 成分は

$$\begin{aligned} \Sigma_{ab} = \frac{1}{N} \left( \bm x^{(0) \top} \bm x^{(0)} - \bm x^{(0) \top} \bm x^{(a)} - \bm x^{(0) \top} \bm x^{(b)} + \bm x^{(a) \top} \bm x^{(b)} \right). \end{aligned}$$

$\bm T$ の確率分布は最終的に次式のようになる。

$$\begin{aligned} P(\bm T \mid \{\bm x^{(a)}\}, \bm x^{(0)}) &= \prod_{i=1}^M \mathcal N( \bm t^{(:)}_i \mid \bm 0, \bm \Sigma ) \\ &= \prod_{i=1}^M \frac{1}{\sqrt{(2\pi)^n \det(\bm \Sigma)}} \exp\left( -\frac{1}{2} \bm t^{(:) \top}_i \bm \Sigma^{-1} \bm t^{(:)}_i \right). \end{aligned}$$

秩序変数の導入によるTの切り離し

このままでは $\bm T$ は復元結果 $\{\bm x^{(a)}\}$ と原信号 $\bm x^{(0)}$ に依存しており、分配関数の積分の順序が変更できなかったり、項を分割できなかったりするなど都合が悪い。これを解消するために以下の量を導入する。

$$\begin{aligned} q^{00} &= \frac{1}{N} \bm x^{(0) \top} \bm x^{(0)}, \\ q^{0a} &= \frac{1}{N} \bm x^{(0) \top} \bm x^{(a)}, \\ q^{ab} &= \frac{1}{N} \bm x^{(a) \top} \bm x^{(b)}. \end{aligned}$$

$q^{ab} = \bm x^{(a) \top} \bm x^{(b)} / N$ という情報はデルタ関数

$$\begin{aligned} P(\bm Q \mid \{ \bm x^{(a)} \}, \bm x^{(0)}) = \prod_{a=0}^n \prod_{b=0}^n \delta \left( q^{ab} - \frac{1}{N} \bm x^{(a) \top} \bm x^{(b)} \right) \end{aligned}$$

で表現し、$q^{(ab)}$ は $P(\bm Q \mid \{ \bm x^{(a)} \}, \bm x^{(0)})$ に従う新たな確率変数として導入する。こうすると $P(\bm T \mid \{\bm x^{(a)}\}, \bm x^{(0)})$ は $P(\bm T \mid \bm Q)$ に書き換えられる。分配関数は以下のように書ける。

$$\begin{aligned} & \mathbb E_{\bm A, \bm y}~ Z^n(\bm A, \bm y) \\ &= \mathbb E_{\bm x^{(0)}}~ \int \prod_{a=1}^n d \bm x^{(a)}~ \mathbb E_{\bm T \mid \{\bm x^{(a)}\}, \bm x^{(0)}}~ \exp\left( -\beta \sum_{a=1}^n \|\bm x^{(a)}\|_1 - \frac{\lambda}{2} \| \bm T \|_F^2 \right) \\ &= \mathbb E_{\bm x^{(0)}}~ \int \prod_{a=1}^n d \bm x^{(a)}~ \mathbb E_{\bm Q \mid \{\bm x^{(a)}\}, \bm x^{(0)}}~ \mathbb E_{\bm T \mid \bm Q}~ \exp\left( -\beta \sum_{a=1}^n \|\bm x^{(a)}\|_1 - \frac{\lambda}{2} \| \bm T \|_F^2 \right) \\ &= \int d \bm Q \mathbb E_{\bm T \mid \bm Q} \left[ \exp\left( - \frac{\lambda}{2} \| \bm T \|_F^2 \right) \right] \\ & \qquad \qquad \times \mathbb E_{\bm x^{(0)}} \left[ \int \prod_{a=1}^n d \bm x^{(a)} \exp\left( -\beta \sum_{a=1}^n \|\bm x^{(a)}\|_1 \right) P(\bm Q \mid \{ \bm x^{(a)} \}, \bm x^{(0)}) \right] . \end{aligned}$$

$q^{ab}$ は2つのレプリカ $a$, $b$ の最適化問題の解 $\bm x^{(a)}$, $\bm x^{(b)}$ の一致度を表しており、レプリカ番号 $a$ に関する秩序の程度を表すことから秩序変数とかオーバーラップとか呼ばれる。更に自己平均性の仮定のもとでは $q^{00} = \bm x^{(0) \top} \bm x^{(0)} / N$ は $N \to \infty$ において期待値 $\mathbb E[\bm x^{(0) \top} \bm x^{(0)} / N]$ に一致するので

$$\begin{aligned} q^{00} &= \frac{1}{N} \mathbb E_{\bm x^{(0)}} \left[ \| \bm x^{(0)} \|_2^2 \right] \\ &= \frac{1}{N} \sum_{i=1}^N \mathbb E_{x^{(0)}_i} \left[ (x^{(0)}_i)^2 \right] \\ &= \rho \end{aligned}$$

である。以上の秩序変数を用いれば

$$\begin{aligned} \Sigma_{ab} &= q^{00} - q^{0a} - q^{0b} + q^{ab} \\ &= \rho - q^{0a} - q^{0b} + q^{ab} \end{aligned}$$

と書き直せる。

内部エネルギーとエントロピーへの分解

統計力学的には、分配関数の前半の $\bm T$ に関する期待値は「内部エネルギー」、後半の $\bm x^{(0)}$ に関する期待値は「エントロピー」と解釈される、らしい。「内部エネルギー」を

$$\begin{aligned} e(\bm Q) &\coloneqq -\frac{1}{N} \log \mathbb E_{\bm T \mid \bm Q}~ \exp \left( -\frac{\lambda}{2} \| \bm T \|_F^2 \right), \end{aligned}$$

「エントロピー」を

$$\begin{aligned} s(\bm Q) &\coloneqq \frac{1}{N} \log \mathbb E_{\bm x^{(0)}} \int \prod_{a=1}^n d \bm x^{(a)}~ \exp\left( -\beta \sum_{a=1}^n \|\bm x^{(a)}\|_1 \right) P(\bm Q \mid \{ \bm x^{(a)} \}, \bm x^{(0)}) \end{aligned}$$

と定義する。これらの量を用いて分配関数は

$$\begin{aligned} \mathbb E_{\bm A, \bm y}[Z^n(\bm A, \bm y)] &= \int d \bm Q \exp\left( - N e(\bm Q) + N s(\bm Q) \right) \end{aligned}$$

と書ける。

ここまでのまとめ

以上、

1. 観測値の分布の変形
2. 観測行列の分布の変形
3. 秩序変数の導入
4. 内部エネルギーとエントロピーへの分解

という段階を踏んで、

$$\begin{aligned} & \mathbb E_{\bm A, \bm y}~ Z^n(\bm A, \bm y) \\ &= \mathbb E_{\bm A, \bm y}~ \int \prod_{a=1}^n d \bm x^{(a)}~ \exp\left( -\beta \sum_{a=1}^n \|\bm x^{(a)}\|_1 - \frac{\lambda}{2} \sum_{a=1}^n \|\bm y - \bm A \bm x^{(a)}\|_2^2 \right) \\ &= \mathbb E_{\bm A, \bm x^{(0)}}~ \int \prod_{a=1}^n d \bm x^{(a)}~ \exp\left( -\beta \sum_{a=1}^n \|\bm x^{(a)}\|_1 - \frac{\lambda}{2} \sum_{a=1}^n \|\bm A (\bm x^{(0)} - \bm x^{(a)})\|_2^2 \right) \\ &= \mathbb E_{\bm x^{(0)}}~ \int \prod_{a=1}^n d \bm x^{(a)}~ \mathbb E_{\bm T \mid \{\bm x^{(a)}\}, \bm x^{(0)}}~ \exp\left( -\beta \sum_{a=1}^n \|\bm x^{(a)}\|_1 - \frac{\lambda}{2} \| \bm T \|_F^2 \right) \\ &= \mathbb E_{\bm x^{(0)}}~ \int \prod_{a=1}^n d \bm x^{(a)}~ \mathbb E_{\bm Q \mid \{\bm x^{(a)}\}, \bm x^{(0)}}~ \mathbb E_{\bm T \mid \bm Q}~ \exp\left( -\beta \sum_{a=1}^n \|\bm x^{(a)}\|_1 - \frac{\lambda}{2} \| \bm T \|_F^2 \right) \\ &= \int d \bm Q~ \mathbb E_{\bm T \mid \bm Q} \left[ \exp\left( - \frac{\lambda}{2} \| \bm T \|_F^2 \right) \right] \times \mathbb E_{\bm x^{(0)}} \left[ \int \prod_{a=1}^n d \bm x^{(a)}~ \exp\left( -\beta \sum_{a=1}^n \|\bm x^{(a)}\|_1 \right) P(\bm Q \mid \{ \bm x^{(a)} \}, \bm x^{(0)}) \right] \\ &= \int d \bm Q~ \exp\left( - N e(\bm Q) + N s(\bm Q) \right) \end{aligned}$$

と変形し、秩序変数に関する積分へと変形してきた。

レプリカ解析に慣れると、ここまでの計算は暗算できるようになるらしい。

内部エネルギーの計算

内部エネルギー

$$\begin{aligned} e(\bm Q) &\coloneqq -\frac{1}{N} \log \mathbb E_{\bm T \mid \bm Q}~ \exp \left( -\frac{\lambda}{2} \| \bm T \|_F^2 \right) \end{aligned}$$

を多次元Gauss積分

$$\begin{aligned} \int d \bm x \exp\left( - \frac{1}{2} \bm x^\top \bm \Sigma^{-1} \bm x \right) &= \sqrt{(2\pi)^n \det \bm \Sigma} \end{aligned}$$

を使って変形する。

$$\begin{aligned} e(\bm Q) &= -\frac{1}{N} \log \int d \bm T P(\bm T) \exp \left( -\frac{\lambda}{2} \| \bm T \|_F^2 \right) \\ &= -\frac{1}{N} \log \prod_{i=1}^M \int d \bm t^{(:)}_i P(\bm t^{(:)}_i) \exp \left( -\frac{\lambda}{2} \| \bm t^{(:)}_i \|_2^2 \right) \\ &= -\frac{M}{N} \log \int d \bm t^{(:)} \frac{1}{\sqrt{(2 \pi)^n \det \bm \Sigma}} \exp \left( - \frac{1}{2} \bm t^{(:) \top} \left( \bm \Sigma^{-1} + \lambda \bm I_n \right) \bm t^{(:)} \right) \\ &= -\alpha \log \sqrt{\frac{ \det \left( \bm \Sigma^{-1} + \lambda \bm I_n \right)^{-1} }{ \det \bm \Sigma} } \\ &= - \frac{\alpha}{2} \log \frac{ 1 }{ \det \bm \Sigma \det (\bm \Sigma^{-1} + \lambda \bm I_n) } \\ &= \frac{\alpha}{2} \log \det \left( \bm I_n + \lambda \bm \Sigma \right). \end{aligned}$$

途中で $\alpha = M / N$ を用いた。

レプリカ対称性の仮定

行列 $\bm I_n + \lambda \bm \Sigma$ は

$$\begin{aligned} \bm I_n + \lambda \bm \Sigma &= \left[ \delta_{a, b} + \lambda \Sigma_{ab} \right]_{a, b} \\ &= \left[ \delta_{a, b} + \lambda \left( \rho - q^{0a} - q^{0b} + q^{ab} \right) \right]_{a, b} \in \mathbb R^{n \times n} \end{aligned}$$

である。ここに、レプリカ対称性と呼ばれる仮定をおく。

$$\begin{alignedat}{2} &(1) \quad& q^{aa} &= Q, \\ &(2) \quad& q^{ab} &= q \quad (a \neq b), \\ &(3) \quad& q^{0a} &= m, \\ \end{alignedat}$$

以下のような意味である。

1. どのレプリカにおいても、復元後の信号のL2ノルムは一定である。
2. どの2つのレプリカを取ってきても、同じくらい似ているし、同じくらい似ていない。
3. どのレプリカにおいても、原信号と復元後の信号のオーバーラップは等しい。

行列 $\bm Q$ と共分散行列 $\bm \Sigma$ は、以下のように対称性をもつことになる。

$$\begin{aligned} \bm Q &= \begin{bmatrix} Q & q & \cdots & q \\ q & Q & \cdots & q \\ \vdots & \vdots & \ddots & \vdots \\ q & q & \cdots & Q \end{bmatrix}, \\ \bm \Sigma &= \begin{bmatrix} \rho - 2m + Q & \rho - 2m + q & \cdots & \rho - 2m + q \\ \rho - 2m + q & \rho - 2m + Q & \cdots & \rho - 2m + q \\ \vdots & \vdots & \ddots & \vdots \\ \rho - 2m + q & \rho - 2m + q & \cdots & \rho - 2m + Q \end{bmatrix} \end{aligned}$$

以降は $q^{ab}$ の代わりに $Q, q, m$ を用いて計算を進める。

行列式の計算

行列は以下のように書き換えられる。

$$\begin{aligned} \bm I_n + \lambda \bm \Sigma &= \left( 1 + \lambda (Q - q) \right) \bm I_n + \lambda \left( \rho - 2m + q \right) \bm 1_{n \times n}. \end{aligned}$$

$\bm 1_{n \times n}$ は $n \times n$ の全要素が $1$ の行列である。各項の固有値分解は

$$\begin{aligned} \bm I_n &= \sum_{i=1}^n 1 \bm u_i \bm u_i^\top, \\ \bm 1_{n \times n} &= n \bm u_1 \bm u_1^\top. \end{aligned}$$

ただし1つ目の固有ベクトルは $\bm u_1 = \bm 1_n / \sqrt{n}$ であり、またすべての固有ベクトル $\{ \bm u_i \}_{i=1}^n$ は正規直交系に選ばれるとする。線形結合もまた固有値分解された状態で簡単に書くことができて、

$$\begin{aligned} a \bm I_n + b \bm 1_{n \times n} &= (a + bn) \bm u_1 \bm u_1^\top + \sum_{i=2}^n a \bm u_i \bm u_i^\top. \end{aligned}$$

ゆえに、この線形結合された行列は $1$ 個の固有値 $a + bn$ と $n-1$ 個の固有値 $a$ を持つので、行列式は

$$\begin{aligned} \det \left( a \bm I_n + b \bm 1_{n \times n} \right) = a^{n-1} (a + bn) = a^n \left( 1 + n \frac{b}{a} \right). \end{aligned}$$

内部エネルギーは

$$\begin{aligned} e(\bm Q) &= \frac{\alpha}{2} \log \det \left( \bm I + \lambda \bm \Sigma \right) \\ &= \frac{\alpha}{2} \left( \log \left( 1 + n \frac{b}{a} \right) + n \log a \right). \end{aligned}$$

$a=1 + \lambda (Q - q)$ および $b = \lambda (\rho - 2m + q)$ を代入して

$$\begin{aligned} e(\bm Q) &= \frac{\alpha}{2} \left( \log \left( 1 + n \frac{\lambda (\rho - 2m + q)}{1 + \lambda(Q - q)} \right) + n \log \left( 1 + \lambda (Q - q) \right) \right). \end{aligned}$$

$\log( 1 + n(a/b) )$ の $n$ に関するMaclaurin展開

$$\begin{aligned} \log\left(1 + n\frac{a}{b}\right) &= n\frac{a}{b} + \mathcal O(n^2) \end{aligned}$$

を用いて近似すると

$$\begin{aligned} e(\bm Q) &\approx \frac{n \alpha}{2} \left( \frac{\lambda (\rho - 2m + q)}{1 + \lambda(Q - q)} + \log \left( 1 + \lambda (Q - q) \right) \right). \end{aligned}$$

エントロピーの計算

次いでエントロピーの部分の計算。エントロピーは次式で与えられるのであった。

$$\begin{aligned} s(\bm Q) &\coloneqq \frac{1}{N} \log \mathbb E_{\bm x^{(0)}} \int \prod_{a=1}^n d \bm x^{(a)}~ \exp\left( -\beta \sum_{a=1}^n \|\bm x^{(a)}\|_1 \right) P(\bm Q \mid \{ \bm x^{(a)} \}, \bm x^{(0)}). \end{aligned}$$

デルタ関数の積分表示

デルタ関数

$$\begin{aligned} P(\bm Q \mid \{ \bm x^{(a)} \}, \bm x^{(0)}) &= \prod_{a=0}^n \prod_{b=0}^n \delta \left( q^{ab} - \frac{1}{N} \bm x^{(a) \top} \bm x^{(b)} \right) \end{aligned}$$

を指数関数で表すために、Fourier変換の結果から導かれる以下の公式を使う。

$$\begin{aligned} \delta(x) &= \frac{|c|}{2\pi} \int d \tilde x \exp(i c \tilde x x). \end{aligned}$$

ただし $c$ は任意の実数である。これを使ってデルタ関数を積分表示する。レプリカ対称性が成り立つ場合、エントロピーの式中に出てくるデルタ関数は3通りある。

$$\begin{aligned} &(a = b)& &\delta \left( Q - \frac{1}{N} \bm x^{(a) \top} \bm x^{(a)} \right), \\ &(a \ne b)& &\delta \left( q - \frac{1}{N} \bm x^{(a) \top} \bm x^{(b)} \right), \\ &(a \ne 0, b = 0)& &\delta \left( m - \frac{1}{N} \bm x^{(a) \top} \bm x^{(0)} \right). \end{aligned}$$

後の便利のために、各デルタ関数の積分表示の引数について以下のように定数倍する。

1. $a=b$ の場合、$N/2$ 倍
2. $a \ne b$ の場合、$-N/2$ 倍
3. $a \ne 0, b = 0$ の場合、$-N$ 倍

以上を用いて

$$\begin{aligned} \delta \left( Q - \frac{1}{N} \bm x^{(a) \top} \bm x^{(a)} \right) &\propto \int d \tilde Q \exp\left( \frac{i \tilde Q}{2} \left( NQ - \bm x^{(a) \top} \bm x^{(a)} \right) \right), \\ \delta \left( q - \frac{1}{N} \bm x^{(a) \top} \bm x^{(b)} \right) &\propto \int d \tilde q \exp\left( - \frac{i \tilde q}{2} \left( Nq - \bm x^{(a) \top} \bm x^{(b)} \right) \right), \\ \delta \left( m - \frac{1}{N} \bm x^{(a) \top} \bm x^{(0)} \right) &\propto \int d \tilde m \exp\left( - i \tilde m \left( Nm - \bm x^{(a) \top} \bm x^{(0)} \right) \right). \end{aligned}$$

デルタ関数は以下のように変形できる。

$$\begin{aligned} & P(\bm Q \mid \{ \bm x^{(a)} \}, \bm x^{(0)}) \\ &= \prod_{a=0}^n \prod_{b=0}^n \delta \left( q^{ab} - \frac{1}{N} \bm x^{(a) \top} \bm x^{(b)} \right) \\ &\propto \int d \tilde Q~ \int d \tilde q~ \int d \tilde m \\&\qquad\qquad \exp\left( \sum_{a = 1}^n \frac{i \tilde Q}{2} \left( NQ - \bm x^{(a) \top} \bm x^{(a)} \right) - \sum_{a \ne b} \frac{i \tilde q}{2} \left( Nq - \bm x^{(a) \top} \bm x^{(b)} \right) - \sum_{a = 1}^n i \tilde m \left( Nm - \bm x^{(a) \top} \bm x^{(0)} \right) \right) \\ &= \int d \tilde Q~ \int d \tilde q~ \int d \tilde m \exp\left( \frac{i N n \tilde Q Q}{2} - \frac{i N n(n-1) \tilde q q}{2} - i N n \tilde m m \right) \\&\qquad\qquad \exp\left( - \sum_{a = 1}^n \frac{i\tilde Q}{2} \| \bm x^{(a)} \|_2^2 + \sum_{a \ne b} \frac{i \tilde q}{2} \bm x^{(a) \top} \bm x^{(b)} + \sum_{a = 1}^n i \tilde m \bm x^{(a) \top} \bm x^{(0)} \right) \\ &= \int d \tilde Q~ \int d \tilde q~ \int d \tilde m \exp\left( \frac{i N n \tilde Q Q}{2} - \frac{i N n(n-1) \tilde q q}{2} - i N n \tilde m m \right) \\&\qquad\qquad \exp\left( \sum_{a = 1}^n \left( - \frac{i (\tilde Q + \tilde q)}{2} \| \bm x^{(a)} \|_2^2 + i \tilde m \bm x^{(a) \top} \bm x^{(0)} \right) + \frac{i}{2} \left\| \sum_{a = 1}^n \tilde q \bm x^{(a)} \right\|_2^2 \right). \end{aligned}$$

最後の行では以下の関係を用いた。

$$\begin{aligned} \sum_{a \ne b} \bm x^{(a) \top} \bm x^{(b)} &= \left\| \sum_{a = 1}^n \bm x^{(a)} \right\|_2^2 - \sum_{a=1}^n \| \bm x^{(a)} \|_2^2. \end{aligned}$$

Stratonovich-Hubbard変換

多次元Gauss積分

$$\begin{aligned} \int d \bm x \frac{1}{\sqrt{ (2 \pi)^n}} \exp\left( - \frac{1}{2} \bm x^\top \bm x + \bm \mu^\top \bm x \right) = \exp\left( \frac{1}{2} \bm \mu^\top \bm \mu \right) \end{aligned}$$

を逆向きに使って得られる関係式

$$\begin{aligned} \exp \left( -\frac{1}{2} \| \bm x \|_2^2 \right) &= \int d \bm z \frac{1}{\sqrt{(2 \pi)^n}}\exp\left( -\frac{1}{2} \bm z^\top \bm z + \bm z^\top \bm x \right) \\ &= \mathbb E_{\bm z} ~ \exp\left( \bm z^\top \bm x \right) \end{aligned}$$

をStratonovich-Hubbard変換と呼ぶ。ただし $P(\bm z)$ は標準正規分布である。これを用いると、新たな積分変数 $\bm \zeta_i \in \mathbb R^N$ を導入することで

$$\begin{aligned} \exp \left( \frac{i}{2} \left\| \sum_{a = 1}^n \tilde q \bm x^{(a)} \right\|_2^2 \right) &= \mathbb E_{\bm z} ~ \exp\left( \bm z^\top \sum_{a = 1}^n \sqrt{i \tilde q} \bm x^{(a)} \right) \\ &= \mathbb E_{\bm z} ~ \prod_{a = 1}^n \exp\left( \sqrt{i \tilde q} \bm z^\top \bm x^{(a)} \right) \end{aligned}$$

が得られる。デルタ関数は

$$\begin{aligned} & P(\bm Q \mid \{ \bm x^{(a)} \}, \bm x^{(0)}) \\ &\propto \int d \tilde Q~ \int d \tilde q~ \int d \tilde m \exp\left( \frac{i N n \tilde Q Q}{2} - \frac{i N n(n-1) \tilde q q}{2} - i N n \tilde m m \right) \\&\qquad\qquad \mathbb E_{\bm z}~ \prod_{a = 1}^n \exp\left( - \frac{i (\tilde Q + \tilde q)}{2} \| \bm x^{(a)} \|_2^2 + \left( \sqrt{i \tilde q} \bm z + i \tilde m \bm x^{(0)} \right)^\top \bm x^{(a)} \right). \end{aligned}$$

エントロピーは

$$\begin{aligned} & s(\bm Q) \\ &= \frac{1}{N} \log \mathbb E_{\bm x^{(0)}} \int \prod_{a=1}^n d \bm x^{(a)}~ \exp\left( -\beta \sum_{a=1}^n \|\bm x^{(a)}\|_1 \right) P(\bm Q \mid \{ \bm x^{(a)} \}, \bm x^{(0)}) \\ &\approx \frac{1}{N} \log \int d \tilde Q~ \int d \tilde q~ \int d \tilde m \exp\left( \frac{i N n \tilde Q Q}{2} - \frac{i N n(n-1) \tilde q q}{2} - i N n \tilde m m \right) \\&\qquad\qquad \mathbb E_{\bm z}~ \mathbb E_{\bm x^{(0)}}~ \int \prod_{a=1}^n d \bm x^{(a)}~ \prod_{a = 1}^n \exp\left( - \frac{i (\tilde Q + \tilde q)}{2} \| \bm x^{(a)} \|_2^2 + \left( \sqrt{i \tilde q} \bm z + i \tilde m \bm x^{(0)} \right)^\top \bm x^{(a)} - \beta \| \bm x^{(a)} \|_1 \right). \end{aligned}$$

後半の $\bm x^{(a)} \in \mathbb R^N$ $(a=1,\dots,n)$, $\bm x^{(0)} \in \mathbb R^N$, $\bm z \in \mathbb R^N$ に関する積分は、$n$ 個のまったく同等な積分が $N$ 個あると解釈して、

$$\begin{aligned} & \mathbb E_{\bm z}~ \mathbb E_{x^{(0)}}~ \int \prod_{a=1}^n d \bm x^{(a)}~ \prod_{a = 1}^n \exp\left( - \frac{i (\tilde Q + \tilde q)}{2} \| \bm x^{(a)} \|_2^2 + \left( \sqrt{i \tilde q} \bm z + i \tilde m \bm x^{(0)} \right)^\top \bm x^{(a)} - \beta \| \bm x^{(a)} \|_1 \right) \\ &= \left( \mathbb E_{z}~ \mathbb E_{x^{(0)}}~ \left( \int d x~ \exp\left( - \frac{i (\tilde Q + \tilde q)}{2} \| x \|_2^2 + \left( \sqrt{i \tilde q} z + i \tilde m x^{(0)} \right) x - \beta | x | \right) \right)^n \right)^N \\ &= \exp\left( N \log \mathbb E_{z}~ \mathbb E_{x^{(0)}}~ \exp\left( n \log \phi(x^{(0)}, z; \tilde Q, \tilde q, \tilde m) \right) \right). \end{aligned}$$

ただし、$x (= x^{(a)}_j)$ に関する積分について

$$\begin{aligned} \phi(x^{(0)}, z; \tilde Q, \tilde q, \tilde m) \coloneqq \int dx \exp\left( - \frac{i (\tilde Q + \tilde q)}{2} x^2 + \left( \sqrt{i \tilde q} z + i \tilde m x^{(0)} \right) x - \beta |x| \right) \end{aligned}$$

とおいた。$\exp(n)$ と $\log(1 + n)$ のMaclaurin展開

$$\begin{aligned} &(1)& \exp(n) &= 1 + n + \mathcal O(n^2), \\ &(2)& \log(1 + n) &= n + \mathcal O(n^2) \end{aligned}$$

を用いれば、更に

$$\begin{aligned} & \exp\left( N \log \mathbb E_{z}~ \mathbb E_{x^{(0)}}~ \underline{ \exp\left( n \log \phi(x^{(0)}, z; \tilde Q, \tilde q, \tilde m) \right) }_{(1)} \right) \\ &\approx \exp\left( N \underline{ \log \mathbb E_{z}~ \mathbb E_{x^{(0)}}~ \left( 1 + n \log \phi(x^{(0)}, z; \tilde Q, \tilde q, \tilde m) \right) }_{(2)} \right) \\ &\approx \exp\left( N n \mathbb E_{z}~ \mathbb E_{x^{(0)}}~ \log \phi(x^{(0)}, z; \tilde Q, \tilde q, \tilde m) \right). \end{aligned}$$

エントロピーは

$$\begin{aligned} s(\bm Q) &\approx \frac{1}{N} \log \int d \tilde Q~ \int d \tilde q~ \int d \tilde m \\ &\qquad\qquad \exp\left( \frac{i N n \tilde Q Q}{2} - \frac{i N n(n-1) \tilde q q}{2} - i N n \tilde m m + N n \mathbb E_{x^{(0)}}~ \mathbb E_{z}~ \log \phi(x^{(0)}, z; \tilde Q, \tilde q, \tilde m) \right). \end{aligned}$$

Stratonovich-Hubbard変換の副作用である $\tilde Q, \tilde q, \tilde m$ に関する積分は、$N \to \infty$ 極限を取ることで鞍点法により評価できて、

$$\begin{aligned} s(\bm Q) &\approx \frac{1}{N} \log \operatorname*{extr}_{\tilde Q, \tilde q, \tilde m} \exp\left( \frac{N n \tilde Q Q}{2} - \frac{N n(n-1) \tilde q q}{2} - N n \tilde m m + N n \mathbb E_{x^{(0)}}~ \mathbb E_{z}~ \log \phi(x^{(0)}, z; \tilde Q, \tilde q, \tilde m) \right) \\ &= n \operatorname*{extr}_{\tilde Q, \tilde q, \tilde m} \left( \frac{\tilde Q Q}{2} - \frac{(n-1) \tilde q q}{2} - \tilde m m + \mathbb E_{x^{(0)}}~ \mathbb E_{z}~ \log \phi(x^{(0)}, z; \tilde Q, \tilde q, \tilde m) \right). \end{aligned}$$

$n \to 0$ とすることから、$n$ の2次以降の項を無視して、最終的に

$$\begin{aligned} s(\bm Q) &\approx n \operatorname*{extr}_{\tilde Q, \tilde q, \tilde m} \left( \frac{\tilde Q Q}{2} + \frac{\tilde q q}{2} - \tilde m m + \mathbb E_{x^{(0)}}~ \mathbb E_{z}~ \log \phi(x^{(0)}, z; \tilde Q, \tilde q, \tilde m) \right). \end{aligned}$$

変数変換

以上の結果をまとめると

$$\begin{aligned} e(\bm Q) &\approx \frac{n \alpha}{2} \left( \frac{\lambda (\rho - 2m + q)}{1 + \lambda(Q - q)} + \log \left( 1 + \lambda (Q - q) \right) \right), \\ s(\bm Q) &\approx n \operatorname*{extr}_{\tilde Q, \tilde q, \tilde m} \left( \frac{\tilde Q Q}{2} + \frac{\tilde q q}{2} - \tilde m m + \mathbb E_{x^{(0)}}~ \mathbb E_{z}~ \log \phi(x^{(0)}, z; \tilde Q, \tilde q, \tilde m) \right), \\ \phi(x^{(0)}, z; \tilde Q, \tilde q, \tilde m) &= \int dx \exp\left( - \frac{i (\tilde Q + \tilde q)}{2} x^2 + \left( \sqrt{i \tilde q} z + i \tilde m x^{(0)} \right) x - \beta |x| \right). \end{aligned}$$

続いて、$\lambda \to \infty$, $\beta \to \infty$ の極限を取っていく。$\beta \to \infty$ 極限を取るにあたって

$$\begin{aligned} &(1) & Q - q &\sim \mathcal O(1/\beta), && (\to 0) \\ &(2) & \tilde Q + \tilde q &\sim \mathcal O(\beta), && (\to \infty) \\ &(3) & \tilde m &\sim \mathcal O(\beta), && (\to \infty) \\ &(4) & \tilde q &\sim \mathcal O(\beta^2) && (\to \infty) \end{aligned}$$

という問題が生じるので、各変数を $\beta$ に関して $\mathcal O(1)$ に調整するために、以下の変数変換を行う。

$$\begin{aligned} &(1) & Q - q &= \frac{\chi}{\beta}, \\ &(2) & \tilde Q + \tilde q &= \beta \hat Q, \\ &(3) & \tilde m &= \beta \hat m, \\ &(4) & \tilde q &= \beta^2 \hat \chi. \end{aligned}$$

然らば、まず $x$ に関する積分が鞍点法により評価できて、

$$\begin{aligned} \phi(x^{(0)}, z; \tilde Q, \tilde q, \tilde m) &= \int dx \exp\left( - \frac{i (\tilde Q + \tilde q)}{2} x^2 + \left( \sqrt{i \tilde q} z + i \tilde m x^{(0)} \right) x - \beta |x| \right) \\ &= \int dx \exp\left( - \beta \frac{i \hat Q}{2} x^2 + \beta \left( \sqrt{i \hat \chi} z + i \hat m x^{(0)} \right) x - \beta |x| \right) \\ &\propto \max_{x} \exp \left( - \beta \left( \frac{\hat Q}{2} x^2 + \left( \sqrt{\hat \chi} z + \hat m x^{(0)} \right) x - |x| \right) \right) . \end{aligned}$$

内部エネルギーとエントロピーはそれぞれ

$$\begin{aligned} e(\bm Q) &\approx \frac{n \alpha}{2} \left( \frac{\lambda (\rho - 2m + q)}{1 + \lambda(Q - q)} + \log \left( 1 + \lambda (Q - q) \right) \right) \\ &= \frac{n \alpha}{2} \left( \frac{\lambda( \rho - 2 m + Q - \chi / \beta )}{1 + \lambda \chi / \beta} + \log \left( 1 + \frac{\lambda \chi}{\beta} \right) \right) \\ &= \frac{n \alpha \beta}{2} \left( \frac{\rho - 2 m + Q - \chi / \beta}{\beta / \lambda + \chi} + \frac{1}{\beta} \log \left( 1 + \frac{\lambda \chi}{\beta} \right) \right) \\ &\to \frac{n \alpha \beta}{2} \frac{\rho - 2 m + Q}{\chi} + \mathcal O(1), \end{aligned}$$

$$\begin{aligned} s(\bm Q) &\approx n \operatorname*{extr}_{\tilde Q, \tilde q, \tilde m} \left( \frac{\tilde Q Q}{2} + \frac{\tilde q q}{2} - \tilde m m + \log \mathbb E_{x^{(0)}}~ \mathbb E_{z}~ \phi(x^{(0)}, z; \tilde Q, \tilde q, \tilde m) \right) \\ &\approx n \operatorname*{extr}_{\hat \chi, \hat Q, \hat q, \hat m} \left( \frac{(\beta \hat Q - \beta^2 \hat \chi) Q}{2} + \frac{\beta^2 \hat \chi q}{2} - \beta \hat m m \right. \\&\qquad\qquad \left. + \log \mathbb E_{x^{(0)}}~ \mathbb E_{z}~ \max_{x} \exp \left( - \beta \left( \frac{\hat Q}{2} x^2 + \left( \sqrt{\hat \chi} z + \hat m x^{(0)} \right) x - |x| \right) \right) \right) \\ &= n \beta \operatorname*{extr}_{\tilde Q, \tilde q, \tilde m} \left( \frac{\hat Q Q}{2} - \frac{\beta \hat \chi (Q - q)}{2} - \hat m m + \mathbb E_{x^{(0)}}~ \mathbb E_{z}~ \max \left( - \frac{\hat Q}{2} x^2 + \left( \sqrt{\hat \chi} z + \hat m x^{(0)} \right) x - |x| \right) \right) \\ &= n \beta \operatorname*{extr}_{\tilde Q, \tilde q, \tilde m} \left( \frac{\hat Q Q}{2} - \frac{\hat \chi \chi}{2} - \hat m m - \mathbb E_{x^{(0)}}~ \mathbb E_{z}~ \min \left( \frac{\hat Q}{2} x^2 - \left( \sqrt{\hat \chi} z + \hat m x^{(0)} \right) x + |x| \right) \right). \end{aligned}$$

再パラメータ化トリック

以下のような新たな確率変数を導入する。

$$\begin{aligned} \eta \coloneqq \sqrt{\hat \chi} z + \hat m x^{(0)}. \end{aligned}$$

すると

$$\begin{aligned} P(x^{(0)}) &= (1 - \rho) \delta(x^{(0)}) + \rho \mathcal N(0, 1), \\ P(z) &= \mathcal N(0, 1) \end{aligned}$$

であることから、それらの線形結合である $\eta$ の確率分布は以下のようになる。

$$\begin{aligned} P(\eta) &= P(\sqrt{\hat \chi} z + \hat m x^{(0)}) \\ &= (1 - \rho) \mathcal N(0, \hat \chi) + \rho \mathcal N(0, \hat \chi + \hat m^2). \end{aligned}$$

よって標準正規分布に従う新たな変数 $\hat z$ を導入して、

$$\begin{aligned} \mathbb E_\eta[f(\eta)] = (1 - \rho) \mathbb E[f(\sqrt{\hat \chi} \hat z)] + \rho \mathbb E[f(\sqrt{\hat \chi + \hat m^2} \hat z)] \end{aligned}$$

という書き換えが可能である。以上をエントロピーに現れる期待値操作に適用すると、

$$\begin{aligned} & \mathbb E_{x^{(0)}}~ \mathbb E_{z}~ \min \left( \frac{\hat Q}{2} x^2 - \left( \sqrt{\hat \chi} z + \hat m x^{(0)} \right) x + |x| \right) \\ &= \mathbb E_{\eta} \min_x \left( \frac{\hat Q}{2} x^2 - \eta x + |x| \right) \\ &= (1 - \rho) \mathbb E_{\hat z} \min_x \left( \frac{\hat Q}{2} x^2 + \sqrt{\hat \chi} \hat z x - |x| \right) + \rho \mathbb E_{\hat z} \min_x \left( \frac{\hat Q}{2} x^2 + \sqrt{\hat \chi + \hat m^2} \hat z x - |x| \right) \\ &= (1 - \rho) \mathbb E_{\hat z} \Phi(\hat z ; \hat Q, \hat \chi, 0) + \rho \mathbb E_{\hat z} \Phi(\hat z ; \hat Q, \hat \chi, \hat m). \end{aligned}$$

ただし以下のように定義した。

$$\begin{aligned} \Phi(\hat z ; \hat Q, \hat \chi, \hat m) &\coloneqq \min_{x} \left( \frac{\hat Q}{2} x^2 + \sqrt{\hat \chi + \hat m^2} \hat z x - |x| \right). \end{aligned}$$

この関数の最小値は場合分けにより計算できる。

最小値の計算

ここで簡単のために $f(x) = ax^2/2 + bx - c|x|$ $(a \gt 0)$ を考察する。

(1) $x = 0$ のとき は明らかに $0$。

(2) $x \gt 0$ のとき は、平方完成すると

$$\begin{aligned} f(x) &= \frac{a}{2} x^2 + (b - c) x \\ &= \frac{a}{2} \left( x + \frac{b - c}{a} \right)^2 - \frac{(b - c)^2}{2a} \end{aligned}$$

を得るので、グラフの頂点は $x=- (b - c) / (2a)$ にある。ところが $x \gt 0$ でなくてはならないので、$a \gt 0$ かつ $b \gt c$ のとき最小値は $- (b - c)^2 / (2a)$ である。

(3) $x \lt 0$ のとき は、平方完成すると

$$\begin{aligned} f(x) &= \frac{a}{2} x^2 + (b + c) x \\ &= \frac{a}{2} \left( x - \frac{b + c}{a} \right)^2 - \frac{(b + c)^2}{2a} \end{aligned}$$

を得るので、グラフの頂点は $x = (b + c) / (2a)$ にある。ところが $x \lt 0$ でなくてはならないので、$a \gt 0$ かつ $b \lt -c$ のとき最小値は $- (b + c)^2 / (2a)$ である。

以上を取りまとめると

$$\begin{aligned} \min_x f(x) &= \left\{ \begin{aligned} & - \frac{(b + c)^2}{2a}, && b \gt c, \\ & 0, && |b| \lt c, \\ & - \frac{(b - c)^2}{2a} && b \lt c \end{aligned} \right. \\ &= - \frac{1}{2a} \left( |b| - c \right)^2 \Theta \left( |b| - c \right) \end{aligned}$$

ただし $\Theta(x)$ はHeavisideのステップ関数であり、

$$\begin{aligned} \Theta(x) &\coloneqq \mathbb I \left[ x \ge 0 \right] \\ &= \left\{ \begin{aligned} & 1, && x \ge 0, \\ & 0, && x \lt 0 \end{aligned} \right. \end{aligned}$$

で定義される。

以上の結果に $a =\hat Q$, $b = \sqrt{\hat \chi + \hat m^2} \hat z$, $c = 1$ を代入すれば、

$$\begin{aligned} \Phi(\hat z ; \hat Q, \hat \chi, \hat m) &= - \frac{1}{2 \hat Q} \left( \left| \sqrt{\hat \chi + \hat m^2} \hat z \right| - 1 \right)^2 \Theta \left( \left| \sqrt{\hat \chi + \hat m^2} \hat z \right| - 1 \right) \\ &= - \frac{\hat \chi + \hat m^2}{2 \hat Q} \left( |\hat z| - \frac{1}{\sqrt{\hat \chi + \hat m^2}} \right)^2 \Theta \left( |\hat z| - \frac{1}{\sqrt{\hat \chi + \hat m^2}} \right). \end{aligned}$$

期待値計算

つぎに $\Phi(\hat z ; \hat Q, \hat \chi, \hat m)$ の $z$ に関する期待値を計算する。

$$\begin{aligned} \mathbb E_{\hat z}~ \Phi(\hat z ; \hat Q, \hat \chi, \hat m) &= - \frac{\hat \chi + \hat m^2}{2 \hat Q} \mathbb E_{\hat z}~ \left( |\hat z| - \frac{1}{\sqrt{\hat \chi + \hat m^2}} \right)^2 \Theta \left( |\hat z| - \frac{1}{\sqrt{\hat \chi + \hat m^2}} \right) \end{aligned}$$

$x = 1 / \sqrt{\hat \chi + \hat m^2}$ とおくと、期待値部分は

$$\begin{aligned} & \frac{1}{2} \mathbb E_{\hat z}~ \left( |\hat z| - \frac{1}{\sqrt{\hat \chi + \hat m^2}} \right)^2 \Theta \left( |\hat z| - \frac{1}{\sqrt{\hat \chi + \hat m^2}} \right) \\ &= \frac{1}{2} \frac{1}{\sqrt{2 \pi}} \int_{-\infty}^\infty d \hat z~ \exp \left( - \frac{\hat z^2}{2} \right) \left( \hat z^2 - 2 |\hat z| x + x^2 \right) \Theta \left( |\hat z| - x \right) \\ &= \frac{1}{\sqrt{2 \pi}} \int_{x}^\infty d\hat z~ \exp \left( - \frac{\hat z^2}{2} \right) \left( \hat z^2 - 2 \hat z x + x^2 \right). \end{aligned}$$

第1項は

$$\begin{aligned} \int_{x}^\infty d\hat z~ \hat z^2 \exp \left( - \frac{\hat z^2}{2} \right) &= \left[ - \hat z \exp \left( - \frac{\hat z^2}{2} \right) \right]_{x}^\infty + \int_{x}^\infty d\hat z~ \exp \left( - \frac{\hat z^2}{2} \right) \\ &= x \exp \left( - \frac{x^2}{2} \right) + \int_{x}^\infty d\hat z~ \exp \left( - \frac{\hat z^2}{2} \right). \end{aligned}$$

第2項は

$$\begin{aligned} -2 x \int_{x}^\infty d\hat z~ \hat z \exp \left( - \frac{\hat z^2}{2} \right) = -2 x \exp \left( - \frac{x^2}{2} \right). \end{aligned}$$

第3項は

$$\begin{aligned} x^2 \int_{x}^\infty d\hat z~ \exp \left( - \frac{\hat z^2}{2} \right). \end{aligned}$$

以上3項を足し合わせて

$$\begin{aligned} & \frac{1}{2} \mathbb E_{\hat z}~ \left( |\hat z| - \frac{1}{\sqrt{\hat \chi + \hat m^2}} \right)^2 \Theta \left( |\hat z| - \frac{1}{\sqrt{\hat \chi + \hat m^2}} \right) \\ &= (x^2 + 1) \frac{1}{\sqrt{2 \pi}} \int_{x}^\infty d\hat z~ \exp \left( - \frac{\hat z^2}{2} \right) - x \frac{1}{\sqrt{2 \pi}} \exp \left( - \frac{x^2}{2} \right). \end{aligned}$$

これを $x^2 = 1 / (\hat \chi + \hat m^2)$ で割ったものを $G(x)$ と書く、すなわち

$$\begin{aligned} G(x) \coloneqq \left( 1 + \frac{1}{x^2} \right) \frac{1}{\sqrt{2 \pi}} \int_{x}^\infty d\hat z~ \exp \left( - \frac{\hat z^2}{2} \right) - \frac{1}{x} \frac{1}{\sqrt{2 \pi}} \exp \left( - \frac{x^2}{2} \right) \end{aligned}$$

と定めれば

$$\begin{aligned} \mathbb E_{\hat z}~ \Phi(\hat z ; \hat Q, \hat \chi, \hat m) &= - \frac{1}{\hat Q} G\left( \frac{1}{\sqrt{\hat \chi + \hat m^2}} \right). \end{aligned}$$

こうして最小値全体の期待値は

$$\begin{aligned} & \mathbb E_{x^{(0)}}~ \mathbb E_{z}~ \min \left( \frac{\hat Q}{2} x^2 - \left( \sqrt{\hat \chi} z + \hat m x^{(0)} \right) x + |x| \right) \\ &= (1 - \rho) \mathbb E_{\hat z}~ \Phi(\hat z ; \hat Q, \hat \chi, 0) + \rho \mathbb E_{\hat z}~ \Phi(\hat z ; \hat Q, \hat \chi, \hat m) \\ &= - \frac{1 - \rho}{\hat Q} G\left( \frac{1}{\sqrt{\hat \chi}} \right) - \frac{\rho}{\hat Q} G\left( \frac{1}{\sqrt{\hat \chi + \hat m^2}} \right). \end{aligned}$$

エントロピーは

$$\begin{aligned} s(\bm Q) &\approx n \beta \operatorname*{extr}_{\tilde Q, \tilde q, \tilde m} \left( \frac{\hat Q Q}{2} - \frac{\hat \chi \chi}{2} - \hat m m - \mathbb E_{x^{(0)}}~ \mathbb E_{z}~ \min \left( \frac{\hat Q}{2} x^2 - \left( \sqrt{\hat \chi} z + \hat m x^{(0)} \right) x + |x| \right) \right) \\ &= n \beta \operatorname*{extr}_{\hat Q, \hat \chi, \hat m} \left( \frac{\hat Q Q}{2} - \frac{\hat \chi \chi}{2} - \hat m m + \frac{1 - \rho}{\hat Q} G\left( \frac{1}{\sqrt{\hat \chi}} \right) + \frac{\rho}{\hat Q} G\left( \frac{1}{\sqrt{\hat \chi + \hat m^2}} \right) \right). \end{aligned}$$

鞍点法の適用

最終的な分配関数の式は $N \to \infty$ より鞍点法により評価できて、

$$\begin{aligned} & \mathbb E_{\bm A, \bm y}[Z^n(\bm A, \bm y)] \\ &= \int d \bm Q \exp\left( - N e(\bm Q) + N s(\bm Q) \right) \\ &\propto \operatorname*{extr}_{Q, q, m} \exp\left( - N e(\bm Q) + N s(\bm Q) \right) \\ &\propto \operatorname*{extr}_{Q, \chi, m, \hat Q, \hat \chi, \hat m} \exp\left\{ N n \beta \left( - \frac{\alpha}{2} \frac{\rho - 2 m + Q}{\chi} + \frac{\hat Q Q}{2} - \frac{\hat \chi \chi}{2} - \hat m m \right.\right. \\ & \qquad \qquad \qquad \qquad \left. \left. + \frac{1 - \rho}{\hat Q} G\left( \frac{1}{\sqrt{\hat \chi}} \right) + \frac{\rho}{\hat Q} G\left( \frac{1}{\sqrt{\hat \chi + \hat m^2}} \right) \right) \right\} \end{aligned}$$

であり、自由エネルギーは

$$\begin{aligned} f &\approx \operatorname*{extr}_{Q, \chi, m, \hat Q, \hat \chi, \hat m} \left\{ - \frac{\alpha}{2} \frac{\rho - 2 m + Q}{\chi} + \frac{\hat Q Q}{2} - \frac{\hat \chi \chi}{2} - \hat m m \right. \\ & \qquad \qquad \qquad \qquad \left. + \frac{1 - \rho}{\hat Q} G\left( \frac{1}{\sqrt{\hat \chi}} \right) + \frac{\rho}{\hat Q} G\left( \frac{1}{\sqrt{\hat \chi + \hat m^2}} \right) \right\}, \end{aligned}$$

$$\begin{aligned} G(x) &= \left( 1 + \frac{1}{x^2} \right) H(x) - \frac{1}{x} \frac{1}{\sqrt{2 \pi}} \exp \left( - \frac{x^2}{2} \right), \\ H(x) &= \frac{1}{\sqrt{2 \pi}} \int_{x}^\infty d\hat z~ \exp \left( - \frac{\hat z^2}{2} \right). \end{aligned}$$

鞍点方程式

自由エネルギーの括弧の中を各秩序変数に関して偏微分して零となる方程式 (鞍点方程式) を導出する。$Q$, $\chi$, $m$, $\hat Q$, $\hat \chi$, $\hat m$ についてそれぞれ偏微分して零にすると

$$\begin{aligned} 0 &= - \frac{\alpha}{2 \chi} + \frac{\hat Q}{2} , \\ 0 &= \frac{\alpha (\rho - 2m + Q)}{2 \chi^2} - \frac{\hat \chi}{2} , \\ 0 &= \frac{\alpha}{\chi} - \hat m , \\ 0 &= \frac{Q}{2} - \frac{1 - \rho}{\hat Q^2} G\left( \frac{1}{\sqrt{\hat \chi}} \right) - \frac{\rho}{\hat Q^2} G\left( \frac{1}{\sqrt{\hat \chi + \hat m^2}} \right) , \\ 0 &= - \frac{\chi}{2} - \frac{1 - \rho}{\hat Q} \frac{1}{2 \sqrt{\hat \chi^3}} \left. \frac{\partial G (x)}{\partial x} \right|_{x = 1 / \sqrt{\hat \chi}} - \frac{\rho}{\hat Q} \frac{1}{2 \sqrt{(\hat \chi + \hat m^2)^3}} \left. \frac{\partial G (x)}{\partial x} \right|_{x = 1 / \sqrt{\hat \chi + \hat m^2}} , \\ 0 &= - m - \frac{\rho}{\hat Q} \frac{m}{\sqrt{(\hat \chi + \hat m^2)^3}} \left. \frac{\partial G (x)}{\partial x} \right|_{x = 1 / \sqrt{\hat \chi + \hat m^2}} . \end{aligned}$$

ここで $G$ の偏微分は

$$\begin{aligned} \frac{\partial G(x)}{\partial x} &= - \frac{2}{x^3} H(x) \end{aligned}$$

であるので、鞍点方程式はさらに

$$\begin{aligned} 0 &= - \frac{\alpha}{2 \chi} + \frac{\hat Q}{2} , \\ 0 &= \frac{\alpha (\rho - 2m + Q)}{2 \chi^2} - \frac{\hat \chi}{2} , \\ 0 &= \frac{\alpha}{\chi} - \hat m , \\ 0 &= \frac{Q}{2} - \frac{1 - \rho}{\hat Q^2} G\left( \frac{1}{\sqrt{\hat \chi}} \right) - \frac{\rho}{\hat Q^2} G\left( \frac{1}{\sqrt{\hat \chi + \hat m^2}} \right) , \\ 0 &= - \frac{\chi}{2} + \frac{1 - \rho}{\hat Q} H \left( \frac{1}{\sqrt{\hat \chi}} \right) + \frac{\rho}{\hat Q} H \left( \frac{1}{\sqrt{\hat \chi + \hat m^2}} \right) , \\ 0 &= - m + \frac{2 \rho \hat m}{\hat Q} H \left( \frac{1}{\sqrt{\hat \chi + \hat m^2}} \right) . \end{aligned}$$

都合の良いことに、各変数について1次の項が各式に現れる。それを左辺に持ってきて整理すれば以下の自己無撞着方程式に帰着する。

$$\begin{aligned} \hat Q &= \frac{\alpha}{\chi}, \\ \hat \chi &= \frac{\alpha (\rho - 2m + Q)}{\chi^2}, \\ \hat m &= \frac{\alpha}{\chi}, \\ Q &= \frac{2(1 - \rho)}{\hat Q} G \left( \frac{1}{\sqrt{\hat \chi}} \right) + \frac{2\rho}{\hat Q} G \left( \frac{1}{\sqrt{\hat \chi + \hat m^2}} \right), \\ \chi &= \frac{2(1 - \rho)}{\hat Q} H \left( \frac{1}{\sqrt{\hat \chi}} \right) + \frac{2\rho}{\hat Q} H \left( \frac{1}{\sqrt{\hat \chi + \hat m^2}} \right), \\ m &= \frac{2 \rho \hat m}{\hat Q} H \left( \frac{1}{\sqrt{\hat \chi + \hat m^2}} \right). \end{aligned}$$

ただし

$$\begin{aligned} G(x) &= \left( 1 + \frac{1}{x^2} \right) H(x) - \frac{1}{x} \frac{1}{\sqrt{2 \pi}} \exp \left( - \frac{x^2}{2} \right), \\ H(x) &= \frac{1}{\sqrt{2 \pi}} \int_{x}^\infty d\hat z~ \exp \left( - \frac{\hat z^2}{2} \right). \end{aligned}$$

自己無撞着方程式が収束したときの解 $Q$, $\chi, m$, $\hat Q$, $\hat \chi$, $\hat m$ が熱力学極限 $N \to \infty$ における自由エネルギー $f$ の極値を与える。そのときの平均二乗誤差 (MSE) の値

$$\begin{aligned} \mathrm{MSE} &= \lim_{N \to \infty} \lim_{\beta \to \infty} \mathbb E_{\bm A}~ \mathbb E_{\bm x^{(0)}}~ \mathbb E_{\bm x \mid \bm A, \bm x^{(0)}}~ \frac{1}{N} \sum_{i=1}^N (x^{(0)}_i - x_i)^2 \\ &= \rho - 2m + Q \end{aligned}$$

の様子を調べることで、L1ノルム最小化のレプリカ対称性条件下の性能限界を見積もることができる。

鞍点の探索

簡易的にPythonで鞍点を探索するコードを書いてシミュレーションする。計算自体はそれなりに高速なのだが、わずかにでも係数がミスっているととんでもない結果を与えてくるので注意する。

また $H(x)$ は補誤差関数 $\mathrm{erfc}(x)$ を用いて

$$\begin{aligned} H(x) &= \frac{1}{2} \operatorname{erfc}\left( \frac{x}{\sqrt{2}} \right) \end{aligned}$$

と表せることにも注意しておく。

import numpy as np
from scipy.special import erfc
import matplotlib.pyplot as plt
from matplotlib import colors as mcolors


def calc_H(x):
    if x < -7:
        return 1
    if x > 7:
        return 0
    return erfc(x / np.sqrt(2)) / 2


def calc_G(x):
    result = (1 + 1 / x**2) * calc_H(x)
    result -= 1 / x / np.sqrt(2 * np.pi) * np.exp(-(x**2) / 2)
    return result


def calc_hat_Q(chi, alpha):
    hat_Q = alpha / chi
    return hat_Q


def calc_hat_chi(Q, m, chi, alpha, rho):
    hat_chi = alpha * (rho - 2 * m + Q) / chi**2
    if hat_chi < 0:
        return 0
    return hat_chi


def calc_hat_m(chi, alpha):
    hat_m = alpha / chi
    return hat_m


def calc_Q(hat_Q, hat_chi, hat_m, rho):
    result = 2 * (1 - rho) / hat_Q**2 * calc_G(1 / np.sqrt(hat_chi))
    result += 2 * rho / hat_Q**2 * calc_G(1 / np.sqrt(hat_chi + hat_m**2))
    return result


def calc_chi(hat_Q, hat_chi, hat_m, rho):
    result = 2 * (1 - rho) / hat_Q * calc_H(1 / np.sqrt(hat_chi))
    result += 2 * rho / hat_Q * calc_H(1 / np.sqrt(hat_chi + hat_m**2))
    return result


def calc_m(hat_Q, hat_chi, hat_m, rho):
    return 2 * rho * hat_m / hat_Q * calc_H(1 / np.sqrt(hat_chi + hat_m**2))


def saddle_point_loop(
    Q,
    chi,
    m,
    hat_Q,
    hat_chi,
    hat_m,
    alpha,
    rho,
    max_iter,
    damping=0.1,
    tol=1e-10,
):
    MSE_hist = []
    is_converged = False
    hist = []
    for i in range(max_iter):
        hat_Q = (1 - damping) * hat_Q + damping * calc_hat_Q(chi, alpha)
        hat_chi = (1 - damping) * hat_chi + damping * calc_hat_chi(
            Q, m, chi, alpha, rho
        )
        hat_m = (1 - damping) * hat_m + damping * calc_hat_m(chi, alpha)
        Q = (1 - damping) * Q + damping * calc_Q(hat_Q, hat_chi, hat_m, rho)
        chi = (1 - damping) * chi + damping * calc_chi(hat_Q, hat_chi, hat_m, rho)
        m = (1 - damping) * m + damping * calc_m(hat_Q, hat_chi, hat_m, rho)

        hist.append([hat_Q, hat_chi, hat_m, Q, chi, m])

        MSE = rho - 2 * m + Q
        MSE_hist.append(MSE)

        if MSE < tol:
            is_converged = True
            break
    return hist, MSE_hist, is_converged


def compute_converged_matrix(I, J):
    converged_matrix = np.zeros((I, J))
    mse_matrix = np.zeros((I, J))
    for i, alpha in enumerate(np.linspace(0, 1 + 1 / I, I)):
        for j, rho in enumerate(np.linspace(0, 1 + 1 / J, J)):
            print(i, j)
            Q = 1
            chi = 1
            m = 0
            hat_Q = 1
            hat_chi = 1
            hat_m = 0

            hist, mse_hist, is_converged = saddle_point_loop(
                Q, chi, m, hat_Q, hat_chi, hat_m, alpha, rho, max_iter=2000
            )

            mse_matrix[i, j] = mse_hist[-1]
            converged_matrix[i, j] = is_converged

    return mse_matrix, converged_matrix


mse_matrix, converged_matrix = compute_converged_matrix(50, 50)


fig, ax = plt.subplots(1, 1, figsize=(3, 3), dpi=200)
im = ax.imshow(mse_matrix[::-1], cmap="viridis", extent=[0, 1, 0, 1], norm=mcolors.LogNorm())
cbar = fig.colorbar(im, orientation="vertical")
axpos = ax.get_position()
caxpos = cbar.ax.get_position()
cbar.ax.set_position([caxpos.x0, axpos.y0, caxpos.width, axpos.height])

ax.set_xlabel(r"$\rho = \dfrac{K}{N}$")
ax.set_ylabel(r"$\alpha = \dfrac{M}{N}$")
plt.show()

ここで3つほど注意点を挙げておく。

1. $H(x)$ 関数は $x \to \pm \infty$ で指数関数的に収束するので、十分に大きい $A$ に対して $x \gt A$ で $H(x) \approx 0$、$x \lt -A$ で $H(x) \approx 1$ として計算してよい。
2. 計算を安定化させるために、damping因子 $0 \lt \mathrm{damping} \lt 1$ を導入して、前回の値と計算された新たな値の中間の値へ遷移している。damping因子の値が小さいほど前回の値を重視するため、収束が安定する。
3. MSEが零に収束する場合、MSEは急激に減衰していくので計算が不安定になる。計算の不安定さを解消する方法は色々あるが、今回はMSEが一定の値 tol 以下になったら再構成に成功したとみなしてループを抜けるようにしている。この方法は簡単だが、計算結果から MSE < tol の箇所の情報が欠落してしまう。

このコードは1-2分程度で以下の図を生成する。

図の明るい箇所はMSEが大きい箇所で、再構成に失敗する箇所を、暗い箇所はMSEが小さい箇所で、再構成に成功する箇所を示している。

図の右上は変化がやや滑らかであるのに対し、左下は変化が急峻である。$\alpha \approx 0$ 付近では1次転移が生じていることを示唆している。

プロットの下限を $10^{-3}$ とする場合は以下のようになる。

この図は、よく知られたL1ノルム最小化による再構成の典型的な性能限界とよく整合している。左上が成功層、右下が失敗層である。

参考文献

1. Y. Kabashima, T. Wadayama, and T. Tanaka, A typical reconstruction limit for compressed sensing based on Lp-norm minimization, J. Stat. Mech.: Theory and Experiment, L09003 (2009).
2. 黒木 学, 清水 昌平, 湊 真一, 石畠 正和, 樺島 祥介, 田中 和之, 木村 陽一, 玉田 嘉紀, 確率的グラフィカルモデル, 共立出版, p. 280 (2016).
3. 片岡 駿, 大関 真之, 安田 宗樹, 田中 和之, 画像処理の統計モデリング: 確率的グラフィカルモデルとスパースモデリングからのアプローチ, 共立出版, p. 174 (2018).

脚注

1. レプリカ数 $n$ を逆温度の比 $\tilde \beta / \beta$ で定義する方法は「パーシャルアニーリング」の文脈で使う。クエンチ系を考察する場合はむしろ、$Z^n \approx 1 + n \log Z$ というテイラー展開より $\log Z = \lim_{n \to +0} (Z^n - 1) / n = \lim_{\nu \to +0} \partial_n Z^n |_{n = \nu}$ という式変形の結果として分配関数の $n$ 乗が出てくるという解釈が一般的だと思う。