Paper deep dive
A user's guide to PINNs in geometric analysis: lessons from the asymptotic Plateau problem
Tancredi Schettini Gherardini
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 91%
Last extracted: 8/3/2026, 2:08:26 AM
Summary
This paper presents a methodological guide for applying Physics-Informed Neural Networks (PINNs) to geometric analysis, specifically the asymptotic Plateau problem in hyperbolic space. It details two key innovations: encoding geometry directly into the neural network architecture to satisfy boundary conditions exactly, and optimizing the evaluation of PDE residuals using forward-mode automatic differentiation and computational graph compilation. These techniques significantly reduce training time while providing numerical evidence for Joel Fine's conjecture relating minimal surfaces in H^4 to HOMFLY polynomial coefficients.
Entities (9)
Relation Signals (6)
PINNs → appliedto → Asymptotic Plateau Problem
confidence 95% · introduced a machine learning framework based on physics-informed neural networks (PINNs), aimed at constructing near-minimal discs in hyperbolic space asymptotic to a prescribed knot
Asymptotic Plateau Problem → definedin → Hyperbolic Space H^4
confidence 95% · The asymptotic Plateau problem asks whether a closed submanifold... of the sphere at infinity of hyperbolic space Hn+1 bounds a complete, properly immersed minimal submanifold
PINNs → providesevidencefor → Fine's Conjecture
confidence 92% · We used this method to provide numerical evidence for a conjecture of Joel Fine relating minimal surfaces in H4 to the coefficients of the HOMFLY polynomial
Forward-Mode Automatic Differentiation → optimizes → PINNs
confidence 90% · replacing nested reverse-mode automatic differentiation with the forward propagation of second-order jets... reduce the cost of a training step by a factor of roughly forty to fifty
Computational Graph Compilation → optimizes → PINNs
confidence 90% · compiling the computational graph of the residual once instead of rebuilding it at every optimisation step... reduce the cost of a training step by a factor of roughly forty to fifty
PINNs → usedby → Marco Usula
confidence 90% · a joint work with Marco Usula, where we introduced a machine learning framework based on physics-informed neural networks (PINNs)
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:This proceedings contribution elaborates on the findings of arXiv:2605.26234v2: a joint work with Marco Usula, where we introduced a machine learning framework based on physics-informed neural networks (PINNs), aimed at constructing near-minimal discs in hyperbolic space asymptotic to a prescribed knot at infinity. We used this method to provide numerical evidence for a conjecture of Joel Fine relating minimal surfaces in $H^{4}$ to the coefficients of the HOMFLY polynomial. This is a methodological companion to that paper, based on a presentation given at the 2026 edition of the workshop "DANGER: Data, Numbers, and Geometry". Rather than reviewing the results, which are presented extensively in the preprint above, we discuss the two aspects of the framework which, in our experience, determined whether the method worked at all. First, the geometry of the problem must be encoded in the architecture of the model, so that the boundary condition and asymptotics at infinity hold exactly for every value of the learnable parameters - leaving us with a single-component loss function; second, the evaluation of the PDE residual must be engineered with care to ensure that complete trainings can be performed in a reasonable time. On the latter point, we describe two implementation techniques which are not spelled out in detail in the original paper: replacing nested reverse-mode automatic differentiation with the forward propagation of second-order jets, and compiling the computational graph of the residual once instead of rebuilding it at every optimisation step. Together, on identical hardware, these two changes reduce the cost of a training step by a factor of roughly forty to fifty. We hope these methodological discussions can be useful for researchers in differential geometry and geometric analysis who wish to deploy PINNs on problems of their own.
Tags
Links
- Source: https://arxiv.org/abs/2607.28733v1
- Canonical: https://arxiv.org/abs/2607.28733v1
Trouble viewing inline? Open PDF directly →
Full Text
44,245 characters extracted from source content.
Expand or collapse full text
A user’s guide to PINNs in geometric analysis: lessons from the asymptotic Plateau problem Tancredi Schettini Gherardini Abstract. This proceedings contribution elaborates on the findings of arXiv:2605.26234v2: a joint work with Marco Usula, where we introduced a machine learning framework based on physics-informed neural networks (PINNs), aimed at constructing near-minimal discs in hyperbolic space asymptotic to a prescribed knot at infinity. We used this method to provide numerical evidence for a conjecture of Joel Fine relating minimal surfaces in H4 H^4 to the coefficients of the HOMFLY polynomial. This is a methodological companion to that paper, based on a presentation given at the 2026 edition of the workshop “DANGER: Data, Numbers, and Geometry”. Rather than reviewing the results, which are presented extensively in the preprint above, we discuss the two aspects of the framework which, in our experience, determined whether the method worked at all. First, the geometry of the problem must be encoded in the architecture of the model, so that the boundary condition and asymptotics at infinity hold exactly for every value of the learnable parameters – leaving us with a single-component loss function; second, the evaluation of the PDE residual must be engineered with care to ensure that complete trainings can be performed in a reasonable time. On the latter point, we describe two implementation techniques which are not spelled out in detail in the original paper: replacing nested reverse-mode automatic differentiation with the forward propagation of second-order jets, and compiling the computational graph of the residual once instead of rebuilding it at every optimisation step. Together, on identical hardware, these two changes reduce the cost of a training step by a factor of roughly forty to fifty. We hope these methodological discussions can be useful for researchers in differential geometry and geometric analysis who wish to deploy PINNs on problems of their own. University of Bonn and Max Planck Institute for Mathematics, Bonn, Germany. E-mail: tsg@math.uni-bonn.de. 1. Introduction The asymptotic Plateau problem asks whether a closed submanifold Mk⊂SnM^k⊂ S^n of the sphere at infinity of hyperbolic space Hn+1 H^n+1 bounds a complete, properly immersed minimal submanifold of Hn+1 H^n+1. Anderson’s classical theorem [1] guarantees the existence of an area-minimising solution, but says nothing about the structure of the full space of solutions. For surfaces in H4 H^4 bounding knots K⊂S3K⊂ S^3, Fine [5] proved that (for generic K) the minimal surfaces of fixed genus g and self-intersection number d bounding K form a compact oriented 0-dimensional moduli space when g=0g=0, whose signed count is a knot invariant, and explained why these counts should be closely related to the coefficients of the HOMFLY polynomial [6] of K; this precise conjectural correspondence was numerically tested and corroborated in [18]. The paper presented a numerical framework, based on physics-informed neural networks (PINNs) [17], which produces near-minimal discs in H4 H^4 bounding a prescribed knot, locates their self-intersections, and computes the sign of each double point. Testing the resulting self-intersection numbers against the HOMFLY coefficients – for the unknot, the torus knots T(3,2)T(3,2), T(5,2)T(5,2), T(4,3)T(4,3) and T(5,3)T(5,3), the figure-eight, three-twist, Stevedore and square knots, and various mirror images – we found agreement with Fine’s conjecture in every case. We refer to [18] for the geometric background, the precise statements, and the results; the trained models and the full implementation are publicly available111https://github.com/Tancredi-Schettini-Gherardini/deep_plateau. The present note, prepared for the proceedings of the “DANGER: Data, Numbers, and Geometry” 2026 workshop, where this work was presented, deliberately takes a different angle. It is written for geometers and analysts who are considering PINNs as a tool, and its subject is not what we computed but how. Neural networks are being used across geometric analysis and adjacent fields [4, 16, 12, 9, 11, 22, 10], including as a component of computer-assisted proofs [7, 21], and in our experience the difference between a PINN that converges to 10−610^-6 residuals in under an hour on a laptop and one that stalls, or takes days, is not primarily the network or the optimiser: it lies in choices that are usually left undocumented. Two of these choices form the theses of this note. (1) The geometry of the problem should be encoded in the model, not learned. Our model satisfies the boundary condition at infinity exactly, for every value of the learnable parameters, and is asymptotically minimal by construction: the analytically known leading behaviour of solutions near the boundary is built into the architecture, and the network only learns the interior. This removes the boundary term from the loss entirely, and with it the notoriously delicate balancing of boundary against interior penalties (§3); it yields a single-component loss encoding only one geometric property, which in our case is minimality. (2) The evaluation of the PDE residual should be engineered, not delegated. A PINN loss for a second-order equation needs first and second derivatives of the model with respect to its inputs at every collocation point, inside an outer optimisation loop that differentiates with respect to the parameters. The default implementation – nested reverse-mode automatic differentiation, rebuilt at every step – is wasteful in a way that compounds multiplicatively. Choosing the direction of differentiation by counting dimensions (forward for inputs, reverse for parameters), propagating second-order jets through the network in closed form, and compiling the resulting computational graph once, sped up our training by a factor of 4040–5050 on identical hardware (§4); neither technique appears in the companion paper. Section 2 recalls the minimal amount of geometric setup needed to make this note self-contained and explains where the standard PINN recipe struggles on this problem. Sections 3 and 4 develop the two theses. Section 5 distils the discussion into a short list of guidelines which, we believe, transfer to many geometric variational problems beyond the asymptotic Plateau problem. Acknowledgements The results surveyed here were obtained jointly with Marco Usula, whom the author thanks for the collaboration and for many conversations reflected in this note. We also thank Ed Hirst for the invitation to the DANGER workshop, and all the organisers for putting together such a great event. Finally, we also thank once more Joel Fine for his feedback and encouragement during the preparation of the main preprint. We acknowledge the contribution of general-purpose AI agents (Claude Opus 4.8 and Fable 5) in the implementation of the code and in the preparation of the first draft of this paper. We acknowledge the support of the 2024 Max Planck-Humboldt Research Award, bestowed on Geordie Williamson by the Max Planck Society and the Alexander von Humboldt Foundation and hosted by Catharina Stroppel at the University of Bonn. 2. The problem as a PINN task 2.1. Minimal discs in half-space coordinates We work in the upper half-space model of hyperbolic space: on ℝ>0×ℝnR_>0×R^n with coordinates (X,)=(X,Y1,…,Yn)(X, Y)=(X,Y_1,…,Y_n), the hyperbolic metric is ghyp=X−2(dX2+|d|2)g_ hyp=X^-2 (dX^2+|d Y|^2 ), and the boundary at infinity of the closure is the locus X=0≅ℝn\X=0\ ^n, a stereographic chart of SnS^n minus a point. The object we seek is a map u=(X,):D2⟶H¯n+1u=(X, Y):D^2 H^n+1 from the closed unit disc, mapping the interior to the interior, whose boundary restriction u|∂D2u_|∂ D^2 parametrises a prescribed embedded curve γ:S1→ℝnγ:S^1 ^n (for n=3n=3: a knot) in the boundary at infinity, and whose interior image is minimal. In the language of [18, 19, 20], u should be a minimal p-immersion. Minimality is expressed by the vanishing of the tension field τ(u)=tru∗ghyp∇duτ(u)= tr_u^*g_ hyp∇ du, computed with respect to the pull-back metric; in half-space coordinates, writing g=X−2J⊤Jg=X^-2J J for the pull-back metric (J the Jacobian of u in the disc coordinates), the components of τ(u)τ(u) in the orthonormal frame (X∂X,X∂Yk)(X _X,X _Y_k) read (2.1) τX(u) τ^X(u) =1X[ΔgX+1X(∑k=1n|dYk|g2−|dX|g2)], = 1X [ _gX+ 1X ( _k=1^n|dY_k|_g^2-|dX|_g^2 ) ], (2.2) τYk(u) τ^Y_k(u) =1X[ΔgYk−2X⟨dX,dYk⟩g],k=1,…,n, = 1X [ _gY_k- 2X\, dX,dY_k _g ], k=1,…,n, where Δg _g is the Laplace–Beltrami operator of g, acting on a scalar f by (2.3) Δgf=gab∂ab2f+(∂agab)∂bf+gab(∂alogdetg)∂bf. _gf=g^ab _ab^2f+( _ag^ab) _bf+g^ab( _a g) _bf. The equation τ(u)=0τ(u)=0 is a quasi-linear second-order elliptic system, degenerate at the boundary, with no general method available to solve it explicitly. Training minimises the Monte Carlo L2L^2 norm of the residual, (2.4) ℒ(θ)=1N∑i=1N|τ(uθ)|2(pi),L(θ)= 1N _i=1^N |τ(u_θ) |^2(p_i), over a sample pi⊂D2\p_i\⊂ D^2 of interior collocation points; here θ denotes the learnable parameters of the model uθu_θ, and the pointwise norm is taken with respect to the hyperbolic metric, i.e. |τ|2=(τX)2+∑k(τYk)2|τ|^2=(τ^X)^2+ _k(τ^Y_k)^2 in the frame above. 2.2. The standard recipe, and where it strains The textbook PINN prescription for a boundary value problem ℱ[u]=0F[u]=0 in Ω , ℬ[u]=0B[u]=0 on ∂Ω∂ , is to take uθu_θ to be a multi-layer perceptron (MLP) and to minimise a weighted sum of an interior and a boundary penalty, ℒ(θ)=1NΩ∑i|ℱ[uθ](xi)|2+λN∂Ω∑j|ℬ[uθ](xj′)|2,L(θ)= 1N_ _i |F[u_θ](x_i) |^2+ λN_∂ _j |B[u_θ](x_j ) |^2, the derivatives in ℱF being supplied by automatic differentiation. This recipe is general, and its very generality is its weakness on a problem like ours, for three reasons. First, the weighting problem: λ trades boundary accuracy against interior accuracy, the optimal trade-off drifts during training, and a wrong choice produces either a surface that misses its boundary curve or one that fits the boundary while violating the equation; a sizeable literature on adaptive weighting schemes exists precisely because this balance is fragile. Second, the boundary is at infinity: the ambient metric blows up as X→0X→ 0, so uniform penalties on the map do not translate into uniform geometric control, and a network that is free to move the boundary even slightly produces maps that are not asymptotic to the prescribed knot at all. Third, the regularity theory is known and is not exploited: by a theorem of Marx-Kuo [15], a minimal p-submanifold of H¯n+1 H^n+1 with smooth boundary at infinity meets the boundary orthogonally, and admits a precise polyhomogeneous expansion there. A generic soft-constrained network knows nothing of this and must spend its capacity discovering it. The framework of [18] addresses these three points, as we now explain. 3. Encoding the geometry in the model 3.1. A hard-constrained ansatz Instead of asking an MLP to represent u directly, the model wraps the network inside a composite map whose structure carries the geometry. Fix an embedding γ:S1→ℝ3γ:S^1 ^3 (the knot), a boundary defining function ρ for the disc (a smooth function on D2D^2, positive in the interior, vanishing on ∂D2∂ D^2, but whose exterior derivative is non-zero on ∂D2∂ D^2), an extension operator ext ext producing a map ext(γ):D2→ℝ3 ext(γ):D^2 ^3 with ext(γ)|∂D2=γ ext(γ)_|∂ D^2=γ, and a decay exponent k∈1,2k∈\1,2\. The model is (3.1) uθ=(ρeNNX,ext(γ)+ρkN),u_θ= (\;ρ\,e^N^X,\; ext(γ)+ρ^k\,N Y ), where N=(NNX,N):ℝ2→ℝ4N=(N^X,N Y):R^2 ^4 is the only learnable ingredient — in all our experiments a plain MLP with four hidden layers of width 6464 and tanh activations, 12 93212\,932 parameters in total. Three structural facts hold for every value of θ: (1) the first component is positive in the interior and vanishes exactly on ∂D2∂ D^2, so the image lies in H4 H^4 and reaches the boundary at infinity precisely over ∂D2∂ D^2; (2) the boundary restriction of uθu_θ is exactly γ: the factor ρkρ^k kills the network correction at the boundary, where ext(γ) ext(γ) restricts to γ; (3) consequently the loss (2.4) contains no boundary term at all: the weighting problem of §2.2 disappears, because the constraint manifold has been parametrised away rather than penalised. A fourth fact concerns initialisation, and we single it out because it guided every later refinement: at θ=0θ=0 the model reduces to the explicit map (ρ,ext(γ))(ρ, ext(γ)). The network is therefore best understood as a learnable perturbation of a built-in starting point, and choosing ρ and ext ext amounts to choosing the best explicit approximation to a solution that one can write down. For the round unknot γ(φ)=(cosφ,sinφ,0)γ( )=( , ,0) (we write (r,φ)(r, ) for polar coordinates on D2D^2, reserving θ for the learnable parameters) the best choice is obvious: the totally geodesic copy of H2 H^2, (3.2) uH2(r,φ)=(1−r21+r2,2rcosφ1+r2,2rsinφ1+r2, 0),u_ H^2(r, )= ( 1-r^21+r^2,\, 2r 1+r^2,\, 2r 1+r^2,\,0 ), is an exact solution. Reading off its first component motivates the stereographic boundary defining function ρst=(1−r2)/(1+r2) _st=(1-r^2)/(1+r^2), which we use throughout: with it (and the extension below), θ=0θ=0 reproduces (3.2) exactly when γ is the round unknot. 3.2. Choosing the extension: a regularity and asymptotics story The extension operator is the heart of the construction, and finding the right one took several iterations over the life of the project. We record the sequence, because each step encodes a transferable lesson. Radial fattening The earliest implementations extended the knot as a cone, ext(γ)(r,φ)=c(r)γ(φ) ext(γ)(r, )=c(r)\,γ( ) with c(r)=rc(r)=r or 2r1+r2 2r1+r^2, and smoothed variants. The latter restricts to γ at r=1r=1 and, paired with ρst _st, is again exact for the round unknot; but for a generic curve it is merely continuous at the origin, where the cone has a genuine singularity. Since the network correction in (3.1) is smooth, no value of θ can repair a non-smooth baked-in ingredient: in experiments, training stalls with the error concentrated near the origin, and ad hoc smoothings of the cone tip remove the singularity only at the cost of the boundary asymptotics. The simple lesson: the regularity of every hard-coded ingredient is a ceiling on the regularity of the model – hard constraints transfer their defects to the solution just as reliably as their virtues. Harmonic extension Every smooth γ:S1→ℝ3γ:S^1 ^3 has a unique harmonic extension Γ to the disc, which is smooth up to the boundary and computable to machine precision from the Fourier coefficients of γ: writing z=x+iyz=x+iy and γ∼A0+∑m≥1(Amcosmφ+Bmsinmφ)γ A_0+ _m≥ 1(A_m m +B_m m ), one has Γ=A0+∑m≥1(UmAm+VmBm) =A_0+ _m≥ 1 (U_mA_m+V_mB_m ) with Um+iVm=zmU_m+iV_m=z^m. The stereoharmonic extension ext(γ)=2Γ/(1+r2) ext(γ)=2 /(1+r^2) again reproduces (3.2) for the round unknot, and is now smooth for every γ. This already trains well. Biharmonic extension and orthogonality at infinity The remaining defect concerns asymptotics rather than smoothness. As recalled in §2.2, an actual minimal p-submanifold meets the boundary at infinity orthogonally; equivalently, it is asymptotically minimal, i.e. |τ|2→0|τ|^2→ 0 at the boundary. For a map of the form (3.1) one computes that orthogonality holds if and only if the normal derivative ∂r(ext(γ)+ρkN) _r ( ext(γ)+ρ^kN Y ) is tangent to the curve along ∂D2∂ D^2 [18, §3.2]. If k=1k=1 the network can, in principle, learn to correct a non-orthogonal extension; if k=2k=2 its correction decays too fast to do so, and the extension alone must carry the right asymptotics. This suggests imposing a Neumann condition on the extension in addition to the Dirichlet one — which is exactly what the bi-Laplacian allows. The stereobiharmonic extension is ext(γ)=2Γ/(1+r2) ext(γ)=2 /(1+r^2), where Γ is the unique solution of Δ2Γ=0,Γ|S1=γ,∂rΓ|S1=γ; ^2 =0, _|S^1=γ, _r _|S^1=γ; mode by mode, Γ=A01+r22+∑m≥1Wm(r2)(UmAm+VmBm),Wm(r2)=m+12+1−m2r2, =A_0\, 1+r^22+ _m≥ 1W_m(r^2) (U_mA_m+V_mB_m ), W_m(r^2)= m+12+ 1-m2\,r^2, so it is exactly as computable as the harmonic one. Paired with k=2k=2, the resulting model is orthogonal to the boundary — hence asymptotically minimal — for every value of the learnable parameters, while still reducing to (3.2) at θ=0θ=0 for the round unknot. This is the production configuration used for all results in [18]. Figure 3.1 shows the practical meaning of these choices: the pointwise residual |τ(uθ)|2|τ(u_θ)|^2 of an untrained model already vanishes along the boundary; training only needs to remove the error away from the degenerate boundary region, on compact subsets of the interior where the equation is uniformly elliptic. In our experience this design is the single largest contributor to the reliability of the method: with it, plain Adam [13] followed by L-BFGS [14] reaches expected residuals typically in the 10−710^-7–10−410^-4 range across the knots we tested (up to 5×10−45× 10^-4 for the ten-crossing 1012410_124), with essentially no problem-specific tuning. Figure 3.1. Pointwise squared residual |τ(uθ)|2|τ(u_θ)|^2 for a (perturbed) unknot boundary, before training (left) and after training (right). Note the colour scales. The residual vanishes along ∂D2∂ D^2 even before training: the model is asymptotically minimal by construction, and the optimiser only works in the interior. Figure reproduced from [18]. Remark 3.1. Two comments are needed, for the sake of completeness. First, the boundary condition is satisfied exactly with respect to the degree-1515 Fourier truncation of γ used by the extension; for the explicit trigonometric parametrisations we use, the truncation error is negligible, but the statement “exact boundary condition” should be understood in this sense. Second, hard constraints are not free: they require the constraint manifold (the space of maps with the prescribed boundary and asymptotics, in this case) to be explicitly parametrised, which is possible here because the domain is a disc and the boundary data is a curve. When no such parametrisation is available, penalty methods remain the fallback. Remark 3.2 (What did not work). Since negative results rarely get reported, we list explorations that were tried and abandoned. Periodic-activation (SIREN-type) correctors diverged: their large second derivatives, amplified by the X−2X^-2 factors in (2.1)–(2.2), destabilised training. Feature embeddings in solid harmonics create a gradient blind spot at the origin, where all features vanish to high order. Residual-driven adaptive resampling of collocation points and self-adaptive pointwise loss weights brought no measurable benefit – plausibly because the hard-constrained model leaves no boundary layer for adaptivity to chase – and both interact badly with the static-shape compilation strategy of §4.3. A learnable reparametrisation of the disc, i.e. geometric preconditioning, helped in some specific cases but was not needed in the final pipeline. 4. Evaluating the residual efficiently The second half of this note concerns a purely computational question: what does it cost to evaluate (2.4) and its parameter gradient, and how should that computation be organised? The companion paper records only that second-order derivatives are computed by composing forward- and reverse-mode automatic differentiation [18, §3.3], without entering the details about the cost of the residual evaluation or how to organise it; yet the difference between the naive and the final organisation is a factor of 4040–5050 in wall-clock time (Table 1), on the same hardware, for results identical up to round-off. We describe the two ingredients in turn. 4.1. What a training step must compute Fix a mini-batch of N collocation points. Evaluating (2.1)–(2.3) at each point requires the value u∈ℝ4u ^4, the Jacobian J∈ℝ4×2J ^4× 2 and the Hessian H∈ℝ4×2×2H ^4× 2× 2 of the model with respect to the disc coordinates; everything else – the pull-back metric g=X−2J⊤Jg=X^-2J J, its explicit 2×22× 2 inverse, the derivative terms in (2.3), the residual itself – is plain tensor algebra in (u,J,H)(u,J,H): the full residual is an algebraic expression in the 22-jet of u, with no third derivatives appearing. The optimiser then needs the gradient of the scalar loss with respect to the ∼1.3×104 1.3× 10^4 network parameters. Automatic differentiation (AD) [3, 8] offers two elementary modes, which a geometer will recognise immediately: forward mode computes pushforwards, mapping an input tangent vector to an output tangent vector at the cost of one augmented evaluation; reverse mode computes pullbacks of output covectors. The cost of assembling a full Jacobian is therefore one sweep per input dimension in forward mode, and one sweep per output dimension in reverse. For the parameter gradient the choice is forced and classical: the loss is a single scalar depending on ∼104 10^4 parameters, so reverse mode, i.e. backpropagation, is optimal, and every deep-learning framework provides it. For the spatial derivatives the dimension count points the other way: the map goes from ℝ2R^2 to ℝ4R^4, so forward mode needs 22 sweeps against reverse mode’s 44; and for the Hessian the disparity compounds. Yet the naive PINN implementation ignores the count: one calls the framework’s reverse-mode routine once per output component (with the option that keeps the derivative computation itself differentiable), then differentiates the resulting expressions again for second derivatives. Our first draft of the code did exactly this, in the divergence form of (2.3): per residual evaluation it issued sixteen reverse sweeps – four for J, then twelve more inside the two divergence-form Laplacians ΔgX _gX and Δg _g Y – each sweep materialising a new computational graph on top of the previous one, and the whole tower being differentiated once more, in reverse, for the parameter gradient. This is the pattern that Table 1, row 1, prices at 427427 ms per training step. 4.2. Forward propagation of second-order jets The final implementation computes (u,J,H)(u,J,H) in a single forward pass, with no spatial AD at all, by propagating second-order jets through the network – forward-mode differentiation carried out explicitly, order two at a time (in AD terminology, Taylor-mode [8]). An MLP is a composition of affine maps and coordinatewise activations, and both have trivial jet-transport rules. Writing (a,Ja,Ha)(a,J_a,H_a) for the value, Jacobian and Hessian of the signal at the current layer (batched over collocation points): (4.1) affine a↦Wa+b: a Wa+b: (a,Ja,Ha)⟼(Wa+b,WJa,WHa), (a,J_a,H_a) (Wa+b,\;WJ_a,\;WH_a), (4.2) activation a↦σ(a): a σ(a): (a,Ja,Ha)⟼(σ(a),σ′(a)⊙Ja,σ′(a)⊙Ja⊗Ja+σ′(a)⊙Ha), (a,J_a,H_a) (σ(a),\;σ (a) J_a,\;σ (a) J_a J_a+σ (a) H_a ), with ⊙ the coordinatewise product and (Ja⊗Ja)ij=∂ia∂ja(J_a J_a)_ij= _ia\, _ja per neuron; for σ=tanhσ= one has σ′=1−σ2σ =1-σ^2 and σ′=−2σσ′σ =-2σ , computable from the already-evaluated activation. Initialising with Ja=Id2J_a=Id_2, Ha=0H_a=0 at the input and applying (4.1)–(4.2) layer by layer yields the exact 22-jet of the network at every collocation point in one pass whose cost is a small constant multiple of a plain forward evaluation – roughly the number of jet components carried per neuron, six or seven depending on whether the symmetry of the Hessian is exploited. The composite model (3.1) multiplies the network by explicit functions, and jets compose by the Leibniz rule. The 22-jets of the two geometric ingredients are available in closed form: for ρst _st by elementary calculus, and for the stereobiharmonic extension because differentiating the harmonic polynomials is algebraic (∂xUm=mUm−1 _xU_m=mU_m-1, ∂x2Um=m(m−1)Um−2 _x^2U_m=m(m-1)U_m-2, and so on), so the jet of the extension costs one more pass over the same Fourier recursion. The 22-jet of (3.1) is then assembled by the product rule – for instance the first component X=ρeνX=ρ e^ν, ν=NNXν=N^X, has dX=eν(dρ+ρdν),HessX=eν(Hessρ+dρ⊗dν+dν⊗dρ+ρdν⊗dν+ρHessν),dX=e^ν(dρ+ρ\,dν), \,X=e^ν (Hess\,ρ+dρ dν+dν dρ+ρ\,dν dν+ρ\,Hess\,ν ), and similarly for Yj=ext(γ)j+ρkwjY_j= ext(γ)_j+ρ^kw_j (with k the fixed decay exponent of (3.1)) – after which the residual is evaluated algebraically as in §4.1. Three properties of this organisation matter in practice. (1) It is exact: all derivatives follow exact rules, and the result agrees with the nested-autograd implementation to 5×10−145× 10^-14 in double precision (§4.4), i.e. to accumulated round-off. (2) It separates the two differentiations: spatial derivatives no longer pass through the AD engine at all – collocation points are plain data, and no graph of graphs is ever built; the single remaining reverse pass, for the parameter gradient, differentiates an ordinary composition of tensor operations. (3) It mirrors the mathematics: the implementation is literally the chain rule, the Leibniz rule and the coordinate formula (2.3), in the order in which one would write them on paper. This makes it auditable line by line – a non-trivial virtue when the output feeds a conjecture test – and it makes the whole residual a flat, branch-free tensor program, which is precisely what the compiler discussed next requires. On its own, this reorganisation accounts for a factor ≈22≈ 22–2424 (Table 1, row 3). A remark on precision, complementing [18, §3.3]: we run everything in double precision, since in single precision the rounding incurred in second-derivative computations makes the residual saturate near 10−410^-4 even at a true solution, defeating both the L-BFGS refinement and the double-point analysis. This choice interacts with hardware: consumer GPUs execute float64 at a small fraction of their nominal throughput (and Apple-silicon GPUs not at all), so the pipeline was designed CPU-first – which raises the stakes for the compilation strategy below, as the CPU has no raw-power headroom to hide overheads in. 4.3. Compiling the computational graph once A PyTorch program normally runs in eager mode: every tensor operation is dispatched to a pre-compiled kernel as the Python interpreter encounters it, and the computational graph exists only implicitly, rebuilt at every step for the benefit of the backward pass. For a residual like ours, i.e. hundreds of small operations on modest batches, the per-operation overhead (Python dispatch, kernel launch, allocation of intermediates) dominates the arithmetic. The remedy, available in modern frameworks under names such as torch.compile [2] or jax.jit, is to trace the computation once into an explicit graph, optimise it (fusing chains of elementwise operations into single kernels, eliminating dead code and redundant intermediates), generate machine code – C++ kernels, in our CPU case – and thereafter re-run the compiled artefact, forward and backward, at every step. The catch is that tracing specialises: the compiled artefact is valid for a fixed operation sequence and (in the static regime we use) fixed tensor shapes, and anything that varies from step to step triggers either a fallback to eager execution or an expensive silent recompilation. Concretely, the final implementation: (1) packages the entire residual evaluation — jet propagation (4.1)–(4.2), closed-form jets of ρ and ext(γ) ext(γ), Leibniz assembly, metric algebra, residual — as one flat module, compiled once as a single graph with static shapes (fullgraph=True, dynamic=False); the Fourier coefficients of the extension enter as constant buffers baked into the graph; (2) makes the training loop shape-static to match: the collocation pool (2142^14 points) is drawn once and only re-shuffled between epochs, mini-batches have a fixed size (2102^10; a remainder that would produce a ragged final batch is dropped), and the points carry no AD metadata; (3) compiles one further instance for the L-BFGS refinement phase, whose full-batch closure has a different (but again fixed) input shape; each compiled graph is reused for the entire phase – up to 1.6×1051.6× 10^5 Adam steps and 10410^4 full-batch L-BFGS iterations in a production run. The one-off cost of compilation is 1414–1818 seconds in our setting; amortised across a training run it is negligible, but it is per shape: an innocently varying batch size would pay it again and again, which is the practical reason for the static sampling design in item (2). It also explains a symptom worth knowing: a shape change mid-training manifests as a multi-second, log-silent stall that is easily mistaken for a hang. We add one experience report. Before settling on hand-rolled jets we used the framework’s composable functional transforms (torch.func’s jacfwd/jacrev) to obtain (u,J,H)(u,J,H); these already realise the forward-mode dimension count and run an order of magnitude faster than nested autograd (Table 1, row 2); for some time, they were the production pipeline. But composing them with the graph compiler proved fragile: in our environment (PyTorch 2.9), higher-order functional transforms over a function containing non-trivial input arithmetic could not be traced. The workaround was structural: replace transform-generated derivatives by the explicit jets of §4.2, so that the program presented to the compiler contains only elementary tensor operations. Relative to this eager transform pipeline (the state of the code immediately before compilation was introduced), the compiled jet evaluator is 2.92.9–4.9×4.9× faster (Table 1, rows 2 and 4): this, rather than the smaller margin over the eager jets, is the improvement that was actually experienced when the compiled graph entered the pipeline. The general lesson is that the intersection of advanced features (higher-order AD, functional transforms, compilers) is the least-charted part of any framework, and computations organised as plain tensor programs age better than ones that lean on that intersection. 4.4. Measured impact Table 1 isolates the two contributions on the production task of [18]: identical model (perturbed-trefoil boundary, default architecture, fixed random initialisation – the timings are insensitive to the weights), identical collocation points, identical hardware, double precision throughout; each timed step evaluates the residual on the batch, forms the mean-square loss, and computes the full parameter gradient, i.e. the work of one optimisation step, apart from the parameter update itself, for the two batch sizes used by the two training phases. The baseline row runs the project’s original nested reverse-mode implementation; before timing, we verified that all four pipelines return the same residual field, with maximum pairwise deviation 5.0×10−145.0× 10^-14 on a residual field of maximum modulus ≈53≈ 53. N=210N=2^10 (Adam) N=214N=2^14 (L-BFGS) Residual pipeline ms/step speedup ms/step speedup nested reverse-mode AD (§4.1) 427.5427.5 11 3260.73260.7 11 torch.func transforms, eager (§4.3) 42.042.0 10.2×10.2× 223.7223.7 14.6×14.6× forward 22-jets, eager (§4.2) 17.617.6 24.3×24.3× 143.9143.9 22.7×22.7× forward 22-jets, compiled (§4.3) 8.68.6 49.7×49.7× 78.278.2 41.7×41.7× Table 1. Median wall-clock time per training step (residual + loss + parameter gradient) over 2020 steps after warm-up, on a laptop CPU (Apple M4 Pro, 1010 threads, PyTorch 2.9.1, Python 3.12, float64); speedups are computed from the unrounded timings. The torch.func row reports the faster of the forward-over-reverse and forward-over-forward compositions at each batch size. One-off graph compilation costs 1414–1818 s per input shape (excluded: a one-off cost at the start of each training phase). The four pipelines agree pairwise to 5×10−145× 10^-14 on the same inputs. Let us conclude with three observations. First, the factors compose as claimed: moving from nested reverse-mode AD to the forward-mode functional transforms is worth 1010–15×15×, replacing the transforms by explicit jets a further 1.61.6–2.4×2.4×, and compiling the resulting graph a further 1.81.8–2.1×2.1×, for 4242–50×50× end to end. Second, the gains are broadly stable across a 1616-fold change in batch size, so they reflect the algorithm, not a small-batch overhead artefact. Third, the absolute numbers set the scale of experiments one can contemplate; at the production schedule (10410^4 Adam epochs of 1616 mini-batches plus up to 10410^4 full-batch L-BFGS iterations), the compiled pipeline completes the optimisation in roughly 4040 minutes of compute on a laptop (consistent with the roughly one hour reported in the companion paper) whereas the same schedule through the baseline pipeline extrapolates to more than a day. The entire study of [18], some dozens of trained models plus failed runs and ablations, was carried out on ordinary laptops; without these two optimisations it would have required a cluster, or would not have been done. 5. Guidelines We summarise the discussion as advice to a geometric analyst or differential geometer starting a PINN project. None of it is specific to minimal surfaces. (1) Parametrise constraints away; penalise only what you cannot parametrise. Build boundary conditions, positivity, decay and symmetry into the model: every constraint moved from the loss into the architecture removes a weighting hyperparameter, shrinks the search space, and holds identically throughout training. (2) Spend your regularity theory before you spend your compute. Whatever is known about solutions a priori (boundary asymptotics, decay rates, symmetries) is structure the network would otherwise have to learn, imperfectly, from the residual signal alone. (3) Make θ=0θ=0 meaningful. Arrange the architecture so the zero network is the best explicit approximation you can write down, and check that this initial map is smooth – the network will not fix a singular ansatz. (4) Count dimensions before differentiating. Reverse mode for the parameter gradient; forward mode for spatial derivatives whenever the domain dimension is small, as it is for curves and surfaces. For second-order residuals, propagating the 22-jet explicitly through the network — a few dozen lines of code implementing equations (4.1)–(4.2) – is exact, fast, and auditable. (5) Compile once; design for it. Keep shapes static (fixed collocation pool, fixed batch size), keep the residual a flat branch-free tensor program, and expect one compiled artefact per input shape. Treat a mysterious mid-run stall as a recompilation until proven otherwise. (6) Use double precision, and measure. Second derivatives in single precision floor the residual near 10−410^-4; float64 on a compiled CPU pipeline is a perfectly viable regime. Verify optimised pipelines against a naive reference implementation to round-off, and time the actual training step, not a proxy. 6. Concluding remarks The asymptotic Plateau problem is, in retrospect, an ideal proving ground for geometry-aware PINNs: the boundary condition sits at infinity, the boundary asymptotics of solutions are theorems, an exact model solution exists to anchor the ansatz, and the probed conjecture (signed counts of double points confronted with HOMFLY coefficients [5]) demands a precision that soft constraints and single precision cannot deliver. But nothing in the two theses of this note is specific to it: hard-constrained architectures apply whenever the constraint set can be parametrised, and the jet-propagation and compilation techniques of §4 apply verbatim to any PINN whose domain has low dimension – the typical situation for parametrised curves, surfaces and maps in geometric analysis. The efficiency margin of Table 1 is also what makes the natural next steps thinkable at laptop-to-workstation scale: higher-genus domains, systematic knot families, and, most importantly, pushing residuals towards the regime where the numerical solutions could serve as the starting point of computer-assisted existence proofs, in the spirit of [7, 21, 16]. We hope this note lowers the activation energy for colleagues in geometric analysis and differential geometry to attempt such computations themselves. References [1] M. T. Anderson (1982) Complete minimal varieties in hyperbolic space. Invent. Math. 69, p. 477–494 (English). External Links: ISSN 0020-9910, Document, Link Cited by: §1. [2] J. Ansel, E. Yang, H. He, N. Gimelshein, A. Jain, M. Voznesensky, B. Bao, P. Bell, D. Berard, E. Burovski, et al. (2024) PyTorch 2: faster machine learning through dynamic Python bytecode transformation and graph compilation. In Proceedings of the 29th ACM International Conference on Architectural Support for Programming Languages and Operating Systems, Volume 2 (ASPLOS ’24), p. 929–947. Cited by: §4.3. [3] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind (2018) Automatic differentiation in machine learning: a survey. Journal of Machine Learning Research 18 (153), p. 1–43. Cited by: §4.1. [4] G. Cortes, M. Esteban-Casadevall, Y. Feng, J. Henkel, E. Hirst, T. S. Gherardini, and A. G. Stapleton (2026) A Machine Learning Approach to the Nirenberg Problem. arXiv preprint arXiv:2602.12368. Cited by: §1. [5] J. Fine (2021) Knots, minimal surfaces and J-holomorphic curves. arXiv preprint arXiv:2112.07713. Cited by: §1, §6. [6] P. Freyd, D. Yetter, J. Hoste, W. B. R. Lickorish, K. Millett, and A. Ocneanu (1985) A new polynomial invariant of knots and links. Bull. Am. Math. Soc., New Ser. 12, p. 239–246 (English). External Links: ISSN 0273-0979, Document Cited by: §1. [7] J. Gomez-Serrano (2019) Computer-assisted proofs in PDE: a survey. SeMA Journal 76 (3), p. 459–484. Cited by: §1, §6. [8] A. Griewank and A. Walther (2008) Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation. 2nd edition, SIAM, Philadelphia, PA. External Links: ISBN 978-0-89871-659-7 Cited by: §4.1, §4.2. [9] J. Halverson and F. Ruehle (2024-10) Metric flows with neural networks. Machine Learning: Science and Technology 5 (4), p. 045020. External Links: Document, Link Cited by: §1. [10] K. Hashimoto, K. Kyo, M. Murata, G. Ogiwara, and N. Tanahashi (2026-01) Physics-informed neural network solves minimal surfaces in curved spacetime. Machine Learning: Science and Technology 7 (1), p. 015013. External Links: Document, Link Cited by: §1. [11] E. Hirst, H. N. S. Earp, and T. S. R. Silva (2026) Minimising Willmore Energy via Neural Flow. arXiv preprint arXiv:2604.04321. Cited by: §1. [12] E. Hirst, T. Schettini Gherardini, and A. G. Stapleton (2025) AInstein: Numerical Einstein Metrics via Machine Learning. AI Sci. 1 (2), p. 025001. External Links: 2502.13043, Document Cited by: §1. [13] D. P. Kingma and J. Ba (2015) Adam: A Method for Stochastic Optimization. Proceedings of ICLR. External Links: 1412.6980 Cited by: §3.2. [14] D. C. Liu and J. Nocedal (1989) On the limited memory BFGS method for large scale optimization. Mathematical Programming 45 (3), p. 503–528. External Links: Document Cited by: §3.2. [15] J. Marx-Kuo (2025) Variations of renormalized volume for minimal submanifolds of Poincaré-Einstein manifolds. Commun. Anal. Geom. 33 (1), p. 17–29 (English). External Links: ISSN 1019-8385, Document Cited by: §2.2. [16] D. Platt (2026) Non-uniqueness and symmetries for the Nirenberg problem using computer assistance. arXiv preprint arXiv:2603.29544. Cited by: §1, §6. [17] M. Raissi, P. Perdikaris, and G. E. Karniadakis (2019) Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J. Comput. Phys. 378, p. 686–707 (English). External Links: ISSN 0021-9991, Document, Link Cited by: §1. [18] T. Schettini Gherardini and M. Usula (2026) Minimal surfaces, knots, and neural networks. Note: arXiv preprint arXiv:2605.26234 External Links: 2605.26234, Link Cited by: §1, §1, §2.1, §2.2, Figure 3.1, Figure 3.1, §3.2, §3.2, §4.2, §4.4, §4.4, §4. [19] M. Usula (2026) Biharmonic maps between conformally compact manifolds. Calc. Var. Partial Differ. Equ. 65 (2), p. 24 (English). Note: Id/No 59 External Links: ISSN 0944-2669, Document Cited by: §2.1. [20] M. Usula (2026-05) Isometric embeddings of conformally compact manifolds into hyperbolic spaces. International Mathematics Research Notices 2026 (10), p. rnag098. External Links: ISSN 1073-7928, Document, Link, https://academic.oup.com/imrn/article-pdf/2026/10/rnag098/68350564/rnag098.pdf Cited by: §2.1. [21] Y. Wang, M. Bennani, J. Martens, S. Racanière, S. Blackwell, A. Matthews, S. Nikolov, G. Cao-Labora, D. S. Park, M. Arjovsky, D. Worrall, C. Qin, F. Alet, B. Kozlovskii, N. Tomašev, A. Davies, P. Kohli, T. Buckmaster, B. Georgiev, J. Gómez-Serrano, R. Jiang, and C. Lai (2025) Discovery of Unstable Singularities. arXiv preprint arXiv:2509.14185. Cited by: §1, §6. [22] S. Zhou and X. Ye (2023) Approximating High-Dimensional Minimal Surfaces with Physics-Informed Neural Networks. arXiv preprint arXiv:2309.02589. Cited by: §1.