構造正則化② : メジャライザとISTA、FISTA・その2

前回、凸関数 のメジャライザ

  を導き、これが

を満たすことを説明しました。ここで は点 における凸関数 f の曲率の上限です。

ところでメジャライザのノルム中の定ベクトルを

式(1)

と書けば

となります。

更にL1正則化問題の場合

なのでL1正則化問題の目的関数

に対するメジャライザは

となります。よってL1正則化問題の場合も上のメジャライザを最小化、その点でメジャライザを構築し直し最小化…を繰り返すことで最適解が得られます。ちなみにその際定数は何の寄与も持たないので

を最小化することになります。

ところで、この解は解析的に求められ

式(2)

となります。ここで は以下のsoft-threshold関数です。

以上からISTAのアルゴリズムは以下のように求められます。

===== ISTAのアルゴリズム =======

1.  を適当な値で初期化

2. 現在の について式(1)の z を構築

3. 式(2)で  を更新

4. 収束判定クリアなら終了、そうでなければ2. へ

これがISTAの基本的な解析フローになります。

=======================================

これに基づいて、ISTAの1サイクルを実行する機能を持ったクラスをpythonで作成してみました(なお実効するには予め numpy (あとmatplotlib.pyplot)をインポートしておく必要が有ります)。

ここでは、データ y、伝達行列 X に対し以下を満たす β を求めるL1正則化問題を考えます。

このため目的関数

を最小化することを考えます。

この場合、これまで f(β) と表記していた損失関数とその微分は具体的に

式(3)

となります。

class ISTA:
    def __init__(self):
        self.lambda_val = 0. # 正則化パラメータ
        self.Lipschitz = 0.  # 曲率上限（リプシッツ定数）
        self.X = None    # 伝達行列
        self.y = None    # 観測データ
        self.beta = None # 解ベクトル
        
        self.z = None  # ベクトル z
        self.r = None  # 残差ベクトル: r = y - X * beta

    def set(self, X, y, lambda_val, K):
        '''
        解くべき連立方程式を設定
        '''
        self.lambda_val = lambda_val
        self.X = X
        self.y = y
        self.Lipschitz = K

        self.beta = np.zeros(self.X.shape[1])

    def soft_threshold(self, g, l):
        '''
        ソフト閾値関数
        '''
        return np.sign(g) * np.maximum(np.abs(g) - l, 0.)

    def update_r(self):
        '''
        残差ベクトル r を更新
        '''
        self.r = self.y - self.X.dot(self.beta)

    def update_z(self):
        '''
        ベクトル z を更新
        '''
        self.z = self.beta + self.X.T.dot(self.r) / self.Lipschitz

    def update(self):
        '''
        ISTAの1サイクルの更新:
        r 更新 → z 更新 → ソフト閾値関数適用で beta 更新
        '''
        self.update_r()
        self.update_z()
        self.beta = self.soft_threshold(self.z, self.lambda_val / self.Lipschitz)

これをテストデータで動かしてみましょう。

n x n (下の場合 n = 100) のランダムな行列を X0 とし、真の解を

40 - 70番目の要素が 1、その他が 0

であるベクトル beta0 とします。

また入力データ y0 には X0.dot(beta0) を用います。

またISTAのアルゴリズムに乗せるため X0 は標準化、 y0 は中心化し、それぞれ X、yとしています。

### テストデータ作成 ###
# n x n のランダム行列 X0
n = 100
X0 = np.random.rand(n, n)

# 解ベクトル
# 40 - 70 番目要素が1、それ以外0
beta0 = np.zeros(n)
beta0[40:70] = 1.

# テストデータ: y0 = X0 * beta0
y0 = X0.dot(beta0)

# X0 の標準化: 中心化 + 規格化
X = X0 - np.mean(X0, axis=0)
w = np.linalg.norm(X, axis=0)
X = X / w

# y0 の中心化
y = y0 - y0.mean()

こうして作成した X と y を入力としてISTAで beta0 の再現を試みました。

ここで定数 K をどうやって求めるか、ということが問題になります。

これまで説明してきたように K は 点βにおけるf(β)の曲率の上限値でした。

f が式(3)でかける場合、その曲率は

なのでこの行列の最大固有値（より少し大きな値）が曲率の上限になります。

そこで下では np.linalg.eig で X.T.dot(X) を固有値分解して最大値を求め、その ceil を K としています。

また正則化パラメータ λ は 0.01 として計算しました。

# 固有値の上限を計算
v, _ = np.linalg.eig(X.T.dot(X))
K = np.ceil(v.max())
print(v.max(), K)

# ISTAオブジェクト作成
obj = ISTA()
obj.set(X, y, 0.01, K)

# ISTAのイテレーション
# 解ベクトルの修正量が閾値(epsilon)を下回る
# 或いはイテレーション回数が規定値(maxniter)を上回る
# なら終了
epsilon = 1.e-5
maxniter = 1000

niter = 0
while True:
    dbeta = obj.beta.copy()
    obj.update()
    
    dbeta -= obj.beta
    delta = np.linalg.norm(dbeta)
    if delta < epsilon:
        break
    niter += 1
    if niter >= maxniter:
        break

plt.plot(beta0, '.')
plt.plot(obj.beta / w, '.')

上の計算では epsilon = 1.e-5、maxniter = 1000 としています。

この計算により下の結果が得られました。

ここで横軸が β の要素番号、縦軸が β の各要素の値です。

青い点が真のモデル beta0、オレンジが得られた解(を、X の標準化の際 X の各列を規格化した効果を取り除いたもの) beta / w ( wは X の各列ベクトルのノルムを格納したもの)を表します。

オレンジが青を良く再現できているのが分かります。

なおこの計算でのイテレーション回数 (niter) は 620 回でした（乱数のseedを設定していないので毎回結果は異なりますが…）。

=======================================

さて、前回の記事に続き今回でISTAの概念的な解析フローとサンプルの愚直なpythonクラスを提示しました。

ところで、実用上では定数 K をどう求めるか、という事が問題になってきます。

上で述べたように今回は固有値分解のソルバーで直接 K を計算していますが、大規模問題ではここで大きな計算コストがかかる、或いは問題が大きすぎてそもそも求められない、といった場合も出てきます。

またこれまで軽く流していますが、この定数はLipschitz （リプシッツ）定数と呼ばれていて、本来 Lipschitz 連続(リプシッツ連続)という、関数の連続性の拡大解釈に基づき説明されます。

この辺りについての説明も、あまり理論に偏らない範囲で次回以降のどこかで説明したいと思っています（自分自身への備忘録的なドキュメントとして！）。