arXiv CSopen access

Barzilai-Borwein Fails Superlinear Convergence on an Open Set of Quadratics for Every Dimension $n\geq 4$

Barzilai–Borwein Fails Superlinear Convergence on an Open Set of Quadratics for Every Dimension n ≥ 4 Dawei Li

Xiaotian Jiang Mingyi Hong University of Minnesota {li004678,jian0851,mhong}@umn.edu July 14, 2026

arXiv:2607.21579v1 [math.OC] 23 Jul 2026

Abstract Barzilai–Borwein (BB) method has shown strong practical performance in continuous optimization, yet its convergence dynamics remains poorly understood. In particular, a central unresolved question is whether BB converges superlinearly for almost every strictly convex quadratic problem and initialization. We provide a negative answer to this question. Specifically, for every finite dimension n ≥ 4, we construct a nonempty open, hence positive-Lebesgue-measure, family of strictly convex quadratic problems and initial points for which the long Barzilai–Borwein method (BB1) converges but cannot converge root-superlinearly. More precisely, with the explicit constants ρ = 10−6 , ρ = 0.61, every spectral component of the gradient is bounded above and below by the corresponding geometric sequence. Consequently, the gradient norm and the energy norm of the error satisfy twosided geometric estimates with the same rates, while the objective gap satisfies the corresponding estimates with squared rates. In particular, all three quantities are bounded below by geometric sequences, ruling out superlinear convergence. The construction is highly nontrivial, based on a computer-assisted proof of a nonresonant, attracting seven-cycle of the projectivized BB dynamics in dimension four.

1

Introduction

Gradient descent is one of the oldest and most widely used methods in continuous optimization. Each iteration requires only a gradient and a few vector operations, making it attractive in large-scale problems. Its performance, however, can deteriorate on ill-conditioned settings: exact line search may produce zigzag iterates, whereas a good constant step requires spectral information that is rarely available in advance. The Barzilai–Borwein (BB) method (Barzilai and Borwein, 1988) addresses this difficulty without changing the negative-gradient direction. Instead, it extracts a scalar curvature estimate from the two most recent iterates. This change of step-size rule, rather than search direction, can substantially improve the practical behavior of a basic gradient iteration (Fletcher, 2005; Birgin et al., 2014). For a differentiable objective f : Rn → R, write gk = ∇f (xk ), sk−1 = xk − xk−1 , and yk−1 = gk − gk−1 . The two classical BB step sizes are (Barzilai and Borwein, 1988) αkBB1 =

sT k−1 sk−1 , T sk−1 yk−1

αkBB2 =

sT k−1 yk−1 , T y yk−1 k−1

1

xk+1 = xk − αk gk .

(1)

The first choice makes the scalar matrix (αkBB1 )−1 I the least-squares solution of the secant equation; the second makes αkBB2 I a least-squares approximation of the inverse Hessian. They are therefore often called the long and short BB steps, respectively. The formulation of the two BB steps are different, but they possess similar properties in training. BB runs at essentially the cost of gradient descent while carrying curvature information without forming a matrix and, in its basic form, without performing a line search. On a quadratic, the reciprocal BB step is a Rayleigh quotient of the Hessian; this is also the origin of the name spectral gradient method. Despite this simple construction, BB generates nonlinear, history-dependent dynamics even on quadratic objectives. A longstanding open question is whether BB converges superlinearly for almost every strictly convex quadratic problem and initialization, or whether failure of root-superlinear convergence can persist on a set of positive measure rather than only on exceptional instances. More precisely, the question is whether lim ∥gk ∥1/k = 0 k→∞

for almost every strictly convex quadratic problem and initialization. We aim to resolve this question in this paper. Significance and applications. The BB step has become an algorithmic primitive rather than a single method. Raydan combined it with a non-monotone globalization strategy for large-scale unconstrained optimization (Raydan, 1997); related non-monotone frameworks were developed in Grippo and Sciandrone (2002) and Zhang and Hager (2004), and the spectral projected gradient method brought the same step selection to optimization over convex sets (Birgin et al., 2000, 2001; Dai and Fletcher, 2005b; Andreani et al., 2005). BB-type spectral steps now appear in support-vector-machine training (Serafini et al., 2005), sparse reconstruction and compressed sensing (Figueiredo et al., 2007; Wright et al., 2009; Van Den Berg and Friedlander, 2009; Wen et al., 2010), image restoration (Wang and Ma, 2007; Bonettini et al., 2009), nonnegative matrix factorization (Huang et al., 2015), optimization with orthogonality constraints (Wen and Yin, 2013), Riemannian optimization (Iannazzo and Porcelli, 2018), stochastic and variance-reduced methods (Tan et al., 2016), distributed optimization with locally computed step sizes (Gao et al., 2022), and PDEconstrained optimization in Hilbert space (Azmi and Kunisch, 2020, 2022); more recent safeguarded and adaptive constructions seek global guarantees on broader function classes (Burdakov et al., 2019; Zhou et al., 2026). These descendants alter the pure iteration in essential ways, such as projection, regularization, line searches, stochastic averaging, safeguarding, but all of them inherit its step-size mechanism. Understanding the unmodified deterministic iteration is therefore valuable both in its own right and as a foundation for interpreting its many descendants. Earlier convergence and rate theory. The known convergence theory of BB mostly focuses on low-dimensional problems. Barzilai and Borwein (1988) proved R-superlinear convergence for two-dimensional strictly convex quadratics. Raydan (1993) established global convergence for arbitrary finite-dimensional strictly convex quadratics, and Dai and Liao (2002) proved an R-linear rate in arbitrary dimension, together with a local consequence for sufficiently smooth nonquadratic objectives near a nondegenerate minimizer. In dimension two, Dai (2013) related the rate to the initial spectral balance and showed that superlinear behavior holds for almost every initialization, while linear behavior occurs on an exceptional set. Dai and Fletcher (2005a) developed an asymptotic framework for gradient methods with retarded spectral information and, by computing the asymptotics of BB and its relatives, observed a transition from superlinear to linear behavior at some dimension n ≥ 4 depending on the method. Their statements are asymptotic

2

computations and numerical observations, and no rigorous lower-rate mechanism in dimension n ≥ 4 was identified. For general dimension, rate bounds have been substantially less informative about the observed fast behavior. A recent refinement by Li and Sun (2021) proves, for BB1 on an arbitrary finitedimensional strongly convex quadratic with condition number κ = λn /λ1 , an R-linear upper factor 1 − 1/κ. This is comparable in condition-number scaling to the exact-line search steepest-descent factor. The same work gives a lower-rate example supported only on the two extreme eigenspaces, with exact factor (κ − 1)/(κ + 1), ruling out superlinear convergence. However, for n ≥ 2 its initial point lies in a Lebesgue-null set, so an arbitrarily small perturbation could possibly accelerate the convergence; the same paper explicitly raises the convergence rate under generic initialization as an open question. Parallel lines of work have sought either to globalize BB or to design spectral rules with more transparent asymptotics. These include gradient methods with retards (Friedlander et al., 1998), alternate and cyclic BB schemes (Dai, 2003; Dai et al., 2006), convex combinations of the long and short BB steps (Dai et al., 2018), and steps engineered for two-dimensional termination or favorable spectral limits (Yuan, 2006; Huang et al., 2021, 2022). This literature gives considerable evidence that low-dimensional spectral mechanisms govern the performance of many gradient methods. What has been missing for the original BB1 method is an answer to the following question: In dimension n ≥ 4, is every linear-rate BB1 example an unstable exception, or does genuinely full-dimensional linear-rate behavior persist on a robust family? Contribution of the present lower bound. This paper answers this question for BB1. For every finite dimension n ≥ 4 and every linear term b ∈ Rn , we construct a nonempty open set An ⊂ Sn++ of simple-spectrum matrices such that, for every A ∈ An , there is a nonempty open set Xn (A) of initial points for which the admissibly initialized BB1 iteration is well defined and converges to x∗ = A−1 b. More precisely, if Πi (A) denotes the orthogonal spectral projector associated with the ith eigenvalue, then 10−6k ∥Πi (A)g0 ∥ ≤ ∥Πi (A)gk ∥ ≤ 0.61k ∥Πi (A)g0 ∥,

i = 1, . . . , n,

k ≥ 0.

(2)

Here “admissibly initialized” means that the first step is the inverse gradient Rayleigh quotient, α0 =

g0T g0 , g0T Ag0

x1 = x0 − α0 g0 ,

so the orbit is an actual BB1 trajectory rather than an arbitrarily prescribed pair of projective states. The joint set Ωn = {(A, x0 ) : A ∈ An , x0 ∈ Xn (A)} is open and nonempty in Sn++ × Rn . It therefore has positive Lebesgue measure, as does every initial-point fiber Xn (A). Summing (2) over the orthogonal spectral components gives the same two-sided geometric estimate for the gradient norm and the energy norm of the error; the objective gap satisfies the corresponding squared estimate. In particular, lim inf ∥gk ∥1/k ≥ 10−6 > 0, k→∞

so none of these trajectories is root-superlinearly or Q-superlinearly convergent. An almosteverywhere superlinearity statement on the joint problem-data space is therefore false in every dimension n ≥ 4; in particular, the question of Li and Sun (2021) about generic initialization is answered in the negative on an open set of matrices. The result strengthens the earlier lower-rate picture in three ways. 3

(i) Robustness. The slow behavior is not confined to a specially balanced endpoint subspace or to a zero-measure initialization. It persists under simultaneous perturbations of the spectrum and the initial point. (ii) A dynamical mechanism. The proof identifies a nonresonant, attracting period-seven orbit of the projectivized four-dimensional BB1 dynamics. A certified fifteen-step connection reaches its basin from the admissible-initialization manifold, and a spectral translation makes every radial multiplier contractive. Attraction of the projective cycle, rather than invariance of a coordinate subspace, creates the open basin and the lower geometric rate. (iii) Full-dimensional persistence. Four spectral modes remain active along the asymptotic cycle. A uniform transverse Floquet estimate then embeds the construction into each finite dimension n > 4 while keeping every additional component nonzero and geometrically controlled. The normalized spectral measure consequently approaches a seven-periodic measure rather than collapsing to the two endpoint modes familiar from exact steepest descent (Akaike, 1959; Forsythe, 1968; Pronzato et al., 2006). The periodic orbit is certified by interval arithmetic and a contraction argument; it is not inferred from a long floating-point simulation. This computer-assisted component establishes exact existence, nonresonance, and attraction inside explicit rational boxes; openness, positive measure, and persistence in higher dimension then follow from analytic stability arguments. The theorem is an existence result, not a universal rate formula. It concerns BB1, not every BB variant. The set An is a small open neighborhood of one explicitly constructed spectrum—three separated eigenvalues together with, for n > 4, the remaining eigenvalues clustered in a short interval near the top—so the phrase “every finite dimension” asserts that such a neighborhood exists for each n, not that its size is uniform in n. The constants 10−6 and 0.61 are conservative separation margins rather than claimed optimal asymptotic factors. The contribution is structural: it proves that the dimension-four transition from low-dimensional superlinear behavior to linear-rate behavior can be generated by a stable recurrent attractor and can occupy positive measure. In this sense, the lower bound complements the classical and recent R-linear upper bounds. The upper theory controls how slow BB1 can be; the present result shows that a positive linear root factor is actually realized robustly, and explains the orbit geometry that sustains it. AI-Assisted Development. The sketch of the proof of the main result of this paper was obtained by GPT 5.6 Sol. GPT was given an equivalent problem formulation (component-wise quotient formulation) and asked to construct the lower bound example. The prompt is provided in Appendix B. The construction by GPT is highly non-trivial. Human verification and polishing were done afterwards.

2

The BB1 method and the main theorem

Let Sn++ denote the open cone of real symmetric positive-definite n × n matrices. It is an open subset of the Euclidean space Sn , which has dimension n(n + 1)/2. Fix A ∈ Sn++ and b ∈ Rn , and consider 1 fA,b (x) = xT Ax − bT x, 2

x∗ := A−1 b,

The unique minimizer is x∗ because A ≻ 0.

4

g(x) := ∇fA,b (x) = Ax − b.

(3)

We study the long Barzilai–Borwein method, usually called BB1 Barzilai and Borwein (1988); classical convergence analyses on strictly convex quadratics include Raydan (1993); Dai and Liao (2002). The initialization is part of the theorem and is therefore stated explicitly. Given x0 = ̸ x∗ , put g0 = Ax0 − b and g T g0 α0 := T0 , x1 := x0 − α0 g0 . (4) g0 Ag0 For k ≥ 1, provided the quantities are defined, set sk−1 := xk − xk−1 , yk−1 := gk − gk−1 , αk :=

sTk−1 sk−1 , sTk−1 yk−1

xk+1 := xk − αk gk ,

gk := Axk − b.

(5)

The choice (4) is the inverse Rayleigh quotient of the initial gradient. It makes the first two reciprocal step sizes equal. Main Result.

For a simple-spectrum matrix A, let 0 < λ1 (A) < · · · < λn (A)

be its ordered eigenvalues and let Πi (A) be the orthogonal projector onto the one-dimensional eigenspace corresponding to λi (A). The orientation of an eigenvector is irrelevant because all conclusions below use ∥Πi (A)gk ∥. Theorem 2.1 (Positive-measure lower rate for BB1). Fix an integer n ≥ 4 and a vector b ∈ Rn . There exist a nonempty open set An ⊂ Sn++ consisting of simple-spectrum matrices and, for every A ∈ An , a nonempty open set Xn (A) ⊂ Rn \ {x∗ } with the following properties. For every x0 ∈ Xn (A), the BB1 iteration (4)–(5) is well-defined for every k ≥ 0, converges to x∗ , and satisfies, for every i ∈ {1, . . . , n} and every k ≥ 0, ρk ∥Πi (A)g0 ∥ ≤ ∥Πi (A)gk ∥ ≤ ρk ∥Πi (A)g0 ∥ ,

ρ = 10−6 ,

ρ = 0.61.

In particular, none of the initial spectral components is zero. The union Ωn := {(A, x) : A ∈ An , x ∈ Xn (A)} ⊂ Sn++ × Rn

(6)

(7)

has positive Lebesgue measure in Sn × Rn , and each fiber

is nonempty and open. Consequently Ωn Xn (A) has positive Lebesgue measure in Rn . Moreover, writing ek = xk − x∗ and ∥e∥A = (eT Ae)1/2 , one has ρk ∥g0 ∥ ≤ ∥gk ∥ ≤ ρk ∥g0 ∥ , k

(8)

k

ρ ∥e0 ∥A ≤ ∥ek ∥A ≤ ρ ∥e0 ∥A ,

(9)

ρ2k f (x0 ) − f (x∗ ) ≤ f (xk ) − f (x∗ ) ≤ ρ2k f (x0 ) − f (x∗ ) .

(10)

Thus 1/k

lim inf ∥gk ∥1/k ≥ 10−6 ,

lim inf ∥ek ∥A ≥ 10−6 ,

k→∞

k→∞

(11)

and these BB1 orbits are not root-superlinearly or Q-superlinearly convergent. The proof occupies Sections 4–9. The constants are conservative and are not claimed to be optimal. The result is an existence theorem for an open family; it does not assert the same lower bound for every positive spectrum. 5

Why the BB dynamics are difficult. Even for the quadratic objective (3), BB1 generates a nonlinear delayed dynamical system. Let 0 < λ1 ≤ · · · ≤ λn be the eigenvalues of A, let v1 , . . . , vn P be an orthonormal eigenbasis, and write gk = i dik vi . Then for BB1 on this quadratic, j 2 j=1 (dk−1 ) , j 2 j=1 λj (dk−1 )

Pn

dik+1 = (1 − λi αk )dik ,

αk = Pn

(12)

equivalently, j 2 j=1 (λj − λi )(dk−1 ) . Pn j 2 j=1 λj (dk−1 )

Pn

dik+1 = dik

(13)

Equation (13) displays the first source of difficulty. Diagonalization completely decouples fixedstep gradient descent, but it does not decouple BB1. The update of the ith component di from time k to time k +1 depends on every component at time k −1. Thus, although the objective is diagonal in the eigenbasis, its spectral modes remain coupled through a delayed Rayleigh quotient. This coupled and delayed feedback produces nonlinear behavior. Depending on the Rayleigh quotient selected by the preceding gradient, the multiplier of a spectral component can be negative, close to zero, or greater than one in absolute value. Consequently, components can change sign, high-eigenvalue components can grow, and different eigendirections may take turns dominating the gradient, so neither the objective value nor the gradient norm is monotone. (On nonquadratic objectives, the situation is more delicate – sT k−1 yk−1 may fail to remain positive, and the unsafeguarded method can take excessively long steps or fail to converge even under strong convexity (Fletcher, 2005; Burdakov et al., 2019) – but they are not the subject of this paper.) The second difficulty is determining the actual asymptotic rate. As a component multiplier can be close to zero, a spectral component may become arbitrarily small in a single iteration. An upper estimate ∥gk ∥ ≤ Cq k therefore cannot distinguish genuinely linear convergence from faster convergence, such as superlinear. Such an upper bound also does not determine which spectral modes survive, whether the normalized state approaches a fixed point or a cycle, or what lower asymptotic rate is actually attained. Ruling out superlinear convergence needs a lower bound, or an equally precise description of the asymptotic orbit. There is a third difficulty in answering an almost-everywhere question. A single specially constructed lower bound example could be an isolated exception and hence belong to a set of Lebesgue measure zero (Dai, 2013; Li and Sun, 2021). In practice, such an example may depart from the theoretical trajectory and lead to faster convergence rate due to numerical instability He et al. (2025). A robust lower bound example should additionally prove that the convergence rate persists under perturbations of both the spectrum and the initialization. This is the aspect of BB dynamics for which the classical theory has been least complete.

3

Outline of the proof

The lower bound is not obtained by estimating a generic BB1 trajectory. The proof constructs one rigorously certified orbit of a low-dimensional dynamical system, verifies that this orbit stays uniformly away from every resonance, and then lets stability and continuity enlarge the single orbit into an open family. Figure 2 displays the chain of implications; this section walks through it once. Reduction to a factor bound (§4). On a quadratic, BB1 acts on each spectral component
of the gradient through the exact multiplier identity dik+1 = (ak − λi )/ak dik , where ak = αk−1 is

6

(a) one-step multipliers

(b) spectral components of the gradient

|ak − λi |/ak

0 10 20 30 40 50

i=1

i=2

k 10−

0

10

i=3 i=4 0.61 kkg0 k

observed factor β 1/7 ≈ 0.144

6k g 0k

10 1 10 2 10 3 10 4 10 5 connection k < 15 lower bound 10−6 i=1 i=2 i=3 i=4 10 6 0 10 20 30 40 50 iteration k

log 10 |dki |

upper bound 0.61

100

20

30

40 iteration k

50

60

Figure 1: Numerical realization of the factor bounds along the certified orbit at the central shifted spectrum λ∗ , computed in 200-digit arithmetic. (a) The one-step multipliers |ak − λi | /ak : after the fifteen-step connection (shaded) they are seven-periodic and remain inside the certified band [10−6 , 0.61]; the transient minimum, approximately 1.7 × 10−5 , occurs at k = 3. (b) The four spectral components of the gradient decay in parallel between the two geometric envelopes; the observed asymptotic per-step factor is β 1/7 ≈ 0.144, where β ≈ 1.29 × 10−6 is the contraction of the raw coordinates over one seven-cycle. The figure is an illustration; the rigorous statements are those of Certificate 5.1. the delayed Rayleigh quotient, so that dik = di0 t<k |at − λi | /at . All conclusions of Theorem 2.1 therefore follow from a single uniform one-step inequality, Q

10−6 ≤

|ak − λi | ≤ 0.61 ak

for every i and every k ≥ 0:

(14)

the lower bound is a nonresonance condition (ak never comes too close to an eigenvalue), the upper bound forces convergence. Since the gradient itself tends to zero, its scale is removed: the pair (ak , pk ) of the reciprocal step size and the normalized squared spectral weights evolves under a rational map Fλ on R × ∆n−1 , with the admissible first step encoded by an initialization map Dλ . If this projective state approaches a periodic orbit that is separated from the spectrum, all late multipliers in (14) lie in a fixed compact subinterval of (0, 1), and the finite initial segment is handled by continuity. The certified periodic object (§5). The construction is carried out in dimension four. The unknowns are a free eigenvalue λ̄3 , a projective state z, and an admissible initial weight r; they satisfy eight rational equations, namely Fλ̄7 (z) = z (a seven-cycle) together with Fλ̄15 (Dλ̄ (r)) = z (an exact fifteen-step landing). The landing block is forced by a dimension count: the projective state space is four-dimensional while the admissibly initialized states form a three-dimensional hypersurface, so an attracting cycle need not meet the forward orbit of any genuine BB1 initialization. A Newton-like map H(u) = u − RG(u) with an exact rational preconditioner R is then evaluated in directed interval arithmetic on a box X of radius 10−55 ; the certified inclusions H(X) ⊂ int X and supX ∥DH∥∞ < 1 make Banach’s theorem produce an exact root u∗ (Lemma 5.2). The same run certifies that every normalization denominator (namely Zt ) is positive, that the cycle and the connection stay separated from all eigenvalues, and that an explicit Lyapunov pair (P, Q7 ) is positive 7

definite. Two logically distinct contractions appear here: H only certifies the root; attraction is proved separately. Attraction, rates, and openness in dimension four (§6–§7). The certified inequality P − AT7 P A7 ≻ 0 gives ∥A7 ∥P < 1, so the seven-cycle is Schur stable and possesses an open basin (Lemma 6.3). The rate constants come from two independent mechanisms. A translation λ 7→ λ + 41 commutes with the projective dynamics (Lemma 7.1) and moves the spectrum into (4.99, 8.01) with width below 3.02; since every ak is a convex combination of eigenvalues, this alone yields the upper factor 3.02/4.99 < 0.61. The certified separation yields the lower factor 10−5 /8.01 > 10−6 . Because 1 is not an eigenvalue of A7 , the implicit-function theorem continues the cycle to all nearby spectra, and finite-time continuity of the fifteen-step map carries all nearby admissible initializations into the trapping tube: one trajectory becomes an open product L4 × P4 (Proposition 7.2). Every dimension, and the lift to BB1 problems (§8–§9). For n > 4, the extra eigenvalues are placed in (7.49, 7.51) with zero initial weight, so the four-dimensional cycle survives on an invariant simplex face. The seven-step derivative is block upper triangular there, and the certificate bounds every transverse multiplier by τ (ν) < 0.049, so the boundary cycle attracts in the new directions as well; tilting the extra weights to small positive values and continuing in all n eigenvalues gives open sets Ln × Pn for each fixed finite n (Proposition 8.1), with the rate constants—though not the neighborhood sizes—uniform in n. Finally, three changes of variables lift the projective family √ to problem data: the orthant polar map di0 = σi r p0,i , the continuity of the spectral-projector weights p(A, g) on simple-spectrum matrices, and the diffeomorphism (A, x) 7→ (A, Ax − b). The resulting set Ωn of problems and initial points is open and nonempty, hence of positive Lebesgue measure, and Theorem 2.1 follows.

4

Exact reduction of BB1 to a projective rational map

We first derive the scalar recurrence without suppressing the exceptional index k = 0. Lemma 4.1 (The BB1 step on a quadratic). Suppose A ≻ 0 and gk−1 ̸= 0. If xk = xk−1 − αk−1 gk−1 with αk−1 ̸= 0, then the BB1 quotient in (5) is positive and equals αk =

T g gk−1 k−1 . T gk−1 Agk−1

(15)

In particular, the initialization (4) gives α1 = α0 . Proof. Because gk = Axk − b and gk−1 = Axk−1 − b, subtraction gives yk−1 = gk − gk−1 = A(xk − xk−1 ) = Ask−1 . The preceding gradient step gives sk−1 = −αk−1 gk−1 . Therefore 2 T sTk−1 sk−1 = (−αk−1 gk−1 )T (−αk−1 gk−1 ) = αk−1 gk−1 gk−1 , 2 T sTk−1 yk−1 = sTk−1 Ask−1 = αk−1 gk−1 Agk−1 . T Ag 2 Since A ≻ 0 and gk−1 ̸= 0, the denominator gk−1 k−1 is strictly positive. Since αk−1 > 0, cancellation yields (15); both its numerator and denominator are positive. Taking k = 1 and using (4) proves α1 = α0 .

8

Exact reduction (§4)
a − λi i k dk , dik+1 = (ak , pk ) = Fλk Dλ (p0 ) ak goal: 10−6 ≤ |ak − λi | /ak ≤ 0.61 for all i, k interval Newton + Banach (Lem. 5.2)

Certified periodic object (§5, dimension 4) eight equations G(u∗ ) = 0: seven-cycle F 7 (z∗ ) = z∗ ⊕ admissible landing F 15 (D(r∗ )) = z∗ ; certified Zt > 0, separations at − λ̄i , Lyapunov pair (P, Q7 ), transverse τ (ν̄) discrete Lyapunov inequality (Lem. 6.2, 6.3)

Attraction (§6) P ≻ 0, P − AT7 P A7 ≻ 0 =⇒ ∥A7 ∥P < 1 =⇒ the cycle has an open basin translation conjugacy + IFT (Lem. 7.1, Prop. 7.2)

Rates and openness in dimension 4 (§7) translation λ 7→ λ + 41 ⇒ upper factor 0.61; certified separation ⇒ lower factor 10−6 ; hyperbolicity ⇒ open family L4 × P4 transverse Floquet bound (Prop. 8.1)

Every finite dimension (§8) boundary cycle on an invariant simplex face; transverse multipliers τ (ν) < 0.049 ⇒ attraction persists ⇒ open Ln × Pn for each n ≥ 4 three changes of variables (§9)

Lift and measure (§9) Θσ , projector weights p(A, g), (A, x) 7→ (A, Ax − b) diffeomorphisms ⇒ Ωn open and nonempty ⇒ positive measure ⇒ Theorem 2.1

Figure 2: Map of the proof of Theorem 2.1. The shaded box is the computer-assisted step (Certificate 5.1); every other stage is a pen-and-paper argument that consumes its certified strict inequalities. Arrow labels name the tool carrying each implication. Fix an orthogonal diagonalization A = QΛQT ,

Λ = diag(λ1 , . . . , λn ),

0 < λ1 < · · · < λn ,

(16)

and define the spectral gradient coordinates dk := QT gk ,

dik = (QT gk )i .

(17)

Because Q is orthogonal, gkT gk =

n X

(djk )2 ,

gkT Agk = dTk Λdk =

j=1

n X

λj (djk )2 .

j=1

The gradient update is gk+1 = A(xk − αk gk ) − b = (Axk − b) − αk Agk = (I − αk A)gk . 9

(18)

Multiplication by QT and use of QT AQ = Λ give dik+1 = (1 − αk λi )dik .

(19)

For k = 0, substitution of (4) and (18) into (19) yields j 2 j=1 (λj − λi )(d0 ) . Pn j 2 j=1 λj (d0 )

Pn

di1 = di0

(20)

For k ≥ 1, Lemma 4.1 instead uses gk−1 , and hence j 2 j=1 (λj − λi )(dk−1 ) . Pn j 2 j=1 λj (dk−1 )

Pn

dik+1 = dik

(21)

Equations (20)–(21) are exactly the delayed recurrence certified below. We now remove the common magnitude of dk . Whenever dk ̸= 0, define xk,i := (dik )2 ,

Sk :=

n X

xk,j ,

pk,i :=

j=1

xk,i . Sk

(22)

Then pk belongs to the simplex (

∆n−1 :=

n

p ∈ R : pi ≥ 0,

n X

)

pi = 1 .

i=1

For p ∈ ∆n−1 and a ∈ R, define µλ (p) := Zλ (a, p) :=

n X i=1 n X

λi pi ,

(23)

pi (a − λi )2 ,

(24)

i=1

pi (a − λi )2 , Zλ (a, p)
Fλ (a, p) := µλ (p), Tλ,a (p) ,

(Tλ,a (p))i :=

(25)

Dλ (p) := µλ (p), p .

(26)

The map Tλ,a is used only if Zλ (a, p) > 0. Lemma 4.2 (Exact projective state and multiplier identities). Suppose d0 = ̸ 0, and let p0 be defined by (22). Set a0 := µλ (p0 ), ak := µλ (pk−1 ) (k ≥ 1). (27) As long as the recurrence is defined, (ak , pk ) = Fλk (Dλ (p0 ))

(k ≥ 0),

(28)

and, for every i and k ≥ 0, dik+1 = dik

ak − λ i , ak

dik = di0

k−1 Y

|at − λi | . at t=0

The empty product for k = 0 is one. In particular, a1 = a0 = µλ (p0 ). 10

(29) (30)

Proof. For k ≥ 1, the scalar quotient in (21) is P P j 2 j (λj − λi )(dk−1 ) j λj xk−1,j − λi j xk−1,j P = P j 2 j λj xk−1,j j λj (dk−1 )

P

= =

Sk−1

j λj pk−1,j − λi Sk−1

P

Sk−1

P

j λj pk−1,j

ak − λi µλ (pk−1 ) − λi = . µλ (pk−1 ) ak

The same computation applied to (20), with p0 in place of pk−1 , gives the formula for k = 0. This proves (29). Squaring (29) gives (ak − λi )2 xk+1,i = xk,i . (31) a2k Summing (31) over i and using xk,i = Sk pk,i yields Sk+1 =

X

xk,i

i

(ak − λi )2 Sk X Sk = 2 pk,i (ak − λi )2 = 2 Zλ (ak , pk ). 2 ak ak i ak

If Zλ (ak , pk ) > 0, division of (31) by this identity gives pk+1,i =

pk,i (ak − λi )2 = (Tλ,ak (pk ))i . Zλ (ak , pk )

Definition (27) also gives ak+1 = µλ (pk ). Hence Fλ (ak , pk ) = (ak+1 , pk+1 ). Since (a0 , p0 ) = Dλ (p0 ), induction on k proves (28). Finally, taking absolute values in (29) and multiplying the identities for t = 0, . . . , k − 1 proves (30). Because all λi > 0 and p is a probability vector, µλ (p) ≥ λ1 > 0. Thus ak never causes division by zero. The only remaining possible singularity is Zλ (ak , pk ) = 0; the construction below keeps every orbit in a region where Z is strictly positive.

5

The certified four-dimensional periodic object

The only computer-assisted part of the proof is the validation of one exact zero of eight rational equations and of several strict inequalities at that zero. The validation is performed by the supplied script interval_cert.py, provided in Appendix A. The script uses only Python’s standard decimal module and contains all entries of the rational preconditioner and Lyapunov matrices.

5.1

The eight equations

Begin with the unshifted four-point spectrum

λ̄ = 1, 1.8786699041860466, λ̄3 , 4 .

(32)

The displayed second eigenvalue is a terminating decimal and is interpreted as the exact corresponding rational number. Use the direct state chart z = (a, p1 , p2 , p3 ),

p4 = 1 − p1 − p2 − p3 , 11

and the initialization chart r4 = 1 − r1 − r2 − r3 .

r = (r1 , r2 , r3 ), The unknown is

u = (λ̄3 , a, p1 , p2 , p3 , r1 , r2 , r3 ) ∈ R8 . Define, wherever the iterates are well-defined, 

Fλ̄7 (z) − z



 ∈ R8 . G(u) :=  15 Fλ̄ (Dλ̄ (r)) − z

(33)

The first four equations require a point fixed by the seven-step map. The last four require an exact fifteen-step connection from an admissible BB initial state Dλ̄ (r). The landing condition is essential: a periodic state in the four-dimensional (a, p) space need not lie in the forward orbit of the three-dimensional admissible-initialization hypersurface a = µ(p). Let c = (c1 , . . . , c8 ) be the exact terminating-decimal vector c1 = 2.71278149106533722827430219171662596122131200973311943581550524557629575811, c2 = 1.33344967252388129918980034968384933800281196626470897628663867165111253656, c3 = 0.828108598362557830284706279933642675447754907755381063536594022726890227454, c4 = 0.167043051370226631389809462759502319781038286313814177038878187525592790929, c5 = 0.004717957350367412952812533939501994168806798103523267952088050658941184644, c6 = 0.0000134113750442893933116720965327494222816004773939806475307653666721805974245, c7 = 0.989171705769078149047269711515135535847128223427916068474163743526926233731, c8 = 0.00000349260028439145589133467903469364875523517774851271154130219629904643379924.

Set X := c + [−10−55 , 10−55 ]8 .

(34)

At the center, and uniformly throughout X, the eliminated probabilities satisfy p4 ∈ [0.0001303929168481253726717233673530106024000078272814911, 0.0001303929168481253726717233673530106024000078272814918],

(35)

r4 ∈ [0.0108113902555931701035272817092970210818349409169414378, 0.0108113902555931701035272817092970210818349409169414385].

(36)

Thus p and r are in the relative interior of ∆3 , and 1 < 1.8786699041860466 < λ̄3 < 4 throughout X.

5.2

What the interval calculation proves

Computer-assisted certificate 5.1 (Validated root, stability, and separation). With G, c, and X defined above, the supplied directed-rounding calculation proves all of the following statements. (C1) There is an explicitly recorded exact rational matrix R ∈ R8×8 such that, with H(u) = u − RG(u),
H(X) ⊂ c − RG(c) + I − R[DG(X)] (X − c) ⊂ int X. (37) The largest coordinate ratio of the final box radius to 10−55 is less than 3.876379 × 10−20 . 12

(C2) The same derivative enclosure satisfies sup ∥DH(u)∥∞ ≤ ∥I − R[DG(X)]∥∞ < 7.200691 × 10−27 < 1.

(38)

u∈X

(C3) Every normalization denominator occurring in the seven cycle steps and the fifteen connection steps is positive. More explicitly, the same outward evaluation gives min Ztconn > 0.0009510722601452239.

min Zt > 0.1516357118221540,

0≤t<7

0≤t<15

(39)

(C4) For every exact u ∈ X, the finite connection and cycle obey the uniform separation bounds min at − λ̄i > 1.3886903675 × 10−4 ,

(40)

min at − λ̄i > 7.5691389311 × 10−2 .

(41)

0≤t≤15 1≤i≤4 0≤t<7 1≤i≤4

Both statements remain true, with the same displayed lower bounds, when an additional test eigenvalue ν̄ ∈ [3.49, 3.51] is included in the minimum. (C5) Let A7 = Dz Fλ̄7 (z), with z ranging over the z-part of X. The script records an exact symmetric rational matrix P and proves, for the exact root obtained below, P ≻ 0,

Q7 := P − AT7 P A7 ≻ 0.

(42)

The outward-rounded LDLT pivots for P have the positive lower bounds 977656.6634223271,

65683881.36341604,

9.217413243679623,

323.9590422006289, (43)

and the row diagonal-dominance margins for Q7 have the lower bounds 0.9999999999999999999976374,

0.9999999999999999994270558,

0.9999999999999999994222700,

0.9999999999999999993948258.

(44)

(C6) For every ν̄ ∈ [3.49, 3.51], the seven-step transverse multiplier defined in (67) below satisfies 0.0230006289035345 < τ (ν̄) < 0.0482036188628888 < 1.

(45)

(C7) For the first two phases of the cycle, a0 − a1 > 0.1782019618055700634.

(46)

These machine statements are rigorous real inequalities by the standard directed-rounding validation principle (Krawczyk, 1969; Rump, 2010). A decimal literal is first converted to a Decimal exactly. For each addition, subtraction, multiplication, and reciprocal, the script evaluates once with rounding toward −∞ and once with rounding toward +∞; if the input intervals contain the exact real inputs, the output interval therefore contains the exact result of that elementary operation. Structural induction over an arithmetic expression then shows that its interval evaluation encloses its entire real range over the input box. The automatic-differentiation class stores a pair (v, d) of an interval value and a vector of interval partial derivatives and propagates them by the exact 13

(b) distance to resonance

8

λ4

7

λ3

6

λ2

5

λ1

connection k < 15 seven-periodic for k ≥ 15

0

10

20

30

iteration k

40

100 min |a − λi | i k

reciprocal step size ak

(a) ak against the spectrum 10 1

certified cycle bound 7.57 × 10−2

10 2 10 3

certified transient bound 1.39 × 10−4

10 4

50

0

10

20

30

iteration k

40

50

Figure 3: The certified projective orbit at the central parameter. (a) The reciprocal step size ak against the eigenvalues: the admissible initialization lands on the seven-cycle after exactly fifteen steps and then repeats. (b) The distance from ak to the nearest eigenvalue. The certified separation bounds of Certificate 5.1(C4) are sharp for this orbit: the transient minimum at k = 3 and the cycle minimum match the certified lower bounds to the displayed digits. identities D(u + v) = Du + Dv, D(uv) = vDu + uDv, D(u−1 ) = −u−2 Du; the same induction, applied simultaneously to values and derivatives, shows that the computed matrix [DG(X)] contains DG(u) entrywise for every u ∈ X. Each reciprocal operation asserts that its denominator interval excludes zero, and the explicit positive bounds (39) strengthen this check. All magnitudes used are far from the decimal exponent limits, so no overflow or underflow is involved. The centered enclosure in (37) follows from the identity H(u) − c = H(c) − c +

Z 1

DH(c + t(u − c))(u − c) dt

0

= −RG(c) +

Z 1

I − RDG(c + t(u − c)) (u − c) dt.

0

valid because X is convex and H is C 1 on a neighborhood of X. Every derivative matrix in the integral belongs entrywise to I − R[DG(X)], and interval matrix–box multiplication encloses every product with u − c ∈ X − c as well as its average over t ∈ [0, 1]. The second inclusion in (37) is the directed-rounding output. We next turn the interval statements into exact analytic conclusions. Lemma 5.2 (Existence and uniqueness of an exact solution). There is a unique u∗ ∈ X satisfying G(u∗ ) = 0. Proof. The box X is a nonempty closed subset of (R8 , ∥·∥∞ ) and hence a complete metric space. Certificate 5.1(C1) gives H(X) ⊂ int X ⊂ X, so H maps X into itself, and by Certificate 5.1(C3) every denominator is nonzero on an open neighborhood of X, so H is C 1 there. For u, v ∈ X the segment [v, u] lies in the convex set X, so the mean-value inequality together with Certificate 5.1(C2) gives ∥H(u) − H(v)∥∞ ≤ (7.200691 × 10−27 ) ∥u − v∥∞ . Banach’s fixed-point theorem yields a unique u∗ ∈ X with H(u∗ ) = u∗ . 14

The bound (38) also forces R to be invertible: if y T R = 0 for some y ̸= 0, then y T (I − RDG(u)) = y T , so 1 would be an eigenvalue of (I − RDG(u))T , contradicting ρ(I − RDG(u)) ≤ ∥I − RDG(u)∥∞ < 1. Hence H(u∗ ) = u∗ reads RG(u∗ ) = 0 and gives G(u∗ ) = 0. Conversely, every zero of G in X is a fixed point of H, so u∗ is the only such zero. Write the exact components as u∗ = (λ̄∗3 , z∗ , r∗ ),

z∗ = (a∗ , p∗ ).

(47)

Unpacking G(u∗ ) = 0 gives the exact, rather than numerical, identities Fλ̄7 (z∗ ) = z∗ ,

Fλ̄15 (Dλ̄ (r∗ )) = z∗ .

(48)

The third eigenvalue lies in the radius-10−55 interval centered at c1 , so in particular λ̄∗3 = 2.71278149106533722827430219171662596122 . . . . Equation (46) shows Fλ̄ (z∗ ) ̸= z∗ . The least period of z∗ divides 7; since 7 is prime and the period is not 1, it equals 7.

6

Floquet stability and nonlinear attraction

Let A7 := Dz Fλ̄7 (z∗ ).

(49)

The derivative exists because all seven denominators are positive. We now derive every stability implication of Certificate 5.1(C5). Lemma 6.1 (The certified matrices are positive definite). The exact matrices P and Q7 = P − AT7 P A7 satisfy P ≻ 0 and Q7 ≻ 0. Proof. The interval LDLT recursion is performed without pivoting; every division is valid because each preceding pivot interval is positive. It encloses an exact factorization P = LDLT with L unit lower triangular and the four diagonal entries of D bounded below by the positive numbers in (43); hence P ≻ 0. The matrix Q7 is symmetric, and the interval calculation proves strict row diagonal dominance, X (Q7 )ii − |(Q7 )ij | ≥ mi > 0, (50) j̸=i

with the margins mi bounded below by (44); in particular every diagonal entry is positive. By Gershgorin’s circle theorem, every eigenvalue of Q7 is positive, so Q7 ≻ 0. Lemma 6.2 (Schur stability of the seven-step derivative). Every complex eigenvalue of A7 has modulus strictly smaller than one. Proof. Let v ∈ C4 \ {0} and ζ ∈ C satisfy A7 v = ζv. A real symmetric positive-definite matrix defines a positive Hermitian form on C4 , so v ∗ P v > 0 and v ∗ Q7 v > 0. Using that A7 and P are real, v ∗ Q7 v = v ∗ P v − (A7 v)∗ P (A7 v) = (1 − |ζ|2 )v ∗ P v. Both sides are strictly positive, so |ζ| < 1.

15

Lemma 6.3 (Nonlinear attraction in an adapted norm). The point z∗ is a locally attracting fixed point of Φ := Fλ̄7 . Proof. Define ∥x∥P = (xT P x)1/2 . Since P ≻ 0, this is a norm. On the compact P -unit sphere SP = {x : xT P x = 1} the continuous function x 7→ xT Q7 x is strictly positive by Lemma 6.1 and therefore has a positive minimum δ > 0; by homogeneity, xT Q7 x ≥ δxT P x

for every x ∈ R4 .

(51)

Since Q7 = P − AT7 P A7 , (51) gives ∥A7 x∥2P = xT (P − Q7 )x ≤ (1 − δ) ∥x∥2P , and the left side is nonnegative, so 0 < δ ≤ 1 and ∥A7 ∥P ≤ q0 :=

√

1 − δ < 1.

(52)

All normalization denominators at the cycle phases are positive by (39) and remain positive on a small Euclidean neighborhood of each phase, so Φ is C 1 near z∗ with continuous derivative and DΦ(z∗ ) = A7 . Choose any q with q0 < q < 1. Continuity of DΦ in the induced P -norm gives ε > 0 such that ∥DΦ(z)∥P ≤ q whenever ∥z − z∗ ∥P ≤ ε. (53) The closed P -ball in (53) is convex, so for any z in this ball the mean-value inequality gives ∥Φ(z) − z∗ ∥P = ∥Φ(z) − Φ(z∗ )∥P ≤ q ∥z − z∗ ∥P . Thus the ball is forward invariant under Φ, and induction gives ∥Φm (z) − z∗ ∥P ≤ q m ∥z − z∗ ∥P → 0.

7

Translation and an open four-dimensional family

The projective dynamics depends only on differences a − λi , whereas the original BB multiplier also contains the denominator a. A spectral translation exploits this distinction. Lemma 7.1 (Translation conjugacy). For s ∈ R, let λ′ = λ + s1 and Cs (a, p) = (a + s, p). Wherever the maps are defined, Fλ′ ◦ Cs = Cs ◦ Fλ , Dλ′ = Cs ◦ Dλ . (54) Proof. Because

P

i pi = 1,

µλ′ (p) =

X

(λi + s)pi = µλ (p) + s.

i

Also (a + s) − (λi + s) = a − λi , so Zλ′ (a + s, p) =

X

pi ((a + s) − (λi + s))2 = Zλ (a, p),

i

and every normalized weight in (25) is unchanged. Substitution into (26) proves both identities. Take s = 4. The shifted central spectrum is λ∗ = 5, 5.8786699041860466, 6.712781491065337228 . . . , 8 .

(55)

Translation leaves the projective weights, all differences, all Z values, and all state derivatives unchanged; it adds four only to the scalar a. Let z∗+ = C4 (z∗ ). Then Fλ7∗ (z∗+ ) = z∗+ ,

Fλ15∗ (Dλ∗ (r∗ )) = z∗+ . 16

(56)

Proposition 7.2 (Uniform factor bounds in dimension four). There are nonempty open neighborhoods L4 ⊂ {λ ∈ R4 : 0 < λ1 < · · · < λ4 } and P4 ⊂ int ∆3 of λ∗ and r∗ , respectively, such that for every (λ, p0 ) ∈ L4 × P4 , the orbit (ak , pk ) = Fλk (Dλ (p0 )) is defined for all k and |ak − λi | 10−6 ≤ ≤ 0.61 (1 ≤ i ≤ 4, k ≥ 0). (57) ak Proof. We divide the argument into continuation, uniform attraction, finite-time entry, and the numerical factor bounds. Step 1: continuation of the seven-cycle. Define Ψ(λ, z) := Fλ7 (z) − z. All seven Z denominators are positive at the cycle by (39) and, by continuity, positive nearby, so Ψ is C 1 near (λ∗ , z∗+ ). Identity (56) gives Ψ(λ∗ , z∗+ ) = 0, and Dz Ψ(λ∗ , z∗+ ) = A7 − I is invertible because 1 is not an eigenvalue of A7 (Lemma 6.2). The implicit-function theorem gives an open neighborhood Uλ of λ∗ and a C 1 periodic-point function z(λ) satisfying Fλ7 (z(λ)) = z(λ),

z(λ∗ ) = z∗+ .

(58)

Step 2: uniform local attraction. At the central point, P − AT7 P A7 ≻ 0. Positive definiteness is open in the entries: if Q7 ≻ 0, its minimum on the Euclidean unit sphere is some m > 0, and every symmetric perturbation E with ∥E∥2 < m/2 satisfies xT (Q7 + E)x ≥ m/2 for every Euclidean unit x. The matrix Dz Fλ7 (z(λ)) depends continuously on λ, so after shrinking Uλ , the strict inequality P − Dz Fλ7 (z(λ))T P Dz Fλ7 (z(λ)) ≻ 0 holds for every λ ∈ Uλ . Repeating the compact-unit-sphere argument from Lemma 6.3, now over the closure of a sufficiently small parameter ball, gives a common contraction constant q < 1 and a common radius ε > 0 such that the P -ball Bλ := {z : ∥z − z(λ)∥P < ε}

(59)

is mapped into itself by Fλ7 for every parameter in that smaller ball. Shrinking ε if necessary, all seven intermediate images Fλj (Bλ ), 0 ≤ j < 7, stay in neighborhoods where Z > 0 and where every cycle-phase separation remains strict. Step 3: admissible initial conditions enter the trapping tube. At the central parameter, (56) gives exact entry after fifteen steps. The map (λ, r) 7−→ Fλ15 (Dλ (r))

(60)

is continuous near (λ∗ , r∗ ): it is a finite composition of rational maps, and every denominator on the central connection is positive by (39). Since Bλ∗ is an open neighborhood of z∗+ , continuity gives a joint open neighborhood W of (λ∗ , r∗ ) whose image under (60) lies in the corresponding parameter-dependent trapping ball; W contains a product L4 × P4 of open neighborhoods, and (35) and (36) allow P4 to be chosen inside int ∆3 . For any (λ, p0 ) in this product, the first fifteen iterates are defined by finite-time continuity, the fifteenth state lies in Bλ , and thereafter the seven-step iterates remain in Bλ with the intermediate states in its seven-phase trapping tube. Hence the orbit is defined for all time. 17

Step 4: lower factor bound. At the central parameter and initial weight, the finite connection has absolute separation greater than 1.3886903675 × 10−4 by (40); the cycle has the much larger separation (41). All relevant functions are continuous, so we may shrink the trapping tube and L4 × P4 so that |ak − λi | > 10−5 for all i and all k ≥ 0. (61) This is a finite continuity assertion for 0 ≤ k ≤ 15 and a uniform seven-phase assertion in the forward-invariant trapping tube for k ≥ 15. Shrink L4 further so that λ4 < 8.01. Every scalar ak along an admissible orbit is a convex combination of the eigenvalues, by a0 = µλ (p0 ) and ak = µλ (pk−1 ) for k ≥ 1; thus ak ≤ λ4 < 8.01, and with (61), |ak − λi | 10−5 > > 10−6 . ak 8.01 Step 5: upper factor bound. Shrink L4 once more so that λ1 > 4.99,

λ4 < 8.01,

λ4 − λ1 < 3.02.

(62)

For every convex combination ak ∈ [λ1 , λ4 ] and every λi in the same interval, |ak − λi | ≤ λ4 − λ1 < 3.02, so

ak ≥ λ1 > 4.99,

3.02 |ak − λi | < < 0.61. ak 4.99

This completes all parts of (57).

8

Extension to every finite dimension

Fix n > 4. We embed the four active eigenvalues at indices 1, 2, 3, n and choose any distinct central extra eigenvalues 7.49 < ν4 < ν5 < · · · < νn−1 < 7.51. (63) Such a choice exists for every finite n: for example, equally spaced points in the interval (7.49, 7.51) suffice. The resulting central ordered spectrum is λ∗,n = (5, 5.8786699041860466, 6.712781491065337228 . . . , ν4 , . . . , νn−1 , 8).

(64)

Give every extra coordinate zero projective weight. The corresponding face of ∆n−1 is invariant because each update has the form pi (a − λi )2 p′i = ; Z(a, p) if pi = 0, then p′i = 0. Thus the four-dimensional seven-cycle persists as a boundary cycle of the n-dimensional projective map. Use independent local coordinates (a, p1 , p2 , p3 , q4 , . . . , qn−1 ), where qℓ is the weight at νℓ , and eliminate the active weight at the eigenvalue 8: pn = 1 − p 1 − p 2 − p 3 −

n−1 X ℓ=4

18

qℓ .

At a boundary-cycle phase, q = 0 and Zt =

X

pt,j (at − λj )2 > 0.

j∈{1,2,3,n}

For an extra weight, qℓ′ =

qℓ (a − νℓ )2 . Z(a, p, q)

(65)

Differentiate (65) at q = 0 by the quotient rule gives (a − νℓ )2 Z − qℓ (a − νℓ )2 (∂Z/∂qℓ ) (a − νℓ )2 ∂qℓ′ = = , ∂qℓ q=0 Z2 Z q=0 qℓ (a − νℓ )2 (∂Z/∂qm ) ∂qℓ′ =− =0 ∂qm q=0 Z2 q=0

(m ̸= ℓ).

Every derivative of qℓ′ with respect to a base variable a, p1 , p2 , p3 also contains the prefactor qℓ and is zero at q = 0. Thus the one-step derivative at phase t is block upper triangular: !

Bt Ct , 0 diag(γt (ν4 ), . . . , γt (νn−1 ))

γt (ν) =

(at − ν)2 . Zt

(66)

Products of block upper-triangular matrices are block upper triangular, and their diagonal blocks are the products of the diagonal blocks. Therefore the full seven-step derivative Mn has the form !

Mn =

A7 ∗ , 0 diag(τ (ν4 ), . . . , τ (νn−1 ))

τ (ν) :=

6 Y (at − ν)2 t=0

Zt

.

(67)

The certificate evaluates the unshifted interval ν̄ ∈ [3.49, 3.51]. Translation by four leaves every difference and Zt unchanged, so (45) holds for every ν ∈ [7.49, 7.51], including every extra eigenvalue in (63). The characteristic polynomial of a block upper-triangular matrix is the product of the characteristic polynomials of its diagonal blocks. Thus, the eigenvalues of Mn are those of A7 together with the scalars τ (νℓ ). Lemma 6.2 and (45) imply ρ(Mn ) < 1. Since ρ(Mn ) < 1, we have Mnk

2

(68)

≤ Cγ k for any γ ∈ (ρ(Mn ), 1) and some finite C, so the series Pen :=

∞ X

(MnT )k Mnk

(69)

k=0

converges, satisfies Pen ⪰ I ≻ 0, and obeys the exact identity Pen − MnT Pen Mn = I. It therefore supplies an adapted Lyapunov norm in the full state space, even though the four-dimensional matrix P does not include the new coordinates, and thus the nonlinear attraction proof of Lemma 6.3 applies verbatim in the full state space. The boundary cycle lies on a simplex face, but the rational map is smooth on an ordinary ambient Euclidean neighborhood because every cycle denominator is positive, so its attracting neighborhood intersects the simplex in a relative neighborhood containing interior points. To see this explicitly, start from the boundary initialization r∗ , assign sufficiently small positive values 19

to all extra weights, and subtract their sum from the active weight r∗,4 at eigenvalue 8, which is possible because r∗,4 > 0.0108 by (36). The resulting weight lies in int ∆n−1 and can be arbitrarily close to the boundary weight. The boundary admissible state lands on the cycle after fifteen steps; by continuity of the finite iterate, every sufficiently close interior weight lands in the attracting (n) neighborhood. Choose one such interior weight and call it r∗ . We must also continue the cycle under perturbations of all n eigenvalues. The map Ψn (λ, z) = Fλ7 (z) − z is C 1 near the central boundary cycle because the denominators are positive, and its state derivative Mn − I is invertible by (68). The implicit-function theorem therefore continues the periodic point for every sufficiently small perturbation of all n eigenvalues, and the strict Lyapunov inequality, the finite-time entry property, and all separation inequalities persist after the parameter and initial neighborhoods are shrunk. The certificate’s separation test includes every ν̄ ∈ [3.49, 3.51], hence after translation every ν ∈ [7.49, 7.51]. The factor-bound calculation in Steps 4–5 of Proposition 7.2 therefore applies to active and extra coordinates alike. We have proved the following. Proposition 8.1 (Uniform factor bounds for every finite n). For every finite n ≥ 4, there are nonempty open sets Ln ⊂ {λ ∈ Rn : 0 < λ1 < · · · < λn },

Pn ⊂ int ∆n−1 ,

(70)

such that every admissible projective orbit initialized by (λ, p0 ) ∈ Ln × Pn is defined for all time and satisfies |ak − λi | 10−6 ≤ ≤ 0.61 (1 ≤ i ≤ n, k ≥ 0). (71) ak There is no hidden uniformity claim as n → ∞: for each fixed finite n, finitely many distinct extra eigenvalues are selected and then a possibly n-dependent open neighborhood is taken. The rate constants 10−6 and 0.61, however, are the same for every finite n.

9

Lifting the projective construction to BB problems and measure

We first obtain componentwise gradient bounds for a diagonal quadratic. Fix λ ∈ Ln and choose any p0 ∈ Pn , any radius r > 0, and any sign vector σ ∈ {±1}n . Set √ di0 := σi r p0,i . Then

(72)

i 2 2 i 2 P (dj )2 = p . Conversely, on any fixed open orthant, the map 0,i i (d0 ) = r and (d0 ) / j 0

P

Θσ : (0, ∞) × int ∆n−1 → {d ∈ Rn : sign(di ) = σi },

√ Θσ (r, p)i = σi r pi ,

(73)

has the smooth inverse r = ∥d∥2 , pi = d2i / ∥d∥22 , and is therefore a diffeomorphism. In particular, the image of (1, 2) × Pn is a nonempty open set of initial spectral gradients. Combining Proposition 8.1 with Lemma 4.2 gives, coordinate by coordinate, di0 10−6k ≤ dik = di0

k−1 Y

|at − λi | ≤ di0 0.61k . a t t=0

This is the spectral version of (6). Because Pn ⊂ int ∆n−1 , every di0 is nonzero. The lower bound then makes every dik , and hence every gk , nonzero. Lemma 4.1 shows inductively that every BB denominator sTk−1 yk−1 = sTk−1 Ask−1 is strictly positive. Thus the BB1 iteration is well-defined for 20

all k, rather than merely being a formal solution of the scalar recurrence. The upper bound tends to zero, so gk → 0 and, because A is invertible, xk − x∗ = A−1 gk → 0. We next allow arbitrary eigenvectors and prove openness in the full matrix space. Define An := {A ∈ Sn++ : λ1 (A) < · · · < λn (A), (λ1 (A), . . . , λn (A)) ∈ Ln }.

(74)

This set is nonempty: it contains the diagonal matrix with any spectrum in Ln . It is open because the ordered eigenvalues of a real symmetric matrix are continuous functions of the matrix entries (Weyl’s inequality gives |λi (A) − λi (B)| ≤ ∥A − B∥2 ) and Ln is open. For A ∈ An and g ̸= 0, define the orientation-free spectral weights pi (A, g) :=

g T Πi (A)g ∥Πi (A)g∥2 = . gT g ∥g∥2

(75)

They are nonnegative and sum to one because the spectral projectors are orthogonal and sum to I. The map (A, g) 7→ p(A, g) is continuous on simple-spectrum matrices with g ̸= 0: on the simple-spectrum set, Y A − λj (A)I . (76) Πi (A) = λ (A) − λj (A) j̸=i i To verify (76), apply the polynomial on its right to an eigenvector of A: it equals one on the ith eigenspace and zero on every other eigenspace. Every denominator is nonzero because the spectrum is simple. Therefore, continuity of the eigenvalues makes the right side locally continuous in A, and (75) is a quotient of continuous functions with positive denominator g T g, proving the claim. Define the open set of matrix–gradient pairs Γn := {(A, g) : A ∈ An , 1 < ∥g∥ < 2, p(A, g) ∈ Pn }.

(77)

It is open by the continuity just established. It is nonempty: take a diagonal A with spectrum in √ Ln , any p ∈ Pn , a sign vector σ, and gi = (3/2)σi pi ; then ∥g∥ = 3/2 and p(A, g) = p. For the fixed vector b, the map Jb : Sn++ × Rn → Sn++ × Rn ,

Jb (A, x) = (A, Ax − b),

(78)

is a smooth bijection with smooth inverse Jb−1 (A, g) = (A, A−1 (g + b)), matrix inversion being smooth on Sn++ ; it is therefore a diffeomorphism. Set Ωn := Jb−1 (Γn ),

Xn (A) := {x : (A, x) ∈ Ωn }.

(79)

The set Ωn is nonempty and open, and every fiber Xn (A) is nonempty and open: for a fixed A ∈ An , one can choose a spectral-coordinate gradient from (72) and then put x = A−1 (g + b). Every nonempty open subset of a finite-dimensional Euclidean space has positive Lebesgue measure; applied in dimensions n(n+1)/2+n and n, this proves the positive-measure assertions of Theorem 2.1. It remains to prove the norm, objective, and non-superlinear conclusions. Because Q is orthogonal, (6) gives 10−12k ∥g0 ∥2 ≤ ∥gk ∥2 =

X

dik

2

≤ 0.612k ∥g0 ∥2 ,

i

and taking nonnegative square roots proves (8). Since gk = Aek , one has, in the eigenbasis, ∥ek ∥2A = gkT A−1 gk =

2 n X di k

i=1

21

λi

,

(80)

and applying the component bounds term by term proves (9). Completing the square in (3) gives f (x)−f (x∗ ) = 12 ∥x − x∗ ∥2A , so squaring (9) and multiplying by 1/2 proves (10). The lower bounds in (11) follow by taking kth roots, for example ∥gk ∥1/k ≥ 10−6 ∥g0 ∥1/k → 10−6 . Finally, Q-superlinear 1/k convergence of a positive null sequence (uk+1 /uk → 0) implies uk → 0, so (11) rules out both root-superlinear and Q-superlinear convergence. All assertions of Theorem 2.1 are proved.

10

Scope and limitations

The theorem establishes a robust lower-geometric rate on a constructed open family. Four qualifications are important. 1. It is a result for BB1, not BB2. The initial step is the explicit inverse gradient Rayleigh quotient (4). An arbitrary first step has a related projective description but is not the admissible-initialization map certified here. 2. The conclusion holds on a nonempty open, positive-measure family of spectra and initial points. It is not asserted for every positive spectrum or for almost every spectrum. 3. The computer-assisted step proves existence of one exact orbit inside a tiny rational box. The displayed decimal for λ∗3 is a locator for that exact root, not a claim that the eigenvalue is itself a terminating decimal or has a closed form. 4. The constants 10−6 and 0.61 are safety margins. The proof uses a minimum transient separation of approximately 1.39 × 10−4 and an upper multiplier bound close to 0.6 at the central translated spectrum. No optimality is claimed.

References Hirotugu Akaike. On a successive transformation of probability distribution and its application to the analysis of the optimum gradient method. Annals of the Institute of Statistical Mathematics, 11(1):1–16, 1959. Roberto Andreani, Ernesto G Birgin, José Mario Martínez, and Jinyun Yuan. Spectral projected gradient and variable metric methods for optimization with linear inequalities. IMA Journal of Numerical Analysis, 25(2):221–252, 2005. Behzad Azmi and Karl Kunisch. Analysis of the barzilai-borwein step-sizes for problems in hilbert spaces. Journal of Optimization Theory and Applications, 185(3):819–844, 2020. Behzad Azmi and Karl Kunisch. On the convergence and mesh-independent property of the barzilai– borwein method for pde-constrained optimization. IMA Journal of Numerical Analysis, 42(4): 2984–3021, 2022. Jonathan Barzilai and Jonathan M Borwein. Two-point step size gradient methods. IMA journal of numerical analysis, 8(1):141–148, 1988. Ernesto G Birgin, José Mario Martínez, and Marcos Raydan. Nonmonotone spectral projected gradient methods on convex sets. SIAM Journal on Optimization, 10(4):1196–1211, 2000.

22

Ernesto G Birgin, José Mario Martínez, and Marcos Raydan. Algorithm 813: Spg—software for convex-constrained optimization. ACM Transactions on Mathematical Software (TOMS), 27(3): 340–349, 2001. Ernesto G Birgin, Jose Mario Martínez, and Marcos Raydan. Spectral projected gradient methods: review and perspectives. Journal of Statistical Software, 60:1–21, 2014. Silvia Bonettini, Riccardo Zanella, and Luca Zanni. A scaled gradient projection method for constrained image deblurring. Inverse problems, 25(1):015002, 2009. Oleg Burdakov, Yuhong Dai, and Na Huang. Stabilized barzilai-borwein method. Journal of Computational Mathematics, pages 916–936, 2019. Yu-Hong Dai. Alternate step gradient method. Optimization, 52(4-5):395–415, 2003. Yu-Hong Dai. A new analysis on the barzilai-borwein gradient method. Journal of the operations Research Society of China, 1(2):187–198, 2013. Yu-Hong Dai and Roger Fletcher. On the asymptotic behaviour of some new gradient methods. Mathematical Programming, 103(3):541–559, 2005a. Yu-Hong Dai and Roger Fletcher. Projected barzilai-borwein methods for large-scale box-constrained quadratic programming. Numerische Mathematik, 100(1):21–47, 2005b. Yu-Hong Dai and Li-Zhi Liao. R-linear convergence of the barzilai and borwein gradient method. IMA Journal of numerical analysis, 22(1):1–10, 2002. Yu-Hong Dai, William W Hager, Klaus Schittkowski, and Hongchao Zhang. The cyclic barzilai— borwein method for unconstrained optimization. IMA Journal of Numerical Analysis, 26(3): 604–627, 2006. Yu-Hong Dai, Yakui Huang, and Xin-Wei Liu. A family of spectral gradient methods for optimization. arXiv preprint arXiv:1812.02974, 2018. Mário AT Figueiredo, Robert D Nowak, and Stephen J Wright. Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems. IEEE Journal of selected topics in signal processing, 1(4):586–597, 2007. Roger Fletcher. On the barzilai-borwein method. In Optimization and control with applications, pages 235–256. Springer, 2005. George E Forsythe. On the asymptotic directions of the s-dimensional optimum gradient method. Numerische Mathematik, 11(1):57–76, 1968. Ana Friedlander, José Mario Martínez, Brigida Molina, and Marcus Raydan. Gradient method with retards and generalizations. SIAM Journal on Numerical Analysis, 36(1):275–289, 1998. Juan Gao, Xin-Wei Liu, Yu-Hong Dai, Yakui Huang, and Peng Yang. Achieving geometric convergence for distributed optimization with barzilai-borwein step sizes. Science China. Information Sciences, 65(4):149204, 2022. Luigi Grippo and Marco Sciandrone. Nonmonotone globalization techniques for the barzilai-borwein gradient method. Computational Optimization and Applications, 23(2):143–169, 2002. 23

Chang He, Wenzhi Gao, Bo Jiang, Madeleine Udell, and Shuzhong Zhang. New results on the polyak stepsize: Tight convergence analysis and universal function classes. arXiv preprint arXiv:2512.06231, 2025. Ya-Kui Huang, Yu-Hong Dai, and Xin-Wei Liu. Equipping the barzilai–borwein method with the two dimensional quadratic termination property. SIAM Journal on Optimization, 31(4):3068–3096, 2021. Yakui Huang, Hongwei Liu, and Shuisheng Zhou. Quadratic regularization projected barzilai– borwein method for nonnegative matrix factorization. Data mining and knowledge discovery, 29 (6):1665–1684, 2015. Yakui Huang, Yu-Hong Dai, Xin-Wei Liu, and Hongchao Zhang. On the asymptotic convergence and acceleration of gradient methods. Journal of Scientific Computing, 90(1):7, 2022. Bruno Iannazzo and Margherita Porcelli. The riemannian barzilai–borwein method with nonmonotone line search and the matrix geometric mean computation. IMA Journal of Numerical Analysis, 38(1):495–517, 2018. Rudolf Krawczyk. Newton-algorithmen zur bestimmung von nullstellen mit fehlerschranken. Computing, 4(3):187–201, 1969. Dawei Li and Ruoyu Sun. On a faster R-linear convergence rate of the barzilai-borwein method, 2021. Luc Pronzato, Henry P Wynn, and Anatoly A Zhigljavsky. Asymptotic behaviour of a family of gradient algorithms in Rd and Hilbert spaces. Mathematical programming, 107(3):409–438, 2006. Marcos Raydan. On the barzilai and borwein choice of steplength for the gradient method. IMA Journal of Numerical Analysis, 13(3):321–326, 1993. Marcos Raydan. The barzilai and borwein gradient method for the large scale unconstrained minimization problem. SIAM Journal on Optimization, 7(1):26–33, 1997. Siegfried M Rump. Verification methods: Rigorous results using floating-point arithmetic. In Proceedings of the 2010 International Symposium on Symbolic and Algebraic Computation, pages 3–4, 2010. Thomas Serafini, Gaetano Zanghirati, and Luca Zanni. Gradient projection methods for quadratic programs and applications in training support vector machines. Optimization Methods and Software, 20(2-3):353–378, 2005. Conghui Tan, Shiqian Ma, Yu-Hong Dai, and Yuqiu Qian. Barzilai-borwein step size for stochastic gradient descent. volume 29, 2016. Ewout Van Den Berg and Michael P Friedlander. Probing the pareto frontier for basis pursuit solutions. Siam journal on scientific computing, 31(2):890–912, 2009. Yanfei Wang and Shiqian Ma. Projected barzilai–borwein method for large-scale nonnegative image restoration. Inverse Problems in Science and Engineering, 15(6):559–583, 2007. Zaiwen Wen and Wotao Yin. A feasible method for optimization with orthogonality constraints. Mathematical Programming, 142(1):397–434, 2013. 24

Zaiwen Wen, Wotao Yin, Donald Goldfarb, and Yin Zhang. A fast algorithm for sparse reconstruction based on shrinkage, subspace optimization, and continuation. SIAM Journal on Scientific Computing, 32(4):1832–1857, 2010. Stephen J Wright, Robert D Nowak, and Mário AT Figueiredo. Sparse reconstruction by separable approximation. IEEE Transactions on signal processing, 57(7):2479–2493, 2009. Ya-xiang Yuan. A new stepsize for the steepest descent method. Journal of Computational Mathematics, pages 149–156, 2006. Hongchao Zhang and William W Hager. A nonmonotone line search technique and its application to unconstrained optimization. SIAM journal on Optimization, 14(4):1043–1056, 2004. Danqing Zhou, Shiqian Ma, and Junfeng Yang. Adabb: Adaptive barzilai-borwein method for convex optimization. Mathematics of Operations Research, 51(1):715–745, 2026.

25

A

Computer-Assisted Certificate

The computer-assisted component of the proof is implemented by the following Python script. It uses only Python’s standard decimal module. 1 from

decimal import Decimal, getcontext, localcontext, ROUND_FLOOR, ROUND_CEILING

2 3 getcontext().prec 4 D

= Decimal

= 120

5 6 def

opdown(fn): with localcontext() as c: 8 c.prec=120; c.rounding=ROUND_FLOOR 9 return +fn() 10 def opup(fn): 11 with localcontext() as c: 12 c.prec=120; c.rounding=ROUND_CEILING 13 return +fn() 7

14 15 class 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

IV: __slots__=('lo','hi') def __init__(self,lo,hi=None): self.lo=D(str(lo)) if not isinstance(lo,D) else lo self.hi=self.lo if hi is None else (D(str(hi)) if not isinstance(hi,D) else hi) assert self.lo <= self.hi def __add__(a,b): b=asiv(b); return IV(opdown(lambda:a.lo+b.lo),opup(lambda:a.hi+b.hi)) __radd__=__add__ def __neg__(a): return IV(-a.hi,-a.lo) def __sub__(a,b): return a+(-asiv(b)) def __rsub__(a,b): return asiv(b)+(-a) def __mul__(a,b): b=asiv(b) los=[opdown(lambda x=x,y=y:x*y) for x in (a.lo,a.hi) for y in (b.lo,b.hi)] his=[opup(lambda x=x,y=y:x*y) for x in (a.lo,a.hi) for y in (b.lo,b.hi)] return IV(min(los),max(his)) __rmul__=__mul__ def recip(a): assert not (a.lo <= 0 <= a.hi) valslo=[opdown(lambda x=x:D(1)/x) for x in (a.lo,a.hi)] valshi=[opup(lambda x=x:D(1)/x) for x in (a.lo,a.hi)] return IV(min(valslo),max(valshi)) def __truediv__(a,b): return a*asiv(b).recip() def __rtruediv__(a,b): return asiv(b)*a.recip() def __repr__(a): return f'[{a.lo},{a.hi}]'

41 42 def 43 def 44 45 46

asiv(x): return x if isinstance(x,IV) else IV(x) sumup(vals): s=D(0) for v in vals: s=opup(lambda s=s,v=v:s+v) return s

47 48 class 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63

AD: __slots__=('v','d') def __init__(self,v,d=None): self.v=asiv(v); self.d=[IV(0) for _ in range(8)] if d is None else d def __add__(a,b): b=asad(b); return AD(a.v+b.v,[x+y for x,y in zip(a.d,b.d)]) __radd__=__add__ def __neg__(a): return AD(-a.v,[-x for x in a.d]) def __sub__(a,b): return a+(-asad(b)) def __rsub__(a,b): return asad(b)+(-a) def __mul__(a,b): b=asad(b); return AD(a.v*b.v,[x*b.v+a.v*y for x,y in zip(a.d,b.d)]) __rmul__=__mul__ def recip(a): return AD(a.v.recip(),[-x/(a.v*a.v) for x in a.d]) def __truediv__(a,b): return a*asad(b).recip() def __rtruediv__(a,b): return asad(b)*a.recip()

64

26

65 def 66 def 67

asad(x): return x if isinstance(x,AD) else AD(x) var(box,j): der=[IV(0) for _ in range(8)];der[j]=IV(1);return AD(box,der)

68 69 L2=AD('1.8786699041860466') 70 ONE=AD(1);

FOUR=AD(4)

71 72 def 73 74 75 76 77

fmap_with_den(z,la): m=z[0]; p=[z[1],z[2],z[3],ONE-z[1]-z[2]-z[3]] terms=[p[i]*(m-la[i])*(m-la[i]) for i in range(4)] den=sum(terms,AD(0)) mean=sum((la[i]*p[i] for i in range(4)),AD(0)) return [mean]+[terms[i]/den for i in range(3)],den

78 79 def 80 81

fmap(z,la): out,_=fmap_with_den(z,la) return out

82 83 def 84 85

fpow(z,la,n): for _ in range(n): z=fmap(z,la) return z

86 87 def 88 89 90 91 92 93 94 95 96

f_and_j(boxes): x=[var(boxes[i],i) for i in range(8)] la=[ONE,L2,x[0],FOUR] z=x[1:5] a=fpow(z,la,7) p=x[5:8]; p4=ONE-sum(p,AD(0)) m0=sum((la[i]*([p[0],p[1],p[2],p4][i]) for i in range(4)),AD(0)) b=fpow([m0]+p,la,15) out=[a[i]-z[i] for i in range(4)]+[b[i]-z[i] for i in range(4)] return [o.v for o in out],[o.d for o in out]

97 98 C=[D(s)

for s in [

99 '2.71278149106533722827430219171662596122131200973311943581550524557629575811',

100 '1.33344967252388129918980034968384933800281196626470897628663867165111253656',

101 '0.828108598362557830284706279933642675447754907755381063536594022726890227454', 102 '0.167043051370226631389809462759502319781038286313814177038878187525592790929', 103 '0.004717957350367412952812533939501994168806798103523267952088050658941184644',

104 '0.0000134113750442893933116720965327494222816004773939806475307653666721805974245', 105 '0.989171705769078149047269711515135535847128223427916068474163743526926233731',

106 '0.00000349260028439145589133467903469364875523517774851271154130219629904643379924']] 107 108 RSTR='''

109 -0.810526348848592888901084652460863262098889862944445915997922,307.24267338133438527153922854692102527957829 ⌋

4613703728839309,307.741025629901050999029082768910651032865241495000968766894,312.7827425715368916260963 ⌋ 23271366901386778871589158283516632,-0.0137368376192333684238175487126343869576306541180311710900804,-421 ⌋ ,→ .401532972233614339176178157325006069741606160158648367929,-420.72062833282751211095253495839339095174156 ⌋ ,→ 8565011272813875,-417.510650350147764867647073770427680507322497423363849322798 110 -2.53381240743594641550513665431248219325091959956406240279772,219.146844709590120107670663557533796892315501 ⌋ ,→ 194180693740921,221.314903815662600782780829104204280970210726156970789478092,238.43306113888209660651090 ⌋ ,→ 7554404630970281465079411306474734,-0.0066156137536975764499876413028288031850150714975420122664002,-202. ⌋ ,→ 94552899548185110702795449214798069667985467445722312406,-202.6176076629107747123467952188369797982499185 ⌋ ,→ 28977551841755,-201.071693306206914696638563015692889394595107198327734687951 111 1.21923182241183846078543106866408237572050290347089038684799,-283.490269055458601705948958623004565413293682 ⌋ ,→ 800367378712469,-284.215637538925680074135617356259178167657374408793138161663,-294.814697516715266368336 ⌋ ,→ 657507713365927445251871977984552338,0.00510430037289669853951692869951609336796693289217752688876014,156 ⌋ ,→ .583346292001534198852820229253143092303058113580534397193,156.330337418986871290621900890772993639761529 ⌋ ,→ 519065013423077,155.137581686689234202670907269428105368197596348629742204233 112 -1.26108044272173568710731141145266481650790341533600008992776,313.942924362985217723054543998906451668776939 ⌋ ,→ 851163166266309,314.691544869619096433878614662418685401023724831284472124436,326.12764070745284415061672 ⌋ ,→ 026001723134524840696137617430052,-0.00557213932631546351238737770698205496104361101774514176877326,-170. ⌋ ,→ 935124890502247958801827793221206576836705555576496966117,-170.658926275954081329267222464116649321236505 ⌋ ,→ 862472492490618,-169.356847511562687722423100405437534964897323010662551564416 113 0.043078195604279001914097963327737890401646235452894255110503,-31.297922833392128192888211392587198270527115 ⌋ ,→ 7312435670381131,-31.3219714710519066420310359824923010174212504037635542365452,-32.166217200933550156209 ⌋ ,→ 2089499137302982952674704895954239567,0.000478099922934842210827835744935335196661930745438524440498953,1 ⌋ ,→ 4.6665517947567543301073202404399084503649597758063756119065,14.64285343967317710227853739404576405908733 ⌋ ,→ 5950989354927027,14.531132658756835739832753728374049724120361151869572985717 ,→ ,→

27

114 0.000226520043874099667313273699369756542509219592106852411297113,-0.0225135215590217304102381601094316680754 ⌋

744631503520593497794,-0.0226850497610954124630365608985981365370252310505830687908151,-0.025225307271391 ⌋ 448839622398185040987284482193153534946168519,-0.00006464162903863897657234459873986315211483969963726214 ⌋ ,→ 61687481,-0.0407002504446717796356496112619123261021708978107613923529244,-0.0408037515085730062610317248 ⌋ ,→ 805338145956192819108221297272651,-0.0395655532419233021327540805034448469654922061236616605213023 115 -0.0591009778675161000065292747250136710540662107925102493119043,4.525428637470154004681198684575136861825405 ⌋ ,→ 10968085333767212,4.5708376271796282978337702766393546845959286011088156970192,5.252779805830252696321420 ⌋ ,→ 81451745114823486084850385324393674,0.0182267787018476751330387775829333443438138046050174382552172,13.72 ⌋ ,→ 94182492651964487088572696615033597888320575657007499638,13.754769271794975636070159437239004779426662361 ⌋ ,→ 8339207810693,13.39537325223256904440264858605473726414587831319750699597 116 0.0000756712444366593089630182277257919839520913409911366332434105,-0.008198976004242307532640264918658378766 ⌋ ,→ 18103810542438832290775,-0.00825475135727810105199639024201679021048816336552888167047913,-0.009157358086 ⌋ ,→ 25744896227380169097910909397341703431881848992986,-0.000019729063909076274555279269010657568336641436089 ⌋ ,→ 8273244501419,-0.0159791300873645197383620839265853550901431071050090962014869,-0.01601853297980077370645 ⌋ ,→ 68268873709203041685752673075101630795,-0.0155338627526389510516621091918715132796619616194482426763043''' 117 R=[[IV(s) for s in line.split(',')] for line in RSTR.strip().splitlines()] ,→ ,→

118 119 PSTR='''

120 977656.6634223271235507758960608141,227877753.5872023336353569232362594,226076347.85889685492460652333685333, ⌋

213421596.35978306381023106533985262

,→

121 227877753.5872023336353569232362594,53180721626.999672123361693796538579,52760940596.820424133577893141389299 ⌋

,49812202173.83420304678138035243017

,→

122 226076347.85889685492460652333685333,52760940596.820424133577893141389299,52344478955.77919231215377678997530 ⌋

3,49419059318.10867586670666295706405

,→

123 213421596.35978306381023106533985262,49812202173.83420304678138035243017,49419059318.10867586670666295706405, ⌋

46657455589.35772965711287496916893''' for s in line.split(',')] for line in PSTR.strip().splitlines()]

,→

124 P=[[IV(s) 125 126 def

mm(A,B): return [[sum((A[i][k]*B[k][j] for k in range(len(B))),IV(0)) for j in range(len(B[0]))] for i in ,→ range(len(A))] 128 def mv(A,v): return [sum((A[i][j]*v[j] for j in range(len(v))),IV(0)) for i in range(len(A))] 129 def tr(A): return [list(x) for x in zip(*A)] 127

130 131 def 132 133 134 135 136 137 138 139 140

ldl_pivots(A): n=len(A); L=[[IV(0) for _ in range(n)] for _ in range(n)]; piv=[] for i in range(n): L[i][i]=IV(1) di=A[i][i]-sum((L[i][k]*L[i][k]*piv[k] for k in range(i)),IV(0)) piv.append(di) for j in range(i+1,n): num=A[j][i]-sum((L[j][k]*L[i][k]*piv[k] for k in range(i)),IV(0)) L[j][i]=num/di return piv

141 142 if 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164

__name__=='__main__': rad=D('1e-55') boxes=[IV(c-rad,c+rad) for c in C] _,J=f_and_j(boxes) fx,_=f_and_j([IV(c) for c in C]) rf=mv(R,fx) RJ=mm(R,J) B=[[IV(1 if i==j else 0)-RJ[i][j] for j in range(8)] for i in range(8)] dx=[IV(-rad,rad) for _ in range(8)] bd=mv(B,dx) kc=[IV(C[i])-rf[i]+bd[i] for i in range(8)] print('Krawczyk ratios |K-c|/radius:') for i in range(8): assert C[i]-rad < kc[i].lo and kc[i].hi < C[i]+rad ratio=max(abs(kc[i].lo-C[i]),abs(kc[i].hi-C[i]))/rad print(i,ratio) bnorm=max(sumup(max(abs(v.lo),abs(v.hi)) for v in row) for row in B) assert bnorm < 1 print('contraction infinity row-sum bound:',bnorm) # Period derivative A = derivative of first four outputs wrt z plus I. A=[[J[i][1+j]+(1 if i==j else 0) for j in range(4)] for i in range(4)] print('A intervals widths / centers:') for row in A: print(' '.join(str((v.lo+v.hi)/2) for v in row))

28

print('P LDL lower pivots:') ppiv=ldl_pivots(P); assert all(x.lo>0 for x in ppiv) print(*(x.lo for x in ppiv)) Q=[[P[i][j]-mm(mm(tr(A),P),A)[i][j] for j in range(4)] for i in range(4)] print('Q diagonal-dominance lower margins:') for i in range(4): off=sumup(max(abs(Q[i][j].lo),abs(Q[i][j].hi)) for j in range(4) if j!=i) margin=opdown(lambda:Q[i][i].lo-off); assert margin>0 print(margin) # Uniform transverse Floquet multiplier for any inserted eigenvalue nu in [3.49,3.51]. xv=[var(boxes[i],i) for i in range(8)] la=[ONE,L2,xv[0],FOUR]; zz=xv[1:5]; nu=AD(IV('3.49','3.51')); tau=AD(1)

165 166 167 168 169 170 171 172 173 174 175 176 177

cycle_state=xv[1:5]; cycle_z_min=D('Infinity') for _ in range(7): cycle_state,den=fmap_with_den(cycle_state,la) cycle_z_min=min(cycle_z_min,den.v.lo) p=xv[5:8]; p4=ONE-sum(p,AD(0)) conn_state=[sum((la[i]*[p[0],p[1],p[2],p4][i] for i in range(4)),AD(0))]+p conn_z_min=D('Infinity') for _ in range(15): conn_state,den=fmap_with_den(conn_state,la) conn_z_min=min(conn_z_min,den.v.lo) print('minimum cycle normalization denominator:',cycle_z_min) print('minimum connection normalization denominator:',conn_z_min) assert cycle_z_min>D('0.1516357118221540') assert conn_z_min>D('0.0009510722601452239')

178 179 180 181 182 183 184 185 186 187 188 189 190 191 192

cycle_next,_=fmap_with_den(xv[1:5],la) phase_gap=(xv[1]-cycle_next[0]).v print('first cycle phase difference a0-a1:',phase_gap) assert phase_gap.lo>D('0.1782019618055700634')

193 194 195 196 197

for _ in range(7): m=zz[0]; pp=[zz[1],zz[2],zz[3],ONE-zz[1]-zz[2]-zz[3]] den=sum((pp[i]*(m-la[i])*(m-la[i]) for i in range(4)),AD(0)) tau=tau*(m-nu)*(m-nu)/den zz=fmap(zz,la) print('transverse tau [3.49,3.51]:',tau.v) assert D(0)<tau.v.lo and tau.v.hi<1 # No multiplier vanishes during the exceptional initial segment or on the cycle. def sep(a,b): a=a.v if isinstance(a,AD) else asiv(a); b=b.v if isinstance(b,AD) else asiv(b) if a.hi < b.lo: return opdown(lambda:b.lo-a.hi) if b.hi < a.lo: return opdown(lambda:a.lo-b.hi) return D(0) xx=[sum((la[i]*[p[0],p[1],p[2],p4][i] for i in range(4)),AD(0))]+p transient_sep=D('Infinity') for _ in range(16): transient_sep=min(transient_sep,*(sep(xx[0],ell) for ell in la),sep(xx[0],nu)) xx=fmap(xx,la) cycle_sep=D('Infinity'); xx=xv[1:5] for _ in range(7): cycle_sep=min(cycle_sep,*(sep(xx[0],ell) for ell in la),sep(xx[0],nu)) xx=fmap(xx,la) print('minimum initial-segment m-lambda separation:',transient_sep) print('minimum cycle separation (including nu box):',cycle_sep) assert transient_sep>0 and cycle_sep>0

198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222

B

Prompt Used for Proof Verification

The following prompt was submitted to GPT 5.6 Sol. Current task statement

29

For $n\geq4$, given $n$ positive real numbers $\lambda_1\leq ... \leq \lambda_n$. Consider the ,→ sequence $\{d^1_k, ..., d^n_k\} (k\geq 0)$ where $d^i_k$ are real numbers for $1\leq i\leq n$, ,→ satisfying $$\begin{equation} d^i_{k+1}=d^i_k\cdot\left( \frac{ \sum_{j=1}^n(\lambda_j-\lambda_i)(d^j_{k-1})^2 } { \sum_{j=1}^n\lambda_j(d^j_{k-1})^2 } \right), ~~i=1, \cdots, n. \end{equation}$$ Specifically, in the first iteration, the relation becomes $$\begin{equation} d^i_1=d^i_0\cdot\left(\frac{\sum_{j=1}^n(\lambda_j-\lambda_i)(d^j_0)^2} ,→ {\sum_{j=1}^n\lambda_j(d^j_0)^2}\right), ~~i=1, \cdots, n. \end{equation}$$ Now, given $2n$ real numbers $\lambda_1, ... \lambda_n, d_0^1, ..., d_0^n$, characterize the ,→ decreasing rate of the sequence $|d_k^i|$. In particular, prove or disprove the following: For a positive measure of $(\lambda_1, ..., \lambda_n)\in\mathbb{R}^n$ and a positive measure of ,→ $(d_0^1, ..., d_0^n)\in\mathbb{R}^n$, $|d_k^i|$ converges at a linear rate, i.e., there exists ,→ positive constants $C_i$ and $0<\rho_i<1$ only related to the $2n$ initial real numbers, such ,→ that $|d_k^i|\geq C_i\rho_i^k$. Assume for purposes of this task that a complete affirmative proof or counterexample exists. ,→ Partial progress does not count unless it implies exactly the resolution above. In particular, ,→ proofs for zero measure of $(\lambda_1, ..., \lambda_n)\in\mathbb{R}^n$ or zero measure of ,→ $(d_0^1, ..., d_0^n)\in\mathbb{R}^n$, $|d_k^i|$, constructions of counter-examples with zero ,→ measure of $(\lambda_1, ..., \lambda_n)\in\mathbb{R}^n$ or zero measure of $(d_0^1, ..., ,→ d_0^n)\in\mathbb{R}^n$, $|d_k^i|$, reductions to another unproved conjecture, computational verification through any fixed $n$, and candidate counterexamples without a complete certificate are insufficient. Use multiagent v2 aggressively and dynamically. You have up to 64 concurrent agents available. Do ,→ not use a fixed assignment such as "N agents for strategy X." Instead, manage the search using ,→ the following heuristics: - Begin with a genuinely diverse portfolio of approaches. Agents should explore substantially ,→ different formulations, invariants, reductions, algebraic viewpoints, structural inductions, ,→ decompositions, flow formulations, transition systems, embeddings, extremal arguments, and ,→ computational sanity checks. - Do not tell most agents the currently favored approach. Preserve independence during early rounds ,→ so that agents do not all converge to the same attractive but incomplete reduction. - Maintain an explicit registry of approach families. Group agents by the mathematical idea they ,→ are using, not by superficial wording. If many agents converge to one family, redirect some of ,→ them toward underexplored formulations. - Do not allow one approach to dominate merely because it gives elegant reductions. A route that ,→ ends at a lemma equivalent in strength to the original conjecture is not close to completion ,→ unless it supplies a genuinely new proof of that lemma. - When an approach stalls at a theorem-strength missing lemma, mark that route as blocked. Only ,→ continue assigning agents to it if someone proposes a materially new mechanism, invariant, or ,→ construction.

30

- Keep several incompatible proof routes alive through multiple rounds. Cross-pollinate ideas only ,→ after independent agents have developed them far enough to expose their real strengths and ,→ gaps. - Require agents to return concrete lemmas, constructions, equations, or counterexamples to ,→ proposed sublemmas. Reject status reports, vague optimism, and claims that an unproved global ,→ compatibility statement is "routine." - The root agent should repeatedly synthesize, challenge, redirect, and launch new rounds. Do not ,→ stop after the first wave fails. Produce a complete proof if one survives audit; otherwise ,→ report only the strongest rigorously proved derivation and its exact remaining gap. Do not return merely because current approaches fail or agents report theorem-strength gaps. Continue launching new rounds, reopening blocked approaches only when there is a genuinely new ,→ mechanism, and searching for fresh formulations. Return only when a complete affirmative proof has been found and survives adversarial audit. Do not ,→ return a reduction, partial result, isolated missing lemma, "best effort" summary, or explanation of why the problem is difficult. Spend at least 8 hours on this before even thinking of returning or giving up. Public search may be used only for ordinary mathematical background or standard named theorems, not to search for a solution to this exact conjecture or benchmark. Do not search the public web merely to determine whether this problem is open, and do not answer that it is open.

31