理論ばっかり追っていて悶々してきたので、IRLSでL1残差最小化が解けないか実験してみる。
- 第5章 厳密解から近似解へ に『スパースモデリング』5章のPython実装あり。
- スパースモデリング：第3章 追跡アルゴリズム は『スパースモデリング』3章のPython実装。
IRLSの実装は カレル大学卒論 を参考に。Pythonで簡単にできた。
import numpy
# IRLS法によりPhi @ x = yのスパース解を求める
def irls_update(Phi, x, y, order):
EPSILON = 10 ** (-8.0)
# 重みの計算
weight = numpy.abs(y - Phi @ x).flatten()
# 小さくなりすぎた重みは打ち切る
weight[weight < EPSILON] = EPSILON
# 対角行列に展開
W = numpy.diag(weight ** (order - 2))
# 更新後の係数: Phi.T @ W @ Phi @ x = Phi.T @ W @ y の解
return numpy.linalg.solve(Phi.T @ W @ Phi, Phi.T @ W @ y)
if __name__ == "__main__":
DIMENSION = 2
NUM_SAMPLES = 100
NUM_ITERATION = 50
# 解ベクトル
X_ANSWER = numpy.array([0.5, 0.5]).reshape((DIMENSION, 1))
x = numpy.zeros((DIMENSION, 1))
xhistory = numpy.zeros((DIMENSION, NUM_ITERATION))
# 観測を生成
Phi = numpy.random.rand(NUM_SAMPLES, DIMENSION)
y = Phi @ X_ANSWER
# 加法的雑音を重畳
# yrand = y + numpy.random.normal(0, 0.3, (NUM_SAMPLES, 1))
yrand = y + numpy.random.laplace(0, 0.3, (NUM_SAMPLES, 1))
error = numpy.zeros(NUM_ITERATION)
emp_error = numpy.zeros(NUM_ITERATION)
# IRLSを繰り返し適用
for count in range(NUM_ITERATION):
x = irls_update(Phi, x, yrand, 1)
xhistory[:, count] = x.reshape(2)
error[count] = numpy.linalg.norm(y - Phi @ x, ord=1) / NUM_SAMPLES
emp_error[count] = numpy.linalg.norm(yrand - Phi @ x, ord=1) / NUM_SAMPLES
実装は楽だったけど、誤差解析が沼。
- 誤差を重畳してみると、真の誤差と経験誤差が当然一致しない。
- 経験誤差的には局所解に入っている印象。
- サンプル数が少ないと大域最小解に入らないケースあり（経験誤差曲面の最小値が真の誤差の曲面の最小値に不一致）
- 経験誤差の曲面は二次曲線に見える。（2次式の最小化を考えているから当然のはず。）
- 最小二乗解よりも誤差が悪い時がある。最小二乗解はorder=2とすれば良くて、その時重み行列Wは単位行列になり、普通の最小二乗法と一致。
思いつき:
- IRLSは評価関数の最小化を考える時閉形式で求まるので何も考えない。パラメータに関してもう一度微分できるのでニュートン法使えそう。
- フィルタのときのように逐次的に求められない？
- パラメータ全てではなく1こずつ。サンプルについても1こずつ。更新していく。評価関数の最小化は平均値の最小化に見受けられるので、逐次的に更新しても良いように見える。
今日は遅いのでもう寝る。