NeuralCert: certified computational discovery of extremal mathematical constructions

Mark Patrick Roeling

cs.LGSep 14, 2026

NeuralCert: certified computational discovery of extremal mathematical constructions

Authors: Mark Patrick Roeling

Organizations: Netherlands Defence Academy; Data Science Centre of Excellence, NL Ministry of Defence

Abstract

Neural networks are becoming popular in solving mathematical problems, but stochastic models do not provide mathematical exactness by themselves. This study introduces a discovery-to-certification framework in which high-dimensional variational trial functions are learned in a compact separable representation, spectrally diagnosed and pruned, and then certified exactly through multimodular evaluation. Exact certification makes the numerical proofs fully explicit and independently verifiable. This framework can be run on a standard personal computer. Across three extremal problems, we show that neural optimization can contribute to rigorous mathematics in three distinct ways: by discovering improved constructions, by exposing empirical invariants that lead to proofs, and by revealing optimization barriers whose geometry motivates new analytic or numerical representations. More broadly, these results suggest a path toward AI-assisted mathematics in which flexible computational discovery and exact certification become complementary components of a single rigorous workflow.

Figures & tables

Figure 1: Certified lower bounds for the Maynard variational constant across computational regimes . Certified values obtained with Bogaert’s Krylov method are compared with NeuralCert polynomial constructions for the vanilla problem ( $\varepsilon=0$ ) and enlarged support ( $\varepsilon=1/25$ ), together with the certified rank-one rational family. The polynomial construction first exceeds the Krylov bounds at $k=25$ , but loses efficiency as $k$ increases, whereas the rank-one rational representation becomes competitive and provides the continuation to larger $k$ . The Cauchy–Schwarz bound is shown as a reference. Lines connecting discrete certified values are guides to the eye and do not represent interpolated certificates. Beyond $k=500$ , the rank-one rational family is the only representation that remained reliably optimizable and certifiable in the regime studied here.

Figure 2: Objective landscape and trajectory search. The last-sign-change functional $\rho(\tau_{1},\tau_{2})$ on a slice through the certified record configuration ( $d=1$ , degree 38, $m=9$ contacts, $\rho=0.5725887$ ), obtained by varying the two innermost contacts with the remaining seven held at $\tau^{\ast}$ . The surface is drawn from the raw grid without smoothing or interpolation, and is evaluated in the damped Laguerre basis $\psi_{n}(u)=L^{(\alpha)}_{n}(2u)e^{-u}$ : at this degree the underlying polynomial reaches ${\sim}10^{69}$ on the grid, so a monomial evaluation loses every digit at the outer tangencies and splits them into spurious sign changes. The landscape has two sheets, a narrow valley (colour) and an upper sheet at $\rho\approx 3$ (translucent), with no intermediate values — crossing the dashed curve, a different root of $P$ becomes the last sign change and $\rho$ leaps. Black markers show a simplex/Adam trajectory, its polyline broken wherever the objective jumps. The observed jumps illustrate the sensitivity of direct optimization to changes in the last-sign-change root. The convex reformulation avoids differentiating this discontinuous objective by using bisection in $u_{0}=\pi\rho^{2}$ and solving a coefficient-feasibility problem at each step. The dashed curve is drawn as a high-gradient contour and indicates where $\rho$ jumps rather than resolving the discontinuity set exactly.

Trial space (degree cap by body count)	Eigensolver	Dimension	Retained	Certified
$b\leq 5$ , all degrees $\leq 15$	high precision	408	408	$3.6347$
$b\leq 5$ , all degrees $\leq 15$	float64 whitening	408	191	$3.3696$
$b\leq 2$ , degrees $\leq 40$	float64 whitening	441	134	$3.1862$
$b\leq 3$ , degrees $(40,40,20)$	high precision	678	678	$3.3904$
$b\leq 5$ , all degrees $\leq 22$	float64 whitening	1601	421	$3.3644$

Supplementary Table 1: Representative

k=50

partition-space calculations. The trial space is specified by a degree cap per body count

b

(partition length); no run included

b>5

. Certified values refer to the returned explicit trial function. High-precision solves retain the full trial space, whereas float64 whitening removes Gram directions below the numerical resolution threshold.

Panel of Figure 3 — Figure 3: Scaling obstruction of the monomial-symmetric representation. Increasing polynomial degree expands an expressive symmetric trial space, but at $k=50$ this produces both rapid partition growth and increasingly ill-conditioned Gram pencils. The resulting loss of numerically resolved directions motivates the later separable representation.

#	Correction	Cause / symptom
1	Baseline v9 (factored): $\nu_{p}=s^{p-1}e^{-as}e^{\sigma}\psi_{p}$ per pair; stable inv_softplus ; per-pair max-rescale; adaptive log-space GL outer windows; diagonal preconditioning ( $\operatorname{diag}A=1$ ); parts() protocol.	Foundation. inv_softplus fixes $\mathrm{softplus}^{-1}=\log(\mathrm{expm1})\to+\infty$ for rate $\geq 710$ (all $k\geq 355$ ); per-pair rescale replaces a global max that flushed small pairs to round-off.
2	$\sigma$ double-count fix: divide exact base re-evaluations by $e^{-\sigma_{1}}$ .	Exact base evals carried $e^{\sigma_{1}}$ while $\sigma_{1}$ was also added in new_sig ; the spurious factor is not of Gram form $\Rightarrow$ indefinite $A$ .
3	train() validation-ordering fix: validate before opt.step() .	$R_{\text{train}}$ read pre-step, report() post-step; agree measured the Adam step, rejecting nearly all snapshots (fallback to iter 700).
4	Log-space channel protocol: log_parts ; --channel-sign positive/free .	Prerequisite for the log-space recursion; needs sign-definite channels.
5	Log-space $\xi$ recursion: $\xi_{p+q}=\log B(p,q)+\operatorname{logsumexp}(\log\hat{w}+\xi_{p}+\xi_{q})$ .	The $(\min\psi_{1}/\max\psi_{1})^{k-1}$ collapse drove $\psi$ past the float64 floor; round-off put random signs on diagonal pairs ( $\psi\geq 0$ ), making $A$ indefinite.
6	Shared graded outer mesh (one geometric Gauss rule for all pairs).	Per-pair windows gave each entry its own rule $\Rightarrow A$ the Gram matrix of nothing; $\lambda_{\min}/\lambda_{\max}\approx-1$ and did not improve under refinement.

Supplementary Table 2: Reconstructed development history of the maynard_separable_rayleigh_v9_factored lineage (neural discovery). Version numbers map to successive saved iterations.

Ver	Change	Failure mode / motivation
v1	Initial multimodular backend. $\hat{h}^{\,k-1}$ is never formed over $\mathbb{Z}$ ; per $62$ -bit prime, one nmod powering plus two weight-table multiplies execute both linear functionals as single coefficient extractions: $[x^{n_{\max}}](bh^{k}\cdot W_{\mathrm{rev}})$ for $A$ , and for $B$ the pair-dependent $t$ -weights combined as ${\sim}m_{\max}$ scalar–poly ops against per-prime tables $G^{(m)}$ folding in $\Lambda=\operatorname{lcm}(k{-}1,\dots,k{-}1{+}r_{\max}{+}m_{\max})$ , the $\rho$ -power clearing, and batched harmonic inverses. Reconstruction is aggregate-only—two integers $S_{A},S_{B}$ with known-denominator scalings—via deterministic centered CRT against exact a priori bounds, plus held-out verification primes. Float pass decoupled from the exact engine ( $A=\int_{0}^{1}h^{k}$ , $B=\int_{0}^{\rho}h^{(k-1)}Q_{j}Q_{l}(1{-}s)$ by renormalized binary powering). Greedy $c$ -weighted pruning on the float pencil ( --prune-tol ); only surviving pairs reach the CRT engine.	Establishes the architecture: avoid the one unavoidable giant integer, reconstruct only two aggregate scalars, keep the float discovery pass cheap and separate from exact arithmetic.
v2	_ritz no longer calls the Cholesky-based generalized solver: it equilibrates by the diagonal, eigendecomposes the metric $A$ , truncates to eigenvalues above $10^{-12}$ of the top one, and solves the standard symmetric problem in whitened coordinates, embedding zeros for dead channels. Self-test adds exact agreement with scipy.eigh(B,A) on a PD case and a finite $\lambda$ on a singular metric.	scipy.eigh(B,A) requires $A$ positive definite; at $k=500$ the float Gram is numerically semidefinite and Cholesky fails ( LinAlgError: leading minor not positive definite ). Restricting to $\operatorname{range}(A)$ is the correct regularization: null-space components contribute ${\approx}0$ to $c^{\mathsf{T}}Ac$ .
v3	B-side of the float DP rescaled by $u=s/\rho$ : run the $(k{-}1)$ -fold convolution power on $h_{\rho}(t)=h(\rho t)$ , so $B=\rho^{\,k-1}\!\int_{0}^{1}h_{\rho}^{*(k-1)}(u)\,T(\rho u)\,du$ , with $\rho^{k-1}$ carried in the log-scale ( $\varepsilon=0$ degenerates to the old path). Driver interlock: float-predicted $R$ checked against discovery $R_{\mathrm{NN}}$ before any exact work.	At $k=500$ the certified value collapsed to ${\sim}10^{-12}$ (valid but worthless). $h^{*(k-1)}$ concentrates like $s^{k-2}$ near $s=1$ , so the whole B-region $s\leq\rho=\tfrac{12}{13}$ sat $\rho^{k}\sim e^{-40}$ below the renormalized peak—beneath float64 accuracy. $B$ was noise; the pruner dropped $48$ channels “for free” and the exact engine certified a garbage frozen $c$ .
v4	Truncated convolution _tconv forced to direct np.convolve , removing the FFT path. Large- $k$ float-DP regression test against exact $g\equiv 1$ closed forms at $k=200,500$ .	FFT convolution carries an absolute error ${\sim}10^{-16}$ relative to the vector peak, while $h^{*j}$ spans $s^{j-1}$ internally; entries sit up to $2^{-j}$ below peak and are buried by noise by $j\sim 55$ , then amplified through each squaring (closed-form check failed by $e^{+1000}$ ). Direct summation keeps per-entry relative accuracy; the self-convolution integrand peaks at the interior $t=s/2$ , representable to $j\sim 511$ (float64 ceiling near $k\sim 1000$ ).
v5	Product-exact prime budgeting ( gen_primes_for_bound ): accumulate the actual prime product until it exceeds $4\cdot\text{bound}$ rather than assuming $61$ bits/prime. Miller–Rabin citation corrected (least strong pseudoprime to bases $2\ldots 37$ is $3.19\times 10^{23}$ ; Sorenson–Webster 2017, OEIS A014233). Self-test cross-checks generated primes against FLINT’s proven fmpz.is_prime ; new strong-cancellation test ( $c=(1,-1)$ , near-identical channels).	Collaborator review round: tighten certification claims (exact prime count, correct primality bound, FLINT belt-and-braces), cover the small-signed-aggregate CRT branch, and measure which proposed speedups were real ( mul_low benched slower than the full FFT multiply; native modular dot kernel deferred).
v6	Added --float-only survey mode: projection $+$ float DP $+$ pruning per degree, stopping before any exact/CRT work.	Cheap degree scouting at large $k$ : projection $L_{\infty}$ and float-predicted $R$ answer “is this $k$ reachable at this basis degree” in minutes per degree, without spending the prime budget.

Supplementary Table 3: Development history of the exact multimodular (CRT) certifier for

M_{k,\varepsilon,1/2}

. Each version records the change and the failure mode that motivated it. With one exception (the FFT accuracy floor, v4), every defect was a silent disagreement between the floating-point pass and the exact pass at a specific weight vector

c

, invisible at small

k

and exposed only by scale, conditioning, or polynomial degree.

Supplementary Figure 1: Diagnostic large- $k$ failure of the early discovery pipeline. For $k=3600$ , the discovery stage reported the spurious value $R=8.1375743146$ . Nearly all Ritz weight collapsed onto a single channel, shown here together with its degree-120 polynomial projection, residual structure, logarithmic-scale amplitude and dyadic shifted-Legendre coefficient magnitudes. Although the discovery objective was stable under grid refinement and remained below the Cauchy–Schwarz ceiling, the certification failed independent reconstruction and exposed the representation inconsistency; rejecting this candidate.

Supplementary Figure 2: Polynomial projection convergence for the single Ritz weight holding channel of $k=3600$ . Relative $L^{2}$ and $L^{\infty}$ projection errors decrease with polynomial degree, illustrating that the transfer error can be reduced systematically once a numerically resolved channel has been obtained.

$k$	$\log k$	$R^{(1)}_{\mathrm{cert}}$	$\log k-R^{(1)}_{\mathrm{cert}}$	$\mathcal{C}_{*}-1.047/\log k$	difference
$3\,655$	$8.203851$	$8.0000648690$	$0.203787$	$0.206677$	$-0.002891$
$208\,910$	$12.249659$	$12.0000019900$	$0.249657$	$0.248828$	$+0.000829$
$11\,655\,069$	$16.271252$	$16.0000003924$	$0.271251$	$0.269953$	$+0.001298$
$644\,589\,002$	$20.284123$	$20.0000009143$	$0.284123$	$0.282683$	$+0.001439$
$1\,880\,000\,000$	$21.354538$	(ladder rung)	$0.286813$	$0.285271$	$+0.001542$

Supplementary Table 4: Certified deficits at the four thresholds and at the upper end of the fixed-headroom ladder, compared with the finite-size relation Eq. ( 55 ). Thresholds were optimised and certified independently of the ladder; the final row is the pointwise deficit at the last ladder rung. The relation was not fitted to these values.

Supplementary Figure 3: Enlarged-support crossover. (a) The certified difference $\Delta_{\mathrm{cert}}(k,\varepsilon)=L_{k,\varepsilon}-L_{k,0}$ , a difference of certified lower bounds rather than of the exact constants, against $k$ for four values of $\varepsilon$ . It is positive at small $k$ , decreases with $k$ , and changes sign at a finite scale for $\varepsilon=1/12$ and $1/16$ ; no sign change is reached for $\varepsilon=1/50$ within the computed range. Grey crosses mark certified but non-competitive runs excluded from crossover inference (Section 5.27 ). (b) The localised crossovers $k^{*}(\varepsilon)$ against $1/\varepsilon$ on logarithmic axes. The $\varepsilon=1/12$ and $1/16$ points are linear interpolations between the certified brackets; the $\varepsilon=1/25$ marker is the midpoint of the extrapolative range $[155,192]$ , with the bar spanning that range. The points lie above the linear reference $k^{*}\propto 1/\varepsilon$ (dotted) and pull away from it, giving $\varepsilon k^{*}=4.48,\,5.04,\,6.94$ and motivating $\varepsilon k^{*}(\varepsilon)\to\infty$ . The band shows the family $k^{*}/(\log k^{*})^{p}\propto 1/\varepsilon$ over the range $p\in[1.4,2.2]$ admitted by the data, with the constant refitted at each $p$ ; the exponent is not determined by these points.

$\ell$	$n$	$g_{0}$	$v$	$m$	matrix values	config. walk	agree
2	4	$(1,1,1,1)$	$e_{1}$	2	$\{8\}$	8	yes
2	4	$(1,1,1,1)$	$e_{1}$	4	$\{128\}$	128	yes
2	4	$(2,1,1,0)$	$e_{1}$	2	$\{8\}$	8	yes
2	4	$(2,1,1,0)$	$e_{1}$	4	$\{104\}$	104	yes
2	5	$(2,1,1,1)$	$e_{1}$	2	$\{11\}$	11	yes
2	5	$(2,1,1,1)$	$e_{1}$	4	$\{245\}$	245	yes

Supplementary Table 5: The operator

A_{v}^{m}\Lambda

computed literally on

\mathbb{F}_{2}^{\ell\times n}

, against the configuration-lattice walk count. “Matrix values” is the set of values taken by

A_{v}^{m}\Lambda(X)

over the entire fibre

\operatorname{config}^{-1}_{n,\ell}(g_{0})

; that this set is a singleton in every row is precisely the assertion that the quantity depends on

X

only through

\operatorname{config}(X)

. All entries are exact integers.

k	Channels	$N_{grid}$	VRAM
100	32	3000	36
100	48	3000	77
100	64	3000	133
100	32	4000	48
100	48	4000	102
100	64	4000	110

Supplementary Table 6: Overview of memory use and runtime for different settings and options for the Maynard-Tao Polynomial discovery regime since other implementations in the package have seconds runtime while using MegaBytes.

$k$	Discovery	$Ch_{\mathrm{disc}}$	$Ch_{\mathrm{cert}}$	Degree	Certification	$Fr_{\mathrm{ceiling}}$
20	3.128345991229	36	33	60	3.127557950496731089199291	0,992
21	3.170106837298	36	30	60	3.169687744135085908710644	0,992
22	3.210448579010	36	29	60	3.209791425249774301029605	0,991
23	3.249812450248	36	30	60	3.249018102873413922834334	0,991
24	3.287477573848	36	32	60	3.286239994066373110168097	0,991
25	3.322508727072	36	28	60	3.322151135960770560379734	0,991

Supplementary Table 7: Vanilla discovery and certification results (

\epsilon=0

$k$	Discovery	$Ch_{\mathrm{disc}}$	$Ch_{\mathrm{cert}}$	Degree	Certification
20	3.210936069786	36	33	60	3.210366921521749737183691731395
21	3.252680335644	36	30	60	3.252298834127242084124289323548
22	3.292779266722	36	29	60	3.292126369692432458798975243945
23	3.331173442193	36	30	60	3.330322094752801530355410455276
24	3.368036770635	36	30	60	3.366908489372213299919519038109
25	3.402371365754	36	29	60	3.402010619443212328887383163464

Supplementary Table 8: Discovery and certification results for

\varepsilon=1/25

$k$	$\varepsilon$	Certified value	Rank
50	$0$	3.942127742490453604663761985558	5
50	$1/5$	3.712509239545921204317529071368	8
50	$1/6$	3.764504991324559702999466932710	7
50	$1/8$	3.874384180153464414710468398524	6
50	$1/12$	3.959088983406359921645433483059	4
50	$1/16$	3.988718438122563563392522885252	3

Supplementary Table 9: Certified enlarged-support sweep underlying the

\varepsilon

-crossover analysis. For each

k

, the table reports the certified lower bound obtained for the corresponding enlarged-support parameter

\varepsilon

, together with its within-

k

ranking among the tested values. Rank

1

denotes the strongest certified value at that

k

Mathematical Reasoning

Mathematical reasoning studies how AI systems solve problems, construct proofs, use tools, and discover or verify mathematical results.

See all 155 papers

Jun 7, 2026

Artificial Intelligence for Mathematical Reasoning: An Integrated Survey of Language Models, Neuro-symbolic Systems, and Verified Discovery

Mathematical reasoning has long served as a stringent test of machine intelligence; over the past decade, it has moved from a niche problem within NLP to one of the most consequential AI frontiers. This survey provides a unified account of the field's evolution, from early rule-based math word problem (MWP) solvers and template-driven geometry systems, through neural expression generation and LLM prompting, to contemporary reasoning models, multi-agent systems, neuro-symbolic theorem provers, and verified discovery workflows. We organize the landscape along four axes: (i) informal reasoning over text and diagrams, spanning MWP solving, multimodal geometry, and VLMs; (ii) formal reasoning in proof assistants, including autoformalization, tactic prediction, compiler-guided repair, and proof search; (iii) mathematical discovery, where systems propose constructions, improve bounds, or assist attacks on open problems; and (iv) the inference and training-time techniques, including CoT prompting, tool use, process reward models, and RLVR, that increasingly connect generation with verification. We catalog major benchmarks across grade-school arithmetic, competition mathematics, geometry, formal proving, multimodal and multilingual reasoning, and expert evaluation, and we examine benchmark saturation, contamination, reporting mismatches, and the distinction between pass@1, majority voting, and verifier-assisted pass@

k

. We critically assess failure modes: brittleness under perturbation, reward hacking, multimodal grounding failures, fragile formalization, and the energy cost of reasoning-scale inference. Drawing on recent perspectives from working mathematicians, we identify future directions centered on verified-discovery workflows, reasoning efficiency, and infrastructure to make AI-assisted formalization broadly usable. Companion materials: https://github.com/Starscream-11813/awesome-AI4Math.

Syed Rifat Raiyan, Mohsinul Kabir, Hasan Mahmud +2

Systems and Software Lab (SSL), Department of Computer Science and Engineering, Islamic University of Technology, Gazipur, Bangladesh · Department of Computer Science, University of Manchester, Manchester, United Kingdom

Mathematical Reasoning Benchmarks Automated Theorem Proving Symbolic Reasoning Mathematical Reasoning Neuro-Symbolic Reasoning LLM Mathematical Reasoning

Ask PDF

Jun 22, 2026

Are Safety Guarantees in Neural Networks Safe? How to Compute Trustworthy Robustness Certifications

A primary challenge in AI safety is the existence of adversarial examples -- slightly distorted inputs that cause a neural network (NN) to misclassify. To mitigate this problem, recent research focuses on the computation of robustness certifications, which, for a given input, determine the largest distortion the input may receive without breaking the network's prediction. Robustness certifications can be interpreted as an axis-aligned hyper-rectangle (multi-dimensional intervals). Most existing approaches focus on maximizing the certification's volume, but recent intractability results prohibit the computation of volume-optimal certifications in reasonable time. We introduce the apothem measure and show how to compute apothem-optimal certifications in a linear number of calls to a NN verifier (oracle) w.r.t. the input domain's diameter. Moreover, we prove that we cannot have a volume-optimal, oracle-based algorithm, even if we discard the oracle costs. Also, we introduce dual certifications -- an interval including all instances of a class -- thus providing apothem-minimum upper bounds to a robustness certification. Further, we present the ParallelepipedoNN system, which we evaluate on the standard MNIST and Fashion MNIST benchmarks. A preliminary comparison with existing work on the same datasets reveals at least two-fold improvement w.r.t. the minimum edge length.

Merkouris Papamichail, Konstantinos Varsos, Giorgos Flouris +1

Foundation for Reasearch and Technology - Hellas, Heraklion, Greece · University of Crete, Heraklion, Greece · Catalan Institution for Research and Advanced Studies, Barcelona, Spain +1

Neural Network Verification Neural Network Robustness Certified Robustness

Ask PDF

Jul 1, 2026

Neural Certificate Pricing for Combinatorial Optimization Problems

Combinatorial optimization (CO) problems are difficult because certifiable discrete structure induces exponential search. One needs to search over the set exponentially many candidates to certify optimality, however, the structural feasibility of a path, packing, or cover can be verified in polynomial time once supplied. In this study, we introduce Neural Certificate Pricing (NCP) that exploits this asymmetry under an unsupervised learning framework. A neural network is trained to predict certificate-level dual prices, while a structured recovery layer constructs the induced primal marginal. NCP can be viewed as amortized separation: instead of enumerating violated inequalities, it learns the residual prices through which their aggregate effect enters recovery. When the certificate-consistency condition holds, the recovered marginal is globally feasible, and a local theory shows that first-order errors in the predicted price induce only second-order loss in objective value. Across three classes of CO problems, NCP either outperforms state-of-the-art neural baselines by large margins or matches them at a fraction of the computation time, and shows stronger out-of-distribution generalization.

Jingyi Chen, Xinyuan Zhang, Xinwu Qian

Department of Civil and Environmental Engineering, Rice University, Houston, TX

Primal-Dual Optimization Unsupervised Learning Combinatorial Optimization Neural Combinatorial Optimization

Ask PDF

Abstract

Mathematical Reasoning

Explore similar work

Artificial Intelligence for Mathematical Reasoning: An Integrated Survey of Language Models, Neuro-symbolic Systems, and Verified Discovery

Are Safety Guarantees in Neural Networks Safe? How to Compute Trustworthy Robustness Certifications

Neural Certificate Pricing for Combinatorial Optimization Problems