ADMMについてまとめてみる


こんにちは、sakura__mochiです。 授業でADMMによる画像のノイズ除去について学んだので、まとめたいと思います。

## 1. どんなもの？
ADMM(Alternating Direction Method of Multipliers)とは、大域的に難しい凸最適化問題を「分割して交互に解く」手法です。特に以下のような制約付き凸問題を想定します：

```math
\begin{aligned}
&\min_{x, z} \quad f(x) + g(z) \\
&\text{subject to} \quad A x + B z = c.
\end{aligned}\tag{1.1}
```
$f, g$ は凸関数であり、$A, B$ は行列、$c$ はベクトルです。ADMMはラグランジュ関数に二乗ペナルティ項（ペナルティ係数ρ）を加えた拡張ラグランジュ (Augmented Lagrangian) に基づく手法です。拡張ラグランジアン（augmented Lagrangian）を
```math
Lρ(x,z,u) = f(x) + g(z) + u^T (Ax + Bz − c)+(ρ/2)\|Ax + Bz − c\|_2^2.\tag{1.2}
```
として構成するとき、ADMM は次の反復から成ります。
```math
x^{k+1} := \arg\min_{x}\,
L_\rho(x, z^k, u^k),
\tag{1.3}
```
```math
z^{k+1} := \arg\min_{z}\, L_\rho(x^{k+1}, z, u^k),
\tag{1.4}
```
```math
u^{k+1} := u^k + \rho \bigl(Ax^{k+1} + Bz^{k+1} - c\bigr),
\tag{1.5}
```
このようにADMMでは$x$と$z$を交互に最適化し（$x$更新では最新の$z$を用い、その後$z$更新で最新の$x$を用いる）、最後に双対変数$u$を更新します。


## 2. 先行研究と比べてどこがすごいの？

ADMMは、FISTAやSplit Bergman法に比べ分散処理や柔軟性で群を抜いています。収束速度は並程度でエルゴード平均$O(1/k)$と言われています。

基盤となる理論提案から応用比較研究まで数多くの先行研究がありますが、総じてADMMは**古くから知られる手法が現代の大規模問題に再評価され普及した**という歴史を持ち、FISTAやChambolle-Pock法は**新規に提案された理論的改良**として位置づけられます。

以下に関連研究を載せます：


- [Boyd et al., 2011 – Distributed Optimization and Statistical Learning via the ADMM](https://web.stanford.edu/~boyd/papers/pdf/admm_distr_stats.pdf)
- [Glowinski & Marroco, 1975](https://www.ijcai.org/proceedings/2019/0629.pdf) / [Gabay & Mercier, 1976](https://www.sciencedirect.com/science/article/pii/0898122176900031?via%3Dihub)
- [Goldstein & Osher, 2009 – The Split Bregman Method for L1 Regularized Problems](https://ww3.math.ucla.edu/camreport/cam08-29.pdf)
- [Beck & Teboulle, 2009 – A Fast Iterative Shrinkage-Thresholding Algorithm](https://ieeexplore.ieee.org/document/4959678)
- [Chambolle & Pock, 2011 – A First-Order Primal-Dual Algorithm for Convex Problems](https://link.springer.com/article/10.1007/s11425-015-5079-0)

## 3. 技術や手法の"キモ"はどこにある？
ADMMの"キモ"は、大域的な問題を解くために局所的な部分問題へ分解し、それらを協調させる(漸次的に解く)ことです。

今日のデータ分析は様々な分野で応用されていますが、扱われるデータには共通の特徴が見られます：
1. データセットがとてつもなく大きい
2. 高次元である
3. データの収集や保管が分散的であること

その結果、データの複雑性を捉えつつ、並列or分散的にデータを処理できる拡張性が求められます。これら多くの問題は凸最適化問題として提起することができ、ADMMはこのような分散的に凸最適化問題を解くことに適しています。


## 4. どうやって有効だと検証した？

LASSO, SVM,sparse logistic regression, basis pursuit, covariance selectionなどの問題で応用事例がありますが、今回は画像のノイズ除去問題を通して有効性を検証したいと思います。
問題は次のように定義されます：
> 「ノイズが付加された画像しか観測できない状況で、元のきれいな画像を推定する」

たとえば、下記の例では左がノイズなしの元画像、右がガウシアンノイズ（正規分布に従うノイズ）を付加した画像です。目標は「右の観測画像から左の元画像をなるべく忠実に復元すること」になります。

| ノイズなし | ノイズあり
| -- | -- |
|![スクリーンショット 2025-04-02 20.22.03.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3933005/057fba73-5bcb-4206-9a54-d29a05d74817.png)|![スクリーンショット 2025-04-02 20.23.20.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3933005/a332ee04-6f5e-4b9e-97a7-1068de60ef2b.png)|

<br>
この問題を凸最適化問題の枠組みで解く方法を考えていきます。

一般に“自然な画像”では、隣り合うピクセル間の画素値は急激には変化しにくく、画像のほとんどの領域は比較的滑らかであるという特徴があります。一方で、物体の輪郭やエッジの部分では値の変化が大きくなる、という局所的なパターンが存在します。

この「隣接ピクセル間の滑らかさ」を定量化するために、画像に差分フィルターをかけてみます。以下は元画像・ノイズ付き画像に対して垂直方向の差分フィルターを掛けた結果です。

| 差分フィルタリング(ノイズなし) | 差分フィルタリング(ノイズあり) |
| -- | -- |
|![スクリーンショット 2025-04-02 21.13.11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3933005/2b16356f-f6b5-478f-bdf7-6a0e1aaf2671.png)|![スクリーンショット 2025-04-02 21.13.36.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3933005/959c385d-6ca6-4bc3-a571-29f84c42a91e.png)|

画像から分かるように、ノイズが乗ることで高周波成分（細かい乱れ）が顕著になり，差分画像にも大きなばらつきが生じている様子が分かります。これを定式化すると、垂直フィルター$D_v$、水平フィルター$D_h$、求める画像$x$とした時、

```math
\|D_vx\|_2^2 + \|D_hx\|_2^2\tag{4.1}
```

が最小値となるxを求めれば良いことになります。ただし、「とにかく差分が小さい画像」を求めると、物体の輪郭やテクスチャが犠牲になる可能性があります。そのため、元の画像にもある程度近い必要があります。これを式にすると、

```math
\frac{1}{2}\|x-\tilde{x}\|_2^2 +  \lambda\|D_vx\|_2^2 + \lambda\|D_hx\|_2^2\tag{4.2}
```
ここで、$D = \begin{bmatrix}D_v \\\D_h \end{bmatrix}$と置くと、(4.2)は
```math
\frac{1}{2}\|x-\tilde{x}\|_2^2 +  \lambda\|Dx\|_2^2 \tag{4.3}
```
となる。こちらを式(1.1)に当てはめると以下のようになりまs。
```math
f(x)=\frac{1}{2}\|x-\tilde{x}\|_2^2 
```
```math
g(x)=\lambda\|Dx\|_2^2
```
```math
z = Dx, A = D, B = -I, c = 0  \tag{Iは単位行列 }
```
したがって、式(1.3)、(1.4)、(1.5)は、
```math
x^{(k+1)}= arg\min_{x}(\frac{1}{2}\|x-\tilde{x}\|^2_2 + \frac{\rho}{2}\|Dx-z^{(k)}+u^{(k)}\|^2_2 )
```
```math
= (I + \rho D^TD)^{-1}(\tilde{x}+\rho D^T((z)^{(k)}-u^{(k)}))\tag{4.4}
```

```math
z^{(k+1)}= arg\min_{z}(\lambda\ \|z\|^2_2 + \frac{\rho}{2}\|Dx^{(k+1)}-z+u^{(k)}\|^2_2)
```
```math
= \frac{\rho}{2\lambda+\rho}(Dx^{(k+1)}+u^{(k)})\tag{4.5}

```
```math
u^{(k+1)}+u^{(k)}+(x^{(k+1)}+z^{(k+1)})\tag{4.6}
```
となります。さて、最後にmatlabのコードを書いてこの問題を解いてみます。

```matlab:ADMM
clear; clc; close all;

img = imread('Pepper.bmp');
if size(img,3) == 3
    img = rgb2gray(img);
end
img = double(img);

sigma = 20.0;
noise = sigma * randn(size(img));
img_noisy = img + noise;

[height, width] = size(img_noisy);
N = height * width;

lambda = 1.0;
rho = 1.0;
numIter = 50;

I_h = speye(height);
D0_v = -I_h + circshift(I_h, -1, 1); 
I_w = speye(width);
D0_w = -I_w + circshift(I_w, -1, 1);

Dv = kron(I_w, D0_v);

Dh = kron(D0_w, I_h);
D = [Dv; Dh];


tilde_x = reshape(img_noisy, [N, 1]);

x = tilde_x;
z = zeros(size(D,1), 1);
u = zeros(size(D,1), 1);

I_N = speye(N);
Dt = D';
DtD = Dt * D;
A = I_N + rho * DtD;  

for k = 1:numIter
    b = tilde_x + rho * (Dt * (z - u));
    x = A \ b; 

    Dx = D * x; 
    z = (rho / (2*lambda + rho)) * (Dx + u);

    u = u + (Dx - z);
end

img_denoised = reshape(x, [height, width]);

figure(1);
imshow(uint8(img), [0 255]);

figure(2)
imshow(uint8(img_noisy), [0 255]);

figure(3)
imshow(uint8(img_denoised), [0 255]);
```
| BEFORE | AFTER |
| -- | -- |
| ![スクリーンショット 2025-04-06 13.38.32.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3933005/eb190ae1-899d-4551-a851-54fcd09cde70.png)|![スクリーンショット 2025-04-06 13.38.36.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3933005/5eb874c2-a657-4e38-99be-04e6a15676a3.png)| 


## 5. 議論はあるか？
- **収束速度の改良:** ADMM自体をネステロフ加速する試み（Accelerated ADMM）や、**過緩和 (Over-relaxation)** や**線形化 (Linearized ADMM)** によるステップ高速化などが提案されています。
    - [Ouyang et al. (2015) Accelerated ADMM](https://arxiv.org/abs/1401.6607)
- **マルチブロックADMM:** 3つ以上の部分に対するADMM拡張は収束保証が難題でしたが、近年では**ランダムな更新順序**や**近似修正**によって安定化する方法が研究されています。
    - [A Survey of ADMM Variants for Distributed Optimization: Problems, Algorithms and Features](https://arxiv.org/abs/2208.03700)
    - [Eckstein, Bertsekas Generalized ADMM](https://arxiv.org/abs/1705.06191)
- **非凸問題への適用:** 深層学習の重み最適化など非凸問題にADMMを使う試みも増えています。理論保証は難しいものの、**経験的な成功例**（圧縮センシングの非凸ペナルティやディープラーニングの正則化最適化など）が報告されており、ADMMの適用範囲拡大が図られています。
- **ADMMと機械学習の融合:** ADMMの計算ステップをニューラルネットワークで展開した**深層展開 (Deep Unrolling)が注目されています。これは反復アルゴリズムを有限ステップのニューラルネットとして解釈し、データ駆動で学習させる手法です。ADMM-Netや学習型ADMM**が提案され、画像再構成や通信最適化で従来より高速な収束を実現しています。

## 参考文献
- [Boyd et al., 2011 – Distributed Optimization and Statistical Learning via the ADMM](https://web.stanford.edu/~boyd/papers/pdf/admm_distr_stats.pdf)
- [Glowinski & Marroco, 1975](https://www.ijcai.org/proceedings/2019/0629.pdf) / [Gabay & Mercier, 1976](https://www.sciencedirect.com/science/article/pii/0898122176900031?via%3Dihub)
- [Goldstein & Osher, 2009 – The Split Bregman Method for L1 Regularized Problems]()
- [Beck & Teboulle, 2009 – A Fast Iterative Shrinkage-Thresholding Algorithm](https://ieeexplore.ieee.org/document/4959678)
- [Chambolle & Pock, 2011 – A First-Order Primal-Dual Algorithm for Convex Problems](https://link.springer.com/article/10.1007/s11425-015-5079-0)
- https://mirror2image.wordpress.com/2013/01/02/split-bregman-for-total-variation-parameter-choice/#:~:text=Interesting,a%20panacea%2C%20I%20had%20great
- [[2208.03700] A Survey of ADMM Variants for Distributed Optimization: Problems, Algorithms and Features](https://arxiv.org/abs/2208.03700)
- [Ouyang et al. (2015)の**Accelerated ADMM**](https://arxiv.org/abs/1401.6607)
- [EcksteinとBertsekasの**Generalized ADMM (GADMM)](https://arxiv.org/abs/1705.06191)