How to use scipy `minimize` on a difference between two vectors?

Question

I have two vectors w1 and w2 (each of length 100), and I want to minimize the sum of their absolute difference i.e.

import numpy as np


def diff(w: np.ndarray) -> float:
    """Get the sum of absolute differences in the vector w.

    Args:
        w: A flattened vector of length 200, with the first 100 elements
            pertaining to w1, and the last 100 elements pertaining to w2.

    Returns:
        sum of absolute differences.
    """
    return np.sum(np.absolute(w[:100] - w[-100:]))

I need to write diff() as only taking one argument since scipy.opyimize.minimize requires the array passed to the x0 argument to be 1 dimensional.

As for constraints, I have

1. w1 is fixed and not allowed to change
2. w2 is allowed to change
3. The sum of absolute values w2 is between 0.1 and 1.1: 0.1 <= sum(abs(w2)) <= 1.1
4. |w2_i| < 0.01 for any element i in w2

I am confused as to how we code these constraints using the Bounds and LinearConstraints objects. What I've tried so far is the following

from scipy.optimize import minimize, Bounds, LinearConstraint

bounds = Bounds(lb=[-0.01] * 200, ub=[0.01] * 200)  # constraint #4
lc = LinearConstraint([[1] * 200], [0.1], [1.1])  # constraint #3
res = minimize(
    fun=diff,
    method='trust-constr',
    x0=w,  # my flattened vector containing w1 first 100 elements, and w2 in last 100 elements
    bounds=bounds,
    constraints=(lc)
)

My logic for the bounds variable is from constrain #4, and for the lc variable comes from constrain #3. However I know I've coded this wrong because because the lower and upper bounds are of length 200 which seems to indicate they are applied to both w1 and w2 whereas I only wan't to apply the constrains to w2 (I get an error ValueError: operands could not be broadcast together with shapes (200,) (100,) if I try to change the length of the array in Bounds from 200 to 100).

The shapes and argument types for LinearConstraint are especially confusing to me, but I did try to follow the scipy example.

This current implementation never seems to finish, it just hangs forever.

How do I properly implement bounds and LinearConstraint so that it satisfies my constraints list above, if that is even possible?

Answer 1

Your problem can easily be formulated as a linear optimization problem (LP). You only need to reformulate all absolute values of the optimization variables.

Changing the notation slightly (x is now the optimization variable w2 and w is just your given vector w1), your problem reads as

min |w_1 - x_1| + .... + |w_N - x_N|

s.t. lb <= |x1| + ... + |xN| <= ub         (3)
                       |x_i| <= 0.01 - eps (4) (models the strict inequality)                  

where eps is just a sufficiently small number in order to model the strict inequality.

Let's consider the constraint (3). Here, we add additional positive variables z and define z_i = |x_i|. Then, we replace all absolute values |x_i| by z_i and impose the constraints -x_i <= z_i <= x_i which model the relationship z_i = |x_i|. Similarly, you can proceed with the objective and the constraint (4). The latter is by the way trivial and equivalent to -(0.01 - eps) <= x_i <= 0.01 - eps.

In the end, your optimization problem should read (assuming that all your w_i are positive):

min u1 + .... + uN

s.t.      lb <= z1 + ... + zN <= ub
          -x <=      z        <= x
 -0.01 + eps <=      x        <= 0.01 - eps
    -(w-x)   <=      u        <= w - x
          0  <= z
          0  <= u

with 3*N optimization variables x1, ..., xN, u1, ..., uN, z1, ..., zN. It isn't hard to write these constraints as an matrix-vector product A_ineq * x <= b_ineq. Then, you can solve it by scipy.optimize.linprog as follows:

import numpy as np
from scipy.optimize import minimize, linprog

n = 100
w = np.abs(np.random.randn(n))
eps = 1e-10
lb = 0.1
ub = 1.1

# linear constraints: A_ub * (x, z, u)^T <= b_ub
A_ineq = np.block([
    [np.zeros(n), np.ones(n), np.zeros(n)],
    [np.zeros(n), -np.ones(n), np.zeros(n)],
    [-np.eye(n),  np.eye(n), np.zeros((n, n))],
    [-np.eye(n), -np.eye(n), np.zeros((n, n))],
    [ np.eye(n), np.zeros((n, n)), -np.eye(n)],
    [ np.eye(n), np.zeros((n, n)),  np.eye(n)],
])

b_ineq = np.hstack((ub, -lb, np.zeros(n), np.zeros(n), w, w))

# bounds: lower <= (x, z, u)^T <= upper 
lower = np.hstack(((-0.01 + eps) * np.ones(n), np.zeros(n), np.zeros(n)))
upper = np.hstack((( 0.01 - eps) * np.ones(n), np.inf*np.ones(n), np.inf*np.ones(n)))
bounds = [(l, u) for (l, u) in zip(lower, upper)]

# objective: c^T * (x, z, u)
c = np.hstack((np.zeros(n), np.zeros(n), np.ones(n)))

# solve the problem
res = linprog(c, A_ub=A_ineq, b_ub=b_ineq, method="highs")

# your solution
x = res.x[:n]
print(res.message)
print(x)

Some notes in arbitrary order:

• It's highly recommended to solve linear optimization problems with linprog instead of minimize. The former provides an interface to HiGHS, a high-performance LP solver HiGHs that outperforms all algorithms under the hood of minimize. However, it's also worth mentioning that minimize is meant to be used for nonlinear optimization problems.

• In case your values w are not all positive, we need to change the formulation.

Answer 2

You can (and perhaps should, for clarity), use the args argument in minimize, and provide the fixed vector as an extra argument to your function.

If you set up your equation as follows:

def diff(w2, w1):
    return np.sum(np.absolute(w1 - w2))

and your constraints with

bounds = Bounds(lb=[-0.01] * 100, ub=[0.01] * 100)  # constraint #4
lc = LinearConstraint([[1] * 100], [0.1], [1.1])  # constraint #3

and then do

res = minimize(
    fun=diff,
    method='trust-constr',
    x0=w1,
    args=(w2,),
    bounds=bounds,
    constraints=[lc]
)

Then:

print(res.success, res.status, res.nit, np.abs(res.x).sum(), all(np.abs(res.x) < 0.01))

yields (for me at least)

(True, 1, 17, 0.9841520351691752, True)

which seems to be what you want.

Note that my test inputs are:

w1 = (np.arange(100) - 50) / 1000
w2 = np.ones(100, dtype=float)

which may or may not be favourable to the fitting procedure.