2018-04-12

近接勾配法応用編その1 ～スパースコーディング、辞書学習からの超解像～

はじめに

前回の記事で近接勾配法の基本と実装を見てきました。今回はこの近接勾配法を利用して、スパースコーディングの応用の１つである超解像技術を調べPythonで実装してみました*1。画像処理関係の実装はほぼ未経験だったので、その分野の人が見れば前処理等が足りなかったり、アンチパターン的なことをやってたりするかもしれませんが、その点は指摘頂けると幸いです。

スパースコーディングとは

スパースコーディングは観測した信号 ${\bf y}$ を、「辞書」と呼ばれる $m$ 個の基底ベクトルの集合　 ${\bf D} = \left[{\bf d}_1, {\bf d}_2, \cdots, {\bf d}_m \right]$ のスパースな線形結合で表現する符号化手法です。
今、この線形結合の係数を ${\bf w}$ とすると、スパースコーディングは下記のように定式化されます。
$\displaystyle {\bf y }\simeq {\bf D} {\bf w} \tag{1}$ 　
ただし、 ${\bf w}$ の要素の多くは $0$ （つまり ${\bf w}$ はスパース）です。

このようなスパースコーディングを使うことによりデータの圧縮を行えるだけでなく、画像からのノイズ除去や今回取り上げる超解像など様々な技術へ応用が可能です。
スパースコーディングを成功させる上で重要となるのは辞書 ${\bf D}$ です。 ${\bf D}$ がテキトーだとうまくいきません*2。ということで、辞書 ${\bf D}$ は実データからの学習により構築します。ここで、学習データの個数を $n$ 、学習データをまとめた行列 ${\bf Y} = \left[{\bf y}_1 \cdots, {\bf y}_n\right]$ 、各線形結合の係数をまとめた行列 ${\bf W} = \left[{\bf w}_1, \cdots, {\bf w}_n\right]$ を定義すると、辞書学習の目的関数は下記のようになります。
$\displaystyle {\rm arg} \min_{{\bf D}, {\bf W}} \| {\bf Y} - {\bf D}{\bf W} \|^2_2 \tag{2}$ 　
ただし、 ${\bf W}$ はスパース。

この方法で辞書 ${\bf D}$ が得られたら、次はその辞書を用いて入力画像の符号化します。そして、その係数と辞書から画像を再構成し超解像画像を得るという流れになります。
ということで、スパースコーディングによる超解像生成は（1）辞書を学習するフェーズと（2）その辞書を元に画像を再構成するフェーズの2フェーズで実現できます。

画像再構成のための辞書学習

実装した手法の説明

では、今回実装した学習用の画像データを元に辞書を学習する方法を説明します。
まず、最終的に得たい解像度の画像データ（高解像画像データ）複数枚を学習データとして用意します。この時、超解像を行いたい画像のカテゴリが決まっているなら、そのカテゴリに合った画像を用意するのがよさそうです。
この学習データを用いて、辞書 ${\bf D}$ を学習する流れを下記に示します。

学習用画像を小領域を取り出した「画像パッチ」を $n$ 個用意し、それぞれ縦に並べたベクトル ${\bf y}_1, {\bf y}_2, \cdots, {\bf y}_n$ および、これらのベクトルを横に並べた行列を ${\bf Y} = \left[{\bf y}_1 \cdots, {\bf y}_n\right]$ を作る。
${\bf Y}$ の各行が平均 $0$ 、分散 $1$ になるように標準化する。
${\hat {\bf D}}$ をランダム行列で初期化する。
${\hat {\bf D}}$ を固定して、次の最適化問題を $i = 1, 2, \cdots, n$ について解き、行列 ${\hat {\bf W}} = \left[{\hat {\bf w}_1}, \cdots, {\hat {\bf w}_n} \right]$ を得る。

$\displaystyle {{\hat {\bf w}_i}} = {\rm arg} \min_{{\bf w}_i} \| {\bf y}_i - {\hat {\bf D}}{\bf w}_i \|^2_2 + \lambda \| {\bf w}_i \|_1\tag{3}$

$\hat{{\bf W}}$ を固定して、次の最適化問題を解き $\hat{{\bf D}}$ を得る。

$\displaystyle \hat{{\bf D}}= \arg \min_{{\bf D}} \| {\bf Y} - {\bf D} {\hat{ \bf W}} \|^2_2 \tag{4}$

手順4,手順5の操作を収束するまで、もしくは一定回数繰り返す。

ここで、式（3）はlassoの式であるため解がスパースになることが期待され、更に近接勾配法で解くことができます。なお、今回、近接勾配法の応用編ということで手順4の方法でスパースな係数 $\bf{w}_i$ を求めていますが、この方法以外にもOrthogonal Matching PursuitやIterative Reweighted Least Squaresなどのアルゴリズムがあります。
また、手順5で式（4）の解を得る方法としてはK-SVDという手法がよく使われますが、今回は簡単のために下記のように解析的に解いた結果を使用します。
$\displaystyle {\hat {\bf D}} = {\bf Y} {\hat {\bf W}}^T ( {\hat {\bf W}} {\hat {\bf W}}^T )^{-1} \tag{5}$

この手順4、手順5を繰り返して ${\bf W}, {\bf D}$ の同時最適化を行うわけですが、この両変数に対しては凸関数ではないため（式（3）、式（4）はそれぞれ凸関数ですが）、求まる結果は局所最適解です。
なお、この学習を収束まで待っているとかなり時間がかかるため、今回は適当な回数ループさせたら学習を完了させています。

ちなみに、最初は手順2の標準化を省いた実装をしていましたが、これだとうまく超解像ができませんでした。やはり、学習データの標準化は重要のようです。

実装したコード

ここまでの手順を実現するpythonコードです。

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

#dorectory以下のファイルを再帰的に検索する関数
def find_all_files(directory):
    for root, dirs, files in os.walk(directory):
        yield root
        for file in files:
            yield os.path.join(root, file)

#おれおれ画像・ベクトル・行列変換関数群
def color2gray_img(img):
    return img.convert('L')

def im2vec(img):
    img_arr = np.array(img.convert('L'),'float')
    return img_arr.reshape(-1)

def gray_im2float_vec(gray_img):
    return im_scale2float((im2vec(gray_img)))

def color_im2float_vec(color_img):
    return gray_im2float_vec(color2gray_img(color_img))

def vec2image(vec, H, W):
    tmp_vec = copy.deepcopy(vec)
    tmp_vec = tmp_vec.reshape((H, W))
    return Image.fromarray(np.uint8(tmp_vec))

def float_vec2image(vec, H, W):
    return vec2image(float2im_scale(vec), H, W)

def arr2image(arr):
    arr = copy.deepcopy(arr)
    return Image.fromarray(np.uint8(arr))

def float_arr2image(arr):
    return arr2image(float2im_scale(arr))

def im_scale2float(vec):
    return vec.astype(np.float)

def float2im_scale(vec):
    vec[vec > 255.0] = 255.0
    vec[vec < 0.0] = 0.0
    return vec

#画像データからパッチを切り出す関数
#shift_numで切り出す際のずらし量を指定
def gray_im2patch_vec(im, patch_hight, patch_width, shift_num):
    patchs = np.zeros((patch_hight*patch_width,
                        (1 + (im.size[0] - patch_hight)/shift_num + (1 if (im.size[0] - patch_hight) % shift_num  != 0 else 0)) *
                        (1 + (im.size[1] - patch_width)/shift_num + (1 if (im.size[1] - patch_width) % shift_num  != 0 else 0))))
    im_arr = im_scale2float(np.array(im))
    k = 0
    for i in range(1 + (im.size[0] - patch_hight)/shift_num + (1 if (im.size[0] - patch_hight) % shift_num  != 0 else 0)):
        h_over_shift = min(0, im.size[0] - ((i* shift_num) + patch_hight)) #ずらしたことによって画像からはみ出る量(縦方向)
        for j in range(1 + (im.size[1] - patch_width)/shift_num + (1 if (im.size[1] - patch_width) % shift_num  != 0 else 0)):
            w_over_shift = min(0, im.size[1] - ((j* shift_num) + patch_width)) #ずらしたことによって画像からはみ出る量(横方向)
            #パッチの切り出し　画像からはみ出たら、そのはみ出た分を戻して切り出し
            patchs[:, k] = im_arr[(i* shift_num) + h_over_shift:((i* shift_num) + patch_hight) + h_over_shift,
                                    (j* shift_num) + w_over_shift:((j* shift_num) + patch_width) + w_over_shift].reshape(-1)
            k += 1
    return patchs


#式（5）を計算する関数
def dictionary_learning_with_mod(W, Y):
    D = np.dot(np.dot(Y, W.T), np.linalg.inv(np.dot(W, W.T)))
    result = np.linalg.norm(Y - np.dot(D, W)) ** 2
    return D, result

#辞書クラス
class dictionary_with_learning():
    def __init__(self, weight_learner, dictionary_learner, base_num, max_loop_num = 8):
        self.Dic = np.array([])
        self.weight_learner = weight_learner
        self.dictionary_learner = dictionary_learner
        self.base_num = base_num

    #学習データ行列Yから辞書を学習するメソッド
    def learn(self, Y, eps = 1e-8):
        #手順3:ランダム行列で辞書の初期化
        D = np.random.randn(Y.shape[0], self.base_num)

        result = -1
        loop_num = 0
        while 1:
            loop_num += 1
            W = np.zeros((self.base_num, Y.shape[1]))
            #手順4:Dを固定して係数を求める
            for i in range(Y.shape[1]):
                [W[:, i], tmp] = self.weight_learner(D, Y[:, i])

            #手順5:Wを固定してDを求める
            [D, new_result] = self.dictionary_learner(W, Y)

            if result < 0 or np.abs(new_result - result) > eps or loop_num < max_loop_num:
                result = new_result
            else: #手順4、5を収束するまで繰り返す。
                break;
        self.Dic = D
        return D

    def get_dic_array(self):
        return self.Dic



#画像パッチから学習される辞書クラス
class image_patch_dictionary():
    def __init__(self, dictionary):
        self.dictonary = dictionary
        self.patch_hight = 0
        self.patch_width = 0

    #学習画像の集合から辞書を学習するメソッド
    def learn_images(self, image_list, patch_hight, patch_width, eps=2e-2, shift_num=1):
        self.patch_hight = patch_hight
        self.patch_width = patch_width
        Y = np.empty((patch_hight * patch_width, 0))
        #手順1:画像パッチの切り出し
        for img in image_list:
            gray_img = color2gray_img(img)
            Ytmp = gray_im2patch_vec(gray_img, patch_hight, patch_width,shift_num) 
            Y = np.concatenate((Y, Ytmp), axis=1)

        #手順2:学習データ行列の標準化
        Y_ = (Y - np.mean(Y, axis=1, keepdims = True))/np.std(Y, axis = 1, keepdims=True)
        self.dictonary.learn(Y_, eps)

    def get_dic(self):
        return self.dictonary.get_dic_array()

    def get_patch_hight(self):
        return self.patch_hight

    def get_patch_width(self):
        return self.patch_width


def main():
    lam = 5
    base_num = 144
    patch_hight = 12
    patch_width = 12
    train_img_list = []
    for file in find_all_files('C:\\data'):
        if os.path.isfile(file):
            try:
                train_img_list.append(Image.open(file))
                if len(train_img_list) >= 50:
                    break;
            except:
                print file

    weight_learning_method = lambda W, y:ISTA(W, y, lam) #近接勾配法を行う関数を使用
    dic_learning_method = lambda W, Y:dictionary_learning_with_mod(W, Y)#解析的に辞書Dを求める方法を使用
    im_dic = image_patch_dictionary(dictionary_with_learning(weight_learning_method, dic_learning_method, base_num))
    im_dic.learn_images(train_img_list, patch_hight, patch_width, shift_num = patch_hight)

    np.save('dic_array.npy', im_dic.dictonary.Dic)

本実装は、画像⇔行列・ベクトル変換を行う関数群、パッチ切り出しを行う関数、画像パッチ辞書クラスimage_patch_dictionary、辞書クラスdictionary_with_learningから構成されています。
image_patch_dictionaryは実際の辞書学習オブジェクトであるdictionary_with_learningをインスタンス変数として持ち、与えられた学習画像群からパッチの切り出しを行い、辞書を学習させる機能を持ちます。
そのため、辞書学習を行うのはdictionary_with_learningになりますが、係数の学習を行うweight_learnerと辞書の学習を行うdictionary_learnerをインスタンスとして持ち、生成時にどのアルゴリズムで係数推定、辞書推定を行うかを指定します。今回、係数の計算には近接勾配法を使用しますので、前回の記事記載のISTA関数をweight_learnerとして指定し、dictionary_learnerには新たに作成したdictionary_learning_with_mod関数を指定します（main関数参照）。もし係数推定として別のアルゴリズムを使いたければweight_learnerを変えればよいし、辞書学習もK-SVDなどを使いたければdictionary_learnerを変更すればよいです。

画像再構成による超解像

実装した手法の説明

ここからは、スパースコーディングによる超解像の手法の説明です。手法の大まかな考え方は「学習した辞書 ${\bf D}$ をダウンサンプリングした辞書 ${\bf D}_l$ で入力画像を符号化する係数を求めた後、その係数と元の辞書 $\bf D$ から画像を再構成する」というものです。以下、順を追って説明します。

画像パッチの切り出し

${\bar{\bf Y}}$

${\bar{\bf y}_i}$

$h \times w$

$s$

$\frac{h}{s} \times \frac{w}{s}$

${\bar {\bf y}_i}$

$\frac{hw}{s^2}$

アーティファクト

辞書のダウンサンプリング

${\bf D}$

作用素

${\bf D}_l$

${\bf D}$

$hw$

$\times$

$m$

${\bf d}$

$h \times w$

$s$

${\bf D}_l$

$\frac{hw}{s^2} \times m$

符号化係数の推定

${\bf D}_l$

最適化問題

$\displaystyle {{\hat {\bf w}_i}} = {\rm arg} \min_{{\bf w}_i} \| {\bar {\bf y}_i} - {\bf D}_l {\bf w}_i \|^2_2 + \lambda \| {\bf w}_i \|_1\tag{6}$

高解像画像パッチの生成

${\bf x}_i$

${\bf D}$

${\hat{\bf w}_i}$

$\displaystyle {\bf x}_i = {\bf D}{\hat {\bf w}_i} \tag{7}$

パッチをつなぎ合わせて高解像画像を再構成

${\bf X}_0$

$p$

ピクセル

$ps$

ピクセル

バックプロジェクションによる画像の補正

${\bf X}_0$

${\bar {\bf Y}}$

${\bf X}$

作用素

${\bf L}$

${\bar {\bf Y}} = {\bf L} {\bf X}$

${\bf X}_0$

最適化問題

$\displaystyle {\bf X}^* = {\rm arg} \min_{{\bf X}} \|{\bar {\bf Y}} - {\bf L} {\bf X}\|^2_2 + \mu \|{\bf X} - {\bf X_0} \|\tag{8}$

$\bf{X}_0, {\bar{\bf Y}}, {\bf X}$

${\bf L}$

逆行列

実装したコード

以上の内容の実装が下記です。モジュールのimportやおれおれ画像変換関数群は辞書学習時のコードと同様のものを使っています。

#バックプロジェクションで使うダウンサンプラー行列を求める関数
#ダウンサンプラーはスパースかつサイズが大きい行列なのでscipy.sparseのスパース行列型として返す
def get_downsampler_arr(from_H, from_W, to_H, to_W):
    S = sci.sparse.lil_matrix((to_H * to_W, from_H * from_W))
    k = 0
    for i in range(0, from_H, from_H/to_H):
        for j in range(0, from_W, from_W/to_W):
            s = sci.sparse.lil_matrix((from_H, from_W))
            s[i : i + from_H/to_H, j:j + from_W/to_W] = (float(to_W)/float(from_W)) * (float(to_H)/float(from_H))
            S[k, :] = s.reshape((1, from_H*from_W))
            k += 1
    return  S.tocsr()

#ベクトルデータのダウンサンプリング関数
def vec_downsampling(vec, from_H, from_W, to_H, to_W):
    arr = vec.reshape((from_H, from_W))
    ret_arr = arr_downsampling(arr, from_H, from_W, to_H, to_W)
    return ret_arr.reshape(-1)

#行列データのダウンサンプリング関数
def arr_downsampling(arr, from_H, from_W, to_H, to_W):
    ret_arr = np.zeros((to_H, to_W))
    k = 0
    for i in range(0, from_H, from_H/to_H):
        for j in range(0, from_W, from_W/to_W):
            ret_arr[i/(from_H/to_H), j/(from_W/to_W)] = (np.average(arr[i : i + from_H/to_H, j:j + from_W/to_W]))
            k += 1
    return ret_arr

#画像パッチをつなぎ合わせて画像を再構成する関数
#shift_numがずらし量、パッチがオーバーラップしている個所は足し合わせて平均化
def patch_vecs2gray_im(patchs, im_hight, im_width, patch_hight, patch_width, shift_num):
    im_arr = np.zeros((im_hight, im_width))
    sum_count_arr = np.zeros((im_hight, im_width)) #各領域何回パッチを足したか表す行列
    i = 0
    j = 0
    w_over_shift = 0
    h_over_shift = 0
    for k in range(patchs.shape[1]):
        P = patchs[:, k].reshape((patch_hight, patch_width))
        im_arr[(i* shift_num) + h_over_shift:((i* shift_num) + patch_hight) + h_over_shift,
            (j* shift_num) + w_over_shift:((j* shift_num) + patch_width) + w_over_shift] += P #指定の領域にパッチの足しこみ
        sum_count_arr[(i* shift_num) + h_over_shift :(i* shift_num) + patch_hight + h_over_shift,
            (j* shift_num) + w_over_shift:(j* shift_num) + patch_width + w_over_shift] += 1 #当該領域のパッチを足しこんだ回数をカウントアップ
        if j < ((im_width - patch_width)/shift_num + (1 if (im_width - patch_width) % shift_num  != 0 else 0)):
            j += 1
            w_over_shift = min(0, im_width - ((j* shift_num) + patch_width)) #パッチを足しこむ領域が画像からはみ出た時、そのはみ出た分を戻すための変数（横方向）
        else:
            j = 0
            i += 1
            h_over_shift = min(0, im_hight - ((i* shift_num) + patch_hight)) #パッチを足しこむ領域が画像からはみ出た時、そのはみ出た分を戻すための変数（縦方向）
            w_over_shift = 0
    im_arr /= sum_count_arr
    return float_arr2image(im_arr), im_arr

#バックプロジェクションを行う関数
def back_projection(x0, y, L, c, mu, x_init = [], eps = 1e-7):
    if x_init == []:
        xt = sci.matrix(np.random.rand(x0.shape[0], x0.shape[1])) *255.0
    else:
        xt = x_init
    result = -1
    while 1:
        x_new = xt - mu * (L.T.tocsr()*(L * xt - y) + c* (xt - x0))
        result_new = sci.linalg.norm(L * x_new - y, 2) + c * sci.linalg.norm(x_new - x0, 2)
        if np.abs(result_new - result) < eps:
            break;
        xt = x_new
        result = result_new
    return x_new, result_new

#超解像計算クラス
class super_resolution:
    def __init__(self, dictionary, weight_learner, back_projection):
        self.Dh = dictionary.get_dic()
        self.patch_hight = dictionary.get_patch_hight()
        self.patch_width = dictionary.get_patch_width()
        self.weight_learner = weight_learner
        self.back_projection_method = back_projection

    def get_super_resolution(self, image, scaling, shift_num = 1):
        #手順1:画像パッチの切り出し
        image = color2gray_img(image)
        Y = gray_im2patch_vec(image, self.patch_hight/scaling, self.patch_width/scaling, shift_num)

        #手順2:辞書のダウンサンプリング
        Dl = np.zeros((self.Dh.shape[0]/(scaling**2), self.Dh.shape[1]))
        for i in range(self.Dh.shape[1]):
            Dl[:, i] = vec_downsampling(self.Dh[:, i], self.patch_hight, self.patch_width,
                                        self.patch_hight/scaling, self.patch_width/scaling)
       
        im_arr = np.zeros(image.size)
        i = 0
        j = 0
        Patchs = np.zeros(((self.patch_hight) * (self.patch_width), Y.shape[1]))
        W = np.zeros((self.Dh.shape[1], Y.shape[1]))
        for k in range(Y.shape[1]):
            #手順3:符号化係数の推定
            [w, result] = self.weight_learner(Dl, Y[:, k])

            #手順4:高解像画像パッチの生成
            Patchs[:, k] = np.dot(self.Dh, w)

        #手順5:パッチをつなぎ合わせて高解像画像を再構成
        [reconstructed_im, X0] = patch_vecs2gray_im(Patchs, image.size[0] * scaling , image.size[1] * scaling , self.patch_hight, self.patch_width, shift_num* scaling)
        L = get_downsampler_arr(image.size[0] * scaling, image.size[1] * scaling, image.size[0], image.size[1])

        #手順6:バックプロジェクションによる画像の補正
        x = np.array(self.back_projection_method(sci.matrix(X0.reshape(-1)).T, sci.matrix(np.array(image).reshape(-1)).T, L)[0])
        return float_arr2image(x.reshape((image.size[0] * scaling, image.size[1] * scaling)))


def main():
    test_img = Image.open('Lenna.png')
    scale = 3
    lam = 50
    base_num = 144
    patch_hight = 12
    patch_width = 12
    shift_num = 1
    c = 0.8
    mu = 0.9
    weight_learning_method = lambda W, y:ISTA(W, y, lam)
    back_projection_method = lambda x0, y, L:back_projection(x0, y, L, c, mu)
    dic_learning_method = lambda W, Y:dictionary_learning_with_mod(W, Y)
    im_dic = image_patch_dictionary(dictionary_with_learning(weight_learning_method, dic_learning_method, base_num))
    im_dic.dictonary.Dic = np.load('dic_array.npy')
    im_dic.patch_hight = 12
    im_dic.patch_width = 12
    SR = super_resolution(im_dic, weight_learning_method, back_projection_method)
    result_img = SR.get_super_resolution(test_img, scale, shift_num)
    result_img.save("Lenna_SR.png")

超解像計算クラスsuper_resolutionは、画像パッチ辞書オブジェクト、係数推定関数オブジェクト、バックプロジェクション関数オブジェクトを生成時に指定し、get_super_resolutionメソッドで指定した辞書、関数に基づいた超解像画像の生成を行います。推定方法は上記で解説したとおりです。
画像パッチをつなぎ合わせて画像の再構成を行う関数patch_vecs2gray_imは、パッチのオーバーラップを許し、オーバーラップした部分は足し合わせた回数分割ることで平均をとっています。また、バックプロジェクションを行う関数back_projectionは、サイズが大きくスパースな行列Lを扱うので、scipy.matrixクラスで各引数を受け取る関数としています。ダウンサンプラー行列はget_downsampler_arr関数で求めています。

実験

まず辞書学習ですが、サイズ250×250の顔画像50枚から12×12のパッチ22000枚程度切り出して学習を行いました。また、辞書の基底数はパッチベクトルの次元と同じ144、学習時の $\lambda$ は $\lambda = 5$ にしています。なお、辞書の学習繰り返し回数は8で打ち切っています（それ以上繰り返した学習辞書も作りましたが、生成される超解像画像の精度が逆に悪くなるという結果に・・・）。
超解像画像生成時は、 $\lambda=50$ 、バックプロジェクション時のパラメータ $c = 0.8$ 、 $\mu = 0.9$ 、パッチのずらし量 $1$ で実験しています。実験に使用する画像データはおなじみのレナさんです。

83×83の画像を249×249の画像へ（拡大率3倍）

入力画像

f:id:YamagenSakam:20180411231300p:plain

左：通常の拡大　右：超解像アルゴリズムの出力

f:id:YamagenSakam:20180411232337p:plain 　 f:id:YamagenSakam:20180411231626p:plain

左の通常の拡大はPILのresizeで得られる画像です。少しわかりにくいですが、通常の拡大に比べて右の超解像画像は滑らかになっているのがわかると思います。

62×62の画像を248×248の画像へ（拡大率4倍）

入力画像

f:id:YamagenSakam:20180411231846p:plain

左：通常の拡大　右：超解像アルゴリズムの出力

f:id:YamagenSakam:20180411231342p:plain 　 f:id:YamagenSakam:20180411231925p:plain

もう少しわかりやすくするため、更に拡大率を上げた画像を生成してみました。通常の拡大に比べると滑らかになってはいますが、ブロック的なノイズが乗ってしまっています。アルゴリズム的な限界なのか、辞書やパラメータを洗練すればもっとうまくいくのか、その辺りの検証は今後の課題です。

まとめ

今回は近接勾配法の応用例として辞書学習および更にその応用である超解像画像生成をやってみました。
実際やってみると、パッチのつなぎ目が汚くなったりノイズまみれになったりと結構苦労しました。その過程で色々わかったこと・思ったことを下記にまとめます（その分野では当たり前のことばかりかも知れませんが・・・）。また、やり残したこと、今後やりたいことも併せて記します。

やってみてわかったこと、思ったこと
- 辞書学習時はパッチベクトルを平均 ${\bf 0}$ 、分散 $1$ へ標準化した方がよい。
- 画像再構成時はパッチ間の境界を滑らかにするため、元画像から切り出すパッチは領域を重複させて切り出し、再構成後のパッチを重ね合わせることで画像全体を再構成する方がよい。
- 学習させすぎると精度が落ちるので、途中で学習を打ち切った方が超解像の結果がよい。
- $\lambda$ によって結構結果が変わる。係数推定結果の非零要素数がパッチベクトルの要素数と大体同じになる程度が妥当か。
- バックプロジェクションはあまり効かない？（ダウンサンプリング以外にボケなどを考慮した ${\bf L}$ すると有効か？）

やり残したこと、今後やりたいこと
- 学習用画像データのバリエーションを増やして実験。
- 学習回数が増えると精度が落ちるのはなぜか数理的に考察。
- 推定した画像の一部の要素が0~255の範囲を超えることがある。制約条件をつけられないか検討。
- K-SVDなどの他の辞書学習手法の適用。
- 独立成分分析やウェーブレットを用いた辞書行列の適用。
- 他のダウンサンプリング辞書生成方法の適用。
- 辞書学習以外の超解像手法との比較。

参考文献

スパース表現の数理とその応用：スパースコーディング全般の解説。
スパースコーディングを用いたマルチフレーム超解像：超解像手法の論文。
近傍パッチの情報に基づく重みを利用したスパース再構成による一枚超解像：画像パッチ再構成時近傍のパッチも利用する手法を提案。
Image Super-Resolution via Sparse Representation：かなり有名な論文。今回使用した手法とは異なり、ダウンサンプリング辞書は辞書学習時に一緒に学習してる。
辞書学習のアルゴリズム：スライド資料。K-SVDなどの解説がわかりやすい。
scikit-learnでSparse CodingとDictionary Learning -実践編：scikit-learnを用いた辞書学習の実装が書かれている。

*1:実務等で利用するときはscikit-learn辺りのライブラリを使うべきですが、今回は理解を深めるためnumpy、scipy、PILを使った実装を行っています。

*2:実は ${\bf D}$ が正規分布に従ったランダムな行列でも、符号化としてはある程度うまくいくという研究もあります。ただ超解像やその他画像再構成への応用という観点ではどうなんでしょう？

2018-02-14

近接勾配法とproximal operator

機械学習連続最適化スパース学習近接勾配法

はじめに

前回、前々回とl1ノルム正則化項をつけた離散フーリエ変換により、スパースな解が得られること、そして基底ベクトルを並べた行列 ${\bf W}$ がユニタリ行列（直交行列）のため解析的に解けることを見てきました。それでは、 ${\bf W}$ が一般の行列の場合どうするか？今度こそ解析的に解けないので、反復法により最適解を求めます。しかし、目的関数にはl1ノルムという微分できない項があるため、最急降下法などのアルゴリズムは使えません。
このような微分できない項を含む最適化問題を効率的に解く手法の1つに近接勾配法（Proximal Gradient Method）があります。今回はこの近接勾配法と近接勾配法において重要な役割を果たすproximal operatorについて書いていきます。また、pythonで近接勾配法の実装を行いましたので、そのソースコードも載せておきます。

近接勾配法、proximal operatorとは

まずは、近接勾配法で解ける最適化問題の定義です。

$\displaystyle\min_{{\bf x}} f({\bf x} )+ g({\bf x}) \tag{1}$

ただし、 $f({\bf x})$ は微分可能な凸関数、 $g({\bf x})$ は微分可能とは限らない凸関数*1です。
l1ノルム正則化項付きの問題を式（1）に当てはめると、 $g({\bf x}) = \lambda \|{\bf x}\|_1$ となります。

次にproximal operatorの説明です。とりあえず定義を見てみましょう。関数 $g({\bf x})$ のproximal operatorは以下のように定義されます。

$\displaystyle {\rm prox}_{\gamma g}({\bf v}) = \arg \min_{{\bf x}} g({\bf x}) + \frac{1}{2 \gamma} \|{\bf x} - {\bf v} \|^2_2 \tag{2}$

ここで $\gamma$ は後述の更新式のステップ幅です。
このproximal operatorは何者か？というのは、 $g({\bf x})$ として以下のような指示関数を考えると見えてきます。

$\begin{eqnarray} g({\bf x}) = I_c(\bf{x}) &=& \left\{ \begin{array}{ll} 0 & ({\bf x} \in C) \\ \infty & ({\bf x} \notin C) \\ \end{array} \right. \end{eqnarray} \tag{3}$

ここで $C$ は閉凸集合です。この関数は ${\bf x}$ が集合 $C$ の要素なら0、 $C$ の要素でないなら $\infty$ に吹っ飛ばすという関数です。そして、この指示関数のproximal operatorは以下のようになります。

$\displaystyle \begin{eqnarray} {\rm prox}_{\gamma g}({\bf v}) &=& \arg \min_{{\bf x}} I_c(\bf{x}) + \frac{1}{2 \gamma} \|{\bf x} - {\bf v} \|^2_2 \\ &=& \arg \min_{{\bf x} \in C} \|{\bf x} - {\bf v} \|^2_2 \\ &=& \Pi_C({\bf v})\\ \end{eqnarray} \tag{4}$

${\bf v}$ が集合 $C$ の要素の場合、そのまま ${\bf x} = {\bf v}$ とすれば $\|{\bf x} - {\bf v} \|^2_2$ を最小化できます。一方で、 ${\bf v}$ が集合 $C$ の要素ではない場合、集合 $C$ の要素で ${\bf v}$ とのユークリッド距離を最小とするものが選ばれます。すなわち、 ${\bf v}$ の集合 $C$ への射影です（このような射影作用素を $\Pi_C({\bf v})$ として定義しています）。
よって、proximal operatorは射影作用素の関数 $g$ への一般化と考えることができ、 ${\bf v}$ が $g$ によって定められる領域からはみ出していたら、その領域に入るように射影する作用素と見なせます。
proximal operatorの具体例を挙げると、 $g({\bf x})$ が $g({\bf x}) = \lambda\| \bf{x}\|_1$ というl1ノルムであれば、前回述べたソフトしきい値作用素がproximal operatorになります。つまり、

$\begin{eqnarray} {\rm prox}_{\gamma g}({\bf v}) = S_{\gamma \lambda} ({\bf v}) &=& \left\{ \begin{array}{ll} v_i - \gamma \lambda & (v_i \geq \gamma \lambda) \\ 0 & (- \gamma \lambda < v_i < \gamma \lambda) \\ v_i + \gamma \lambda & (v_i \leq -\gamma \lambda) \\ \end{array} \right. \end{eqnarray} \tag{5}$

です。そして、近接勾配法の更新式はこのproximal operatorを使って以下のようになります。

$\displaystyle{\bf x}_{t+1} = {\rm prox}_{\gamma g}({\bf x}_t - \gamma \nabla f({\bf x}_t)) \tag{6}$

この更新式の収束性や正しさなどの説明は他の文献を参照いただくとして、直観的には次のような解釈ができると思います。
まず、 ${\bf x}_t - \gamma \nabla f({\bf x}_t)$ は最急降下法の更新式です。また、上述のようにproximal operatorは $g$ によって定められる領域に射影する作用素です。従って、最急降下法の1ステップで得られる値が $g$ の領域からはみ出してたら、その領域へ射影するという操作を繰り返して最適解を探すアルゴリズムと見ることができます。
このように、関数 $f({\bf x})$ の勾配と関数 $g({\bf x})$ のproximal operatorがわかれば、近接勾配法により最適解が得られます。
なお、 $\gamma$ の決め方も収束のためには重要ですが、こちらについても別の文献を参照してください。

Moreau decompositionと双対ノルム

さて、ここまで見たように近接勾配法を使うためにはproximal operatorを求める必要があります。しかし、関数 $g({\bf x})$ によっては、式（2）から直接的にproximal operatorを求めるのが難しいことがあります。例えば、 $g({\bf x})=\lambda \| {\bf x} \|_2$ とl2ノルムとした場合（2乗が付いていないことに注意）、式（2）を解こうとして ${\bf x}$ で微分しても、 $1/\| {\bf x} \|_2$ がかかる項が出てきてしまうので厄介です。
そこで便利なのが以下に示すMoreau decompositionが成り立つというproximal operatorの性質です。

$\displaystyle {\bf v} = {\rm prox}_{g}({\bf v}) +{\rm prox}_{g^*}({\bf v}) \tag{7}$

ここで、 ${\rm prox}_{g^*}({\bf v})$ は $g({\bf x})$ の共役関数 $g^*({\bf x})$ のproximal operatorです。この共役関数とは下記の式で表される関数のことを言います。

$\displaystyle g^*({\bf z}) = \sup_{{\bf x}} {\bf z}^T {\bf x} - g({\bf x}) \tag{8}$

よって、この共役関数のproximal operatorさえ分かれば、元の関数のproximal operatorも求められます。
では、 $g({\bf x}) = \lambda\| \bf{x}\|$ というノルム関数の場合どうなるかを見ていきます。まずノルム関数の共役関数は以下の指示関数となります。

$\begin{eqnarray} g^*({\bf z}) = I_{\it B} \left({\bf z} \right) &=& \left\{ \begin{array}{ll} 0 & (\|z \|_{*} \leq \lambda) \\ \infty & (\|z \|_{*} > \lambda) \\ \end{array} \right. \end{eqnarray}\tag{9}$

ここで、 $\| \bf{x}\|_{*}$ は双対ノルムと呼ばれるもので、以下の式で定義されます。

$\displaystyle\| {\bf z} \|_{*} = \sup_{{\bf x}} {\bf z}^T {\bf x} \quad ({\rm s.t.} \quad \| {\bf x} \| \leq 1) \tag{10}$

ということで、共役関数のproximal operatorは式（4）でみたように下記の射影作用素になります。

${\rm prox}_{g^*}({\bf v}) = \Pi_{\it B} ({\bf v}) \tag{11}$

ただし、 ${\it B}$ は ${\it B}=\left\{ {\bf x} \, | \, \| \bf{x} \|_* \leq \lambda \right\}$ で表される集合です。
つまり、 ${\bf v}$ を $\| {\bf v} \|_* \leq \lambda$ の範囲に射影する作用素がノルム関数の共役proximal operatorとなります。
おそらく式だけをザッと書いているので何故ノルムの共役関数が指示関数なのか？双対ノルムってそもそも何？どうやって求めるの？などの疑問が湧くかと思います。その辺りを少し補足しておきます。

双対ノルムについての補足

説明が前後しますが、後の説明のためにまずここで双対ノルムについて補足を書きます。そのための準備としてヘルダーの不等式という不等式を紹介します。
$1 \leq p,q \leq \infty$ という2つの整数に、 $1/p + 1/q = 1$ という関係があれば、

${\bf z}^T{\bf x} \leq \| {\bf x}\|_p \| {\bf z} \|_q \tag{12}$

が成り立つというのがヘルダーの不等式です。実は、 $\| \cdot \|_p$ の双対ノルムは $\| \cdot \|_q$ であり、その逆も成り立ちます。そのため、l2ノルムの双対ノルムはそのままl2ノルムですし、l1ノルムの双対ノルムはl $\infty$ ノルムです。
このことはヘルダーの不等式から示すことができます。まずヘルダーの不等式より、任意の $\|{\bf x}\|_p \leq 1$ に対して $\| {\bf z} \|_q \geq {\bf z}^T {\bf x}$ が成り立ちます。lpノルムの双対ノルム $\|{\bf z} \|_*$ は $\|{\bf x}\|_p \leq 1$ という条件下での ${\bf z}^T {\bf x}$ の最大値です。ここで ${\bf x}$ をうまくとることにより、等式を成立させることができます。つまり、

$\displaystyle \begin{eqnarray} \| {\bf z} \|_{*} &=& \sup_{{\bf x}} {\bf z}^T {\bf x} \quad ({\rm s.t.} \quad \| {\bf x} \|_p \leq 1) \\ & =& \| {\bf z} \|_q \\ \end{eqnarray} \tag{13}$

が成立します。

何故、ノルム関数の共役関数が指示関数なのか？

まず、 $\| {\bf z} \|_* > \lambda$ の時を考えます。式（8）の定義より $g^*({\bf z}) = \sup_{{\bf x}} {\bf z}^T {\bf x} - \lambda \| {\bf x} \|$ なわけですが、双対ノルムの定義式（10）と前述の前提より、 ${\bf z}^T {\bf x} > \lambda$ かつ $\|{\bf x}\| \leq 1$ となる ${\bf x}$ が存在することになります。これを満たす ${\bf x}$ を ${\hat{\bf x}}$ とした場合、その ${\hat{\bf x}}$ に定数 $t$ をかけた $t {\hat{\bf x}}$ も式（8）の最適解の候補として入ってきます。これはもちろん、 $t \rightarrow \infty$ とした場合も同様です。これらのことから、

$t{\bf z}^T {\hat {\bf x}} - \lambda \| t{\hat {\bf x}} \| = t({\bf z}^T {\hat {\bf x}} - \lambda \| {\hat {\bf x}} \| ) \rightarrow \infty \tag{14}$

となるため、 $\| {\bf z} \|_* > \lambda$ のときは $g^*({\bf z}) = \infty$ となります。
一方で、 $\| {\bf z} \|_* \leq \lambda$ の時は、式（12）のヘルダーの不等式より ${\bf z}^T{\bf x} \leq \| {\bf x}\| \| {\bf z} \|_*$ なので、

${\bf z}^T{\bf x} - \lambda \| {\bf x}\| \leq \| {\bf x}\| \| {\bf z} \|_* - \lambda \| {\bf x} \| = \| {\bf x}\| ( \| {\bf z} \|_* - \lambda) \tag{15}$

です。前提より、 $\| {\bf z} \|_* \leq \lambda$ なので、左辺（つまりは $g^*({\bf z})$ ）が0以上になることはありません。ここで、 ${\bf x}= {\bf 0}$ という取り方もできるので、式（15）の左辺を0にすることができます。よって、 $\| {\bf z} \|_* \leq \lambda$ のときは $g^*({\bf z}) = 0$ となります。

結局のところproimal operatorはどうやって求めればよいか？

2つ方法があると思います。1つ目は式（2）から直接求める方法です。例えば、前回記したようの方法でl1ノルム正則化項のproximal operator式（5）を求めるというものです。
もう1つの方法としては、式（7）のMoreau decompositionの性質を使う方法です。特に正則化項がノルムなら、共役関数のproximal operatorは ${\it B}=\left\{ {\bf x} \, | \, \| \bf{x} \|_* \leq \lambda \right\}$ という超球への射影なので簡単に求められます。例えば、l2ノルムの双対ノルムはl2ノルムなので、

$\displaystyle\begin{eqnarray} {\rm prox}_{g^{*}} ({\bf v})= \Pi_{\it B}({\bf v}) = \left\{ \begin{array}{ll} \lambda \frac{\bf{v}}{\| \bf {v} \|_2} \quad (\| {\bf v} \|_2 > \lambda) \\ {\bf v} \quad (\| {\bf v} \|_2 \leq \lambda) \\ \end{array} \right. \end{eqnarray} \tag{15}$

であり、求めたいproximal operatorは、

$\displaystyle\begin{eqnarray} {\rm prox}_{g} ({\bf v})= {\bf v} - \rm {prox}_{g^*} ({\bf v})= \left\{ \begin{array}{ll} {\bf v} -\lambda \frac{{\bf v}}{\| {\bf v} \|_2} \quad (\| {\bf v} \|_2 > \lambda) \\ 0 \quad (\| {\bf v} \|_2 \leq \lambda) \\ \end{array} \right. \end{eqnarray} \tag{16}$

となります。

なお、式（2）にはステップ幅の $\gamma$ が付いています。 $\gamma$ が付いていても基本的には上式の $\lambda$ を $\gamma \lambda$ に置き換えればよいだけですが、 $\gamma$ を含めたMoreau decompositionの一般形を書いておきましょう。

$\displaystyle {\bf v} = {\rm prox}_{\gamma g}({\bf v}) +\gamma{\rm prox}_{g^*/\gamma} \left(\frac{{\bf v}}{\gamma} \right) \tag{17}$

という形式になります。したがって、 $g$ がノルム関数の場合は、

$\displaystyle {\rm prox}_{\gamma g}( {\bf v}) = {\bf v} - \gamma{\rm prox}_{g^*/\gamma} \left(\frac{{\bf v}}{\gamma} \right) = {\bf v} - \gamma \Pi_{\it B} \left(\frac{{\bf v}}{\gamma} \right) \tag{18}$

になります。ちなみに、多くの文献では $g({\bf x })=\| {\bf x}\|$ としているので、 ${\it B}=\left\{ {\bf x} \, | \, \| \bf{x} \|_* \leq 1 \right\}$ として理論を展開している点に注意してください。

様々な正則化項のproximal operator

ここまでノルム関数が正則化項の時のproximal operator導出方法を見てきました。ここで、いくつかの正則化項に対するproximal operatorを記載しておきます。

l2ノルムの2乗正則化： $g({\bf v}) = \lambda \| {\bf v} \|^2_2$

リッジ回帰等でよく使われる正則化項です。普通にproximal operatorの式を微分して0とおけば求められます。2乗しているので、以下のように式（16）とは異なる結果が得られます。

$\displaystyle prox_{\gamma g} ({\bf v})= \frac{1}{1+2\gamma\lambda} {\bf v} \tag{19}$

Elastic Net： $g({\bf v}) = \lambda_1 \| {\bf v} \|_1 + \lambda_2 \| {\bf v} \|^2_2$

l1ノルム正則化とl2ノルム2乗正則化の混合です。l1ノルムは2つ変数が同程度の寄与度だったとしてもどちらかしか選ばれない（一方が0になってしまう）ケースが多いですが、Elastic Netではこの点を改善できるようです。このproximal operatorはl1ノルム正則化の時と同様、変数ごとに分離して求めることができます。

$\displaystyle prox_{\gamma g} ({\bf v})= \frac{1}{1+2\gamma\lambda_2} S_{\gamma \lambda_1} ({\bf v}) \tag{20}$

グループl1ノルム（グループに重複がない場合）： $\displaystyle g({\bf v}) = \lambda \sum_{g \in G} \| {\bf v}_g \|_p$

変数がいくつかのグループに分けられるとき、それぞれのグループのlpノルムを取って、更にその和を取るノルムです（ベクトル ${\bf v}$ のうち、グループ $g$ に所属する変数を並べたのが ${\bf v}_g$ ）。これを正則化項に用いるとグループごとにスパースになる、ならないが決まります。このグループl1ノルムを正則化項として使った回帰をgroup lassoと言います。
これはグループに重複がなく完全に分離ができるため、グループごとにlpノルムのproximal operatorを求めればよいです。p=2の時のグループ $g$ のproximal operatorは下記の通りです。（式（16）とほぼ同じです。）

$\displaystyle\begin{eqnarray} {\rm prox}_{\gamma g} ({\bf v}_g)= \left\{ \begin{array}{ll} {\bf v}_g -\gamma \lambda \frac{{\bf v}_g}{\| {\bf v}_g \|_2} \quad (\| {\bf v}_g \|_2 > \gamma \lambda) \\ 0 \quad (\| {\bf v}_g \|_2 \leq \gamma \lambda) \\ \end{array} \right. \end{eqnarray} \tag{21}$

ちなみに、グループに重複がある場合も提案されていますが、ここでは割愛します（まだ追いきれてない）。

トレースノルム： $g({\bf Z}) = \lambda \| {\bf Z} \|_{{\rm trace}} = {\rm Tr}\left(\sqrt{{\bf Z}^T {\bf Z}}\right)$

最後は行列関数に関する正則化項です。トレースノルムは $\sqrt{{\bf Z}^T {\bf Z}}$ のトレースを正則化項とするものなので、 ${\bf Z}$ の特異値を $\sigma({\bf Z})=\left[\sigma_1, \sigma_2, \cdots, \sigma_m \right]$ とすると、 $\| {\bf Z} \|_{{\rm trace}} = \sum_{i=1}^{m} \sigma_i$ です。そのため、このトレースノルムを正則化項として使うと特異値が0になりやすく、結果低ランクな解が得られるようになります。
それではproximal operatorを見ていきます。 ${\bf Z}$ の特異値分解を ${\bf Z} = {\bf U} {\rm diag}(\sigma ({\bf Z})) {\bf V}^T$ とすると、proximal operatorは以下のようになります。

$\displaystyle {\rm prox}_{\gamma g} ({\bf Z}) = {\bf U} {\rm diag}(S_{\gamma \lambda}(\sigma ({\bf Z}))) {\bf V}^T \tag{22}$

特異値分解して特異値に対してソフトしきい値作要素を施した後に再構築すればよいわけです。

近接勾配法を実装

ということで式（1）を解く近接勾配法を実装しました。

近接勾配法を行う関数

ということで、まずは近接勾配法そのものの実装です。

def proximal_gradient(grad_f, prox, gamma, objective, init_x, tol = 1e-9):
    x = init_x
    result = objective(x)
    while 1:
        x_new = prox(x - gamma * grad_f(x), gamma)
        result_new = objective(x_new)
        if (np.abs(result - result_new)/np.abs(result) < tol) == True :
            break;
        x = x_new
        result = result_new
    return x_new, result_new

引数について説明します。
grad_fは $\nabla f({\bf x})$ にあたります。つまり、grad_fはxを引数として、その勾配を返す関数です。
proxはproximal operatorです。proxも関数で、入力変数vとステップ幅のgammaの2引数をとり、proximal operatorの計算結果を返します。
objectiveは目的関数です。objectiveはxを引数として、目的関数の計算結果を返す関数です。
init_xは初期値、tolは収束の許容誤差です。
やっていることは、式（6）の更新式をそのまま素直に実装しているだけです。

l1ノルム正則化項のproximal operator

次にproximal operatorの実装です。

def prox_norm1(v, gamma, lam):
    return soft_thresh(v, lam * gamma)

def soft_thresh(b, lam):
    x_hat = np.zeros(b.shape[0])
    x_hat[b >= lam] = b[b >= lam] - lam
    x_hat[b <= -lam] = b[b <= -lam] + lam
    return x_hat

prox_norm1がl1ノルムのproximal operatorを計算する関数です。補助としてsoft_threshという関数を作っています。

lasso回帰問題を解く関数

以下は、 $\min_{{\bf x}}\| {\bf y} - {\bf W}{\bf x} \|_2^2 + \lambda \| {\bf x}\|_1$ というlasso回帰問題を解く関数です。なお、l1ノルム正則化項に対する近接勾配法はISTA（Iterative Shrinkage Thresholding Algorithm）と呼ばれるため、関数名をISTAとしています。

def ISTA(W, y, lam):
    objective = lambda x:np.sum(np.power(y - np.dot(W, x), 2))/2.0 + lam * np.sum(np.abs(x)) #(1)
    grad_f = lambda x: np.dot(W.T, np.dot(W, x) - y)  #(2)
    (u, l, v) = np.linalg.svd(W)   #(3)
    gamma = 1/max(l.real*l.real)   #(4)
    prox = lambda v, gamma:prox_norm1(v, gamma, lam)  #(5)
    x_init = randn(W.shape[1])  #(6)
    (x_hat, result) = proximal_gradient(grad_f, prox, gamma, objective, x_init, 1e-5)  #(7)
    return x_hat, result

(1)ではラムダ式を使ってobjectiveをxを引数としてとり、上記の目的関数を計算する関数として定義しています。
また、 $\nabla f({\bf x})={\bf W}^T ({\bf W} {\bf x} - {\bf y})$ であり、(2)はこの計算を行うgrad_fをxに引数とる関数として定義しています。
(3)、(4)は $\gamma$ を定めています。最大特異値の逆数より小さい値を $\gamma$ として定めると収束するようですが、その辺りはまだ理解ができておりません。とりあえず、文献の通りに実装します。
(5)がproximal operatorです。prox_norm1関数はv,gamma,lamの3引数を受け取りますが、proximal_gradient関数に渡すproxはvとgammaの2引数のみ取る関数なので、ラムダ式を使って再定義しています。
(6)はランダムに初期値を決めて、(7)でproximal_gradient関数にて問題を解きます。

proximal operatorを別なものにする

上記ISTA関数のproximal operatorとobjectiveの部分だけ変えれば他の正則化にも対応可能です。例としてグループl1ノルム正則化（group lasso）の場合を示します。
まずは、グループl1ノルムのproximal operatorの実装です。

def prox_group_norm(v, groups, gamma, lam):
    u = np.zeros(v.shape[0])
    for i in range(1, max(groups) + 1):
        u[groups == i] = block_soft_thresh(v[groups == i] , gamma * lam)
    return u

def block_soft_thresh(b, lam):
    return max(0, 1 - lam/np.linalg.norm(b)) * b

groupsは1～G（Gはグループ数）を要素に持つリストで、各変数がどのグループに属しているかを表します。後はグループごとにblock_soft_thresh関数にて、l2ノルムのproximal operatorを計算しています。
そして、上記のISTA関数内の(5)でproxとして、prox_group_normを指定すればgroup lassoになります。

    prox = lambda v, gamma:prox_group_norm(v, groups, gamma, lam) #(5)'

もし、Elastic Netを使いたかったらここにElastic Netのproximal operatorを持ってくればよいし、他のノルムの場合も同様です（本当はobjectiveも変更する必要がありますが上記では省略しています）。
また、 $f({\bf x})$ を別なものにしたかったら $\nabla f$ を別途定義して、grad_fとすればよいです。
行列関数の場合も同様で、grad_fとproximal operator (あと、objective)さえ定義してあげたら機能します。そのため、proximal operatorとしてトレースノルムを指定してあげると、低ランク解を得ることができます。

まとめ

今回は近接勾配法とproximal operatorについて見てきました。上記の実装については次の機会に応用事例で実験したいと思います（ちなみに、疑似データで実験したところそれっぽく動いていました）。
また、今回は取り上げませんでしたが、この界隈ではISTAの高速版FISTAや拡張ラグランジュ法、その改良であるAlternating Direction Method of Multipliers (ADMM) など様々なアルゴリズムが提案されており、非常に奥が深いです。

応用例（2018/6/15追記）

近接勾配法の応用例の記事をいくつか書いたのでリンクを貼っておきます。

参考文献

最後にこの記事を書くにあたって参考にした文献です（基本的に前回とほぼ同じです）。

スパースモデリングのための凸最適化—近接勾配法による高速アルゴリズム：今回の件を調べるきっかけとなった解説論文
Proximal Algorithms：共役関数、双対ノルム、Moreau decompositionについて非常に参考になった
Optimization with Sparsity-Inducing Penalties：同上
Convex Optimization：UCバークレーの講義資料？上２つの内容を補間するのに役にたった

また日本語の書籍で、proximal operator界隈が詳しく書かれているのは以下2冊かなと思います（他にもあったら教えてください）。

スパース性に基づく機械学習 (機械学習プロフェッショナルシリーズ)

作者: 冨岡亮太
出版社/メーカー: 講談社
発売日: 2016/12/23
メディア: Kindle版
この商品を含むブログを見る

機械学習のための連続最適化 (機械学習プロフェッショナルシリーズ)

作者: 金森敬文,鈴木大慈,竹内一郎,佐藤一誠
出版社/メーカー: 講談社
発売日: 2017/09/08
メディア: Kindle版
この商品を含むブログを見る

*1:もうひとつ、「単純」という条件も付きます。

2018-02-04

l1ノルム正則化フーリエ変換を複素数のまま解く

機械学習スパース学習

はじめに

前回の記事では、離散フーリエ変換を係数をLasso回帰として求めスパースな解を得る方法を書きましたが、この際に複素数を取り扱いたくないため実数のみの問題に変換して、ソルバに解かせるということを行いました。今回は、複素数のまま問題を解くことを考え、複素数の最適解を得る方針で考えます。この方針で考察することで、以下のことが見えてきたので記事に残そうと思った次第です。

前回の記事では「この問題は解析的には解けず」などと書いているが、実は ${\bf W}$ が離散フーリエの基底ベクトルを並べた行列であれば解析的に解けること
結局は、離散フーリエ変換してしきい値未満の要素を0にしていることと（ほぼ）等価であること
複素数としてのl1ノルムを考えると、別な解が得られること

ちなみに、自分は「この方法だとなくうまくいきそうだな」と思ったことはまず試してみて、後から理屈を考える（特にうまくいかなかったとき）というスタンスなので、理論的、数学的な厳密性はあまり考慮していません。したがって、本来理論的にはあり得ない式展開、考え方をしている場合もあるかもしれませんがご了承ください。

定式化

前回と同様、正則化項をつけた回帰問題として定式化します。

$\displaystyle \min_{\bf x} \| {\bf y} - {\bf W}{\bf x}\|^2_2 + \lambda \Omega({\bf x}) \tag{1}$

${\bf y}$ は元の時間軸の信号で、 ${\bf y} = \left[ y(0), y(1), \cdots, y(T - 1) \right]^T$ 、 ${\bf x}$ が求めたい離散フーリエ係数であり、 ${\bf x} = \left[ x(0), x(1), \cdots, x(F - 1) \right]^T$ です。ここで ${\bf x}$ の各要素は複素数です。 ${\bf y}$ は実信号なので実数ベクトルですが、虚部がすべて0の複素数と見なしてもよいでしょう。また、行列 ${\bf W}$ の各列ベクトル ${\bf w}_k ( k = 0, 1, \cdots, T-1)$ は
${\bf w}_k = \frac{1}{T}\left[\exp(i\frac{2\pi k \times (0)}{T}) ,\exp(i\frac{2\pi k \times (1)}{T}), \cdots, \exp(i\frac{2\pi k \times (T - 1)}{T} ) \right]^T$
とします。前回の定義とほぼ一緒ですが、各要素を要素数 $T$ で割っています。こうすると、 ${\bf W}$ はユニタリ行列になります。つまり、 ${\bf W}^* {\bf W} = {\bar {\bf W}}^T {\bf W} = {\bf I}$ です( ${\bf W}^*$ は随伴行列、 ${\bar {\bf W}}$ は複素共役を表します)。
後は正則化項 $\lambda \Omega({\bf x})$ です。普通に考えると
$\Omega({\bf x}) = \| {\bf x} \|_1 = | x(0)| + |x(1)|+ \cdots + |x(T-1)| \tag{2}$

（ただし、 $|x|$ は $x$ の絶対値を表し $|x| = \sqrt{{\bar x} x}$ ）としたいところですが、実はこの定義だと前回と同じ結果にはなりません。このことは後ほど説明するとして、今は
$\Omega({\bf x}) = | {\rm Re}( x(0))| + | {\rm Im}( x(0))| + \cdots + | {\rm Re}( x(T - 1))| + | {\rm Im}( x(T - 1))| \tag{3}$

とします（ ${\rm Re(x)}$ は $x$ の実部、 ${\rm Im(x)}$ は $x$ の虚部を表す）。式（2）と式（3）の大きな違いは、実部と虚部が関連し合っているか、完全に独立なものとして扱うか、です。そのため、式（3）を正則化項として用いることは、実部と虚部を分けて解を求める前回の考え方と一致します。

解析解の導出と実装

当初は「スパースモデリングのための凸最適化」の解説論文を読んで近接勾配法で解こうなんて考えていましたが、解説を読んでいくと実は、 ${\bf W}$ がユニタリ行列であるため解析的に解が得られることがわかりました。
ということでその導出です。まず、解説に従うと、式（1）は以下のように式を変形できます。
$\begin{eqnarray} \displaystyle (1) &=& \min_{\bf x} \| {\bf y} - {\bf W}{\bf x}\|^2_2 + \lambda \Omega({\bf x}) \\ \displaystyle &=& \min_{\bf x} \| {\bf W}^*{\bf y} - {\bf W}^*{\bf W}{\bf x}\|^2_2 + \lambda \Omega({\bf x}) \\ \displaystyle &=& \min_{\bf x} \| {\bf b} - {\bf x}\|^2_2 + \lambda \Omega({\bf x}) \ \ \ (\because {\bf b} = {\bf W}^*{\bf y})\\ \displaystyle &=& \min_{\bf x} \sum_{i = 0}^{T - 1} \left\{ |(b(i) - x(i))|^2+ \lambda( | {\rm Re}( x(i))| + | {\rm Im}( x(i))| ) \right\} \\ \displaystyle &=& \min_{\bf x} \sum_{i = 0}^{T - 1} \left[ \left\{ {\rm Re} (b(i)) - {\rm Re} ( x(i) ) \right\} ^2 + \lambda | {\rm Re}( x(i))| + \left\{ {\rm Im} (b(i)) - {\rm Im} ( x(i) ) \right\} ^2 + \lambda | {\rm Im}( x(i))| \right] \\ \displaystyle &=&\sum_{i = 0}^{T - 1} \min_{x(i)} \left[ \left\{ {\rm Re} (b(i)) - {\rm Re} ( x(i) ) \right\} ^2 + \lambda | {\rm Re}( x(i))| + \left\{ {\rm Im} (b(i)) - {\rm Im} ( x(i) ) \right\} ^2 + \lambda | {\rm Im}( x(i))| \right] \end{eqnarray}$
このように式（1）の問題は各要素独立に最適解を選ぶ問題に変形されます。しかも、実部と虚部も分離しているので、こちらについても独立して最適解を選ぶことができます。ようするに、以下の最適化問題を各要素の実部、虚部それぞれ解けばよいわけです。
$\displaystyle\min_{x} (b- x)^2 + \lambda | x| \tag{4}$
この解は次のような式で得られます。
$\begin{eqnarray} {\hat x} = S_{\lambda} (b) &=& \left\{ \begin{array}{ll} b - \frac{\lambda}{2} & (b \geq \frac{\lambda}{2}) \\ 0 & (- \frac{\lambda}{2} < b < \frac{\lambda}{2}) \\ b + \frac{\lambda}{2} & (b \leq - \frac{\lambda}{2}) \\ \end{array} \right. \end{eqnarray} \tag{5}$
解説論文と異なり式（1）および式（4）の第1項に $1/2$ がかかっていないため、最適解も $\lambda$ に $1/2$ かけていることに注意してください。ここで式（5）の $S_{\lambda}(b)$ はソフトしきい値作用素と呼ばれる作用素です。
このソフトしきい値作用素を用いると、今回得たい式（1）の最適解は下記のように表せます。
${\hat{\bf x}} = S_{\lambda}({\rm Re}({\bf W}^* {\bf y}))+ i S_{\lambda}({\rm Im}({\bf W}^* {\bf y})) \tag{6}$

さて、ここで通常の離散フーリエ変換の係数は ${\bf W}^* {\bf y}$ で得られることを思い出すと、式（6）は通常の離散フーリエ変換の結果の実部、虚部それぞれに対してソフトしきい値作用素を作用させていることになります。
これは式（5）により離散フーリエ係数の絶対値が $\lambda/2$ 未満であれば0にすることに対応します。l1ノルム正則化項の離散フーリエ変換なんて言っていますが、紐解いていくと離散フーリエの結果に対ししきい値未満ならば0にするという割りとよくやる操作とほぼ同じになるわけですね（しきい値を超えている要素には $\pm \lambda/2$ を加える点は異なりますが）。
ということで、以下は上式によるスパース離散フーリエ変換を行う関数の実装です。

def sparse_dft_with_thresh(y, T, lam):
    Z = np.meshgrid(np.arange(T),np.zeros((T,1)))[0]
    u = np.arange(T)
    V = (Z.T * u) * (2 * np.pi /T)
    W= np.exp( 1j * V)/T
    return (soft_thresh((np.dot(W.T.conj(), y)).real, lam/2)
            +1j * soft_thresh((np.dot(W.T.conj(), y)).imag, lam/2))

def soft_thresh(b, lam):
    x_hat = np.zeros(b.shape[0])
    x_hat[b >= lam] = b[b >= lam] - lam
    x_hat[b <= -lam] = b[b <= -lam] + lam
    return x_hat

sparse_dft_with_threshがスパースな離散フーリエ変換を行う関数で、soft_threshがソフトしきい値作用素を施す関数です。soft_threshを解説と合わせるために、引数で渡す時にlam/2をしています。

複素数l1ノルム正則化

ここで1つ疑問が残ります。正則化として式（2）のノルムを使うとどうなるか？です。これまでの議論から、式（2）のノルムを使った場合、式（1）は次のように変形できます。
$\displaystyle \sum_{i = 0}^{T - 1} \min_{x(i)} \left\{ |(b(i) - x(i))|^2+ \lambda| ( x(i)| \right\} \tag{7}$
この式を見ると要素同士は独立して最適解を選ぶことができますが、実部と虚部は独立させることができません（ $|x(i)| = \sqrt{{\bar x(i)} x(i)}$ なので）。
この最適解を得るために、まず以下のような式を考えます。
$\displaystyle {\rm prox}_{\lambda}({\bf b})= \min_{\bf x} \frac{1}{2} \| {\bf b} - {\bf x} \|_2^2+ \lambda \Omega({\bf x}) \tag{8}$
この式（8）は近接作用素（proximal operator）と呼ばれるもので、近接勾配法などで重要な役割を持つ作用素です。今、式（7）について複素数を実部と虚部を並べた2次元ベクトルとして考えると（つまり、 ${\bf x} = \left[{\rm Re}(x) , {\rm Im}(x)\right]^T$ ）、式（7）と式（8）の問題は（第1項の $1/2$ を除いて）等価になります（ $\Omega({\bf x}) = \| x \|_2$ ）。
ということで、式（8）を解くことを考えますが、実はこれはgroup lassoにおけるproximal operatorを求めることと等価です。group lassoのproximal operatorの導出やその基礎となる劣微分、凸解析、双対理論などは「ptimization with sparsity-inducing penalties」や「Proximal Algorithms」にまとまっています。正直、まだ理解できていない部分が多いですが、これらの資料に記載されているgroup lassoにおけるproximal operatorを使うと式（7）の解は下記のようになります（たぶん）。
$\begin{eqnarray} {\hat x} = S_{\lambda} (b) &=& \left\{ \begin{array}{ll} b - \frac{\lambda b}{2|b|}& (|b| > \frac{\lambda}{2}) \\ 0 & (|b| \leq - \frac{\lambda}{2}) \\ \end{array} \right. \end{eqnarray} \tag{9}$

式（6）では、実部、虚部それぞれでしきい値未満だと0になりましたが、こちらは複素数の絶対値がしきい値未満だと0になります。つまり、実部と虚部は運命共同体で、しきい値未満となる要素は実部も虚部一緒に0になり、その逆もまたしかりです。ということで実装です。

def group_lasso_prox(b, lam):
    x_hat = np.zeros(b.shape[0], "complex")
    x_hat[np.abs(b) >= lam] = (b[np.abs(b) >= lam] -
                                lam *  b[np.abs(b) >= lam]/np.abs(b)[np.abs(b) >= lam])
    return x_hat

def sparse_dft_with_complex_l1_norm(y, T, lam):
    Z = np.meshgrid(np.arange(T),np.zeros((T,1)))[0]
    u = np.arange(T)
    V = (Z.T * u) * (2 * np.pi /T)
    W= np.exp( 1j * V)/T
    return (group_lasso_prox(np.dot(W.T.conj(), y), lam/2))

一応、実験した結果も載せておきます。実験には前回同様sin波とcos波とホワイトノイズを重ね合わせた疑似データを用います。

def main():
    T = 2**10
    y = (3 * np.cos(28* 2*np.pi * np.arange(T)/T) +2 * np.sin(28* 2*np.pi * np.arange(T)/T) +
            np.cos(400* 2*np.pi * np.arange(T)/T) + np.sin(125 * 2*np.pi * np.arange(T)/T) + 5 * randn(T))

    x1 = sparse_dft_with_thresh(y, T, 0.5)
    plt.figure()
    plt.xlim((0, T))
    plt.xticks(np.arange(0, T , 100))
    plt.ylim((-2, 2))
    plt.plot(x1.real, "r")
    plt.hold(True)
    plt.plot(x1.imag, "b")

    x2 = sparse_dft_with_complex_l1_norm(y, T, 0.5)
    plt.figure()
    plt.xlim((0, T))
    plt.xticks(np.arange(0, T , 100))
    plt.ylim((-2, 2))
    plt.plot(x2.real, "r")
    plt.hold(True)
    plt.plot(x2.imag, "b")

    plt.show()

実部と虚部を独立の変数として扱った正則化の結果
f:id:YamagenSakam:20180204165656p:plain

複素数l1ノルムによる正則化の結果
f:id:YamagenSakam:20180204165749p:plain

こうやって見ると、前者の方がスパース性という意味だと高いです。実部、虚部を運命共同体にする理由がない限り、前者のように独立して扱った方がよいように思えます。

最後に

ここまで見て感じたかと思いますが、結局のところ単なるしきい値での枝刈りに帰着するので、実装まで落としこんでしまうと正直あまり面白いものではありませんし、応用としてもあまり意味のあるものでもありません（これやるくらいなら、普通にFFTしてしきい値未満を0にする方が断然早いです）。
ただ、ここまで導出する過程で、ソフトしきい値作用素やproximal operator、劣微分、凸解析などなど最適化理論の色々なトピックスに触れることができました（まだ全然理解できていないですが）。
冒頭に書いた思いついたことはまず試してみるというスタンスは、試すために知らないトピックスを調べるきっかけになるというメリットがあります。やってみてダメなら何がダメか更に深追いするモチベーションにもなります。後はこういう試行錯誤は何より楽しいですね。ということで、今後は近接勾配法や凸解析の理論などの知識を深めていければと思います。

2018-01-29

l1ノルム正則化でスパースな離散フーリエ変換

機械学習スパース学習

はじめに

離散フーリエ変換は離散化された時間軸の信号を周波数軸に写像する基底変換の一種で、周波数解析などによく使われます。この離散フーリエ変換の係数は、基底ベクトルの重み付き線形和で表したときの誤差の2乗和を最小にする重みと見なせます。つまり、離散フーリエ係数は一種の線形回帰問題の解として得られます。ということは、l1ノルム正則化項つけたLasso回帰として解けばスパースな離散フーリエ変換ができるのでは？ということでやってみました。

離散フーリエ変換の線形回帰的定式化

元の信号のデータ点数を $T$ 個、基底ベクトルの個数を $F$ 個とし、上述の通り回帰問題として離散フーリエ係数を求める問題を定式化すると下記のようになります。

$\displaystyle \min_{\bf x} \| {\bf y} - {\bf W}{\bf x}\|^2_2 \tag{1}$

ここで、 ${\bf y}$ は元の時間軸の信号で、 ${\bf y} = \left[ y(0), y(1), \cdots, y(T - 1) \right]^T$ です。
${\bf W} = \left[{\bf w}_0, {\bf w}_1, \cdots, {\bf w}_{F - 1} \right]$ は $F$ 個の基底ベクトルを劣ベクトルとして並べた $T \times F$ の行列です。
一般的に離散フーリエ変換では ${\bf W}$ は ${\bf w}_k = \left[\exp(i\frac{2\pi k \times (0)}{T}) ,\exp(i\frac{2\pi k \times (1)}{T}), \cdots, \exp(i\frac{2\pi k \times (T - 1)}{T} ) \right]^T$ という複素数の要素を持つ行列で、各列ベクトル同士は直交します( ${{\bf w}_i^{*} {\bf w}_j} = 0 \quad i \neq j$ )。
そして、 $\bf{x}$ が求めたい離散フーリエ係数であり、 ${\bf x} = \left[ x(0), x(1), \cdots, x(F - 1) \right]^T$ です。
通常、 $F = T$ です。そのため、(1)の最適解は ${\hat{\bf x}} = {\bf W}^{-1} {\bf y} = {\bf W}^{*} {\bf y}$ で表され、最小2乗誤差は0になります。

l1ノルムによるスパース化

ということで、回帰問題としてみなせるならl1ノルムの正則化項をつける（つまりは、Lasso回帰を行う）とスパースな解が得られるだろうという考えです。(1)の目的関数に正則化項を付け加えます。

$\displaystyle \min_{\bf x} \| {\bf y} - {\bf W}{\bf x}\|^2_2 + \lambda \| {\bf x} \|_1 \tag{2}$

こうなると最早この問題は解析的には解けず、正則化項は $\|{\bf x} \|_1 \leq \gamma$ という制約条件を付けるのと等価なので、最適解による最小2乗誤差も一般的には0にはなりません。ということで、cvxpyというソルバを用いて上記の最適化問題を解きます。
なお、この定式化は入口こそ離散フーリエ変換を回帰問題と見なしてLasso回帰を行うという考えですが、結局のところこれは、この資料等に記載されいるノイズあり圧縮センシングの一種です。　

複素数の最適化問題から実数のみの最適化問題へ

離散フーリエ変換は上述のように、基底ベクトルが複素数で得られる解も一般的に複素数となります。ところが、このままだと実装上面倒くさい（というより、ソルバが動いてくれるかわからない）ので実数のみの式に変換します。といっても、 $\exp({i \theta}) = \cos\theta + i\sin\theta$ より実部と虚部に分けてそれぞれ個別の変数として考えるだけです。
具体的には、 ${\bf W}$ を ${\bf W} = \left[ {\bf w}_0^c, {\bf w}_1^c, \cdots, {\bf w}_{T - 1}^c, {\bf w}_0^s, {\bf w}_1^s \cdots, {\bf w}_{T - 1}^s \right]$ という $T \times 2T$ の行列として再定義します。ここで、
${\bf w}_k^c = \left[\cos(\frac{2\pi k \times (0)}{T}) ,\cos(\frac{2\pi k \times (1)}{T}), \cdots, \cos(\frac{2\pi k \times (T - 1)}{T} )\right]^T$
${\bf w}_k^s = \left[\sin(\frac{2\pi k \times (0)}{T}) ,\sin(\frac{2\pi k \times (1)}{T}), \cdots, \sin(\frac{2\pi k \times (T - 1)}{T} )\right]^T$
です。あわせて説明変数 ${\bf x}$ を要素数 $2T$ のベクトルとして式（２）を解けばよいです。この変換で基底ベクトルの数も説明変数の要素数も倍になったように見えますが、元々も複素数でそれぞれが実部と虚部を持っていたため本質的には変わりません。
また、離散フーリエ変換の場合、 $x(k) + {\bar x(k + \frac{T}{2}})$ という関係があるため、実質 $0\sim\frac{T}{2}-1$ 番目の要素までを求めれば、残りの $\frac{T}{2}\sim T-1$ 番目は自動で決定されます。これは、l1ノルム正則化項をつけても変わりません（実際に試しました）。そのため、説明変数の要素数を $2T\rightarrow T$ へとシュリンクすることも可能ですが、今回は理解を深めるためにもそのまま解くようにしました。

cvxpyを用いた実装

上述の通り、今回は最適化問題を解くためにcvxpyというツールを利用しました。このcvxpyの使い方はこちらの記事を参考にしました。cvxpyの使い方は他の解説記事に譲るとして、早速式（２）を解く関数のソースコードを示します。

from cvxpy import *
import numpy as np
from numpy.random import *
import matplotlib.pyplot as plt

#y:原信号、T:信号のデータ点数、lam:ハイパーパラメータλ
def sparse_dft(y, T, lam):
    #行列Wを生成
    V = np.meshgrid(np.arange(T),np.zeros((T,1)))[0]
    u = np.arange(T)
    V = (V.T * u) * (2 * np.pi /T)
    W = np.c_[np.cos(V), np.sin(V)]

    #変数xを定義
    x = Variable(T * 2)

    #ハイパーパラメータを定義
    lamda = Parameter(sign="positive")
    lamda.value =lam

    #目的関数を定義
    objective = Minimize(sum_squares((y - W * x))/T + lamda*norm(x, 1))
    p = Problem(objective)

    #最適化計算
    result = p.solve()

    return result, x

特に特別なことはしておらず、最初に ${\bf W}$ を生成して、後はcvxpyの解説に従って変数、ハイパーパラメータ、目的関数を定義し、p.solve()で問題を解きます。

実験

以下のようなコードで疑似データを発生させて実験してみました。

    T = 2**10
    y = (3 * np.cos(28* 2*np.pi * np.arange(T)/T) +2 * np.sin(28* 2*np.pi * np.arange(T)/T) +
            np.cos(400* 2*np.pi * np.arange(T)/T) + np.sin(125 * 2*np.pi * np.arange(T)/T) + 5 * randn(T) )

いくつかの周波数が異なるsin波とcos波を重み付きで重ね合わせて、更に正規分布に従うノイズを乗せた1024点のデータです。疑似データの波形は以下の通り。
f:id:YamagenSakam:20180129000732p:plain

通常の離散フーリエ変換

まずは上記の疑似データにFFTにより離散フーリエ変換を行いました。

    Y = np.fft.fft(y)/T #スケーリングを合わすため要素数で割る
    plt.figure()
    plt.xlim((0, T))
    plt.xticks(np.arange(0, T , 100))
    plt.ylim((-1.5, 1.5))
    plt.plot(Y.real , "r")
    plt.hold(True)
    plt.plot(Y.imag , "b")

結果は下記の図です。
f:id:YamagenSakam:20180129000744p:plain
重ねたsin波、cos波に対応する周波数成分が大きな値になるのは当たり前ですが、ノイズが重畳しているせいで全周波数成分まんべんなく値が入ります。

Lassoによる離散フーリエ変換

次は今回記した手法での離散フーリエ変換（ $\lambda = 0.5$ で実験）。上述のsparse_dft関数で計算します。

    (result, x) = sparse_dft(y, T, 0.5)

    plt.figure()
    plt.xlim((0, T))
    plt.xticks(np.arange(0, T , 100))
    plt.ylim((-1.5, 1.5))
    plt.plot(x.value[0:T,0], "r")
    plt.hold(True)
    plt.plot(x.value[T:T*2,0], "b")

f:id:YamagenSakam:20180129000758p:plain
通常の離散フーリエ変換とは異なりスパースな解が得られ、余分な周波数成分は0になっていることがわかります。

逆変換

得られたフーリエ係数から逆変換を行い元の信号を再構築してみました。こうすることで、余計なノイズ成分が除去された信号が得られるはずです。
なお、逆変換のため、得られた ${\bf x}$ から $\cos$ にかかる係数を実部、 $\sin$ にかかる係数を虚部として複素数に変換してから、ifftメソッドで逆変換を行っています。

    X=(x.value[0:T] + (x.value[T:T*2] * 1j))*T
    y_hat = np.fft.ifft(X.T)

    plt.figure()
    plt.xlim((0, T))
    plt.xticks(np.arange(0, T , 100))
    plt.ylim((-20, 20))
    plt.plot(y_hat.real.T, "r")

以下が原波形から今回の方法でノイズを除去したデータです。
f:id:YamagenSakam:20180129000809p:plain
また、以下が上述の疑似データ生成のコードからノイズ項を取り除いたデータの波形です（つまり、真に得たい波形）。
f:id:YamagenSakam:20180129002851p:plain
確かにノイズが除去された波形が得られています。実際のノイズ項を除いた信号と比べると若干スケーリングが小さくなっていますが、悪くない結果だと思います。

まとめ

離散フーリエ変換を回帰問題と見なし、目的関数にl1ノルムの正則化項をつけて誤差の最小化を行うことで、スパースな離散フーリエ係数が得られました。またそれを逆変換することで原信号からノイズを除去した信号を得ることもできました。また、最適化問題を解くためにcvxpyというツールを利用しました。
今回は何も考えずcvxpyに問題を突っ込んで解いただけでしたが、今後は、Lasso回帰の最適化等で使われる近接勾配法やcvxpyの詳しい使い方、凸解析など最適化の理論、実装両面で理解を深めていきたいと思います。

2018-01-20

はてなダイアリーからはてなブログへ移行

このたび、本ブログをはてなダイアリーからはてなブログへ移行しました。旧はてなダイアリーの記事にアクセスすると、対応する本ブログの記事へのリダイレクトされるようなので、ブックマーク登録等の変更は不要です。
さて、この移行作業に当たって、元のダイアリーで書いていた数式がうまく表示されず修正に手間取った箇所がいくつかあったので、備忘録の意味でもここに記しておきます。（大体は元の書き方がよろしくなかっただけですが）

小中大括弧でエラー

元々の書き方

[tex:\(x + y \)]
[tex:\{x + y \}]
[tex:\left{x + y \right}]
[tex:\[x + y \]]
[tex:\left[x + y \right]]

見え方
$x + y$ ：おかしい。
$\{x + y \}$ ：うまく表示される。
$\left{x + y \right}$ ：おかしい。
$[x + y ]$ ：うまく表示される。
$\left[x + y \right$ ]：おかしい。

小括弧は\エスケープだけではエラーになり、中大括弧は\left, \rightをつけるとエラーになります。そのため、小括弧には\left、\rightをつける、中大括弧は\left\{ \right}\}のように\left、\rightの後、更に\でエスケープする方法で対策しました。

修正後の書き方：

[tex:\left(x + y \right)]
[tex:\{x + y \}]
[tex:\left\{x + y \right\}]
[tex:\[x + y \]]
[tex:\left\[x + y \right\]]

見え方
$\left(x + y \right)$ ：\left, \right追加でOKになった。
$\{x + y \}$ ：変更なし。
$\left\{x + y \right\}$ ：\エスケープ追加でOKになった。
$[x + y ]$ ：変更なし。
$\left[x + y \right]$ ：\エスケープ追加でOKになった。

sum、min、max等の上置き,下置きが上付き、下付きになる

元々の書き方

[tex:\sum_{i = 1}^{N}x_i]
[tex:\min_{x} f(x)]

見え方
$\sum_{i = 1}^{N}x_i$
$\min_{x} f(x)$

本当は $\Sigma$ や $\min$ の真上、真下に添え字が来るようにしたいのに、右上、右下に来ています。これは\displaystyleにすることで解決しました。

修正後の書き方

[tex:\displaystyle \sum_{i = 1}^{N}x_i]
[tex:\displaystyle \min_{x} f(x)]

見え方
$\displaystyle \sum_{i = 1}^{N}x_i$
$\displaystyle \min_{x} f(x)$

\bfにより意図しないところまでも太字化される

元々の書き方

[tex:(\bf{x}_N)]

見え方
$(\bf{x}_N)$ ：xの先にある添え字のNや)まで太字化される。

わかりにくいですが、xだけ太字にしたいのに、それ以降の文字まで太字になっています。これは元々の書き方がそもそも間違っており、中括弧の位置を\bf{x}→{\bf x}と変更すると修正できます。

修正後の書き方

[tex:({\bf x}_N)]

見え方
$({\bf x}_N)$ ：xのみが太字。

ということで、今回は以上。せっかくはてなブログに移行しましたが、次回更新は遠い未来になるかも知れません。

2018-01-03

NIPS2017論文メモ

機械学習

お久しぶりです。2014年9月以来、実に3年ぶりの更新です。しばらく機械学習関連の勉強は滞っていたのですが、思うところがあって再開。とりあえず、昨年12月に開催されたNIPS2017のプロシーディングス中から面白そうな記事をピックアップして斜め読みしました。備忘録＆リハビリのためにも、記事を更新してみよう思った次第です。

Deep Subspace Clustering Networks

部分空間クラスタリング+オートエンコーダで非線形なクラスタリングを行う手法を提案。

部分空間クラスタリングは、同じクラスなら同じ（低次元）部分空間埋め込まれている（つまり、同じクラスに所属している別のサンプルの線形結合で表せる）という仮定の下クラスタリングを行う枠組みで、SSCやLRRが有名。サンプル集合を $\bf{X}$ として、いずれの手法も $\bf{X}=\bf{X}\bf{C}$ という制約条件のもと、線形結合の行列 $\bf{C}$ のスパース化や低ランク化を行う（前者は $l1$ ノルムの最小化、後者はトレースノルムの最小化）。 $\bf{C}$ を求めたら $|\bf{C}|+|\bf{C}^T|$ をaffinity matrixとしてスペクトラルクラスタリングを行うのが一般的（他にも色々やり方があるらしい）。

これらの方法では非線形なクラスタリングは行えない。そこで、オートエンコーダを利用しようというのがこの論文の提案。オートエンコーダはニューラルネットを入力と出力を同じデータとして学習させることで、中間層の出力 $\bf{z}$ が非線形な特徴抽出となるという手法。入力から特徴を抽出する中間層まで層をエンコーダ層、特徴を抽出する中間層から出力までの層をデコーダ層と呼ぶ。この論文では、このエンコーダ層とデコーダ層の間に $\bf{Z}\bf{C}$ とするSelf-Expressive層を新たに設けることを提案( $\bf{Z}$ はエンコーダ層で特徴抽出されたサンプルの集合)。これは、線形の全結合となる。この $\bf{C}$ と他の層の重みを求めるため、下記の損失関数を最小化。
$L(\Theta, {\bf C} ) = \frac{1}{2} \|{\bf X} - {\bf \hat{X}_{\Theta}} \|^2_{F}+\lambda_1 \|{\bf C} \|_p + \frac{\lambda_2}{2} \|{\bf Z_{\Theta_e}} - {\bf Z_{\Theta_e}} {\bf C} \|^2_{F}$ $\quad {\rm s.t.} \quad {\rm diag}( {\bf C} )={\bf 0}$
$\Theta$ はニューラルネットの重み係数集合、 ${\bf Z_{\Theta_e}}$ はエンコーダ層の出力集合（バッチ）を表す。ノルム項は $l1$ と $l2$ の両方を実験。

上述のように、Self-Expressive層は線形の全結合層とみなせるので、バックプロパゲーションにより一括で学習できる。この方法で $\bf{C}$ を求めたら後は従来通りにスペクトラルクラスタリング。結局のところ、オートエンコーダによって特徴抽出したデータに対して部分空間クラスタリングを行うが、その部分空間クラスタリング用の $\bf{C}$ もニューラルネット上に組み込んで、一緒に学習してしまおうという手法。

顔画像クラスタリングのベンチーマークデータセット（Extended Yale BとORL）や物体画像クラスタリングのベンチマークデータセット（COIL20とCOIL100）で実験。結果、他の手法（カーネル部分空間クラスタリングなど非線形クラスタリング手法含む）よりも高精度なクラスタリングが可能。ノルム項ついては、何故か $l2$ ノルムを使った方が良い結果に。

Mixture-Rank Matrix Approximation for Collaborative Filtering

協調フィルタリングでよく使われるLow-rank matrix approximationにおいて、低ランクと高ランクの近似行列を混合することにより精度を上げた手法。
実際、レイティング数が少ないアイテムやユーザに対する近似行列は低ランクにした方が推定結果は良く、逆にレイティング数に対する近似行列は高ランクにした方が推定結果は良い。

従来の手法はユーザ・アイテムの実際のレイティングマトリックスを $\bf{R}$ として、 $\bf{R}\approx \bf{\hat{R}} = \bf{U}\bf{V}^T$ （ $\bf{U}, \bf{V}$ ともにk列の行列）となる $\bf{U}, \bf{V}$ を求め、i番目のユーザのj番目のアイテムに対する嗜好スコアを $\bf{U}_i\bf{V}^T_j$ （それぞれ、 $\bf{U}$ のi行目、 $\bf{V}$ のj行目のベクトル）として推定する。PMFという手法では、この問題を確率的にモデル化してMAP推定により $\bf{U}, \bf{V}$ を求める手法を提供している。ここでkは固定の値とするため、ユーザおよびアイテムのレイティング数にかかわらず、ランクkの近似行列を推定する。

本論文では、このPMFをベースに複数のランクkを重み付きで混合させたモデルを提案。具体的には下記の確率モデルを構築している。
$\displaystyle\Pr({\bf R}|{\bf U},{\bf V},{\bf \alpha},{\bf \beta},\sigma^2)=\prod_{i=1}^m \prod_{j=1}^n \left[\sum_{k=1}^K \alpha_i^k \beta_j^k N(R_{i,j}|{\bf U}_i^k{{\bf V}_j^k}^T,\sigma^2) \right]^{ \mathbb{1}_{i,j}}$
ここで、 ${\textbf{U}_i^k}$ はランクkの行列 $\textbf{U}^k$ のi行目（ユーザiの特徴ベクトル）、 ${\textbf{V}_j^k}$ はランクkの行列 $\textbf{V}^k$ のj行目（アイテムjの特徴ベクトル）を表す。また、 $\alpha_i^k$ はユーザiにおけるランクkモデルの重み、 ${\beta_j^k}$ はアイテムjにおけるランクkモデルの重み。これにより、ランクkのモデルがユーザi、アイテムjにマッチするなら $\alpha_i^k,\beta_i^k$ は大きくなるはずであり、マッチしないなら $\alpha_i^k,\beta_i^k$ は小さくなるはず（期待的には0になる）。これらの $\bf{U},\bf{V},\bf{\alpha},{\bf \beta}$ をMAP推定により求めようというのが提案手法。

MAP推定のため事前分布を与える。 $\textbf{U}_i^k$ 、 $\textbf{V}_j^k$ の事前分布はそれぞれ、平均 $\textbf{0}$ 、共分散 $\sigma_U^2 \textbf{I}$ 、 $\sigma_V^2 \textbf{I}$ のガウス分布としている（行列 $\textbf{U},\textbf{V}$ の事前分布はこれらの総乗）。また、 ${\bf \alpha}^k, {\bf \beta}^k$ はそれぞれスパースなベクトルにしたいので、 $\alpha_i^k,\beta_i^k$ の事前分布はラプラス分布とする（それぞれ、平均 $\mu_{\alpha},\mu_{\beta}$ 、分散 $b_{\alpha},b_{\beta}$ というパラメータを持つ）。これらを積算することによって導き出される事後分布のlogにマイナスをつけた損失関数を最小化する。実際には、そのままの損失関数を最小化するのは難しいので、イェンセンの不等式で導き出される下限を最小化する。

$\textbf{U},\textbf{V},{\bf \alpha},{\bf \beta}$ は $\sigma,\sigma_U^2,\sigma_V^2,\mu_{\alpha},\mu_{\beta},b_{\alpha},b_{\beta}$ を固定して確率的勾配法で解く。一方、　 $\sigma,\sigma_U^2,\sigma_V^2,\mu_{\alpha},\mu_{\beta},b_{\alpha},b_{\beta}$ も学習したいので、今度は $\textbf{U},\textbf{V},{\bf \alpha},{\bf \beta}$ を固定して最尤推定する（こちらは、解析的に求まる）。これらの操作を収束するまで繰り返す。

MovieLensのデータセットやNetflix Prizeのデータセットで実験。結果、kを固定するPMFに比べて高い推定精度を得られるし、他の手法と比べても有意な差が得られる。また、混合するモデルの個数（k個数）はたくさんかつ細かくした方が良いようだ（計算時間はかかるけど）。

以下2つの論文については、まだほとんど理解できていない。とりあえず、備忘録として残す。

Kernel Feature Selection via Conditional Covariance Minimization

カーネル法を使って特徴とラベル間の非線形関係をとらえて特徴抽出する手法。条件付き共分散作要素のトレースを基準にすることを提案。（相互）共分散作要素は2つのヒルベルト空間間の共分散を表しており、そこから条件付き共分散作要素が導き出される。
論文では条件付き共分散作要素のトレースがラベルとの独立性を評価する指標となると主張し、その理由を説明している。また、一般に特徴選択の最適化をまともに行おうとするとNP困難な問題になってしまうので、ここでは離散最適化から連続最適化へ問題を緩和している。また、それでも計算量的にしんどいので、いくつか計算量を抑えるための目的関数、制約条件の変換も行っている。

共分散作要素をもう少し理解してから再チャレンジが必要。今のところ共分散行列のヒルベルト空間への拡張と理解。日本語だとこの辺りが良さげ。

Unsupervised Image-to-Image Translation Networks

2つのドメインの画像セットがあり、一方のドメインの画像をもう一方のドメインの画像に変換する手法を提案した論文。例えば、低解像画像から高解像画像へ変換するタスクや、朝の風景画像を夜の風景画像に変換するタスクなどがあげられる。Unspervisedなので、学習用に2つのドメイン間で対応する画像のセットが与えられるわけではなく、ただ単にドメインAの画像セットとドメインBの画像セットが与えられるのみ。
これを実現するために2つのドメインのデータは同じlatent spaceを共有するという仮定を置く。この仮定をもとに、variational autoencoder（VAE）とgenerative adversarial network（GAN）によりドメイン変換を行う。

まず、VAEでは2つのドメインのインプットが同じlatent spaceに変換されるように、エンコーダ層を学習。この同じlatent spaceへ変換というのを実現するために、エンコーダ層の最後の数層とデコーダ層の最初の数層はドメイン間で重みを共有するという制約をつける（なぜこれで、うまくいくのかがイマイチわからない）。
次にGANは、adversarial discriminatorが画像生成（ドメイン変換）のネットワークによって生成された画像か実際の画像かを高精度に判定できるように学習させる。画像生成ネットワーク側はそれに負けじと高精度なドメイン変換を行う。
これら2つを両立させる目的関数（VAEとGANと制約項の和）を最適化することにより、ドメイン変換ネットワークを学習する。

VAEやGANは名前しか聞いたことなかったが、面白い枠組みだし応用範囲も広そう。まだまだブームは続く？

所感

久々にNIPSの論文を読んでみましたが、やはりニューラルネット関連が多いという印象。それでもニューラルネットばかりだけではなく、未だにカーネル法に関する論文も結構あるので（「kernel」で論文タイトルの検索かけると20件ヒット）まだまだこの分野も廃れていない、今後盛り返す可能性も十分ありと感じました。
あと、VAEやGANなど生成モデル+ニューラルネットの枠組みは、今後いろいろな分野に波及しそうな気がします（GAN+トピックモデルとかVAE+HMMとか、もう誰かやってる？）。

2014-09-26

【Python】生存報告がてら昔書いたFFTのコードを晒す

お久しぶりです。ちゃんと生きてます。色々と立て込んでたり何だりでモチベーションが低下していました。このまま、フェードアウトしていくのもよろしくないと思いつつもネタがない・・・なんて思っていたところ、PCのデータを整理してたら、以前作成したFFT（高速フーリエ変換）のプログラムを発見。てことで、このソースコードを今回のネタにしようかなと思います。
なお、フーリエ変換やDFT（離散フーリエ変換）についての解説は割愛します。また、FFTについても多少は説明しますが、DFTからFFT アルゴリズムの導出過程など詳しい説明は省きます。また、データのサンプル数は2のべき乗個に限定します。

FFTのソースコード

ということで早速ソースコードです。まずは、モジュールのインポートとFFTの関数です。

import numpy as np
import matplotlib.pyplot as plt
import time

def fft(x):
    n = x.shape[0]
    X = np.zeros(n ,  dtype = 'complex')
    if n > 2:
        y = x[0:n:2] #(1)
        z = x[1:n:2]

        Y = fft(y) #(2)
        Z = fft(z)

        Wn = np.exp(-1j * (2 * np.pi * np.arange(0 , n/2))/float(n)) #(3)
        X[0:n/2] = Y + Wn * Z
        X[n/2:n] = Y - Wn * Z
        return X
    else:
        X[0] = x[0] + x[1] #(4)
        X[1] = x[0] - x[1]
        return X

モジュールはいつも通り、numpyをインポート。後は、後々の動作確認で使うモジュールです。
関数fftが今回の肝であるFFTを行う関数ですが、見ての通り、再帰を使って書いてます。通常、スタックオーバーフローとかを避けるために、ビットリバース処理なんかをして繰り返し構文で書くのですが、今回はそこまでやってません。こっちの方が簡単だし直感的なので再帰で書く方法を採用しています。
ソースコード内のコメントにあるように、ポイントは4つ。
まず、(1)で元の信号xから偶数番目だけ取り出した信号yと奇数番目だけ取り出した信号zを作ります。
(2)では、fft関数を再帰呼び出ししてyとzのFFTの結果Y,Zを得ています。このY,Zを利用すると元の信号xのFFTの結果Xは次の式で、表されます。
$\left\{ \begin{array}{l} X(i) = Y(i) + \exp \left( -j \frac{2 \pi i}{n} \right) Z(i) \: (0 \leq i < \frac{n}{2}) \\ X(\frac{n}{2} + i) = Y(i) - \exp \left( -j \frac{2 \pi i}{n} \right) Z(i) \: (0 \leq i < \frac{n}{2}) \end{array} \right.$

(3)はこの処理を行っています。n点のDFTをこの式のようにn/2点のDFTの和に書きなおすことができるってのが、FFTの勘所です。
そして、(1)、(2)の再帰呼び出し処理を繰り返していくと、最終的にn = 2に行きつきます（元のデータサンプル数が2のべき乗の場合）。ここが再帰の終着点で、2点のDFTは(4)のように凄く簡単な計算で求められます。
これだけです。numpy使えば簡単にFFTのプログラムが書けます。ただし、numpyには元々FFTを行うメソッドが用意されているので、普通はわざわざ書く必要はありません。処理時間もこっちの方が断然速いです。

動作確認

てことで、一応動作確認。

def main():
    #テスト用の信号生成
    n =2**10
    x = (3 * np.cos(250* 2*np.pi * np.arange(0,n)/n) +
         2 * np.sin(250* 2*np.pi * np.arange(0,n)/n) +
         np.cos(30* 2*np.pi * np.arange(0,n)/n) )

    tic = time.clock()
    X = fft(x) #FFTの計算
    toc = time.clock()
    print toc - tic #処理時間表示

    #以下、結果グラフ表示
    plt.figure()
    plt.plot(X.real , "r")
    plt.xlim((0, n))
    plt.xticks(np.arange(0, 1024 , 100))
    plt.hold(True)
    plt.plot(X.imag , "b")
    plt.show()

テスト用の信号は、データのサンプル数が1024個で、周波数250のcos波、sin波、周波数30のcos波をそれぞれ適当なゲインをかけて足し合わせたものを使います。で、結果はこちら。

赤が実数部、青が虚数部です。うん、それっぽい。numpyに用意されているfftメソッドも同様の結果を吐き出します。また処理時間は、1024点なら30[ms]程度で処理できます。しかし、既存のメソッドの方は0.2[ms]程度。速い・・・

まとめ

今回はFFTのプログラム（再帰版）について話しました。上述のように、FFTはループ処理で書くことの方が多いようで、僕が知っている参考書では、大抵そちらの方法のサンプルプログラムが記述されています。しかし、それらは全部C言語で書かれているものなので、他の言語のサンプルコードが載っているものだと、どうなのか気になるところです。
あと、今回はサンプル点数が2のべき乗個である場合という最も基本的なFFTのみ記述しましたが、サンプル点数が素数の場合のFFT アルゴリズムが提案されていたりと応用的なFFT アルゴリズムも数々あります。これらを調べてみるのも、面白いかもしれません。

関係ないですが、転職により愛知から大阪に引っ越ししてきました。新天地からの更新でした。

はじめに

スパースコーディングとは

画像再構成のための辞書学習

実装した手法の説明

実装したコード

画像再構成による超解像

実装した手法の説明

実装したコード

実験

83×83の画像を249×249の画像へ（拡大率3倍）

62×62の画像を248×248の画像へ（拡大率4倍）

まとめ

はじめに

近接勾配法、proximal operatorとは

Moreau decompositionと双対ノルム

双対ノルムについての補足

何故、ノルム関数の共役関数が指示関数なのか？

結局のところproimal operatorはどうやって求めればよいか？

様々な正則化項のproximal operator

l2ノルムの2乗正則化：

Elastic Net：

グループl1ノルム（グループに重複がない場合）：

トレースノルム ：

近接勾配法を実装

近接勾配法を行う関数

l1ノルム正則化項のproximal operator

lasso回帰問題を解く関数

proximal operatorを別なものにする

まとめ

応用例（2018/6/15追記）

参考文献

はじめに

定式化

解析解の導出と実装

複素数l1ノルム正則化

最後に

はじめに

離散フーリエ変換の線形回帰的定式化

l1ノルムによるスパース化

複素数の最適化問題から実数のみの最適化問題へ

cvxpyを用いた実装

実験

通常の離散フーリエ変換

Lassoによる離散フーリエ変換

逆変換

まとめ

小中大括弧でエラー

sum、min、max等の上置き,下置きが上付き、下付きになる

\bfにより意図しないところまでも太字化される

Deep Subspace Clustering Networks

Mixture-Rank Matrix Approximation for Collaborative Filtering

Kernel Feature Selection via Conditional Covariance Minimization

Unsupervised Image-to-Image Translation Networks

所感

FFTのソースコード

動作確認

まとめ

l2ノルムの2乗正則化： $g({\bf v}) = \lambda \| {\bf v} \|^2_2$

Elastic Net： $g({\bf v}) = \lambda_1 \| {\bf v} \|_1 + \lambda_2 \| {\bf v} \|^2_2$

グループl1ノルム（グループに重複がない場合）： $\displaystyle g({\bf v}) = \lambda \sum_{g \in G} \| {\bf v}_g \|_p$

トレースノルム： $g({\bf Z}) = \lambda \| {\bf Z} \|_{{\rm trace}} = {\rm Tr}\left(\sqrt{{\bf Z}^T {\bf Z}}\right)$