Background theory

Background theory

In one sentence: the operations-research ideas behind the Julia side — which location problem this is, why the LP relaxation certifies a lower bound, what a λ-sweep can and cannot reach, why the solver is a warm-started simplex, and the three algorithms inside the rounding — with the proofs the neighbouring pages only cite, checked against lp_run.jl.

Where this sits: Julia builds one LP from the GeoDMS OD matrix (The facility location model explained) → the λ-sweep solves it along a price grid (The lambda sweep and the Pareto frontier) → rounding turns fractions into pharmacies (From LP fractions to real pharmacies) → scenarios (Scenarios S1 S2 S3). This page owns the theory; the neighbours own pipeline detail and numbers; Maarten's formal statement is Lambda sweep method. Written against commit e25416b, September 2026.

1. Which family of location problem this is

Facility location chooses which candidate sites to open and which client each serves. Three classic families:

• Uncapacitated facility location (UFL) (Glossary) — minimise travel plus a fixed cost per open site; the number of sites is an output.
• p-median (Glossary; Hakimi 1964, integer-LP form by ReVelle & Swain 1970) — minimise demand-weighted travel with exactly p sites; the count is an input.
• Covering models — fewest sites so every client is within a standard D; travel counts only through the yes/no "within D".

build_lp_warmstart in lp_run.jl (around line 397) is a UFL in which every site has the same fixed cost λ, so one price traces the whole count-versus-travel trade-off. Two features pull it towards the other families. Soft coverage (Glossary): the client constraint is Σ_j x_ij ≤ 1, not = 1, and an unserved client pays the stranding price BIG — the covering idea in disguise, introduced (commit 4ec47a9, July 2026) so that the frontier runs past today's count and S1 and S2 can be bracketed everywhere (The facility location model explained; BIG per cost function on Travel cost functions). And S1 — least travel at today's pharmacy count (Scenarios S1 S2 S3) — is a p-median instance: the sweep does not solve it, it hopes to pass through it; §5 says when it cannot.

2. Integer programme, LP relaxation, and why a relaxation bounds

The honest problem has y_j ∈ {0, 1}: a pharmacy is open or not — an integer programme (IP; Glossary). The LP relaxation (LP for linear programme: linear objective and constraints, variables free to take any real value within their bounds; Glossary) replaces {0, 1} by 0 ≤ y_j ≤ 1 (around line 406), allowing fractions of a pharmacy. An LP with a million OD rows solves in seconds to minutes; the IP would need branch-and-bound (Glossary), the exact method that repeatedly splits on a fractional y and re-solves LPs — hours per point at this size. Why no exact solve is run per λ is not recorded; the likely reason is cost across ~88 sweeps while the LP already supplies the bound (inference). The roadmap keeps one fixed-count exact solve to pin S1/S2 (doc/todo.md item 5).

Why a relaxation is a bound: every integer solution is still allowed in the relaxed problem, so the relaxation minimises over a larger set and its optimum can only be lower or equal. And any concrete set of pharmacies you exhibit — the multistart set — is allowed in the IP, so its cost is higher or equal:

LP relaxation  ≤  integer optimum  ≤  any feasible integer solution

This chain is the basis of the "certified" language on the results pages (§8).

3. The integrality gap, and why the strong formulation is nearly integral

The integrality gap (Glossary) is the distance between the LP optimum and the integer optimum. Often it is zero. The three-clients-on-a-ring toy on The facility location model explained (one resident per client, every client 1 minute from two of the three candidates and 3 minutes from the third, λ = 2) has a real one: the best integer answer costs 7, the LP's y = (½, ½, ½) costs 6 — gap 1/7 ≈ 14 %. What that page does not show is why the LP cannot do better than 6. The reason is LP duality: every LP has a twin "pricing" problem whose best value equals the LP's optimum, so any valid set of prices is a floor under every LP solution. Here the prices are easy to exhibit: charge every client 2 and let it pass 1 of that on to each of its two 1-minute sites. Each site then collects exactly its λ = 2 — none is over-subscribed — and what a client keeps towards any site (2 minus what it passes on) never exceeds its travel time there (1 ≤ 1 at a near site, 2 ≤ 3 at the far one). Those prices are valid and worth 3 × 2 = 6, so the LP cannot go below 6 (HiGHS confirms it) and the gap of 1 is genuine: the geography resists a clean integer covering.

Why "strong". The linking constraint is written once per OD pair, x_ij ≤ y_j (around line 426) — the strong formulation (Glossary) — not once per site as Σ_i x_ij ≤ |I|·y_j, the weak one that lets the LP open a hundredth of a pharmacy anywhere and gives a useless bound (size cost and measured frac_y: The facility location model explained). In practice the strong LP comes out integral almost everywhere; the ring toy, where each client's cheap options chain into a loop of odd length, is the canonical exception. ReVelle & Swain reported the regularity in 1970 and Lambda sweep method §1 calls it "known to be nearly integral; the sweeps confirm it" — empirical, not a theorem.

4. Supporting lines: the two facts the code relies on

Put every achievable open set in the plane N = number of open pharmacies, T = total population travel cost. The objective T + λ·N is constant along lines of slope −λ, so minimising it slides such a line down until it last touches the cloud of achievable points: a supporting line (Glossary) touching the Pareto frontier (Glossary), the lower-left outline of the cloud. The picture, a worked three-client sweep and the charts are on The lambda sweep and the Pareto frontier; here are the two consequences the code relies on, with their proofs.

• The count never rises when λ rises. If S₁ is optimal at λ₁ and S₂ at λ₂ > λ₁, then S₁ is at least as good as S₂ at price λ₁ and the reverse at λ₂; adding the two inequalities gives (λ₂ − λ₁)(N₂ − N₁) ≤ 0 — a higher price never buys more sites. The same holds for the LP with N = Σy* (y* denotes the LP's optimal y values), which is why the sweep can bisect S1 on sum_y; multistart travel is monotone only "up to rounding noise" (lambda_sweep_simplex.jl, comment around line 318), so the S2 bisection on cost_c, the multistart travel column, can wobble.
• The LP curve has no bumps. Any 50/50 mix of two allowed LP solutions is again an allowed LP solution, with the average count and the average travel. So the cloud of (count, travel) points the LP can produce has no dents, the sliding line traces a curve with no bumps, and λ only ever increases along it. A bump in a plotted LP lower-bound curve, or λ not monotone along it, means the points did not come from one and the same LP. Issue #54 applied this test (Maarten confirmed the LP bound must be convex) to the FRI LINEAR curve, where w (the grid value; λ = w × €100 000 — Glossary) ran 0.109 → 0.115 → 0.119 → 0.123 along the bound because a July sweep on a candidate subsample had been joined to September points on the full set (The lambda sweep and the Pareto frontier).

5. Lagrangian relaxation: why sweeping λ reaches only the lower convex envelope

Write V(p) for the least travel with exactly p sites open. The sweep never imposes Σy = p; it prices it. Moving a constraint into the objective with a multiplier is a Lagrangian relaxation (Glossary), and Cornuéjols, Fisher & Nemhauser (1977) did exactly this for the location problem:

L(λ) = min_S [ T(S) + λ·|S| ] − λ·p   ≤   V(p)   for every λ

In words: L(λ) is what the sweep's LP computes at price λ — the best travel-plus-facility total over every open set S, |S| being the number of sites in S — minus the λ·p you would have paid at the target count. It can never exceed V(p), because the best p-site set is one of the sets the minimum ranges over, and for that set the two λ terms cancel. Keeping the best of these bounds over all λ gives the highest curve with no bumps that fits under all the points (p, V(p)) — a line of slope −λ pushed up from below until it touches. That curve is the lower convex envelope (Glossary), and it only touches points with no neighbour-chord above them. Consequently:

a count p whose point (p, V(p)) lies above the chord joining its neighbours — a concave dent — is never produced, whatever λ.

The smallest example — three clients, four candidates, V(1) = 30, V(2) = 20, V(3) = 0, so the chord from (1, 30) to (3, 0) passes through (2, 15) < 20 — is worked through on The lambda sweep and the Pareto frontier: the count jumps 3 → 1 at λ = 15 and a bisection aimed at "today's count = 2" bounces between the endpoints. The LP relaxation of that toy behaves the same (solved with HiGHS while writing this page; the script is not in the repository): Σy* = 3 for λ < 15, 1 for λ > 15 (a tie at 15), never 2 — the dent is a property of the geography, not of relaxing integrality.

Chord interpolation, the in-sweep bisection and the roadmap's fixed-count exact solve are the pipeline's answers (Scenarios S1 S2 S3). Dents are shallow in practice: after issue #52, S1 was within 1 % or one facility of today's count in 87 of 88 sweeps (median 0.13 %, measured on the unpushed commit 216ac77; the bisection itself stops at REFINE_TOL = 0.2 %), and the exception, country-level Poland LINEAR, is a solve-time limit, not a dent.

6. Solving the LPs: warm-started dual simplex versus interior point

Simplex (Glossary → dual simplex) walks from corner to corner of the LP's feasible region. A corner is a vertex (Glossary), described by a basis (Glossary) — the set of variables allowed to be non-zero — and each step, a pivot (Glossary), swaps one variable in and one out.

solve_at_w! (around line 460) keeps one JuMP/HiGHS model for the whole sweep and, between two λ values, changes only the objective coefficients of the y_j. So the old optimal vertex is still a feasible vertex of the same region — only its reduced costs (Glossary), the "is it worth leaving this vertex along variable v?" numbers, changed, and for a small Δλ few of them change sign, so the old vertex is close to the new optimum. HiGHS detects the reusable basis itself: that is the warm start (Glossary), why the driver is sequential — "warm-start gives the speedup, not threads" (commit 011cafb) — and why the next grid point costs a small fraction of the cold solve (pivot counts from the Netherlands log: The lambda sweep and the Pareto frontier).

Why the log says dual simplex (Glossary): simplex has two variants. The primal one keeps every constraint satisfied and fixes wrong-signed reduced costs step by step; the dual one keeps the reduced costs right and fixes violated constraints. After a pure price change the old vertex still satisfies every constraint but some reduced costs have the wrong sign — the textbook case for the primal variant. HiGHS runs the dual simplex because that is its default (the code sets no simplex_strategy), so it first restores the reduced costs in a short phase 1 — the Ph1/Du counts on the iteration-0 row of every warm solve — and then iterates; the primal/dual choice is incidental to the saving. Warm start fails when the step is large: at the top of the grid a 1-2-5 step moves every y-coefficient by Δλ ≈ 10⁵–10⁶, thousands of reduced costs flip, and the LP is degenerate (Glossary) — many basic variables sit at zero, so many pivots make no progress.

Interior-point methods (IPM; Glossary) do not walk corners; they move through the inside of the region, each iteration solving one large linear system. The iteration count is roughly constant whatever the starting basis — "basis-distance-immune" (comment, around line 479) — so a single far solve is safe, but every solve costs the same, and when the optimum is not unique the answer is a blend of optimal vertices: many small fractional y_j. Crossover (Glossary) pushes that solution to a vertex and leaves a basis a later simplex can warm-start from; the code turns it on whenever SOLVER=ipm. The solver history (IPM in July 2026 and back, the IPM OTHER_ERROR on LOGISTIC, the time-out cascade) is on The lambda sweep and the Pareto frontier and Decision log; why the LINEAR soft LP survives IPM is not recorded.

7. Turning fractions into a set: the matheuristic

A matheuristic (Glossary) is a heuristic that uses a mathematical-programming solution as raw material. Here the LP's y* fixes the count p = round(Σy*) (around line 531; Julia rounds halves to even, so 1.5 → 2), splits candidates into must-open (y* ≥ 1 − 10⁻⁶), fractional and closed, and seeds a local search — so that "the discrete pattern stays an apples-to-apples point on the same Pareto curve" (issue #47). The step-by-step walkthrough, worked toy and MS_* knobs are on From LP fractions to real pharmacies; here are the three ideas inside it and why each is licensed.

Submodularity and the lazy greedy (CELF). A set function is submodular (Glossary) when adding an element to a smaller set helps at least as much as adding it to a larger one — diminishing returns. The travel saving of an open set is monotone submodular: once a client is served closely, no further site can save as much for it. Cornuéjols, Fisher & Nemhauser (1977) showed that plain greedy — add the site with the largest marginal gain, p times — reaches at least 1 − 1/e ≈ 63 % of the best possible saving for this location problem (Nemhauser, Wolsey & Fisher 1978 generalised it to all monotone submodular functions). That is the licence for a greedy seed. The CELF trick (Glossary; Leskovec et al. 2007): because gains only ever decrease as the set grows, a gain computed earlier is an upper bound on the current gain — so keep the stale gains, recompute only the largest, and if it is still the largest it is the true best. greedy_round (lines 23–88) does this and after each opening invalidates only candidates sharing a client whose best cost just improved — exact, because a gain can only change through such a client. If greedy stalls before p sites it tops up by y* descending so the count matches top-p, the seed that simply opens the p candidates with the largest y*.

Weighted reservoir sampling (A-Res). Efraimidis & Spirakis (2006) draw k items without replacement with probability proportional to weights w_j in one pass (A-Res, Glossary): give each item the key u_j^{1/w_j} with u_j uniform on (0, 1) and keep the k largest keys. In words: u is below 1, so raising it to the power 1/w with a large w pushes it towards 1 — heavy items get keys near the top, light items keys near zero. randomized_seed (lines 148–156) uses log(rand)/y_j; log is increasing, so the ordering is the same. Issue #47 has Maarten's proof that the draw lands on j with probability y_j/Σy — "the relaxed solution is read directly as selection probabilities"; the generator has a fixed seed, so the randomness only diversifies the starting points.

Vertex substitution and the fast interchange. The classic p-median local search (vertex substitution, Glossary; Teitz & Bart 1968, not cited in the repository) swaps one open site out and one closed site in, accepting the swap if travel falls. Resende & Werneck (2007) showed the exact change from swapping f in and r out decomposes into three cheap parts — what f saves on its own, minus what closing r costs, plus a correction for the clients who lose r but are caught by f:

profit(f, r) = gain(f) − loss(r) + extra(f, r)

swap_round! (lines 161–209) computes all three in one pass over the OD, from each client's nearest and second-nearest open cost, and applies the best profit above 10⁻⁶, at most MS_ROUNDS = 12 times per seed; since every accepted swap strictly lowers travel, the polish terminates. The worked swap, and why only fractional-y* sites may enter or leave, are on From LP fractions to real pharmacies.

8. Certified optimality gaps

Per λ the sweep has travel_relax — LP travel including the unserved fraction priced at BIG (around lines 506–513; that is what makes the bound coverage-honest) — at count Σy*, and travel_c_multi at count n_open = round(Σy*). The rigorous statement is on the λ-priced total (Glossary → certified gap):

travel_relax + λ·Σy*   ≤   integer optimum at λ   ≤   travel_multi + λ·n_open

The deck's "multi vs relax" column and the README's λ-axis chart compare the two travel numbers only; since |n_open − Σy*| ≤ ½ the facility terms differ by up to λ/2, so that comparison can be off by that much and even negative wherever the count was rounded up. The worked toy, the measured gaps and the Netherlands LINEAR w = 0.02 row are on From LP fractions to real pharmacies. That is the precise content of the results page's "the certified bound is on total cost, not on the count". Under LOGISTIC the bound is on Σ pop·c(t), not on minutes; mean_t, the mean travel time of served residents, carries no bound. The gap mixes the integrality gap with whatever the heuristic leaves on the table; only an exact solve could separate them.

In the code

Idea	File → function	Lines (around)
UFL with soft coverage, strong linking	lp_run.jl → build_lp_warmstart	397–430
Objective-only update, warm start, solver choice, travel_relax	lp_run.jl → solve_at_w!	460–513
CELF greedy · A-Res draw · fast interchange · multistart	lp_run.jl → greedy_round, randomized_seed, swap_round!, multistart_round	23–246
Monotone bisection on sum_y / cost_c	lambda_sweep_simplex.jl → bisect!	343–414
Stranding price BIG, c(t)	settings.jl → c, big_cost	73–97

Permalink pinned to e25416b: solve_at_w!. Function by function: Code walkthrough lp_run and Code walkthrough lambda_sweep_simplex.

Pitfalls and open questions

• Convexity is a test, not a guarantee (§4): a bump in a frontier chart, or λ not monotone along it, means two different LPs were joined (issue #54) — Scope rules and data gaps.
• A concave dent is unreachable by any λ (§5). How often round(Σy*) lands in one in real areas is not measured; the 87/88 pins of issue #52 suggest rarely.
• The certified gap is on total cost per λ (§8); the "multi vs relax" travel column can be negative.
• Warm re-solve behaviour is confirmed only from the pre-soft-coverage Netherlands log; the IPM failure on LOGISTIC is stated in the solve_at_w! comment, not shown in any log (unverified beyond that comment).
• Inf gains in the greedy. best_cost starts at Inf for clients no must-open site reaches (line 31), so the first candidate reaching such a client wins with gain Inf (strict >, line 54): the greedy covers somebody first and weighs travel only after. Documented nowhere; harmless for the reported result (the polish re-scores every set with BIG) and the likely reason — inference, not recorded — why travel_greedy can be far worse than travel_multi at high λ.

References

One sentence each on the role in this pipeline; the formal list is in Lambda sweep method.

• ReVelle & Swain (1970), Central facilities location, Geographical Analysis 2 — the p-median as an integer LP with assignment (x) and openness (y) variables (notation adopted in commit 4cee238); its LP relaxation usually integral.
• Cornuéjols, Fisher & Nemhauser (1977), Location of bank accounts to optimize float, Management Science 23 — the Lagrangian relaxation (§5) and the greedy (1 − 1/e) guarantee (§7).
• Leskovec et al. (2007), Cost-effective outbreak detection in networks, KDD — CELF, the lazy greedy of greedy_round.
• Efraimidis & Spirakis (2006), Weighted random sampling with a reservoir, Information Processing Letters 97 — algorithm A-Res, the key used in randomized_seed.
• Resende & Werneck (2007), A fast swap-based local search procedure for location problems, Annals of Operations Research 150 — the gain/loss/extra decomposition of swap_round!.
• Albuquerque, Figueiredo & Genre-Grandpierre (2026), SSRN 7133060 — RSSV: random candidate subsets, spatial voting, then an exact solve on the survivors; the planned replacement for the stride subsample, not implemented as of September 2026 (known only from the deck and README — unverified).

See also

• The facility location model explained · The lambda sweep and the Pareto frontier · From LP fractions to real pharmacies · Scenarios S1 S2 S3
• Lambda sweep method · Code walkthrough lp_run · Code walkthrough lambda_sweep_simplex · Decision log · Alternatives not pursued · Glossary