My linear constraint for scipy.optimize.minimize is
ones = np.ones_like(x)
np.outer(x, ones) - np.outer(ones, x) > something
where something is a given matrix.
(Mathematically, a_ij < x_i - x_j < a_ji).
How do I express this in terms of scipy.optimize.LinearConstraint which wouldn't accept a 3d "matrix" as a 1st argument?
I'd suggest writing a loop that forms your constraint one row at a time.
If you want to add a constraint corresponding to a_ij < x_i - x_j < a_ji, you can do that by adding a vector to your constraint matrix which is all zeros, except for a one in the i'th position and a negative one in the j'th position. When SciPy takes the dot product of this vector with x, the result will be x_i - x_j.
You probably also want to avoid duplicate constraints. If you constrain both the value x_i - x_j and the value x_j - x_i, then those values are the same, except multiplied by -1. Having duplicate constraints can cause issues with some minimize methods due to low rank in the constraint matrix. In the below example, I have arbitrarily chosen that i < j for all constraints.
Here's an example.
import numpy as np
import scipy
something = np.arange(9).reshape(3, 3).T
constraints_list = []
lb = []
ub = []
N = 3
for i in range(N):
for j in range(N):
if i == j:
# If i == j, then x_i - x_j = 0, and constraint is always satisfied or always violated
continue
if j > i:
# If j > i, then this constraint duplicates another constraint.
continue
new_cons = np.zeros(N)
new_cons[i] = 1
new_cons[j] = -1
constraints_list.append(new_cons)
lb.append(something[i, j])
ub.append(something[j, i])
constraints_array = np.array(constraints_list)
constraint = scipy.optimize.LinearConstraint(constraints_array, lb, ub)
# Compute constraint values in old way
x = np.array([1.1, 1.2, 1.3])
print("Old constraint values")
ones = np.ones_like(x)
print(np.outer(x, ones) - np.outer(ones, x))
# Compute constraint values in new way
print("New constraint values (same as lower triangle of above array)")
print(constraint.A @ x)
print("Old constraint lower bound")
print(something)
print("New constraint lower bound (same as lower triangle of above array)")
print(constraint.lb)
Output:
Old constraint values
[[ 0. -0.1 -0.2]
[ 0.1 0. -0.1]
[ 0.2 0.1 0. ]]
New constraint values (same as lower triangle of above array)
[0.1 0.2 0.1]
Old constraint lower bound
[[0 3 6]
[1 4 7]
[2 5 8]]
New constraint lower bound (same as lower triangle of above array)
[1 2 5]
This essentially "flattens" your constraint into a 2D linear matrix.
itertool.combinations - but how would you avoid it altogether?You'll want your constraint to scale if the vector x is long. The proper way to do this is feed LinearConstraint sparse matrices:
import numpy as np
import scipy.sparse
def make_constraint_array(
n: int, # The number of decision variables
dtype: np.dtype = np.int8,
) -> scipy.sparse.sparray:
# The number of unique constraints
p = n*(n - 1)//2
# Sparse index type. Scipy requires this for certain calls.
itype = np.int32
# Index (column) pointers for the left-hand CSC array
decreasing = np.arange(n - 1, -1, -1, dtype=itype)
sums = decreasing.cumsum(dtype=itype)
indptr = np.zeros(n + 1, dtype=itype)
indptr[1:] = sums
# Left (minuend) terms: a column of n-1 ones, then n-2 ones, [...] then 0 ones.
# "standard" constructor form of https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csc_array.html
left = scipy.sparse.csc_array((
np.ones(shape=p, dtype=dtype), # data
np.arange(p, dtype=itype), # indices
indptr,
))
# Offsets (columns) for the right-hand COO array
repeat_offsets = np.repeat(a=n - sums[:-1], repeats=decreasing[:-1])
# Right (subtrahend) diagonal terms
# scipy.sparse.dia_array is one possible form, but half efficiency of the COO style.
# "existing data" constructor form of https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.coo_array.html
right = scipy.sparse.coo_array((
left.data,
(left.indices, left.indices + repeat_offsets), # y, x coords
))
return left - right
def demo() -> None:
print(make_constraint_array(n=6).toarray())
if __name__ == '__main__':
demo()
[[ 1 -1 0 0 0 0]
[ 1 0 -1 0 0 0]
[ 1 0 0 -1 0 0]
[ 1 0 0 0 -1 0]
[ 1 0 0 0 0 -1]
[ 0 1 -1 0 0 0]
[ 0 1 0 -1 0 0]
[ 0 1 0 0 -1 0]
[ 0 1 0 0 0 -1]
[ 0 0 1 -1 0 0]
[ 0 0 1 0 -1 0]
[ 0 0 1 0 0 -1]
[ 0 0 0 1 -1 0]
[ 0 0 0 1 0 -1]
[ 0 0 0 0 1 -1]]
onescan be simplified tox[:,None]-x[None,:](a broadcasted difference).x? If it's a two-dimensional matrix,np.outer(x, ones) - np.outer(ones, x)does not produce a 3D matrix, but a 2D matrix.np.subtract.outer(x.ravel(), x.ravel())or equivalent.xis a vector (cf.a_ij < x_i - x_j < a_ji)