<p>为了让您了解如何使用诸如scipy的<a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linprog.html#scipy.optimize.linprog" rel="nofollow noreferrer">linprog</a><strong>这样的低级工具来实现这一点,我们必须模仿它们的标准形式</strong>:</p>
<p>其基本思想是:</p>
<ul>
<li>对于一个辅助变量,建议的平均值</li>
<li>添加n个辅助变量来处理绝对值<a href="http://lpsolve.sourceforge.net/5.5/absolute.htm" rel="nofollow noreferrer">like documented in lpsolve's docs</a></li>
<li>(我使用n个额外的aux-vars作为<code>x-mean</code>的临时变量)</li>
</ul>
<p>现在这个例子起作用了。在</p>
<p>对于你的数据,你必须小心界限。确保这个问题是可行的。这个问题本身是非常不稳定的,因为硬约束用于<code>Ax=b</code>;通常您会最小化一些范数/最小二乘法(不再是LP;QP/SOCP),并将此错误添加到目标中)!在</p>
<p>可能需要在某个时刻将解算器从<code>method='simplex'</code>切换到{<cd4>}(仅在scipy1.0之后可用)。在</p>
<p>备选方案:</p>
<p>使用<a href="http://www.cvxpy.org/en/latest/index.html" rel="nofollow noreferrer">cvxpy</a>,公式化要容易得多(两个变体都提到过),而且您可以免费获得一个相当好的解算器(ECOS)。在</p>
<p>代码:</p>
<pre><code>import numpy as np
import scipy.optimize as spo
np.set_printoptions(linewidth=120)
np.random.seed(1)
""" Create random data """
M, N = 2, 3
A = np.random.random(size=(M,N))
x_hidden = np.random.random(size=N)
b = A.dot(x_hidden) # target
n = N
print('Original A')
print(A)
print('hidden x: ', x_hidden)
print('target: ', b)
""" Optimize as LP """
def solve(A, b, n):
print('Reformulation')
am, an = A.shape
# Introduce aux-vars
# 1: y = mean
# n: z = x - mean
# n: abs(z)
n_plus_aux_vars = 3*n + 1
# Equality constraint: y = mean
eq_mean_A = np.zeros(n_plus_aux_vars)
eq_mean_A[:n] = 1. / n
eq_mean_A[n] = -1.
eq_mean_b = np.array([0])
print('y = mean A:')
print(eq_mean_A)
print('y = mean b:')
print(eq_mean_b)
# Equality constraints: Ax = b
eq_A = np.zeros((am, n_plus_aux_vars))
eq_A[:, :n] = A[:, :n]
eq_b = np.copy(b)
print('Ax=b A:')
print(eq_A)
print('Ax=b b:')
print(eq_b)
# Equality constraints: z = x - mean
eq_mean_A_z = np.hstack([-np.eye(n), np.ones((n, 1)), + np.eye(n), np.zeros((n, n))])
eq_mean_b_z = np.zeros(n)
print('z = x - mean A:')
print(eq_mean_A_z)
print('z = x - mean b:')
print(eq_mean_b_z)
# Inequality constraints: absolute values -> x <= x' ; -x <= x'
ineq_abs_0_A = np.hstack([np.zeros((n, n)), np.zeros((n, 1)), np.eye(n), -np.eye(n)])
ineq_abs_0_b = np.zeros(n)
ineq_abs_1_A = np.hstack([np.zeros((n, n)), np.zeros((n, 1)), -np.eye(n), -np.eye(n)])
ineq_abs_1_b = np.zeros(n)
# Bounds
# REMARK: Do not touch anything besides the first bounds-row!
bounds = [(0., 1.) for i in range(n)] + \
[(None, None)] + \
[(None, None) for i in range(n)] + \
[(0, None) for i in range(n)]
# Objective
c = np.zeros(n_plus_aux_vars)
c[-n:] = 1
A_eq = np.vstack((eq_mean_A, eq_A, eq_mean_A_z))
b_eq = np.hstack([eq_mean_b, eq_b, eq_mean_b_z])
A_ineq = np.vstack((ineq_abs_0_A, ineq_abs_1_A))
b_ineq = np.hstack([ineq_abs_0_b, ineq_abs_1_b])
print('solve...')
result = spo.linprog(c, A_ineq, b_ineq, A_eq, b_eq, bounds=bounds, method='simplex')
print(result)
x = result.x[:n]
print('x: ', x)
print('residual Ax-b: ', A.dot(x) - b)
print('mean: ', result.x[n])
print('x - mean: ', x - result.x[n])
print('l1-norm(x - mean) / objective: ', np.linalg.norm(x - result.x[n], 1))
solve(A, b, n)
</code></pre>
<p>输出:</p>
^{pr2}$
<p>现在<strong>对于您的原始数据,</strong>,事情变得很棘手!在</p>
<p>您需要:</p>
<ul>
<li>缩放数据,因为变量的大小会损害解算器
<ul>
<li>根据您的任务,您可能需要分析缩放/反转的效果</li>
</ul></li>
<li>使用内部点解算器</li>
<li>小心界限</li>
</ul>
<p>示例:</p>
<pre><code>A = np.array([[3582.000000, 3394.000000, 3356.000000, 3256.000000, 3415.000000,
3505.000000, 3442.000000, 3381.000000, 3392.000000],
[233445193.000000, 217344811.000000, 237746918.000000,
219204892.000000, 225272825.000000, 242819442.000000,
215258725.000000, 227681178.000000, 215189377.000000],
[559090945.000000, 496500751.000000, 493639029.000000,
461547877.000000, 501057960.000000, 505073223.000000,
490458632.000000, 503102998.000000, 487026744.000000]])
b = np.array([8531., 1079115532., 429386951.])
A /= 10000. # scaling
b /= 10000. # scaling
bounds = [(-50., 50.) for i in range(n)] + \
...
result = spo.linprog(c, A_ineq, b_ineq, A_eq, b_eq, bounds=bounds, method='interior-point')
</code></pre>
<p>输出:</p>
<pre><code>solve...
con: array([ 3.98410760e-09, 1.18067724e-08, 8.12879938e-04, 1.75969041e-03, -3.93853838e-09, -3.96305566e-09,
-4.10043555e-09, -3.94957667e-09, -3.88362764e-09, -3.89452381e-09, -3.95134592e-09, -3.92182287e-09,
-3.85762178e-09])
fun: 52.742900697626389
message: 'Optimization terminated successfully.'
nit: 8
slack: array([ 5.13245265e+01, 1.89309145e-08, 1.83429094e-09, 4.28687782e-09, 1.03726911e-08, 2.77000474e-09,
1.41837413e+00, 6.75769654e-09, 8.65285462e-10, 2.78501844e-09, 3.09591539e-09, 5.27429006e+01,
1.30944103e-08, 5.32994799e-09, 3.15369669e-08, 2.51943821e-09, 7.54848797e-09, 3.22510447e-09])
status: 0
success: True
x: array([ -2.51938304e+01, 4.68432810e-01, 2.68398831e+01, 4.68432822e-01, 4.68432815e-01, 4.68432832e-01,
-2.40754247e-01, 4.68432818e-01, 4.68432819e-01, 4.68432822e-01, -2.56622633e+01, -7.91749954e-09,
2.63714503e+01, 4.40376624e-09, -2.52137156e-09, 1.43834811e-08, -7.09187065e-01, 3.95395716e-10,
1.17990950e-09, 2.56622633e+01, 1.10134149e-08, 2.63714503e+01, 8.69064406e-09, 7.85131955e-09,
1.71534858e-08, 7.09187068e-01, 7.15309226e-09, 2.04519496e-09])
x: [-25.19383044 0.46843281 26.83988313 0.46843282 0.46843282 0.46843283 -0.24075425 0.46843282 0.46843282]
residual Ax-b: [ -1.18067724e-08 -8.12879938e-04 -1.75969041e-03]
mean: 0.468432821891
x - mean: [ -2.56622633e+01 -1.18805552e-08 2.63714503e+01 4.54189575e-10 -6.40499920e-09 1.04889573e-08 -7.09187069e-01
-3.52642715e-09 -2.67771227e-09]
l1-norm(x - mean) / objective: 52.7429006758
</code></pre>
<p><strong>编辑</strong></p>
<p>这里有一个基于SOCP的最小二乘(软约束)方法,我建议在数值稳定性方面!这种方法可以而且应该根据您的需要进行调整。它是使用前面提到的使用ECOS solver的<a href="http://www.cvxpy.org/en/latest/" rel="nofollow noreferrer">cvxpy</a>建模工具实现的。在</p>
<p>基本思路:</p>
<ul>
<li><strong>而不是</strong>:<code>min(l1-norm(x - mean(x)) st. Ax=b</code></li>
<li><strong>求解</strong>:<code>min(l2-norm(Ax-b) + c * l1-norm(x - mean(x)))</code>
<ul>
<li>其中<code>c</code>是非负权衡参数</li>
<li>小<code>c</code>:<code>Ax=b</code>更重要</li>
<li>大<code>c</code>:<code>x - mean(x)</code>更重要</li>
</ul></li>
</ul>
<p>数据和<code>-50, 50</code>和<code>c=1e-3</code>边界的示例代码和输出:</p>
<pre><code>import numpy as np
import cvxpy as cvx
""" DATA """
A = np.array([[3582.000000, 3394.000000, 3356.000000, 3256.000000, 3415.000000,
3505.000000, 3442.000000, 3381.000000, 3392.000000],
[233445193.000000, 217344811.000000, 237746918.000000,
219204892.000000, 225272825.000000, 242819442.000000,
215258725.000000, 227681178.000000, 215189377.000000],
[559090945.000000, 496500751.000000, 493639029.000000,
461547877.000000, 501057960.000000, 505073223.000000,
490458632.000000, 503102998.000000, 487026744.000000]])
b = np.array([8531., 1079115532., 429386951.])
n = 9
# A /= 10000. scaling would be a good idea
# b /= 10000. """
""" SOCP-based least-squares approach """
def solve(A, b, n, c=1e-1):
x = cvx.Variable(n)
y = cvx.Variable(1) # mean
lower_bounds = np.zeros(n) - 50 # -50
upper_bounds = np.zeros(n) + 50 # 50
constraints = []
constraints.append(x >= lower_bounds)
constraints.append(x <= upper_bounds)
constraints.append(y == cvx.sum_entries(x) / n)
objective = cvx.Minimize(cvx.norm(A*x-b, 2) + c * cvx.norm(x - y, 1))
problem = cvx.Problem(objective, constraints)
problem.solve(solver=cvx.ECOS, verbose=True)
print('Objective: ', problem.value)
print('x: ', x.T.value)
print('mean: ', y.value)
print('Ax-b: ')
print((A*x - b).value)
print('x - mean: ', (x - y).T.value)
solve(A, b, n)
</code></pre>
<p>输出:</p>
<pre><code> ECOS 2.0.4 - (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 +2.637e-17 -1.550e+06 +7e+08 1e-01 2e-04 1e+00 2e+07 - - 1 1 - | - -
1 -8.613e+04 -1.014e+05 +8e+06 1e-03 2e-06 2e+03 2e+05 0.9890 1e-04 2 1 1 | 0 0
2 -1.287e+03 -1.464e+03 +1e+05 1e-05 9e-08 4e+01 3e+03 0.9872 1e-04 3 1 1 | 0 0
3 +1.794e+02 +1.900e+02 +2e+03 2e-07 1e-07 1e+01 5e+01 0.9890 5e-03 5 3 4 | 0 0
4 -1.388e+00 -6.826e-01 +1e+02 1e-08 7e-08 9e-01 3e+00 0.9458 4e-03 7 6 6 | 0 0
5 +5.491e+00 +5.683e+00 +1e+01 1e-09 8e-09 2e-01 3e-01 0.9617 6e-02 1 1 1 | 0 0
6 +6.480e+00 +6.505e+00 +1e+00 2e-10 5e-10 3e-02 4e-02 0.8928 2e-02 1 1 1 | 0 0
7 +6.746e+00 +6.746e+00 +2e-02 3e-12 5e-10 5e-04 6e-04 0.9890 5e-03 1 0 0 | 0 0
8 +6.759e+00 +6.759e+00 +3e-04 2e-12 2e-10 6e-06 7e-06 0.9890 1e-04 1 0 0 | 0 0
9 +6.759e+00 +6.759e+00 +3e-06 2e-13 2e-10 6e-08 8e-08 0.9890 1e-04 2 0 0 | 0 0
10 +6.758e+00 +6.758e+00 +3e-08 5e-14 2e-10 7e-10 9e-10 0.9890 1e-04 1 0 0 | 0 0
OPTIMAL (within feastol=2.0e-10, reltol=4.7e-09, abstol=3.2e-08).
Runtime: 0.002901 seconds.
Objective: 6.757722879805085
x: [[-18.09169736 -5.55768047 11.12130645 11.48355878 -1.13982006
12.4290884 -3.00165819 -1.05158589 -2.4468432 ]]
mean: 0.416074272576
Ax-b:
[[ 2.17051777e-03]
[ 1.90734863e-06]
[ -5.72204590e-06]]
x - mean: [[-18.50777164 -5.97375474 10.70523218 11.0674845 -1.55589434
12.01301413 -3.41773246 -1.46766016 -2.86291747]]
</code></pre>
<p>这种方法将始终输出一个可行的解决方案(对于我们的任务),然后您可以决定观察到的残差,如果它对您有用的话。在</p>
<p>正如您所观察到的,在所有公式中,0的下限都是致命的(看看您的数据之间的巨大差异!)。在</p>
<p>在这里,下限为0将得到一个具有高残余误差的解。在</p>
<p>例如:</p>
<ul>
<li><code>c=1e-7</code></li>
<li>边界=<code>0 / 15</code></li>
</ul>
<p>输出:</p>
<pre><code>Objective: 785913288.2410747
x: [[ -5.57966858e-08 -4.74997454e-08 1.56066068e+00 1.68021234e-07
-3.55602958e-08 1.75340641e-06 -4.69609562e-08 -3.10216680e-08
-4.39482554e-08]]
mean: 0.173406926909
Ax-b:
[[ -3.29341696e+03]
[ -7.08072860e+08]
[ 3.41016903e+08]]
x - mean: [[-0.17340698 -0.17340697 1.38725375 -0.17340676 -0.17340696 -0.17340517
-0.17340697 -0.17340696 -0.17340697]]
</code></pre>