基于Scipy minimize的组合权重约束优化问题咨询

问题核心 #

你要解决的是**混合整数二次规划（MIQP）**问题：在最小化二次拟合误差（( |A\mathbf{w} - \mathbf{b}|_2^2 )）的前提下，要求每个权重要么为0，要么绝对值不小于阈值( T )（即( |w_i| \geq T )或( w_i=0 )）。这种非连续约束无法用传统连续优化工具直接处理，必须转化为整数+连续变量的组合优化问题。

原启发式方案的不足 #

你提到的“截断小权重后重新优化”是一种简单近似，但存在明显缺陷：

• 每次截断会缩小可行域，容易陷入局部最优，无法保证全局最优解。
• 迭代过程依赖手动调整，效率低且结果不稳定。

严谨的MIQP建模与实现 #

通过引入二进制变量将非连续约束线性化：
对每个权重( w_i )，定义两个二进制变量( b_i^+ )和( b_i^- )（均属于( {0,1} )），分别标记权重为正非零、负非零的状态，同时添加以下约束：

1. 若( b_i^+=1 )，则( w_i \geq T )；若( b_i^-=1 )，则( w_i \leq -T )；若两者都为0，则( w_i=0 )。
2. ( b_i^+ + b_i^- \leq 1 )（同一权重不能同时为正非零和负非零）。

用大M法将这些约束转化为线性不等式（M取权重可能的最大绝对值，比如100），即可用MIQP求解器处理。

代码实现（基于cvxpy） #

cvxpy支持整数变量和二次目标，是实现这类问题的理想工具：

import cvxpy as cp
import numpy as np
import pandas as pd

class TestMinimizer:
    def __init__(self, X, target, min_size=None):
        self.X = X
        self.target = target
        self.min_size = min_size  # 对应阈值T

    def solve(self) -> pd.Series:
        n = self.target.shape[0]
        # 定义变量：连续权重w，二进制变量标记正负非零状态
        w = cp.Variable(n)
        b_plus = cp.Variable(n, boolean=True)
        b_minus = cp.Variable(n, boolean=True)
        
        # 目标函数：最小化二次拟合误差
        objective = cp.Minimize(cp.sum_squares(self.X @ w - self.target))
        
        constraints = []
        if self.min_size is not None:
            T = self.min_size
            M = 100  # 根据实际场景调整，设为权重的最大可能绝对值
            for i in range(n):
                # 正非零约束：b_plus[i]=1时，w[i] >= T
                constraints.append(w[i] >= T - M * (1 - b_plus[i]))
                constraints.append(w[i] <= M * b_plus[i])
                # 负非零约束：b_minus[i]=1时，w[i] <= -T
                constraints.append(w[i] <= -T + M * (1 - b_minus[i]))
                constraints.append(w[i] >= -M * b_minus[i])
                # 互斥约束：同一权重不能同时为正负非零
                constraints.append(b_plus[i] + b_minus[i] <= 1)
        
        # 求解MIQP问题，需安装支持整数规划的求解器（如GUROBI、SCIP）
        prob = cp.Problem(objective, constraints)
        prob.solve(solver=cp.GUROBI)  # 若用SCIP，替换为cp.SCIP
        
        # 返回结果
        index = self.target.index if hasattr(self.target, 'index') else range(n)
        return pd.Series(w.value, index=index)

求解器选择 #

• 商业求解器如GUROBI、CPLEX在MIQP问题上速度快、精度高，提供学术免费许可证。
• 开源求解器如SCIP也可通过cvxpy调用，适合无商业授权的场景，但性能略逊。

启发式替代方案（快速近似） #

如果无法使用MIQP求解器，可迭代优化截断小权重：

1. 先求解无约束最小二乘得到初始权重。
2. 将绝对值小于T的权重强制设为0。
3. 固定这些0权重，重新优化剩余权重。
4. 重复步骤2-3直到没有权重被截断。

代码示例：

def solve_heuristic(self) -> pd.Series:
    # 初始无约束解
    w = np.linalg.lstsq(self.X, self.target, rcond=None)[0]
    T = self.min_size
    if T is None:
        return pd.Series(w, index=self.target.index if hasattr(self.target, 'index') else range(len(w)))
    
    changed = True
    while changed:
        # 筛选出需要保留的权重（绝对值>=T）
        mask = np.abs(w) >= T
        if np.all(mask):
            changed = False
            break
        # 固定小权重为0，优化剩余权重
        free_idx = np.where(mask)[0]
        X_free = self.X[:, free_idx]
        # 计算剩余权重的最优解
        w_free = np.linalg.lstsq(X_free, self.target - self.X[:, ~mask] @ w[~mask], rcond=None)[0]
        w[free_idx] = w_free
        w[~mask] = 0
    
    index = self.target.index if hasattr(self.target, 'index') else range(len(w))
    return pd.Series(w, index=index)

总结 #

• 追求全局最优解时，优先选择MIQP方法，搭配专业整数规划求解器。
• 启发式方法实现简单、速度快，但仅能得到近似解，适合对精度要求不高的场景。

内容的提问来源于stack exchange，提问作者Branck