【machine learning】线性模型（二）Lasso代码实例

稀疏信号上的LASSO与弹性网

使用估计Lasso和弹性网络回归对 人工破坏并加噪的稀疏信号进行建模。然后将估计系数与真实值作了比较。

创建一个随机生成的标准正态分布的二维数组X，数据集的维度50*100.

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)

使用 numpy.arange() 函数生成一个从 0 到 n_features - 1 的整数数组，表示特征的索引。

(-1) ** idx 计算权重向量的符号，使得奇数索引对应的权重为负，偶数索引对应的权重为正。np.exp() 给每个权重附加一个衰减因子，衰减因子随特征索引增大而减小（指数函数值越小，衰减越快）。这一步骤通常用于模拟某些特征对目标变量影响的强弱和方向。

将 coef 向量中索引大于等于10的部分设置为0，实现权重向量的稀疏化。这意味着在生成的目标变量 y 中，仅前10个特征对 y 有非零贡献，其余特征的影响被忽略。这种操作通常用于模拟实际问题中部分特征对目标变量影响显著，而其他特征影响较弱甚至无关的情形。

# 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:]

alpha 是Lasso回归中的正则化参数。模型系数惩罚程度，值越大，对系数的收缩力度越大，模型会更倾向于产生稀疏解（即更多系数为零）。
创建Lasso模型，它在最小化残差平方和的同时，引入L1范数惩罚项来促使某些特征的系数接近于零，实现特征选择和防止过拟合。
训练模型并进行预测：略
计算R²分数
R²值介于0到1之间，其中1表示完美拟合（预测值与真实值完全一致），0表示模型的预测能力与仅用均值预测的效果相当，负值则表明模型的预测效果还不如仅用均值预测。r2_score() 函数接收实际的测试集标签 y_test 和模型预测值 y_pred_lasso 作为参数，返回R²分数 r2_score_lasso，用来量化模型在测试集上的表现。

# 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)

压缩感知：基于L1(Lasso)先验层析重建

简介

压缩感知（Compressed Sensing, CS）是一种信号处理和采样理论，它允许在信号采样数远低于奈奎斯特定理要求的采样率（即信号的Nyquist采样率）的情况下，通过有效的数学方法精确重构原始信号。基于L1范数（Lasso先验）的层析重建是压缩感知在图像重建领域的一种应用，特别是在医学影像（如CT、MRI）和遥感成像等场景中。

基本原理

利用信号的稀疏性（或可压缩性）来实现低采样率下的有效重构。其主要包含三个步骤：

稀疏表示：原始信号可以通过一个过完备字典（Overcomplete Dictionary）以稀疏形式表示，即信号在字典基下的系数大部分为零或接近零。
欠采样：通过一个非适应性（Non-adaptive）测量矩阵对信号进行远低于Nyquist采样率的采样，得到观测值（或测量值）。
稀疏重建：利用优化算法从观测值出发，寻找与观测值最为匹配的、在字典基下具有最小稀疏度的信号表示。这是一个通常具有非凸、非线性特性的优化问题。

Lasso先验与L1范数

在压缩感知的稀疏重建阶段，Lasso方法是一种常用的正则化回归技术，它引入了L1范数作为先验约束，鼓励解的稀疏性。Lasso回归的目标函数可以表示为：
$$[\arg \underset{\mathbf{\theta }}{\min }\frac{1}{2}||\mathbf{y}-\mathbf{A}\mathbf{\theta }|{|}_2^2+\alpha ||\mathbf{\theta }|{|}_1]$$
其中：$\mathbf{y}$ 是观测值向量。$\mathbf{A}$是测量矩阵。$\mathbf{\theta }$ 是待求的稀疏信号表示（或系数向量）。$||\cdot |{|}_2$表示L2范数（欧几里得范数），用于衡量模型拟合误差。$||\cdot |{|}_1$ 表示L1范数，对应Lasso中的L1正则化项，它对$\mathbf{\theta }$中的绝对值求和，鼓励系数向量的稀疏性。$\alpha$是正则化参数。
Lasso通过求解上述优化问题，得到一个稀疏的(\boldsymbol{\theta})，进而重构出原始信号。

基于L1（Lasso）先验的层析重建过程如下：

稀疏表示：将待重建图像视为一个稀疏信号，通常选择一个合适的过完备字典（如小波基、离散余弦变换基、稀疏基等），使其在该字典下具有稀疏表示。
欠采样：在实际扫描过程中，通过非均匀、非适应性采样方案（如Radon变换、傅立叶采样等）采集远低于Nyquist采样率的投影数据（即(\mathbf{y})）。
L1范数优化：利用Lasso回归或其他L1正则化优化算法（如坐标下降法、梯度投影法、ADMM等），以采集到的投影数据$\mathbf{y}$和测量矩阵$\mathbf{A}$为输入，求解上述Lasso目标函数，得到稀疏的图像系数$\mathbf{\theta }$。
图像重构：将得到的稀疏系数$\mathbf{\theta }$通过所选字典基进行逆变换，重构出原始图像。

使用np.ravel()将x转化为一维数组
计算整数索引 floor_x： 对于 x 中的每个元素，计算 (x - orig) / dx 后向下取整得到整数部分，即该元素最接近的较小整数索引。使用 np.floor() 函数完成取整操作，并将结果转换为 np.int64 类型以确保整数精度。
计算权重因子 alpha： 计算 (x - orig) / dx 后的小数部分，即元素距离其整数索引的相对位置（在 [0, 1) 区间内）。这个小数部分代表了该元素与相邻整数索引之间的权重分配比例。

_generate_center_coordinates 生成一个以原点为中心、边长为 l_x 的正方形网格中每个单元格的中心坐标。

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

build_projection_operator(l_x, n_dir) 用于计算断层成像的设计矩阵。设计矩阵描述了图像的线积分（投影）与其自身的对应关系，在断层成像场景中起到核心作用，用于从不同角度获取的投影数据中重建图像。

函数步骤：

1. 生成网格坐标： 调用 _generate_center_coordinates(l_x) 函数，返回边长为 l_x 的正方形，图像中每个像素的中心坐标 X 和 Y。
2. 定义投影角度： 创建包含 n_dir 个均匀分布于 0 到 π 弧度（不包括π）之间的角度数组 angles。这些角度代表了采集投影的不同方向。
3. 初始化列表： 设置空列表 data_inds、weights 和 camera_inds，用于存储设计矩阵中非零条目的索引和权重。
4. 创建数据展开索引： 生成正方形图像中所有像素的扁平索引序列 data_unravel_indices。将该序列复制一次，以考虑到每个像素对两个互补角度的投影均有贡献。
5. 遍历角度： 对于 angles 中的每个角度 angle：
   旋转坐标： 通过变换 Xrot = np.cos(angle) * X - np.sin(angle) * Y 将像素中心坐标 X 和 Y 旋转至当前角度。模拟图像沿给定角度方向投射至一条垂直线的过程。
   计算权重和索引： 调用 _weights 函数，传入旋转后的坐标 Xrot 和固定步长 dx=1（假设像素宽度为1）。该函数返回一个包含两个数组的元组：inds（投影线上的像素中心索引）和 w（对应的权重）。orig 参数设为 X.min() 以确保正确索引。
   应用掩码并追加至列表： 应用逻辑掩码，仅保留索引在范围 [0, l_x) 内的条目及其对应的权重。将掩码后的 inds（加上 i * l_x 以考虑设计矩阵的不同行）、w 和 data_unravel_indices 分别追加到各自列表。
6. 构造稀疏矩阵： 使用收集到的 weights、camera_inds 和 data_inds 列表，创建一个稀疏矩阵的坐标列表（COO）格式，使用 sparse.coo_matrix。这代表了断层成像设计矩阵 proj_operator，其中每个非零条目表示特定像素对特定投影测量的贡献。

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

使用 scipy.ndimage.gaussian_filter 对二值掩模应用高斯滤波。滤波器核的标准差 sigma 设为 l / n_pts。这会模糊前景点的边缘，使得前景与背景之间过渡更平滑。
创建一个新的二值掩模 res，通过执行两个条件的逻辑与运算：
mask > mask.mean()：模糊掩模中值大于掩模均值的元素（本质上是一种阈值操作）。
mask_outer：之前定义的外侧圆形掩模，确保只考虑圆内的点。

返回通过逻辑异或运算得到的最终合成二值图像。运算的双方分别是 res 和 res 的腐蚀版本。ndimage.binary_erosion 通过移除边界像素减小前景区域的大小，实质上缩小了前景对象。与腐蚀图像进行异或运算引入了额外的复杂性，创造出更为复杂的合成图像。

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，其中 True 值表示圆内的点，False 值表示圆外的点'''
    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))

使用线性回归模型（Ridge 和 Lasso）对合成数据进行断层成像重建

调用 build_projection_operator 函数，生成一个投影算子（proj_operator），用于从图像（数据）中获取特定角度的投影。使用 l 作为图像边长，l // 7 作为投影角度的数量（约等于 18.3 度间隔）。
用 generate_synthetic_data 函数，生成一个大小为 l × l 的二值合成图像数据（data）。
proj = proj_operator * data.ravel()[:, np.newaxis]：将 data 展平为一维数组，然后增加一维以适应矩阵乘法，最后与投影算子 proj_operator 相乘，得到断层成像的投影数据（proj）。
proj += 0.15 * np.random.randn(*proj.shape)：向投影数据中添加高斯噪声。乘以 0.15 控制噪声幅度，np.random.randn(*proj.shape) 生成与 proj 形状相同的随机高斯分布噪声。

使用L2正则化重建

创建一个 Ridge 回归模型实例，设置正则化参数 alpha 为 0.2，表示使用 L2 范数惩罚项。
rgr_ridge.fit(proj_operator, proj.ravel())：使用投影算子 proj_operator 和投影数据 proj.ravel() 训练 Ridge 回归模型。
rec_l2 = rgr_ridge.coef_.reshape(l, l)：将模型系数（即重建的图像）重塑为 l × l 的二维数组，得到使用 L2 正则化的重建图像 rec_l2。

使用 L1（Lasso）正则化进行重建

# 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()