Multi-dimensional scipy.optimize.LinearConstraint?

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.