Optimization

Motivation

Second-order cone programming (SOCP) minimizes a linear function over the intersection of an affine subspace with a product of Lorentz cones. It contains convex quadratically constrained programs, robust linear programs with ellipsoidal uncertainty, and many norm-minimization problems, and it is solved in practice by primal–dual interior-point methods. At each iteration such a method computes a search direction from a linearized, symmetrized form of the complementarity conditions; the symmetrization depends on a scaling matrix $G$, and different choices of $G$ give different directions. The two directions used in practice are the HKM direction (Helmberg–Rendl–Vanderbei–Wolkowicz, Kojima–Shindoh–Hara, Monteiro) and the NT direction (Nesterov–Todd).

A desirable property of a search direction is scale invariance: if the problem is transformed by an automorphism of the cone, the direction computed for the transformed problem is the image of the original direction. A scale-invariant method behaves identically on all problems related by such a change of coordinates, so its behaviour does not depend on an arbitrary choice of coordinates within each cone.

Timeline. Tsuchiya (1999, Optim. Methods Softw. 11 & 12, pp. 141–182; reference [27] of the paper) proved that the HKM and NT directions for SOCP are scale invariant, with a two-page argument for each. Tütüncü, Toh and Todd (2003), in the paper that describes the solver SDPT3 version 3.0, derived closed forms for the blocks $(\mathcal E_i)^{-1}\mathcal F_i$ of the Schur complement matrix ((23) for HKM, (29) for NT) and used them to give a short proof of Tsuchiya's result (their Proposition 1).

Setting

There are $n_q$ second-order cone blocks. Block $i$ has dimension $q_i \ge 1$, and a vector $x_i \in \mathbb R^{q_i}$ is split into its first coordinate $x_i^0$ and the rest $x_i^1$. The cone is $K_q^{q_i} = \{x : x^0 \ge \|x^1\|_2\}$, and with $\bar J = \operatorname{diag}(1, -I)$ and $J = -\bar J$, the function $\gamma(x) = \sqrt{(x^0)^2 - \langle x^1, x^1 \rangle} = \sqrt{x^\top \bar J x}$ is defined on the cone.

The data are matrices $A_i \in \mathbb R^{q_i \times m}$, a vector $b \in \mathbb R^m$ and vectors $c_i \in \mathbb R^{q_i}$. The primal problem minimizes $\sum_i \langle c_i, x_i \rangle$ subject to $\sum_i A_i^\top x_i = b$ and $x_i \in K_q^{q_i}$; the dual maximizes $b^\top y$ subject to $A_i y + z_i = c_i$ and $z_i \in K_q^{q_i}$.

The arrow matrix of $v$ is $\operatorname{Arw}(v) = \begin{pmatrix} v^0 & (v^1)^\top \\ v^1 & v^0 I \end{pmatrix}$. Given an interior iterate $(x, y, z)$ and a scaling block $G_i$, set

$$\mathcal E_i = \operatorname{Arw}(G_i^{-1} z_i)\, G_i, \qquad \mathcal F_i = \operatorname{Arw}(G_i x_i)\, G_i^{-1}, \qquad T_{G_i}(x_i, z_i) = \operatorname{Arw}(G_i x_i)(G_i^{-1} z_i).$$

With $\mu = \langle x, z \rangle / n_q$ and a centering parameter $\sigma$, a search direction $(\Delta x, \Delta y, \Delta z)$ solves the Newton system

$$A_i \Delta y + \Delta z_i = c_i - z_i - A_i y, \qquad \sum_i A_i^\top \Delta x_i = b - \sum_i A_i^\top x_i, \qquad \mathcal E_i \Delta x_i + \mathcal F_i \Delta z_i = \sigma\mu\, e_1 - T_{G_i}(x_i, z_i).$$

The HKM direction takes $G_i$ from the explicit formula (19) in $z_i$; the NT direction takes $G_i$ from (20)–(21) in $x_i$ and $z_i$.

The automorphism group $\mathcal G_i$ of $K_q^{q_i}$ consists of the invertible linear maps $F$ with $F K_q^{q_i} = K_q^{q_i}$; these are exactly the matrices with $F^\top \bar J F = \lambda^2 \bar J$ for some $\lambda > 0$ and $F_{00} > 0$. Scaling by $F = \operatorname{diag}(F_1, \dots, F_{n_q})$ sends the data and iterate to $\hat A_i = F_i^{-1} A_i$, $\hat b = b$, $\hat c_i = F_i^{-1} c_i$, $\hat x_i = F_i^\top x_i$, $\hat y = y$, $\hat z_i = F_i^{-1} z_i$.

Formalization targets

Goal: Proposition 1 (p. 11)

For either direction (HKM or NT), every interior iterate, every $\sigma$, and every choice of $F_i \in \mathcal G_i$,

$$(\Delta x, \Delta y, \Delta z) \text{ solves the Newton system for } (A, b, c; x, y, z) \iff (F^\top \Delta x,\ \Delta y,\ F^{-1} \Delta z) \text{ solves it for } (\hat A, \hat b, \hat c; \hat x, \hat y, \hat z).$$

Milestones

In attack order: the identities (31) for $F$ with $F^\top \bar J F = \bar J$; the Schur complement reduction (10)–(15); $\operatorname{Arw}(G_i^{-1} z_i) = I$ for HKM; the HKM closed forms (23) and (24)–(25); the right-hand-side formula (28) for both directions; the HKM invariance step of the proof; the NT identities $G_i x_i = G_i^{-1} z_i$ and (29); the NT block formula (30); and the NT invariance step of the proof.

Significance

The proposition shows that SDPT3's two second-order-cone directions do not depend on how each cone is coordinatized. The closed forms (23), (24), (29) and (30) are what SDPT3 uses to assemble the Schur complement matrix: they show that each block is a diagonal matrix plus a low-rank term, which is the basis for its sparse and Sherman–Morrison–Woodbury handling of large cones.

The result is proved twice in the literature (Tsuchiya 1999 and this paper). As far as is known, neither proof has been machine-checked, and no formal library contains the HKM or NT scalings, the arrow matrix of the Jordan algebra of the Lorentz cone, or its automorphism group. This mission produces those definitions and a checked proof. The printed proof assumes $F_i^\top \bar J F_i = \bar J$, that is $\lambda = 1$; the goal covers the whole group $\mathcal G_i$, so a complete formalization also closes the case $\lambda \ne 1$, where (31) and the remark that $\omega_i$ is unchanged no longer hold as printed.

Difficulty

Every step is linear algebra, but none is a one-line rewrite. The inverse of $\mathcal E_i$ appears everywhere, so one first needs invertibility of $G_i$ and $\mathcal E_i$ for interior iterates, and the HKM and NT scalings must be shown to satisfy $G_i e_1 = z_i$ and $G_i x_i = G_i^{-1} z_i$ respectively, which are identities involving $\gamma$, square roots and the denominators $\gamma(z) + z^0$ and $1 + t^0$. Transporting the Newton system directly through $F$ fails: $\mathcal E_i$ and $\mathcal F_i$ are built from $G_i$, which is not equivariant in an obvious way, and the HKM scaling at $\hat z_i$ is not $F_i^{-1} G_i$. For $\lambda \ne 1$, the quantities $\gamma$, $\omega$ and $t$ scale by powers of $\lambda$, and these factors must be tracked until they cancel.

Formalization scope

A block of dimension $q_i$ is Fin (d i + 1) → ℝ, so $q_i = d_i + 1 \ge 1$; $x_i^0$ is index 0 and $x_i^1$ the indices Fin.succ j. Mathlib's norm on this type is the sup norm, so every Euclidean quantity is written through $\bar J$. Interior points are encoded by $x^0 > 0$ and $(x^0)^2 - \langle x^1, x^1 \rangle > 0$. Matrix inverses are Matrix.inv, which is $0$ on singular input; the milestones that invert a matrix state its invertibility.

Standing assumptions: the problem is a pure SOCP ($n_s = n_l = 0$); the iterate is interior; $\mu = \langle x, z \rangle / n_q$. The paper's assumption that $A$ has full row rank is dropped, because the goal is stated for solution sets of the Newton system and needs no uniqueness. $\mathcal G_i$ is encoded algebraically (positive multiples of orthochronous Lorentz transformations) rather than as $\{F : F K = K\}$, so that the goal does not require the classification of cone automorphisms. The goal assumes nothing about the scaled iterate; its interiority is a consequence.

The Newton system is defined by its raw rows with $\mathcal E_i$, $\mathcal F_i$ from (8) and $G_i$ from (19)/(21). A formalization that defined the direction through the closed forms (23)/(29), or defined $G_i$ implicitly by $G_i e_1 = z_i$ or $G_i x_i = G_i^{-1} z_i$, would make the goal nearly definitional and is ruled out.

Reusable infrastructure: the arrow matrix and $\gamma$ of the Lorentz cone, the HKM and NT scalings, and the algebra of Lorentz transformations. Contributions toward any milestone, and alternative proofs of the goal that bypass the milestones, are welcome.

Selected references

• R. H. Tütüncü, K. C. Toh, M. J. Todd, Solving semidefinite-quadratic-linear programs using SDPT3, Mathematical Programming 95 (2003) 189–217. https://doi.org/10.1007/s10107-002-0347-5
• T. Tsuchiya, A convergence analysis of the scaling-invariant primal-dual path-following algorithms for second-order cone programming, Optimization Methods and Software 11 & 12 (1999) 141–182 (reference [27] of the SDPT3 paper; no stable link recorded).
• F. Alizadeh, D. Goldfarb, Second-order cone programming, Mathematical Programming 95 (2003) 3–51. https://doi.org/10.1007/s10107-002-0339-5
• Yu. Nesterov, M. J. Todd, Self-scaled barriers and interior-point methods in convex programming, Mathematics of Operations Research 22 (1997) 1–42. https://doi.org/10.1287/moor.22.1.1

Motivation

Obstacle problems, contact problems and control-constrained optimal control problems lead, after discretization, to complementarity systems of the form "find $y \le \psi$ and $\lambda \ge 0$ with $Ay + \lambda = f$ and $\lambda_i(y_i - \psi_i) = 0$". The primal-dual active set strategy solves such systems by guessing which constraints are active, solving one linear system, and updating the guess. Each iteration costs a single linear solve, and in practice the method stops after few iterations.

Hintermüller, Ito and Kunisch (SIAM J. Optim. 13(3), 2002; authors' version hal-01660511v1) showed that the strategy is a semismooth Newton method for a reformulation of the complementarity system through a max-function, which gives local superlinear convergence. For convergence from arbitrary initial data the paper gives three sufficient conditions on the matrix. The first is that $A$ is an M-matrix (Theorem 3.2), a class that many discretized second-order elliptic operators belong to. The M-matrix property is not stable under small perturbations: off-diagonal entries may become positive (p. 10). Theorem 3.4, the subject of this mission, shows that global convergence nevertheless survives small perturbations of an M-matrix. The paper notes that this also makes it relevant for numerical implementations.

Setting

Fix $n \ge 0$, vectors $f, \psi \in \mathbb{R}^n$ and a parameter $c > 0$. Vectors are ordered componentwise. For a matrix $A \in \mathbb{R}^{n\times n}$ the problem is

$$Ay + \lambda = f, \qquad \lambda - \max\big(0, \lambda + c(y - \psi)\big) = 0, \tag{3.1}$$

with the maximum taken componentwise. The paper writes $A$ for the matrix and calligraphic $\mathcal{A}$, $\mathcal{I}$ for the active and inactive index sets.

The primal-dual active set algorithm starts from an arbitrary pair $(y^0, \lambda^0)$. Given $(y^k, \lambda^k)$ it sets

$$\mathcal{A}_k = \{ i : \lambda^k_i + c(y^k_i - \psi_i) > 0 \}, \qquad \mathcal{I}_k = \{ i : \lambda^k_i + c(y^k_i - \psi_i) \le 0 \},$$

and determines $(y^{k+1}, \lambda^{k+1})$ from

$$Ay^{k+1} + \lambda^{k+1} = f, \qquad y^{k+1} = \psi \text{ on } \mathcal{A}_k, \qquad \lambda^{k+1} = 0 \text{ on } \mathcal{I}_k .$$

For an index set $\mathcal{I} \subseteq \{1,\dots,n\}$ with complement $\mathcal{A}$, $A_{\mathcal{I}}$ is the principal submatrix of $A$ on $\mathcal{I}$ and $A_{\mathcal{I}\mathcal{A}}$ the block with rows in $\mathcal{I}$ and columns in $\mathcal{A}$. The step above is uniquely solvable exactly when $A_{\mathcal{I}_k}$ is invertible.

A square matrix $M$ is an M-matrix if it is nonsingular, $m_{ij} \le 0$ for $i \ne j$, and $M^{-1} \ge 0$ entrywise. For a rectangular matrix $B$, $\|B\|_1$ is the matrix norm subordinate to the one-norms on both spaces, i.e. the maximal absolute column sum $\|B\|_1 = \max_j \sum_i |b_{ij}|$.

Formalization targets

Goal: Theorem 3.4

Let $M$ be an M-matrix. There is $\varepsilon > 0$ such that for every $K$ with $\|K\|_1 < \varepsilon$, every $f, \psi$ and every $c > 0$, with $A = M + K$:

1. (3.1) has exactly one solution $x^* = (y^*, \lambda^*)$;
2. the algorithm is well-defined: from every $(y^k, \lambda^k)$ the step has exactly one solution;
3. for arbitrary initial data,

$$\lim_{k\to\infty} (y^k, \lambda^k) = (y^*, \lambda^*).$$

The threshold $\varepsilon$ is left unspecified; the theorem asserts only its existence.

Milestones (the steps of the paper's proof)

• Principal submatrices $M_{\mathcal{I}}$ of an M-matrix are nonsingular M-matrices.
• With $\rho \ge \|M_{\mathcal{I}}^{-1}K_{\mathcal{I}}\|_1$ for every $\mathcal{I}$ and $\rho < \tfrac12$, every $A_{\mathcal{I}}$ is invertible and

$$A_{\mathcal{I}}^{-1} = \Big(I_{\mathcal{I}} + \sum_{i=1}^{\infty}\big(-M_{\mathcal{I}}^{-1}K_{\mathcal{I}}\big)^i\Big) M_{\mathcal{I}}^{-1}.$$

• Under the same hypothesis, for $v \le 0$, $\sum_{i\in\mathcal{I}} \big(-A_{\mathcal{I}}^{-1}v\big)_i \ge \frac{1-2\rho}{1-\rho}\,\|M_{\mathcal{I}}^{-1}v\|_1$.
• Whenever $A_{\mathcal{I}}$ is invertible, entrywise $A_{\mathcal{I}}^{-1}A_{\mathcal{I}\mathcal{A}} \le M_{\mathcal{I}}^{-1}K_{\mathcal{I}\mathcal{A}} - M_{\mathcal{I}}^{-1}K_{\mathcal{I}}(M+K)_{\mathcal{I}}^{-1}A_{\mathcal{I}\mathcal{A}}$.
• For $\|K\|_1$ small (uniformly in $f, \psi, c$) and $k \ge 1$: $\sum_i y^{k+1}_i < \sum_i y^k_i$ unless $y^{k+1} = y^k$, so $\mathcal{M}(y) = \sum_i y_i$ is a merit function.

A companion statement records the paper's remark after the theorem: for $\|K\|_1$ small, $M + K$ is a P-matrix.

Significance

The result. Theorem 3.2 makes global convergence depend on a sign pattern, which a perturbation of any size can destroy. Theorem 3.4 replaces the sign pattern by a neighbourhood: every matrix close enough to an M-matrix in $\|\cdot\|_1$ inherits unique solvability of (3.1), well-posedness of every iteration, and convergence from every starting point. As the paper observes, solvability of (3.1) for every $f$ is equivalent to $A$ being a P-matrix (Berman–Plemmons, Thm. 10.2.15), so the theorem also yields a sufficient condition for the P-matrix property.

Formalizing it. The paper's argument is a short sketch that refers back to the proof of Theorem 3.3 and asserts that "$\|K\|_1$ can be chosen sufficiently small". A formalization makes the uniformity of that choice explicit (it depends only on $M$), checks the strict decrease of the merit function including its equality case, and fixes the role of the first iterate, from which the sign facts used in (3.8) hold. To our knowledge none of the three global convergence theorems of the paper has a machine-checked proof.

Difficulty

The obvious argument is to transfer Theorem 3.2 by continuity. It fails because the M-matrix proof is monotonicity-based: it shows $y^* \le y^{k+1} \le y^k$, and that ordering uses $A_{\mathcal{I}}^{-1} \ge 0$ and $A_{\mathcal{I}\mathcal{A}} \le 0$, sign conditions that a perturbation breaks immediately. What survives is weaker: a scalar merit function $\sum_i y_i$ whose decrease must be quantified. The perturbation has to be controlled uniformly over all $2^n$ index sets, both in the multiplier term (where the decrease must remain strict whenever $\lambda^k_{\mathcal{I}} \ne 0$) and in the coupling term $A_{\mathcal{I}}^{-1}A_{\mathcal{I}\mathcal{A}}(y^k-\psi)_{\mathcal{A}}$, whose positive part must be absorbed by $-\sum_{\mathcal{A}}(y^k_i - \psi_i)$. Finally, strict decrease of a merit function does not by itself give convergence; the finiteness of the possible active sets must be combined with it.

Formalization scope

Vectors are Fin n → ℝ and matrices Matrix (Fin n) (Fin n) ℝ; index sets are Finset (Fin n) and blocks are submatrices along subtype inclusions. Conventions committed to:

• $\|B\|_1$ is the maximal absolute column sum, computed as a supremum in $\mathbb{R}_{\ge0}$, equal to $0$ for a matrix with no columns; it is not the entrywise sum or another operator norm;
• the M-matrix predicate requires nonsingularity explicitly (IsUnit of the determinant), because Lean's inverse of a singular matrix is $0$; it is stated for any finite index type so that it applies to principal blocks;
• "$\|K\|_1$ sufficiently small" is one $\varepsilon > 0$ depending only on $M$ (and $n$), quantified before $f$, $\psi$, $c$ and the initial data;
• a step is a relation and a run is an infinite sequence of steps from arbitrary initial data; step (iv)'s "Stop" is not modelled;
• the convergence conclusion is stated for every solution of (3.1), which by the first conclusion is the solution; convergence is in the product topology;
• the merit-function milestone holds from $k = 1$ on, as in the paper.

The goal does not assume that $M + K$ is a P-matrix, that its principal blocks are invertible, that $\rho < \tfrac12$, or that the merit function decreases: all of these are consequences. A version with $\varepsilon$ chosen after $f$, $\psi$, $c$ or the run, or with any of these as hypotheses, would be trivial or different and is ruled out. Well-definedness is a conclusion, so runs from every initial pair exist and the convergence conclusion is not vacuous.

In the g-estimate, the page's final $\|M^{-1}\lambda^k_{\mathcal{I}}\|_1$ is a misprint for $\|M_{\mathcal{I}}^{-1}\lambda^k_{\mathcal{I}}\|_1$; the statement uses $M_{\mathcal{I}}^{-1}$.

A complete development needs: principal submatrices of M-matrices, Neumann series for matrices with the column-sum norm, the block form of one step and the identity (3.8), the sign facts $\lambda^k \le 0$ on $\mathcal{I}_k$ and $y^k \ge \psi$ on $\mathcal{A}_k$ for $k \ge 1$, and continuity of principal minors for the P-matrix statement. The M-matrix and Neumann-series lemmas are reusable beyond this mission. Proofs of any milestone, of the companion, or of the goal through another argument are welcome.

Selected references

• M. Hintermüller, K. Ito, K. Kunisch, The primal-dual active set strategy as a semismooth Newton method, SIAM J. Optim. 13(3):865–888, 2002. https://doi.org/10.1137/S1052623401383558 ; authors' version https://hal.science/hal-01660511v1
• A. Berman, R. J. Plemmons, Nonnegative Matrices in the Mathematical Sciences, Academic Press, 1979 (SIAM Classics reprint 1994). https://doi.org/10.1137/1.9781611971262
• R. H. W. Hoppe, Multigrid algorithms for variational inequalities, SIAM J. Numer. Anal. 24:1046–1065, 1987. https://doi.org/10.1137/0724069

Motivation

Airlines, hotels and car-rental firms sell reservations ahead of service, and a fraction of those reservations cancel or do not show up. Overbooking — accepting more reservations than there is capacity — recovers revenue that cancellations would otherwise destroy, at the risk of denying service to some customers. In practice inventory is often substitutable: a passenger booked in economy can be upgraded to business, a guest booked in a standard room can be given a suite. With substitution, the right overbooking level for one reservation class depends on the levels of all the others.

Karaesmen and van Ryzin (Operations Research 52(1):83–104, 2004) model this as a two-period problem: overbooking decisions in a reservation period, followed by an optimal assignment of the customers who actually show up to inventory classes in a service period. Their structural result, Theorem 2, establishes the shape of the expected net revenue as a function of the overbooking levels. That shape is what justifies critical-level booking policies and the stochastic-gradient method of their §4, and it is the textbook account of joint overbooking with substitution (Talluri and van Ryzin, The Theory and Practice of Revenue Management, 2004, Proposition 4.4).

Setting

There are $n$ reservation classes $i = 1, \dots, n$ and $m + 1$ inventory classes $j = 0, 1, \dots, m$. Class $j = 0$ is virtual: assigning a customer to it means the customer is denied service.

Reservation period. The firm holds $x_i$ reservations of class $i$ and chooses overbooking levels $u_i \ge x_i$. Accepting a class-$i$ reservation earns $r_i$; a cancellation of class $i$ is refunded $q_i$, with $q_i \le r_i$. The number of class-$i$ customers who survive to the service period is a random variable $Z_i(u_i)$ with values in $\{0, 1, 2, \dots\}$. The family $\{Z_i(u_i)\}$ has the semigroup property: if $Y_1$ and $Y_2$ are independent with the laws of $Z_i(u_i)$ and $Z_i(s_i)$, then $Y_1 + Y_2$ has the law of $Z_i(u_i + s_i)$. The binomial model ($Z_i(u_i) \sim \mathrm{Bin}(u_i, p_i)$ for integer $u_i$) and the Poisson model ($Z_i(u_i) \sim \mathrm{Poisson}(u_i p_i)$) both have it. The survivals of different classes are independent.

Service period. Given the survivals $z = (z_1, \dots, z_n) \ge 0$, net benefits $a_{ij}$ (of either sign) and capacities $c_j$, the firm solves the transportation problem

$$V_0(z, c) = \max \sum_{i=1}^n \sum_{j=0}^m a_{ij} y_{ij} \quad \text{s.t.} \quad \sum_{j=0}^m y_{ij} = z_i, \ \ \sum_{i=1}^n y_{ij} \le c_j, \ \ y_{ij} \ge 0.$$

Expected net revenue. The objective of the reservation period is

$$G(u) = \sum_{i=1}^n r_i (u_i - x_i) - \mathbb E\Big[\sum_{i=1}^n q_i\big(u_i - Z_i(u_i)\big)\Big] + \mathbb E\big[V_0(Z(u), c)\big],$$

to be maximized over $u \ge x$.

A function $f$ on $\mathbb R^n$ is submodular on $S$ if $f(s \vee s') + f(s \wedge s') \le f(s) + f(s')$ for $s, s' \in S$, with $\vee, \wedge$ the componentwise maximum and minimum.

Formalization targets

Goal: Theorem 2

For every $u$ with levels in the parameter set $P$, every class $i$ and all $\varepsilon, \alpha \in P$,

$$G(u + (\varepsilon + \alpha) e_i) - G(u + \alpha e_i) \le G(u + \varepsilon e_i) - G(u),$$

and $G$ is submodular on $P^n$:

$$G(u \vee u') + G(u \wedge u') \le G(u) + G(u').$$

The first statement is the paper's componentwise concavity (6); the second is its submodularity, proved in the paper through the pairwise inequality (7) and Theorem 1.

Milestones

1. LP duality for the transportation problem (§2.2, p. 86): for $z \ge 0$ the primal attains its optimum and equals the dual optimum.
2. Lemma 1: $V_0(\cdot, c)$ is submodular on $\mathbb R^n_{\ge 0}$.
3. Lemma 2: $V_0$ is concave in $(z_1, \dots, z_n)$ and concave in $(c_1, \dots, c_m)$.
4. Theorem 1: on a product set, submodularity is equivalent to decreasing differences.
5. Under the semigroup property, $\mathbb E[Z_i(u_i)]$ is additive in $u_i$.
6. The semigroup representation $\mathbb E[V_0(Z(u + \varepsilon e_i))] = \mathbb E[V_0(Z(u) + Y e_i)]$ with $Y$ independent of $Z(u)$ and distributed as $Z_i(\varepsilon)$.
7. The sample-path inequalities of the proof: nonincreasing increments of $V_0$ along one coordinate, and the pairwise submodularity inequality for two coordinates.

Significance

Componentwise concavity of $G$ implies that for each reservation class there is a critical booking level beyond which expected revenue does not increase, other levels held fixed. Submodularity implies that these critical levels are nonincreasing in the bookings accepted for any other class: heavier overbooking in one class leaves less spare substitute capacity for the others. These two properties underpin the heuristics and the stochastic-gradient optimization of §4–§5 of the paper.

The result is proved in the paper (for $n = 2$, with the general case described as straightforward). A machine-checked version would cover the general semigroup case, which includes the binomial model the paper calls realistic, real-valued levels and arbitrary step sizes. The platform already holds a proved special case, RevenueManagement.joint_overbooking_concave_submodular (Talluri–van Ryzin, Proposition 4.4): Poisson survivals, integer levels and unit steps only. This mission generalizes it in all three respects and also formalizes the transportation-problem facts (LP duality, submodularity and concavity of $V_0$) as reusable statements.

Difficulty

The obvious route — differentiate $G$ and use concavity of $V_0$ — fails twice. First, $V_0$ is piecewise linear and the survivals are discrete, so there is nothing to differentiate; under the binomial model $G$ is only defined at integer levels. Second, concavity of $V_0$ alone does not survive the expectation: $\mathbb E[V_0(Z(u))]$ is a mixture over laws that change with $u$, and a concave function of a random vector whose law moves with the parameter need not be concave in the parameter. The semigroup property is what replaces a change of law by the addition of an independent block of survivors, and that coupling must be built for every coordinate and every pair of coordinates at once. On the LP side, submodularity of a transportation-problem value in its demands is a lattice-programming fact (Topkis, Theorem 3.4.1; Shapley 1962), not a consequence of concavity.

Formalization scope

All declarations are in the namespace SubstOverbooking.Structure. Reservation classes are Fin n (0-based), inventory classes Fin (m + 1) with 0 the virtual denied-service class. Conventions committed to:

• Uncapacitated virtual class. The paper gives class 0 "finite but very high capacity" and then needs it to absorb every survivor; with Poisson survivals no finite capacity does. Constraint (4) is imposed only for $j \ne 0$, $c_0$ is unused, and the dual fixes $\lambda_0 = 0$. Capacities satisfy $c_j \ge 0$ for $j \ne 0$.
• $V_0$ on real demands. $V_0$ is the supremum of the TP objective over feasible assignments, for real $z$; every statement assumes $z \ge 0$, where the supremum is attained and finite.
• Survival laws. Each class has a family $u \mapsto L_i(u)$ of probability mass functions on $\mathbb N$, with finite means. Levels range over a parameter set $P \subseteq [0, \infty)$ that contains $0$ and is closed under addition and under nonnegative differences: $P = \mathbb N$ covers the binomial model, $P = [0, \infty)$ the Poisson model. The semigroup property is the identity of laws $L_i(u + s) = (L_i(u) \otimes L_i(s)) \circ (+)^{-1}$ for $u, s \in P$, and independence across classes is the product law.
• Concavity is stated as the increment inequality (6), which is meaningful on a discrete $P$; submodularity of $G$ is lattice submodularity on $P^n$ via the published SupermodularOn applied to $-G$.
• Theorem 1 is stated for product sets $S = T_1 \times \dots \times T_n$; as printed for arbitrary $S \subset \mathbb R^n$ it is false, and every set the paper applies it to is a product.
• The refund condition $q_i \le r_i$ is kept as a hypothesis of the goal, as on the page.

A formalization that fixes Poisson laws, integer levels or unit steps, or that defines $G$ or $V_0$ through the inequalities to be proved, is the already-proved special case or an empty statement, and is not the goal.

A complete development needs finite LP duality for the transportation problem, the lattice-programming argument for Lemma 1, product measures and image measures on $\mathbb N^n$, and the coupling of the semigroup property. The LP and lattice parts are reusable beyond this mission. Proofs of any milestone are welcome, as are alternative arguments for Lemma 1 (for instance via monotonicity of optimal duals, after Erlenkotter 1970).

Selected references

• I. Karaesmen and G. van Ryzin, Overbooking with Substitutable Inventory Classes, Operations Research 52(1):83–104, 2004. https://doi.org/10.1287/opre.1030.0079
• K. T. Talluri and G. J. van Ryzin, The Theory and Practice of Revenue Management, Springer, 2004, §4.4. https://doi.org/10.1007/b139000
• D. M. Topkis, Supermodularity and Complementarity, Princeton University Press, 1998, Theorem 3.4.1. https://doi.org/10.1515/9781400822539
• R. K. Sundaram, A First Course in Optimization Theory, Cambridge University Press, 1996. https://doi.org/10.1017/CBO9780511804526
• L. S. Shapley, Complements and substitutes in the optimal assignment problem, Naval Research Logistics Quarterly 9(1):45–48, 1962. https://doi.org/10.1002/nav.3800090106
• M. Shaked and J. G. Shanthikumar, Stochastic convexity and its applications, Advances in Applied Probability 20(2):427–446, 1988. https://doi.org/10.2307/1427398

Why the Nesterov–Todd scaling of a second-order cone matters

Conic quadratic (second-order cone) programs are solved in practice by primal-dual interior-point methods. Each iteration solves a Newton system for the central-path equations, and for non-polyhedral cones that system is only guaranteed to be well posed after the problem is scaled: the primal point $x$ and the dual slack $s$ are transformed by a matrix $\Theta W$ chosen from the current iterate. Among the scalings analysed by Monteiro and Tsuchiya (2000), the best complexity bounds are obtained with the Nesterov–Todd (NT) scaling of Nesterov and Todd (1997), which maps $x$ and $s$ to the same point.

Andersen, Roos and Terlaky (2003) describe the conic quadratic interior-point solver of MOSEK. Their Lemma 4.2 gives the NT scaling in closed form for each of the three cone types they use. It is what makes the NT direction cheap: each cone block of $W$ is determined by one vector $w$, and $W$, $W^2$ and their inverses can be applied in time linear in the dimension. For the quadratic cone the formula goes back to Tsuchiya (1999); the paper re-derives it and extends it to the rotated quadratic cone.

Setting

Fix one cone block of dimension $d$. Coordinates are $x = (x_1,\dots,x_d)$, $x_{k:d}$ denotes the subvector of coordinates $k,\dots,d$, and $\|\cdot\|$ is the Euclidean norm. The block is one of

• the positive half-line $\mathbb R_+$ ($d=1$), with $T = Q = 1$;
• the quadratic cone $K^q = \{x : x_1 \ge \|x_{2:d}\|\}$ ($d\ge1$), with $T = I$ and $Q = \operatorname{diag}(1,-1,\dots,-1)$;
• the rotated quadratic cone $K^r = \{x : 2x_1x_2 \ge \|x_{3:d}\|^2,\ x_1,x_2\ge0\}$ ($d\ge2$), with $T$ the identity except for the leading block $\tfrac1{\sqrt2}\begin{bmatrix}1&1\\1&-1\end{bmatrix}$, and $Q$ the matrix with $Q_{12} = Q_{21} = 1$, $Q_{jj} = -1$ for $j\ge3$, and zeros elsewhere.

A scaling matrix (Definition 3.3) is a symmetric positive definite $W$ with $WQW = Q$. For $\theta>0$ the scaled points are $\bar x = \theta W x$ and $\bar s = (\theta W)^{-1}s$. The pair $(\theta, W)$ is an NT scaling of $(x, s)$ when $\bar x = \bar s$, equation (28); equivalently $s = \theta^2W^2x$, equation (29).

Formalization targets

Goal: Lemma 4.2 (p. 15)

For $x, s$ in the interior of the cone:

1. every NT scaling satisfies

$$\theta^2 = \sqrt{\frac{s^TQs}{x^TQx}};\qquad(30)$$

2. with this $\theta$ and

$$w = \frac{\theta^{-1}s + \theta Qx}{\sqrt2\sqrt{x^Ts + \sqrt{x^TQx\,s^TQs}}},\qquad(32)$$

the following matrices $W$ make $(\theta, W)$ an NT scaling: $W = \theta^{-1}(X^{-1}S)^{1/2}$ for $\mathbb R_+$;

$$W = \begin{bmatrix} w_1 & w_{2:d}^T\\ w_{2:d} & I + \frac{w_{2:d}w_{2:d}^T}{1+w_1}\end{bmatrix} = -Q + \frac{(e_1+w)(e_1+w)^T}{1+e_1^Tw}\qquad(31)$$

for $K^q$; and

$$W = -Q + \frac{(Te_1+w)(Te_1+w)^T}{1+e_1^TTw}\qquad(34)$$

for $K^r$. In the last two cases $W^2 = -Q + 2ww^T$, equations (33) and (35).

Milestones

Lemma 3.3 i)–iii) (scaling preserves $x^Ts$ and rescales $x^TQx$, $s^TQs$ by $\theta^{\pm2}$); (30) on its own; the normalization $w_1^2 - \|w_{2:d}\|^2 = 1$ and the square (62) of the arrow form (31); positivity of the radicand of (32); the identity $2\theta^2(w^Tx)^2 = x^Ts + \sqrt{x^TQx\,s^TQs}$ and the formula (32); for the rotated cone, the transport (63) by $T$, the identities $\hat\theta = \theta$ and $\hat w = Tw$, and $TWT = \hat W$. Two companions: Lemma 4.3, $(\theta W)^{-2} = \theta^{-2}QW^2Q$, and Lemma 3.3 iv), that a scaling maps the cone and its interior onto themselves.

Significance

The result. Lemma 4.2 turns the NT scaling, defined implicitly by (28), into explicit data: one scalar $\theta$ and one vector $w$ per cone. Combined with Lemma 4.3 it shows that every product with $W$, $W^2$ or their inverses in the scaled Newton system (27) is a rank-one correction of a signature matrix, which is what the paper's implementation relies on. The rotated-cone formula (34) is the paper's own addition and is what lets rotated cones be handled natively instead of through a change of variables.

Formalizing it. The result is proved on paper; no machine-checked version is known to exist, and the platform has no formalization of second-order-cone scalings. The work is the matrix algebra of the Lorentz group — the signature matrix $Q$, the matrices with $WQW = Q$, and their action on $K^q$ and $K^r$ — which is reusable for every second-order-cone algorithm.

Difficulty

Most of the items are finite-dimensional identities, but they are identities between block matrices of arbitrary size with entries that are nested square roots, and they hold only on the interior of the cone, where the radicands are positive. The derivation in the paper runs backwards: it assumes $W$ has the shape (31) and solves for $w$, taking a square root without fixing its sign. The goal asks for the forward direction — that the displayed $\theta$, $w$, $W$ do satisfy (28) and $W\succ0$ — which needs the normalization $w^TQw = 1$, positive definiteness of $-Q + (e_1+w)(e_1+w)^T/(1+w_1)$, and, for the rotated cone, the reduction to the quadratic cone through $T$ with $TT = I$ and $TQ^rT = Q^q$.

Formalization scope

• Per block. Every statement is made for one cone block. $\Theta$ and $W$ are block diagonal (pp. 11–12), so this is the paper's statement block by block. The cone type is the parameter kind : ConeKind, and vectors are Fin d → ℝ.
• Index base and norms. The paper's $x_1$ is coordinate 0. $\|x_{2:d}\|^2$ is written as a sum of squares (tailSq x 1); Mathlib's norm on Fin d → ℝ, which is the sup norm, is never used.
• Interiors. The interiors are given by strict inequalities. The paper never defines $\operatorname{int}(K)$ concretely; these coincide with the topological interiors.
• Added hypotheses. All are disclosed in the items:
  • the dimension conventions ($d=1$ for $\mathbb R_+$, $d\ge1$ for $K^q$, $d\ge2$ for $K^r$), which the paper leaves implicit;
  • $\theta>0$ in Definition 3.3, which the paper calls positive definite on p. 24;
  • $1+w_1>0$ wherever the page divides by $1+w_1$;
  • $w^Tx>0$ for the last step of the derivation of (32), where the page silently fixes this sign.
• Inverses. The matrix inverse is Mathlib's total inverse. The NT predicate always contains $W\succ0$, so (28) never reads the junk value of a singular inverse.
• Reading of Lemma 4.2. The lemma names $\theta$ and $W$ without quantifying them. The goal states (30) for every NT scaling and asserts that the displayed formulas are an NT scaling. It does not claim that the NT scaling is unique, which the paper does not prove.
• Ruling out a trivial formalization. Stating only (30) would allow a version that holds vacuously if no NT scaling existed. Part 2 of the goal is what excludes this, and a sanity check exhibits $\theta = 1$, $W = I$ at $x = s = e_1$ in $K^q$ for $d=2$.

A complete development needs:

• the cone definitions and the matrices $T$, $Q$;
• Mathlib's Matrix.PosDef and Matrix.inv;
• elementary real square-root algebra.

Lemmas about $\{W : W\succ0,\ WQW = Q\}$ (closure under products and inverses, invariance of $K^q$) are reusable beyond this mission. Proofs of any milestone are welcome, as are alternative arguments for positive definiteness of (31).

Selected references

• E. D. Andersen, C. Roos, T. Terlaky, On implementing a primal-dual interior-point method for conic quadratic optimization, Mathematical Programming 95 (2003) 249–277; authors' preprint of 18 Dec 2000. https://doi.org/10.1007/s10107-002-0349-3
• Yu. E. Nesterov, M. J. Todd, Self-scaled barriers and interior-point methods for convex programming, Mathematics of Operations Research 22(1) (1997) 1–42. https://doi.org/10.1287/moor.22.1.1
• R. D. C. Monteiro, T. Tsuchiya, Polynomial convergence of primal-dual algorithms for the second-order cone program based on the MZ-family of directions, Mathematical Programming 88(1) (2000) 61–83. https://doi.org/10.1007/PL00011378
• T. Tsuchiya, A convergence analysis of the scaling-invariant primal-dual path-following algorithms for second-order cone programming, Optimization Methods and Software 11–12 (1999) 141–182. https://doi.org/10.1080/10556789908805750
• F. Alizadeh, D. Goldfarb, Second-order cone programming, Mathematical Programming 95 (2003) 3–51. https://doi.org/10.1007/s10107-002-0339-5