sklearn学习（1）-线性回归(2)

  

3.Lasso（拟合系数稀疏的线性模型）

3.1 Lasso基础

Lasso 是拟合稀疏系数的线性模型。 它在一些情况下是有用的，因为它倾向于使用具有较少参数值的情况，有效地减少给定解决方案所依赖变量的数量。 因此，Lasso 及其变体是压缩感知领域的基础。 在一定条件下，它可以恢复一组非零权重的精确集（见压缩感知_断层重建）。

在数学公式表达上，它由一个带有 ${\ell }_1$先验的正则项的线性模型组成。 其最小化的目标函数是:

$$\begin{array}{c}\underset{w}{\min }\frac{1}{2n_{samples}}||Xw-y|{|}_2^2+\alpha ||w|{|}_1\end{array}\text{\hspace{1em}(3)}$$

lasso estimate 解决了加上罚项 $\alpha ||w|{|}_1$ 的最小二乘法的最小化，其中， $\alpha$ 是一个常数， $||w|{|}_1$ 是参数向量的 ${\ell }_1-norm$ 范数。

Lasso 类的实现使用了 coordinate descent （坐标下降算法）来拟合系数。

reg = linear_model.Lasso(alpha = 0.1)
reg.fit([[0, 0], [1, 1]], [0, 1])
print("reg.predict([[1, 1]])",reg.predict([[1, 1]]))
# reg.predict([[1, 1]]) [0.8]

对于较简单的任务，同样有用的是函数 lasso_path 。它能够通过搜索所有可能的路径上的值来计算系数。

3.2 Lasso的正则化参数设置

3.2.1 使用交叉验证

scikit-learn 通过交叉验证来公开设置 Lasso alpha 参数的对象: LassoCV 和 LassoLarsCV。 LassoLarsCV 是基于下面将要提到的 最小角回归 算法。 （LassoLarsCV 常用一些）

对于具有许多线性回归的高维数据集， LassoCV 最常见。 然而，LassoLarsCV 在寻找 alpha参数值上更具有优势，而且如果样本数量比特征数量少得多时，通常 LassoLarsCV 比 LassoCV 要快。

3.2.1 基于信息标准的模型选择

有多种选择时，估计器 LassoLarsIC 建议使用 Akaike information criterion （Akaike 信息判据）（AIC）或 Bayes Information criterion （贝叶斯信息判据）（BIC）。 当使用 k-fold 交叉验证时，正则化路径只计算一次而不是 k + 1 次，所以找到 α 的最优值是一种计算上更经济的替代方法。 然而，这样的判据需要对解决方案的*度进行适当的估计，它会假设模型是正确的，对大样本（渐近结果）进行导出，即数据实际上是由该模型生成的。 当问题严重受限（比样本更多的特征）时，它们也容易崩溃。

3.2.1 与 SVM（支持向量机） 的正则化参数的比较

alpha 和 SVM 的正则化参数C 之间的等式关系是 $alpha=1/C$ 或者 $alpha=1/(n_{s}amples\ast C)$ ，并依赖于估计器和模型优化的确切的目标函数。

3.3 Lasso的补充

• 由于 Lasso 回归产生稀疏模型，因此可以用于执行特征选择，详见 基于 L1 的特征选取 。

• 断层重建-这个还听有趣的，更换了一个距离函数但效果迥然不同。

运行如下代码：

# Author: Emmanuelle Gouillart <emmanuelle.gouillart@nsup.org>
# License: BSD 3 clause

import matplotlib.pyplot as plt
import numpy as np
from scipy import ndimage, sparse

from sklearn.linear_model import Lasso, Ridge


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.0
    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.0) ** 2 + (y - l / 2.0) ** 2 < (l / 2.0) ** 2
    mask = np.zeros((l, l))
    points = l * rs.rand(2, n_pts)
    mask[(points[0]).astype(int), (points[1]).astype(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()

输出结果：

图上可以看出 $\ell 1$ 的效果远优于$\ell 2$ .

4.多任务Lasso

4.1 基本概念及公式

MultiTaskLasso 是一个估计多元回归稀疏系数的线性模型： y 是一个形状为(n_samples, n_tasks) 的二维数组，其约束条件和其他回归问题（也称为任务）是一样的，都是所选的特征值。

在数学上，它由一个线性模型组成，以混合的$ \ell_1$ ${\ell }_2$ 作为正则化器进行训练。目标函数最小化是：

$$\begin{array}{c}\underset{w}{\min }\frac{1}{2n_{samples}}||XW-Y|{|}_{Fro}^2+\alpha ||W|{|}_{21}\end{array}\text{\hspace{1em}(4)}$$

其中 $Fro$ 表示 Frobenius 标准（弗罗贝尼乌斯准则）：

$$\begin{array}{c}||A|{|}_{Fro}=\sqrt{\sum _{ij}{a}_{ij}^2}\end{array}\text{\hspace{1em}(5)}$$

并且${\ell }_1$ ${\ell }_2$ 读取为:

$$\begin{array}{c}||A|{|}_{21}=\sum _i\sqrt{\sum _j{a}_{ij}^2}\end{array}\text{\hspace{1em}(6)}$$

MultiTaskLasso 类的实现使用了坐标下降作为拟合系数的算法。

4.2 Lasso与MultiTaskLasso的对比

• 下图比较了通过使用简单的 Lasso 或 MultiTaskLasso 得到的 W 中非零的位置。 Lasso 估计产生分散的非零值，而 MultiTaskLasso 的一整列都是非零的。

import numpy as np
from sklearn.linear_model import Lasso, MultiTaskLasso
import matplotlib.pyplot as plt

# 生成二维数据
rng = np.random.RandomState(42)

# Generate some 2D coefficients with sine waves with random frequency and phase
n_samples, n_features, n_tasks = 100, 30, 40
n_relevant_features = 5
coef = np.zeros((n_tasks, n_features))
times = np.linspace(0, 2 * np.pi, n_tasks)
for k in range(n_relevant_features):
    coef[:, k] = np.sin((1.0 + rng.randn(1)) * times + 3 * rng.randn(1))

X = rng.randn(n_samples, n_features)
Y = np.dot(X, coef.T) + rng.randn(n_samples, n_tasks)

# 创建线性模型
# Lasso
coef_lasso_ = np.array([Lasso(alpha=0.5).fit(X, y).coef_ for y in Y.T])
# MultiTaskLasso
coef_multi_task_lasso_ = MultiTaskLasso(alpha=1.0).fit(X, Y).coef_

# 进行多任务Lasso模型训练
fig = plt.figure(figsize=(8, 5))
plt.subplot(1, 2, 1)
plt.spy(coef_lasso_)
plt.xlabel("Feature")
plt.ylabel("Time (or Task)")
plt.text(10, 5, "Lasso")
plt.subplot(1, 2, 2)
plt.spy(coef_multi_task_lasso_)
plt.xlabel("Feature")
plt.ylabel("Time (or Task)")
plt.text(10, 5, "MultiTaskLasso")
fig.suptitle("Coefficient non-zero location")

feature_to_plot = 0
plt.figure()
lw = 2
plt.plot(coef[:, feature_to_plot], color="seagreen", linewidth=lw, label="Ground truth")
plt.plot(
    coef_lasso_[:, feature_to_plot], color="cornflowerblue", linewidth=lw, label="Lasso"
)
plt.plot(
    coef_multi_task_lasso_[:, feature_to_plot],
    color="gold",
    linewidth=lw,
    label="MultiTaskLasso",
)
plt.legend(loc="upper center")
plt.axis("tight")
plt.ylim([-1.1, 1.1])
plt.show()

• 从上图可以看出，MultiTaskLasso的拟合结果对比Lasso的结果，明显要更为贴合真实数据。