交互方向乗数法による最適化と画像ノイズ除去への応用

はじめに

これまでの記事で近接勾配法と、それによるスパース解や低ランク解に導く正則化項を付随した最適化問題の解法、そしてその応用を見てきました。正則化項に変数間の絡みがなく各変数が独立に扱える場合は、正則化項のproximal operatorが解析的に求まるため近接勾配法は有効な手段です。ところが、各変数間に絡みがあり分離できない場合、一般的にproximal operatorの計算は容易ではありません（解析的に求まらない）。そこで今回は、変数間に絡みがある正則化項が付いていても最適解を導出することができる交互方向乗数法（Alternating Direction Method of Multipliers : ADMM）のアルゴリズムを見ていきます。また、それを用いた画像のノイズ除去をPythonで実装したので紹介します。

交互方向乗数法（ADMM）とは

いきなりですがADMMのアルゴリズムを記載します。まずADMMの対象となる最適化問題は下記のようなものになります。

$\displaystyle \min_{{\bf x},{\bf y}} f({\bf x}) + g({\bf y}) \ \ {\rm s.t.} \ \ {\bf A}{\bf x} + {\bf B}{\bf y} = {\bf 0} \tag{1}$

この問題を解くために、下記に示す「拡張ラグランジュ関数」というものを定義します。

$\displaystyle L_{\rho}({\bf x},{\bf y}, {\bf \lambda})= f({\bf x}) + g({\bf y}) + {\bf \lambda}^T ({\bf A}{\bf x} + {\bf B}{\bf y}) + \frac{\rho}{2} \| ({\bf A}{\bf x} + {\bf B}{\bf y}) \|^2_2 \tag{2}$

この拡張ラグランジュ関数に関する詳細は補足の項や参考文献をご参照ください。とりあえず、ここでは通常のラグランジュ関数に等式制約が満たされないことに対する罰則項（第3項）がついたものと考えて問題ありません。この拡張ラグランジュ関数を用いると式（1）を解くADMMアルゴリズムは下記のようになります。

${\bf x}$ 、 ${\bf y}$ 、 ${\bf \lambda}$ を適当な値で初期化する
以下の操作を収束するまで繰り返す

$\displaystyle {\bf x}_k \leftarrow {\rm arg}\min_{\bf x} L_{\rho} ({\bf x}, {\bf y}_{k-1}, {\bf \lambda}_{k-1} ) \tag{3}$
$\displaystyle {\bf y}_k \leftarrow {\rm arg}\min_{\bf y} L_{\rho} ({\bf x}_k, {\bf y}, {\bf \lambda}_{k-1} ) \tag{4}$
$\displaystyle {\bf \lambda}_k \leftarrow \lambda_{k-1} + \rho({\bf A}{\bf x}_k + {\bf B}{\bf y}_k) \tag{5}$

これだけです。アルゴリズムとしては非常に簡単です。特に式（3）と式（4）では ${\bf x}$ と ${\bf y}$ それぞれに関する最適化問題を独立で解いているので、元々の最適化問題が関数同士の和が含まれており解くのが難しい場合でも、式（1）の形式に変換できればADMMにより簡単に解くことができます。このことを次の章で見ていきます。が、その前に補足として、上記アルゴリズムがどのように出てきたかを理解するために、双対上昇法、拡張ラグランジュ関数について簡単に説明します（細かいことはかなり省略するので、詳しい説明は参考文献をご参照ください）。

【補足】双対上昇法、拡張ラグランジュ関数について

ADMMは最適化問題の強双対性という性質を利用しています。ここでは、

$\displaystyle \min_{{\bf x}} f({\bf x}) \ \ {\rm s.t.} \ \ {\bf A}{\bf x} = {\bf b} \tag{6}$

という最適化問題を考えます。この問題のラグランジュ関数は、

$\displaystyle L({\bf x}, {\bf \lambda})= f({\bf x}) + {\bf \lambda}^T( {\bf A}{\bf x} - {\bf b} )\tag{7}$

と定義されます。 $f({\bf x})$ が真凸関数である場合は、 $\nabla L({\bf x}^*, {\bf \lambda}^*) = 0$ となる ${\bf x}^*$ と ${\bf \lambda}^*$ が存在すれば、 ${\bf x}^*$ が式（6）の最適解です。
更に $f({\bf x})$ が真凸関数であるとき式（7）について下記も成り立ちます。

$\displaystyle \min_{{\bf x}} \max_{{\bf \lambda}} L({\bf x}, {\bf \lambda}) = \max_{{\bf \lambda}} \min_{{\bf x}} L({\bf x}, {\bf \lambda})\tag{8}$

これは強双対性と呼ばれる性質で、左辺は式（6）と同値の問題です。この強双対性を利用して最適化問題を解くのが双対上昇法や拡張ラグランジュ関数法、そして今回の主題であるADMMです。

双対上昇法

双対上昇法は式（6）の問題を強双対性から双対問題という問題に変換し、そちらを解こうという方針の手法です。まず、式（8）の右辺に問題を当てはめ変形すると、

$\displaystyle \begin{eqnarray} \max_{{\bf \lambda}} \min_{{\bf x}} L({\bf x}, {\bf \lambda}) &=& \max_{{\bf \lambda}} -f^*(-{\bf A}^T {\bf \lambda}) - {\bf \lambda}^T {\bf b} \\ & =& \min_{{\bf \lambda}} f^*(-{\bf A}^T {\bf \lambda}) + {\bf \lambda}^T {\bf b}\\ & =& \min_{{\bf \lambda}} \phi({\bf \lambda}) \end{eqnarray} \tag{9}$

という問題になります。ここで、 $f^*$ は以前説明した共役関数で $f^*({\bf z})= \sup_{{\bf x}} {\bf z}^T {\bf x} - f({\bf x})$ です。この式（9）の問題は双対問題と呼ばれ、これを解く方法を双対上昇法と言います。更にこれを勾配法に基づいて解く場合は宇沢の方法といい、

$\displaystyle {\bf \lambda}_k = {\bf \lambda}_{k - 1} - \alpha \nabla \phi({\bf \lambda}_k) \tag{10}$

という更新式になります。で、この勾配を求めるためには、

$\displaystyle \nabla f^*(-{\bf A}^T {\bf \lambda})= {\rm arg} \min_{\bf x}f(\bf x) + {\bf \lambda}_k^T {\bf A} {\bf x}\tag{11}$

を計算する必要があります。結局のところこれは、式（7）にあるラグランジュ関数 $L({\bf x}, {\bf \lambda}_k)$ の ${\bf x}$ に関する最小化です。つまり、主問題である ${\bf x}$ に関する最適化と双対問題である ${\bf \lambda}$ に関する最適化を交互に行うアルゴリズムとなります。
この双対上昇法では共役関数 $f^*$ の勾配を求めるため、共役関数が微分可能でなければ使えません。しかし、一般的に共役関数は微分可能ではなく、しかも関数 $f$ が強凸という性質*1を持っていなければアルゴリズムの収束も保証されません*2。そこで、次に説明する拡張ラグランジュの登場です。

拡張ラグランジュ関数法

まず式（6）の問題について、

$\displaystyle \min_{{\bf x}} f({\bf x}) + \frac{\rho}{2} \| {\bf A}{\bf x} - {\bf b} \|_2^2\ \ {\rm s.t.} \ \ {\bf A}{\bf x} = {\bf b} \tag{12}$

という式に書きなおします。目的関数の第２項として $\frac{\rho}{2} \| {\bf A}{\bf x} - {\bf b} \|_2^2$ を付け加えましたが制約条件より $0$ になるため、式（12）は式（6）と等価です。次に式（12）のラグランジュ関数は

$\displaystyle L_{\rho}({\bf x}, {\bf \lambda})= f({\bf x}) + {\bf \lambda}^T ( {\bf A}{\bf x} - {\bf b}) + \frac{\rho}{2} \| {\bf A}{\bf x} - {\bf b} \|_2^2 \tag{13}$

となります。これは拡張ラグランジュ関数と呼ばれる関数で、ラグランジュ関数に対し制約を満たさないことに対する罰則項 $\frac{\rho}{2} \| {\bf A}{\bf x} - {\bf b} \|_2^2$ を加えたものとなっています。実はこの罰則項があることにより主問題は強凸となり、たとえ $f$ の共役関数 $f^*$ が微分不可能だったとしても双対問題が滑らかになるため双対上昇法が適用できます（その理由については、参考文献をご参照ください）。このように、ラグランジュ関数に制約条件を満たさないことに対する罰則項 $\frac{\rho}{2} h({\bf x} )$ を加えることで問題を解きやすくした上で、 ${\bf x}$ と ${\bf \lambda}$ の更新を交互していくのが拡張ラグランジュ関数法です。

ADMM再考

さて、式（1）は式（6）において変数が ${\bf x}$ と ${\bf y}$ に分解可能という特別なケースです。式（1）を拡張ラグランジュ関数法で解こうとした時に ${\bf x}$ と ${\bf y}$ に関する同時最適化が困難なことが多いです。そこで1回の更新で ${\bf x}$ と ${\bf y}$ の同時最適化はあきらめ、個別に最適化していこうというのがADMMです。

変数間に絡みのある正則化項付き最適化問題のADMMによる解法

L1ノルム正則化は変数間が独立しており、proximal operatorが解析的に求まります。L2ノルム正則化（２乗しない）はルートを取る段階で変数間の絡みが出てきますがMoreau decompositionをうまく使ってproximal operator計算することができます（近接勾配法の記事参照）。p=2の重複なしのグループ正則化については、重複がないため各グループ独立したL2ノルムとしてproximal operatorとして計算できます。
しかし、上述のように変数間に絡みがあるとproximal operatorは容易には計算できません。例えば、重複ありのグループ正則化（p=2）の場合、同じ変数が複数のグループにまたがって現れるため、各グループ独立してproximal operatorを計算できなくなってしまいます。そこで、そのような正則化がつく最適化問題をうまく式（1）の形式かつ $g({\bf y})$ がproximal operatorが計算しやすい形に変形し、ADMMで解くことを考えます。重複ありグループ正則化を例にとって考えてみると、

$\displaystyle \min_{{\bf x}} f({\bf x}) +C\sum_{g \in G} \| {\bf x}_g \|_2 \tag{14}\ \ \ \ \ \ \ \ \ G=\left\{g_1, g_2, \cdots, g_k\right\}$

という問題を解くことになるのですが、これを下記のような問題に変形します。

$\displaystyle\min_{{\bf x}} f({\bf x}) +C \sum_{g' \in G'} \| {\bf y}_{g'} \|_2 \ \ {\rm s.t.} \ \ {\bf y} = {\bf B}{\bf x}\tag{15}$

ここで ${\bf B}$ をうまく定義し $G'$ がグループ間の重複がないようにします。具体的には ${\bf B}_{g_1},{\bf B}_{g_2}, \cdots, {\bf B}_{g_k}$ を

$\displaystyle {\bf B}_{g_1} = \left[ \begin{array}{ccccccc} 1 & 0 & 0 & 0 & 0 & \ldots & 0 \\ 0 & 1 & 0 & 0 & 0 & \ldots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & 0 & \ldots & 1 \\ \end{array} \right] ,\cdots, \displaystyle {\bf B}_{g_k} = \left[ \begin{array}{ccccccc} 1 & 0 & 0 & 0 & 0 & \ldots & 0 \\ 0 & 0 & 0 & 1 & 0 & \ldots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & 0 & \ldots & 1 \\ \end{array} \right]$

などのように、それぞれ ${\bf x}_g={\bf B}_g {\bf x}$ と左からかけると、グループ $g$ に所属する要素のみを抽出する行列とします。そして、 ${\bf B}=\left[{\bf B}_{g_1}^T, {\bf B}_{g_2}^T, \cdots, {\bf B}_{g_k}^T \right]^T$ とすると、式（15）より ${\bf y}={\bf B}{\bf x}=\left[{\bf x}_{g_1}^T, {\bf x}_{g_2}^T, \cdots, {\bf x}_{g_k}^T\right]^T$ と ${\bf y}$ は各グループに対応する部分ベクトルを並べたベクトルとなります。こうすることで、正則化項 $C \sum_{g' \in G'} \| {\bf y}_{g'} \|_2$ は重複なしのグループ正則化になります。
ということで、この式（15）に対してADMMを適用することを考えます。この中で ${\bf y}$ に関する更新式は、

$\displaystyle \begin{eqnarray} {\bf y}_k &\leftarrow& {\rm arg}\min_{\bf y} L_{\rho} ({\bf x}_k, {\bf y}, {\bf \lambda}_{k-1} ) \\ & = & {\rm arg}\min_{\bf y} C \sum_{g' \in G'} \| {\bf y}_{g'} \|_2 + {\bf \lambda}_{k-1}^T ( {\bf y} - {\bf B}{\bf x}_k) + \frac{\rho}{2}\| {\bf y} - {\bf B}{\bf x}_k \|_2^2 \\ &=& {\rm arg}\min_{\bf y} \frac{1}{\rho} C \sum_{g' \in G'} \| {\bf y}_{g'} \|_2 + \frac{1}{2} \|{\bf y} -({\bf B} {\bf x}_k + {\bf \lambda}_{k-1}/\rho ) \|_2^2 \\ &=& {\rm prox}_{ \psi/\rho}({\bf B} {\bf x}_k + {\bf \lambda}_{k-1}/\rho) \end{eqnarray} \tag{16}$

となります。ここで、 ${\rm prox}_{ \psi}(\cdot)$ は、グループ正則化項 $C \sum_{g' \in G'} \| {\bf y}_{g'} \|_2$ のproximal operatorです（こちらの式（21）参照）。このように、元々の問題が正則化項のproximal operaorを求めるのが困難であり近接勾配法を利用できない場合でも、うまく式変形を行いADMMを適用すれば解けるようになることもあるため、ADMMはかなり強力なアルゴリズムだと言えます。

TV正則化によるノイズ除去への応用

冒頭でも述べたように、今回はこのADMMを応用して画像のノイズ除去を行ってみます。今、入力画像を ${\bf v}$ として、この ${\bf v}$ からのノイズ除去は、全変動正則化（TV正則化）と呼ばれる正則化項 $\|\cdot\|_{TV}$ を用いて以下のような定式化が提案されています。

$\displaystyle{\rm arg} \min_{\bf x} \|{\bf v} - {\bf x} \|_2^2 + C \| {\bf x}\|_{TV} \tag{17}$

式（17）の第１項は出力画像 ${\bf x}$ がなるべく入力画像 ${\bf v}$ に近する働きがあります。第２項がTV正則化項で下記のような式となります。

$\displaystyle \| {\bf x}\|_{TV} = \sum_{(i,j)}\sqrt{({\bf x}_{(i,j)} - {\bf x}_{(i + 1, j)})^2 + ({\bf x}_{(i,j)} - {\bf x}_{(i, j + 1)})^2} \tag{18}$

ここで $(i,j)$ はピクセル位置を表します*3。基本的に画像は滑らかな性質を有しているはずなので、隣り合うピクセルの差は小さいはずです。そのため、このTV正則化は隣り合うピクセルの差が大きい場合の罰則項となり、画像を滑らかにしノイズ成分を除去する働きがあります。
さて、この式（17）をADMMで解くことを考えます。TV正則化項も重複ありグループ正則化と同じく変数間に絡みがあり、直接proximal operatorを計算するのは困難です。よって、重複ありグループ正則化を考えたときと同じく、 ${\bf y} = {\bf B}{\bf x}$ として、重複なしグループ正則化項に変換することを考えます。これは以下のように、

$\displaystyle {\bf B}_{g_1} = \left[ \begin{array}{ccccccc} 1 & -1 & \cdots & 0 & 0 & \cdots & 0 \\ 1 & 0 & \cdots & -1 & 0 & \cdots & 0 \end{array} \right] , \displaystyle {\bf B}_{g_2} = \left[ \begin{array}{ccccccc} 0 & 1 & -1 & \cdots & 0 & \cdots & 0 \\ 0 & 1 & 0& \cdots & -1 & \cdots & 0 \end{array} \right] , \cdots$

と、 ${\bf B}_g {\bf x} = \left[{\bf x}_{(i,j)} - {\bf x}_{(i + 1, j)}, {\bf x}_{(i,j)} - {\bf x}_{(i, j+1)}\right]$ となるように定義すれば、式（18）は ${\bf B}_g {\bf x}$ のL2ノルムの和なので、 ${\bf B}=\left[{\bf B}_{g_1}^T, {\bf B}_{g_2}^T, \cdots, {\bf B}_{g_k}^T \right]^T$ とすればTV正則化を重複なしグループ正則化に変換することができます。よって、この ${\bf B}$ を用いれば式（16）の更新式で ${\bf y}$ を更新することができます。一方、 ${\bf x}$ についての更新式は、式（3）について ${\bf x}$ で微分して ${\bf 0}$ となる点を求めればよいので解析的に求まり、以下のようになります。

$\displaystyle \begin{eqnarray} {\bf x}_k &\leftarrow& {\rm arg}\min_{\bf x} L_{\rho} ({\bf x}, {\bf y}_{k-1}, {\bf \lambda}_{k-1} ) \\ & = & (2 {\bf I} + \rho {\bf B}^T {\bf B})^{-1} (2 {\bf v} - {\bf B}^T ({\bf \lambda}_{k-1} - \rho {\bf y}_{k-1})) \end{eqnarray} \tag{19}$

実装と実験

ノイズ除去アルゴリズムの実装

ということで、TV正則化によるノイズ除去を実装し、実験してみました。ここで、全てのピクセルを１つのベクトルとして扱ってしまうと画像サイズによってはメモリが足りなくなってしまうので、超解像を行ったときのように画像パッチを切り出し、それぞれのパッチに対してノイズ除去を行った後、そのパッチをつなぎ合わせて画像を再構築するという方法を取りました。なお、この画像切り出しやパッチのつなぎ合わせなどは、超解像の記事で記した「おれおれ画像・ベクトル・行列変換関数群」も利用してるので、そちらもご参照ください。

import os
import sys
import numpy as np
import scipy as sci
from scipy import sparse
from numpy.random import *
from PIL import Image

#ブロックソフト閾値計算
def block_soft_thresh(b, lam):
    return max(0, 1 - lam/np.linalg.norm(b, 2)) * b

#グループ正則化のproximal operator
def prox_group_norm(v, groups, gamma, lam):
    u = np.zeros(v.shape[0])
    for i in range(1, np.max(groups) + 1):
        u[groups == i] = block_soft_thresh(v[groups == i] , gamma * lam)
    return u

#グループ正則化計算
def gloup_l1_norm(x, groups, p):
    s = 0
    for i in range(1, np.max(groups) + 1):
        s += np.linalg.norm(x[groups == i], p)
    return s

#関数をメモ化する関数
def memoize(f):
    table = {}
    def func(*args):
        if not args in table:
            table[args] = f(*args)
        return table[args]
    return func

#ADMMを行う関数
#argmin_x, argmin_y, update_lam, objectiveはそれぞれ変数名に準ずる計算を行う関数
def ADMM(argmin_x, argmin_y, update_lam, p, objective, x_init, y_init, lam_init, tol= 1e-9):
    x = x_init
    y = y_init
    lam = lam_init
    result = objective(x)
    while 1:
        x_new = argmin_x(y, lam, p)
        y_new = argmin_y(x_new, lam, p)
        lam_new = update_lam(x_new, y_new, lam, p)
        result_new = objective(x_new)
        if result_new < tol or (np.abs(result - result_new)/np.abs(result) < tol) == True :
            break
        x = x_new
        y = y_new
        lam = lam_new
        result = result_new
    return x_new, result_new


#TV正則化付きのスムージング
#N,M:画像のピクセル行数、列数
#v:入力画像のベクトル
#B:正則化の変換行列（スパースなので、scipyのsparse行列を渡す）
#groups:変換後ベクトルの各要素の所属グループ
#C:正則化の係数
#p：拡張ラグランジュの係数
def TV_reg_smoothing(N, M, v, B, groups, C, p) :
    #inv(2I + pB^T B)を計算する関数。アルゴリズムによってはpが可変なので、pを引数として受け取る。
    #関数をメモ化して同一のpが入力された場合、再計算不要としている。
    inv_H = memoize(lambda p:np.array(np.linalg.inv( 2.0 * np.eye(B.shape[1], B.shape[1]) + p * (B.T * B))))

    argmin_x = lambda y, lam, p:np.dot(inv_H(p), 2.0 * v - np.array(B.T * (-p * y + lam)))
    argmin_y = lambda x, lam, p:prox_group_norm((B * x) + lam/p, groups, 1.0/p, C)
    update_lam = lambda x, y, lam, p:lam + p*((B * x) - y)
    objective = lambda x:np.linalg.norm(v - x, 2)**2 + C * gloup_l1_norm((B * x), groups, 2)

    x_init = np.random.randn(B.shape[1])
    y_init = (B * x_init)
    lam_init = np.zeros(B.shape[0])

    (x, result) = ADMM(argmin_x, argmin_y, update_lam, p, objective, x_init, y_init, lam_init, 1e-9)

    return x, result

#グループ変換行列を計算する関数
#N,M:画像のピクセル行数、列数
def calc_group_trans_mat(N, M):
    B = sci.sparse.lil_matrix((2 * (M - 1) * (N - 1) + (M - 1) + (N - 1), M * N))
    groups = np.zeros(B.shape[0], 'int')
    k = 0
    for i in range(N):
        for j in range(M):
            base = i * M + j
            if i < N -1 and j < M -1:
                B[k, base] = 1
                B[k, base + 1] = -1
                B[k + 1, base] = 1
                B[k + 1, base + M] = -1
                groups[k] = int(k/2)  + int(k % 2) + 1
                groups[k + 1] = int(k/2)  + int(k % 2) + 1
                k += 2
            #一番下の行のピクセルは右隣のピクセルとの差分のみ計算
            elif i >= N - 1 and j < M - 1:
                B[k, base] = 1
                B[k, base + 1] = -1
                groups[k] = int(k/2) + int(k % 2) + 1
                k += 1
            #一番右の劣のピクセルは下のピクセルとの差分のみ計算
            elif i < N - 1 and j >= M - 1:
                B[k, base] = 1
                B[k, base + M] = -1
                groups[k] = int(k/2) + int(k % 2) + 1
                k += 1
    return B, groups

#グレースケール画像に対するノイズ除去
#img:画像データ（PILのimageクラス）
#patch_hight:切り出す画像パッチの高さ
#patch_width:切り出す画像パッチの幅
#shift_num:画像を切り出す際のずらし量（重複して切り出してもOK）
#C:正則化の係数
#p:拡張ラグランジュの係数
def denoise_gray_img(img, patch_hight, patch_width, shift_num, C, p):
    #グループ変換行列の計算
    [B, groups] = calc_group_trans_mat(patch_hight, patch_width)

    #画像パッチの切り出し
    patchs = gray_im2patch_vec(img, patch_hight, patch_width, shift_num)

    new_patchs = np.zeros((patch_hight * patch_width, patchs.shape[1]))
    for i in range(patchs.shape[1]):
        #各パッチに対してノイズ除去を施す
        new_patchs[:, i] = TV_reg_smoothing(patch_hight, patch_width, patchs[:, i], B, groups, C, p)[0]

    #パッチをつなぎ合わせて画像の再構成
    [new_img, img_arr] = patch_vecs2gray_im(new_patchs, img.size[0], img.size[1], patch_hight, patch_width, shift_num)
    return new_img

細かい点はソースコードとそのコメントをご参照ください。とりあえず、denoise_gray_imgにPILのimageクラスのオブジェクトと切り出すパッチのサイズ、パッチ切り出し時のずらし量を渡せばノイズ除去されたimageクラスのオブジェクトを返してくれるように実装しています。超解像の時と同様、パッチを切り出す際に領域の重複を許し、再構成時に重ね合わせて平均を取るようにしています。また、画像の最後の行および列は右隣のピクセル（および下のピクセル）とのみ差分を取るようにBを工夫しています（calc_group_trans_mat関数内）。
注意する点としては、TV_reg_smoothingが受け取るBはスパースな行列なので、scipyのsparse行列のクラス渡します。また、同関数内にあるinv_Hはpを引数にとり、式（19）の逆行列部分を計算する関数ですが、同一のpに対して再度同じ計算を行うと時間がかかるので、memoizeによりメモ化し同じpが入力されたときに再計算しないようにしています。
なお、更新式はADMM関数の中に直接記述するのではなく、ADMM関数はそれぞれの更新式の計算を行う関数を引数として受け取るようにし、目的関数の出力が収束するまでそれらの更新式関数を呼び出すという設計にしているため、更新式の定義もTV_reg_smoothing関数の中で行っています。

実験と結果

いつものように256×256のレナさん画像で実験を行います。
f:id:YamagenSakam:20180630185527p:plain

これに対し、標準偏差30のガウスノイズを乗せた画像に対してノイズ除去を行ってみました。

def main():
    test_img = Image.open("ADMM_input.png")
    img_arr = np.array(test_img.convert('L'),'float')
    img_arr = (img_arr  + 30 * (np.random.randn(test_img.size[0], test_img.size[1]) ))
    img_arr[img_arr >255] = 255
    img_arr[img_arr <0] = 0
    test_img = Image.fromarray(np.uint8(img_arr))
    patch_hight = 20
    patch_width = 20
    shift_num = 10
    C = 50
    p = 20
    img2 = denoise_gray_img(test_img, patch_hight, patch_width, shift_num, C, p)
    img2.save("ADMM_result.png")

画像パッチのサイズを20×20ずらし量10として切り出し、C=50、p=20で実験を行った結果は以下のようになります。

　
左が入力画像（標準偏差30のガウスノイズを重畳）、右がノイズ除去後の画像

元画像レベルとまではいきませんが、それなりにノイズ除去ができているかと思います。ただ、パッチを切り貼りしている影響か、少しゆがんでしまっているようにも見えます。更に追実験として入力画像のノイズレベルを標準偏差50と上げ、その時のノイズ除去もやってみました（C=80、p=1で実験）。

左が入力画像（標準偏差50のガウスノイズを重畳）、右がノイズ除去後の画像

ノイズレベルの割にはある程度ノイズ除去できているかと思いますが、TV正則化項の罰則に引っ張られ画像がノッペリとしてしまいました。Cをもっと調整すればもう少し改善できるかも知れません。

まとめ

今回は近接勾配法よりも幅広い問題を解くことができるADMMについて説明し、それを用いて画像のノイズ除去をやってみました。今回紹介した重複ありグループ正則化やTV正則化以外にも、特徴間のグラフ構造や階層構造なんかを表す正則化項がついた最適化問題もADMMを使えば解けるようになり、加えて別途制約条件がある場合でも式変形をうまくやればADMMで解けるようになるケースもあるなど応用範囲が広い手法です。更に収束も割りと速く、 $\rho$ の値も適当に決めても $0$ より大きければ大体収束するので、使い勝手もなかなか良いです。