All missions

Motivation

Markov chain Monte Carlo (MCMC) estimates an expectation $\int f\,d\pi$ by the average of $f$ along a Markov chain whose invariant distribution is $\pi$. Many chains share the same $\pi$: every Metropolis–Hastings acceptance rule that satisfies detailed balance, every mixture of such kernels, every choice of proposal. Practitioners need a criterion for preferring one of them. The standard yardstick is the asymptotic variance of the ergodic average, the constant in the Markov chain central limit theorem. A smaller asymptotic variance means fewer iterations for the same Monte Carlo error.

Peskun (1973) compared chains on a finite state space through a partial order on transition matrices: if one reversible matrix moves off the diagonal at least as much as another, entry by entry, its asymptotic variances are no larger, for every function. That result justifies the Metropolis–Hastings acceptance probability as the best possible among reversible acceptance rules. It covers only finite state spaces, while MCMC is used almost exclusively on continuous or mixed ones.

Tierney (1998) extended Peskun's theorem to general state spaces, using the spectral approach of Kipnis and Varadhan (1986) for reversible chains. The theorem is the one usually cited when an MCMC paper argues that one sampler dominates another; Mira (2001) surveys orderings built on it.

Setting

Let $(E, \mathcal E)$ be a measurable space in which singletons are measurable, and let $\pi$ be a probability measure on $E$. A Markov kernel $H$ assigns to each $x \in E$ a probability measure $H(x, \cdot)$, measurably in $x$. It acts on functions by $(Hf)(x) = \int f(y)\,H(x,dy)$. The measure $\pi$ is invariant for $H$ if $\int H(x, A)\,\pi(dx) = \pi(A)$ for every $A \in \mathcal E$. The kernel $H$ is reversible with respect to $\pi$ (satisfies detailed balance) if

$$\pi(dx)\,H(x,dy) = \pi(dy)\,H(y,dx),$$

that is, $\int_A H(x,B)\,\pi(dx) = \int_B H(x,A)\,\pi(dx)$ for all $A, B \in \mathcal E$. Reversibility implies invariance.

Write $\langle f, g\rangle = \int fg\,d\pi$, $L^2(\pi)$ for the square-integrable functions and $L^2_0(\pi) = \{g \in L^2(\pi) : \int g\,d\pi = 0\}$.

Off-diagonal domination. For kernels $P_1, P_2$, $P_1 \succeq P_2$ (OffDiagDominates π P₁ P₂) if for $\pi$-almost every $x$,

$$P_1(x, A\setminus\{x\}) \ge P_2(x, A \setminus\{x\}) \quad\text{for all } A \in \mathcal E.$$

So from almost every state $P_1$ moves to every region at least as readily as $P_2$, and the kernels differ only in the probability of staying put.

The chain and its asymptotic variance. For a Markov kernel $H$, let $X_0, X_1, \dots$ be the Markov chain with initial distribution $\pi$ and transition kernel $H$ (chainMeasure π H, a measure on paths $\mathbb N \to E$). For $f \in L^2_0(\pi)$ put $S_n = \sum_{i=1}^n f(X_i)$ (pathSum f n) and

$$v(f, H) = \lim_{n \to \infty} \frac1n \operatorname{Var}_H(S_n) \in [0, \infty].$$

The lag inner products are $\langle f, H^k f\rangle = \int f(x) \int f(y)\,H^k(x,dy)\,\pi(dx)$ (lagInner π H f k), and for $0 \le \lambda < 1$ the regularized variance is $v_\lambda(f,H) = \langle f,f\rangle + 2 \sum_{k\ge1} \lambda^k \langle f, H^k f\rangle$ (vLam π H f lam).

Formalization targets

Goal: Theorem 4 (p. 5)

Let $P_1, P_2$ be Markov kernels reversible with respect to $\pi$, $f \in L^2_0(\pi)$, and $P_1 \succeq P_2$. Then both asymptotic variances exist in $[0,\infty]$ and

$$v(f, P_1) \le v(f, P_2).$$

No rate, constant or regularity of the kernels is fixed. The statement is the ordering itself, valid for every reversible pair and every $f \in L^2_0(\pi)$.

Milestones, in attack order

1. Lemma 3 (p. 5): if $P_1, P_2$ have invariant distribution $\pi$ and $P_1 \succeq P_2$, then $P_2 - P_1$ is a positive operator on $L^2(\pi)$:

$$\iint f(x)f(y)\,\bigl(P_2(x,dy) - P_1(x,dy)\bigr)\,\pi(dx) \ge 0 \qquad (f \in L^2(\pi)).$$

2. A reversible kernel is a self-adjoint contraction on $L^2_0(\pi)$ (p. 5): $\langle Hf, g\rangle = \langle f, Hg\rangle$ and $\|Hf\| \le \|f\|$.
3. Finite-$n$ variance identity (p. 5), for $n \ge 1$:

$$\frac1n \operatorname{Var}_H(S_n) = \langle f,f\rangle + 2\sum_{i=1}^n \frac{n-i}{n}\,\langle f, H^i f\rangle.$$

4. Existence of $v(f,H)$ in $[0,\infty]$ (p. 6).
5. $v_\lambda(f,H) \to v(f,H)$ as $\lambda \uparrow 1$, finite or infinite (p. 6).
6. $v_\lambda(f,P_1) \le v_\lambda(f,P_2)$ for $0 \le \lambda < 1$ when $P_1 \succeq P_2$ (p. 6).

Significance

The result. Theorem 4 turns a pointwise, one-step comparison of kernels, which is easy to check, into a comparison of the quantity that governs Monte Carlo error. Its main consequence, drawn in §3 of the paper, is that the Metropolis–Hastings acceptance probability $\alpha_{MH}(x,y) = \min\{1, r(y,x)\}$ gives the maximal kernel in the off-diagonal order among reversible Metropolis–Hastings kernels with a given proposal. It is therefore optimal in asymptotic variance, on arbitrary state spaces. Proposition 5 of the same paper (a separate mission in this series) combines with it to show that a single Metropolis–Hastings kernel built on a mixture proposal beats the mixture of the component kernels. Later orderings of samplers (Mira 2001; Andrieu and Livingstone 2021) take this theorem as their base case.

Formalizing it. The theorem has been proved since 1998. No machine-checked version exists for general state spaces, and none of its milestones is on the platform. The formalization produces reusable infrastructure: the asymptotic variance of a stationary chain as an extended-real limit on Mathlib's Ionescu–Tulcea path measure, the $L^2$ facts for reversible kernels (self-adjointness, contraction, the covariance formula for path sums), and the positivity of $P_2 - P_1$ under off-diagonal domination. Each of these is used again in any formal treatment of MCMC efficiency or the Markov chain central limit theorem.

Difficulty

The direct approach compares the two finite-$n$ variances. This fails, and not just technically: the paper exhibits two doubly stochastic, symmetric $4\times4$ matrices with $P_1 \succeq P_2$ for which the variance of $f(X_0)+f(X_1)+f(X_2)$ is 15.4 under $P_1$ and 14.8 under $P_2$ (p. 7). Off-diagonal domination orders the lag-one covariances, but higher-order correlations "need not be ordered" (p. 5). The ordering appears only in the limit, and only for reversible kernels. The comparison must pass through an object that sees all lags at once and is monotone along the segment $P_1 + \beta(P_2 - P_1)$, and that object involves resolvents of operators on $L^2_0(\pi)$. The limit may be infinite, so every comparison must be made in $[0,\infty]$. Mathlib has neither the spectral measure of a self-adjoint operator nor the Kipnis–Varadhan theory.

Formalization scope

The state space is {E : Type*} [MeasurableSpace E] with [MeasurableSingletonClass E] wherever off-diagonal domination appears. This is an assumption the paper leaves implicit: $A \setminus \{x\}$ must be an event. $\pi$ is a probability measure and all kernels are Markov kernels. Reversibility is Mathlib's Kernel.IsReversible, invariance is Kernel.Invariant. The function $f$ is measurable with MemLp f 2 π and, for $L^2_0$, ∫ f ∂π = 0; measurability picks a representative of the $L^2$ class and costs nothing. The chain is Kernel.trajMeasure started from $\pi$. The sum runs over $X_1,\dots,X_n$, not $X_0$. Variances are Mathlib's evariance in $[0,\infty]$, and $v(f,H)$ is a Tendsto limit in ℝ≥0∞, so an infinite asymptotic variance is represented. The goal asserts the existence of both limits rather than assuming it, so it cannot hold vacuously. Neither it nor any milestone specializes to finite $E$, to Metropolis–Hastings kernels, or to a chain started from a point. Lemma 3 assumes invariance only, and the theorem requires reversibility, as printed.

Two statements depart in form from the page. The finite-$n$ variance identity and $v_\lambda$ are written through the moments $\langle f, H^k f\rangle$ (a Neumann series) instead of through the spectral measure $e_{f,H}$ and the resolvent $(I-\lambda H)^{-1}$. The two forms agree for a self-adjoint contraction, and this is noted in each item. The paper's appeal to the spectral theorem and to Kipnis and Varadhan (1986) is not restated as an item: a complete development needs it, or an equivalent argument, as part of the proof. Proofs of any milestone, and reusable lemmas on the path measure (stationarity and the marginal laws of $(X_i, X_j)$), are welcome.

Selected references

• L. Tierney, A Note on Metropolis–Hastings Kernels for General State Spaces, Ann. Appl. Probab. 8(1), 1–9, 1998. https://doi.org/10.1214/aoap/1027961031
• P. H. Peskun, Optimum Monte-Carlo sampling using Markov chains, Biometrika 60(3), 607–612, 1973. https://doi.org/10.1093/biomet/60.3.607
• C. Kipnis and S. R. S. Varadhan, Central limit theorem for additive functionals of reversible Markov processes and applications to simple exclusions, Comm. Math. Phys. 104, 1–19, 1986. https://doi.org/10.1007/BF01210789
• A. Mira, Ordering and improving the performance of Monte Carlo Markov chains, Statist. Sci. 16(4), 340–350, 2001. https://doi.org/10.1214/ss/1015346318
• C. Andrieu and S. Livingstone, Peskun–Tierney ordering for Markovian Monte Carlo: beyond the reversible scenario, Ann. Statist. 49(4), 1958–1981, 2021. https://doi.org/10.1214/20-AOS2008

Motivation

The Metropolis–Hastings algorithm (Metropolis et al. 1953; Hastings 1970) is the basic construction of Markov chain Monte Carlo. It turns a target probability distribution $\pi$, known only up to a constant, into a Markov chain that has $\pi$ as its invariant distribution. Bayesian computation depends on it, and a sampler is usually designed by proving one property: reversibility, or detailed balance, with respect to $\pi$.

For discrete state spaces, or when every measure involved has a density with respect to a common reference measure, the condition on the acceptance probability is the familiar $\pi(x)q(x,y)\alpha(x,y)=\pi(y)q(y,x)\alpha(y,x)$. Samplers used in practice often do not fit this setting: deterministic involutive proposals, mixtures of proposals with different supports, and Green's dimension-changing moves for model selection (Green 1995) have no common density. Each of these was treated separately in the literature. Tierney (1998) gives one necessary and sufficient condition that covers all of them, on an arbitrary measurable state space. This mission formalizes that condition.

Setting

Let $(E,\mathcal E)$ be a measurable space, with no topological or countability assumption. Let $\pi$ be a probability measure on $E$, the target. Let $Q(x,dy)$ be a Markov transition kernel on $E$, the proposal, and let $\alpha:E\times E\to[0,1]$ be a measurable function, the acceptance probability. From the current state $x$, a candidate $y$ is drawn from $Q(x,\cdot)$ and accepted with probability $\alpha(x,y)$; otherwise the chain stays at $x$. The resulting Metropolis–Hastings kernel (mhKernel Q α) is

$$P(x,dy)=Q(x,dy)\,\alpha(x,y)+\delta_x(dy)\int\bigl(1-\alpha(x,u)\bigr)\,Q(x,du).\tag{1}$$

A kernel $P$ satisfies detailed balance with respect to $\pi$ if the two measures $\pi(dx)P(x,dy)$ and $\pi(dy)P(y,dx)$ on $\mathcal E\otimes\mathcal E$ are equal (Eq. (2)). Equivalently, $\int_A P(x,B)\,\pi(dx)=\int_B P(x,A)\,\pi(dx)$ for all measurable $A,B$, which is Mathlib's Kernel.IsReversible P π.

Put $\mu(dx,dy)=\pi(dx)Q(x,dy)$ (π ⊗ₘ Q), and let $\mu^T(dx,dy)=\mu(dy,dx)$ be its image under the swap $(x,y)\mapsto(y,x)$. A symmetric split (IsSymmetricSplit) is a measurable set $R\subseteq E\times E$ with $(x,y)\in R\iff(y,x)\in R$, such that $\mu$ and $\mu^T$ are mutually absolutely continuous on $R$ and mutually singular on its complement $R^c$. Informally, $R$ consists of the pairs between which the proposal can move in both directions. A ratio version (IsRatioVersion) is a measurable $r:E\times E\to(0,\infty)$ that is a density of $\mu_R$ (the restriction of $\mu$ to $R$) with respect to $\mu^T_R$, and satisfies $r(x,y)=1/r(y,x)$ at every point.

Formalization targets

Goal: Theorem 2 (p. 3)

For every $\pi$, $Q$, $\alpha$ as above and every symmetric split $R$ and ratio version $r$ for $\mu=\pi\otimes Q$:

$$P\ \text{satisfies detailed balance w.r.t. }\pi\iff \begin{cases}\text{(i)}\ \ \alpha=0\quad\mu\text{-a.e. on }R^c,\\[2pt] \text{(ii)}\ \alpha(x,y)\,r(x,y)=\alpha(y,x)\quad\mu\text{-a.e. on }R.\end{cases}$$

The paper phrases the left side as condition (4), $\mu(dx,dy)\alpha(x,y)=\mu^T(dx,dy)\alpha(y,x)$, which its text identifies with (2) for the kernel (1). The goal states it for the kernel itself.

Milestones

1. Proposition 1 (p. 2). For every $\sigma$-finite measure $\mu$ on $E\times E$: a symmetric split $R$ exists; any two symmetric splits differ by a set null for both $\mu$ and $\mu^T$; and every symmetric split admits a ratio version.
2. §2, Eqs. (2)–(3) (p. 2). The kernel (1) satisfies (2) if and only if

$$\pi(dx)Q(x,dy)\alpha(x,y)=\pi(dy)Q(y,dx)\alpha(y,x),\tag{3}$$

that is, the rejection mass on the diagonal does not affect reversibility.

Companion item

§2, special case 1 (pp. 3–4). If $\pi(dx)=\pi(x)\nu(dx)$ and $Q(x,dy)=q(x,y)\nu(dy)$ for a $\sigma$-finite $\nu$, then $R=\{\pi(x)q(x,y)>0,\ \pi(y)q(y,x)>0\}$ is a symmetric split, $r=\pi(x)q(x,y)/(\pi(y)q(y,x))$ is a density of $\mu_R$ with respect to $\mu^T_R$, and detailed balance is equivalent to the two conditions holding $\nu\times\nu$-almost everywhere.

Significance

Theorem 2 lets reversibility be checked in the same way for every Metropolis–Hastings variant: compute $R$ and $r$ for the proposal, then verify (i) and (ii). The paper derives from it the reversibility of the standard acceptance probability $\alpha_{MH}=\min\{1,r(y,x)\}$ on $R$ (and $0$ off $R$), and the three special cases of §2 are instances. Together with Mathlib's Kernel.IsReversible.invariant, it yields that $\pi$ is invariant for the sampler. This is the correctness statement of every MCMC method built on the Metropolis–Hastings kernel. Missions II and III of this series (Peskun ordering; mixture proposals) take reversible Metropolis–Hastings kernels as their objects.

The result is proved in the paper. The mission adds a machine-checked proof at the paper's full generality: no densities, no dominating measure, no countability of $\mathcal E$. The platform currently has only finite-state statements (MarkovMixing.metropolis_stationary, a sufficiency direction on a Fintype state space with a matrix proposal), so neither the general kernel nor the converse direction is formalized there.

Difficulty

The obvious argument works with densities: write both sides of (3) as densities with respect to one reference measure and compare them pointwise. On a general space no such reference is given for $\mu$ and $\mu^T$ together, and even $\mu+\mu^T$ yields densities only up to null sets. Pointwise comparison of densities is therefore not available, and the statement mixes three kinds of almost-everywhere claim ($\mu$-a.e., $\mu^T$-a.e., and a.e. for the restrictions to $R$ and $R^c$). The Lean statement also has to hold for every version of $R$ and $r$, not one convenient choice. Milestone 2 has its own content: the diagonal part of $P$ is a measure concentrated on the diagonal, which need not be a measurable set, and its symmetry has to be shown without that measurability.

Formalization scope

• Space. {E : Type*} [MeasurableSpace E] with nothing else: no measurable singletons, no topology, no countable generation. π : Measure E with [IsProbabilityMeasure π] and Q : Kernel E E with [IsMarkovKernel Q].
• Acceptance probability. α : E × E → ℝ≥0∞ with the hypotheses Measurable α and ∀ p, α p ≤ 1, which is the paper's measurable $\alpha:E\times E\to[0,1]$. The kernel is Q.withDensity (fun x y => α (x, y)) + Kernel.withDensity Kernel.id (fun x _ => ∫⁻ u, (1 - α (x, u)) ∂(Q x)). The measurability hypothesis rules out the junk zero kernel that Kernel.withDensity returns for a non-measurable density.
• Detailed balance is Kernel.IsReversible. It agrees with the measure identity (2) because rectangles determine a finite measure on $\mathcal E\otimes\mathcal E$.
• (i) and (ii) are almost-everywhere statements for the restrictions of $\mu$ to $R^c$ and to $R$ respectively; neither is required pointwise. Condition (ii) off $R$ would be false in general.
• $R$ and $r$ are universally quantified in the goal. A formalization that fixes one specific Radon–Nikodym derivative, or drops the everywhere conditions $0<r<\infty$, $r(x,y)=1/r(y,x)$, proves a different statement. Proposition 1's existence clause shows that the goal's hypotheses can be met, so the goal is not vacuous. A sorry-free check in the workspace confirms this for $Q=\pi$ with $R=E\times E$, $r\equiv1$.
• Added hypothesis. The companion item assumes $\nu$ is $\sigma$-finite, which the paper leaves implicit in "$\nu\times\nu$-almost all".
• Infrastructure. The definitions IsSymmetricSplit and IsRatioVersion (a symmetric Lebesgue-type decomposition of a measure against its transpose) are reusable for any reversibility argument on product spaces. Lemmas on Measure.map Prod.swap of compProd and withDensity, and on the symmetry of measures carried by the diagonal, are also welcome contributions.

Selected references

• L. Tierney, A Note on Metropolis–Hastings Kernels for General State Spaces, Ann. Appl. Probab. 8(1) (1998) 1–9. https://doi.org/10.1214/aoap/1027961031
• N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E. Teller, Equation of State Calculations by Fast Computing Machines, J. Chem. Phys. 21 (1953) 1087–1092. https://doi.org/10.1063/1.1699114
• W. K. Hastings, Monte Carlo Sampling Methods Using Markov Chains and Their Applications, Biometrika 57 (1970) 97–109. https://doi.org/10.1093/biomet/57.1.97
• L. Tierney, Markov Chains for Exploring Posterior Distributions, Ann. Statist. 22 (1994) 1701–1762. https://doi.org/10.1214/aos/1176325750
• P. J. Green, Reversible Jump Markov Chain Monte Carlo Computation and Bayesian Model Determination, Biometrika 82 (1995) 711–732. https://doi.org/10.1093/biomet/82.4.711

Motivation

A variational inequality asks for a point of a convex set at which a vector field points "inward" against every feasible direction. The format covers the first-order optimality conditions of constrained optimization, nonlinear complementarity problems, traffic and economic equilibria (Wardrop, Walrasian, Nash–Cournot), and systems of nonlinear equations; see Harker and Pang's survey (Math. Programming 48, 1990) and Facchinei and Pang's monograph (Springer, 2003).

When the map has no special structure (not strongly monotone, not Lipschitz with known constant, not affine) and the feasible set is a general closed convex set, the practical algorithms are projection methods. The oldest is Korpelevich's extragradient method (1976). Without a known Lipschitz constant, extragradient-type methods need a linesearch in which every trial point costs one projection onto the feasible set, and projection onto a general convex set is itself an optimization problem.

Solodov and Svaiter (SIAM J. Control Optim. 37 (1999) 765–776) proposed a method that spends exactly two projections per iteration, whatever the linesearch does, and proved global convergence under only continuity of the map and a generalized monotonicity condition weaker than pseudomonotonicity. The method, often called the hyperplane projection method, is a standard reference point for later projection and extragradient-type algorithms.

Timeline:

• 1976: Korpelevich, extragradient method, Lipschitz monotone maps.
• 1987–1994: Khobotov (1987), Iusem (1994) and others: extragradient variants with Armijo-type stepsize rules, which need one projection per trial step.
• 1997: Iusem and Svaiter, a separating-hyperplane variant of extragradient for monotone maps (reference [9] of the paper).
• 1999: Solodov and Svaiter, Algorithm 2.1: two projections per iteration, convergence under condition (1.2) below.

Setting

Work in $\mathbb{R}^n$ with the Euclidean inner product $\langle\cdot,\cdot\rangle$ and norm $\|\cdot\|$. Let $C \subseteq \mathbb{R}^n$ be closed and convex and $F : \mathbb{R}^n \to \mathbb{R}^n$ continuous. The problem $\mathrm{VI}(F, C)$ is to find $x^*$ with

$$x^* \in C, \qquad \langle F(x^*), x - x^*\rangle \ge 0 \quad \text{for all } x \in C. \tag{1.1}$$

Its solution set is $S$. The projection onto a nonempty closed convex set $K$ is $P_K[x] := \arg\min_{y \in K}\|y - x\|$. The projected residual is $r(x) := x - P_C[x - F(x)]$; its zeros are exactly the points of $S$.

Condition (1.2) requires, for every $x^* \in S$,

$$\langle F(x), x - x^*\rangle \ge 0 \qquad \text{for all } x \in C. \tag{1.2}$$

It holds when $F$ is monotone or pseudomonotone, and in cases where $F$ is neither.

Algorithm 2.1. Fix $\gamma, \sigma \in (0,1)$ and $x^0 \in C$. Given $x^i$: if $r(x^i) = 0$, stop. Otherwise let $k_i$ be the smallest nonnegative integer $k$ with

$$\langle F(x^i - \gamma^k r(x^i)), r(x^i)\rangle \ge \sigma\|r(x^i)\|^2, \tag{2.1}$$

set $\eta_i = \gamma^{k_i}$, $z^i = x^i - \eta_i r(x^i)$, $H_i = \{x \mid \langle F(z^i), x - z^i\rangle \le 0\}$, and

$$x^{i+1} = P_{C \cap H_i}[x^i].$$

The hyperplane $\partial H_i$ separates $x^i$ from $S$.

Formalization targets

Goal: Theorem 2.1

If $C$ is closed and convex, $F$ is continuous, $S \ne \emptyset$ and (1.2) holds, then every sequence generated by Algorithm 2.1 converges to a single point of $S$:

$$\exists\, \hat x \in S:\quad x^i \to \hat x \quad (i \to \infty).$$

The theorem fixes no rate and no constant; it asserts convergence of the whole sequence, not only of a subsequence.

Milestones (in attack order)

• Lemma 2.1 (p. 768): for nonempty closed convex $B$, $\langle x - P_B[x], z - P_B[x]\rangle \le 0$ for $z \in B$, and $\|P_B[x]-P_B[y]\|^2 \le \|x-y\|^2 - \|P_B[x]-x+y-P_B[y]\|^2$.
• Residual characterization (p. 767): $x \in S \iff r(x) = 0$.
• (2.5) (p. 769): $\langle F(x), r(x)\rangle \ge \|r(x)\|^2$ for $x \in C$.
• Linesearch well-definedness (p. 769): for $x \in C$ with $r(x) \ne 0$, some $k$ satisfies (2.1).
• Lemma 2.2 (p. 768): $x^{i+1} = P_{C\cap H_i}[\bar x^i]$ with $\bar x^i = P_{H_i}[x^i]$.
• (2.6) (pp. 769–770): $\|x^{i+1}-x^*\|^2 \le \|x^i-x^*\|^2 - \|x^{i+1}-\bar x^i\|^2 - \big(\eta_i\sigma/\|F(z^i)\|\big)^2\|r(x^i)\|^4$ for every $x^* \in S$.
• (2.8) (p. 770): $\eta_i\|r(x^i)\| \to 0$.

Significance

Theorem 2.1 gives global convergence of a projection method for variational inequalities with no Lipschitz constant, no monotonicity and no knowledge of the problem beyond continuity and (1.2), at a fixed cost of two projections per iteration. Condition (1.2) covers pseudomonotone maps, which arise as gradients of pseudoconvex functions and in equilibrium models where monotonicity fails. The separating-hyperplane-and-project template of the proof is reused throughout the later literature on projection, proximal and hybrid methods for monotone inclusions.

The result has been proved since 1999. To the best of available knowledge no machine-checked proof of it, or of any convergence theorem for a projection method for variational inequalities, exists in Lean or Mathlib. This mission produces the statement and the supporting layer: a Euclidean projection onto closed convex sets with its standard inequalities, variational inequality solution sets, the projected residual, and a formal model of an Armijo-type linesearch algorithm with termination.

Difficulty

The Fejér-type inequality (2.6) quickly gives bounded iterates and $\eta_i\|r(x^i)\| \to 0$. The obvious next step, concluding $r(x^i) \to 0$, fails: nothing prevents the stepsizes $\eta_i$ from tending to zero, and in that regime the product going to zero says nothing about the residual. This regime is where the minimality of $k_i$ and the continuity of $F$ enter, and it is the step a naive formalization (for instance one that drops minimality, or fixes the stepsize) cannot reach. A second subtlety is that (1.2) is needed at an accumulation point that is only known to lie in $S$ at the end of the argument, which is why the condition must hold for every $x^* \in S$. Finally, subsequential convergence must be upgraded to convergence of the whole sequence to one solution; convergence of a subsequence, or of the distance to $S$, is strictly weaker.

Formalization scope

• $\mathbb{R}^n$ is EuclideanSpace ℝ (Fin n) (not Fin n → ℝ, whose norm is the sup norm). The accumulation-point step needs finite dimension; no Hilbert-space generalization is intended.
• Projection encoding. projOnto K x is a nearest point of $K$ to $x$ when one exists, chosen by Classical.choose, and the junk value $x$ otherwise. On nonempty closed convex sets it is exactly $P_K[x]$; the paper only projects onto such sets ($C$, $H_i$, $C \cap H_i$), so the junk value is never reached under the hypotheses.
• Stopping-rule encoding. A run is a sequence x : ℕ → ℝⁿ with Armijo indices k : ℕ → ℕ (predicate IsAlg21Run). If $r(x^i) = 0$ the method has stopped and the run stalls, $x^{i+1} = x^i$; otherwise $k_i$ is the least index satisfying (2.1) and $x^{i+1} = P_{C \cap H_i}[x^i]$. A stalled point is a solution, so finitely terminating runs are included in the goal.
• Parameters $\gamma, \sigma$ are real with $0 < \gamma < 1$, $0 < \sigma < 1$, universally quantified; $n$, $C$, $F$ and $x^0 \in C$ are arbitrary.
• (2.6) is stated for one generic step ($x \in C$, $r(x)\ne 0$, $k$ satisfying (2.1)) rather than along a run; it is the same inequality with $x^i, k_i$ abstracted.
• Trivializing formalizations are ruled out: condition (1.2) is quantified over every solution and every $x \in C$ (not replaced by monotonicity or an existential), $k_i$ is the least index satisfying (2.1), the update projects $x^i$ onto $C \cap H_i$ (not onto $C$ alone), the stopped case is pinned down by the stall encoding, and the conclusion is convergence of the whole sequence to one solution, not $r(x^i) \to 0$ or $\operatorname{dist}(x^i, S) \to 0$.
• Needed infrastructure: existence, uniqueness and variational characterization of the projection (Mathlib has exists_norm_eq_iInf_of_complete_convex and norm_eq_iInf_iff_real_inner_le_zero), firm nonexpansiveness, the explicit projection onto a halfspace, and a bounded-sequence subsequence argument in $\mathbb{R}^n$. The projection lemmas are reusable for any projection-type method; contributions proving them as standalone lemmas are welcome.

Selected references

• M. V. Solodov and B. F. Svaiter, A New Projection Method for Variational Inequality Problems, SIAM J. Control Optim. 37(3), 765–776, 1999. https://doi.org/10.1137/S0363012997317475
• G. M. Korpelevich, The extragradient method for finding saddle points and other problems, Matecon 12, 747–756, 1976.
• A. N. Iusem and B. F. Svaiter, A variant of Korpelevich's method for variational inequalities with a new search strategy, Optimization 42, 309–321, 1997. https://doi.org/10.1080/02331939708844365
• P. T. Harker and J.-S. Pang, Finite-dimensional variational inequality and nonlinear complementarity problems: a survey of theory, algorithms and applications, Math. Programming 48, 161–220, 1990. https://doi.org/10.1007/BF01582255
• F. Facchinei and J.-S. Pang, Finite-Dimensional Variational Inequalities and Complementarity Problems, Springer, 2003. https://doi.org/10.1007/b97543

Motivation

Shor's 1997 paper (SIAM J. Comput. 26(5), arXiv:quant-ph/9508027) gives a polynomial-time quantum algorithm for factoring integers. The quantum computer does not factor directly: it finds the multiplicative order of an element modulo $n$. The step from order finding to factoring is classical and randomized, and goes back to Miller's 1976 work on primality testing (G. L. Miller, Riemann's hypothesis and tests for primality, J. Comput. System Sci. 13 (1976)). Every account of Shor's algorithm, and every resource estimate for breaking RSA with a quantum computer, depends on this reduction succeeding with a constant probability per trial. This mission formalizes that probability bound as Shor states it on p. 1498 of the published paper.

Setting

Let $n > 1$ be an odd integer with prime factorization

$$n = \prod_{i=1}^{k} p_i^{\alpha_i},$$

so $k$ is the number of distinct prime factors of $n$, all odd. The unit group $(\mathbb{Z}/n\mathbb{Z})^\times$ consists of the residues coprime to $n$; it has $\varphi(n)$ elements, where $\varphi$ is Euler's totient function.

For a unit $x$ the order $r = \operatorname{ord}_n(x)$ is the least positive integer with $x^r \equiv 1 \pmod n$. For each $i$ the local order $r_i$ is the order of $x \bmod p_i^{\alpha_i}$, taken modulo the full prime power, not modulo $p_i$. For a positive integer $m$, $\nu_2(m)$ denotes the exponent of the largest power of $2$ dividing $m$.

The reduction is: choose $x$ uniformly at random from $(\mathbb{Z}/n\mathbb{Z})^\times$, obtain its order $r$ (from the quantum subroutine), and compute

$$g(x) = \gcd\bigl(x^{r/2} - 1,\ n\bigr).$$

The procedure yields a nontrivial factor at $x$ when $r$ is even and $1 < g(x) < n$. In Lean this event is ShorAlgorithms.Reduction.successEvent n u for u : (ZMod n)ˣ, and $r_i$ is localOrder n u p for p ∈ n.primeFactors.

Formalization targets

Goal: the success probability

$$\Pr_{x \in (\mathbb{Z}/n)^\times}\bigl[r \text{ even and } 1 < \gcd(x^{r/2}-1, n) < n\bigr] \;\ge\; 1 - \frac{1}{2^{k-1}}.$$

It is stated for every odd $n > 1$. For a prime power ($k = 1$) the bound is $0$, so the statement says nothing there; it is informative exactly when $n$ is not a prime power, as the paper remarks. The constant is sharp: for $n = 21$ exactly $6$ of the $12$ units succeed, so $1 - 1/2^{k}$ in place of $1 - 1/2^{k-1}$ would be false.

Milestones, in the order the page uses them

1. Success criterion. If $r$ is even and $x^{r/2} \not\equiv -1 \pmod n$, then $1 < \gcd(x^{r/2}-1, n) < n$.
2. Order is the lcm. $r = \operatorname{lcm}(r_1, \dots, r_k)$.
3. Failure forces agreement. For odd $n$, if the procedure fails at $x$, then $\nu_2(r_1) = \cdots = \nu_2(r_k)$.
4. At most half per odd prime power. For an odd prime $p$ and $\alpha \ge 1$, at most $\varphi(p^\alpha)/2$ units modulo $p^\alpha$ have order with a prescribed 2-adic valuation.
5. All agree rarely. The units for which $\nu_2(r_1) = \cdots = \nu_2(r_k)$ number at most $\varphi(n)/2^{k-1}$.

Significance

The result. The bound turns an order-finding oracle into a factoring algorithm: when $n$ is odd and not a prime power, each trial succeeds with probability at least $1/2$, so $t$ independent trials all fail with probability at most $2^{-t}$. Even numbers and prime powers are split classically, as the paper notes, so the bound completes the reduction from factoring to order finding. The same criterion — a square root of $1$ other than $\pm 1$ splits $n$ — underlies the Miller–Rabin test and several classical factoring methods.

Formalizing it. The mathematics is classical and proved; the paper gives a sketch of one paragraph. This mission writes out the sketch as machine-checked statements over Mathlib's ZMod, including the probabilistic step, which in the paper is an informal appeal to the Chinese remainder theorem and "50% probability of agreeing with the previous ones". Mathlib already has the needed ingredients (cyclicity of $(\mathbb{Z}/p^\alpha)^\times$ for odd $p$, ZMod.chineseRemainder, ZMod.card_units_eq_totient), but not the reduction or its probability bound.

Difficulty

The success criterion (milestone 1) is elementary. The substance is the counting. The obvious route — treating the $\nu_2(r_i)$ as independent and each "equal to the previous one with probability $1/2$" — needs both a precise product decomposition of the unit group modulo $n$ into the unit groups modulo $p_i^{\alpha_i}$, compatible with the local orders, and the count in a cyclic group of even order of the elements whose order has a given 2-adic valuation. The informal phrase "at most a 50% probability of agreeing with the previous ones" hides a conditioning argument over $k - 1$ coordinates that has to be done by an explicit cardinality bound. A second pitfall is milestone 3: its converse direction and its forward direction use oddness of $n$ in different places, and modulo a power of $2$ the argument breaks because $-1 \equiv 1 \pmod 2$.

Formalization scope

• Sample space. Uniform on (ZMod n)ˣ; probabilities are stated in cleared-denominator form, $(1 - 2^{-(k-1)})\,\varphi(n) \le \#\{\text{successes}\}$ in $\mathbb{R}$, with the count as Nat.card of a subtype. Non-units have no multiplicative order and are not sampled.
• The gcd. $x^{r/2}$ is represented by its least nonnegative residue .val, which is at least $1$ for a unit when $n > 1$, so the natural-number subtraction in val - 1 never truncates. $r/2$ is natural-number division, used only under Even r.
• $k$. n.primeFactors.card, at least $1$ for $n > 1$, so k - 1 does not truncate. Since $n$ is odd this equals the page's "number of distinct odd prime factors".
• Local orders. The order of the image of $x$ in ZMod (p ^ n.factorization p) under the reduction homomorphism.
• Hypotheses. The goal assumes exactly $n$ odd and $n > 1$. It does not assume "not a prime power": that clause in the paper describes when the bound is useful. Milestones 1 and 2 do not assume $n$ odd, because they do not need it; milestones 3 and 5 do.
• No trivialization. Counting over all of ZMod n instead of the units would put non-units (with junk order $0$) into the denominator; the goal counts over (ZMod n)ˣ and divides by $\varphi(n)$. The goal's constant is the paper's $1 - 1/2^{k-1}$, which is attained, so it cannot be weakened into a triviality without changing the theorem.
• Welcome contributions. A reusable counting lemma for elements of prescribed 2-adic order in a finite cyclic group; the transfer of ZMod.chineseRemainder to unit groups and to local orders; and proofs of the milestones in any order.

Selected references

• P. W. Shor, Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer, SIAM J. Comput. 26(5):1484–1509, 1997. https://doi.org/10.1137/S0097539795293172 (preprint arXiv:quant-ph/9508027, https://arxiv.org/abs/quant-ph/9508027)
• G. L. Miller, Riemann's hypothesis and tests for primality, J. Comput. System Sci. 13(3):300–317, 1976. https://doi.org/10.1016/S0022-0000(76)80043-8
• D. E. Knuth, The Art of Computer Programming, Vol. 2: Seminumerical Algorithms, 2nd ed., Addison-Wesley, 1981.
• G. H. Hardy and E. M. Wright, An Introduction to the Theory of Numbers, 5th ed., Oxford University Press, 1979 (Theorem 121, Chinese remainder theorem).

Motivation

Shor's factoring and discrete logarithm algorithms (Shor 1997) reduce both problems to sampling from the output of a quantum Fourier transform: a register holding a superposition with a hidden period is transformed, then measured, and the measured value carries information about the period. The whole speed-up rests on one engineering fact: for $q = 2^l$ the $q \times q$ Fourier matrix, which has $q^2 = 4^l$ entries, can be applied by a quantum circuit of only $O(l^2)$ elementary gates, each acting on one or two bits.

That circuit was found independently by Coppersmith (IBM RC 19642, 1994) and Deutsch, and Shor's §4 presents it following Ekert and Jozsa (Rev. Mod. Phys. 68, 1996). It is the quantum analogue of the radix-2 fast Fourier transform, and it reappears in phase estimation, in the hidden subgroup algorithms for abelian groups, and in every textbook account of quantum computation. This mission formalizes Shor's statement that the circuit computes the Fourier matrix, up to a reversal of the output bits, together with the two displayed steps of its verification.

Setting

A register of $l$ bits has one basis state $|a\rangle = |a_{l-1} a_{l-2} \dots a_0\rangle$ for every bit string, with $a_0$ the least significant bit; the string encodes the integer

$$a = \sum_{j=0}^{l-1} 2^j a_j, \qquad 0 \le a < q = 2^l .$$

A state is a vector of complex amplitudes, one per basis state. A gate is a matrix whose rows are indexed by input basis vectors and whose columns are indexed by output basis vectors (§2, p. 1489). A gate on some of the bits acts on those bits through its matrix and leaves the others alone (§2, p. 1490): the output amplitude at a string $b$ is the sum, over the possible input values $u$ of the acted-on bits, of the input amplitude at $b$ with those bits replaced by $u$, times the matrix entry from $u$ to the corresponding bits of $b$.

The Fourier matrix $A_q$ (eq. (4.1)) is the $q \times q$ matrix with $(a, c)$ entry $q^{-1/2}\exp(2\pi i\,ac/q)$; it takes $|a\rangle$ to $q^{-1/2}\sum_{c=0}^{q-1}\exp(2\pi i\,ac/q)\,|c\rangle$.

The circuit uses two gates (eqs. (4.2), (4.3)):

• $R_j$ acts on bit $j$ with matrix $\frac{1}{\sqrt 2}\begin{pmatrix} 1 & 1 \\ 1 & -1\end{pmatrix}$;
• $S_{j,k}$, for $j < k$, acts on bits $j$ and $k$ with matrix $\mathrm{diag}(1, 1, 1, e^{i\theta_{k-j}})$, where $\theta_{k-j} = \pi / 2^{k-j}$; it multiplies the amplitude of a basis state by $e^{i\theta_{k-j}}$ when bits $j$ and $k$ are both $1$.

The circuit is the gate sequence (4.4), applied from left to right:

$$R_{l-1}\, S_{l-2,l-1}\, R_{l-2}\, S_{l-3,l-1}\, S_{l-3,l-2}\, R_{l-3} \cdots R_1\, S_{0,l-1}\, S_{0,l-2} \cdots S_{0,2}\, S_{0,1}\, R_0 ,$$

that is, for $j = l-1, \dots, 0$ it applies $S_{j,l-1}, \dots, S_{j,j+1}$ and then $R_j$. On three bits it is $R_2 S_{1,2} R_1 S_{0,2} S_{0,1} R_0$.

The bit reversal of a string $b$ is the string $c$ with $c_k = b_{l-1-k}$.

Formalization targets

Goal: the circuit computes $A_q$ up to bit reversal (§4, p. 1496)

For every $l \ge 0$ and every basis state $|a\rangle$, with $q = 2^l$,

$$\text{circuit}\,|a\rangle = \frac{1}{q^{1/2}}\sum_{b}\exp(2\pi i\,ac/q)\,|b\rangle, \qquad c = \text{bit reversal of } b .$$

The sum runs over all $l$-bit strings $b$. Reading the output register in reverse order therefore yields $A_q|a\rangle$.

Milestone 1: the amplitude along the circuit (§4, eq. (4.5), p. 1496)

The amplitude of $|b\rangle$ in $\text{circuit}\,|a\rangle$ is

$$2^{-l/2}\exp\Big(i\Big(\sum_{0\le j<l}\pi a_jb_j + \sum_{0\le j<k<l}\frac{\pi}{2^{k-j}}a_jb_k\Big)\Big).$$

Milestone 2: the phase identity (§4, eqs. (4.6)–(4.10), pp. 1496–1497)

For bit strings $a, b$, with $c$ the bit reversal of $b$ and $a, c$ their values,

$$\exp\Big(i\Big(\sum_{0\le j<l}\pi a_jb_j + \sum_{0\le j<k<l}\frac{\pi}{2^{k-j}}a_jb_k\Big)\Big) = \exp(2\pi i\,ac/q).$$

Significance

The result. The goal says that $A_q$, a dense unitary on $2^l$ amplitudes, is realized by $l$ one-bit gates and $l(l-1)/2$ two-bit gates, followed by a relabelling of the output. This is what makes the Fourier sampling step of the factoring algorithm (§5) and of the discrete logarithm algorithm (§6) polynomial in the number of bits. Without it, the analyses of those sections describe measurements of states that no efficient circuit is known to prepare. The bit reversal is a real part of the statement: the circuit does not compute $A_q$ itself for $l \ge 2$, and an implementation must either permute the output bits or read them in reverse order.

Formalizing it. The identity is classical and its proof is short on paper; its content lies in bookkeeping that is easy to get wrong: which bit a gate touches, which of the input or output value a bit holds when a phase gate acts, and the order in which the gates are applied. A machine-checked version fixes all of these conventions explicitly and yields a reusable model of gate-level circuits on bit strings. To the drafter's knowledge there is no Lean formalization of this circuit on Prove2Me; the platform's FastFourierTransform rows state the classical Cooley–Tukey recursion for the unnormalized transform with the opposite sign, which is a different object.

Difficulty

Multiplying out the gate matrices is not an argument beyond tiny $l$: the difficulty is to control the whole product of $l(l+1)/2$ gates symbolically. The paper's argument (§4, p. 1496) is informal at exactly the points a formal proof must make precise: that only one sequence of intermediate basis states from $|a\rangle$ to $|b\rangle$ carries nonzero amplitude, and that at the moment $S_{j,k}$ acts, bit $k$ already holds its output value $b_k$ while bit $j$ still holds its input value $a_j$. Both facts depend on the order (4.4); applying the same gates in the reverse order gives a different unitary, whose phase pairs $b_j$ with $a_k$. The phase identity then holds only modulo $2\pi$, not as an equality of real numbers, so it cannot be closed by rearranging sums alone.

Formalization scope

• States. A bit string on $l$ bits is Fin l → Fin 2, with bit j equal to $a_j$ and value $\sum_j 2^j a_j$ (least significant bit first). A state is a function from bit strings to ℂ. The case $l = 0$ ($q = 1$, empty circuit) is included.
• Gate application follows the row = input convention: a one-bit gate $M$ on bit $j$ sends $\psi$ to $b \mapsto \sum_u \psi(b[j\mapsto u])\,M_{u, b_j}$, and a two-bit gate acts analogously on the pair (bit $j$, bit $k$).
• The circuit is the gate list (4.4) — an explicit list of constructors R j and S j k — run by a left fold, so the leftmost gate acts first. It is not defined as the matrix $A_q$ or by its entries; a definition of that kind would make the goal true by unfolding and is ruled out.
• Angles. $\theta_{k-j} = \pi/2^{k-j}$ uses natural-number subtraction, which is exact because the list only contains $S_{j,k}$ with $j < k$.
• Normalization. The prefactor $q^{-1/2}$ is written $(\sqrt{2^l})^{-1}$.
• Basis form. The goal is stated on basis states, as on the page; by linearity it determines the circuit on every state.
• Not stated. The gate count: the paper's sentence "we thus need to use $l(l-1)/2$ quantum gates" (p. 1496) counts only the gates $S_{j,k}$; the sequence (4.4) also contains the $l$ gates $R_j$, $l(l+1)/2$ gates in all. The approximate transform of Coppersmith and the one-bit construction of Griffiths and Niu (p. 1497) are out of scope, as is the polynomial-time claim.

Welcome contributions: proofs of the two milestones and the goal; general lemmas about one- and two-bit gate actions on Fin l → Fin 2 states (commutation of gates on disjoint bits, linearity, action on basis states), which are reusable for any gate-level circuit; and a proof that the circuit is unitary.

Selected references

• P. W. Shor, Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer, SIAM J. Comput. 26(5):1484–1509, 1997. https://doi.org/10.1137/S0097539795293172
• D. Coppersmith, An approximate Fourier transform useful in quantum factoring, IBM Research Report RC 19642, 1994. https://arxiv.org/abs/quant-ph/0201067
• A. Ekert and R. Jozsa, Quantum computation and Shor's factoring algorithm, Rev. Mod. Phys. 68:733–753, 1996. https://doi.org/10.1103/RevModPhys.68.733
• R. B. Griffiths and C.-S. Niu, Semiclassical Fourier transform for quantum computation, Phys. Rev. Lett. 76:3228–3231, 1996. https://doi.org/10.1103/PhysRevLett.76.3228

Motivation

Conjugate gradient methods minimize a smooth function $f:\mathbb R^n\to\mathbb R$ using only gradients and a handful of stored vectors. This makes them the standard choice when $n$ is too large for Newton or quasi-Newton methods, which store an $n\times n$ matrix. On a strictly convex quadratic with exact line searches the classical method of Hestenes and Stiefel terminates in at most $n$ steps. On general functions, and with the inexact line searches used in practice, its behaviour is much less clear.

D. F. Shanno's 1978 paper in Mathematics of Operations Research (doi:10.1287/moor.3.3.244) links conjugate gradient methods to quasi-Newton methods. It writes the search direction as $-\hat H g$, where $\hat H$ is a positive definite approximation of the inverse Hessian that is never stored. The resulting "memoryless" BFGS directions give descent without exact line searches. The paper's new algorithm uses two BFGS updates: one from the last restart and one from the current step. Its first update is scaled by the Oren–Spedicato factor $\gamma_t$. Shanno and Phua's CONMIN code implements the algorithm, and the memoryless BFGS direction is the one-pair case of the later limited-memory BFGS methods.

• 1952: Hestenes and Stiefel, linear conjugate gradients.
• 1964: Fletcher and Reeves, nonlinear conjugate gradients.
• 1969: Polak and Ribière, a second nonlinear variant.
• 1972: Beale, a restart procedure that keeps the computed direction $d_t$.
• 1977: Powell's restart criterion (Powell 1977).
• 1978: Shanno's reformulation as memoryless and two-update quasi-Newton methods (this paper).

Setting

Vectors are columns in $\mathbb R^n$. A prime denotes transpose: $u'v$ is the inner product and $uv'$ the outer product. An iterative method produces points $x_k$, steps $p_k = x_{k+1}-x_k = \alpha_k d_k$ along search directions $d_k$, gradients $g_k = \nabla f(x_k)$, and gradient changes $y_k = g_{k+1}-g_k$. A line search is exact when $p_k'g_{k+1} = 0$.

The BFGS update of a matrix $H$ with the pair $(p,y)$ is

$$H^+ = H - \frac{p\,y'H + H y\,p'}{p'y} + \left(1+\frac{y'Hy}{p'y}\right)\frac{pp'}{p'y}.$$

A restart cycle begins at iteration $t$. At a later iteration $k>t$, Shanno's self-scaled restart matrix is

$$\hat H_k = \gamma_t\left(I - \frac{p_ty_t'+y_tp_t'}{p_t'y_t} + \frac{y_t'y_t}{p_t'y_t}\frac{p_tp_t'}{p_t'y_t}\right) + \frac{p_tp_t'}{p_t'y_t}, \qquad \gamma_t = \frac{p_t'y_t}{y_t'y_t}.$$

The matrix $\hat H_{k+1}$ is its BFGS update with $(p_k,y_k)$, and the self-scaled two-update direction is $d_{k+1} = -\hat H_{k+1}g_{k+1}$. The unscaled variant uses the BFGS update of $I$ in place of the first matrix.

Beale's direction is

$$d_{k+1} = -g_{k+1} + \frac{y_k'g_{k+1}}{d_k'y_k}d_k + \frac{y_t'g_{k+1}}{d_t'y_t}d_t.$$

The quadratic case has gradient $g(x) = Ax + c$ with $A$ symmetric positive definite.

Formalization targets

Goal: reduction of the self-scaled method to Beale's method

Let $A$ be symmetric positive definite, $g_i = Ax_i + c$ and $t<k$. Assume that for $t\le i\le k$ we have $x_{i+1} = x_i + p_i$, $p_i = \alpha_i d_i$ and $p_i'g_{i+1}=0$, that $p_t'Ap_i = 0$ for $t<i\le k$, and that $p_t, p_k \ne 0$. Then

$$-\hat H_{k+1}g_{k+1} = \gamma_t\left(-g_{k+1} + \frac{y_k'g_{k+1}}{d_k'y_k}d_k + \frac{y_t'g_{k+1}}{d_t'y_t}d_t\right).$$

This is the paper's claim that "for $f(x)$ quadratic with exact searches each of the above methods reduces exactly to Beale's method defined by (28)", with the conclusion (44). The scale is exactly $\gamma_t$.

Companion statements

• The unscaled two-update direction equals Beale's direction exactly.
• Both two-update directions are descent directions, $g_{k+1}'d_{k+1} < 0$, whenever $p_t'y_t > 0$ and $p_k'y_k > 0$. No exact search is needed.

Milestones on the path

• (34): the expansion of $-\hat H_{k+1}g_{k+1}$.
• (40): its form under an exact search.
• (41): gradients along a run on a quadratic.
• $p_t'g_{k+1} = 0$.
• (38), corrected by a factor 2: the action of the self-scaled restart matrix.
• (42): its form when $p_t'g_{k+1} = 0$.
• (43): the direction after substitution.

Significance

The result. The reduction says that the new algorithm reproduces Beale's restarted conjugate gradient directions on a quadratic with exact line searches. So it keeps the finite-termination and rate-of-convergence properties behind Beale's restart. Away from that setting it behaves as a quasi-Newton method, whose directions are descent directions under any line search with $p'y>0$. The two regimes are what justify relaxing the line search, which the paper's computations exploit. The scale $\gamma_t$ changes only the length of the step, not its direction.

Formalizing it. The claim is proved in the paper by a short computation, and no machine-checked version is known. This mission produces:

• a checked statement of the claim with every hypothesis explicit, including the conjugacy the proof takes as known;
• a corrected version of display (38), which is misprinted;
• reusable definitions of the additive BFGS update and of Beale's direction.

Difficulty

The obvious attempt is to expand both BFGS updates symbolically and compare with Beale's formula. This fails without two facts that are not algebraic identities. The first is that the restart step stays orthogonal to every later gradient, $p_t'g_{k+1}=0$. It needs the affine gradient of a quadratic, the exact search at the restart step, and conjugacy of $p_t$ with all later steps. The second is $y_k'p_t = 0$, which again comes from conjugacy. Beale's formula also has to be matched in its $d$-form: the coefficient $\frac{y'g}{d'y}d$ equals $\frac{y'g}{p'y}p$ only when the step length is nonzero. The descent statements need a different argument: the BFGS update of a positive definite matrix with $p'y>0$ must be shown to remain positive definite, and this has to be done twice.

Formalization scope

Vectors are Fin n → ℝ and matrices Matrix (Fin n) (Fin n) ℝ. The inner product $u'v$ is u ⬝ᵥ v, the outer product $uv'$ is vecMulVec u v, and $Hv$ is H *ᵥ v. The quadratic enters only through its gradient A *ᵥ x + c with A.PosDef. The paper's (4) is the case $c = -A\hat x$. Iterates, steps, directions and gradients are sequences indexed by ℕ. Division is Lean's total division. Every statement that divides therefore carries hypotheses making its denominators nonzero: $p_t \ne 0$ and $p_k\ne0$ in the quadratic statements, and $p_t'y_t\ne 0$ or $p'y>0$ in the generic ones.

The conjugacy $p_t'Ap_i=0$ for $t<i\le k$ is a hypothesis, exactly as the paper's proof uses it. It is not derived from a full run of Beale's algorithm. The range $k\le t+n-1$ of Beale's formula is not assumed.

Several formalizations would make the claim easier than the paper's, and none of them is used:

• defining the direction by the expanded formula (34), or the restart matrix by (38);
• adding orthogonality or conjugacy hypotheses beyond those listed;
• concluding only that the two directions are parallel;
• dropping the restart term of Beale's direction;
• allowing a zero denominator.

The generic milestones — (34), (40), (38), (42), (43) and the descent statements — are statements about arbitrary vectors and matrices and are reusable for any BFGS-based method. Proofs of any milestone, and alternative derivations of the goal, are welcome.

Selected references

• D. F. Shanno, Conjugate Gradient Methods with Inexact Searches, Mathematics of Operations Research 3(3) (1978) 244–256. https://doi.org/10.1287/moor.3.3.244
• E. M. L. Beale, A derivation of conjugate gradients, in F. A. Lootsma (ed.), Numerical Methods for Nonlinear Optimization, Academic Press, 1972, 39–43.
• M. J. D. Powell, Restart procedures for the conjugate gradient method, Mathematical Programming 12 (1977) 241–254. https://doi.org/10.1007/BF01593790
• M. R. Hestenes and E. Stiefel, Methods of conjugate gradients for solving linear systems, J. Res. Nat. Bur. Standards 49 (1952) 409–436. https://doi.org/10.6028/jres.049.044
• S. S. Oren and E. Spedicato, Optimal conditioning of self-scaling variable metric algorithms, Mathematical Programming 10 (1976) 70–90. https://doi.org/10.1007/BF01580654