## Wednesday, August 10, 2022

### Large sparse transportation model with CVXPY,CVXR

In [1] I was trying out different formulations of a large, sparse (but very easy) transportation model using different modeling tools and solvers. The conclusion was:

• Exploiting sparsity leads to the best performance
• GAMS/Cplex is about 10 times as fast as Python/PuLP/CBC
• GAMS/Cplex is about 1000 times as fast as R/OMPR/GLPK

Don't overinterpret these numbers: this is a single model, which may not be at all representable for your models. Let's see how CVXPY[2] and CVXR[3] are doing. To recap we want to solve a large (but easy) transportation model:

Dense Transportation Model
\begin{align}\min&\sum_{i,j}\color{darkblue}c_{i,j}\cdot\color{darkred}x_{i,j}\\ & \sum_j \color{darkred}x_{i,j} \le \color{darkblue}a_i && \forall i \\& \sum_i \color{darkred}x_{i,j} \ge \color{darkblue}b_j && \forall j \\ & \color{darkred}x_{i,j} \ge 0 \end{align}

The additional wrinkle is that only a small fraction of the links $$i \rightarrow j$$ exist.

There are (at least) three ways to model this:

1. Use a large cost coefficient for forbidden links.
2. Fix $$\color{darkred}x_{i,j}=0$$ for non-existing links.
3. Exploit sparsity directly by only summing over existing links.
CVXPY/CVXR are matrix-based modeling tools. That gives some extra things to think about. A transportation model in matrix-notation can look like:

Transportation Model in Matrix Notation
\begin{align}\min\>&{\bf tr}(\color{darkblue}C^T\color{darkred}X)\\&\color{darkred}X\color{darkblue}e \le \color{darkblue}a \\ & \color{darkred}X^t\color{darkblue}e \ge \color{darkblue}b \\ & \color{darkred}X \ge 0 \end{align}

The vectors $$\color{darkblue}e$$ are all ones, appropriately sized, and $$\bf tr$$ is the trace function. Not the most intuitive formulation, if you ask me. This model can be easily transcribed into a CVXPY Python-based model as follows:

#-------------------------------------------------
# dense transportation model
# matrix formulation
#-------------------------------------------------

# X has same shape as our cost matrix
X = cp.Variable(cost.shape,"X")

# two different e vectors
ei = np.ones(capacity.shape)
ej = np.ones(demand.shape)

prob = cp.Problem(
cp.Minimize(cp.trace(cost.T@X)),
[X.T@ei >= demand,
X@ej <= capacity,
X>=0])
prob.solve(verbose=True)
print("status:",prob.status)
print("objective:",prob.value)

Here @ indicates matrix multiplication. Arguably, a more readable model is using elementwise multiplication and summations:

#-------------------------------------------------
# dense transportation model
# sum formulation
#-------------------------------------------------

# X has same shape as our cost matrix
X = cp.Variable(cost.shape,"X")

prob = cp.Problem(
cp.Minimize(cp.sum(cp.multiply(cost,X))),
[cp.sum(X,axis=0) >= demand,
cp.sum(X,axis=1) <= capacity,
X >= 0])
prob.solve(verbose=True)
print("status:",prob.status)
print("objective:",prob.value)

Here cp.multiply() indicates elementwise multiplication.

The same model in R, using cvxr can look like:

X <- Variable(rows=n,cols=m,name='X',nonneg=T)

prob <- Problem(Minimize(sum_entries(cost*X)),
constraints = list(
sum_entries(X,axis=1) >= dem,
sum_entries(X,axis=2) <= cap))

res <- solve(prob,verbose=T)
res$status res$value

#### Solvers

I use the default solvers for this exercise, and default settings. For CVXPY the default solver is ECOS[4], and CVXR uses OSQP[5] by default. These are not standard Simplex-based solvers.

#### Formulation 1: Dense Transportation Model with Large Cost Coefficients.

Here we apply the costs:

where link is a 0/1 constant matrix with values: 0 if the link does not exist, and 1 if the link exists.

CVXPY uses the conic algorithm ECOS for solving LPs. This solver cannot solve the problem using default settings:

98  -4.660e+04  -4.660e+04  +5e+05  1e-04  7e-15  8e-02  2e+00  0.0232  6e-01   1  2  2 |  0  0
99  -4.655e+04  -4.655e+04  +4e+05  1e-04  9e-15  1e-01  2e+00  0.0272  7e-01   1  2  2 |  0  0
100  -4.654e+04  -4.654e+04  +4e+05  1e-04  1e-14  1e-01  2e+00  0.0407  7e-01   1  2  1 |  0  0
Maximum number of iterations reached, recovering best iterate (88) and stopping.

RAN OUT OF ITERATIONS (reached feastol=1.2e-05, reltol=4.7e-02, abstol=1.1e+03).
Runtime: 13.731298 seconds.

We can add a larger iteration limit. This gives:

257  +1.202e+04  +1.202e+04  +2e-01  3e-08  9e-17  4e-07  7e-07  0.9890  7e-02   1  1  1 |  0  0
258  +1.202e+04  +1.202e+04  +5e-03  1e-09  6e-17  1e-08  2e-08  0.9713  1e-04   2  1  1 |  0  0
259  +1.202e+04  +1.202e+04  +5e-05  1e-11  8e-17  1e-10  2e-10  0.9890  1e-04   2  1  1 |  0  0

OPTIMAL (within feastol=1.4e-11, reltol=4.4e-09, abstol=5.3e-05).
Runtime: 34.450007 seconds.

-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
(CVXPY) Aug 10 12:44:50 PM: Problem status: optimal
(CVXPY) Aug 10 12:44:50 PM: Optimal value: 1.202e+04
(CVXPY) Aug 10 12:44:50 PM: Compilation took 5.605e+00 seconds
(CVXPY) Aug 10 12:44:50 PM: Solver (including time spent in interface) took 3.446e+01 seconds
status: optimal
objective: 12015.21812055369

Similarly, CVXR using the default OSQP solver has troubles with this model:

9600   2.5164e+05   5.37e-04   6.25e+00   1.01e-01   6.13e+01s
9800   2.5184e+05   5.37e-04   5.97e+00   1.01e-01   6.25e+01s
10000   2.5167e+05   5.37e-04   5.98e+00   1.01e-01   6.38e+01s

status:               solved inaccurate
number of iterations: 10000
optimal objective:    251671.7459
run time:             6.38e+01s
optimal rho estimate: 1.53e-01

Warning messages:
1: In size(const) : NAs introduced by coercion to integer range
2: In size(const) : NAs introduced by coercion to integer range

I am not sure what the message about NAs is about. But in any case: this model was not solved to optimality. We can again increase the iteration limit. This gives only slightly better results with a reported objective of

optimal objective:    94120.8210
we are actually not close to the optimal solution.

This is a disappointing result.

#### Formulation 2: Dense Transportation Model with Fixed Variables.

The second formulation is using:

x[i,j] = 0 if link[i,j] = 0
I don't really know how to do this in a vectorized way in CVXPY/CVXR. Very strange that this is not supported, as this is an operation that we see quite often. We could form a long series of scalar equalities. I chose to formulate:
or in matrix notation: $\color{darkred}X \circ (\color{darkred}e\cdot\color{darkred}e^T-\color{darkblue}{\mathit link}) = 0$ where $$\circ$$ indicates the Hadamar or element-wise product. In Python notation this becomes: cp.multiply(1-link,X)==0.

When running this with CVXPY/ECOS we see success!

C:\tmp>py trnsport_cvx.py
CHECK: sum capacity = 10000.000000000016, sum demand = 10000.000000000002
===============================================================================
CVXPY
v1.2.1
===============================================================================
(CVXPY) Aug 10 01:16:27 AM: Your problem has 250000 variables, 4 constraints, and 0 parameters.
(CVXPY) Aug 10 01:16:27 AM: It is compliant with the following grammars: DCP, DQCP
(CVXPY) Aug 10 01:16:27 AM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)
(CVXPY) Aug 10 01:16:27 AM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
(CVXPY) Aug 10 01:16:27 AM: Compiling problem (target solver=ECOS).
(CVXPY) Aug 10 01:16:27 AM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> ECOS
(CVXPY) Aug 10 01:16:27 AM: Applying reduction Dcp2Cone
(CVXPY) Aug 10 01:16:27 AM: Applying reduction CvxAttr2Constr
(CVXPY) Aug 10 01:16:27 AM: Applying reduction ConeMatrixStuffing
(CVXPY) Aug 10 01:16:27 AM: Applying reduction ECOS
(CVXPY) Aug 10 01:16:27 AM: Finished problem compilation (took 4.100e-01 seconds).
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
(CVXPY) Aug 10 01:16:27 AM: Invoking solver ECOS  to obtain a solution.

ECOS 2.0.10 - (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  +5.471e+04  +5.471e+04  +2e+06  5e-01  5e-01  1e+00  1e+01    ---    ---    1  1  - |  -  -
1  +3.188e+04  +3.188e+04  +1e+06  3e-01  2e-01  2e+00  5e+00  0.6991  3e-01   1  1  1 |  0  0
2  +2.670e+04  +2.670e+04  +1e+06  2e-01  2e-01  3e+00  4e+00  0.3511  6e-01   1  1  0 |  0  0
3  +1.558e+04  +1.558e+04  +4e+05  8e-02  6e-02  2e+00  2e+00  0.9518  3e-01   1  1  1 |  0  0
4  +1.381e+04  +1.381e+04  +2e+05  3e-02  2e-02  8e-01  7e-01  0.7207  2e-01   1  1  0 |  0  0
5  +1.318e+04  +1.318e+04  +1e+05  2e-02  2e-02  4e-01  4e-01  0.6765  5e-01   1  1  0 |  0  0
6  +1.261e+04  +1.261e+04  +5e+04  8e-03  7e-03  2e-01  2e-01  0.6174  1e-01   1  1  1 |  0  0
7  +1.235e+04  +1.235e+04  +3e+04  4e-03  3e-03  9e-02  1e-01  0.6661  3e-01   1  1  0 |  0  0
8  +1.218e+04  +1.218e+04  +1e+04  2e-03  1e-03  3e-02  5e-02  0.7733  3e-01   1  1  1 |  0  0
9  +1.211e+04  +1.211e+04  +6e+03  8e-04  7e-04  2e-02  2e-02  0.7199  3e-01   1  1  1 |  0  0
10  +1.206e+04  +1.206e+04  +3e+03  4e-04  4e-04  8e-03  1e-02  0.7449  3e-01   1  1  1 |  0  0
11  +1.204e+04  +1.204e+04  +1e+03  2e-04  2e-04  3e-03  5e-03  0.6062  9e-02   1  1  1 |  0  0
12  +1.202e+04  +1.202e+04  +5e+02  7e-05  7e-05  1e-03  2e-03  0.6782  1e-01   1  1  1 |  0  0
13  +1.202e+04  +1.202e+04  +2e+02  2e-05  2e-05  4e-04  6e-04  0.7878  1e-01   1  1  1 |  0  0
14  +1.202e+04  +1.202e+04  +9e+01  1e-05  1e-05  2e-04  4e-04  0.4928  2e-01   2  1  1 |  0  0
15  +1.202e+04  +1.202e+04  +3e+01  3e-06  3e-06  6e-05  1e-04  0.8573  2e-01   2  1  1 |  0  0
16  +1.202e+04  +1.202e+04  +2e+01  2e-06  2e-06  4e-05  6e-05  0.6635  4e-01   3  1  2 |  0  0
17  +1.202e+04  +1.202e+04  +2e+00  2e-07  2e-07  5e-06  7e-06  0.8947  2e-02   4  1  2 |  0  0
18  +1.202e+04  +1.202e+04  +2e-01  3e-08  2e-08  5e-07  8e-07  0.9343  4e-02   6  2  3 |  0  0
19  +1.202e+04  +1.202e+04  +2e-03  3e-10  3e-10  6e-09  9e-09  0.9890  1e-03   1  1  1 |  0  0
20  +1.202e+04  +1.202e+04  +3e-05  3e-12  3e-12  6e-11  1e-10  0.9890  1e-04   2  0  0 |  0  0

OPTIMAL (within feastol=3.4e-12, reltol=2.2e-09, abstol=2.6e-05).
Runtime: 4.309761 seconds.

-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
(CVXPY) Aug 10 01:16:31 AM: Problem status: optimal
(CVXPY) Aug 10 01:16:31 AM: Optimal value: 1.202e+04
(CVXPY) Aug 10 01:16:31 AM: Compilation took 4.100e-01 seconds
(CVXPY) Aug 10 01:16:31 AM: Solver (including time spent in interface) took 4.322e+00 seconds
status: optimal
objective: 12015.21812552263

C:\tmp>

The first thing we notice is that counting is difficult: Your problem has 250000 variables, 4 constraints. We seem to count here individual variables and constraint blocks. I guarantee you there are more than 4 constraints. That sounds strange. I mean, why in the world would you do that? Choose the same units!

Now let's see if we also have some success with the CVXR/OSQP setup. Unfortunately, this still does not solve:

-----------------------------------------------------------------
OSQP v0.6.0  -  Operator Splitting QP Solver
(c) Bartolomeo Stellato,  Goran Banjac
University of Oxford  -  Stanford University 2019
-----------------------------------------------------------------
problem:  variables n = 250000, constraints m = 501000
nnz(P) + nnz(A) = 950334
settings: linear system solver = qdldl,
eps_abs = 1.0e-05, eps_rel = 1.0e-05,
eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
sigma = 1.00e-06, alpha = 1.60, max_iter = 10000
check_termination: on (interval 25),
scaling: on, scaled_termination: off
warm start: on, polish: on, time_limit: off

iter  objective    pri res    dua res    rho        time
1  -5.4672e+05   4.33e+01   4.22e+01   1.00e-01   8.30e-01s
200   1.2118e+04   5.96e-02   1.12e-01   1.00e-01   3.08e+00s
400   1.2047e+04   1.95e-02   9.76e-02   1.00e-01   5.46e+00s
600   1.2031e+04   7.57e-03   3.55e-02   1.00e-01   7.73e+00s
800   1.2043e+04   4.77e-03   2.27e-02   1.00e-01   9.95e+00s
1000   1.2023e+04   1.63e-02   1.40e-02   1.87e-02   1.23e+01s
1200   1.2055e+04   7.39e-03   8.88e-03   1.87e-02   1.45e+01s
1400   1.2034e+04   6.94e-03   7.34e-03   1.87e-02   1.67e+01s
1600   1.2048e+04   3.92e-03   4.51e-03   1.87e-02   1.89e+01s
1800   1.2047e+04   3.34e-03   4.63e-03   1.87e-02   2.12e+01s
2000   1.2038e+04   3.35e-03   2.90e-03   1.87e-02   2.34e+01s
2200   1.2041e+04   2.34e-03   2.50e-03   1.87e-02   2.56e+01s
2400   1.2040e+04   2.32e-03   2.47e-03   1.87e-02   2.78e+01s
2600   1.2046e+04   1.75e-03   1.82e-03   1.87e-02   3.00e+01s
2800   1.2044e+04   1.86e-03   1.45e-03   1.87e-02   3.22e+01s
3000   1.2044e+04   1.39e-03   1.47e-03   1.87e-02   3.45e+01s
3200   1.2041e+04   1.33e-03   1.34e-03   1.87e-02   3.67e+01s
3400   1.2045e+04   1.02e-03   1.62e-03   1.87e-02   3.89e+01s
3600   1.2043e+04   9.00e-04   1.20e-03   1.87e-02   4.12e+01s
3800   1.2043e+04   1.04e-03   1.12e-03   1.87e-02   4.34e+01s
4000   1.2040e+04   1.32e-03   1.17e-03   1.87e-02   4.56e+01s
4200   1.2044e+04   7.46e-04   9.03e-04   1.87e-02   4.78e+01s
4400   1.2043e+04   6.78e-04   8.80e-04   1.87e-02   5.00e+01s
4600   1.2042e+04   6.31e-04   1.19e-03   1.87e-02   5.22e+01s
4800   1.2042e+04   7.80e-04   8.06e-04   1.87e-02   5.44e+01s
5000   1.2042e+04   6.57e-04   8.05e-04   1.87e-02   5.67e+01s
5200   1.2044e+04   6.17e-04   9.06e-04   1.87e-02   5.89e+01s
5400   1.2042e+04   6.12e-04   6.06e-04   1.87e-02   6.11e+01s
5600   1.2041e+04   6.23e-04   6.60e-04   1.87e-02   6.33e+01s
5800   1.2041e+04   3.80e-04   7.29e-04   1.87e-02   6.55e+01s
6000   1.2043e+04   5.39e-04   6.00e-04   1.87e-02   6.78e+01s
6200   1.2043e+04   4.63e-04   6.70e-04   1.87e-02   7.00e+01s
6400   1.2042e+04   4.04e-04   6.56e-04   1.87e-02   7.22e+01s
6600   1.2041e+04   5.45e-04   4.48e-04   1.87e-02   7.44e+01s
6800   1.2042e+04   6.47e-04   5.94e-04   1.87e-02   7.66e+01s
7000   1.2043e+04   3.81e-04   5.65e-04   1.87e-02   7.88e+01s
7200   1.2043e+04   7.66e-04   4.04e-04   1.87e-02   8.11e+01s
7400   1.2042e+04   2.75e-04   5.74e-04   1.87e-02   8.33e+01s
7600   1.2043e+04   5.49e-04   4.90e-04   1.87e-02   8.55e+01s
7800   1.2042e+04   3.88e-04   3.64e-04   1.87e-02   8.77e+01s
8000   1.2042e+04   3.03e-04   5.45e-04   1.87e-02   8.99e+01s
8200   1.2043e+04   3.03e-04   5.13e-04   1.87e-02   9.21e+01s
8400   1.2043e+04   3.23e-04   3.85e-04   1.87e-02   9.43e+01s
8600   1.2042e+04   3.14e-04   4.93e-04   1.87e-02   9.65e+01s
8800   1.2043e+04   2.98e-04   3.68e-04   1.87e-02   9.87e+01s
9000   1.2042e+04   3.10e-04   2.46e-04   1.87e-02   1.01e+02s
9200   1.2043e+04   2.24e-04   4.20e-04   1.87e-02   1.03e+02s
9400   1.2042e+04   2.69e-04   3.78e-04   1.87e-02   1.05e+02s
9600   1.2042e+04   2.80e-04   2.56e-04   1.87e-02   1.08e+02s
9800   1.2043e+04   2.09e-04   4.05e-04   1.87e-02   1.10e+02s
10000   1.2042e+04   2.78e-04   3.05e-04   1.87e-02   1.12e+02s

status:               solved inaccurate
number of iterations: 10000
optimal objective:    12042.4795
run time:             1.12e+02s
optimal rho estimate: 9.00e-03

We still see some worrisome warning messages:

Warning messages:
1: In size(const) : NAs introduced by coercion to integer range
2: In size(const) : NAs introduced by coercion to integer range
3: In size(const) : NAs introduced by coercion to integer range

#### Formulation 3: exploit sparsity

This is a more difficult exercise.  CVXPY and CVXR do not have a conditional generation facility. Or in other words, we cannot generate a sparse variable. In practical modeling, this is rather essential. So how can we generate only variables that are related to existing links? Well, we can do everything ourselves and basically do all the work that CVXPY/CVXR is supposed to do. We can also fudge things a bit. We still generate all variables, but we multiply each variable in each constraint by 0 if it should not exist. This is only a partial solution, as we still generate all variables. However, a variable with all zeros in its column should be removed by the solver. We hope that is happening, but we have little feedback that is actually happening.

In CVXPY I did:

#----------------------------------------------------
# model 3
# exploit sparsity
#----------------------------------------------------

# X has same shape as our cost matrix
X = cp.Variable(cost.shape,"X")

prob = cp.Problem(
)

prob.solve(verbose=True)
print("status:",prob.status)
print("objective:",prob.value)

This actually seems to help a bit:

C:\tmp>py trnsport_cvx.py
CHECK: sum capacity = 10000.000000000016, sum demand = 10000.000000000002
===============================================================================
CVXPY
v1.2.1
===============================================================================
(CVXPY) Aug 10 01:57:25 AM: Your problem has 250000 variables, 3 constraints, and 0 parameters.
(CVXPY) Aug 10 01:57:25 AM: It is compliant with the following grammars: DCP, DQCP
(CVXPY) Aug 10 01:57:25 AM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)
(CVXPY) Aug 10 01:57:25 AM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
(CVXPY) Aug 10 01:57:25 AM: Compiling problem (target solver=ECOS).
(CVXPY) Aug 10 01:57:25 AM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> ECOS
(CVXPY) Aug 10 01:57:25 AM: Applying reduction Dcp2Cone
(CVXPY) Aug 10 01:57:25 AM: Applying reduction CvxAttr2Constr
(CVXPY) Aug 10 01:57:25 AM: Applying reduction ConeMatrixStuffing
(CVXPY) Aug 10 01:57:25 AM: Applying reduction ECOS
(CVXPY) Aug 10 01:57:25 AM: Finished problem compilation (took 1.830e-01 seconds).
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
(CVXPY) Aug 10 01:57:25 AM: Invoking solver ECOS  to obtain a solution.

ECOS 2.0.10 - (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  +5.471e+04  +5.471e+04  +2e+06  5e-01  3e-01  1e+00  1e+01    ---    ---    1  1  - |  -  -
1  +3.188e+04  +3.188e+04  +1e+06  3e-01  1e-01  2e+00  5e+00  0.6991  3e-01   1  1  1 |  0  0
2  +2.670e+04  +2.670e+04  +1e+06  2e-01  1e-01  3e+00  4e+00  0.3511  6e-01   1  1  0 |  0  0
3  +1.558e+04  +1.558e+04  +4e+05  8e-02  4e-02  2e+00  2e+00  0.9519  3e-01   1  1  0 |  0  0
4  +1.381e+04  +1.381e+04  +2e+05  3e-02  2e-02  8e-01  7e-01  0.7207  2e-01   1  0  0 |  0  0
5  +1.318e+04  +1.318e+04  +1e+05  2e-02  1e-02  4e-01  4e-01  0.6765  5e-01   1  0  0 |  0  0
6  +1.261e+04  +1.261e+04  +5e+04  8e-03  5e-03  2e-01  2e-01  0.6174  1e-01   1  0  0 |  0  0
7  +1.235e+04  +1.235e+04  +3e+04  4e-03  3e-03  9e-02  1e-01  0.6661  3e-01   1  0  0 |  0  0
8  +1.218e+04  +1.218e+04  +1e+04  2e-03  1e-03  3e-02  5e-02  0.7733  3e-01   1  0  1 |  0  0
9  +1.211e+04  +1.211e+04  +6e+03  8e-04  6e-04  2e-02  2e-02  0.7198  3e-01   1  1  1 |  0  0
10  +1.206e+04  +1.206e+04  +3e+03  4e-04  3e-04  8e-03  1e-02  0.7449  3e-01   1  1  1 |  0  0
11  +1.204e+04  +1.204e+04  +1e+03  2e-04  1e-04  3e-03  5e-03  0.6063  9e-02   1  1  1 |  0  0
12  +1.202e+04  +1.202e+04  +5e+02  7e-05  5e-05  1e-03  2e-03  0.6782  1e-01   1  1  1 |  0  0
13  +1.202e+04  +1.202e+04  +2e+02  2e-05  2e-05  4e-04  6e-04  0.7878  1e-01   1  1  1 |  0  0
14  +1.202e+04  +1.202e+04  +9e+01  1e-05  9e-06  2e-04  4e-04  0.4929  2e-01   1  1  1 |  0  0
15  +1.202e+04  +1.202e+04  +3e+01  3e-06  3e-06  6e-05  1e-04  0.8568  2e-01   1  1  1 |  0  0
16  +1.202e+04  +1.202e+04  +2e+01  2e-06  2e-06  4e-05  6e-05  0.6638  4e-01   1  1  1 |  0  0
17  +1.202e+04  +1.202e+04  +2e+00  2e-07  2e-07  5e-06  7e-06  0.8945  2e-02   1  1  1 |  0  0
18  +1.202e+04  +1.202e+04  +2e-01  3e-08  2e-08  5e-07  8e-07  0.9341  4e-02   2  1  1 |  0  0
19  +1.202e+04  +1.202e+04  +2e-03  3e-10  2e-10  6e-09  9e-09  0.9890  1e-03   3  0  0 |  0  0
20  +1.202e+04  +1.202e+04  +3e-05  3e-12  3e-12  6e-11  1e-10  0.9890  1e-04   5  0  0 |  0  0

OPTIMAL (within feastol=3.4e-12, reltol=2.2e-09, abstol=2.6e-05).
Runtime: 2.398611 seconds.

-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
(CVXPY) Aug 10 01:57:27 AM: Problem status: optimal
(CVXPY) Aug 10 01:57:27 AM: Optimal value: 1.202e+04
(CVXPY) Aug 10 01:57:27 AM: Compilation took 1.830e-01 seconds
(CVXPY) Aug 10 01:57:27 AM: Solver (including time spent in interface) took 2.408e+00 seconds
status: optimal
objective: 12015.218125523525

C:\tmp>

This is about twice as fast as method 2.

I also tried this with CVXR/OSQP, but the result was still a non-optimal solution: solved inaccurate. After increasing the iteration limit the status became solved. However, the objective value was not really close to the optimal objective.

#### Conclusion

The conclusion is not encouraging for OSQP. That solver has difficulties finding the optimal solution.

CVXPY/ECOSCVXR/OSQP
Large CostFixed VariablesSparse FormulationLarge CostFixed VariablesSparse Formulation
Variables250,000250,000250,000250,000250,000250,000
Constraintsnananananana
StatusOptimalOptimalOptimalsolved solvedsolved
Objective12015.21812015.21812015.21894120.821012042.462912042.4477
Compilation Time5.60.40.1
Solver Time34.54.42.369.316993.8
Iterations2592020107751445016100

CVXPY and CVXR do not allow variables with more than 2 dimensions. Also, it does not support sparse variables. ECOS delivers correct optimal solutions, but OSQP gives suboptimal solutions when solving this very large but very easy LP. Note that this is with default tolerances eps_rel and eps_abs. We can tighten these to get better solutions, but at a considerable computational cost. OSQP is a bit of a special solver based on a so-called first-order method (allowing relatively large problems to be approximated quickly using limited amounts of memory). It may not be the most suited as a default solver.

It is interesting to see that our haphazard sparse formulation works in the sense that it gives the correct solution, but in less time (compilation time and solver time).

The sparse formulation in CVXPY/ECOS is about 10 times as slow as GAMS/CPLEX [1].

#### References

1. Transportation model with some non-existing links, https://yetanothermathprogrammingconsultant.blogspot.com/2022/07/transportation-model-with-some-non.html
2. Welcome to CVXPY 1.2, https://www.cvxpy.org/
3. Convex Optimization in R, https://cvxr.rbind.io/
4. Alexander Domahidi, Eric Chu, Stephen Boyd, ECOS: An SOCP solver for embedded systems,  Proceedings European Control Conference, pages 3071-3076, Zurich, July 2013.
5. OSQP, https://osqp.org/

#### Appendix: Python code

import numpy as np
import pickle
import pandas as pd
import cvxpy as cp

#
# get the GAMS data
# stored in a pickle file
#

datafile = r"c:\tmp\assignment.pkl"
capacity = data['capacity']
demand = data['demand']
cost = data['cost']

#
# we should have a balanced problem
#
print(f"CHECK: sum capacity = {sum(capacity)}, sum demand = {sum(demand)}")

#-------------------------------------------------
# problem 1: large cost
#-------------------------------------------------

#-------------------------------------------------
# dense transportation model
# matrix formulation
#-------------------------------------------------

# X has same shape as our cost matrix
X = cp.Variable(cost.shape,"X")

# two different e vectors
ei = np.ones(capacity.shape)
ej = np.ones(demand.shape)

prob = cp.Problem(
cp.Minimize(cp.trace(largecost.T@X)),
[X.T@ei >= demand,
X@ej <= capacity,
X>=0])
prob.solve(verbose=True,max_iters=1000)
print("status:",prob.status)
print("objective:",prob.value)

#-------------------------------------------------
# dense transportation model
# sum formulation
#-------------------------------------------------

# X has same shape as our cost matrix
X = cp.Variable(cost.shape,"X",nonneg=True)

prob = cp.Problem(
cp.Minimize(cp.sum(cp.multiply(largecost,X))),
[cp.sum(X,axis=0) >= demand,
cp.sum(X,axis=1) <= capacity])
prob.solve(verbose=True,max_iters=1000)
print("status:",prob.status)
print("objective:",prob.value)

#----------------------------------------------------
# model 2
# fix variables
#----------------------------------------------------

# X has same shape as our cost matrix
X = cp.Variable(cost.shape,"X",nonneg=True)

prob = cp.Problem(
cp.Minimize(cp.sum(cp.multiply(cost,X))),
[cp.sum(X,axis=0) >= demand,
cp.sum(X,axis=1) <= capacity,

prob.solve(verbose=True)
print("status:",prob.status)
print("objective:",prob.value)

#----------------------------------------------------
# model 3
# exploit sparsity (somewhat)
#----------------------------------------------------

# X has same shape as our cost matrix
X = cp.Variable(cost.shape,"X")

prob = cp.Problem(
)

prob.solve(verbose=True)
print("status:",prob.status)
print("objective:",prob.value)

#### Appendix: R code

#-----------------------------------------------------------
# libraries
#-----------------------------------------------------------

library(CVXR)

#-----------------------------------------------------------
#-----------------------------------------------------------

#-----------------------------------------------------------
# convert from dataframes to matrices/vectors
#-----------------------------------------------------------

convSup <- 1:nrow(capacity)
names(convSup) <- capacity$i convDem <- 1:nrow(demand) names(convDem) <- demand$j

cost<-matrix(data=0,nrow=length(convSup),ncol=length(convDem))
# loop is fast enough
for (r in 1:nrow(c)) {
cost[convSup[c$i[r]],convDem[c$j[r]]] = c$value[r] } islink<-matrix(data=0,nrow=length(convSup),ncol=length(convDem)) # loop is fast enough for (r in 1:nrow(plink)) { islink[convSup[plink$i[r]],convDem[plink$j[r]]] = plink$value[r]
}

cap <- capacity$value dem <- demand$value

paste(sum(cap),sum(dem))

n <- length(cap)
m <- length(dem)

#-----------------------------------------------------------
# Model 1: large cost
#-----------------------------------------------------------

X <- Variable(rows=n,cols=m,name='X',nonneg=T)

prob <- Problem(Minimize(sum_entries(largecost*X)),
constraints = list(
sum_entries(X,axis=1) >= dem,
sum_entries(X,axis=2) <= cap))

res <- solve(prob,verbose=T,max_iter=1000000)
res$status res$value

#-----------------------------------------------------------
# Model 2: Fixed Variables
#-----------------------------------------------------------

X <- Variable(rows=n,cols=m,name='X',nonneg=T)

prob <- Problem(Minimize(sum_entries(cost*X)),
constraints = list(
sum_entries(X,axis=1) >= dem,
sum_entries(X,axis=2) <= cap,

res <- solve(prob,verbose=T,max_iter=1000000)
res$status res$value

#-----------------------------------------------------------
# Model 3: Exploit sparsity
#-----------------------------------------------------------

X <- Variable(rows=n,cols=m,name='X')

res$status res$value