如何在Python中建模求解Dantzig Selector对应的线性规划问题

前置依赖 #

使用cvxpy凸优化库实现，无需手动拆解复杂线性规划参数，安装命令：

pip install cvxpy numpy

完整可运行代码 #

在你现有预处理代码的基础上，补充线性规划求解逻辑即可：

import numpy as np
import cvxpy as cp

# 原有的预处理代码
n = 5
k = 23

y = np.array([[1],
              [0],
              [0],
              [1],
              [0]])
X = np.array([[1.1, 0, 0.7, 0.8, 0.9, 0.2, 0.3, 0.5, 0.2, 0.2, 1.2, 1.1, 0.5, 0.5, 0.7, 1.2, 1.3, 0.8, 0.9, 1.7, 1.2, 1.9, 0.9],
              [0.3, 0.1, 0.7, 0.4, 0.9, 0.9, 0.1, 0.8, 0.1, 0.2, 1.1, 0, 0.9, 0.4, 1.4, 1.4, 0.1, 0.5, 1.8, 1.6, 1.2, 1.8, 0.3],
              [0.1, 0.1, 0.3, 0.9, 0.7, 0.8, 0, 0.7, 0.8, 0.2, 1.1, 1.1, 0.5, 0.5, 0.8, 1.5, 0.2, 0.5, 1.6, 1.5, 1.2, 1.7, 0.5],
              [1.2, 0.2, 0.9, 0.8, 0.6, 0.2, 0.3, 0.5, 0.3, 0.2, 1.2, 1.1, 0.5, 0, 0.7, 1.2, 1.3, 0.8, 0.9, 1.7, 1.2, 1.9, 0.9],
              [0.2, 0.1, 0.6, 0, 0.5, 1.1, 0.2, 0.5, 0.9, 0.2, 1.2, 1.1, 0.8, 1.6, 0.5, 1.3, 0.2, 0.5, 1.7, 1.2, 1.2, 1.9, 0.1]])


# 缺失值填充
X_row_minima = np.where(X>0,X,X.max()).min(1)
X[X==0] = X_row_minima/2

# X按行归一化
X = X/np.linalg.norm(X, ord=2, axis=1, keepdims=True)

# y标准化
y = y - np.mean(y) / np.std(y)

Xt = X.T
Xty = Xt @ y
d0 = max(abs(Xty))
delta_list = np.linspace(0, d0, 100)

# 新增：Dantzig Selector求解逻辑
def solve_dantzig_selector(X, y, delta):
    n, k = X.shape
    # 定义优化变量
    beta = cp.Variable((k, 1))
    # L1范数转化为线性规划的辅助变量
    u = cp.Variable((k, 1))
    # 目标函数：最小化L1范数
    objective = cp.Minimize(cp.sum(u))
    # 约束条件
    constraints = [
        beta >= -u,
        beta <= u,
        X.T @ (y - X @ beta) <= delta,
        X.T @ (y - X @ beta) >= -delta
    ]
    # 求解
    prob = cp.Problem(objective, constraints)
    prob.solve()
    return beta.value

# 遍历所有delta值求解
beta_results = []
for delta in delta_list:
    beta_hat = solve_dantzig_selector(X, y, delta)
    beta_results.append(beta_hat)

说明 #

• 上述代码直接对应Dantzig Selector的原始问题定义，无需手动拆解线性规划的标准形式，降低了实现门槛
• 求解器默认使用ECOS，是轻量高效的开源线性规划求解器，适配你的问题规模
• 最终beta_results列表存储了所有delta值对应的β估计结果，可根据需求选取合适的参数

内容的提问来源于stack exchange，提问作者Yuyu F.