簡述

凸優化會很詳細地講解這三個演算法，這個學期剛好有這門課。
這裡以期末的大作業的專案中的一個題目作為講解。

題目

考慮線性測量b=Ax+e，其中b為50維的測量值，A為50*100維的測量矩陣，x為100維的未知稀疏向量且稀疏度為5，e為50維的測量噪聲。從b和A中恢復x的一範數規範化最小二乘法模型（任務!!）
$min\mathrm{\mid }$

∣ A x − b ∣ ∣ 2 2

+ ( p / 2 ) ∣ ∣ x ∣ ∣ 1 min||Ax-b||_2^2 +(p/2)||x||_1

• p為非負的正則化引數。
• 提示：
  在本實驗中，設x的真值中的非零元素 服從標準正態分佈，A中的元素服從標準正態分佈，e中的元素服從均值為0，方差為0.1的高斯分佈。

要求使用的演算法：

• 鄰近點梯度下降法
• 交替方向乘子法
• 次梯度法

實驗部分

生成資料

• generate-data.py
• 儲存資料到二進位制檔案中

import numpy as np
import random
ASize = (50, 100)
XSize = 100
A = np.random.normal(0, 1, ASize)
X = np.zeros(XSize)
e = np.random.normal(0, 0.1, 50)
XIndex = random.sample(list(range(XSize)), 5)  # 5 稀疏度
for xi in XIndex:
    X[xi] = np.random.randn()

b = np.dot(A, X) + e

np.save("A.npy", A)
np.save("X.npy", X)
np.save("b.npy", b)

鄰近點梯度下降法

$min f_0(x) = s(x) + r(x)$

• s為光滑項
• r為非光滑項

演算法過程：

$x^{k+\frac{1}{2} } = x^k - \alpha * \nabla{s(x^k)} \\ x^{k+1} = argmin_x{r(x) +\frac{1}{2*\alpha} || x-x^{k+\frac{1}{2}}||^2}$

解析：

這一問本質上就是Lasso的鄰近點梯度下降問題。
代入之前的演算法過程，得到下面的結果(相比於原來的Lasso簡單的變下就好了~)

$x^{k+\frac{1}{2} } = x^k - 2* \alpha * A^T(Ax^k-b) \\ x^{k+1} = argmin_x{ \{(\frac{p}{2})||x||_1+\frac{1}{2*\alpha} || x-x^{k+\frac{1}{2}}||^2} \}$

• 求解argmin的方法–軟門限法

關於光滑的部分，直接求導，在不光滑的部分，就求次梯度。
然後關於每個分量部分不進行。可以參照程式碼中的片段，在中間就為0，在區間之外就是對應的一個正數~

實現領近點梯度法

import numpy as np

A = np.load('A.npy')
b = np.load('b.npy')
X = np.load('X.npy')

ASize = (50, 100)
BSize = 50
XSize = 100
alpha = 0.001
P_half = 0.01
Xk = np.zeros(XSize)
zero = np.zeros(XSize)
while True:
    Xk_half = Xk - alpha * np.dot(A.T, np.dot(A, Xk) - b)
    # 軟門限運算元
    Xk_new = zero.copy()
    for i in range(XSize):
        if Xk_half[i] < - alpha * P_half:
            Xk_new[i] = Xk_half[i] + alpha * P_half
        elif Xk_half[i] > alpha * P_half:
            Xk_new[i] = Xk_half[i] - alpha * P_half
    if np.linalg.norm(Xk_new - Xk, ord=2) < 1e-5:
        break
    else:
        Xk = Xk_new.copy()

print(Xk)
print(X)

附加上畫圖的程式碼

• 藍色的是演算法最優值的距離（最終的收斂點的距離）
• 黃色的是預測的真實值的距離（一開始生成的資料）
• 距離都用的是二範數

import matplotlib.pyplot as plt
import numpy as np

A = np.load('A.npy')
b = np.load('b.npy')
X = np.load('X.npy')

ASize = (50, 100)
BSize = 50
XSize = 100
alpha = 0.005
P_half = 0.01
Xk = np.zeros(XSize)
zero = np.zeros(XSize)

X_opt_dst_steps = []
X_dst_steps = []
while True:
    Xk_half = Xk - alpha * np.dot(A.T, np.dot(A, Xk) - b)
    # 軟門限運算元
    Xk_new = zero.copy()
    for i in range(XSize):
        if Xk_half[i] < - alpha * P_half:
            Xk_new[i] = Xk_half[i] + alpha * P_half
        elif Xk_half[i] > alpha * P_half:
            Xk_new[i] = Xk_half[i] - alpha * P_half
    X_dst_steps.append(np.linalg.norm(Xk_new - X, ord=2))
    X_opt_dst_steps.append(Xk_new)
    if np.linalg.norm(Xk_new - Xk, ord=2) < 1e-5:
        break
    else:
        Xk = Xk_new.copy()

print(Xk)
print(X)

X_opt = X_opt_dst_steps[-1]

for i, data in enumerate(X_opt_dst_steps):
    X_opt_dst_steps[i] = np.linalg.norm(data - X_opt, ord=2)
plt.title("Distance")
plt.plot(X_opt_dst_steps, label='X-opt-distance')
plt.plot(X_dst_steps, label='X-real-distance')
plt.legend()
plt.show()

交替方向乘子法

$min f_1(x)+f_2(y) \\ s.t. Ax+By=d$

演算法過程：

$(x^{k+1}, y^{k+1}) = argmin_{x, y} \{L_c(x, y, v^k)\} \\ v^{k+1} =v^k +c(Ax^{k+1} + B y^{k+1})$

相關推薦