Large Linear Programs

A linear program with many variables and few constraints is a common shape in applied optimization: allocate a budget across thousands of options, subject to a handful of requirements. General LP solvers see only a sparse matrix, which is the right view when nothing more is known. Here more is known: the constraint matrix has only p dense rows, and an interior-point method can exploit that to make every iteration cost O(p² M).

This guide works one example from start to finish. We state a media plan as a linear program, solve it with scipy’s linprog (HiGHS) to establish ground truth and a time to beat, explain the structure Cvxium exploits, and then solve it with Cvxium: about 28 times faster than linprog’s default method, and about 6 times faster than the fastest method HiGHS offers. Three settings decide whether that speedup appears at all: the choice of solver, the scale of the objective, and the starting point. We cover each, read the shadow prices off the solution, and close with how the speedup grows with problem size. Every python block runs, in order, as one script.

The media plan

Example 1. An advertiser buys impressions from M = 20,000 placements (sites, apps, time slots). Placement i sells up to u_i thousand impressions at c_i dollars per thousand. The advertiser cares about p = 50 audience segments, and a fraction R_ki of placement i’s impressions reaches segment k. Buying x_i thousand impressions from each placement delivers (R x)_k thousand impressions to segment k, and the campaign must deliver at least r_k to each. The cheapest plan solves

\[\begin{split} \begin{array}{ll} \mbox{minimize} & c^T x \\ \mbox{subject to} & R x \succeq r \\ & 0 \preceq x \preceq u. \end{array} \end{split}\]

There are M variables, p dense inequality rows, and 2M bounds. Each placement’s audience concentrates on a few segments, so we draw each column of R from a Dirichlet distribution with parameter 0.5, and we ask for 30% of the impressions available to each segment.

import time

import numpy as np
import scipy.optimize

M, p = 20_000, 50
rng = np.random.default_rng(1)
R = rng.dirichlet(np.full(p, 0.5), size=M).T  # shape (p, M); columns sum to 1
u = rng.uniform(50.0, 150.0, M)  # capacity, thousands of impressions
c = rng.lognormal(np.log(5.0), 0.5, M)  # dollars per thousand impressions
r = 0.3 * (R @ u)  # thousands of impressions required per segment

A baseline with scipy

scipy.optimize.linprog takes inequalities as A_ub x ≤ b_ub, so we negate the reach constraints. Its default method, "highs", lets HiGHS choose an algorithm, and on this problem it chooses the dual simplex. HiGHS’s interior-point method, "highs-ipm", is faster here, so to steel-man the alternative we time both.

bounds = np.column_stack([np.zeros(M), u])
c_millions = c / 1e6  # see "Scale the objective" below

start = time.perf_counter()
reference = scipy.optimize.linprog(c_millions, A_ub=-R, b_ub=-r, bounds=bounds)
scipy_seconds = time.perf_counter() - start
assert reference.status == 0

start = time.perf_counter()
reference_ipm = scipy.optimize.linprog(
    c_millions, A_ub=-R, b_ub=-r, bounds=bounds, method="highs-ipm"
)
scipy_ipm_seconds = time.perf_counter() - start
assert reference_ipm.status == 0
print(
    f"HiGHS: ${reference.fun:.4f}M in {scipy_seconds:.2f} s (default), "
    f"{scipy_ipm_seconds:.2f} s (interior point)"
)

The cheapest plan costs about $1.74 million. The default method finds it in about 4.3 seconds, the interior-point method in about 0.9. (All timings in this guide come from an Apple M2 Pro with scipy 1.17 and NumPy 2.4; expect different absolute numbers elsewhere, but similar ratios.)

The structure Cvxium exploits

An interior-point method replaces the inequalities with a logarithmic barrier and takes Newton steps on

\[ t\, c^T x - \sum_{i=1}^M \log x_i - \sum_{i=1}^M \log(u_i - x_i) - \sum_{k=1}^p \log\big((R x)_k - r_k\big). \]

The objective is linear, so it contributes nothing to the Hessian. The bounds touch one variable each, so they contribute a diagonal matrix D. The reach constraints contribute R^T S^{-2} R, where S = diag(R x - r), a matrix of rank p. The Newton system is therefore

\[ \big(D + R^T S^{-2} R\big)\, \Delta x = -g, \]

a diagonal plus a rank-p update. A dense solve would cost O(M³), about 8 × 10¹² operations at M = 20,000. The Woodbury identity instead reduces it to a p × p system:

\[ \big(D + R^T S^{-2} R\big)^{-1} = D^{-1} - D^{-1} R^T \big(S^2 + R D^{-1} R^T\big)^{-1} R D^{-1}. \]

Forming R D^{-1} R^T costs O(p² M), factoring it O(p³), and everything else O(p M). The primal-dual method, which we use below, solves a system with the same pattern, only with D and S weighted by the dual variables. Interior-point methods need a few dozen iterations, nearly independent of M, so the whole solve costs O(p² M) times a small constant.

Cvxium assembles this structure without our help. LinearObjective reports a ZeroSystem as its Hessian, BoundsConstraintLayer a DiagonalSystem, and LinearInequalityConstraintLayer a low-rank update; summing them yields a LowRankUpdatedSystem over a diagonal base, whose solve applies the Woodbury identity.

Solving with Cvxium

Example 1, cont’d. The problem is a composition of three library components:

from cvxium import (
    BoundsConstraintLayer,
    LinearInequalityConstraintLayer,
    LinearObjective,
    PrimalDualSolver,
    Problem,
    SolveStatus,
)

problem = Problem.minimize(LinearObjective(c_millions)).subject_to(
    ineq=[
        BoundsConstraintLayer(lb=0.0, ub=u),
        LinearInequalityConstraintLayer(A=R, cl=r),
    ]
)

x0 = 0.5 * u  # buys half of every placement: strictly feasible
start = time.perf_counter()
result = problem.solve(solver=PrimalDualSolver(), x0=x0)
cvxium_seconds = time.perf_counter() - start

x = result.solution
assert result.status == SolveStatus.OPTIMAL
assert np.all(R @ x >= r) and np.all(x >= 0.0) and np.all(x <= u)
assert abs(result.objective_value - reference.fun) <= 1e-5
print(
    f"Cvxium: ${result.objective_value:.4f}M in {cvxium_seconds:.2f} s "
    f"({result.n_iterations} iterations), "
    f"{scipy_seconds / cvxium_seconds:.0f}x faster than linprog's default, "
    f"{scipy_ipm_seconds / cvxium_seconds:.1f}x faster than HiGHS's IPM"
)

Cvxium takes 19 iterations and about 0.15 seconds: roughly 28 times faster than linprog’s default and 6 times faster than HiGHS’s interior-point method, and it agrees with both on the optimal cost to better than one part in a million. The solution reports a duality-gap bound, result.suboptimality, of about 6 × 10⁻⁷: the plan costs at most sixty cents more than the cheapest possible plan.

The plans themselves agree less closely than their costs. HiGHS’s simplex returns a vertex: with 50 binding reach constraints, at most 50 placements are bought partially, and every other placement is bought fully or not at all. An interior-point method stops just inside the feasible region, so Cvxium’s plan sits a duality gap’s worth from that vertex, with each placement’s purchase near, but not exactly at, its bound. Counting a placement as fractional only when it sits more than 0.1% of its capacity from either bound finds 51, nearly the vertex’s 50; where the optimum is not unique, the two methods may also settle on different optimal plans of equal cost.

Shadow prices

The multipliers on the reach constraints answer the question the advertiser will ask next: what does one more thousand impressions in segment k cost? inequality_multipliers stacks the layers’ multipliers in order, so the last p belong to the reach layer. HiGHS reports the same quantities as ineqlin.marginals, negated because we negated the constraints.

shadow_price = 1e6 * result.inequality_multipliers[-p:]  # $ per thousand
np.testing.assert_allclose(
    shadow_price, -1e6 * reference.ineqlin.marginals, rtol=1e-3
)
cheapest, dearest = np.argmin(shadow_price), np.argmax(shadow_price)
print(
    f"An extra thousand impressions costs ${shadow_price[cheapest]:.2f} in "
    f"segment {cheapest} and ${shadow_price[dearest]:.2f} in segment {dearest}"
)

Every segment’s constraint binds, and the marginal costs range from about \(1.70 to about \)5.50 per thousand impressions. These numbers often matter more than the plan itself: they tell the advertiser which targets are expensive to raise and which are nearly free.

Getting the settings right

The speedup above depends on three choices. Each default is reasonable for some problems, and each is wrong for this one.

Use the primal-dual solver

BarrierSolver, the more general default elsewhere in Cvxium, solves a sequence of centering problems, one per value of the barrier parameter t, each to high accuracy. On a linear program that costs far more Newton steps than the primal-dual method, which updates t at every step:

from cvxium import BarrierSolver

start = time.perf_counter()
barrier_result = problem.solve(solver=BarrierSolver(), x0=x0)
barrier_seconds = time.perf_counter() - start
newton_steps = sum(barrier_result.barrier_diagnostics.inner_nits)
assert abs(barrier_result.objective_value - reference.fun) <= 1e-5
print(f"Barrier: {newton_steps} Newton steps in {barrier_seconds:.2f} s")

The barrier method takes 88 Newton steps over 11 centering problems, about 0.7 seconds: faster than either HiGHS method, but four to five times slower than the primal-dual method. For an LP, start with PrimalDualSolver.

Scale the objective

Cvxium’s stopping rule is absolute: OuterSettings.tolerance, 10⁻⁶ by default, bounds the duality gap in the objective’s own units. We expressed costs in millions of dollars, so the default asks for a plan within a dollar of optimal. Had we left them in dollars, the default would ask for a millionth of a dollar on a $1.74 million plan, a relative accuracy of 10⁻¹², beyond what double precision can deliver through the barrier. The primal-dual solver then grinds to its iteration cap: in our measurement it ran all 1,000 iterations, took about 12 seconds, and returned SolveStatus.ACCEPTABLE.

Either rescale the objective so its optimal value is of order one, or state the tolerance in the objective’s units:

from cvxium import OptimizationSettings, OuterSettings

dollars = Problem.minimize(LinearObjective(c)).subject_to(
    ineq=[
        BoundsConstraintLayer(lb=0.0, ub=u),
        LinearInequalityConstraintLayer(A=R, cl=r),
    ]
)
within_a_dollar = OptimizationSettings(outer=OuterSettings(tolerance=1.0))
dollar_result = dollars.solve(solver=PrimalDualSolver(within_a_dollar), x0=x0)
assert dollar_result.status == SolveStatus.OPTIMAL
assert abs(dollar_result.objective_value - 1e6 * reference.fun) <= 10.0

This solve takes the same 19 or so iterations as before.

Start from a strictly feasible point

An interior-point method needs a starting point strictly inside every constraint. Without one, Cvxium first solves a Phase I problem to find it, which here adds about a quarter to the solve time. Problems like this one usually have an obvious interior point: buying half of every placement delivers 50% of each segment’s available impressions, comfortably above the 30% target, and stays strictly inside the bounds. Passing it as x0 skips Phase I. (Cvxium checks x0 and falls back to Phase I if it is not strictly feasible.)

start = time.perf_counter()
phase1_result = problem.solve(solver=PrimalDualSolver())
phase1_seconds = time.perf_counter() - start
assert abs(phase1_result.objective_value - reference.fun) <= 1e-5
print(f"Without x0: {phase1_seconds:.2f} s; with x0: {cvxium_seconds:.2f} s")

How the speedup scales

The table compares linprog’s default method, HiGHS’s interior-point method, and Cvxium (PrimalDualSolver from x0 = u/2) on media plans generated as above, timing the fastest of three runs of each. The last two columns are Cvxium’s speedups.

M	p	default (s)	HiGHS IPM (s)	Cvxium (s)	iterations	vs. default	vs. IPM
2,000	50	0.094	0.076	0.019	18	5×	4.0×
5,000	50	0.43	0.19	0.036	19	12×	5.4×
20,000	50	4.2	0.92	0.15	19	28×	6.0×
50,000	50	19.2	2.44	0.43	20	45×	5.6×
20,000	10	0.17	0.19	0.084	17	2×	2.2×
20,000	100	28.1	2.14	0.31	22	89×	6.8×
20,000	200	265	4.43	0.88	28	300×	5.0×

Cvxium’s iteration count stays between 17 and 28 across every size, and its time grows linearly in M, as the analysis above predicts. The default simplex grows much faster, in both M and p, because the number of pivots grows with the problem as well as the cost of each one; at p = 200 it takes over four minutes. HiGHS’s interior-point method scales like Cvxium’s, so against it the advantage holds steady at a factor of 4 to 7, less when p is small and fixed costs dominate.

The default deserves the attention because it is what linprog runs unless told otherwise, and nothing about this problem warns that the simplex will be slow. Whichever solver we compare with, the lesson is the same: time the alternatives on the problem at hand.

The comparison has limits worth naming. HiGHS solves problems Cvxium does not, including sparse LPs with many constraints, where a diagonal-plus-low-rank Hessian is no longer cheap, and mixed-integer programs; it also returns a vertex solution, which some applications need. The advantage here comes from the shape of the problem: many variables, few dense constraints, and bounds.

Summary

A linear program with M variables, p dense constraints, and bounds has a barrier Hessian that is a diagonal plus a rank-p update, so each interior-point iteration costs O(p² M). Cvxium assembles that structure from library components, LinearObjective, BoundsConstraintLayer, and LinearInequalityConstraintLayer, with no code of our own, and on the media plan solves the problem about 28 times faster than linprog’s default and 6 times faster than HiGHS’s interior-point method, while reporting the same optimal cost and the same shadow prices.

To get that speed on a problem of your own:

1. Express the problem with library layers: bounds as a BoundsConstraintLayer, dense rows as a LinearInequalityConstraintLayer (or an EqualityConstraintLayer for equalities).

2. Solve it with PrimalDualSolver.

3. Scale the objective so its optimal value is of order one, or set OuterSettings(tolerance=...) in the objective’s units.

4. Pass a strictly feasible x0 when one is obvious.

5. Check the answer against linprog once, and read the shadow prices from inequality_multipliers.