All missions

Motivation

Online retailers show their catalogue on a sequence of pages, and most visitors never reach the last one. Which products appear on the first page, and at what prices, therefore decides what a typical consumer can buy at all. Gallego, Li, Truong and Wang (Operations Research, 2020) model this as the product framing problem: the consideration set of a consumer is the set of products on the pages that consumer is willing to view, and the number of pages viewed is random. Their Section 7 adds prices to the problem and asks how a retailer should frame and price jointly when consumers choose by the multinomial logit (MNL) model.

The question connects two literatures. In assortment pricing under MNL with a common price sensitivity, the optimal prices of all offered products are known to be equal (Anderson, de Palma and Thisse 1992; Gallego and Stefanescu, as cited on p. 16 of the paper), and the nested-logit extension of Gallego and Wang (Operations Research, 2014) keeps an adjusted markup constant across nests. In search and ranking models, the order in which products are shown changes which ones are bought. The paper's structural results say what survives of the equal-price rule once consideration sets are nested by page, and which products belong on the early pages.

Setting

There are $n$ products, $i \in [n] = \{1,\dots,n\}$, and $m$ pages, each holding at most $p$ products. A framing places each product on one page $x(i)\in[m]$ or leaves it undisplayed. A consumer views the first $X$ pages, where $X\in[m]$ is random with law $\lambda(x)=\mathbf P[X=x]$, independent of the framing. A consumer who views $x$ pages considers $S(x)$, the products on pages $1,\dots,x$, so $S(1)\subseteq S(2)\subseteq\cdots\subseteq S(m)$.

Product $i$ has a quality $a_i\in\mathbb R$ and a price $r_i\in\mathbb R$, and the price sensitivity is $\beta>0$. The mean utility of product $i$ is $u_i=a_i-\beta r_i$, and the outside option has utility $u_0=0$. Under the MNL model a consumer with consideration set $S$ buys $i\in S$ with probability

$$P(i,S)=\frac{e^{u_i}}{1+\sum_{k\in S}e^{u_k}},$$

and $P(i,S)=0$ for $i\notin S$. The expected revenue from that consumer is $R(r\mid S)=\sum_{i\in S}r_iP(i,S)$, and the total expected revenue of a framing and a price vector is

$$\mathbf E\big[R(r\mid S(X))\big]=\sum_{x=1}^m\lambda(x)\,R\big(r\mid S(x)\big).$$

For a fixed framing, $R(a)$ denotes the optimal value of this quantity over all price vectors, as a function of the quality vector $a$.

Formalization targets

Goal: Theorem 6 (p. 17)

The goal is the existence of an optimal joint solution with the paper's structure: there are a feasible framing and prices $r\in\mathbb R^n$ such that no feasible framing with any prices earns more, and

$$|S(x)|=\min(n,\;x\,p)\quad\text{for all }x\in[m],\qquad a_k\le a_i\ \text{ whenever } i \text{ is displayed and } k \text{ is undisplayed or on a later page.}$$

Pages are filled in order until all products are displayed, and the products appear in descending order of quality.

Milestones

1. (15), p. 43. The partial derivative of the total expected revenue in the price of a displayed product:

$$\frac{\partial\,\mathbf E[R(r\mid S(X))]}{\partial r_i}=\beta\sum_{l=x(i)}^m\lambda(l)P\big(i,S(l)\big)\Big\{\tfrac1\beta+R\big(r\mid S(l)\big)-r_i\Big\}.$$

2. (17), p. 43. At every optimal price vector, each displayed price is $1/\beta$ plus a weighted average of the revenues $R(r\mid S(l))$, $l\ge x(i)$.
3. Theorem 5, p. 16. For a fixed framing, every optimal price vector is constant on each page, $r_i=\theta_{x(i)}$, with $\theta_1\le\theta_2\le\cdots\le\theta_m$.
4. (19), p. 44. The optimal value $R(a)$ is nondecreasing in the quality of any displayed product.

Significance

Theorem 5 reduces the pricing problem of a fixed framing from $n$ prices to $m$ page prices and says that later pages carry higher prices; this is the page-level analogue of the constant-markup property of MNL pricing, and it runs opposite to the ordering found in oligopoly search models (Arbatskaya 2007). Theorem 6 removes the framing decision almost entirely: a joint optimum is obtained by sorting the products by quality and filling pages in that order, leaving only page prices to choose. The paper's approximation algorithm for joint pricing and framing (Theorem 7, NEST-P, with guarantee $6/\pi^2$) starts from this structure.

None of these statements has a machine-checked proof. A complete development would give the first formal treatment of MNL price optimization with nested consideration sets, including the existence of optimal prices, which the paper takes for granted. The paper's argument for the monotonicity $\theta_1\le\cdots\le\theta_m$ and its envelope-theorem step (19) are informal; a formal proof has to supply both.

Difficulty

The expected revenue is not concave in the prices, even with two products on two pages (Example 1, p. 17), so the first-order condition alone does not identify the optimum, and existence of an optimal price vector has to be argued from the behaviour of the revenue as prices tend to $\pm\infty$. The within-page equality of prices follows from the first-order condition, but the ordering of page prices needs a global comparison of revenues across consideration sets. In Theorem 6, the obvious exchange argument (swap a higher-quality product forward) does not obviously keep revenue from decreasing when the swapped products have different prices; the paper splits into two cases, and in one of them the improvement comes from continuously moving qualities, which requires control of the optimal value as a function of $a$ rather than of a fixed price vector.

Formalization scope

Products are Fin n, pages are the natural numbers $1,\dots,m$ with $m\ge1$, and capacities satisfy $p\ge1$. A framing is f : Fin n → ℕ with page f i ∈ [m] for displayed products and f i = 0 for undisplayed ones (the paper's $x(i)=m+1$); feasibility means at most $p$ products per page. The law of $X$ is a nonnegative function on $[m]$ summing to one. Prices are finite reals: the paper's priced-out products ($r_i=+\infty$) coincide with undisplayed ones, which the framing already allows. "Descending order" and "increases" are read weakly.

Theorem 5 and (17) are stated for every maximizing price vector and only for products on pages that some consumer reaches ($\Lambda(x)=\mathbf P[X\ge x]>0$), a restriction the paper never discusses: a product on a page no consumer reaches can carry any price at an optimum. They do not assert that a maximizer exists. Theorem 6 is read existentially ("some optimal solution has this structure"); the universal reading is false when some $\lambda(x)=0$, because exchanging products between pages $x$ and $x+1$ then changes no revenue. (19) is formalized as its stated consequence, the monotonicity of the optimal value $R(a)=\sup_r\mathbf E[R(r\mid S(X))]$, together with the boundedness of the revenues; the envelope identity itself would need a differentiable selection of optimal prices that the paper does not establish.

A trivializing formalization of the goal would take the optimum over framings that are already sorted, over a single fixed framing, or over prices fixed in advance; the goal's optimality clause ranges over all feasible framings and all real price vectors. The development needs elementary calculus of the MNL revenue (HasDerivAt, Real.exp), finite sums over pages, and a compactness or limiting argument for the existence of optimal prices. Pages are those of the authors' accepted manuscript, which differ from the journal typesetting. Proofs of any milestone, and lemmas on MNL pricing reusable beyond this paper (boundedness of $R(r\mid S)$, existence of optimal MNL prices, the constant-price property for a single assortment), are welcome.

Selected references

• G. Gallego, A. Li, V.-A. Truong, X. Wang, Approximation Algorithms for Product Framing and Pricing, Operations Research 68(1), 2020. https://doi.org/10.1287/opre.2019.1875
• S. P. Anderson, A. de Palma, J.-F. Thisse, Discrete Choice Theory of Product Differentiation, MIT Press, 1992 (book; cited by the paper for the equal-price property of MNL pricing).
• G. Gallego, R. Wang, Multiproduct Price Optimization and Competition under the Nested Logit Model with Product-Differentiated Price Sensitivities, Operations Research 62(2), 2014. https://doi.org/10.1287/opre.2013.1249
• M. Arbatskaya, Ordered Search, The RAND Journal of Economics 38(1), 2007. http://www.jstor.org/stable/25046295

Motivation

Minimizing a convex function $g$ on a real Hilbert space $\mathcal H$ by first-order methods has a continuous-time counterpart: second-order differential equations whose trajectories descend $g$. Su, Boyd and Candès (20) showed that Nesterov's accelerated gradient method is a discretization of $\ddot x+\frac{\alpha}{t}\dot x+\nabla g(x)=0$ with $\alpha=3$, along which $g(x(t))-\min g=O(1/t^2)$. Two modifications of this equation have been studied separately because each adds a property the plain system lacks.

• A Hessian-driven damping term $\beta\nabla^2g(x(t))\dot x(t)$ (Attouch, Peypouquet, Redont, J. Differential Equations 2016) damps oscillations and, for $\beta>0$, makes $t^2\|\nabla g(x(t))\|^2$ integrable, while keeping fast values and weak convergence of trajectories.
• A Tikhonov regularization term $\epsilon(t)x(t)$ with $\epsilon(t)\downarrow0$ (Attouch, Chbani, Riahi, J. Math. Anal. Appl. 2018) can steer trajectories to the minimum-norm minimizer.

Boţ, Csetnek and László (arXiv:1911.12845v2, 2020) combine both terms and ask which properties of the two parent systems survive. This mission covers the first of their answers: when $\epsilon$ decays fast enough, the fast rates and the weak convergence of trajectories are preserved.

Timeline. 2016 (published 2018): weak convergence of trajectories of $\ddot x+\frac\alpha t\dot x+\nabla g(x)=0$ for $\alpha>3$ (Attouch, Chbani, Peypouquet, Redont, Math. Program. 168). 2016: the same with Hessian-driven damping, $\beta>0$ (Attouch, Peypouquet, Redont). 2018: Tikhonov regularization of that system, with the $O(1/t^2)$ rate kept and strong convergence to the minimum-norm minimizer under slow decay of $\epsilon$ (Attouch, Chbani, Riahi). 2020: both terms together (Boţ, Csetnek, László, this paper).

Setting

Fix a real Hilbert space $\mathcal H$, a starting time $t_0>0$, parameters $\alpha\ge3$ and $\beta\ge0$, and initial data $u_0,v_0\in\mathcal H$. The General assumption of the paper is:

• $g:\mathcal H\to\mathbb R$ is convex and twice Fréchet differentiable, $\nabla g$ is Lipschitz continuous on bounded sets, and $\operatorname{argmin} g\ne\emptyset$; $\min g$ denotes its minimal value;
• $\epsilon:[t_0,+\infty)\to[0,+\infty)$ is nonincreasing, of class $C^1$, and $\epsilon(t)\to0$ as $t\to+\infty$.

The dynamical system is

$$\ddot x(t)+\frac{\alpha}{t}\dot x(t)+\beta\nabla^2g(x(t))\dot x(t)+\nabla g(x(t))+\epsilon(t)x(t)=0,\quad t\ge t_0,\quad x(t_0)=u_0,\ \dot x(t_0)=v_0. \tag{5}$$

A global $C^2$-solution is a twice continuously differentiable $x:[t_0,+\infty)\to\mathcal H$ satisfying (5) for every $t\ge t_0$. The decay of $\epsilon$ is controlled by $\int_{t_0}^{+\infty}t\epsilon(t)\,dt<+\infty$ together with one of two conditions:

• (a) there exist $a>1$ and $t_1\ge t_0$ with $\dot\epsilon(t)\le-\frac{a\beta}{2}\epsilon^2(t)$ for all $t\ge t_1$;
• (b) there exist $a>0$ and $t_1\ge t_0$ with $\epsilon(t)\le\frac at$ for all $t\ge t_1$.

A curve $x(t)$ converges weakly to $\bar x$ if $\langle x(t),v\rangle\to\langle\bar x,v\rangle$ for every $v\in\mathcal H$. The Lean development lives in the namespace TikhonovHDD.Weak, with IsSolution, GenAssumptionG, GenAssumptionEps, CondA, CondB, argmin, minVal and WeakTendstoAtTop for these objects.

Formalization targets

Goal: Theorem 3.5 (p. 15)

Under the General assumption, $\int_{t_0}^{+\infty}t\epsilon(t)\,dt<+\infty$ and (a) or (b), if $\alpha>3$ then every global $C^2$-solution of (5) satisfies

$$\exists\,\bar x\in\operatorname{argmin} g:\qquad x(t)\rightharpoonup\bar x\quad(t\to+\infty).$$

Milestones

In the order the argument uses them:

1. Theorem 2.1 (p. 4): existence and uniqueness of a global $C^2$-solution of (5).
2. Lemma A.2 (p. 29): a bounded-below, locally absolutely continuous $F$ with $\dot F\le G\in L^1$ a.e. has a finite limit.
3. Theorem 3.3 (p. 11), split three ways: for $\alpha\ge3$, $g(x(t))-\min g=O(1/t^2)$; for $\alpha>3$, $x$ is bounded and $t(g(x(t))-\min g)$, $t\|\dot x\|^2$, $t\epsilon\|x-x^*\|^2$, $t\epsilon\|x\|^2\in L^1$; and separately $t^2\|\nabla g(x(t))\|^2\in L^1$.
4. Lemma A.3 (p. 30): $u(t)+\frac t\alpha\dot u(t)\to\bar u$ with $\alpha>0$ implies $u(t)\to\bar u$.
5. Theorem 3.4 (p. 12): for $\alpha>3$ and every $x^*\in\operatorname{argmin} g$, $t\langle\nabla g(x),x-x^*\rangle\in L^1$, $\lim\|x(t)-x^*\|$ and $\lim t\langle\dot x+\beta\nabla g(x),x-x^*\rangle$ exist, $g(x(t))-\min g=o(1/t^2)$, $\|\dot x+\beta\nabla g(x)\|=o(1/t)$, and $t^2\epsilon(t)\|x(t)\|^2\to0$.

Significance

The result. Theorem 3.5 shows that the Tikhonov term, when it decays at least as fast as $\int t\epsilon<\infty$ requires, does not disturb the asymptotics of the Hessian-damped system: values converge at rate $o(1/t^2)$ and trajectories converge weakly to a minimizer. With $\beta=0$ it recovers the corresponding results for the regularized system without Hessian damping, and with $\epsilon\equiv0$ those for the Hessian-damped system. It is also the boundary case of the paper's second half, where slower decay of $\epsilon$ is shown to give strong convergence to the minimum-norm minimizer instead.

Formalizing it. The results are proved in the paper; none of them is machine-checked. A formal development would produce reusable pieces: a Lyapunov-energy analysis of second-order systems with time-dependent coefficients, the asymptotic Lemmas A.2 and A.3, and the bridge from "the distance to every minimizer converges and weak cluster points are minimizers" to weak convergence of a curve (the continuous Opial lemma).

Difficulty

The obvious approach, a single Lyapunov function that decreases along trajectories, fails here: the Tikhonov term makes the natural energy $\mathcal E_b$ only almost nonincreasing, with an error $\propto t\epsilon(t)\|x^*\|^2$ that must be integrated, and the Hessian term produces a cross term $-\beta\epsilon(t)t^2\langle\nabla g(x),x\rangle$ that is controlled only by condition (a) or (b). The step from $O(1/t^2)$ to $o(1/t^2)$ and to the existence of $\lim\|x(t)-x^*\|$ cannot be read off a single bounded energy: boundedness of $\mathcal E_b$ gives rates in $O$, not in $o$, and says nothing directly about $\|x(t)-x^*\|$. Weak convergence then requires Opial's lemma, which is not in Mathlib in continuous form. Existence of solutions is itself nontrivial in infinite dimension: the right-hand side involves $\nabla^2g$, so the system has to be rewritten as a first-order system before the Cauchy–Lipschitz theorem applies.

Formalization scope

• $\mathcal H$ is any real Hilbert space (InnerProductSpace ℝ H, CompleteSpace H); nothing is finite-dimensional, so weak convergence is genuinely weaker than norm convergence. It is stated through inner products with every $v\in\mathcal H$; stating it as norm convergence, or testing only some $v$, would change the theorem.
• Trajectories are maps $\mathbb R\to\mathcal H$ with velocity and acceleration maps; derivatives are one-sided at $t_0$ and values before $t_0$ play no role. Every theorem is stated for every $C^2$-solution; Theorem 2.1 shows this class is a singleton, so the hypothesis is not vacuous. A sorry-free check that the hypotheses of the goal can be met ($g=0$, $\epsilon=0$, $x=0$) accompanies the drafts.
• $\nabla^2g(x)v$ is the derivative of $\nabla g$ at $x$ applied to $v$; $\dot\epsilon$ is a derivative map tied to $\epsilon$ on $[t_0,+\infty)$; $\min g=\inf g$, attained since $\operatorname{argmin} g\ne\emptyset$; $L^1$ membership is Lebesgue integrability on $[t_0,+\infty)$ (never a comparison of a Bochner integral with a number); $O$, $o$ are Asymptotics relations at $+\infty$.
• Repair. Theorem 2.1 carries the extra hypothesis "$\beta=0$ or $g\in C^2$". As printed (only twice differentiable $g$) it is false for $\beta>0$: on $\mathbb R$, $g'(x)=2x+x^2\sin(1/x)$ gives a strongly convex $g$ with bounded but discontinuous $g''$, and the solution from $u_0=0$, $v_0=1$ is not $C^2$. The paper's own proof uses continuity of $\nabla^2 g$.
• Flagged, not repaired. The last claim of Theorem 3.3, $t^2\|\nabla g(x(t))\|^2\in L^1$, is its own item: the paper's argument for it uses a term proportional to $\beta$, so for $\beta=0$ the paper gives no proof. No counterexample is known, and the statement is drafted as printed.
• Referenced, not re-posed. Lemma A.4 (continuous Opial lemma, p. 30), the last step of the goal's proof, is the published InertialAVD.Traj.lemma_A_2 (same hypotheses: $S$ nonempty, $\lim\|x(t)-z\|$ exists for every $z\in S$, every weak sequential limit point lies in $S$; same conclusion). It is a reference milestone of this mission; its weak convergence InertialAVD.Traj.WeakTendstoAtTop has the same body as TikhonovHDD.Weak.WeakTendstoAtTop.
• The model is defined from scratch: the published HessianDamping.DINAVD.Setting concerns a different equation (no $\epsilon x$ term) and is not reused.
• Welcome contributions: a general Cauchy–Lipschitz existence theorem for locally Lipschitz first-order systems in Banach spaces on $[t_0,+\infty)$ with a blow-up alternative, the continuous Opial lemma, and Lemma A.2 for absolutely continuous functions.

Selected references

• R.I. Boţ, E.R. Csetnek, S.C. László, Tikhonov regularization of a second order dynamical system with Hessian driven damping, arXiv:1911.12845v2, 2020; published in Math. Program., https://doi.org/10.1007/s10107-020-01528-8. https://arxiv.org/abs/1911.12845v2
• H. Attouch, J. Peypouquet, P. Redont, Fast convex optimization via inertial dynamics with Hessian driven damping, J. Differential Equations 261(10), 5734–5783, 2016. https://doi.org/10.1016/j.jde.2016.08.020
• H. Attouch, Z. Chbani, H. Riahi, Combining fast inertial dynamics for convex optimization with Tikhonov regularization, J. Math. Anal. Appl. 457(2), 1065–1094, 2018. https://doi.org/10.1016/j.jmaa.2016.12.017
• H. Attouch, Z. Chbani, J. Peypouquet, P. Redont, Fast convergence of inertial dynamics and algorithms with asymptotic vanishing viscosity, Math. Program. 168, 123–175, 2018. https://doi.org/10.1007/s10107-016-0992-8
• W. Su, S. Boyd, E.J. Candès, A differential equation for modeling Nesterov's accelerated gradient method: theory and insights, J. Mach. Learn. Res. 17(153), 1–43, 2016. https://jmlr.org/papers/v17/15-084.html

Motivation

Thompson sampling (TS) chooses each action with the posterior probability that it is optimal. Its practical success in online advertising, recommendation and experimentation raised the question of why it works and when it does not. The tutorial of Russo, Van Roy, Kazerouni, Osband and Wen (Found. Trends Mach. Learn. 11(1), 2018) answers part of that question in §8.1 with an information-theoretic regret bound due to Russo and Van Roy (JMLR 17, 2016). The bound applies to every algorithm, not only TS. It reduces the analysis of a particular algorithm to one quantity, the information ratio: the squared expected regret of a period divided by the information the period yields about the optimal action.

Timeline:

• 2016: Russo and Van Roy, An information-theoretic analysis of Thompson sampling (JMLR 17), introduce the information ratio with the regret and the mutual information both conditioned on a realised history, and prove the regret bound and bounds on the ratio for TS (for example $|\mathcal X|/2$ for $|\mathcal X|$-armed bandits and $d/2$ for linear bandits).
• 2018: Russo and Van Roy, Learning to optimize via information-directed sampling (Oper. Res. 66), design algorithms that minimise the information ratio directly.
• 2018: Russo and Van Roy, Satisficing in time-sensitive bandit learning (arXiv:1803.02855), observe that the same bounds hold with the information ratio averaged over histories. This is the form (8.7) the tutorial uses (footnote 1, p. 75).
• 2018: the tutorial gives the bound as (8.8), with a three-line proof on p. 76 whose last step is cited from the 2016 paper.

Setting

An agent acts in periods $t=1,2,\dots$. In each period it selects an action $x_t$ from a finite set $\mathcal X$ and observes an outcome $y_t$ from a finite set $\mathcal Y$. An unknown parameter $\theta$, taking values in a measurable space $\Theta$, has a prior distribution $P$. Given $\theta$ and the action $x$, the outcome is drawn from a probability vector $q_\theta(\cdot\mid x)$ that depends measurably on $\theta$. The agent receives the reward $r(y_t)$ for a known function $r$. The mean reward of action $x$ is

$$\mu(x,\theta)=\sum_{y}q_\theta(y\mid x)\,r(y),$$

and the optimal action $x^*=x^*(\theta)$ is a measurable selector of $\arg\max_x\mu(x,\theta)$. Ties may be broken in any measurable way, and every result holds for every selector.

The history $\mathbb H_t$ is the list $(x_1,y_1),\dots,(x_t,y_t)$. An algorithm $\pi$ assigns to every history a probability vector over $\mathcal X$, from which the next action is drawn. The joint law of $(\theta,\mathbb H_t)$ is $P(d\theta)\prod_{s\le t}\pi(x_s\mid\mathbb H_{s-1})\,q_\theta(y_s\mid x_s)$, and every expectation below is taken under it. The expected cumulative regret, or Bayesian regret, is

$$\mathbb E[\mathrm{Regret}(T)]=\mathbb E\Big[\sum_{t=1}^T\big(\mu(x^*,\theta)-\mu(x_t,\theta)\big)\Big].$$

The entropy of the prior distribution of $x^*$ is $H(x^*)=-\sum_a\mathbb P(x^*=a)\log\mathbb P(x^*=a)$. The conditional mutual information $I(x^*;(x_t,y_t)\mid\mathbb H_{t-1})$ is the expected information that the period-$t$ observation carries about $x^*$ given the past, averaged over the past. The information ratio of period $t$ is (8.7):

$$\Gamma_t=\frac{\big(\mathbb E[\mu(x^*,\theta)-\mu(x_t,\theta)]\big)^2}{I(x^*;(x_t,y_t)\mid\mathbb H_{t-1})}.$$

Formalization targets

Goal: (8.8)

For every model and every algorithm, every horizon $T$ and every $\bar\Gamma\ge0$ with $\big(\mathbb E[\mu(x^*,\theta)-\mu(x_t,\theta)]\big)^2\le\bar\Gamma\,I(x^*;(x_t,y_t)\mid\mathbb H_{t-1})$ for $t=1,\dots,T$,

$$\mathbb E[\mathrm{Regret}(T)]\le\sqrt{\bar\Gamma\,H(x^*)\,T}.$$

When every mutual information is positive, the least such $\bar\Gamma$ is $\max_{t\le T}\Gamma_t$, which is the tutorial's $\bar\Gamma$. A companion theorem states that literal form under the positivity hypothesis.

Milestones, in the order of the p. 76 proof

1. Regret is additive over periods:

$$\mathbb E[\mathrm{Regret}(T)]=\sum_{t=1}^T\mathbb E[\mu(x^*,\theta)-\mu(x_t,\theta)].$$

2. The Jensen (Cauchy–Schwarz) step: if $\delta_t^2\le\bar\Gamma I_t$ for every $t$ and $\sum_tI_t\le H$, then $\sum_t\delta_t\le\sqrt{T\bar\Gamma H}$. This milestone is the published, proved theorem BanditAlgorithm.information_ratio_cumulative_bound_signed. It states the step in the form the proof uses, with the entropy budget as a hypothesis.
3. The entropy budget:

$$\sum_{t=1}^T I\big(x^*;(x_t,y_t)\mid\mathbb H_{t-1}\big)\le H(x^*).$$

Significance

The result. (8.8) separates the regret analysis of an algorithm into two independent problems. One is a per-period bound on the information ratio, which depends on the structure of the problem and the algorithm. The other is the entropy of the optimal action, which depends only on the prior. Combined with the bound $\bar\Gamma\le|\mathcal X|/2$ that Russo and Van Roy (2016) prove for TS on $|\mathcal X|$-armed bandits, and with $H(x^*)\le\log|\mathcal X|$, it gives $\sqrt{|\mathcal X|\,T\log|\mathcal X|/2}$. On the shortest-path problem of p. 76, $\bar\Gamma\le d/2$ makes the regret depend on the number $d$ of edges and on $H(x^*)$, not on the exponentially large number of paths. The same inequality is the design criterion of information-directed sampling, which chooses actions to minimise the ratio.

Formalizing it. The bound and its proof are known. The work here is to make the probabilistic model explicit, with a general parameter space, an outcome kernel and randomized algorithms; to formalize the per-period decomposition of the regret and the chain-rule argument behind the entropy budget, which the tutorial cites rather than proves; and to obtain (8.8) for every algorithm in one statement. The Proved platform results for Lattimore and Szepesvári's Bayesian adversarial bandit (Theorem 36.5 and its chain-rule lemmas) are about a different model: one algorithm (TS), the best arm in hindsight of a reward matrix, and the budget $\log k$. To our knowledge, (8.8) itself has not been formalized.

Difficulty

The Cauchy–Schwarz step is elementary and already proved on the platform. Milestone 1 is bookkeeping: the cumulative regret integrates over histories of length $T$, while the period-$t$ term integrates over histories of length $t$, and later periods must be integrated out.

The substance is the entropy budget. The obvious attempt bounds each period's information by $H(x^*)$, which loses a factor $T$. The bound is uniform in $T$ only because the per-period informations telescope: their sum is the information that the whole history carries about $x^*$, which is at most the entropy of $x^*$. Each step requires conditional distributions of $x^*$ given histories, where $\theta$ lives in a general measurable space and is integrated against the prior. Zero-probability histories must also be handled consistently with the convention $0\log0=0$.

Formalization scope

• $\mathcal X$ and $\mathcal Y$ are finite types; $\Theta$ is an arbitrary measurable space with a probability measure $P$. The tutorial allows infinite action sets and general outcomes. The finite specialization makes $H(x^*)$ and every mutual information a finite sum. The kernel $q_\theta(\cdot\mid x)$ of §4 replaces the noise representation $y_t=g(x_t,\theta,w_t)$ of p. 70; both describe the same laws.
• Standing hypotheses: $q_\theta(y\mid x)\ge0$, $\sum_yq_\theta(y\mid x)=1$, measurability of $\theta\mapsto q_\theta(y\mid x)$, measurability of the level sets of $x^*$, and $\mu(x,\theta)\le\mu(x^*(\theta),\theta)$. They make every integrand bounded and measurable, so no integrability hypothesis is added.
• An algorithm is a behavioural policy, giving a probability vector over actions for each history. Under perfect recall this is equivalent to an algorithm with internal randomization. Periods are 0-based in Lean: Lean period $s$ is tutorial period $t=s+1$.
• Logarithms are natural in both $I$ and $H$. The conditional mutual information is defined directly as an average of Kullback–Leibler divergences, with joint probabilities. It is not defined as a difference of unconditional mutual informations: that would make the chain rule true by definition and empty milestone 3. $\Gamma_t$ is never written as a quotient in the goal, so a period with zero information cannot turn the bound into one with $\bar\Gamma=0$.
• Needed infrastructure: marginalisation over later periods of the history law, the chain rule of conditional mutual information for finite variables averaged against a general prior, and $I(x^*;\mathbb H_T)\le H(x^*)$. The Proved platform results InformationTheory.kernel_average_categorical_kl_le_marginal_entropy and its _general variant (average KL to the mixture is at most the mixture's entropy), and the adversarial-bandit lemmas BanditAlgorithm.bayesian_ts_conditional_information_gain_chain_rule and bayesian_history_mutualInformation_endpoint_bound, are related tools. A chain-rule library for this history model would be reusable beyond this mission, for example for information-directed sampling. Proofs of auxiliary lemmas are welcome.

The tutorial's argument is on pp. 75–76. The tutorial cites Russo and Van Roy (2016) for the original results, including the entropy budget.

Selected references

• D. J. Russo, B. Van Roy, A. Kazerouni, I. Osband, Z. Wen, A Tutorial on Thompson Sampling, Foundations and Trends in Machine Learning 11(1):1–96, 2018. https://doi.org/10.1561/2200000070
• D. Russo, B. Van Roy, An Information-Theoretic Analysis of Thompson Sampling, Journal of Machine Learning Research 17(68):1–30, 2016. https://jmlr.org/papers/v17/14-087.html (arXiv: https://arxiv.org/abs/1403.5341)
• D. Russo, B. Van Roy, Learning to Optimize via Information-Directed Sampling, Operations Research 66(1):230–252, 2018. https://doi.org/10.1287/opre.2017.1663
• D. Russo, B. Van Roy, Satisficing in Time-Sensitive Bandit Learning, arXiv preprint, 2018. https://arxiv.org/abs/1803.02855
• T. Lattimore, C. Szepesvári, Bandit Algorithms, Cambridge University Press, 2020, Chapter 36. https://doi.org/10.1017/9781108571401

Why Sinkhorn's algorithm needs a convergence rate

Sinkhorn's algorithm is the standard method for computing entropically regularized optimal transport between two histograms. It is used throughout machine learning, statistics, imaging and graphics because each iteration is two matrix–vector products followed by entrywise divisions, which parallelize on GPUs and differentiate easily. The same iteration, under the names matrix scaling, RAS or iterative proportional fitting, has been used since the 1930s in transportation planning, economics and contingency-table analysis. Whether the iteration converges, and how fast, determines how many iterations a practitioner must run and which stopping criteria are reliable.

This mission formalizes the global convergence analysis of Sinkhorn's algorithm as presented in Chapter 4 (§4.2, Remarks 4.12–4.14) of G. Peyré and M. Cuturi, Computational Optimal Transport (2019).

Timeline. Sinkhorn (1964) proved that a positive matrix can be scaled to be doubly stochastic, and Sinkhorn and Knopp (1967) extended this to prescribed marginals. Birkhoff (1957) and Samelson (1957) independently showed that a positive linear map is a strict contraction of the positive cone for Hilbert's projective metric, with an explicit ratio. Franklin and Lorenz (1989) combined the two to prove that Sinkhorn's scalings converge linearly in Hilbert's metric at rate $\lambda(K)^2$. Linial, Samorodnitsky and Wigderson (1998) related the iteration to the permanent.

Setting

Fix $n, m \ge 1$, a matrix $K \in \mathbb{R}^{n\times m}$ with positive entries (in the book, the Gibbs kernel $K_{i,j} = e^{-C_{i,j}/\varepsilon}$), and vectors $a \in \mathbb{R}^n$, $b \in \mathbb{R}^m$ with positive entries. Products $\odot$ and quotients of vectors are entrywise. The regularized transport plan has the form $P = \mathrm{diag}(u)K\mathrm{diag}(v)$, where the scalings $u, v$ solve

$$u \odot (Kv) = a, \qquad v \odot (K^\top u) = b. \tag{4.14}$$

Sinkhorn's algorithm alternately enforces each equation:

$$v^{(0)} = \mathbb{1}_m,\qquad u^{(\ell+1)} = \frac{a}{Kv^{(\ell)}},\qquad v^{(\ell+1)} = \frac{b}{K^\top u^{(\ell+1)}}. \tag{4.15}$$

For vectors with positive entries, Hilbert's projective metric is

$$d_{\mathcal H}(u,u') = \log\max_{i,j}\frac{u_i u'_j}{u_j u'_i},$$

which vanishes exactly when $u$ and $u'$ are proportional, and equals the variation seminorm $\|\log u - \log u'\|_{\mathrm{var}}$ with $\|f\|_{\mathrm{var}} = \max_i f_i - \min_i f_i$. Birkhoff's constants of $K$ are

$$\eta(K) = \max_{i,j,k,\ell}\frac{K_{i,k}K_{j,\ell}}{K_{j,k}K_{i,\ell}},\qquad \lambda(K) = \frac{\sqrt{\eta(K)}-1}{\sqrt{\eta(K)}+1}.$$

Write $P^{(\ell)} = \mathrm{diag}(u^{(\ell)})K\mathrm{diag}(v^{(\ell)})$ and $P^\star = \mathrm{diag}(u^\star)K\mathrm{diag}(v^\star)$ for a positive solution $(u^\star, v^\star)$ of (4.14).

Formalization targets

Goal: Theorem 4.2 (Franklin–Lorenz)

For every $\ell$,

$$d_{\mathcal H}(v^{(\ell)}, v^\star) \le \lambda(K)^{2\ell}\, d_{\mathcal H}(v^{(0)}, v^\star),\qquad d_{\mathcal H}(u^{(\ell+1)}, u^\star) \le \lambda(K)^{2\ell+1}\, d_{\mathcal H}(v^{(0)}, v^\star),$$

so both distances tend to $0$; for $\ell \ge 1$,

$$d_{\mathcal H}(u^{(\ell)}, u^\star) \le \frac{d_{\mathcal H}(P^{(\ell)}\mathbb{1}_m, a)}{1-\lambda(K)^2},\qquad \max_{i,j}\big|\log P^{(\ell)}_{i,j} - \log P^\star_{i,j}\big| \le d_{\mathcal H}(u^{(\ell)},u^\star) + d_{\mathcal H}(v^{(\ell)},v^\star),$$

and for $\ell \ge 0$ the mirrored bound for $v^{(\ell)}$ with the half-step coupling $\mathrm{diag}(u^{(\ell+1)})K\mathrm{diag}(v^{(\ell)})$. The book writes the rate as $O(\lambda(K)^{2\ell})$; the constants stated are those its proof yields.

Milestones, in attack order

1. (4.21): $d_{\mathcal H}(u,u') = \|\log u - \log u'\|_{\mathrm{var}}$.
2. Remark 4.12: $d_{\mathcal H}$ is a distance on rays (triangle inequality; zero exactly on proportional vectors).
3. Invariances: $d_{\mathcal H}(v,v') = d_{\mathcal H}(v/v',\mathbb{1}) = d_{\mathcal H}(\mathbb{1}/v,\mathbb{1}/v')$ and invariance under positive diagonal scaling.
4. Theorem 4.1 (Birkhoff–Hopf): $d_{\mathcal H}(Kv,Kv') \le \lambda(K)\,d_{\mathcal H}(v,v')$ and $\lambda(K) < 1$.
5. One half-step contracts by $\lambda(K)$: $d_{\mathcal H}(u^{(\ell+1)},u^\star) \le \lambda(K)\,d_{\mathcal H}(v^{(\ell)},v^\star)$, and the $v$-analogue.
6. The display behind (4.23): $d_{\mathcal H}(u^{(\ell)},u^\star) \le d_{\mathcal H}(a, u^{(\ell)}\odot Kv^{(\ell)}) + \lambda(K)^2 d_{\mathcal H}(u^{(\ell)},u^\star)$, and its mirror.

Significance

The result. Theorem 4.2 gives a global linear rate that holds from the first iteration, independent of the starting point, with a ratio computable from $K$ alone. Through (4.21) the rate transfers to the dual potentials $\varepsilon\log u^{(\ell)}$, and the bounds (4.23) turn the marginal violations $P^{(\ell)}\mathbb{1}_m$ versus $a$ — which are computed anyway — into certified stopping criteria. Theorem 4.1 alone also yields a quantitative Perron–Frobenius theorem (Remark 4.13).

Formalizing it. All results here are classical and proved in the literature (Birkhoff 1957; Franklin–Lorenz 1989); Theorem 4.1 is quoted by the book without proof. To our knowledge neither Hilbert's projective metric on the positive cone, the Birkhoff–Hopf theorem, nor the convergence rate of Sinkhorn's algorithm has a machine-checked proof in Mathlib. The mission produces a reusable development of $d_{\mathcal H}$ and Birkhoff's contraction theorem, and a certified rate and stopping rule for the most widely used optimal-transport solver.

Difficulty

Milestones 1–3 and 5–6 are finite manipulations of maxima of ratios. The central difficulty is Theorem 4.1. The obvious argument — bound $(Kv)_i/(Kv')_i$ entrywise by $\max_k v_k/v'_k$ — only shows that $K$ is non-expansive ($\lambda \le 1$); it loses the strict contraction entirely. The sharp ratio $(\sqrt{\eta}-1)/(\sqrt{\eta}+1)$ depends on how $K$ acts on the two vectors jointly, which the entrywise argument does not see. The book quotes the theorem without proof, so the argument has to be supplied from Birkhoff's original paper or later expositions.

Formalization scope

Lean namespace CompOT.Sinkhorn. Indices are $0$-based (Fin n, Fin m); vectors are functions Fin n → ℝ, matrices Matrix (Fin n) (Fin m) ℝ, and $\odot$ and $/$ are the pointwise operations on functions. $d_{\mathcal H}$, $\eta$ and $\|\cdot\|_{\mathrm{var}}$ use real iSup/iInf over the finite index sets. Every theorem assumes $0 < n$, $0 < m$ and positive entries for $K$, $a$, $b$ and every vector passed to $d_{\mathcal H}$, which excludes Lean's junk values for division by zero and logarithms of non-positive numbers. Sinkhorn's algorithm is a run predicate (IsSinkhornRun) stating the recursion (4.15) with $v^{(0)} = \mathbb{1}_m$; $u^{(0)}$ is not an iterate and is unconstrained. The solution $(u^\star, v^\star)$ is any positive solution of (4.14), taken as a hypothesis (it exists iff $\sum_i a_i = \sum_j b_j$); the book's histograms are probability vectors, but only positivity is used. $P^\star$ in (4.24) is $\mathrm{diag}(u^\star)K\mathrm{diag}(v^\star)$, the solution of the regularized problem (4.2) by Proposition 4.3.

Deviations from the print, all disclosed in the items: the $O(\cdot)$ of (4.22) is replaced by the proof's explicit constants; the limit $(u^{(\ell)}, v^{(\ell)}) \to (u^\star, v^\star)$ is stated as $d_{\mathcal H} \to 0$, since scalings are determined only up to $(ru, v/r)$; the second inequality of (4.23) is false as printed (its right side is $0$ for $\ell \ge 1$) and is stated with the half-step coupling the mirrored proof gives.

A definition of $d_{\mathcal H}$ without the logarithm, or as a maximum of single-entry ratios, would make the statements different theorems; the definition here is the book's, and (4.21) is a milestone checking it against the variation seminorm. A run predicate not tied to $K$, $a$, $b$, or an unsatisfiable one, is ruled out: the predicate is the recursion itself and is inhabited for every $K$, $a$, $b$.

Contributions welcome: the metric properties of $d_{\mathcal H}$ (reusable for Perron–Frobenius and nonlinear consensus), a proof of the Birkhoff–Hopf theorem, and positivity lemmas for Sinkhorn iterates.

Selected references

• G. Peyré, M. Cuturi, Computational Optimal Transport, Foundations and Trends in Machine Learning 11(5–6):355–607, 2019. https://doi.org/10.1561/2200000073 (§4.2, pp. 431–443)
• J. Franklin, J. Lorenz, On the scaling of multidimensional matrices, Linear Algebra and its Applications 114/115:717–735, 1989. https://doi.org/10.1016/0024-3795(89)90490-4
• G. Birkhoff, Extensions of Jentzsch's theorem, Transactions of the AMS 85(1):219–227, 1957. https://doi.org/10.1090/S0002-9947-1957-0087058-6
• H. Samelson, On the Perron–Frobenius theorem, Michigan Mathematical Journal 4(1):57–59, 1957. https://doi.org/10.1307/mmj/1028990177
• R. Sinkhorn, A relationship between arbitrary positive matrices and doubly stochastic matrices, Annals of Mathematical Statistics 35(2):876–879, 1964. https://doi.org/10.1214/aoms/1177703591
• N. Linial, A. Samorodnitsky, A. Wigderson, A deterministic strongly polynomial algorithm for matrix scaling and approximate permanents, STOC 1998. https://doi.org/10.1145/276698.276880

Motivation

The Kalman filter estimates the hidden state of a linear system from noisy observations. At every step it solves a minimum mean square error (MMSE) estimation problem: given a joint normal distribution of a signal $x$ and an observation $y$, find the estimator $\psi(y)$ minimizing $\mathbb E[\|x-\psi(y)\|^2]$. The answer, the conditional mean, is optimal only if the model distribution is correct, and in practice the system matrices and noise covariances are estimated and misspecified.

Shafieezadeh-Abadeh, Nguyen, Kuhn and Mohajerin Esfahani (arXiv:1809.08830, NeurIPS 2018) replace the single model by a set of plausible models: all normal distributions within type-2 Wasserstein distance $\rho$ of a nominal normal distribution. The estimator is then chosen to minimize the worst-case mean square error over this set. Their Theorem 2.5 shows that this robust estimation problem, although posed over an infinite-dimensional space of estimators and a nonconvex set of distributions, has the same value as a finite convex program over covariance matrices, and that the robust estimator is affine. Applied at every time step, this yields the Wasserstein distributionally robust Kalman filter of the paper's §4.

This mission formalizes Theorem 2.5 together with the minimax theorem (Theorem 2.3) and the lemmas of its proof.

Setting

Let $n, m \ge 0$ and $d = n+m$. A point of $\mathbb R^d$ is $z = [x; y]$ with signal $x\in\mathbb R^n$ and observation $y\in\mathbb R^m$. A symmetric matrix $S \in \mathbb R^{d\times d}$ is partitioned as $S = \begin{bmatrix} S_{xx} & S_{xy}\\ S_{yx} & S_{yy}\end{bmatrix}$. $\mathbb S^d_+$ ($\mathbb S^d_{++}$) denotes the positive semidefinite (definite) matrices, $A \succeq B$ means $A-B\in\mathbb S^d_+$, $\langle A,B\rangle = \operatorname{Tr}[A^\top B]$, and $A^{1/2}$ is the positive semidefinite square root of $A\succeq 0$.

Normal distributions. $\mathcal N_d(c,S)$ is the normal distribution on $\mathbb R^d$ with mean $c$ and covariance $S\in\mathbb S^d_+$; singular $S$ is allowed. $\mathcal N_d$ is the set of all of them.

Wasserstein distance. For distributions $\mathbb Q_1,\mathbb Q_2$ on $\mathbb R^d$,

$$W_2(\mathbb Q_1,\mathbb Q_2) = \inf_{\pi} \Bigl(\int \|z_1-z_2\|^2\,\pi(dz_1,dz_2)\Bigr)^{1/2},$$

the infimum ranging over all distributions $\pi$ on $\mathbb R^d\times\mathbb R^d$ with marginals $\mathbb Q_1$ and $\mathbb Q_2$.

Ambiguity set. Fix $\mu = [\mu_x;\mu_y]\in\mathbb R^d$, $\Sigma\in\mathbb S^d_{++}$, the nominal distribution $\mathbb P=\mathcal N_d(\mu,\Sigma)$ and a radius $\rho\ge 0$. The Wasserstein ambiguity set is

$$\mathcal P = \{\mathbb Q\in\mathcal N_d : W_2(\mathbb Q,\mathbb P)\le\rho\}.$$

Robust estimation. Let $\mathcal L$ be the set of all measurable functions $\psi:\mathbb R^m\to\mathbb R^n$. The distributionally robust MMSE problem is

$$\inf_{\psi\in\mathcal L}\ \sup_{\mathbb Q\in\mathcal P}\ \mathbb E^{\mathbb Q}\bigl[\|x-\psi(y)\|^2\bigr], \tag{2}$$

and a distributionally robust MMSE estimator is a $\psi$ attaining the outer infimum. Nature's problem exchanges the two optimizations; a $\mathbb Q^\star\in\mathcal P$ attaining its outer supremum is a least favorable prior.

The finite program. With $\underline\sigma = \lambda_{\min}(\Sigma)$, program (5) maximizes $f(S)=\operatorname{Tr}[S_{xx}-S_{xy}S_{yy}^{-1}S_{yx}]$ over $S\in\mathbb S^d_+$ (with $S_{xx}\in\mathbb S^n_+$, $S_{yy}\in\mathbb S^m_+$) subject to

$$\operatorname{Tr}\bigl[S+\Sigma-2(\Sigma^{1/2}S\Sigma^{1/2})^{1/2}\bigr]\le\rho^2,\qquad S\succeq\underline\sigma I_d .$$

Formalization targets

Goal: Theorem 2.5 (p. 3)

1. The optimal value of (2) equals the optimal value of (5).
2. If $S^\star$ is optimal in (5), then

$$\psi^\star(y) = S^\star_{xy}(S^\star_{yy})^{-1}(y-\mu_y)+\mu_x$$

is a distributionally robust MMSE estimator. 3. For the same $S^\star$, $\mathbb Q^\star=\mathcal N_d(\mu,S^\star)$ lies in $\mathcal P$ and is a least favorable prior.

Existence and uniqueness of an optimal $S^\star$ are not part of the statement.

Milestones, in the order of the proof

• Proposition 2.2 (Gelbrich formula): for $\Sigma_1,\Sigma_2\succeq 0$,

$$W_2\bigl(\mathcal N_d(\mu_1,\Sigma_1),\mathcal N_d(\mu_2,\Sigma_2)\bigr)=\sqrt{\|\mu_1-\mu_2\|^2+\operatorname{Tr}\bigl[\Sigma_1+\Sigma_2-2(\Sigma_2^{1/2}\Sigma_1\Sigma_2^{1/2})^{1/2}\bigr]}.$$

• Lemma A.1: the closed form of $\sup_{S\succeq0}\langle D,S\rangle-\gamma\operatorname{Tr}[S-2(\Sigma^{1/2}S\Sigma^{1/2})^{1/2}]$ and its unique maximizer.
• (A.1b): under a fixed normal distribution, affine estimators $y\mapsto Gy+g$ achieve the minimum mean square error over $\mathcal L$.
• Theorem 2.3 (minimax theorem): $\inf_{\psi\in\mathcal L}\sup_{\mathbb Q\in\mathcal P}\mathbb E^{\mathbb Q}[\|x-\psi(y)\|^2] = \sup_{\mathbb Q\in\mathcal P}\inf_{\psi\in\mathcal L}\mathbb E^{\mathbb Q}[\|x-\psi(y)\|^2]$.
• Lemma A.2: the maximizer of $\langle S,D\rangle$ under the Wasserstein-type trace constraint, and $S^\star\succeq\underline\sigma I_d$.
• (A.6): for $S\succ0$, $\min_G\bigl\langle\begin{bmatrix}I_n&-G\\-G^\top&G^\top G\end{bmatrix},S\bigr\rangle=f(S)$, uniquely at $G^\star=S_{xy}S_{yy}^{-1}$.

Significance

The result. Theorem 2.5 gives the robust estimator in closed form from one optimal covariance matrix, and identifies the least favorable distribution as a normal distribution with the nominal mean. Together with Theorem 2.3 it shows that the pair $(\psi^\star,\mathbb Q^\star)$ is a saddle point: the robust estimator is the Bayesian estimator for the least favorable prior. For $\rho=0$ the program has the single feasible point $\Sigma$ and the estimator reduces to the classical conditional mean, so the theorem contains the Kalman update as a special case. The finite program is what makes the robust filter of §4 computable; the paper's §3 (formalized in the companion mission of this series) solves it with a Frank–Wolfe method.

Formalizing it. The result is proved in the paper; nothing here is open. No part of it is machine-checked yet. The Gelbrich formula for normal distributions (Proposition 2.2) is cited from the literature and has no formal proof in Mathlib; neither does Lemma A.1, cited from earlier work of the same authors. A closely related statement for a Wasserstein ball over all distributions around an elliptical nominal (Kuhn et al. 2019, Theorem 25) is posed on the platform as WassersteinDRO.Shrinkage.distributionally_robust_mmse_estimator; it is a different theorem, since the ambiguity set there is convex and contains non-normal distributions.

Difficulty

The ambiguity set $\mathcal P$ is not convex, because a mixture of normal distributions is generally not normal. The standard route to a minimax equality, Sion's theorem, needs convexity on nature's side and therefore does not apply to (2) directly; the obvious attempt to interchange the infimum and the supremum fails at this step. The space $\mathcal L$ of measurable estimators is infinite-dimensional, and the Wasserstein constraint is an infimum over couplings, so a proof must pass to finite-dimensional descriptions on both sides (means and covariances for $\mathbb Q$, affine maps for $\psi$) and justify each passage. The matrix-analytic parts involve the nonlinear map $S\mapsto(\Sigma^{1/2}S\Sigma^{1/2})^{1/2}$, for which Mathlib has square roots but little calculus.

Formalization scope

• $\mathbb R^d$ is EuclideanSpace ℝ (Fin n ⊕ Fin m); the Sum.inl coordinates are $x$, the Sum.inr coordinates are $y$, and the blocks of a matrix are toBlocks₁₁, toBlocks₁₂, toBlocks₂₁, toBlocks₂₂.
• $\mathcal N_d(c,S)$ is Mathlib's multivariateGaussian c S with S.PosSemidef; only the nominal $\Sigma$ is required to be positive definite.
• $W_2$ is the published coupling-based WassersteinDRO.Shrinkage.wassersteinDistance 2 in $[0,\infty]$. The ambiguity set is defined through it, not through the Gelbrich formula; the formula is Proposition 2.2, a milestone.
• Expectations are lower Lebesgue integrals in $[0,\infty]$, with no integrability guard. A non-integrable estimator therefore has infinite risk and cannot lower the infimum through a junk value.
• The values of (2) and of nature's problem are in $[0,\infty]$; the value of (5) is the published sdpValue in EReal, and the two are compared through the coercion.
• $\underline\sigma$ is passed as a real number that is the least eigenvalue of $\Sigma$.
• "Equivalent" in Theorem 2.5 is read as the three claims above. A version that asserts only equality of values, or optimality of $\psi^\star$ among affine estimators only, is weaker than the theorem and is not the goal.
• Lemma A.2 carries the added hypothesis $\rho>0$, without which its algebraic equation has no admissible root.

A complete development needs the Gelbrich formula for (possibly degenerate) normal distributions, conditional means of jointly normal vectors, a minimax theorem for convex–concave functions on compact sets, Lagrangian duality for a convex semidefinite program, and calculus of the positive semidefinite square root. The first three are reusable well beyond this mission, and contributions of any of them as separate results are welcome.

Selected references

• S. Shafieezadeh-Abadeh, V. A. Nguyen, D. Kuhn, P. Mohajerin Esfahani, Wasserstein Distributionally Robust Kalman Filtering, NeurIPS 2018; arXiv:1809.08830v3. https://arxiv.org/abs/1809.08830
• M. Gelbrich, On a formula for the L² Wasserstein metric between measures on Euclidean and Hilbert spaces, Mathematische Nachrichten 147, 1990. https://doi.org/10.1002/mana.19901470121
• V. A. Nguyen, D. Kuhn, P. Mohajerin Esfahani, Distributionally robust inverse covariance estimation, arXiv:1802.01981, 2018. https://arxiv.org/abs/1802.01981
• D. Kuhn, P. Mohajerin Esfahani, V. A. Nguyen, S. Shafieezadeh-Abadeh, Wasserstein distributionally robust optimization: theory and applications in machine learning, INFORMS TutORials, 2019; arXiv:1908.08729. https://arxiv.org/abs/1908.08729
• D. P. Bertsekas, Convex Optimization Theory, Athena Scientific, 2009.

Motivation

A Kalman filter estimates the hidden state $x$ of a linear system from observations $y$, assuming the joint law of $z = (x, y)$ is a known Gaussian. When that law is only approximately known, Shafieezadeh-Abadeh, Nguyen, Kuhn and Mohajerin Esfahani (arXiv:1809.08830v3, NeurIPS 2018) replace it by the worst Gaussian within a 2-Wasserstein ball of radius $\rho$ around a nominal Gaussian with covariance $\Sigma$. Their Theorem 2.5 shows that the resulting minimax estimation problem is equivalent to a finite nonlinear semidefinite program, program (5), whose objective is $f(S) = \mathrm{Tr}[S_{xx} - S_{xy}S_{yy}^{-1}S_{yx}]$.

Each step of the robust filter solves program (5), so the filter is only practical if (5) can be solved quickly. The paper solves it with a Frank–Wolfe method (Algorithm 2), which avoids projections by replacing the objective with its linearization at the current iterate and maximizing that linear function over the original feasible set. That inner problem, the direction-finding subproblem, is again a nonlinear SDP. The paper's Algorithm 1 solves it by a one-dimensional bisection, and Theorem 3.2 states that this bisection returns a feasible solution whose objective value is within a prescribed tolerance $\varepsilon$ of the optimum. This mission formalizes that guarantee.

Setting

Fix dimensions $n \ge 1$ (state) and $m$ (observation), $d = n + m$, and index coordinates so that the first $n$ belong to $x$. Let $\mathbb{S}^d_+$ and $\mathbb{S}^d_{++}$ denote the positive semidefinite and positive definite $d\times d$ matrices, and $\langle A, B\rangle = \mathrm{Tr}[A^\top B]$ the trace inner product. The data are a nominal covariance $\Sigma \in \mathbb{S}^d_{++}$ with smallest eigenvalue $\underline\sigma = \lambda_{\min}(\Sigma)$, a radius $\rho > 0$, a tolerance $\varepsilon > 0$, and the current Frank–Wolfe iterate $S \in \mathbb{S}^d_+$, written in blocks $S_{xx}, S_{xy}, S_{yx}, S_{yy}$.

The gradient matrix is

$$D = \nabla f(S) = [I_n, -G]^\top [I_n, -G], \qquad G = S_{xy}S_{yy}^{-1},$$

which is positive semidefinite and nonzero. The direction-finding subproblem (7b) is

$$\max_{L \succeq \underline\sigma I_d}\ \langle L, D\rangle \quad \text{s.t.}\quad \mathrm{Tr}\Big[L + \Sigma - 2\big(\Sigma^{1/2} L \Sigma^{1/2}\big)^{1/2}\Big] \le \rho^2 ;$$

its constraint says that the Gaussian with covariance $L$ lies within Wasserstein distance $\rho$ of the nominal Gaussian, and $\mathcal F$ denotes its feasible set. For $\gamma$ with $\gamma I_d \succ D$ the paper defines

$$h(\gamma) = \rho^2 - \big\langle \Sigma, (I_d - \gamma(\gamma I_d - D)^{-1})^2\big\rangle, \qquad L(\gamma) = \gamma^2 (\gamma I_d - D)^{-1}\Sigma(\gamma I_d - D)^{-1},$$ $$\Delta(\gamma) = \gamma(\rho^2 - \mathrm{Tr}[\Sigma]) - \langle L(\gamma), D\rangle + \gamma^2\langle(\gamma I_d - D)^{-1}, \Sigma\rangle .$$

Algorithm 1 takes the largest eigenvalue $\lambda_1$ of $D$ and a unit eigenvector $v_1$, and starts from the bracket $[LB, UB] = [\gamma_{\min}, \gamma_{\max}]$ with

$$\gamma_{\min} = \lambda_1\big(1 + \sqrt{v_1^\top\Sigma v_1}/\rho\big), \qquad \gamma_{\max} = \lambda_1\big(1 + \sqrt{\mathrm{Tr}[\Sigma]}/\rho\big).$$

Each pass sets $\gamma = (UB + LB)/2$ and $L = L(\gamma)$, replaces $LB$ by $\gamma$ if $h(\gamma) < 0$ and $UB$ by $\gamma$ otherwise, and the loop exits with output $L$ as soon as $h(\gamma) > 0$ and $\Delta(\gamma) < \varepsilon$.

Formalization targets

Goal: Theorem 3.2

For all inputs as above:

1. Soundness. At every pass $k$ whose exit test holds, with trial point $\gamma_k$,

$$L(\gamma_k) \in \mathcal F \quad\text{and}\quad \langle L', D\rangle \le \langle L(\gamma_k), D\rangle + \varepsilon \ \text{ for all } L' \in \mathcal F .$$

2. Termination. If $h$ has no zero at a dyadic point $\gamma_{\min} + j(\gamma_{\max} - \gamma_{\min})/2^k$, $0 < j \le 2^k$, then some pass satisfies the exit test.

The termination hypothesis is a correction: the printed theorem asserts that the algorithm outputs a solution, but if a trial point (or $\gamma_{\max}$) coincides exactly with the root of $h$, the upper end of the bracket becomes that root, all later trial points have $h < 0$, and the loop never exits. Soundness holds unconditionally. The normalization $\|v_1\| = 1$ is also added; the page says only "an eigenvector", and $\gamma_{\min}$ is not a lower bound for a rescaled one.

Milestones (App. A.3, p. 13, in the order the proof uses them)

• $L(\gamma) \in \mathcal F$ whenever $\gamma I_d \succ D$ and $h(\gamma) > 0$.
• Weak duality: $\langle L', D\rangle \le \gamma(\rho^2 - \mathrm{Tr}[\Sigma]) + \gamma^2\langle (\gamma I_d - D)^{-1}, \Sigma\rangle$ for all $L' \in \mathcal F$ and $\gamma \ge 0$ with $\gamma I_d \succ D$.
• The suboptimality display: if $h(\gamma^\star) = 0$ then $\langle L(\gamma^\star) - L(\gamma), D\rangle \le \gamma(\rho^2 - \mathrm{Tr}[\Sigma]) + \gamma^2\langle(\gamma I_d - D)^{-1}, \Sigma\rangle - \langle L(\gamma), D\rangle$.
• Lemma A.3: every root $\gamma^\star$ of $h$ with $\gamma^\star I_d \succ D$ lies in $[\gamma_{\min}, \gamma_{\max}]$.

Significance

Theorem 3.2 is what makes Algorithm 1 a valid linear-maximization oracle for the Frank–Wolfe method on program (5): the convergence analysis of Frank–Wolfe with inexact oracles requires exactly a feasible point whose linearized objective is within a known tolerance of the optimum. It also shows that the inner problem, a semidefinite program with a matrix square root in the constraint, reduces to a scalar root-finding problem once $D$ is diagonalized, which is why the robust filter runs at a cost comparable to an eigendecomposition per step.

The result is proved in the paper; no machine-checked version exists. The formal version adds what the paper leaves implicit: the termination edge case, the normalization of $v_1$, and the requirement $D \ne 0$ (here $n \ge 1$). The milestone statements on Lagrangian weak duality for the Wasserstein–Gaussian trust region and on the location of the multiplier are reusable for other covariance-robust programs of the same shape, such as Wasserstein shrinkage estimation.

Difficulty

The algorithm is elementary; the analysis is not. Feasibility of $L(\gamma)$ requires computing the matrix square root $(\Sigma^{1/2}L(\gamma)\Sigma^{1/2})^{1/2}$ in closed form and comparing eigenvalues of non-commuting products to obtain $L(\gamma) \succeq \underline\sigma I_d$. The dual bound is a semidefinite Lagrangian duality statement over a constraint involving a matrix square root, whose inner supremum (Lemma A.1, cited from earlier work) must be evaluated in closed form. Termination needs monotonicity and continuity of $h$ on $(\lambda_1, \infty)$ and a careful account of which side of the root each bisection update lands on; the obvious invariant "the root stays strictly inside the bracket" fails precisely in the excluded dyadic case.

Formalization scope

Matrices are indexed by Fin n ⊕ Fin m with blocks toBlocks₁₁ ($xx$), toBlocks₁₂ ($xy$), toBlocks₂₂ ($yy$). The feasible set $\mathcal F$ is the published definition WassersteinDRO.Shrinkage.sdpFeasibleSet ρ Σ σ, which uses the published positive semidefinite square root psdSqrt. Mathlib's matrix inverse returns $0$ on singular matrices; every statement involving $h$, $L$ or $\Delta$ assumes $\gamma I_d \succ D$ or works at trial points that lie above $\lambda_1$, and $S_{yy}^{-1}$ is allowed to be this junk value because the theorem admits any $S \in \mathbb{S}^d_+$ (then $D = \mathrm{diag}(I_n, 0)$). The constants $\underline\sigma$ and $\lambda_1$ are passed as reals together with hypotheses that pin them to the smallest eigenvalue of $\Sigma$ and the largest eigenvalue of $D$. The algorithm is modelled as its sequence of passes, not as a function returning an output, so a non-terminating run cannot silently produce a value.

A formalization that defines the output with a fuel bound or a default value, that states only optimality of $L(\gamma^\star)$ at the exact root, or that assumes termination as a hypothesis would not be Theorem 3.2; $\varepsilon$-suboptimality is stated against every feasible point, not against a supremum that could be a junk value.

A complete development needs the spectral theorem for symmetric matrices, uniqueness of the positive semidefinite square root, and semidefinite Lagrangian duality; Mathlib supplies the first two. Proofs of the milestones are welcome contributions.

Selected references

• S. Shafieezadeh-Abadeh, V. A. Nguyen, D. Kuhn, P. Mohajerin Esfahani, Wasserstein Distributionally Robust Kalman Filtering, NeurIPS 2018; arXiv:1809.08830v3. https://arxiv.org/abs/1809.08830
• V. A. Nguyen, D. Kuhn, P. Mohajerin Esfahani, Distributionally Robust Inverse Covariance Estimation: The Wasserstein Shrinkage Estimator, Optimization Online, 2018; arXiv:1805.07194 (reference [18] of the paper; source of Lemma A.1). https://arxiv.org/abs/1805.07194
• M. Jaggi, Revisiting Frank–Wolfe: Projection-Free Sparse Convex Optimization, ICML 2013 (reference [13]; Frank–Wolfe with approximate oracles). https://proceedings.mlr.press/v28/jaggi13.html