Why won't SciPy's minimize function work when linprog will?

Question

I have 3,787 nodes and 23,292 arcs. These nodes can be connected by the arcs. I was able to solve this problem with SciPy's linear programming solver, linprog, by using the dot products of utility and the solution as the objective function and stacking two matrices like this:

# nodes and arcs are already set

# Objective coefficients
coefficients = np.array(
    [pair_utility[(node1, node2)] for (node1, node2) in arcs]
)

A = # constraints matrix
b = # upper bounds here

# Bounds for variables
x_bounds = [(0, 1) for _ in arcs]

# Solve the linear programming problem
result = linprog(
    coefficients,
    A_ub=A,
    b_ub=b,
    bounds=x_bounds,
    method="highs",
    options={"disp": True},
)

This solution was lightning fast. However, I now need to add a new constraint which can only be represented as a function.

This means that I had to change the solver to SciPy's minimize function. Without even implementing the new constraint, I am not able to get a solution because the optimization tries to create a matrix with 5 billion plus elements, which my device cannot handle. Here is the code that creates the error:

x0 = [0 for _ in enumerate(arcs)]

# Objective coefficients
coefficients = np.array(
    [pair_utility[(trip1, trip2)] for (trip1, trip2) in arcs]
)
f = lambda x: coefficients.dot(x)

# Inequality constraints (Ax <= b)

A = # constraints matrix
b = # upper bound here
cons = ({"type": "ineq", "fun": lambda x: b - A.dot(x)},)

# Bounds for variables
x_bounds = [(0, 1) for _ in self.arcs]

# Solve the linear programming problem
result = minimize(f, x0, constraints=cons, bounds=x_bounds) 

The matrix of 5 billion elements is created inside the minimize function on the last line of the code snippet.

Am I implementing the new function wrong here? Can SciPy just not handle my large problem without a supercomputer and some time? Where should I go from here?

Here is a way to implement the problem with both the linprog and minimize function:

import uuid
from random import Random

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

rand = Random()

nodes = [uuid.uuid4().hex for _ in range(3787)]
temp = [tuple(rand.sample(range(1, 3787), 2)) for _ in range(23292)]
arcs = [(nodes[arc[0]], nodes[arc[1]]) for arc in temp]

pair_utility = {arc: -rand.random() for arc in arcs}

x0 = [1 for _ in enumerate(arcs)]

# Objective coefficients
coefficients = np.array([pair_utility[arc] for arc in arcs])
f = lambda x: coefficients.dot(x)

# Inequality constraints (Ax <= b)
A_incoming = np.zeros((len(nodes), len(arcs)))
A_outgoing = np.zeros((len(nodes), len(arcs)))

for i, node in enumerate(nodes):
    for j, (trip1_id, trip2_id) in enumerate(arcs):
        if node == trip1_id:
            A_outgoing[i, j] = 1
        elif node == trip2_id:
            A_incoming[i, j] = 1

A = np.vstack((A_incoming, A_outgoing))
b = np.ones(2 * len(nodes))
cons = ({"type": "ineq", "fun": lambda x: b - A.dot(x)},)

# Bounds for variables
x_bounds = [(0, 1) for _ in arcs]

# Solve the linear programming problem
result1 = linprog(
    coefficients,
    A_ub=A,
    b_ub=b,
    bounds=x_bounds,
    method="highs",
    options={"disp": True},
)
result2 = minimize(
    f,
    x0,
    constraints=cons,
    bounds=x_bounds,
    # method="L-BFGS-B",
    options={"disp": True},
)

And here is the error message:

Unable to allocate 38.3 GiB for an array with shape (5141509334,) and data type float64
    result2 = minimize(
              ^^^^^^^^^
numpy.core._exceptions._ArrayMemoryError: Unable to allocate 38.3 GiB for an array with shape (5141509334,) and data type float64```