least-squares - CVXPY 在二次规划优化问题上返回不可行/不准确

Question

我正在尝试使用 CVXPY 来解决非负最小二乘问题（附加约束是解向量中的条目之和必须等于 1）。但是，当我使用 SCS 求解器在这个简单的二次程序上运行 CVXPY 时，我让求解器运行最多 100000 次迭代，最后我遇到了一个错误，说二次程序不可行/不准确。但是，我非常有信心二次规划的数学设置是正确的（因此问题的约束不应该是不可行的）。此外，当我在较小的 A 矩阵和较小的 b 向量（只有前几行）上运行该程序时，求解器运行良好。这是我的代码和当前错误输出的快照：

x = cp.Variable(61)
prob = cp.Problem(cp.Minimize(cp.sum_squares(A @ x - b)), [x >= 0, cp.sum(x) == 1])
print('The problem is QP: ', prob.is_qp())
result = prob.solve(solver = cp.SCS, verbose = True, max_iters = 100000)
print("status: ", prob.status)
print("optimal value: ", prob.value)
print("cvxpy solution:")
print(x.value, np.sum(np.array(x)))

下面是输出：

The problem is QP:  True
----------------------------------------------------------------------------
    SCS v2.1.1 - Splitting Conic Solver
    (c) Brendan O'Donoghue, Stanford University, 2012
----------------------------------------------------------------------------
Lin-sys: sparse-direct, nnz in A = 2446306
eps = 1.00e-04, alpha = 1.50, max_iters = 100000, normalize = 1, scale = 1.00
acceleration_lookback = 10, rho_x = 1.00e-03
Variables n = 62, constraints m = 278545
Cones:  primal zero / dual free vars: 1
    linear vars: 61
    soc vars: 278483, soc blks: 1
Setup time: 1.96e+00s
----------------------------------------------------------------------------
 Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
     0| 5.73e+18  1.10e+20  1.00e+00 -4.09e+22  1.55e+26  1.54e+26  1.47e-01 
   100| 8.00e+16  9.05e+18  1.79e-01  1.60e+25  2.29e+25  7.01e+24  8.82e+00 
   200| 6.83e+16  6.25e+17  8.10e-02  1.66e+25  1.95e+25  2.93e+24  1.81e+01 
   300| 6.62e+16  1.42e+18  1.00e-01  1.60e+25  1.96e+25  3.56e+24  2.61e+01 
   400| 9.49e+16  3.76e+18  1.98e-01  1.64e+25  2.45e+25  8.07e+24  3.39e+01 
   500| 6.24e+16  1.69e+19  3.12e-03  1.74e+25  1.75e+25  9.06e+22  4.20e+01 
   600| 8.47e+16  5.06e+18  1.42e-01  1.65e+25  2.20e+25  5.47e+24  5.04e+01 
   700| 8.57e+16  4.41e+18  1.55e-01  1.64e+25  2.25e+25  6.04e+24  5.79e+01 
   800| 6.03e+16  1.60e+18  1.30e-02  1.72e+25  1.77e+25  4.50e+23  6.51e+01 
   900| 7.35e+17  2.84e+20  3.21e-01  1.56e+25  3.04e+25  1.46e+25  7.21e+01 
  1000| 9.11e+16  1.38e+19  1.40e-01  1.68e+25  2.23e+25  5.46e+24  7.92e+01 
.........many more iterations in between omitted........
99500| 2.34e+15  1.56e+17  1.61e-01  3.30e+24  4.57e+24  1.27e+24  5.32e+03 
 99600| 1.46e+18  1.10e+21  9.12e-01  2.56e+24  5.55e+25  5.26e+25  5.33e+03 
 99700| 1.94e+15  1.41e+16  8.97e-02  3.40e+24  4.07e+24  6.71e+23  5.34e+03 
 99800| 3.31e+17  6.68e+20  5.62e-01  2.71e+24  9.67e+24  6.91e+24  5.35e+03 
 99900| 2.51e+15  9.96e+17  1.03e-01  3.41e+24  4.19e+24  7.81e+23  5.35e+03 
100000| 2.21e+15  3.79e+17  1.40e-01  3.29e+24  4.37e+24  1.07e+24  5.36e+03 
----------------------------------------------------------------------------
Status: Infeasible/Inaccurate
Hit max_iters, solution may be inaccurate
Timing: Solve time: 5.36e+03s
    Lin-sys: nnz in L factor: 2726743, avg solve time: 1.60e-02s
    Cones: avg projection time: 7.99e-04s
    Acceleration: avg step time: 2.72e-02s
----------------------------------------------------------------------------
Certificate of primal infeasibility:
dist(y, K*) = 0.0000e+00
|A'y|_2 * |b|_2 = 2.5778e-01
b'y = -1.0000
============================================================================
status:  infeasible_inaccurate
optimal value:  None
cvxpy solution:
None var59256

有谁知道为什么：(1) CVXPY 可能无法解决我的问题？这只是需要更多求解器迭代的问题吗？如果有，有多少？(2) 列名“pri res”、“dua res”、“rel gap”等是什么意思？我正在尝试按照这些列中的值来了解求解过程中发生的情况。

此外，对于那些想要跟随我正在使用的数据集的人，这里是一个包含我的 A 矩阵（维度 278481x61）的 CSV 文件，这是一个包含我的 b 向量（维度 278481x1）的 CSV 文件。提前致谢！

score 1 · Accepted Answer

一些备注：

介绍

问题：可重现的代码

问题设置很懒惰：您可以添加导入和 csv 读取，而不是让我们复制它......
- 由于解析不同，现在很容易对不同的数据进行实验......

任务+方法

这看起来像以下组合：
- 通用数学选择求解器
- 机器学习规模数据
= 经常达到极限！

解决者：期望

SCS 是基于一阶的，与 ECOS 或其他高阶方法相比，预计其鲁棒性较差

你的模型

观察

求解器：SCS（默认选项）
- 根据残差的进度失败/运行严重（可能由于没有数字问题）
  - 在我身边看起来更疯狂
求解器：ECOS（默认选项）
- 失败（由于数字问题，最初不可行）

分析

您将无法仅通过增加迭代限制来解决这些数值问题
- 更复杂的可行性公差和 co 调整。将需要！

改造/修复

我们能让这些求解器工作吗？我认同。

与其最小化平方和，不如让我们最小化l2 范数。这与我们的解决方案等价，如果我们对这个值感兴趣，我们可能会平方目标值！

这是出于以下动机：

我们强烈鼓励的一种特殊的重新表述是消除二次形式——即像 sum_square、sum(square(.)) 或 quad_form 这样的函数——只要可以使用 norm 构建等效模型。我们的经验告诉我们，二次型通常对 CVX 使用的基础求解器构成数值挑战。

我们承认这个建议违背了传统观点：二次形式是典型的光滑凸函数，而范数是非光滑的，因此难以处理。但是对于 CVX 使用的圆锥求解器，这种智慧完全是倒退的。这是最适合圆锥公式和解决方案的规范。二次形式是通过将它们转换为圆锥形式来处理的——实际上是使用规范！这种转换过程带来了一些有趣的扩展挑战。如果建模者可以消除执行这种转换的需要，那就更好了。

代码

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

A = pd.read_csv('A_matrix.csv').to_numpy()
b = pd.read_csv('b_vector.csv').to_numpy().ravel()

x = cp.Variable(61)
prob = cp.Problem(cp.Minimize(cp.norm(A @ x - b)), [x >= 0, cp.sum(x) == 1])
result = prob.solve(solver = cp.SCS, verbose = True)

print("optimal value: ", prob.value)
print("cvxpy solution:")
print(x.value, np.sum(x.value))

输出求解器=cp.SCS（慢 CPU）

有效的求解器状态 + 慢 + 解决方案看起来不够稳健 -> 对称地在 0 附近波动 -> 关于 x=>0 的原始特征误差较大

可能可以通过调整来改进，但在这里使用不同的求解器可能会更好！这里不做太多分析。

----------------------------------------------------------------------------
        SCS v2.1.2 - Splitting Conic Solver
        (c) Brendan O'Donoghue, Stanford University, 2012
----------------------------------------------------------------------------
Lin-sys: sparse-direct, nnz in A = 2446295
eps = 1.00e-04, alpha = 1.50, max_iters = 5000, normalize = 1, scale = 1.00
acceleration_lookback = 10, rho_x = 1.00e-03
Variables n = 62, constraints m = 278543
Cones:  primal zero / dual free vars: 1
        linear vars: 61
        soc vars: 278481, soc blks: 1
Setup time: 1.63e+00s
----------------------------------------------------------------------------
Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
    0| 9.14e+18  1.39e+20  1.00e+00 -5.16e+20  6.20e+23  6.04e+23  1.28e-01
  100| 8.63e-01  1.90e+02  7.79e-01  5.96e+02  4.81e+03  1.17e-14  1.04e+01
  200| 5.09e-02  3.50e+02  1.00e+00  6.20e+03 -1.16e+02  5.88e-15  2.08e+01
  300| 3.00e-01  3.71e+03  7.64e-01  9.62e+03  7.19e+04  4.05e-15  3.17e+01
  400| 5.19e-02  1.76e+02  1.91e-01  4.71e+03  6.94e+03  3.87e-15  4.25e+01
  500| 4.60e-02  2.66e+02  2.83e-01  5.70e+03  1.02e+04  6.48e-15  5.25e+01
  600| 5.13e-03  1.08e+02  1.24e-01  5.80e+03  7.44e+03  1.72e-14  6.23e+01
  700| 3.35e-03  6.81e+01  9.64e-02  5.39e+03  4.44e+03  5.94e-15  7.15e+01
  800| 1.62e-02  8.52e+01  1.17e-01  5.51e+03  6.97e+03  3.96e-15  8.07e+01
  900| 1.93e-02  1.57e+01  1.89e-02  5.58e+03  5.38e+03  5.04e-15  8.98e+01
  1000| 6.94e-03  6.85e+00  7.97e-03  5.57e+03  5.48e+03  1.75e-15  9.91e+01
  1100| 4.64e-03  7.64e+00  1.42e-02  5.66e+03  5.50e+03  1.91e-15  1.09e+02
  1200| 2.25e-04  3.25e-01  4.00e-04  5.61e+03  5.60e+03  5.33e-15  1.18e+02
  1300| 4.73e-05  9.05e-02  5.78e-05  5.60e+03  5.60e+03  6.16e-15  1.28e+02
  1400| 6.27e-07  4.58e-03  3.22e-06  5.60e+03  5.60e+03  7.17e-15  1.36e+02
  1500| 2.02e-07  5.27e-05  4.58e-08  5.60e+03  5.60e+03  5.61e-15  1.46e+02
----------------------------------------------------------------------------
Status: Solved
Timing: Solve time: 1.46e+02s
        Lin-sys: nnz in L factor: 2726730, avg solve time: 2.54e-02s
        Cones: avg projection time: 1.16e-03s
        Acceleration: avg step time: 5.61e-02s
----------------------------------------------------------------------------
Error metrics:
dist(s, K) = 7.7307e-12, dist(y, K*) = 0.0000e+00, s'y/|s||y| = 3.0820e-18
primal res: |Ax + s - b|_2 / (1 + |b|_2) = 2.0159e-07
dual res:   |A'y + c|_2 / (1 + |c|_2) = 5.2702e-05
rel gap:    |c'x + b'y| / (1 + |c'x| + |b'y|) = 4.5764e-08
----------------------------------------------------------------------------
c'x = 5602.9367, -b'y = 5602.9362
============================================================================
optimal value:  5602.936687635506
cvxpy solution:
[ 1.33143619e-06 -5.20173272e-07 -5.63980428e-08 -9.44340768e-08
  6.07765135e-07  7.55998810e-07  8.45038786e-07  2.65626921e-06
-1.35669263e-07 -4.88286704e-07 -1.09285233e-06  8.63799377e-07
  2.85145903e-07 -1.22240651e-06  2.14526505e-07 -2.40179173e-06
-1.75042884e-07 -1.27680170e-06 -1.40486649e-06 -1.12113037e-06
-2.26601198e-07  1.39878723e-07 -3.19396803e-06 -6.36480154e-07
  2.16005860e-05  1.18205616e-06  2.15620316e-06 -1.94093348e-07
-1.88356275e-06 -7.07687270e-06 -1.99902966e-06 -2.28894738e-06
  1.00000188e+00 -9.95601469e-07 -1.26333877e-06  1.26336565e-06
-5.31474195e-08 -9.81111443e-07  2.22755569e-07 -7.49418940e-07
-4.77882668e-07  6.89785690e-07 -2.46822613e-06 -5.73596077e-08
  5.99307819e-07 -2.57301316e-07 -7.59150986e-07 -1.23753681e-08
-1.39938273e-06  1.48526305e-06 -2.41075790e-06 -3.50224485e-07
  1.79214177e-08  6.71852182e-07 -5.10880844e-06  2.44821668e-07
-2.88655782e-06 -2.45457029e-07 -4.97712502e-07 -1.44497848e-06
-2.22294748e-07] 0.9999895863519757

输出solver=cp.ECOS（慢CPU）

有效的求解器状态 + 更快 + 解决方案看起来不错

ECOS 2.0.7 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS

It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT
0  +0.000e+00  -0.000e+00  +7e+04  9e-01  1e-04  1e+00  1e+03    ---    ---    1  2  - |  -  -
1  -5.108e+01  -4.292e+01  +7e+04  6e-01  1e-04  9e+00  1e+03  0.0218  3e-01   4  5  4 |  0  0
2  +2.187e+02  +2.387e+02  +5e+04  6e-01  8e-05  2e+01  8e+02  0.3427  4e-02   4  5  5 |  0  0
3  +1.109e+03  +1.118e+03  +1e+04  3e-01  2e-05  9e+00  2e+02  0.7403  6e-03   4  5  5 |  0  0
4  +1.873e+03  +1.888e+03  +1e+04  2e-01  2e-05  1e+01  2e+02  0.2332  1e-01   5  6  6 |  0  0
5  +3.534e+03  +3.565e+03  +4e+03  8e-02  8e-06  3e+01  7e+01  0.7060  2e-01   5  6  6 |  0  0
6  +5.452e+03  +5.453e+03  +2e+02  2e-03  2e-07  1e+00  3e+00  0.9752  2e-03   6  8  8 |  0  0
7  +5.584e+03  +5.585e+03  +4e+01  4e-04  4e-08  4e-01  7e-01  0.8069  6e-02   2  2  2 |  0  0
8  +5.602e+03  +5.602e+03  +5e+00  5e-05  6e-09  8e-02  9e-02  0.9250  5e-02   2  2  2 |  0  0
9  +5.603e+03  +5.603e+03  +1e-01  1e-06  1e-10  2e-03  2e-03  0.9798  2e-03   5  5  5 |  0  0
10  +5.603e+03  +5.603e+03  +6e-03  4e-07  6e-12  9e-05  1e-04  0.9498  3e-04   5  5  5 |  0  0
11  +5.603e+03  +5.603e+03  +4e-04  4e-07  3e-13  7e-06  6e-06  0.9890  4e-02   1  2  2 |  0  0
12  +5.603e+03  +5.603e+03  +1e-05  4e-08  8e-15  2e-07  2e-07  0.9816  8e-03   1  2  2 |  0  0
13  +5.603e+03  +5.603e+03  +2e-07  7e-10  1e-16  2e-09  2e-09  0.9890  1e-04   5  3  4 |  0  0

OPTIMAL (within feastol=7.0e-10, reltol=2.7e-11, abstol=1.5e-07).
Runtime: 18.727676 seconds.

optimal value:  5602.936884707248
cvxpy solution:
[7.47985848e-11 3.58238148e-11 4.53994101e-11 3.73056632e-11
3.47224797e-11 3.62895261e-11 3.59367993e-11 4.03642466e-11
3.58643375e-11 3.24886989e-11 3.25080912e-11 3.34866983e-11
3.66296670e-11 3.89612422e-11 3.54489152e-11 7.07301971e-11
3.95949165e-11 3.68235605e-11 3.05918372e-11 3.43890675e-11
3.71817538e-11 3.62561876e-11 3.55281653e-11 3.55800928e-11
4.10876077e-11 4.12877804e-11 4.11174782e-11 3.35519296e-11
3.43716575e-11 3.56588133e-11 3.66118962e-11 3.68789703e-11
9.99999998e-01 3.34857869e-11 3.21984616e-11 5.82577263e-11
2.85751155e-11 3.64710243e-11 3.59930485e-11 5.04742702e-11
3.07026084e-11 3.36507487e-11 4.19786324e-11 8.35032700e-11
3.33575857e-11 3.42732986e-11 3.70599423e-11 4.73856413e-11
3.39708564e-11 3.64354428e-11 2.95022064e-11 3.46315519e-11
3.04124702e-11 4.07870093e-11 3.57782184e-11 3.71824186e-11
3.72394185e-11 4.48194963e-11 4.09635820e-11 6.45638394e-11
4.00297122e-11] 0.9999999999673748

奥特罗

最后的评论

以上重新制定可能足以帮助您解决问题
ECOS 在这个大规模问题上击败 SCS 是不直观的，但总是会发生，我不会分析它（SCS 仍然是一个很好的解决方案，但 ECOS 也是，尽管两者是非常不同的方法！开源社区应该很高兴有那些！）
如果我有时间实现更多定制的东西，我认为我不会使用这些具有 ML 规模数据的求解器
- 这里想到的大规模求解的最简单方法：
  - 投影（加速）梯度（这将是一阶方法，对于您的约束而言非常强大！）
    - “投影到概率单纯形上”（您需要）是一件常见的（经过充分研究的）事情
根据最终的权重/系数：数据看起来很奇怪！
- 似乎有一个非常占主导地位的专栏（信息泄漏）；我不知道
- 四舍五入，两个求解器都将输出解向量：[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]