线性回归模型-Lasso

Lasso是一个能产生系数系数的线性模型. Lasso在某些情况下很有用，因为它倾向于选择具有较少参数值的解决方案，当变量之间是线性相关的, 可以有效地减少变量数量。由于这个原因，Lasso及其变体是压缩感知领域的基础。在某些条件下，它可以返回精确的非零参数集.
数学上, 它有一个带有L1正则的线性模型组成. 需要最小化的损失函数:
$\underset{w}{\min }\frac{1}{2n_{samples}}{\Vert Xw-y\Vert }_2^2+\alpha {\Vert w\Vert }_1$
Lasso算法计算最小二乘法加上$\alpha {\Vert w\Vert }_1$的惩罚项, $\alpha$是一个常数, ${\Vert w\Vert }_1$是参数向量的L1的模.
sklearn中Lasso类中使用坐标下降法作为算法来训练参数.

示例-sklearn Lasso类简单使用

代码:

  1 from sklearn import linear_model
  2 reg = linear_model.Lasso(alpha = 0.1)
  3 reg.fit([[0, 0], [1, 1]], [0, 1])
  4 '''
  5 Lasso(alpha=0.1, copy_X=True, fit_intercept=True, max_iter=1000,
  6    normalize=False, positive=False, precompute=False, random_state=None,
  7    selection='cyclic', tol=0.0001, warm_start=False)
  8 '''
  9 print(reg.predict([[1, 1]]))
 10 print(reg.coef_)
 11 print(reg.intercept_)

运行结果

[0.8]
[0.6 0. ]
0.2

示例-Lasso and Elastic Net for Sparse Signals

在手动产生的带噪声的稀疏信号数据上评估Lasso和Elastic-Net两种回归模型. 得到的系数与ground-truth比较.

代码:

  1 print(__doc__)
  2 
  3 import numpy as np
  4 import matplotlib.pyplot as plt
  5 
  6 from sklearn.metrics import r2_score
  7 
  8 # Generate some sparse data to play with
  9 
 10 np.random.seed(42)
 11 
 12 n_samples, n_features = 50, 200
 13 X = np.random.randn(n_samples, n_features)  # Generate randomly X
 14 coef = 3 * np.random.randn(n_features)  # Generate randomly coefficients
 15 inds = np.arange(n_features) # Generate indexs 
 16 np.random.shuffle(inds) # Indexs shuffle
 17 coef[inds[10:]] = 0 # Sparsify coefficients
 18 y = np.dot(X, coef)
 19 
 20 # Add noise
 21 y += 0.01 * np.random.normal(size = n_samples)
 22 
 23 # Split data in train set and test set
 24 n_samples = X.shape[0]
 25 X_train, y_train = X[:n_samples // 2], y[:n_samples // 2]
 26 X_test, y_test = X[n_samples // 2:], y[n_samples // 2:]
 27 
 28 # Lasso
 29 from sklearn.linear_model import Lasso
 30 
 31 alpha = 0.1
 32 lasso = Lasso(alpha = alpha)
 33 
 34 y_pred_lasso = lasso.fit(X_train, y_train).predict(X_test)
 35 r2_score_lasso = r2_score(y_test, y_pred_lasso)
 36 print(lasso)
 37 print('r^2 on test data: %f' % r2_score_lasso)
 38 
 39 # ElasticNet
 40 from sklearn.linear_model import ElasticNet
 41 
 42 enet = ElasticNet(alpha = alpha, l1_ratio = 0.7)
 43 
 44 y_pred_enet = enet.fit(X_train, y_train).predict(X_test)
 45 r2_score_enet = r2_score(y_test, y_pred_enet)
 46 print(enet)
 47 print('r^2 on test data: %f' % r2_score_enet)
 48 
 49 plt.plot(lasso.coef_, color = 'gold', linewidth = 2, label = 'Lasso coefficients')
 50 plt.plot(enet.coef_, color = 'lightgreen', linewidth = 2, label = 'Elastic net coefficients')
 51 plt.plot(coef, '--', color = 'navy', label = 'original coefficients')
 52 plt.legend(loc = 'best')
 53 plt.title('Lasso R^2: %f, Elastic Net R^2: %f' % (r2_score_lasso, r2_score_enet))
 54 plt.show()

运行结果:

Lasso(alpha=0.1, copy_X=True, fit_intercept=True, max_iter=1000,
   normalize=False, positive=False, precompute=False, random_state=None,
   selection='cyclic', tol=0.0001, warm_start=False)
r^2 on test data: 0.385982
ElasticNet(alpha=0.1, copy_X=True, fit_intercept=True, l1_ratio=0.7,
      max_iter=1000, normalize=False, positive=False, precompute=False,
      random_state=None, selection='cyclic', tol=0.0001, warm_start=False)
r^2 on test data: 0.240498

示例-Compressive sensing: tomography reconstruction with L1 prior (Lasso)

有利于加深对L1和L2正则的区别的理解
This example shows the reconstruction of an image from a set of parallel projections, acquired along different angles. Such a dataset is acquired in computed tomography (CT).

Without any prior information on the sample, the number of projections required to reconstruct the image is of the order of the linear size l of the image (in pixels). For simplicity we consider here a sparse image, where only pixels on the boundary of objects have a non-zero value. Such data could correspond for example to a cellular material. Note however that most images are sparse in a different basis, such as the Haar wavelets. Only l/7 projections are acquired, therefore it is necessary to use prior information available on the sample (its sparsity): this is an example of compressive sensing.

The tomography projection operation is a linear transformation. In addition to the data-fidelity term corresponding to a linear regression, we penalize the L1 norm of the image to account for its sparsity. The resulting optimization problem is called the Lasso. We use the class sklearn.linear_model.Lasso, that uses the coordinate descent algorithm. Importantly, this implementation is more computationally efficient on a sparse matrix, than the projection operator used here.

The reconstruction with L1 penalization gives a result with zero error (all pixels are successfully labeled with 0 or 1), even if noise was added to the projections. In comparison, an L2 penalization (sklearn.linear_model.Ridge) produces a large number of labeling errors for the pixels. Important artifacts are observed on the reconstructed image, contrary to the L1 penalization. Note in particular the circular artifact separating the pixels in the corners, that have contributed to fewer projections than the central disk.

代码

 45     angles = np.linspace(0, np.pi, n_dir, endpoint = False)
 46     data_inds, weights, camera_inds = [], [], []
 47     data_unravel_indices = np.arange(l_x ** 2)
 48     data_unravel_indices = np.hstack((data_unravel_indices, data_unravel_indices))
 49     
 50     for i, angle in enumerate(angles):
 51         Xrot = np.cos(angle) * X - np.sin(angle) * Y
 52         inds, w = _weights(Xrot, dx = 1, orig = X.min())
 53         mask = np.logical_and(inds >= 0, inds < l_x)
 54         weights += list(w[mask])
 55         camera_inds += list(inds[mask] + i * l_x)
 56         data_inds += list(data_unravel_indices[mask])
 57     proj_operator = sparse.coo_matrix((weights, (camera_inds, data_inds)))
 58     return proj_operator
 59 
 60 def generate_synthetic_data():
 61     """ Syntheic binary data"""
 62     rs = np.random.RandomState(0)
 63     n_pts = 36
 64     x, y = np.ogrid[0:l, 0:l]
 65     mask_outer = (x - l / 2.) ** 2 + (y - l / 2.) ** 2 < (l / 2. ) ** 2
 66     mask = np.zeros((l, l))
 67     points = l * rs.rand(2, n_pts)
 68     mask[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1
 69     mask = ndimage.gaussian_filter(mask, sigma = l / n_pts)
 70     res = np.logical_and(mask > mask.mean(), mask_outer)
 71     return np.logical_xor(res, ndimage.binary_erosion(res))
 72 
 73 # Generate synthetic images, and projections
 74 l = 128
 75 proj_operator = build_projection_operator(l, l // 7)
 76 data = generate_synthetic_data()
 77 proj = proj_operator * data.ravel()[:, np.newaxis]
 78 proj += 0.15 * np.random.randn(*proj.shape)
 79 
 80 # Reconstruction with L2(Ridge) penalization
 81 rgr_ridge = Ridge(alpha = 0.2)
 82 rgr_ridge.fit(proj_operator, proj.ravel())
 83 rec_l2 = rgr_ridge.coef_.reshape(l, l)
 84 
 85 # Reconstruction with L1(Lasso) penalization
 86 # the best value of alpha was determined using cross validation with LassoCV
 87 rgr_lasso = Lasso(alpha = 0.001)
 88 rgr_lasso.fit(proj_operator, proj.ravel())
 89 rec_l1 = rgr_lasso.coef_.reshape(l, l)
 90 
 91 plt.figure(figsize = (8, 3.3))
 92 plt.subplot(131)
 93 plt.imshow(data, cmap = plt.cm.gray, interpolation = 'nearest')
 94 plt.axis('off')
 95 plt.title('original image')
 96 plt.subplot(132)
 97 plt.imshow(rec_l2, cmap = plt.cm.gray, interpolation = 'nearest')
 98 plt.title('L2 penalization')
 99 plt.axis('off')
100 plt.subplot(133)
101 plt.imshow(rec_l1, cmap = plt.cm.gray, interpolation = 'nearest')
102 plt.title('L1 penalization')
103 plt.axis('off')
104 
105 plt.subplots_adjust(hspace = 0.01, wspace = 0.01, top = 1, bottom = 0, left = 0, right = 1)
106 plt.show()                                                                                            

运行结果:

Reference

• numpy.random.randn
• numpy.random.normal
• numpy.random.shuffle