解答¶
1¶
∂∂wLρ(w,z,α)=∂∂w(12‖y−Xw‖2L2+λ‖z‖L1+αT(w−z)+ρ2‖w−z‖2L2)=(−XTy+XTXw)+α+ρ(w−z)=0∴(XTX+ρI)w=XTy−α+ρz∴argminwLρ(w,z,α)=(XTX+ρI)−1(XTy−α+ρz)2¶
z の正負で場合分けをすれば、
argminz{c|z|+12(z−z0)2}={z0−c(z>0)z0+c(0>z)0(otherwise.)={z0−c(z0>c)z0+c(−c>z0)0(otherwise.)3¶
argminziLρ(w,z,α)=argminzi{λ‖zi‖L1−αizi+ρ2(−2wizi+z2i)}=argminzi{λ‖zi‖L1+ρ2(z2i−2(wi+αiρ)zi)}=argminzi{λρ‖zi‖L1+12(zi−(wi+αiρ))2}これは、2 の式において
{c⟶λρz0⟶wi+αiρとした場合に対応する。以上より、
argminzLρ(w,z,α)=proxλρ|∗|(wi+αiρ)ここでは、「L1正則化項の下での線形回帰問題」 を考えた。全体の流れは以下
- 二乗和誤差関数 12‖y−Xw‖2 にL1正則化項を加えた誤差関数を定義する。（λ は手動で決定。）
- 解析的に最小解を求めるのが難しいので、新しい変数 z を代入してそれぞれ独立の変数としてみる。
- とはいえ z=w という関係は成り立っているので、ラグランジュ乗数 {αi} を導入して、制約条件を付け加える。 L(w,z,α)=12‖y−Xw‖2L2+λ‖z‖L1+αT(w−z)
- このままでも良いが、凸性を増すために、制約条件を二次の形で加える。この式を拡張ラグランジュ関数(Augmented Lagrangian)と呼ぶ。[参考:知能システム論 第3回] Lρ(w,z,α)=12‖y−Xw‖2L2+λ‖z‖L1+αT(w−z)+ρ‖w−z‖2L2
- 解析的に求めることができないので、拡張ラグランジュ関数の最小化と双対変数の勾配上昇を繰り返す。
- 拡張ラグランジュ関数の最小化 w⟵argminwLρ(w,z,α)with fixed z,αz⟵argminzLρ(w,z,α)with fixed w,α
- 双対変数の勾配上昇 α⟵α+ρ∇ω(α)=α+ρ(w−z)
実装¶
以下では、実際に「線形回帰」「L1正則化項の下での線形回帰」「L2正則化項の下での線形回帰」のそれぞれを実装し、違いや特徴を調べる。
In [1]:
import numpy as np
import matplotlib.pyplot as plt
In [2]:
from kerasy.utils.data_generator import generateSin, generateGausian
sample data¶
In [3]:
N = 10 # data size
xmin = 0
xmax = 1
seed = 0
In [4]:
X_test_ori = np.linspace(xmin,xmax,1000)
Y_test = np.sin(2*np.pi*X_test_ori)
In [5]:
X_train_ori,Y_train = generateSin(N, xmin=xmin, xmax=xmax, seed=seed)
_, Noise = generateGausian(N, x=X_train_ori, seed=seed)
Y_train += Noise
In [6]:
plt.scatter(X_train_ori,Y_train,edgecolors='blue',facecolor="white",label="train data")
plt.plot(X_test_ori,Y_test,color="black",label="$t=\sin(2\pi x)$")
plt.legend(), plt.grid()
plt.show()
Training¶
In [7]:
from kerasy.utils.preprocessing import PolynomialBaseTransformer
In [8]:
Ms = [2,4,8,16]
Linear Regression¶
In [9]:
fig=plt.figure(figsize=(16,6))
for i,M in enumerate(Ms):
axT=fig.add_subplot(2,4,i+1)
#=== Transform ===
phi = PolynomialBaseTransformer(M)
X_train = phi.transform(X_train_ori)
X_test = phi.transform(X_test_ori)
#=== Training ===
w = np.linalg.solve(X_train.T.dot(X_train), X_train.T.dot(Y_train))
Y_pred = X_test.dot(w)
#=== RSS ===
RSS = np.sqrt(np.sum(np.square(X_train.dot(w)-Y_train)))
# Plot.
axT.plot(X_test_ori,Y_test,color="black", label="$t=\sin(2\pi x)$")
axT.plot(X_test_ori,Y_pred,color="red", label=f"prediction ($M={M}$)")
axT.scatter(X_train_ori,Y_train,edgecolors='blue',facecolor="white",label="train data")
axT.set_title(f"Results (RSS=${RSS:.3f}$,M=${M}$)"), axT.set_xlabel("$x$"), axT.set_ylabel("$t$")
axT.set_xlim(0,1), axT.set_ylim(-1.2,1.2), axT.legend()
axB=fig.add_subplot(2,4,i+5)
axB.barh(np.arange(len(w)),w, color="black")
axB.set_ylim(0,max(Ms))
plt.tight_layout()
plt.savefig(f"LinearRegression_seed{seed}.png")
plt.show()
Linear Regression + L2 norm¶
In [10]:
lambda_ = 1e-4
In [11]:
fig=plt.figure(figsize=(16,6))
for i,M in enumerate(Ms):
axT=fig.add_subplot(2,4,i+1)
#=== Transform ===
phi = PolynomialBaseTransformer(M)
X_train = phi.transform(X_train_ori)
X_test = phi.transform(X_test_ori)
#=== Training ===
w = np.linalg.solve(lambda_*np.identity(M+1) + X_train.T.dot(X_train), X_train.T.dot(Y_train))
Y_pred = X_test.dot(w)
#=== RSS ===
RSS = np.sqrt(np.sum(np.square(X_train.dot(w)-Y_train)))
# Plot.
axT.plot(X_test_ori,Y_test,color="black", label="$t=\sin(2\pi x)$")
axT.plot(X_test_ori,Y_pred,color="red", label=f"prediction ($M={M}$)")
axT.scatter(X_train_ori,Y_train,edgecolors='blue',facecolor="white",label="train data")
axT.set_title(f"Results (RSS=${RSS:.3f}$,M=${M}$)"), axT.set_xlabel("$x$"), axT.set_ylabel("$t$")
axT.set_xlim(0,1), axT.set_ylim(-1.2,1.2), axT.legend()
axB=fig.add_subplot(2,4,i+5)
axB.barh(np.arange(len(w)),w, color="black")
axB.set_ylim(0,max(Ms))
plt.tight_layout()
plt.savefig(f"LinearRegression_L2_seed{seed}.png")
plt.show()
Linear Regression + L1 norm¶
In [12]:
lambda_ = 1e-3
rho = 1e-3
In [13]:
def prox(w,alpha,rho,lambda_):
z0 = w + alpha/rho
c = lambda_/rho
return z0-c if c<z0 else z0+c if z0<-c else 0
In [14]:
fig=plt.figure(figsize=(16,6))
for i,M in enumerate(Ms):
axT=fig.add_subplot(2,4,i+1)
#=== Transform ===
phi = PolynomialBaseTransformer(M)
X_train = phi.transform(X_train_ori)
X_test = phi.transform(X_test_ori)
#=== Training ===
alpha = np.ones(shape=(M+1))
z = np.ones(shape=(M+1))
while True:
w = np.linalg.solve(X_train.T.dot(X_train)+rho*np.identity(M+1), X_train.T.dot(Y_train)-alpha+rho*z)
z = np.asarray([prox(w_,alpha_,rho,lambda_) for w_,alpha_ in zip(w,alpha)])
alpha += rho*(w-z)
if np.sqrt(np.sum(np.square(w-z))) < 1e-9:
break
Y_pred = X_test.dot(w)
#=== RSS ===
RSS = np.sqrt(np.sum(np.square(X_train.dot(w)-Y_train)))
# Plot.
axT.plot(X_test_ori,Y_test,color="black", label="$t=\sin(2\pi x)$")
axT.plot(X_test_ori,Y_pred,color="red", label=f"prediction ($M={M}$)")
axT.scatter(X_train_ori,Y_train,edgecolors='blue',facecolor="white",label="train data")
axT.set_title(f"Results (RSS=${RSS:.3f}$,M=${M}$)"), axT.set_xlabel("$x$"), axT.set_ylabel("$t$")
axT.set_xlim(0,1), axT.set_ylim(-1.2,1.2), axT.legend()
axB=fig.add_subplot(2,4,i+5)
axB.barh(np.arange(len(w)),w, color="black")
axB.set_ylim(0,max(Ms))
plt.tight_layout()
plt.savefig(f"LinearRegression_L1_seed{seed}.png")
plt.show()
In [ ]: