如何用Python求解大规模0-1整数优化的最小二乘问题？

我有线性方程组 $Ax=b$，其中：

• $A$ 为 $m×n$ 稀疏矩阵，非零元素约 $n*m^{1/2}$ 个
• $x$ 为 $n×1$ 向量，$b$ 为 $m×1$ 向量
• 规模可达 $n=28680、m=22500$ 或 $62500$

该系统大概率无精确解，需找到仅含0和1的向量 $x_0$，使 $\sum((Ax-b)^2)$ 最小，即求解：
$$\min \sum((Ax-b)^2) \quad s.t. \quad x \in {0,1}^n$$

已尝试的解法 #

• 先求解实数域下的简化问题 $\sum((Ax-b)^2) \to \min$，用TensorFlow快速得到解，但取整后效果极差；
• 使用CVXPY库调用SCIP混合整数二次规划求解器，限制求解时间200秒，仅能处理小规模（如 $n=4950、m=1600$），规模扩大后SCIP返回零向量，但测试发现单个1的向量都比零向量更接近最优解。相关代码如下：

import cvxpy as cp
import numpy as np

# 模拟输入，在m=1600、n=4950时已出现问题
np.random.seed(0)
m, n= 1600, 4950
A = np.round(np.random.rand(m, n)/(1-1/m**(1/2))/2)*255*(m**(1/2)/1800)**1.1
b = 255*(np.random.randn(m))

x = cp.Variable(n, integer=True)
objective = cp.Minimize(cp.sum_squares(A @ x - b))
constraints = []
constraints.extend([x >= 0.0, x <= 1])
prob = cp.Problem(objective,constraints)
prob.solve(solver=cp.SCIP,verbose=True,scip_params={"limits/time": 50})
print("Status: ", prob.status)
print("The optimal value is", prob.value)
print("A solution x is")
print(x.value)
print("Non-zero ",sum([i for i in x.value]))

更新尝试 #

原本想避免导数不连续的函数，但发现绝对误差和 $\sum(|Ax-b|)$ 在示例中表现优于平方和。我的真实矩阵A并非示例中的噪声矩阵，虽无精确解，但零向量绝非最优（测试含单个1的向量中约80%更优）。改用绝对误差和后，SCIP不再返回零向量，但解中1的数量过多，且零向量对应的 $\sum(|Ax-b|)$ 反而更小，对此存疑。

疑问 #

1. 是否有更优的解法？
2. 是否有其他免费求解器可选？
3. 若SCIP为最优选择，如何使其适配大规模m、n？为何规模扩大后失效？

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