機器學習-04 基于sklearn 廣義線性模型- Lasso回歸
- Lasso回歸
- 坐标下降算法
- 官方手冊示例再現
- 重要代碼解釋
Lasso回歸
稀疏系數是指含零較多的系數。這種現象的産生可能是特征值設定的原因，比如性别男性為1女性為0，或者天氣晴天為1陰天為0，這種非黑及白的選擇如果有很多，可能會産生一溜零的情況。百度百科上面這段話寫的特别好，特地摘抄在下面“該方法是一種壓縮估計。它通過構造一個懲罰函數得到一個較為精煉的模型，使得它壓縮一些回歸系數，即強制系數絕對值之和小于某個固定值；同時設定一些回歸系數為零。是以保留了子集收縮的優點，是一種處理具有複共線性資料的有偏估計。”Lasso回歸就是拟合稀疏系數的線性模型。它傾向于使用具有較少參數值的情況，有效地減少給定解決方案所依賴變量的數量。 是以，Lasso 及其變體是壓縮感覺領域的基礎。
其最小化的目标函數是:
min ω 1 2 n s a m e p l e s ∥ X ω − y ∥ 2 2 + α ∥ ω ∥ 1 \underset{\omega }{\mathop{\min }}\,\frac{1}{2{{n}_{sameples}}}{{\left\| X\omega -y \right\|}_{2}}^{2}+\alpha {{\left\| \omega \right\|}_{1}} ωmin2nsameples1∥Xω−y∥22+α∥ω∥1
對其參數的求解有兩種方法坐标下降算法和最小角回歸 。
坐标下降算法
看了很多博文，感覺把坐标下降算法說清楚的百度百科，現摘錄如下：
“坐标下降法（coordinate descent）是一種非梯度優化算法。算法在每次疊代中，在目前點處沿一個坐标方向進行一維搜尋以求得一個函數的局部極小值。在整個過程中循環使用不同的坐标方向。對于不可拆分的函數而言，算法可能無法在較小的疊代步數中求得最優解。為了加速收斂，可以采用一個适當的坐标系，例如通過主成分分析獲得一個坐标間盡可能不互相關聯的新坐标系。”
官方手冊示例再現
示例一: Lasso和Elastic Net(彈性網絡)在稀疏信号上的表現：
import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import r2_score
# #############################################################################
# Generate some sparse data to play with
np.random.seed(42)
n_samples, n_features = 50, 100
X = np.random.randn(n_samples, n_features)
# Decreasing coef w. alternated signs for visualization
idx = np.arange(n_features)
coef = (-1) ** idx * np.exp(-idx / 10)
coef[10:] = 0 # sparsify coef
y = np.dot(X, coef)
# Add noise
y += 0.01 * np.random.normal(size=n_samples)
# Split data in train set and test set
n_samples = X.shape[0]
X_train, y_train = X[:n_samples // 2], y[:n_samples // 2]
X_test, y_test = X[n_samples // 2:], y[n_samples // 2:]
# #############################################################################
# Lasso
from sklearn.linear_model import Lasso
alpha = 0.1
lasso = Lasso(alpha=alpha)
y_pred_lasso = lasso.fit(X_train, y_train).predict(X_test)
r2_score_lasso = r2_score(y_test, y_pred_lasso)
print(lasso)
print("r^2 on test data : %f" % r2_score_lasso)
# #############################################################################
# ElasticNet
from sklearn.linear_model import ElasticNet
enet = ElasticNet(alpha=alpha, l1_ratio=0.7)
y_pred_enet = enet.fit(X_train, y_train).predict(X_test)
r2_score_enet = r2_score(y_test, y_pred_enet)
print(enet)
print("r^2 on test data : %f" % r2_score_enet)
m, s, _ = plt.stem(np.where(enet.coef_)[0], enet.coef_[enet.coef_ != 0],
markerfmt='x', label='Elastic net coefficients')
plt.setp([m, s], color="#2ca02c")
m, s, _ = plt.stem(np.where(lasso.coef_)[0], lasso.coef_[lasso.coef_ != 0],
markerfmt='x', label='Lasso coefficients')
plt.setp([m, s], color='#ff7f0e')
plt.stem(np.where(coef)[0], coef[coef != 0], label='true coefficients', markerfmt='bx')
plt.legend(loc='best')
plt.title("Lasso $R^2$: %.3f, Elastic Net $R^2$: %.3f"
% (r2_score_lasso, r2_score_enet))
plt.show()
示例二: 壓縮感覺_斷層重建：
層析成像投影操作是線性變換。除了對應于線性回歸的資料保真度項之外，我們還對圖像的L1範數進行了懲罰，以考慮其稀疏性。由此産生的優化問題稱為套索。我們使用class sklearn.linear_model.Lasso，它使用坐标下降算法。重要的是，與在此使用的投影運算符相比，該實作在稀疏矩陣上的計算效率更高。
# Author: Emmanuelle Gouillart <[email protected]>
# License: BSD 3 clause
import numpy as np
from scipy import sparse
from scipy import ndimage
from sklearn.linear_model import Lasso
from sklearn.linear_model import Ridge
import matplotlib.pyplot as plt
def _weights(x, dx=1, orig=0):
x = np.ravel(x)
floor_x = np.floor((x - orig) / dx).astype(np.int64)
alpha = (x - orig - floor_x * dx) / dx
return np.hstack((floor_x, floor_x + 1)), np.hstack((1 - alpha, alpha))
def _generate_center_coordinates(l_x):
X, Y = np.mgrid[:l_x, :l_x].astype(np.float64)
center = l_x / 2.
X += 0.5 - center
Y += 0.5 - center
return X, Y
def build_projection_operator(l_x, n_dir):
""" Compute the tomography design matrix.
Parameters
----------
l_x : int
linear size of image array
n_dir : int
number of angles at which projections are acquired.
Returns
-------
p : sparse matrix of shape (n_dir l_x, l_x**2)
"""
X, Y = _generate_center_coordinates(l_x)
angles = np.linspace(0, np.pi, n_dir, endpoint=False)
data_inds, weights, camera_inds = [], [], []
data_unravel_indices = np.arange(l_x ** 2)
data_unravel_indices = np.hstack((data_unravel_indices,
data_unravel_indices))
for i, angle in enumerate(angles):
Xrot = np.cos(angle) * X - np.sin(angle) * Y
inds, w = _weights(Xrot, dx=1, orig=X.min())
mask = np.logical_and(inds >= 0, inds < l_x)
weights += list(w[mask])
camera_inds += list(inds[mask] + i * l_x)
data_inds += list(data_unravel_indices[mask])
proj_operator = sparse.coo_matrix((weights, (camera_inds, data_inds)))
return proj_operator
def generate_synthetic_data():
""" Synthetic binary data """
rs = np.random.RandomState(0)
n_pts = 36
x, y = np.ogrid[0:l, 0:l]
mask_outer = (x - l / 2.) ** 2 + (y - l / 2.) ** 2 < (l / 2.) ** 2
mask = np.zeros((l, l))
points = l * rs.rand(2, n_pts)
mask[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1
mask = ndimage.gaussian_filter(mask, sigma=l / n_pts)
res = np.logical_and(mask > mask.mean(), mask_outer)
return np.logical_xor(res, ndimage.binary_erosion(res))
# Generate synthetic images, and projections
l = 128
proj_operator = build_projection_operator(l, l // 7)
data = generate_synthetic_data()
proj = proj_operator * data.ravel()[:, np.newaxis]
proj += 0.15 * np.random.randn(*proj.shape)
# Reconstruction with L2 (Ridge) penalization
rgr_ridge = Ridge(alpha=0.2)
rgr_ridge.fit(proj_operator, proj.ravel())
rec_l2 = rgr_ridge.coef_.reshape(l, l)
# Reconstruction with L1 (Lasso) penalization
# the best value of alpha was determined using cross validation
# with LassoCV
rgr_lasso = Lasso(alpha=0.001)
rgr_lasso.fit(proj_operator, proj.ravel())
rec_l1 = rgr_lasso.coef_.reshape(l, l)
plt.figure(figsize=(8, 3.3))
plt.subplot(131)
plt.imshow(data, cmap=plt.cm.gray, interpolation='nearest')
plt.axis('off')
plt.title('original image')
plt.subplot(132)
plt.imshow(rec_l2, cmap=plt.cm.gray, interpolation='nearest')
plt.title('L2 penalization')
plt.axis('off')
plt.subplot(133)
plt.imshow(rec_l1, cmap=plt.cm.gray, interpolation='nearest')
plt.title('L1 penalization')
plt.axis('off')
plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0,
right=1)
plt.show()
即使将噪聲添加到投影中，使用L1罰分進行的重建也會産生零誤差的結果（所有像素都成功标記為0或1）。相比之下，L2懲罰（sklearn.linear_model.Ridge）對像素産生大量标記錯誤。與L1懲罰相反，在重建的圖像上觀察到重要的僞影。特别要注意的是，圓形假象将角落中的像素分開，與中央磁盤相比，投影造成的投影更少。
重要代碼解釋
關于電子計算機掃描的輸入與輸出問題：
掃描所得資訊經計算而獲得每個體素的X射線衰減系數或吸收系數，再排列成矩陣，即數字矩陣（digital matrix），數字矩陣可存貯于磁盤或CD光牒中。經數字/模拟轉換器（digital/analog converter）把數字矩陣中的每個數字轉為由黑到白不等灰階的小方塊，即像素（pixel），并按矩陣排列，即構成CT圖像。是以，CT圖像是重建圖像。每個體素的X射線吸收系數可以通過不同的數學方法算出。
alpha = 0.1
lasso = Lasso(alpha=alpha)
y_pred_lasso = lasso.fit(X_train, y_train).predict(X_test)
- 聲明Lasso模型。
- 利用fit()進行模型訓練，使用predict()函數進行模型預測。