scipy.optimize Constrained Minimization SLSQP - Cannot Match the target 100%

Question

I am trying to get the minimized weight that's closer to average weight combined. My current problem is using the SLSQP solver I cannot find the right weight that meet the target 100%. Is there another solver I could use to solve my problem? Or any math suggestions. Please help.

My math right now is

           **min⁡(∑|x-mean(x)|)**

           **s.t.** Aw-b=0, w>=0


           **bound** 0.2'<'x<0.5

Data

  nRow# Variable1   Variable2   Variable3   

  1 3582.00 233445193.00    559090945.00    

  2 3394.00 217344811.00    496500751.00    

  3 3356.00 237746918.00    493639029.00    

  4 3256.00 219204892.00    461547877.00    

  5 3415.00 225272825.00    501057960.00    

  6 3505.00 242819442.00    505073223.00    

  7 3442.00 215258725.00    490458632.00    

  8 3381.00 227681178.00    503102998.00    

  9 3392.00 215189377.00    487026744.00    

    w1          w2                     w3
  Target    8531.00 429386951.00    1079115532.00

Question: Find the minimized weight that are closer to average

Python 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]) 

n=9
def fsolveMin(A,b,n):

    # The constraints: Ax = b

    def cons(x):
        return A.dot(x)-b

    cons = ({'type':'eq','fun':cons},

            {'type':'ineq','fun':lambda x:x[0]})

    # The minimizing constraints: the total absolute difference between the coefficients

    def fn(x):
        return np.sum(np.abs(x-np.mean(x)))       

    # Initialize the coefficients randomly
    z0 = abs(np.random.randn(len(A[1,:])))

    # Set up bound
   # bnds = [(0, None)]*n

    # Solve the problem

    sol = minimize(fn, x0 = z0, constraints = cons, method = 'SLSQP', options={'disp': True})



    #expected 35%
    print(sol.x)
    print(A.dot(sol.x))

    #print(fn(sol.x))

print(str(fsolveMin(A,b,n))+"

")

sascha · Accepted Answer

To give you an idea how bloated this will get with low-level tools like scipy's linprog as we have to mimic their standard-form:

The basic idea is to:

add an auxiliary-variable for the mean, like Erwin proposed
add n-auxiliary-variables to handle the absolute values like documented in lpsolve's docs
(and i'm using n-extra-aux-vars as temp-vars for x-mean)

Now this example works.

For your data you have to be careful in regards to the bounds. Make sure, the problem stays feasible. The problem itself is pretty unstable as hard-constraints are used for Ax=b; usually you would minimize some norm / least-squares here (no LP anymore; QP/SOCP) and add this error to the objective)!

It might be needed to switch the solver from method='simplex' to method='interior-point' at some point (only available since scipy 1.0).

Alternative:

Using cvxpy, formulation is much easier (both variants mentioned) and you get quite a good solver (ECOS) for free.

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)

Output:

Original A
[[  4.17022005e-01   7.20324493e-01   1.14374817e-04]
 [  3.02332573e-01   1.46755891e-01   9.23385948e-02]]
hidden x:  [ 0.18626021  0.34556073  0.39676747]
target:  [ 0.32663584  0.14366255]
Reformulation
y = mean A:
[ 0.33333333  0.33333333  0.33333333 -1.          0.          0.          0.          0.          0.          0.        ]
y = mean b:
[0]
Ax=b A:
[[  4.17022005e-01   7.20324493e-01   1.14374817e-04   0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00
    0.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  3.02332573e-01   1.46755891e-01   9.23385948e-02   0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00
    0.00000000e+00   0.00000000e+00   0.00000000e+00]]
Ax=b b:
[ 0.32663584  0.14366255]
z = x - mean A:
[[-1. -0. -0.  1.  1.  0.  0.  0.  0.  0.]
 [-0. -1. -0.  1.  0.  1.  0.  0.  0.  0.]
 [-0. -0. -1.  1.  0.  0.  1.  0.  0.  0.]]
z = x - mean b:
[ 0.  0.  0.]
solve...
     fun: 0.078779576294411263
 message: 'Optimization terminated successfully.'
     nit: 10
   slack: array([ 0.07877958,  0.        ,  0.        ,  0.        ,  0.07877958,  0.        ,  0.76273076,  0.68395118,
        0.72334097])
  status: 0
 success: True
       x: array([ 0.23726924,  0.31604882,  0.27665903,  0.27665903, -0.03938979,  0.03938979,  0.        ,  0.03938979,
        0.03938979,  0.        ])
x:  [ 0.23726924  0.31604882  0.27665903]
residual Ax-b:  [  5.55111512e-17   0.00000000e+00]
mean:  0.276659030053
x - mean:  [-0.03938979  0.03938979  0.        ]
l1-norm(x - mean) / objective:  0.0787795762944

Now for your original data, things get tough!

You need to:

scale your data as the magnitude of your variables will hurt the solver
- depending on your task you might need to analyze the effects of scaling / invert it
use the interior-point solver
be careful with the bounds

Example:

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')

Output:

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

Edit

Here a SOCP-based least-squares (soft-constrained) approach, which i would recommend in terms of numerical-stability! This approach can and should be tuned for whatever you need. It's implemented using the already mentioned cvxpy modelling-tool using the ECOS solver.

The basic idea:

instead of: min(l1-norm(x - mean(x)) st. Ax=b
solve: min(l2-norm(Ax-b) + c * l1-norm(x - mean(x)))
- where c is the nonnegative trade-off parameter
- small c: Ax=b is more important
- big c: x - mean(x) is more important

Example code & output for your data and bounds of -50, 50 and c=1e-3:

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)

Output:

 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]]

This approach will always output a feasible solution (for our task) and you can then decide on the observed residual if it works for you.

As you observed, a lower-bound of 0 is deadly, in all formulations (look at the magnitude difference in your data!).

Here a lower-bound of 0 will get you a solution with some high residual-error.

E.g.:

c=1e-7
bounds = 0 / 15

Output:

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]]

scipy.optimize Constrained Minimization SLSQP - Cannot Match the target 100%

Answers (2)

Related Questions