问题描述
我正在尝试使用 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”等是什么意思?我正在尝试按照这些列中的值来了解求解过程中发生的情况。
另外,对于那些想要关注我正在使用的数据集的人,here 是一个包含我的 A 矩阵(尺寸为 278481x61)的 CSV 文件,here 是一个 CSV 文件拿着我的 b 向量(尺寸 278481x1)。提前致谢!
解决方法
一些说明:
介绍
问题:可重现的代码
- 问题设置很懒惰:您可以添加导入和 csv 读取,而不是让我们复制它...
- 由于不同的解析,现在很容易对不同的数据进行实验...
任务 + 方法
- 这看起来像是以下内容的组合:
- 通用数学优化求解器
- 机器学习规模数据
- = 经常达到极限!
求解器:期望
- SCS 是基于一阶的,与 ECOS 或其他高阶方法相比,它的鲁棒性预计较低
您的模型
观察
- 求解器:SCS(默认选项)
- 根据残差的进度失败/造成严重破坏(可能是由于没有数字问题)
- 我这边看起来更疯狂
- 根据残差的进度失败/造成严重破坏(可能是由于没有数字问题)
- 求解器:ECOS(默认选项)
- 失败(由于数字问题,最初不可行)
分析
-
仅通过增加迭代限制无法解决这些数值问题
- 对可行性容忍度和 co. 进行更复杂的调整。将需要!
改造/修复
我们能让这些求解器工作吗?我想是的。
不是最小化平方和,而是最小化l2-norm。这对于我们的解决方案来说是等价的,如果我们对这个值感兴趣,我们可能平方这个目标值!
这是由 this 促成的:
我们强烈鼓励的一个特别的重新表述是消除二次形式——即 sum_square、sum(square(.)) 或 quad_form 之类的函数——只要有可能使用范数来构造等效模型。我们的经验告诉我们,二次形式通常会给 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 相关的大型原始 feas 误差
可能可以通过调整来改进,但在这里使用不同的求解器可能会更好!这里没有做太多分析。
----------------------------------------------------------------------------
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
输出求解器=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.]