Numerical Analysis

Motivation

Backward stochastic differential equations (BSDEs) describe the price and the hedge of a contingent claim in models with nonlinear pricing rules (differential interest rates, funding costs, reflected or constrained claims), and give probabilistic representations of semilinear parabolic PDEs (El Karoui, Peng and Quenez, 1997). Their numerical solution in moderate dimension is done by simulation: a backward recursion over a time grid in which each conditional expectation is replaced by a least-squares regression on simulated paths, the same device as the regression method for Bermudan options of Longstaff and Schwartz (2001).

Gobet, Lemor and Warin (2005) split the error of such a scheme into three parts: time discretization (Theorem 1), projection on finite function bases (Theorem 2), and the replacement of $\mathbf L_2$ projections by empirical regressions on $M$ simulated paths (Theorem 3). This mission formalizes the third part, which the authors describe as the major contribution of the paper. Its point is that the error from the simulations is controlled nonasymptotically, step by step, without blowing up as the time step $h$ shrinks, even though every regression of the backward recursion reuses the same simulated paths.

Setting

A model consists of a horizon $T>0$, a drift $b$, a diffusion $\sigma$ satisfying the Lipschitz condition (H1), and a driver $f(t,x,y,z)$ satisfying (H2): $|f(t_2,x_2,y_2,z_2)-f(t_1,x_1,y_1,z_1)|\le C_f(|t_2-t_1|^{1/2}+|x_2-x_1|+|y_2-y_1|+|z_2-z_1|)$. For $N\ge1$, $h=T/N$, $t_k=kh$, the Euler scheme is $S^N_{t_{k+1}}=S^N_{t_k}+b(t_k,S^N_{t_k})h+\sigma(t_k,S^N_{t_k})\Delta W_k$, where the increments $\Delta W_k\sim\mathcal N(0,hI_q)$ are independent of the past. A chain $P^N_{t_k}\in\mathbb R^{d'}$ extends $S^N_{t_k}$, and $\Phi^N(P^N_{t_N})$ is the terminal value.

At each time $t_k$, function bases $p_{0,k}$ (for $Y$) and $p_{l,k}$, $1\le l\le q$ (for the components of $Z$) are fixed, orthonormal in the sense $\mathbb E[p_{l,k}(P^N_{t_k})p_{l,k}(P^N_{t_k})^*]=\mathrm{Id}$. The projection–Picard scheme (Definition 1) computes coefficients $\alpha^{i,I}_k$ by $I$ Picard iterations of an $\mathbf L_2$ least-squares problem, and sets $Y^{N,I,I}_{t_k}=\alpha^{I,I}_{0,k}\cdot p_{0,k}$, $Z^{N,I,I}_{l,t_k}=\alpha^{I,I}_{l,k}\cdot p_{l,k}$.

The empirical scheme (4) replaces the expectation by an average over $M$ independent simulations $(P^{N,m},\Delta W^m)$ of the path: $\alpha^{i,I,M}_k$ minimizes

$$\frac1M\sum_{m=1}^M\Big(Y^{N,I,I,M,m}_{t_{k+1}}-\alpha_0\cdot p^m_{0,k}+hf^m_k(\alpha^{i-1,I,M}_k)-\sum_{l=1}^q\alpha_l\cdot p^m_{l,k}\Delta W^m_{l,k}\Big)^2 .$$

The outputs are truncated: with $\rho^N_{l,k}(x)=\max(1,C_0|p_{l,k}(x)|)$ and a smooth profile $\xi$ equal to the identity on $[-3/2,3/2]$, $Y^{N,I,I,M}_{t_k}=\hat\rho^N_{0,k}(\alpha^{I,I,M}_{0,k}\cdot p_{0,k})$ with $\hat\rho^N_{l,k}(x)=\rho^N_{l,k}(P^N_{t_k})\xi(x/\rho^N_{l,k}(P^N_{t_k}))$. The regression vector is $[v_k]^*=(p_{0,k}^*,p_{1,k}^*\Delta W_{1,k}/\sqrt h,\dots,p_{q,k}^*\Delta W_{q,k}/\sqrt h)$, and the good event $\mathbf A^M_k$ (27) asks that the empirical matrices $V^M_j=\frac1M\sum_mv^m_j[v^m_j]^*$ and $P^M_{l,j}=\frac1M\sum_m p^m_{l,j}[p^m_{l,j}]^*$ be close to the identity for all $j\ge k$.

Formalization targets

Goal: Theorem 3

For $I\ge3$, orthonormal bases with $\mathbb E|p_{l,k}|^4<\infty$, $C_0$ such that the bounds of Proposition 2 hold, and $h$ small enough, for $0\le k\le N-1$,

$$\mathbb E|Y^{N,I,I}_{t_k}-Y^{N,I,I,M}_{t_k}|^2+h\sum_{j=k}^{N-1}\mathbb E|Z^{N,I,I}_{t_j}-Z^{N,I,I,M}_{t_j}|^2\le 9\sum_{j=k}^{N-1}\mathbb E(|\rho^N_j|^2\mathbf 1_{[\mathbf A^M_k]^c})+Ch^{I-1}\sum_{j=k}^{N-1}[1+|S_0|^2+\mathbb E|\rho^N_j|^2]+\frac{C}{hM}\sum_{j=k}^{N-1}\epsilon_j,$$

where $\epsilon_j$ collects second moments of $v_jv_j^*-\mathrm{Id}$, $|v_j|^2|p_{0,j+1}|^2$ and $|v_j|^2(1+|S^N_{t_j}|^2+\dots)$, written out in full in the goal statement. The constant $C$ and the threshold on $h$ depend only on the model and on $\xi$.

Milestones

1. Proposition 2: the a priori bounds $|Y^{N,i,I}_{t_k}|\le\rho^N_{0,k}$, $\sqrt h|Z^{N,i,I}_{l,t_k}|\le\rho^N_{l,k}$ that fix the truncation levels.
2. (28)–(29): the empirical least-squares solution and its contraction inequality $\lambda_{\min}(V^M)|\theta_x|^2\le|\theta_x\cdot v|^2_M\le|x|^2_M$.
3. Lemma 1: on $\mathbf A^M_k$ the empirical Picard iterations contract at rate $Ch$ to a unique fixed point, with error $[Ch]^I$ after $I$ steps.
4. (32): the pathwise bound $|\theta^{i,I,M}_k|^2\le C(\mathcal A^{N,M}_{k+1}+h\mathcal B^{N,M}_k)$ on $\mathbf A^M_k$.
5. (34): the expectation formula $\theta^{\infty,I}_k=\mathbb E(v_k[Y^{N,I,I}_{t_{k+1}}+hf_k(\alpha^{\infty,I}_k)])$.

Significance

Theorem 3 is nonasymptotic: together with Theorems 1 and 2 it lets one compare the three error sources and choose $h$, the bases and $M$ jointly for a target accuracy. The $1/(hM)$ rate shows how many paths a finer time grid requires, and the $h^{I-1}$ term shows that $I=3$ Picard iterations suffice. The term involving $[\mathbf A^M_k]^c$ isolates the event on which the empirical regression matrices are badly conditioned, which the truncation keeps under control.

The result is proved in the paper; nothing in this mission is open mathematically. No machine-checked proof of any part of it is known. A formal proof would check a long chain of estimates whose constants the paper tracks only as a generic $C$, and would settle the two misprints in the printed statement (see below). The definitions layer (Euler scheme, regression schemes, empirical regression matrices) is reusable for other regression Monte Carlo schemes for BSDEs and for optimal stopping.

Difficulty

The obvious approach is to treat each regression as an independent statistical estimation problem and apply a variance bound per time step. This fails because all regressions of the backward recursion use the same $M$ paths: the response at time $t_{k+1}$ is itself a function of the simulations, so the regression at $t_k$ is not a regression of a fixed variable on independent samples. Moreover, the empirical matrix $V^M_k$ may be singular, and on that event the empirical coefficients are unbounded unless truncated. A naive per-step bound also produces a factor $1/h$ at each of the $N=T/h$ steps, which explodes; the proof must keep the accumulated constants of order one.

Formalization scope

• All objects are defined in the namespace RegMCBSDE.Simulation. Continuous time is not used: the increments $\Delta W_k$ are any family that is $\mathcal F_{k+1}$-measurable, independent of $\mathcal F_k$ and $\mathcal N(0,hI_q)$-distributed for some filtration. That covers the Brownian case. (H3), a condition on the terminal functional of the continuous path, is dropped. So is the Markov representation of $P^N$, which no statement uses.
• Expectations of squared quantities are taken in $[0,\infty]$, on both sides of every inequality.
• Euclidean norms are written as sums of squares. $\|A\|\le c$ for symmetric $A$ is written as $|x^*Ax|\le c|x|^2$ for all $x$.
• The schemes are defined as predicates (any minimizer of each least-squares problem), never through a matrix inverse. The empirical coefficients are required to be measurable functions of the simulations, and any such minimizer is allowed, as the paper says the choice is arbitrary.
• The reference path at which the fitted coefficients are evaluated is assumed independent of the $M$ simulations.
• Picard iterations are imposed for every $i\ge1$.
• "For $h$ small enough" is formalized as $T/N<h_0$. Every constant $C$ and every $h_0$ is chosen before $N$, $I$, $M$, $k$, $C_0$, the bases, $S_0$ and the probability space. A formalization in which $C$ is chosen after the scheme data would make every inequality with a positive right-hand side trivially true, and is ruled out.
• "$C_0$ large enough" in Theorem 3 is the hypothesis that the bounds of Proposition 2 hold for $C_0$.
• Proposition 2 is stated for $h$ small enough. The printed statement omits this condition, which its proof needs through the uniform bound (19).
• In (32) and Lemma 1, $\rho^N_{0,N}$, which the paper leaves undefined, is read through the terminal response $\Phi^N(P^{N,m}_{t_N})$.
• Corrections to Theorem 3 (both from the paper's proof, p. 21):
  1. The printed $j=N-1$ summand contains the undefined $p_{0,N}$. It is replaced by $\mathbb E(|v_{N-1}|^2|\Phi^N(P^N_{t_N})|^2)$.
  2. The printed factor $\mathbb E|\rho^N_{0,j}(P^N_{t_j})|^2$ next to $\mathbb E(|v_j|^2|p_{0,j+1}|^2)$ becomes $\mathbb E|\rho^N_{0,j+1}(P^N_{t_{j+1}})|^2$.

Contributions welcome: proofs of the milestones, in particular (28)–(29) (finite-dimensional linear algebra) and (34) (independence and orthonormality), and a construction showing that measurable minimizers of (4) exist.

Selected references

• E. Gobet, J.-P. Lemor, X. Warin, A regression-based Monte Carlo method to solve backward stochastic differential equations, Ann. Appl. Probab. 15(3), 2172–2202, 2005. https://doi.org/10.1214/105051605000000412, preprint https://arxiv.org/abs/math/0508491v1
• N. El Karoui, S. Peng, M. C. Quenez, Backward stochastic differential equations in finance, Math. Finance 7(1), 1–71, 1997. https://doi.org/10.1111/1467-9965.00022
• F. A. Longstaff, E. S. Schwartz, Valuing American options by simulation: a simple least-squares approach, Rev. Financ. Stud. 14(1), 113–147, 2001. https://doi.org/10.1093/rfs/14.1.113
• J. Zhang, A numerical scheme for BSDEs, Ann. Appl. Probab. 14(1), 459–488, 2004. https://doi.org/10.1214/aoap/1075828058
• B. Bouchard, N. Touzi, Discrete-time approximation and Monte-Carlo simulation of backward stochastic differential equations, Stochastic Process. Appl. 111(2), 175–206, 2004. https://doi.org/10.1016/j.spa.2004.01.001

Motivation

Many estimation problems ask for a low-rank matrix from far fewer linear measurements than it has entries: matrix completion, quantum state tomography, recovery of positions from partial distances. The rank minimization problem $\min \operatorname{rank}(X)$ subject to $\mathcal S(X) = \mathcal M$ is NP-hard, so one solves a tractable surrogate instead, either the nuclear-norm minimization $\min \|X\|_*$ subject to $\mathcal S(X) = \mathcal M$ or an iterative scheme such as the iteratively reweighted least squares algorithm IRLS-M of Fornasier, Rauhut and Ward (arXiv:1010.2471).

Guarantees for these surrogates rest on conditions on the measurement map. Two are standard. The restricted isometry property (RIP) asks that $\mathcal S$ nearly preserve the Frobenius norm of every low-rank matrix; it holds with high probability for many random maps, and it is what one can verify in practice. The null space property asks that no nonzero matrix in the kernel of $\mathcal S$ be concentrated on few singular directions; it is what the recovery proofs actually use. This mission formalizes the bridge between them in the paper's form: a small restricted isometry constant implies the strong rank null space property (SRNSP) that drives the convergence analysis of IRLS-M.

Timeline.

• 2005: Candès and Tao introduce restricted isometry constants for sparse vectors and prove exact recovery by $\ell_1$ minimization under a condition on them (IEEE Trans. Inf. Theory).
• 2008: Candès shows that $\delta_{2s} < \sqrt 2 - 1$ suffices in the vector case (C. R. Math.).
• 2010: Recht, Fazel and Parrilo define the rank restricted isometry property and prove that nuclear-norm minimization recovers every rank-$r$ matrix when $\delta_{5r}$ is small enough (SIAM Review).
• 2011: Candès and Plan prove the matrix analogue of the $\sqrt 2 - 1$ threshold and the orthogonality estimate quoted here as Lemma 6.9 (IEEE Trans. Inf. Theory).
• 2011: Fornasier, Rauhut and Ward prove Proposition 6.8, the target of this mission: $\delta_{4k} < \sqrt 2 - 1$ implies the SRNSP of order $k$.

Setting

Fix integers $n, p, m$. A measurement map is a linear map $\mathcal S : \mathbb R^{n\times p} \to \mathbb R^m$; it is given by $m$ matrices $A_1,\dots,A_m$ through $\mathcal S(X)_\ell = \langle A_\ell, X\rangle$, where $\langle X, Y\rangle = \operatorname{Tr}(XY^{\top}) = \sum_{i,j} X_{ij}Y_{ij}$. The Frobenius norm is $\|X\|_F = \langle X, X\rangle^{1/2}$, the nuclear norm $\|X\|_*$ is the sum of the singular values of $X$, and $\|\mathcal S(X)\|_{\ell_2^m}$ is the Euclidean norm of the vector $\mathcal S(X)$. A matrix is $k$-rank if its rank is at most $k$.

The restricted isometry constant $\delta_k = \delta_k(\mathcal S)$ (Definition 1.1) is the smallest $\delta \ge 0$ such that

$$(1 - \delta)\|X\|_F^2 \le \|\mathcal S(X)\|_{\ell_2^m}^2 \le (1+\delta)\|X\|_F^2 \qquad \text{for all $k$-rank } X .$$

The map $\mathcal S$ has the strong rank null space property of order $k$ with constant $\eta \in (0,1)$ (Definition 6.4) if for every $X \in \ker \mathcal S \setminus\{0\}$ and every decomposition $X = X_1 + X_2$ with $X_1$ $k$-rank, there is another decomposition $X = H_1 + H_2$ with $\operatorname{rank} H_1 \le 2k$, $\langle H_1, H_2\rangle = 0$, $X_1H_2^{\top} = 0$, $X_1^{\top}H_2 = 0$, and

$$\|H_1\|_* \le \eta\,\|H_2\|_* .$$

In Lean these are IRLSM.RIP.ripConst A k, IRLSM.RIP.SRNSP A k η and IRLSM.RIP.measSq A X $= \|\mathcal S(X)\|^2_{\ell_2^m}$, built on the published module HighDimStat_MatrixRank_Core (traceInner, frobeniusNorm, nuclearNorm, observationOp).

Formalization targets

Goal: Proposition 6.8

If $\delta_{4k} > 0$ and

$$\delta_{4k} < \sqrt 2 - 1,$$

then $\mathcal S$ satisfies the SRNSP of order $k$ with constant

$$\eta = \sqrt 2\,\frac{\delta_{4k}}{1 - \delta_{3k}} \in (0,1).$$

The constant $\eta$ is the paper's explicit one; the goal asserts the full property for every kernel element and every decomposition, including $\eta \in (0,1)$.

Milestones

1. Lemma 6.9 (Candès–Plan): for $\langle X, Y\rangle = 0$ and $\operatorname{rank}X + \operatorname{rank}Y \le k$, $|\langle \mathcal S(X), \mathcal S(Y)\rangle| \le \delta_k \|X\|_F\|Y\|_F$.
2. Lemma 6.10 (Recht–Fazel–Parrilo): $XZ^{\top} = 0$ and $X^{\top}Z = 0$ imply $\|X + Z\|_* = \|X\|_* + \|Z\|_*$.
3. The block decomposition of the proof: relative to an SVD $X_1 = U(\begin{smallmatrix}\Sigma & 0\\ 0 & 0\end{smallmatrix})V^{\top}$, the matrices $H_0$ and $H_c$ obtained by zeroing blocks of $U^{\top}XV$ satisfy $X = H_0 + H_c$, $\operatorname{rank}H_0 \le 2k$, $X_1H_c^{\top} = 0$, $X_1^{\top}H_c = 0$, $\langle H_0, H_c\rangle = 0$.
4. $\|H\|_* \le \sqrt r\,\|H\|_F$ for $\operatorname{rank}H \le r$, and $\|H\|_F \le \|H + K\|_F$ for $\langle H, K\rangle = 0$.
5. The shelling bound: for a nonincreasing nonnegative sequence cut into blocks of length $\ell$, each block's $\ell_2$ norm is at most $\ell^{-1/2}$ times the $\ell_1$ norm of the previous block.
6. Monotonicity $\delta_k \le \delta_{k'}$ for $k \le k'$, and $\eta \in (0,1)$ under $0 < \delta_{4k} < \sqrt2 - 1$.

Significance

The result. Proposition 6.8 is the step that turns a verifiable, probabilistically typical hypothesis (a small rank restricted isometry constant, which Gaussian and many other random maps satisfy once $m$ is of order $k\max(n,p)$) into the deterministic null-space hypothesis under which the paper proves that IRLS-M converges to the nuclear-norm minimizer and recovers every $k$-rank matrix exactly (Theorem 6.11, Proposition 2.1). Through Lemma 6.6 and Corollary 6.7 of the same paper, the SRNSP also gives stable recovery by nuclear-norm minimization, with error controlled by the best $k$-rank approximation error in the nuclear norm.

Formalizing it. The result is proved on paper; it has no machine-checked proof that we know of. A complete development supplies the rank restricted isometry constant as a usable object (its monotonicity, its polarization estimate Lemma 6.9), the nuclear norm's additivity on orthogonal pieces, its comparison with the Frobenius norm, and the block-splitting of the singular value decomposition. These tools are reusable across low-rank recovery. The sibling mission of this series, IRLS-M Converges to the Nuclear-Norm Minimizer Under the Strong Rank Null Space Property, takes the SRNSP as its hypothesis; together the two give Proposition 2.1 of the paper.

Difficulty

The obvious argument would bound $\|H_1\|_*$ directly by the restricted isometry inequality, but that inequality controls only Frobenius norms of low-rank matrices, while the kernel element $X$ has full rank in general and the property is stated in nuclear norms. The restricted isometry hypothesis therefore cannot be applied to $X$, or to $H_2$, as a whole, and the constant $\eta$ depends on how the high-rank part is handled. In Lean, the singular value decomposition of a rectangular matrix and the behaviour of singular values under block operations are the main missing infrastructure: Mathlib has the spectral theorem for Hermitian matrices but little on singular values of rectangular matrices.

Formalization scope

• Real matrices Matrix (Fin n) (Fin p) ℝ; the paper treats real and complex matrices indifferently. $X^*$ is $X^{\top}$ and $\langle\cdot,\cdot\rangle$ is traceInner.
• Measurement map by matrices: $\mathcal S(X)_\ell = \langle A_\ell, X\rangle$ (observationOp); $\|\mathcal S(X)\|^2_{\ell_2^m}$ is a sum of squares, not Mathlib's sup norm on Fin m → ℝ.
• The squared RIP constant of Definition 1.1, defined as the infimum of the admissible $\delta \ge 0$ (attained, never a junk value), for every $k \in \mathbb N$; the page's range $1 \le k \le n$ is dropped because the goal uses $\delta_{4k}$ with $4k$ possibly larger than $n$. The unsquared constant $(1-\delta)\|X\|_F \le \|\mathcal S(X)\| \le (1+\delta)\|X\|_F$ used by other formalizations is a different quantity and is not used.
• Positivity $\delta_{4k} > 0$ is a hypothesis of the goal: it is Definition 1.1's own "$\delta_k > 0$", and without it $\eta = 0 \notin (0,1)$.
• No n ≤ p assumption: the proof does not use it.
• The block decomposition milestone states the page's explicit construction, with a general (not necessarily diagonal) top-left block, and without the unused hypothesis $X \in \ker\mathcal S$.
• The shelling bound is stated in its vector form for every block; the page writes "for $j \ge 2$" but display (6.7) uses it for every $j \ge 1$.

A formalization that proved only $\eta < 1$, or only the existence of some decomposition without the bound $\|H_1\|_* \le \eta\|H_2\|_*$, or that defined $\delta_k$ as any admissible constant rather than the smallest one, would not be this mission's goal.

Welcome contributions: singular value decomposition for rectangular real matrices in the Core encoding, singular values of block-diagonal and orthogonally-equivalent matrices, and proofs of the milestones in any order.

Selected references

• M. Fornasier, H. Rauhut, R. Ward, Low-rank matrix recovery via iteratively reweighted least squares minimization, SIAM J. Optim. 21(4), 2011; read as arXiv:1010.2471v4. https://arxiv.org/abs/1010.2471
• E. J. Candès, Y. Plan, Tight oracle inequalities for low-rank matrix recovery from a minimal number of noisy random measurements, IEEE Trans. Inf. Theory 57(4), 2011. https://doi.org/10.1109/TIT.2011.2111771
• B. Recht, M. Fazel, P. A. Parrilo, Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization, SIAM Review 52(3), 2010. https://doi.org/10.1137/070697835
• E. J. Candès, T. Tao, Decoding by linear programming, IEEE Trans. Inf. Theory 51(12), 2005. https://doi.org/10.1109/TIT.2005.858979
• E. J. Candès, The restricted isometry property and its implications for compressed sensing, C. R. Math. Acad. Sci. Paris 346, 2008. https://doi.org/10.1016/j.crma.2008.03.014

Motivation

Many data problems ask for a matrix of low rank from far fewer linear measurements than it has entries: matrix completion (recommender systems, where a few ratings of a large user–item table are observed), system identification, and quantum state tomography. Minimizing the rank under the measurement constraints is NP-hard in general. The standard convex surrogate replaces the rank by the nuclear norm, the sum of the singular values, and recovers the low-rank matrix exactly when the measurement map is well conditioned on low-rank matrices (Recht, Fazel, Parrilo 2010; Candès, Recht 2009).

Solving the nuclear norm problem at scale is itself a numerical challenge: interior point methods do not scale, and first-order methods such as singular value thresholding (Cai, Candès, Shen 2010) need many iterations. Fornasier, Rauhut and Ward proposed an alternative, iteratively reweighted least squares for matrices (IRLS-M), which replaces the nonsmooth problem by a sequence of weighted least squares problems whose weights are recomputed from the current iterate. It is the matrix analogue of the IRLS method for sparse vectors of Daubechies, DeVore, Fornasier, Güntürk 2010, and a closely related algorithm was studied at the same time by Mohan and Fazel (Allerton 2010; journal version JMLR 2012).

Timeline. 2008–2010: nuclear norm recovery under rank restricted isometry (Recht–Fazel–Parrilo) and the rank null space property (Recht–Xu–Hassibi). 2010: convergence of IRLS for $\ell_1$ under the vector null space property (Daubechies et al.). 2010–2011: IRLS-M and its convergence under the strong rank null space property (this paper, arXiv:1010.2471, SIAM J. Optim. 21(4), 2011).

Setting

Throughout, $X$ is a real $n\times p$ matrix with $n\le p$, with singular values $\sigma_1(X)\ge\dots\ge\sigma_n(X)\ge0$. The Frobenius norm is $\|X\|_F$, the nuclear norm is $\|X\|_*=\sum_i\sigma_i(X)$, and $\langle X,Y\rangle=\sum_{i,j}X_{ij}Y_{ij}$. A linear measurement map $\mathcal S:\mathbb R^{n\times p}\to\mathbb R^m$ is given by matrices $A_1,\dots,A_m$ through $\mathcal S(X)_l=\langle A_l,X\rangle$; the data are $\mathscr M\in\mathbb R^m$.

The best $k$-rank approximation error is $\rho_k(X)_*=\min_{\operatorname{rank}Z\le k}\|X-Z\|_*$. The map $\mathcal S$ has the strong rank null space property (SRNSP) of order $k$ with constant $\eta\in(0,1)$ if every nonzero $X\in\ker\mathcal S$ and every split $X=X_1+X_2$ with $\operatorname{rank}X_1\le k$ admit another split $X=H_1+H_2$ with $\operatorname{rank}H_1\le2k$, $\langle H_1,H_2\rangle=0$, $X_1H_2^{\mathsf T}=0$, $X_1^{\mathsf T}H_2=0$ and $\|H_1\|_*\le\eta\|H_2\|_*$.

For $\varepsilon>0$ the weight of $X$ is $W=U\Sigma_\varepsilon^{-1}U^{\mathsf T}$, where $XX^{\mathsf T}=U\Sigma^2U^{\mathsf T}$ and $\Sigma_\varepsilon=\operatorname{diag}(\max\{\sigma_j,\varepsilon\})$. The IRLS-M algorithm with parameters $K\in\mathbb N$ and $\gamma>0$ starts from $W^0=I$, $\varepsilon_0=1$ and repeats

$$X^\ell\in\arg\min_{\mathcal S(X)=\mathscr M}\|(W^{\ell-1})^{1/2}X\|_F^2,\qquad\varepsilon_\ell=\min\{\varepsilon_{\ell-1},\gamma\sigma_{K+1}(X^\ell)\},$$

with $W^\ell$ the weight of $X^\ell$ at level $\varepsilon_\ell$; it stops when $\varepsilon_\ell=0$. Two functionals organize the analysis: $\mathcal J(X,W)=\tfrac12(\|W^{1/2}X\|_F^2+\|W^{-1/2}\|_F^2)$, of which IRLS-M is an alternating minimization, and $\mathcal J_\varepsilon(X)=\sum_{i=1}^n j^\varepsilon(\sigma_i(X))$ with $j^\varepsilon(u)=|u|$ for $|u|\ge\varepsilon$ and $(u^2+\varepsilon^2)/(2\varepsilon)$ otherwise.

Formalization targets

Goal: Theorem 6.11

For $\gamma=1/n$, any $K$, a surjective $\mathcal S$ and any run $(X^\ell,\varepsilon_\ell)$:

$$\text{(i)}\quad \mathcal S\ \text{SRNSP of order }K,\ \varepsilon_\ell\to0\ \Longrightarrow\ X^\ell\to\bar X,\ \operatorname{rank}\bar X\le K,\ \bar X=\text{the unique nuclear norm minimizer.}$$

(ii) If $\varepsilon_\ell\to\varepsilon>0$, the iterates are relatively compact, their accumulation points minimize $\mathcal J_\varepsilon$ on the feasible set, and the sequence converges when that minimizer is unique; under the SRNSP with $\eta<1-2/(K-2)$, every accumulation point $\bar X$ satisfies, for every feasible $X$ and $k<K-2\eta/(1-\eta)$,

$$\|X-\bar X\|_*\le\Lambda\rho_k(X)_*,\qquad\Lambda=\frac{4(1+\eta)^2}{(1-\eta)^2((K-k)(1-\eta)-2\eta)}+\frac{2(1+\eta)}{1-\eta}.$$

(iii) Under the same SRNSP, if a feasible matrix of rank at most $k$ exists, then $\varepsilon_\ell\to0$.

Milestones

In attack order: the weighted least squares step (Lemma 5.1, the optimality condition (5.4)); the weight step (Lemma 5.2, Proposition 5.3); monotonicity of $\mathcal J$ along the run (Proposition 6.1); singular value facts (Weyl's Theorem 7.1, attainment of $\rho_k$ at the spectral truncation, Lemma 6.10, Proposition 7.2); null space theory (Theorem 6.3, the inverse triangle inequality Lemma 6.6 with (6.5), Corollary 6.7); and the two steps of the proof of (ii), the optimality condition (6.10) for $\mathcal J_\varepsilon$ and the bound (6.11).

Significance

Theorem 6.11 is the convergence guarantee of IRLS-M. Combined with the fact that a small rank restricted isometry constant implies the SRNSP (Proposition 6.8, the companion mission), it yields Proposition 2.1 of the paper: under $\delta_{4K}<\sqrt2-1$, IRLS-M recovers every rank-$k$ matrix exactly and approximately low-rank matrices stably. Part (iii) says that the algorithm detects exact low-rank data on its own: the smoothing parameter must vanish.

The result is proved on paper; no machine-checked proof exists. A complete formalization would give a machine-checked convergence proof of an IRLS-type method; none is on the platform, for vectors or matrices. It also produces reusable matrix analysis: Weyl's perturbation bound for singular values, the nuclear norm best rank-$k$ approximation (a Mirsky-type theorem), additivity of the nuclear norm under orthogonality, the null space characterization of nuclear norm recovery, and optimality conditions for spectral functions.

Difficulty

Two steps resist the obvious argument. First, the weight update is a constrained minimization over the positive definite cone; identifying its solution (Proposition 5.3) requires a duality or spectral argument rather than setting a gradient to zero, because the constraint $W\preceq\varepsilon^{-1}I$ is active on small singular values. Second, in case (ii) the limit functional $\mathcal J_\varepsilon$ is convex but not strictly convex, so accumulation points cannot be pinned down by uniqueness; the analysis must instead pass the optimality condition (5.4) of the iterates to the limit, using continuity of the weights in the iterate, and derive (6.10). The step relating $\nabla\mathcal J_\varepsilon$ to $W\bar X$ needs the differentiability of unitarily invariant spectral functions, which the paper imports from Lewis and Sendov. Weyl's inequality and the spectral truncation facts, routine on paper, are not yet in Mathlib for rectangular matrices.

Formalization scope

All matrices are real (Matrix (Fin n) (Fin p) ℝ); the paper treats real and complex matrices alike. The paper's standing assumption $n\le p$ is a hypothesis wherever $n\times n$ objects ($XX^{\mathsf T}$, $W$, $\gamma=1/n$, $\mathcal J_\varepsilon$) occur. Nuclear norm, Frobenius norm, trace inner product and singular values come from the published definition HighDimStat_MatrixRank_Core; sv X i is $\sigma_{i+1}(X)$, 0-based. The measurement map is given by measurement matrices (observationOp), which represents every linear map to $\mathbb R^m$. The weight is defined with Mathlib's continuous functional calculus, and $W^{1/2}$ is cfc Real.sqrt W. The algorithm is the relation IsRun: each $X^{\ell+1}$ is some minimizer of (2.9), $X^0$ is free, and $\varepsilon_j=0$ after a stop. Accumulation points are MapClusterPt; uniqueness of a nuclear norm minimizer is a strict inequality against every other feasible matrix.

The goal cannot be satisfied vacuously: IsRun has a run for every surjective $\mathcal S$ and every $\mathscr M$ (Lemma 5.1 supplies the minimizers), it carries no positivity or uniqueness assumption a real run could violate, and the theorem keeps all three parts, including the unique-minimizer clause of (i) and the error bound of (ii) for every feasible $X$ and every admissible $k$.

Two printed slips are corrected and disclosed: (5.5) omits the factor $\tfrac12$ of (5.1), and (6.5) reverses the bracket of (6.4) at $\rho_k=0$. Contributions welcome: Weyl's inequality and the Mirsky-type truncation theorem for rectangular matrices, the spectral calculus behind Propositions 5.3 and (6.10), and proofs of the null space results, which are independent of the algorithm.

Selected references

• M. Fornasier, H. Rauhut, R. Ward, Low-rank matrix recovery via iteratively reweighted least squares minimization, SIAM J. Optim. 21(4), 2011. https://arxiv.org/abs/1010.2471 (v4 is the version formalized)
• B. Recht, M. Fazel, P. A. Parrilo, Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization, SIAM Review 52(3), 2010. https://doi.org/10.1137/070697835
• B. Recht, W. Xu, B. Hassibi, Null space conditions and thresholds for rank minimization, Math. Program. 127, 2011. https://doi.org/10.1007/s10107-010-0422-2
• I. Daubechies, R. DeVore, M. Fornasier, C. S. Güntürk, Iteratively reweighted least squares minimization for sparse recovery, Comm. Pure Appl. Math. 63(1), 2010. https://doi.org/10.1002/cpa.20303
• K. Mohan, M. Fazel, Iterative reweighted least squares for matrix rank minimization, Proc. Allerton Conference, 2010; journal version Iterative reweighted algorithms for matrix rank minimization, JMLR 13, 2012. https://jmlr.org/papers/v13/mohan12a.html
• E. J. Candès, B. Recht, Exact matrix completion via convex optimization, Found. Comput. Math. 9, 2009. https://doi.org/10.1007/s10208-009-9045-5
• J.-F. Cai, E. J. Candès, Z. Shen, A singular value thresholding algorithm for matrix completion, SIAM J. Optim. 20(4), 2010. https://doi.org/10.1137/080738970
• A. S. Lewis, H. S. Sendov, Nonsmooth analysis of singular values. I: Theory, Set-Valued Analysis 13(3), 2005. https://doi.org/10.1007/s11228-004-7197-7

Motivation

Many objective functions in engineering and statistics are nonsmooth and nonconvex: the spectral abscissa of a parametrized matrix, the largest eigenvalue of a symmetric matrix function, the $H_\infty$ norm of a closed-loop system, or a maximum of finitely many smooth functions. Such functions are typically differentiable almost everywhere, yet their minimizers usually sit exactly where the gradient fails to exist. Gradient descent then zig-zags or stalls, and bundle methods, which are designed for convex functions, lose their guarantees.

Burke, Lewis and Overton (SIAM J. Optim. 15 (2005) 751–779) proposed the gradient sampling (GS) algorithm: at each iterate, sample gradients at random nearby points, take the element of least norm in their convex hull, and search along its negative. The method needs nothing beyond gradients at points of differentiability, and it has since become a standard tool for nonsmooth, nonconvex minimization (Burke, Curtis, Lewis, Overton, Simões, 2020). This mission formalizes the paper's convergence theory.

Timeline. Goldstein (1977) introduced the $\epsilon$-subdifferential of a Lipschitz function and a conceptual descent method built on it. Clarke's generalized gradient (Clarke 1983) supplied the stationarity notion. Burke, Lewis and Overton (2005) made Goldstein's idea implementable by random sampling and proved almost-sure convergence to Clarke $\epsilon$-stationary points for a fixed sampling radius. Kiwiel (2007) revisited the analysis and proved convergence for a modified version of the algorithm.

Setting

Let $f:\mathbb R^n\to\mathbb R$ be locally Lipschitz, and let $D\subseteq\mathbb R^n$ be an open dense set on which $f$ is continuously differentiable. Write $\mathbb B$ for the closed Euclidean unit ball, and fix $\tilde x$ such that the level set $\mathcal L=\{x: f(x)\le f(\tilde x)\}$ is compact.

The Clarke subdifferential $\bar\partial f(x)$ is the convex hull of all limits $\lim_i\nabla f(x+h_i)$ with $h_i\to0$ and $f$ differentiable at each $x+h_i$. For $\epsilon>0$ define

$$G_\epsilon(x)=\operatorname{cl}\operatorname{conv}\nabla f\big((x+\epsilon\mathbb B)\cap D\big),\qquad \rho_\epsilon(x)=\operatorname{dist}\big(0\mid G_\epsilon(x)\big),$$

and the Clarke $\epsilon$-subdifferential $\bar\partial_\epsilon f(x)=\operatorname{cl}\operatorname{conv}\bigcup_{\|y-x\|\le\epsilon}\bar\partial f(y)$. A point is Clarke $\epsilon$-stationary if $0\in\bar\partial_\epsilon f(x)$, and Clarke stationary if $0\in\bar\partial f(x)$.

The GS algorithm. Fix $x^0\in\mathcal L\cap D$, $\gamma,\beta\in(0,1)$, $\epsilon_0>0$, $\nu_0\ge0$, $\mu\in(0,1]$, $\theta\in(0,1]$ and a sample size $m\ge n+1$. At iteration $k$:

1. Draw $u^{k1},\dots,u^{km}$ independently and uniformly from $\mathbb B$ and set $x^{kj}=x^k+\epsilon_k u^{kj}$. If some $x^{kj}\notin D$, stop. Otherwise let $G_k=\operatorname{conv}\{\nabla f(x^k),\nabla f(x^{k1}),\dots,\nabla f(x^{km})\}$.
2. Let $g^k$ be the least-norm element of $G_k$. If $\nu_k=\|g^k\|=0$, stop. If $\|g^k\|\le\nu_k$, set $t_k=0$, $\nu_{k+1}=\theta\nu_k$, $\epsilon_{k+1}=\mu\epsilon_k$ and go to step 4. Otherwise keep $\nu_{k+1}=\nu_k$, $\epsilon_{k+1}=\epsilon_k$ and set $d^k=-g^k/\|g^k\|$.
3. Let $t_k$ be the largest $\gamma^s$, $s\in\{0,1,2,\dots\}$, with $f(x^k+\gamma^s d^k)<f(x^k)-\beta\gamma^s\|g^k\|$.
4. If $x^k+t_kd^k\in D$, set $x^{k+1}=x^k+t_kd^k$. Otherwise pick any $\hat x^k\in x^k+\epsilon_k\mathbb B$ with $\hat x^k+t_kd^k\in D$ and $f(\hat x^k+t_kd^k)<f(x^k)-\beta t_k\|g^k\|$, and set $x^{k+1}=\hat x^k+t_kd^k$.

With $\nu_0=0$ and $\mu=1$ the radius stays fixed, $\epsilon_k=\epsilon_0=\epsilon$: this is fixed-radius mode.

Formalization targets

Goal: Theorem 3.4 (p. 760)

In fixed-radius mode, with probability $1$, either the algorithm stops at some iteration $k_0$ with $\rho_\epsilon(x^{k_0})=0$, or it runs forever and there is a subsequence $J$ with

$$\rho_\epsilon(x^k)\xrightarrow[k\in J]{}0,\qquad 0\in\bar\partial_\epsilon f(\bar x)\ \text{ for every cluster point }\bar x\text{ of }\{x^k\}_{k\in J}.$$

The theorem asserts nothing about $\|g^k\|$ along $J$; the authors list $\|g^k\|\to_J0$ as an open question.

Milestones

The milestones follow the paper's argument: the minimax characterization of the search direction (Lemma 2.1), the line-search descent estimate and the descent property of the iterates (§2), the almost-sure absence of a Step 1 stop, upper semicontinuity of $\rho_\epsilon$ and the uniform covering of Lemma 3.2(iv)–(v), Lebourg's mean value theorem (Theorem 3.3), the representation $\bar\partial f(x)=\bigcap_{\epsilon>0}G_\epsilon(x)$, the estimate of Lemma 3.1, the inclusion $G_\epsilon(x)\subseteq\bar\partial_\epsilon f(x)$, and the closed graph of $\bar\partial_\epsilon f$.

Companions

Well-posedness of the algorithm (every sample realization admits a run), Corollary 3.5(1) (the $C^1$ case: $\|g^k\|\to0$), Corollary 3.6 (Clarke $\epsilon_j$-stationary points with $\epsilon_j\downarrow0$ cluster at Clarke stationary points, and with a unique stationary point in $\mathcal L$ the sets $C_j$ converge to it in the Hausdorff sense), Corollary 3.7 (with a unique stationary point in $\mathcal L$, small $\epsilon$ forces the run near it), and Theorem 3.8 (the variable-radius mode with $\nu_0>0$).

Significance

Theorem 3.4 is the first convergence guarantee for a practical method on general locally Lipschitz, nonconvex functions that uses only gradients at points of differentiability. It holds for a fixed sampling radius, the setting where most of the work is, and Corollaries 3.6–3.7 and Theorem 3.8 transfer it to radii decreasing to zero, which is how the method is used in practice. Later gradient sampling variants and their analyses start from these results.

The results are proved on paper; no machine-checked proof of them is known. A formalization would make the probabilistic structure precise, which the paper leaves informal (see Formalization scope), and would certify the hypotheses under which the theorem holds. Already the statements record two such corrections.

Difficulty

The natural argument, "$f(x^k)$ decreases and $\mathcal L$ is compact, so the steps $t_k\|g^k\|$ tend to zero, hence $\|g^k\|\to0$", fails: $t_k$ can go to zero while $\|g^k\|$ stays bounded away from zero. The difficulty is to show that, on the event $\inf_k\rho_\epsilon(x^k)>0$, the random samples land infinitely often in a configuration whose sampled hull nearly realizes $\rho_\epsilon$ at a fixed nearby point, uniformly over a compact set of iterates, and that the failed trial step at such an iteration contradicts the Armijo rule. The iterates depend on all past samples, so independence across iterations cannot be used directly.

Formalization scope

• $\mathbb R^n$ is EuclideanSpace ℝ (Fin n), $\nabla f$ is Mathlib's gradient, the Clarke subdifferential is the published ClarkeGradients.Shared.generalizedGradient, and $\operatorname{dist}(0\mid C)$ is Metric.infDist 0 C. Local Lipschitz continuity is Mathlib's LocallyLipschitz; continuous differentiability on $D$ is ContDiffOn ℝ 1 f D.
• A run is a predicate IsGSRun on the iterates, radii, tolerances, step lengths, least-norm vectors, directions and a stopping index $\tau\in\mathbb N\cup\{\infty\}$, imposing Steps 0–4 up to $\tau$. "$\max\gamma^s$" is the least admissible $s$, and the Armijo and (2) inequalities are strict.
• A random run (IsRandomGSRun) lives on a probability space with a filtration $(\mathcal F_k)$: the samples are uniform on $\mathbb B$, independent within an iteration, and independent of $\mathcal F_k$, while $x^k,\epsilon_k,\nu_k$ are $\mathcal F_k$-measurable. The Step 4 choice may use the past and extra randomness but never future samples. "With probability 1" is ∀ᵐ. The paper's event $\mathcal E$, the realizations that hit every positive-measure subset of $\mathbb B^m$ infinitely often, is not used, because no sequence has that property.
• Added hypothesis $\operatorname{vol}(D^c)=0$. An open dense set can have a complement of positive measure, and then Step 1 stops with positive probability, so Theorem 3.4 would fail. Every probabilistic statement, and the two inclusions $\bar\partial f(x)\subseteq G_\epsilon(x)$, $\bar\partial_{\epsilon_1}f(x)\subseteq G_{\epsilon_2}(x)$ (which fail without it), carry this hypothesis.
• Trivializing encodings are ruled out. The run predicate is satisfiable for every sample realization (a companion statement), so the goal is not vacuous. The Step 4 choice cannot anticipate future samples. The conclusion is almost sure, not sure.
• Reusable infrastructure: the Clarke $\epsilon$-subdifferential and its closed graph, Lebourg's mean value theorem for the published Clarke subdifferential, Lemma 2.1 (a finite-dimensional minimax identity), and a conditional Borel–Cantelli argument for adaptively placed targets. Contributions to any of these are welcome, as are proofs of the companion corollaries.

Selected references

• J. V. Burke, A. S. Lewis, M. L. Overton, A robust gradient sampling algorithm for nonsmooth, nonconvex optimization, SIAM J. Optim. 15(3) (2005) 751–779. https://doi.org/10.1137/030601296
• F. H. Clarke, Optimization and Nonsmooth Analysis, Wiley, 1983; reprinted SIAM Classics in Applied Mathematics 5, 1990. https://doi.org/10.1137/1.9781611971309
• A. A. Goldstein, Optimization of Lipschitz continuous functions, Math. Programming 13 (1977) 14–22. https://doi.org/10.1007/BF01584320
• K. C. Kiwiel, Convergence of the gradient sampling algorithm for nonsmooth nonconvex optimization, SIAM J. Optim. 18(2) (2007) 379–388. https://doi.org/10.1137/050639673
• J. V. Burke, F. E. Curtis, A. S. Lewis, M. L. Overton, L. E. A. Simões, Gradient sampling methods for nonsmooth optimization, in Numerical Nonsmooth Optimization, Springer, 2020. https://arxiv.org/abs/1804.11003

Why discretize a backward SDE

A backward stochastic differential equation (BSDE) prescribes the terminal value of a process instead of its initial value. Pardoux and Peng (1990) proved that a BSDE with Lipschitz driver has a unique adapted solution. Coupled with a forward diffusion it gives a forward–backward SDE (FBSDE), which represents the solution of a semilinear parabolic PDE (nonlinear Feynman–Kac formula) and, in mathematical finance, the price and hedging strategy of a contingent claim, including path-dependent ones such as lookback and Asian options (El Karoui, Peng and Quenez, 1997). In high dimension or with path-dependent payoffs the PDE route is impractical, so one wants a time-discretization of the FBSDE itself, with a proven rate of convergence.

Zhang (2004) gave such a scheme for terminal values that are Lipschitz functionals of the whole path of the forward diffusion, and proved that it converges in mean square at rate $|\pi|\log(1/|\pi|)$, and at rate $|\pi|$ for $L^1$-Lipschitz or Markovian terminal values. The key input is a new $L^2$-regularity estimate for the martingale integrand $Z$. The same regularity estimate underlies the analysis of the Bouchard–Touzi scheme (2004) and of regression-based Monte Carlo methods (Gobet, Lemor and Warin, 2005).

Setting

Fix $T>0$, a probability space carrying a one-dimensional standard Brownian motion $W$, and its natural filtration $\mathbb F=\{\mathcal F_t\}$ augmented by the null sets. For $x\in\mathbb R^d$ the FBSDE is

$$X_t=x+\int_0^t b(s,X_s)\,ds+\int_0^t\sigma(s,X_s)\,dW_s,\qquad Y_t=\Phi(X)+\int_t^T f(s,X_s,Y_s,Z_s)\,ds-\int_t^T Z_s\,dW_s ,\tag{2.1}$$

where $X_t\in\mathbb R^d$, $Y_t,Z_t\in\mathbb R$, $b,\sigma,f$ are deterministic functions and $\Phi$ is a deterministic functional of the path $X=(X_t)_{0\le t\le T}$.

Assumption 2.3 asks that $b,\sigma,f$ be continuous, $\tfrac12$-Hölder in time and Lipschitz in the space variables, that $\Phi$ be $L^\infty$-Lipschitz,

$$|\Phi(x_1)-\Phi(x_2)|\le K\sup_{0\le t\le T}|x_1(t)-x_2(t)|$$

for càdlàg paths $x_1,x_2$, all with one constant $K>0$, and that $\sup_t\{|b(t,0)|+|\sigma(t,0)|+|f(t,0,0,0)|\}+|\Phi(0)|\le K$. $\Phi$ is $L^1$-Lipschitz if the supremum can be replaced by $\int_0^T|x_1(t)-x_2(t)|\,dt$.

A partition is $\pi:0=t_0<\dots<t_n=T$, with $\Delta t_i=t_i-t_{i-1}$ and $|\pi|=\max_i\Delta t_i$; it is $\kappa$-uniform if $\Delta t_i\ge|\pi|/\kappa$ for all $i$. The Euler scheme is $X^\pi_{t_0}=x$, $X^\pi_{t_{i}}=X^\pi_{t_{i-1}}+b(t_{i-1},X^\pi_{t_{i-1}})\Delta t_i+\sigma(t_{i-1},X^\pi_{t_{i-1}})(W_{t_i}-W_{t_{i-1}})$, and its step process is $\hat X^\pi_t=X^\pi_{t_{i-1}}$ on $[t_{i-1},t_i)$. With $\xi^\pi=\Phi(\hat X^\pi)$ the backward scheme is $Y^\pi_{t_n}=\xi^\pi$ and, for $t\in[t_{i-1},t_i)$,

$$Y^\pi_t=Y^\pi_{t_i}+f\big(t_i,X^\pi_{t_i},Y^\pi_{t_i},Z^{\pi,1}_{t_i}\big)\Delta t_i-\int_t^{t_i}Z^\pi_r\,dW_r,\qquad Z^{\pi,1}_{t_i}=\frac1{\Delta t_{i+1}}E\Big\{\int_{t_i}^{t_{i+1}}Z^\pi_r\,dr\,\Big|\,\mathcal F_{t_i}\Big\},$$

with $Z^{\pi,1}_{t_n}=0$. The numerical approximations are the step processes $\hat Y^\pi_t=Y^\pi_{t_{i-1}}$ and $\hat Z^\pi_t=Z^{\pi,1}_{t_{i-1}}$ on $[t_{i-1},t_i)$.

Formalization targets

Goal: Theorem 6.1, (6.5)–(6.6)

Assume Assumption 2.3 and that $Z$ is càdlàg. If $f$ does not depend on $z$, or $\pi$ is $\kappa$-uniform, then

$$\sup_{0\le t\le T}E\{|Y_t-\hat Y^\pi_t|^2\}+E\Big\{\int_0^T|Z_t-\hat Z^\pi_t|^2dt\Big\}\le C(1+|x|^2)\,|\pi|\log\frac1{|\pi|}\tag{6.5}$$

for $|\pi|\le e^{-1/2}$, and, if moreover $\Phi$ is $L^1$-Lipschitz or $\Phi(X)=g(X_T)$, the same left-hand side is at most $C(1+|x|^2)|\pi|$ (6.6). The constant $C$ depends only on $T$, $K$, $\kappa$ and $d$. The goal fixes the rates and leaves $C$ unspecified.

Milestones

1. Lemma 3.2: $\|Z_t\|_p\le C_p(1+|x|)$, and the time-regularity $E\{|X_t-X_{t_{i-1}}|^2+|Y_t-Y_{t_{i-1}}|^2\}\le C(1+|x|^2)|\pi|$.
2. Lemma 3.3 and Corollary 3.4, a weighted martingale inequality over a partition.
3. Theorem 3.1, the $L^2$-regularity of $Z$: $\sum_i E\int_{t_{i-1}}^{t_i}[|Z_t-Z_{t_{i-1}}|^2+|Z_t-Z_{t_i}|^2]\,dt\le C(1+|x|^2)|\pi|$.
4. Lemma 4.1, Theorem 4.2 and Corollary 4.4: errors of the Euler scheme, its step process, and $\Phi(\hat X^\pi)$.
5. Lemma 5.4 (a backward discrete Gronwall inequality), Theorem 5.3 and Remark 5.5 (the scheme at the grid points), and Theorem 5.6 (the step processes).

Significance

Theorem 6.1 gives an explicit, implementable scheme for FBSDEs with path-dependent terminal values and a mean-square rate that is sharp up to the logarithm: by Remark 4.3 the factor $\log(1/|\pi|)$ cannot be removed for $L^\infty$-Lipschitz functionals, since it is already present in the uniform error of the Euler step process. Theorem 3.1 is reusable on its own: any time-discretization of a BSDE whose error is measured against $Z_{t_{i-1}}$ needs it.

All results are proved in the paper. As far as is known, none has a machine-checked proof; the platform has no prior BSDE discretization result. The work is formalizing the known proofs, which needs the $L^2$ Itô integral, Itô's isometry and martingale representation, the standard a priori estimates for SDEs and BSDEs, and Doob's maximal inequality in continuous time.

Difficulty

The obvious argument compares the BSDE and the scheme step by step and closes with a discrete Gronwall inequality (Lemma 5.4). Its error terms include $E\int_{t_{i-1}}^{t_i}|Z_r-Z_{t_{i-1}}|^2dr$, and no pointwise estimate $E|Z_t-Z_s|^2\le C|t-s|$ is available: $Z$ is only square integrable in time, with no continuity modulus in general. The proof therefore needs the summed $L^2$-regularity of Theorem 3.1, whose proof goes through smooth approximations of the coefficients and of $\Phi$, the representation of $Z$ by the variational process $\nabla X$, and the martingale inequality of Lemma 3.3. The $\log(1/|\pi|)$ rate comes from the maximum of $n$ Gaussian increments, which an estimate of $\sup_tE$ cannot see.

Formalization scope

Time is ℝ≥0; the state space is EuclideanSpace ℝ (Fin d); the Brownian motion is the coordinate $0$ of a Fin 1-indexed standard Brownian motion; $\mathbb F$ is the augmented filtration. These, the $L^2$ Itô integral and the class $L^2(\mathbb F)$ are reused from the published modules Peng1990_SMP_Stochastic and ReflectedBSDE_Existence_Setting. Every expectation and time integral of a square is a lower Lebesgue integral in $[0,\infty]$, so a non-integrable error cannot make a bound vacuous; this rules out the trivializing formalization in which a Bochner integral of a non-integrable error is $0$.

The solutions $(X,Y,Z)$ and $(Y^\pi,Z^\pi)$ are quantified as solutions of their equations; the Euler scheme, $\hat X^\pi$, $Z^{\pi,1}$, $\hat Y^\pi$, $\hat Z^\pi$ are constructed. Each constant $C$ is chosen before the probability space, the data, the solution and the partition. Disclosed readings: (i) the log-rate statements assume $|\pi|\le e^{-1/2}$, the paper's "$\pi$ fine enough" (p. 478); (ii) "$Z$ is càdlàg" includes $Z_T=Z_{T-}$; (iii) the uniformity constant of Definition 5.2 is $\kappa$, not $K$, and $C$ may depend on it; (iv) the Hölder constant in time is $K$; (v) $\hat X^\pi_T=X^\pi_T$, $\hat Y^\pi_T=Y^\pi_T$, $\hat Z^\pi_T=0$; (vi) Lemma 3.3 uses $|\Lambda_{t_{i-1}}|$ in place of $\Lambda_{t_{i-1}}$, which strengthens it; (vii) Lemma 5.4 assumes $C\ge0$. Completeness of the probability space is not used.

Not formalized here: Lemma 2.2 (smooth approximation of $\Phi$, cited from Ma–Zhang), the background Lemmas 2.4–2.7, and the representation (6.3)–(6.4) of the scheme by functions of the Euler grid values. Welcome contributions: the SDE and BSDE a priori estimates under Lipschitz coefficients, Itô's isometry and martingale representation for the reused Itô integral, and Doob's $L^2$ inequality in continuous time; all are reusable well beyond this mission.

Selected references

• J. Zhang, A numerical scheme for BSDEs, Ann. Appl. Probab. 14(1), 459–488, 2004. https://doi.org/10.1214/aoap/1075828058
• E. Pardoux and S. Peng, Adapted solution of a backward stochastic differential equation, Systems Control Lett. 14, 55–61, 1990. https://doi.org/10.1016/0167-6911(90)90082-6
• N. El Karoui, S. Peng and M. C. Quenez, Backward stochastic differential equations in finance, Math. Finance 7(1), 1–71, 1997. https://doi.org/10.1111/1467-9965.00022
• J. Ma and J. Zhang, Representation theorems for backward stochastic differential equations, Ann. Appl. Probab. 12(4), 1390–1418, 2002. https://doi.org/10.1214/aoap/1037125868
• B. Bouchard and N. Touzi, Discrete-time approximation and Monte-Carlo simulation of backward stochastic differential equations, Stochastic Process. Appl. 111(2), 175–206, 2004. https://doi.org/10.1016/j.spa.2004.01.001
• E. Gobet, J.-P. Lemor and X. Warin, A regression-based Monte Carlo method to solve backward stochastic differential equations, Ann. Appl. Probab. 15(3), 2172–2202, 2005. https://doi.org/10.1214/105051605000000412