Deep Legendre Transform
Authors: Aleksey Minabutdinov
Patrick Cheridito
Venue: NeurIPS 2025 poster ()
Keywords: Convex Conjugates
Legendre Transformation
Deep Learning
Input Convex Neural Networks
Primary area: optimization
TL;DR: We introduce a deep learning method for calculating convex conjugates in high-dimensional settings, significantly enhancing computational efficiency and providing L2-convergence certificates.

Abstract: We introduce a novel deep learning algorithm for computing convex conjugates of differentiable convex functions, a fundamental operation in convex analysis with various applications in different fields such as optimization, control theory, physics and economics. While traditional numerical methods suffer from the curse of dimensionality and become computationally intractable in high dimensions, more recent neural network-based approaches scale better, but have mostly been studied with the aim of solving optimal transport problems and require the solution of complicated optimization or max-min problems. Using an implicit Fenchel formulation of convex conjugation, our approach facilitates an efficient gradient-based framework for the minimization of approximation errors and, as a byproduct, also provides a posteriori estimates of the approximation accuracy. Numerical experiments demonstrate our method's ability to deliver accurate results across different high-dimensional examples. Moreover, by employing symbolic regression with Kolmogorov–Arnold networks, it is able to obtain the exact convex conjugates of specific convex functions.

Deep Legendre Transform
Aleksey Minabutdinoy Center of Economic Research and RiskLab ETH Zurich, Switzerland aminabut.dinov@eth.cch
Patrick Cheridito
Department of Mathematics and RiskLab ETH Zurich, Switzerland patrickc@ethz.com ch
Abstract
We introduce novel deep learning algorithm for computing convex conjugates differentiable convex functions, a fundamental operation convex analysis with various applications in different fields such optimization, control the- ory, physics and economics. While traditional numerical methods suffer from the curse of dimensionality and become computationally intractable in high di- mensions, more recent neural network-based approaches scale better, but have mostly been studied with the aim of solving optimal transport problems and quire the solution of complicated optimization or max-min problems. Using an implicit Fenchel formulation of convex conjugation, our approach facilitates an efficient gradient-based framework for the minimization of approximation eπrоrs and, as byproduct, also provides a posteriori estimates of the approximation accuracy. Numerical experiments demonstrate our method's ability to deliver accurate results across different high-dimensional examples. Moreover, by employ- ing symbolic regression with Kolmogorov-Amoot networks, it is able to obtain the exact convex conjugates of specific convex functions. Code available https: //gi thub. com http:/eeaaoo//Deeeeeeeeeeee -eraneee--eeTTaraaaorr Transform
Introduction
In this paper, we introduce DLT, a deep learning method for computing convex conjugates, which play an important role in different fields such convex analysis, optimization, control theory, physics, and economics.
The modern formulation of convex conjugation, also called Legendre-Fenchel transformation, was introduced by Fenchel 1949 and defines the convex conjugate function f:C - on nonempty subset C C RR by
(*(v) sup -(=))))(=)), yERd. zec
(1.1)
This definition does not need any special assumptions on the function f or the set C' and makes sense for all y Rdd Moreover, it directly follows from (I.I) that the conjugate function f*; Rd (-00, co] is lower semicontinuous and convex.
*Corresponding author
'The Legendre transform-a fundamental operation for switching between dual formulations and vari- ables has broad utility across thermodynamics, mechanies, optimization, and economics. Key applications include: switching between thermodynamic potentials by exchanging extensive variables (entropy, volume) for intensive conjugates (temperature, pressure) in physics; constructing Moreau envelopes variational analysis: deriving indirect utility and profit functions in economics; and computing convex potentials in optimal transport. Recent work has further demonstrated applications in approximate dynamic programming Shariti Kolarijani and Mohajerin Esfahani 2023) and self-concordant smoothing Adeoye and Bemporad (203). This paper introduces the Deep Legendre Transform algorithm, establishes its theoretical loundations, and demonstrates its application to optimal-control problems governed by Hamilton-Jacohi equations.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

In the special case, where f:C + R is differentiable convex function on an open convex set C' S R. one has
f (()) f((F))
(1.2)
for all y E D VF(C). The right side of (12)). called Legendre transform, was used by Legendre [187] when studying the minimal surface problem. Note that for it to be well-defined, the gradient D does not have to be invertible. In fact, can easily be shown that if f is
map ff:C
differentiable and convex, the right side of (I.]) is maximized by any
C' satisfying ff(z)) y. So the preimage (Vf)) (y) in (.2) can be substituted by any point I E (Vf) (v).
Since S" is not always available in closed form, various numerical methods have been developed to approximate it numerically. Most of the classical approaches use the modern definition (1.1) and
approximate (y) for y in a finite subset y C IR: by the discrete Legendre-Fenchel transform
max (zz y) f(x)) ZEX
(1.3)
based on a finite subset X' C C. Of course, for 1.3 to provide a satisfactory approximation to f. the subsets X' and y need to be sufficiently dense. Typically. X' and y are chosen as grids ... x V' of similar size. For instance, for regular grids of the form
X' = X1 x x Xa and y = N 2 N...., 1]d, conditions under which the discrete Legendre-Fenchee
transform converges to f have been given by Corrias [996]. In this case, naive computation of for all y E y. has complexity O(N2d). For the case where N is a power of 2, Brenier [1989 and Noullez and Vergassola [1994] have derived fast algorithms to compute (I.3), similar in spirit to the fast Fourier transform of Cooley and Tukey| [1655. with complexity O((N log N) By first passing to convex hulls, Lucet [199]] has been able to design an algorithm with complexity -(X
X = -=y-
which for X' = 1% mcans O(N4). But since this is still exponential in d, it essentially limits the 1, 2, and 4.
method to dimensions
A different approach, going back to a tweet by Peyré [2020], uses the fact that if C' S R is bounded, [
f can be approximated with the softmax convex conjugate
(v( = elog (focp)((r))))
(x. y)
f(x)
dr
E
for a small E 0, which can be expressed via Gaussian convolution as
1
(2 = Q:' (Q.(F)=() = Ge (v).
where
G = exp (-11-11 1-11 Qe(f) = exp
1
(()) and Q 1='(F)=111-1)) - F-). 11:12 f()
2E
Approximating the convolution (f) * Ge with a sum over finite set X' C C' yields a GPU-friendly algorithm with complexity O(|x]oog] (: see Blondel and Roulet [2024] for details. But since it is grid-based, it still only gives accurate results in low-dimensional settings.
Lucet et al. 2009 have explored the computation of convex conjugates of bivariate piecewise linear quadratic convex functions, and Haque and Lucet (2018) have proposed an algorithm for such functions with linear worst-case complexity. But the scalability to higher dimensions remains
challenge.
Recently, numerical methods for convex conjugation have also been studied in the machine learning literature as a tool solve the Kantorovich dual of optimal transport problems and develop efficient Wasserstein GANs. While their performance has mostly been tested to see how well they help to solve transport problems, promising results have been obtained by training neural networks to predict an optimal ar in the Legendre- transform (I.I). F.g. [Dam et al. [209) and Taghvaei and Jalali 2019 predict from y/ with a neural network hop (v), where he is maximizing the right side of Makkuva et al. [020 do the same using predictors of the form h Vlg for an ICNN
'The entropic log-sur pog-sum-eae sum-exp smoothing a max can be traced back to the Gibbs variational principle; brought into modern convex-opiimistiion analysis by Nesterov (2055
Input Convex Neural Network; see Amos al. 2017

---

Lg. Korotin et al. 2021a and Amos 2023 have proposed to predict an optimal z in 1.1 with
neural network h by minimizing the first order condition
Σ (VT(()(())) vill2
for sample points yi D,
which, in effect, amounts to learning generalized inverse of the gradient map Vf:CC + D; see
In this paper, we employ the following implicit Legendre formulation of convex conjugation:
('(V[(=)) (a, vs(x)) [(x), x C
(1.4)
for differentiable convex functions f: C:C' - R defined on an open convex set C' C Rd. This allows us to approximate f*: D Vf(C) - I by training machine learning model go: D R parametrized by a parameter 9 from set such that
go Vf(x)) Vf())) f(x) for all C.
DLT implements this by solving
min DE
Σ
([o[[()))++(0 (z. v/())))
(1.5)
for training set Xrrain S C. This combines the following advantages:
(i) On the training points I E Rusins thanks to 1.4). J. the algorithm's target values (x, Vf(x)) fir) exactly agree with the true values f (V)()))
(ii) If is chosen as a neural network, the approach readily scales to tups, and if combined with an ICNN as approximator, is guaranteed to output a convex approximation to the true Legendre transform, which is known to be convex.
se-
(iii) Specifying De as a KAN makes it possible to obtain exact solutions for certain convex functions.
(iv) The implicit Fenchel formulation
1.4 can also be used to derive posteriori error estimates.
(i) is a theoretical advantage compared to existing approaches and, as illustrated in Section helow, yields good numerical results, in fact, limited only by the used NN architecture/computing power available. We demonstrate that DLT achieves state-of-the-an performance by comparing it to the direct learning of the convex conjugate (i.e., using f* target function in cases where it is available in closed form).
(ii) ICNNs have been introduced by Amos al. 2017 and it has been shown by Chen et al. 2019
that they can approximate any convex function arbitrarily well on compacts. With neural networks as approximators, the minimization problem (1.5) numerically be solved in high dimensions with stochastic gradient descent algorithm. The choice of ICNNs ensures that the approximation is convex. This regularizes the learning problem and helps against overfitting. Moreover, is helpful in various downstream tasks, such as optimization, control, or optimal transport, where convex conjugation be utilized powerful tool to find optimal solutions. Alternatively, (0166), which in many multidimensional cases provides a better approximation than ICNN but does not guarantee convexity.
ResNet (He et al.
(iii) KANs have recently been introduced by Liu al. (025). They are particularly powerful for approximation tasks requiring high precision. But in addition, they also be used for symbolic regression to recover functions that have a closed form expressions in terms of special functions such polynomials, exp, log, sin cos. In Section 5 below, we show that KANs can compute Legendre transforms certain convex functions exactly.
(iv) As byproduct, our approach yields posteriori I, appapproxintaioon guarantees. Indeed, by testing a candidate approximating function g: - IR against the implicit Legendre transform (1.4) on test set Xrest. C C sampled from distribution μr, it is possible to derive estimates of the L'-approximation error l|g f 10 (D()).)) with respect to the push-forward
μ (Vf)-1 of under the gradient
KKllmgorov Arnold Network; see Liu 2025

---

map Vf:C + D. This is valuable beyond the framework of this paper since it makes it possible to evaluate the performance of any numerical method for convex conjugation.
The training and test sets Xrain, Nest C used for training and testing a candidate approximator Go of the convex conjugate f can be generated by sampling from arbitrary distribution j on C. But the gradient map Vf:C + D is distorting C' and transforms any finite subset X' S C consist of independent realizations of a desired target distribution on D = Vf(C), X' has to be sampled from a distribution whose push-forward o (Vf))² 1 is similar to v. We do this by training a neural network hg: D - R that approximates generalized inverse of V C + D. Then target training and test sets Verain; Vrest C D can be generated according to a desired distribution v on D, and the sets Xtrain = -ha (JVerin))) and Xeest hop (Nues) can be used for training and testing go-
into fx()) S D. If one wants (X())
The remainder of the paper is organized as follows: we derive estimates of the 1a-approximation error of any candidate approximating function of with respect to given probability measure y on the image = VC(CC) of C' under the gradient map ff. In Section 4 we learn a generalized inverse mapping of V that improves the performance of DLT if f is heavily distorting and is needed to estimate the L L'-approximation error with respect to desired target distribution experiments. Sec 6 concludes. Proofs of theoretical results and additional numerical experiments are given in the Appendix.
In Section 2 we introduce DLT. In Section 3
D. In Section we report the results of different numerical
2 Convex conjugation as machine learning problem
Let f: + IR be a differentiable convex function defined on an open convex subset C C Rd Our goal is to approximate the convex conjugate f* Fd + (-o0, oo] on the image D Vf(C) of C under the gradient map V.-
By the Fenchel- Young inequality one has
f(x) + (w) 2 (x, for all x E and y E RRd
with equality if and only if y = ffx). As a result, one obtains that on D, the convex conjugate f" coincides with the Legendre transform
(v) ((V.)) (v)
EED,
(vv) can be substituted by any C' satisfying Vf(x) y. However, while in many situations, the function ff and its gradient V are readily accessible, the inverse gradient map (V f) might not exist, or if it does, it is often not be available in closed form. Therefore, we use the implicit Legendre representation:
where the preimage
('(V((x)) = (z, ((()) f(x), I E
(2.1)
which makes it possible to formulate the numerical approximation of f as machine learning problem. More precisely, for any machine learning model go: D - IR with parameter 0 set θ DLT approximates I on D C Rd by minimizing the empirical squared error
Σ (go( (v)())())) + f(x) - (z. v)(=)))²
(2.2)
over 8 E O In principle, Xiaain can be any finite set of training data points x E C. But in Section 4 below, describe how it can be assembled such that the image f((rrai)) is approximately distributed according to given target distribution v on D = VS(C).
Any numerical solution E θ of the minimization problem (2.2) yields an approximation gj: D - R of the true conjugate function f D + R. The question is how accurate the approximation is.
3 A posteriori approximation error estimates
Thanks to the implicit form of the Legendre transform 2.1 1 it is possible to estimate the L
approximation error of any candidate approximator g: D R of f - R with standard
seee e.g.. Fenchel 1949 or Rockafellar 1970

---

Monte Carlo average even if the true f is not known. Let &rest be test set of random points z E C independently sampled from a probability distribution a over C. Then, due (2.I). one has
Σ + (z, vf()))) Rtest
Σ ( wEA.
Reest
Z 2
which, by the law of large numbers, for Rrest + oo, converges to
[c ( m(dx),
which, in turn, equals the squared I,2 ..-appooxination error
I (olv) 1" (0) = llg f (
(3.1)
with respect to the push-forward 11 (Vf) of the measure μa. This result lets report Monte Carlo estimates of the error llg ff (' as if f were known in closed form.
Sampling from a given distribution in gradient space
The most straight-foree forward way to generate training and test sets Xtrain Mest C C is to sample independent random points according to convenient distribution on C.
But if one wants to minimize the L a'-apppoxinatiion error (3.1 with respect to given distribution on D = V f(C), has to be chosen so that its push-forward o (Vf) equals v. For instance, the goal is a close approximation everywhere on D for bounded subset D S R e.g. a ball or hypercube in Rd one might want to be the uniform distribution D. If D is unbounded, e.g.
D = RRd, one might want v to be a multivariate normal distribution. But if one interested close approximation of f in neighborhood of particular point yo E D, is better to choose distribution v that is concentrated around yo-
To sample points in C such that the distribution of V f(x) inverse of the gradient map Vf: C + D. This can be done by training neural network ha: D - Rd with parameter R2 such that either
) is close to learn generalized
ha Vf(x) R for E (4. or
V o ho(p)) for E D (4.1b)
We choose standard feed-forward neural network he with hidden layers and GELU-activation. To achieve 4.1a or (4.1b), one can either minimize
Σ (ha Vf(x) x)² (4.2a)
Σ (Vfoho()) v)² (4.2b)
or
neWarain
for suitable training Wirain C Verain C D, respectively. 4.2a and 442b both have advantages and disadvantages. Therefore, we combine the two.
We first pretrain ha by minimizing 4.2a) This works points y from D to Id C. Wusin can e.g. be sampled from uniform distribution over if C bounded or multivariate normal distribution if C' is unbounded. If V is not invertible, (4.2a) can typically not be reduced to zero. But pretraining will force h.s to map most of D into C. After pretraining, minimize convex combination of (4.2a) and (4.2b) for Verain independently sampled from v. In every training step, points y E Veraiu for which he ) (v) is not C are omitted from training. approximate be stopped if (4.2b) is sufficiently close to zero, yielding
if the initialized network hes maps some
The training procedure generalized inverse hg of V.. If Viest. is a test set independently sampled of ) Devain according to v, then Xtrain C h (Verain) and Xiest. C (Juss) can be used train and test given h3
machine learning model da to approximate f".
'Confidence intervals of the estimator provided in Lemma A.l Appendix

---

Numerical experiments
In this section, we report results of numerical approximation of the convex conjugates of differen- tiable convex functions using DLT with MLPs, ResNets, ICNNs, and KANs. As an application, we use DLT approximate solutions of Hamilton-Jase PDEs.
Numerical approximations with MLPs, ResNets and ICNNS
We first approximate f": + R with different neural networks such as standard multi-layer perceptrons (MLPs), residual networks (ResNets) and input convex neural networks (ICNNs). We use MLPs with two hidden layers containing 128 units and GELU activation. Our ResNets have two residual blocks, each containing two dense layers, resulting in 4-layer architecture with skip connections and GELU activation. The employed MLP-ICNNs have the same architecture as the MLPs but constrain hidden weights to non-negative and apply softplus activations to ensure convexity. The ICNNs are specified as in Amos et al. [2017]. using softplus activation and non- negative weights in two hidden layers along with input skip connections. Parameters were initialized according to Gaussian distribution, but for ICNN, the weights were squared to ensure non -ongetiatiiity For training, we used the Adam optimize
As first test, we compare DLT in the three cases of Table 1 where f is known analytically, to directly learning f using pairs (y:, f (1:)) as targets.
Table 1: Functions f. domains, Legendre transforms, and dual images D Vf(C).
Function
f(x)
Domain C
(()
D D= VF(C)
Quadratic
NN 1)a 233(
$ Σ vî
1)
Neg-L.og
- - log(s;)
d Σ Elog(-ws)
{-1/2 x C} R
-0.0,-
Neg-Entropy Σ log
Σ explui 1)
(log + E
- 1.3,3.3)
Table shows that DLT matches the performance of direct learning both in terms of approximation accuracy and training time. This empirical finding is maintained across different network architectures and dimensions, confirming that the performance of DLT is only constrained by the neural networks' our
capacity to approximate multidimensionall functions and not by its implicit formulation. experiments, ResNet provided the best approximations despite not guaranteeing convexity, while MLP-ICNN showed limitations due to the lack of skip connections compared to ICNN
Functions without closed form conjugate Most functions do admit closed form expression for the convex conjugate. However, DLT gives the same a posteriori error estimates as if I were known explicitly and learned directly on D. For illustration, we consider the quadratic over linear function
2
f(x)
31
It +1
not defined everywhere on Rd. But be shown that it is convex on the convex C {a E Rd 11 + +.... Fa > 0): see e.g. Boyd and Vandenberghe (2004). The i-th component of the gradient V [((x)) is
2 + 2x, + Ej (2x, (2x, 2z))
(5.1)
(x)
Fdd
1)2
But the convex conjugate f is not known in closed form, and the image D VS(C) is non-trivial subset of Ed Even though we do not know f explicitly, DLT yields accurate ICNN approximations of f on D. Here, we used an ICNN with hidden layers and ELU activations, we sampled random
7 All computations were performed on NVIDIA T4 GPU and Intel Core Xeon CPUs running Python 3.11.12. More experiments. including the algorithms by Lucet's and Peyré, are available in the supplementary material.
see Kingma and 2017
Additional mulitdmensional experiments confirming that DLT has the same performance as direct learning of known conjugates in terms of approximation accuracy and training time given in the Appendix.

---

Table Comparison of DLT with direct learning (Dir) for d 10 and 50. Errors were computed with Monte Carlo on test set of 4096 i.id. random points sampled from Unif (Vf)
10
50
Function
Model
RMSE RMSE
Aerain (8)
RMSEpur RRSS
ferain (8)
MLP
2.96c-5/2.86e:
86/84
3.55e-3/3.59e.
88/87
MLPPICCN
6.21e+0/3.15e4
87/77
6.46e+2/6.446422
87/77
Quadratic
ResNet
2.43e-5/3.00es 107/101
2.92e-3/3.25c
124/120
ICNN
1.79e-2/1.23e-2
93/93
7.18e-2/6.6664
97/98
MLP
6.75c-4/8.438
94/94
3.81c-1/2.81e-
98/98
MLP-ICNN
4.23c+0/4.22 2,0++0
64/62
2.03c+1/2.03er
67/70
Neg. Log
ResNet
7.34c-4/9.22e-4
109/111
L.05e-11.23ee-1
134/131
ICNN
6.47e+0/4.22e+4
54/65
2.03e+1/2.06+11
74/62
MLP
1.06c-3/12262
96/100
2.71c-1/3.91c-1
100/101
MLP-ICNN
4.50e40/4.01e+0
93/101
7.83e+177.87++1
98/105
Neg. Ent.
ResNet
6.62e-4/7.41e-
110/116
6.93e-217.2ee-2
133/135
ICNN
1.31c-2/1.74e-2
99/107
9.39c+119.39eH1
60/69
d
1190 f'1 L2
Time (mm:ss)
2
6.38e-04
4:14
2.39c-03 5
5:14
10
1.80c-02
6:55
20
3.67e-02
14:28
50
6.00e-02
29:19
100
3.68e-02
38:14
Figure 1: Quadratic- ceee--vvee---neeneee benchmark. Left: Sample points from 2D standard half-normal distribution and their images under the gradient map f (5.I). Right: L errors and training times for various dimensions.
points from standard half- half-normal distribution C. Figure depicts these random points and their
images under the gradient map V given in (5.1) for d = 2. The table on the right shows estimates of the L'-aaproximation error ||9a - 5 LP(R(..) together with the times in minutes and seconds
took to train the approximations Og for different Note that the true f* is unknown here, but due to the results of Section is still possible to provide L2 approximation errors.
In Appendix B we study the following additional non-separable examples: (i) quadratic f(x)
2 Qx: (ii) exponentialmminus-lnes f(x) = e(a.a) - (b, x); and (iii) ICNN-represented functions.
Comparison with classical grid methods
Next, we compare DLT against the LLT algorithm of Lucet 1997]. a state-of-the-art grid method for computing convex conjugates. Rather than directly evaluating max (xr, y) f()) over full d-dimensional grid X' of size Ad (which requires (2) operations), Lucet reformulates the d-dimensional maximization as a sequence of nested 1D problems:
(yI Wa) max DI
max y2 + ... + max Lyaza- (vnza-
This decomposition reduces the computational complexity to OddN++1) while still requiring O(N4) memory. Table comppres the methods for fixed discretization of N = 10 points per dimension for two of the functions from Table [ All errors were evaluated on the same random sample
Durt D VS(C) against the true f For low dimensions (d E {2, 6)), Lucet runs quickly with modest memory, making possible to use much finer grids within a reasonable time budget. However, d grows beyond 8, even modest grids become memory time infeasible due to O

---

Table Comparison of Lucet with DLT for dimensions d E (2,6,8 6, 8, 10). Lucet uses N 10 points dim; DLT uses a ResNet (128, 128). Reported: solution time Eolvee evaluation time Evall active memory (MB) and RMSE.
Function
d
Method
Model
Luoive (s)
Leval (s)
MB(act)
RMSE
Lucet
(N4- Grid, N=10)
0.00
0.12
0.0
3.65e-01
DLT
(ResNet, 128,128)
16.97
0.02
2.14e-02
Neg. Log
Lucet
(N Grid, N=10)
15.17
1.02
15.3
1.83e+000
DLT
(ResNet, 128,128)
17.90
0.02
1.1
8.11c-02
Lucet
(Né_ Grid, (Ni-Grid, N N110)
1.96e++03
3.96
1.53e+ 03
2.93e+001
DLT
(ResNet, 128,128,
26.22
0.02
1.1
1.33c-01
10
DLT
(ResNet, 128,128)
32.09
0.02
1.32c-01 1.1
Lucet
- N N=10) 10)
0.00
0.12
0.0
1.42e 01
DLT
(ResNet, 128,128)
17.13
2,02e
Neg. Entropy
Lucet
(N' Grid, N N=10)
15.11
1.01
15.3
7.32e+001
DLT
(ResNet, 128, 128)
18.74
0.02
1.1
7.99--02
Lucet
(N' -rid, N ==10)
1.94e+003
3.98
1.53e 03
1,08e
DLT
(ResNet, 128,128)
27.16
0.02
1.1
1.40e-01
10
DLT
(ResNet, 128,128)
32.51
0.02
1.1
4.47e-01
storage and Od Nd+11) computational complexity. In contrast, DLT maintains constant memory footprint and scales smoothly to dimensions d beyond 9.
Direct vs. inverse sampling: distorting gradient mapping
When the gradient mapping ff is significantly distorting, the image V (2²) of a uniformly sampled training set C C will highly non-uniformly distributed on This skewed distribution can adversely affect learning dynamics, concentrating approximation accuracy in certain regions while neglecting others. The approximate inverse sampling method introduced in Section 4 addresses this sample from any desired distribution (e.g.. uniform) directly D.
issue by allowing
As an illustration, we consider the multidimensional negative logarithm f(x) = ( - log(a;) on C (10-3 10 10- With V fzx) (-1/21.... //////........aa) the gradient map severely distorts C, transforming it to D Vf(C) (-1000, 10) The known closed-form of the convex conjugate (v) Ei log(-y> d.
First, we train with points sampled uniformly from C (direct sampling). Then, we compare with our inverse sampling approach from Section 4, using Merain (Mraan)))) where Viaain is sampled uniformly from D and h approximates the inverse of V.. Table shows the 2 approximation errors and training times for both methods.
The inverse sampling method consistently produces better approximations (1-2 orders of magnitude improvement for > 10) the cost of increased training time. Figure illustrates this improvement for d = 2, where properly addressing the gradient distortion leads to a more uniform distribution of approximation accuracy across D.
Applications to Hamilton-Jacohi equations
The Hamilton-Jacobi (HJ) equation a comeerstone of classical mechanics, optimal control theory, and mathematical physics, which provides natural application domain for DLT.
Consider the Cauchy problem for the HJ equation:
( H(V ((V,((z,))) =
(z,t) Rd (0, o0).
(5.2)
(u(x, 0) g
x RR
for two convex functions H, RRd R.
The Hopf formula НHpp 1965 provides the following expression of u as a double convex conjugate:
u(z,t)= (g" + +tH) (x),
(5.3)

---

Sampling (V1) in D
(1/)))
550
(y. (b-)((b))
550
Scatter Plot
Numencal and Analytical OUT slong axis
Numenica and Analytical DLT along axts
Tue Conjugete
Tue Conjugete
Without Appr.
V2
Y1
y1
y2
Scatter Plot of sampling
Numerical and Analytical DET along axis
Numnerical and Anelytical DCT along axts
Tue Conjugate
Conjuugate
With Appr-
V1
V1
y2
Figure 2: Negative Logarithm: Numerical estimate of the Legendre transform computed with and without inverse sampling compared to the true Legendre transform in two dimensions.
Table 4: Comparison of direct and inverse sampling for the negative logarithm. Inverse sampling uses an approximate inverse map he whose quality measured by 2/(svd))lV/ he id p., 990 is the sidc-length of D and total training time split between ha and 9A-
Direct Sampling
Inverse Sampling
19a
1a
Time
Inv. Map Error
[1ga
--- 1r.
Time ha
Time
Total
2
3.18e-02
12:50
1.32e-03
5.15e-03
4:06
22:15
26:21
5
2.47e-01
13:12
1.72e-03
3.45e-02
6:42
25:45
32:27
10
3.18c+00
14:41
2.17e-03
4.07c-02
9:16
31:22
40:38
20
4.31e+00
15:24
2.65e-03
7.72e-02
13:05
35:14
48:19
50
9.46e+00
17:51
7.98e-03
3.80e-01
23:31
44:24
67:55
5,44e+01
26:10
3.75e-02
3.06e+00
34:27
62:43
97:10
While this formula is explicit, evaluating it in high dimensions remains computationally challenging. To address this, we propose Time-DLT Time-Parameterized Deep Legendre Transform) which approximates (g* + (H)* using DLT with neural network up (x., t) with input variables and t.
Experimental results. We study HJ equation (.2) with H(x) g(x) - 1 x7.
has the explicit solution u(x, 1) Σd of which we can use a reference. use Time-DLT to compute solution C (-a, a²da T for T 2 and compare against the DGM (Deep Galerkin Method) of Sirignano and Spiliopoulos 2018). both implemented using the same ResNet architectures with two residual blocks. DGM minimizes the residual of the PDE C directly, while Time-DLT approximates the Hopf formula (5.3). (able55sshov that Time-DLT outperforms DGM in approximating the true solution across different dimensions and times (however, we used the explicit expression for g" and just computed the outer Legendre transform).
2 2(++)
Near 0, DGM performs better by directly enforcing the initial condition. But for larger the Time-DI LLt solutions are closer to the correct ones indicated by DGM/DLT I error ratio > 1. Note that DGM aims minimize PDE residuals, but small PDE residuals do not always correspond to small approximation errors.
We repeat the approximations 10 times to statistically significant results and report estimated std. deviations a.

---

Table 5: Comparison of DGM and Time-DLT for different dimensions and times.
d
DGM
Time-DLT
Err. Ratio
L Error
PDE Res.
time (s)
L Error
PDE Res.
time (s)
p
10
0.1
8.90e-2
4.82e-1
82.40
8.46e-2
1.60e40
58.76
0.95
0.09
0.5
2.14e-1
3,04e-1
5.25e-2
9.17e-1
3.69
0,44
1.0
2.26e-1
188e-1
3.68e-2
5.12e-1
7.32
1.30
2.0
1.64e-1
8.05e-2
2.72e-2
2.27e-1
8.64
2.57
30
0.1
3.45e-1
1.18++0
87.18
4.29e-1
2.69e+0
192.62
0.83
0.12
0.5
6,40e-1
8.85e-1
2.60e-1
1.82e+0
2.63
0.34
1.0
1.0l++0
4.68e-1
1.90e-1
9,69e-1
5.15
1.23
2.0
1.21c+0
2.56c-1
1.28c-1
3,98c-1
9.00
2.21
Symbolic regression with KANs
Kolmogorov-Amoll networks (KANs) have recently been introduced by Liu et al. 2025 Inspired by the Kolmogorov- Anrol representation theorem, KANS hold great potential in approximation tasks requiring high precision. Moreover, in cases, where the Legendre transform has closed form expression in terms of special functions, KANs can be employed obtain the exact solution using symbolic regression. This is significant advantage compared to numerical approximation since it uncovers the underlying mathematical structure. To illustrate symbolic regression with KANs, we consider the three functions of Table 1 in two dimensions.
In all three cases we performed a symbolic regression with KANs to minimize the empirical squared error (2.2) for training set Arrain sampled uniformly from the square [1, 10j We used polynomials up the third order, exp and log as basis functions. As can be seen from the results Table KANs
are able to find the correct symbolic expressions for the Legendre transforms of these functions with the right constants (up to rounding errors).
Table 6: DLT with KANs: two-dimensional training regions, L and obtained symbolic expressions.
Function
C
D
1.2 Error
Symbolic Expression
x
10]2
[2, 20)²
0.25yř + 0.25yž
log(z;)
[0]
[-1, -011
5.42e-8
1.0 1o00(-0.000) 0.70)) 1.0 (00g(-))09))) 0.56
D 1== v. log(x;)
[1,10
[1,3.3]
6.78c-05
0.37 exp(1.0p)) + 0.37 exp(1.0y2)
Conclusion, limitation and future work
In this paper, we have introduced DLT (Deep Legendre Transform), a deep learning method for computing convex conjugates of differentiable convex functions f: C' R defined on an open convex set C' Rd Our approach provides posteriori approximation error estimates and scales effectively to high- high-dimensional settings, outperforming traditional numerical methods.
The limitation of DLT is that it needs function C:C + R to be differentiable and convex. Moreover, it only approximates the convex conjugate f*, which can be defined on all of IR. on the image D Vf(C), which is contained in, but in general not equal to, dom f {y E Rd (v) < xc).
The extension of the approach non-differentiable convex functions C IR and approximation of f on dom / will be studied in future work.
Acknowledgments and Disclosure of Funding
This work was supported by the Swiss National Science Foundation under Grant No. 10003723. We grateful to the anonymous reviewers for their constructive comments and suggestions.
10

---

References
A. D. Adeoye and A. Bemporad. Self-concordant smoothing for large-scale convex composite optimization, 2023. arXiv preprint arXiv:2309.01781.
B. Amos. On amortizing convex conjugates for optimal transport. Proceedings of the Eleventh International Conference on Learning Representations (ICLR), 2023.
B. Amos, L.. Xu, and Z. Kolter. Input convex neural networks. In Proceedings of the 34th International Conference on Machine Learning (ICML), volume 70, pages 146-155, PMLR, 2017.
M. Blondel and V. Roulet. The elements of differentiable programming. 2024. URL https
7/ //axiiv... org/ aba/ 2403 14606
S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press, 2004.
Y. Brenier. Un algorithme rapide pour le calcul de transformées de Fenchel discrêtes. C. R. Acad. Sci. Paris. Sér. Math. 308:587-589, 1989.
Y. Chen, Y. Shi, and B. Zhang. Optimal control via neural networks: a convex approach. In Proceedings of the 7th International Conference on Learning Representations (ICML), volume 2019.
W. Cooley and J. W. Tukey, An algorithm for the machine calculation of complex Fourier series. Mathematics of Computation, 90:297-301, 1965.
L. Corrias. Fast Aegendre-Fenchar transform and applications to Hamilton-Jacobi equations and conservation laws. SIAM Journal on Numerical Analysis, 33:1534-1558, 1996.
N. Dam, Q. Hoang, T. Le, T. Nguyen, H. Bui, and D. Phung. Three-player Wasserstein GAN via amortised duality. Proceedings of the Twenty-Fighth International Joint Conference on Artificial Intelligence (AAAI), pages 2202-2208 2019.
W. Fenchel. On conjugate convex functions. Canadian Journal of Mathematics, 1(1):73-77, 1949.
Haque and Y. Lucet. A linear-time algorithm to compute the conjugate of convex piecewise linear-quadratic bivariate functions. Computational Optimization and Applications, 70(2):593-613, 2018.
K. He, X. Zhang, S. Ren, and J. Sun. Deep residual learning for image recognition. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 770-778, 2016.
J.BB. Urruttt Lipschitz r-continuity of the approximative subdifferential convex function. Mathematica Scandinavica, 47:123-133 134, 1980.
E. Hopf. Generalized solutions of nonlinear equations of the first order. Journal of Mathematical Mechanics, 14:951-973, 1965.
D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv: 1412,6980, 2017.
A. Korotin, V. Egiazarian, A. Asadulaev, A. Safin, and E. Burnaev. Wasserstein-2 generative networks. In International Conference on Learning Representations (ICLR), 2021a
A. Korotin, L. Li, A. Genevay, J. Solomon, A. Filippov, and E. Burnaev. Do neural optimal transport solvers work? A continuous Wasserstein-2 benchmark. In Advances in Neural Information Processing Systems (NeurIPS), 2021b.
A.-M. Legendre. Mémoire I'intégration de quelques équations aux différences partielles. Mé- moires de I' l'Académie Royale des Sciences, 1787.
Z. Liu, Y. Wang, S. Vaidya, F. Ruehle, J. Halverson, M. Soljacic, T. Y. Hou, and M. Tegmark. KAN: Kolmogorov-Arnold networks. In The Thirteenth International Conference on Learning Representations (ICLR), 2025.

---

Y. Lucet. Faster than the Fast Legendre Transform, the Linear-Tinee Legendre Transform. Numerical Algorithms, 16(2):171-18, 1997.
Y. Lucet, H. Bauschke, and M. Trienis. The piecewise linear-quadratic model for computational convex analysis. Computational Optimization and Applications, 43(1):95-118, 2009.
A. Makkuva, A. Taghvaei, S. Oh, and J. Lee. Optimal transport mapping via input convex neural networks. Proceedings of the 37th International Conference on Machine Learning (ICML), volume 119 of Proceedings of Machine Learning Research, pages 6672 6681. PMLR, 2020.
Y. Nesterov. Smooth minimization of non-smooth functions. Mathematical Programming, 103(1): 127-152, 2005.
A. Noullez and M. Vergassola. A fast Legendre transform algorithm and applications to the adhesion model. Journal of Scientific Computing, 9:259-281 1994.
G. Peyré. The Legendre transform can be approximated O(n log(n)) using a Gaussian convolution. X, 2020.
R. T. Rockafellar. Convex Analysis. Princeton Mathematical Series. Princeton University Press, Princeton, J., 1970.
M. A. Sharifi Kolarijani and P. Mohajerin Esfahani. Fast approximate dynamic programming for input-aftine dynamics. IEEE Transactions on Automatic Control, 68(10))631-6332,, 2023.
Sirignano and K. Spiliopoulos. DGM: deep learning algorithm for solving partial differential equations. Journal of Computational Physics, 375:1339-1336 2018.
A. Taghvaei and A. Jalali. 2-Wasserstein approximation via restricted convex potentials with applica- tion to improved training for GANs, arXiv preprint arXiv: 1902.07977 2019.
12

---

A Unbiased a posteriori error estimator
Lemma A.1 (Unbiased a posteriori error estimator), Let f:C - R be differentiable convex function on an open convex set C C R and denote VCC).. Fix a probability measure p on C and consider the push-forward y/ = oo (Vf) on D. For a measurable function D - R, denote
(g)
((VV((X)) + (X Vf(x)))
where X1, Xn are i.dd random vectors with distribution Then the following hold:
(i gfg f e L3),,v), one has
E (o) = =f, (o(v)- (()) u(dy) = 1-11--f" ( (Unbiasedness)
and En (g) - 1g --M1EE almost surely for n - oc.
i Ifg f E (D,,,,), then
Var (En(0) o²
for = ( ( VX(I)))
and vn (g) و1g f (B.(())))) N(0, o²) in distribution as co. In particular, if o is an estimate of a and za denotes the a-quantile of (O, 1), then
Eg(g)
2-a022
Vn
is an approximate (1 - -c-oofieeeeece interval for l|g - 8 M'
Proof. (i) By the Fenchel-Young equality (*(VS(x)) (x, V.(x)) (xx), one has
((V](())) f(x)) (X,, Vf(X))) = (o((f(())) - (V
So, taking expectation and using the fush-forward measure == o (Vf) -1 yields
E [(E (o(V)(())) - -f (Vf(()))) ((dr) (/ (9(v) f (v)) v(dy) =lgg- - (
and it follows from the law of large numbers that En(g) + 1[g f (ILAAC,,))) almost surely. This shows (i).
(ii) Denote
Z: (D(VI(X)) (X, VS(X)))
Then (g) 1 Σm Z, and therefore, 71
Var (E)(e)
Var(Z,)
d²
n
So, follows from the central limit theorem that vn distribution. In particular, En Vn 21-/2 an approximate co-ooddeee interval for llg - ( which concludes the proof.
- |gg (1
1 (1: - NN0, in
13

---

B
Technical appendix: additional experiments
Content
Sections B.1 and B.2 describe low-dimensional experiments with grid-based methods: Lucet's s Linear Time Legendre Transform and Peyré's Gaussian convolution approach with Monte Carlo integration.
Section B.3 presents comprehensive benchmarks comparing our DLT approach against direct super- vised learning when analytical conjugates are known. Using identical neural architectures across multiple functions and dimensions d .2..... 200), demonstrate that DLT achieves equivalent L errors without requiring closed-form Legendre transforms.
Section B.4 evaluates DLT on challenging non-separable convex functions. We test quadratic forms with random matrices, aanooentialmminu--iiie plus functions with quadratic pairwise interactions. These experiments show that DLT maintains accuracy across strongly coupled functions, although training time increases due to the more complex computation tree required for gradient evaluation in non-separable functions.
functions, pretrained ICNNs, and coupled soft-
Section B.5 compares DLT with alternative algorithm that uses an approximate inverse gradient mapping instead of exact gradients for training. While DLT trains on exact targets using the Fenchel- Young identity with true gradient values. the alternative method must learn an approximate inverse mapping to construct its training targets. We show that although such an approximation is feasible, it fundamentally depends on the quality of the inverse gradient approximation and provides weaker theoretical guarantees.
Section B.6 describes DLT applications to Ramilton-Jacobl PDEs via the Hopf formula. We introduce DLT that approximates (9"+ +IH)* and compare against Deep Galerkin Methods, showing superior accuracy.
Section B.7 provides additional experiments with approximate inverse sampling for targeting specific dual space regions.
Section B.8 explores KANs integration with DLT for symbolic regression tasks, investigating whether approximate inverse sampling improves symbolic expression recovery in low-dimensional settings.
Grid-based methods: definition-bssed method and Linear Time Legendre Transform
The definition-based method
The definition-based method implements convex conjugation directly according to the Legendre-
Fenchel definition LI.
(*) = sup ((5,z) - f(x)}
=EX
In explicit multi-dimensional form: (81 02..,,......88.. 88) max 21.22
{8121 + 8202 + 8.0a - s(1 22, (2, x2))
Key characteristics of the definition based method
1. Brute-force approach:
The algorithm exhaustively evaluates all possible coordinate combi- nations for each slope combination.
2. Computational complexity: O(N24), where N the number of grid points per dimension and is the dimensionality. This complexity arises from:
(Nd)) operations iterate through all slope combinations
((N)) operations to search through all coordinate combinations for each slope
Memory usage: O(Nd) for storing the output array, with minimal additional temporary storage.
14

---

The definition based method provides a baseline implementation before considering more contempo- rary approaches like the Linear-time Legendre Transform (LLT).
The nested formulation
The nested method exploits the mathematical structure of the Legendre-Fenchel transform by reformulating sequence of nested maximization problems. For function / (() in d dimensions, the transform can be rewritten in nested form
(s1, $2, ())) = max $121 + max $222 + ... max {saTa - f(1 22 ................................ maa }
This decomposition allows dimensiion-b--dimessioo computttion,, starting from the innermost term.
Algorithm structure
ID Linear-time Legendre Transform (LLT) The LLT algorithm by Lucet [1997] efficiently computes (s) in O(n + m) time, where n is the number of data points and m is the number of slopes:
1. Convex hull step: Compute the convex hull of the discrete data points (z,, ff (x:)).
2. Slope computation: Calculate slopes between consecutive vertices of the convex hull:
G f( f(x)))
Fi+1 32
3. Merging step: For each target slope 8j, find the interval [Ca, G+tal containing 8,, then compute:
(*(,) = 83 Is f(x;)
d-dimensional algorithm
The nested approach proceeds in stages:
1. Compute the innermost term:
Va(æ1 Id- 8d) = max (sgad f311,
2d
2. Work outward for d- down to
Vi I; Si. ......,,s) = max (5,24 + Vs+1 (21,..., 3441...,,8))
2i
3. The final result f" (8, 82) = V3((8)...,58a)
Comparative characteristics
Assuming equal grid resolution in primal and dual spaces (n mn):
1. Dimensional decomposition: The nested method decomposes the d-dimensional problem into d sequential lower-dimensional maximizations. The resulting performance (generally) depends on the order in which this decomposition is performed.
2. Computational complexity: (ddv++1) for the nested approach, compared to 2() for the definition-based method. Each stage requires ((N++1) operations, with d total stages.
3. Memory usage: ((Nd) for the nested method.
4. Error dependency: Transform accuracy for both methods depends on grid resolution and therefore identical.
15

---

Results. Each benchmark table reports performance and accuracy metrics for various input dimen- sions d, using uniform grids with n m 10 points per dimension TT
The column lucet (s) reports the runtime of the Nested Lucet method in seconds, while lhet (s) gives the runtime of the brute-force definition-based method (evaluated only for d S 5 due to computational constraints). The columns MB₁ and p indicate the peak memory usage in megabytes (MB) during execution for the Lucet and Def methods, respectively; 0.0 reads that less than 0.05 MB was used. The column max L-I measures the maximum absolute difference between the Lucet and Def outputs across all grid points. Lastly, max err captures the interpolation error of the Lucet method relative to a reference solution obtained via uniform sampling, and RMSE reports the root mean squared error over the same validation sample. Some columns are indicated as since direct computations are too expensive to run.
Generally, both methods struggle with the curse of dimensionality; however, Lucet's nested LLT method is more efficient than the direct computation on the grids.
Table B1: Functions, grid domains (10' total points) and Legendre transforms.
Function
Expression
Legendre Transform
Quadratic
f(x) 0.5 Σ x
[-3, 3)d
f*) = 0.5 Ef
[-3,3] 3j
Neg-Log
[(x) Σ log(z;)
5.0j
(*) -d E log(-s:)
[-5.0, 001]
Quad-Linens
f(x)
Σ xd Eatt
3]
No explicit form
f[[0,
Table B2: Quadratic function benchmark: Lucet vs. definition- method for different dimensions.
d
&uori (s)
fDer (s)
MB,
max
RMSE
1
0.001
0,000
0.0
5.56e
4.11e-02
2
0.002
0.002
0.0
1.11e-01
8.01e-01 02
3
0.012
0.018
0.0
1.66e
1.17e-01
0.158
0.614
0.2
2.11c-01
1.50e-01
5
1,981
83.838
1.5
2.67e-01
1.86e
6
25.979
15.3
3.27e-01
2.20e-01
299.697
152.6
3.38e-01
2.66e-01
8
3419.780
1525.9
3.87c-01
2.90c-01
The nested approach demonstrates substantial computational advantages as dimensionality increases, albeit at the expense of increased memory requirements. However, despite the sparse grid used, it becomes intractable for 2 10.
Softmax-Baser approximation Legendre-Fencher transforms
The softmax-based approach, proposed by Peyré, see Peyré (2020) and Blondel and Roulet (2024), offers an alternative method for approximating Legendre- Fenchel transforms through convolution-based formulation. For a convex function f defined on bounded domain C' C IR the Legendre-Fennchel transform can be approximated using:
Eenechle-Fennhel
How we choose the dual sampling grid. Let G- = (2((0) be the Cartesian primal grid. For a proper,
convex, cocrcive function every optimal slope is sub -obo-rrddentt of at some sampled node, i.c. Of(g.) U. (z(())) contains all slopes that can ever maximise (s, x:) - f f(x). We therefore compute (analytically when possihle, otherwise once on the discrete grid) coordinate-wise bounds 8, S 0z, (x)( S for all a set the dual grid to the Cartesian product G. {s, + he MM-1 h. (3: s,) 1). If Of..) C G. the finite- slope-gre grid theorem d-D extension Lucet's Prop. Lucet 1997 see also Hiriart-Urrtty 1980 guarantees
fā_ (s)
max x² f(x) (*(s)
every E G..
so further enlarging the slope range cannot change the transform at any sampled point-only the interpolation error between nodes depends on how fine m is chosen. Example: for f(x) - - log 2; on x. E [0.1, 5) we have J E |-10,- -0.2]: taking uniform dual grid over [-10, -0.2) already covers all sub graddrnddnnsss and yields the exact conjugate values that grid.
16

---

Table B3: Benchmark results for the negative log function. Lucet and definition based methods
compared across dimensions. MB: active memory (MB); max L L-D|: D/: maximum discrepancy between methods.
d
Lnoe (s)
Ppe (s)
MB,
MB
max L D
max
RMSE
1
0.000
0,000
0.0
0.0
0,00c+000
3.95e-01
2.15c-01
2
0.001
0.001
0.0
0.0
8.88e 16
6.89e-01
3.25e 01
3
0.022
0.031
0.0
0.0
1.78c-10 15
9.38e-00
4.57e-01
0.154
0.659
0.2
0.5
3,55e
1.38c++00
6.36e 01
1.986
85.515
1.5
5.3
3.55e- 15
1.48++00
7.55e-01
6
25.952
15.3
1.99e+000
9.44e-01
_
7
305.878
152.6
2.19e+000
1009++00
8
3449.547
1525.9
2.26c-000
1.17c+000
Table B4: Estimated runtime comparison: definition bastd vs. Lucet methods with N = 10 points dimension. The definition- based method has ((N2d)) complexity while Lucet achieves Odd4+11).
d
Points/dim
Definition- method (est.)
Lucet method (est.)
Speedup
10
50.00
0.62
8.00c++01
6
10
1.39 hrs
7.50
6.67e++02
10
5.79 days
1.46 min
5.71e+03
1.59 yrs
16.67 min
5.00e++04
10
158.55 yrs
3.12 hrs
4.44c+005
10
10
15854,90 yrs
34.72 hrs
4,00e-006
11
10
1.59e-+06 yrs
381.94 hrs
3.64e+007
12
1.59c+-08 yrs
4166.67 hrs
3.33e+008
1.59c++ yrs
45138.89 hrs 3.08c+009
14
10
1.59c+ 12 yrs
4861 11.11 hrs 2.86c+10
(ʃo exp (z, (a,s)- f(x)
dx
where & >> 0 a smoothing parameter that controls the accuracy of the approximation. As 8 - 0, the softmax approximation f² converges to the exact transform f*.
This formulation can rewritten terms of Gaussian convolution:
f2(8) = Q=1
Q+(T) * Ge (s)
where: Gz(x) = exp(- 1 is Gaussian kernel, (z(f(()) exp - and (F)( - || log F(s).
Implementation approaches
Three primary techniques were employed:
Direct integration: For low dimensions (d 2), the integral is directly approximated grid.
FF'T-based convolution: For higher dimensions, the method uses the Fast Fourier Transform to compute the convolution efficiently, with complexity OddN' log N).
(Quasi) Monte Carlo: to approximate integrals. More memory efficient than integration but still expensive to iterate through dual grids or samples {
17

---

Empirical results
In our experiments Monte Carlo approach showed more stable results for performing integration for small values of E, so we focus it. Our MC implementation uses Sobol sequence generate low-discrepancy sample points, achieving better coverage of the integration domain compared to purely random sampling. Table presents results the negative logarithm function using 100,000 samples. The dual grid is the same as in LLT experiments from the previous section, see Table B3
Experiments confirm (E) convergence to the exact Legendre transform - Peyré's approach is suitable when an estimate of the transform is required at only one point in the dual space, but it is still too expensive iterate over dual grids/samples, especially nigh-dimensional setups: d 4.
Table B5: Results Peyré for negative log function using quasi- Carlo (Sobol) with 100,000
samples across dimensions and smoothing parameter
Sample Size
Time (s)
MB
Err
RMSE
0.001
4MC(100k)
1.690
6.43
1.98e-00
4.43e 02
0.010
qMC(100k)
1.714
6.13
2.17e-01
5.33e
0.100
qMC(100k)
1.711
6.14
2.21e-01
1.30e-01
0,500
qMC(100K)
1.793
6.13
6.38e-01
3.55e-01
0.001
qMC(100k)
3.833
9.26
3.95c-01
6.54c-02
0,010
qMC(100K)
3,936
9.26
4.33e-01
8.75e-02
0.100
(MC(100k)
3.361
9.26
4.42e
2.39e-01
0.500
qMC(100K)
3.365
9.29
1.28
5.20e-01
0.001
(MC(100K)
23.794
13.16
6.03e-01
8.73e 02
0.010
qMC(100K)
29.116
12.46
6.49e 01
1.21e-01
0.100
qMC(100K)
16.801
12,46
6.64e-01
3.47e-01
0,500
qMC(100K)
16.960
12.46
1.9e+000
6.56e-01
0.001
9MC(100k)
419.791
18.21
8.29e-01
1.60e
0,010
9MC(100K)
561.630
18.66
8.65e-01
1.93e-01
0.100
(MC()00K)
317.030
20.06
8.85e-01
4.53e-01
0.500
(MC()00K)
317.895
20.06
2.55e+00
7.80e-01
Table B6: Results of Peyré for negative entropy function using quasi-Monte Carlo (Sobol) with 100,000 samples across dimensions and smoothing parameter
d
Sample Size
Time (s)
MB
Max Er
RMSE
0.001
qMC(100K)
4.379
11.61
4.42e-01
1.03e-00 001
0.010
qMC(100k)
4.196
6.15
4.73e-01
1.18c-01
0.100
qMC(100k)
4.145
6.13
5.83e-01
2.19c-01
0.500
qMC(100k)
4.150
6.14
9.02c-01
4.75c-01
0.001
qMC(100k)
8.207
9.26
8.87c-01
1,57c-01
0.010
qMC(100K)
8.206
9.28
9,47e-01
1.87c-01
0.100
qMC(100k)
7,933
9.26
1.1e++00
3.75e-01
0.500
qMC(100k)
7.993
9.26
1.80c+00
7.04c 01
0.001
qMC(100K)
30,088
12.47
1.39++00
2.29e-01
0.010
9MC(100k)
33.914
12.50
1.44e+000
2.60c-01
0.100
qMC(100k)
22.726
13.92
1.75e+00
5.26e-01
0.500
qMC(100k)
22.507
13.22
2.70c+00
8.98c-01
0.001
qMC(100K)
403.168
18.61
2.01c+00
4.10c-01
0.010
qMC(100k)
524.356
17.96
2.00c+00
4.33e-01
0.100
qMC(100k)
301.965
20.07
2.51c++00
7.17c-01
0.500
(MC(100k)
306.520
17.91
3.57e+400
1.07e+000
Our experiments showed that the method performs better for smooth functions such as negative entropy. and less effectively for functions with steep gradients like the negative logarithm.
We tested several MC variants (basic MC/ Sobol- MMobllQMM//// importance adaptive); Sobol- Sobol-QMC was most
stable in this setting.
18

---

DLT: Comparing different NN architectures
Comparison with direct learning when the Legendre transform is known
critical validation of DLT is to demonstrate that it achieves state-of-the-art performance by comparing it to the direct learning of the convex conjugate in cases where the analytic form of f is already known. For this comparison, we implemented a comprehensive benchmark that compares DLT against explicit direct supervised learning across several convex functions with known convex conjugates in multiple dimensions.
Benchmark setup
We compare learning approaches:
DLT: trains neural network go to satisfy 9o (v(((()) = (z, Vf(x)) - f(x).
Direct: directly trains ge so that ge = (*),, which is only possible if f is known explicitly.
The direct approach serves as our "gold standard" since it has direct access to the target function values. However, this advantage is only available in the special cases where f has a known analytical expression- assumption that does not hold in many practical applications. Our goal is to demonstrate that DLT matches the performance of this idealized case.
We implemented both approaches using identical neural network architectures to ensure a fair comparison:
MLP:
Standard layer perceptron with two hidden layers (128 units each) using ReLU activations
Feed- ferd-forward network with enforced convexity using only positive weights (two hidden layers with 128 units each)
MLP-ICNN:
ResNet: Residual network with two residual blocks (128 units each)
ICNN: Input-convex neural network of Amos al. [017 with positive z-weights and direct x-skip connections (two hidden layers with 128 units each)
The implementation uses JAX for efficient computation with identical optimization settings for both approaches: Adam optimizer with learning rate 10-3 and batch size 128 d. We employed an early stopping criterion where training terminates if the Lă error falls below 10-6, with a maximum limit of 50,000 iterations ensure computational feasibility across all experiments. To estimate statistical significance, we use 10 runs and report both the mean p and standard deviation a of residuals achieved with direct and indirect methods.
key consideration in our benchmark design was ensuring that sampling distributions are correctly matched between primal and dual spaces. The gradient map Vʃ C transforms the distribution of points x E C' into a distribution over = Vf(x) E D. For each function, we designed specific sampling strategies:
Quadratic f(x) |=| Since Vf(x) x, we sample x ~ ~( T) for the implicit approach, which naturally gives = x ~ N00 I) in the dual space for the explicit approach.
Negative logarithm f(x) - -'_ log(z,): For the implicit approach, we sample uniformly in log-space from 10) df to adequately cover the positive domain. Since Vz(z) = this transforms to E (-10, -0.1d in dual space. For the explicit approach, we directly sample y from this same distribution.
Negative entropy f(x) Ei_ 2, log(x.): Similarly, sample x uniformly in log-space for the implicit approach. With gradient fz(z) log(z) + 1, this maps to (-1.3, 3.3)d in dual space, which we match in our explicit approach sampling.
PDLT is implicit approach, effectively, defined by L.5
19

---

This distribution matching is essential for a fair comparison, as it ensures both approaches face comparable learning challenges.
Results and analysis We tested these three convex functions with known analytic conjugates across dimensions d E {2, 5, 10, 20}:
Quadratic: f(x) |01
1
(")
2
2
d
Negative logarithm: f(x) Σ
log(zs);
(² - Σ
log(- (og(-Wi) -d
i=1
==1
d
Negative entropy: f(x) Σ J; log(2;),
(" Σ
exp(y): 1)
The quadratic function f(z) = | is particularly elegant for our comparison since it is self- conjugate (f f) and has gradient V f(r I. This means that when using DLT with samples (O, 7), we get y = V (x)) = x ~ (0, T) in the dual space well, making the comparison especially clean.
Table B7 presents L errors and training times for both approaches. These benchmark results reveal key finding: DLT matches the performance of direct learning if the convex conjugate known analytically. The 2 errors for both methods are statistically identical.
When the dual validation samples drawn from the training distribution, DLT learns the convex conjugate just as accurately (and almost as fast) as direct supervised learning of f - even though never sees the true dual values.
These results validate a key advantage of our approach: DLT achieves the same approximation accuracy as methods that require knowing the analytical form of f", without actually needing this knowledge. This makes DLT applicable the wide range of problems where f does not have known closed form, while still maintaining state-of-the-art performance in cases where the Legendre transform is known.
Experiments with non-separable functions
To demonstrate DLT's capability beyond simple separable functions, we evaluate its performance on several challenging non-separable convex functions. Table B10 presents results across four representative function classes:
Quadratic form with random SPD matrix: We test f(x) = f=='Qz where Q is a randomly generated symmetric positive definite matrix with controlled condition number x R 102. The coupling through Q makes this strongly non-separable, requiring the approximator learn the full quadratic form structure.
Exponeentillminnu--ineee::: The function f(x) exp((a,x)) - (b, x) where a, b E g are independently sampled unit vectors. Specifically, we sample Craws braw ~ NO, Ia) and normalize them as Graw/ /lwrl//azwww |l2 and b braw/ (bba// 12-
Pretrained ICNN: We first train 2-layer ICNN (128 hidden units per layer) for 50,000 iterations to approximate target quadratic function, then freeze it and use it as the convex function f. This represents a realistic scenario where the target function itself a learned model from data.
Coupled Soft-plus: The function f(z) Eie, log(1 + exp(x; + +z)) explicitly couples all pairs of variables with (dd) interaction terms. This severely challenges classical methods due to the combinatorial explosion of terms.
Experimental setup. All experiments employ ResNet approximators with 2 blocks of 128 hidden units each, trained using the DLT framework. We train for I 120,000 iterations with learning rate of 10-3 and Huber loss (6 1.0) to ensure robust convergence. For each function class, we apply domain-specific whitening transformations. Training samples are drawn from appropriate domains
20

---

Table B7: Comparison of DLT with direct learning for known Legendre transforms
L Error
Time (s)
Function
d
Model
DLT
Dir
DLT
Dir
MLP
2.96c-05
2.86e-05
85.5
84.4
MLP-ICNN
6,21e+00
3.15c+01
87.4
76.9
10
ResNet
2.43e-05
3.00e-05
106.8
101.2
ICNN
1.79e-02
1.23e-02
92.7
93.0
MLP
2.28c-04
3.41c-04
84.9
84.0
Quadratic
MLP-ICNN
1.20c+02
1.13e+02
71.3
82.6
20
ResNet
2.53e-04
2.32e-04
107.7
104.9
ICNN
3.45e-02
3.48e-02
95.1
93.3
MLP
3.55e-03
3.59e-03
87.5
86.6
MLP-ICNN
6,46e+02
6.46e+02
87.3
76.6
50
ResNet
2.92e-03
3.25e-03
123.8
119.7
ICNN
7.18e-02
6.66e-02
97.4
97.8
MLP
6.75e-04
8.43e-04
94.1
94.4
MLP-ICCN
4.23e+00
4.23e+00
63.7
62.4
10
ResNet
7.34e-04
9.22e-04
108.7
111.0
ICNN
6,47e+00
4.22e+00
54.1
65.4
MLP
4.31e-03
3.60e-03
93.8
95.8
Neg. Log
MLP-ICNN
7.43e+00
7.16c+00
86.2
96.8
20
ResNet
5.36e-03
6.20e-03
114.0
115.2
ICNN
8.75e+00
8.69e+00
71.4
63.5
MLP 3.81e-01
2.81e-01
97.5
97.8
MLP-ICNN
2.03e+01
2.03e+01
67.3
70.4
50
ResNet
1.05e-01
1.23e-01
133.5
130.7
ICNN
2.03e+01
2.06e+01
73.5
61.9
MIP 1.06e-03
1.26e-03
95.6
100.3
MLP-ICNN
4.56e+00
4.01e+00
92.8
101.0
10
ResNet
6.62e-04
7.41e-04
109.8
115.7
ICNN
1.31e-02
1.74e-02
98.8
107.2
MLP
5.88e-03
7.83e-03
94.5
99.8
Neg. Entropy
MLP-ICNN
2.80e+01
2.82e+01
92.7
102.2
20
ResNet
3.88e-03
4.24e-03
113.7
118.6
ICNN
3.68e+01
3.67e+01
68.1
72.3
MLP 2.71e-01
3.91e-01
99.6
100.8
MLP-ICNN
7.83e401
7.87e+01
97.5
104.6
50
ResNet
6.93e-02
7.27e-02
132.8
135.0
ICNN
9.39c401
9,39c+01
59.8
69.1
with truncated normal or uniform distribution (see Table B9), with batch sizes adaptively set to 25%
of the dimension (B = [0.25d We evaluate performance at dimensions E (10, 20, 50) and report results averaged over independent trials with different random seeds.
Results. The results demonstrate that DLT successfully learns all four non-separable function classes, achieving low RMSE emrors across different dimensions. The quadratic SPD function shows the best performance with errors decreasing as dimension increases (from 1.15 x 10 at d = 10 to 2.01 x 10-3 33 at d = 50), demonstrating DLT's effectiveness learning structured covariance patterns. function presents greater challenges due to its unbounded gradient
The poponential-ninu--inenn domain, yet DLT maintains reasonable accuracy. The pretrained ICNN case validates DLT's ability to approximate learned neural network functions.
21

---

Table B8: Comparison of DLT with direct learning for known Legendre transforms
L Error
Time (s)
Function
d
Model
DLT
Dir
DLT
Dir
MLP
3.55e-03
3.59e-03
87.5
86.6
MLP-ICNN
6.46e+02
6.46e+02
87.3
76.6
50
ResNet
2.92e-03
3.25e-03
123.8
119.7
ICNN
7.18e-02
6.66e-02
97.4
97.8
MLP
2.77c-02
2.46c-02
95.0
95.3
Quadratic
MLP-ICNN
2.54c+03
2.59c+03
96.8
87.4
100
ResNet
1.88c-02
1.91c-02
136.3
132.0
ICNN
9.59e-02
9.06e-02
105.7
105.5
MLP
1.43e-01
1.50e-01
113.5
112.7
MLP-ICNN
1.01e+04
1.01e+04
97.6
113.7
200
ResNet
9.53e-02
9.40e-02
163.1
164.4
ICNN
4.01e-01
3.73e-01
124.1
125.8
MLP
3.81c-01
2.81e-01
97.5
97.8
MLP-ICCN
2.03e+01
2.03e+01
67.3
70.4
50
ResNet
1.05e-01
1.23e-01
133.5
130.7
ICNN
2.03e+01
2.06e+01
73.5
61.9
MLP
8.56e+00
6.67c+00
106.2
106.1
Neg. Log
MLP-ICNN
4.00e+01
3,99c+01
73.0
66.2
100
ResNet
1.38e+00
1,93e+00
141.5
142.3
ICNN
3.99e401
3.98e+01
80.3
71.1
MLP
3.29e+01
2.68e+01
124.7
125.1
MLP-ICNN
8.21e+01
8.20e+01
68.9
78.7
200
ResNet
2.92e+01
2.99e+01
168.7
174.1
ICNN
1.70e+01
1.68e+01
131.2
137.0
MLP 2.71e-01
3.91e-01
99.6
100.8
MLP-ICNN
7.83e+01
7.87e+01
97.5
104.6
50
ResNet
6.93e-02
7.27e-02
132.8
135.0
ICNN
9.39e+01
9.39e+01
59.8
69.1
MLP
1.77e+01
1.65e+01
106.5
108.7
Neg. Entropy
MLP-ICNN
1.45++02
1.44++02
103.5
111.2
100
ResNet
1.53e+00
1.32e+00
142.0
144.2
ICNN
1.93e402
1.93e+02
62.9
79.9
MLP
4.92e+01
5.07e+01
124.2
127.3
MLP-ICNN
2.87e+02
2.87e+02
125.3
133.3
200
ResNet
5.08e+01
4.02e+01
94.2
175.4
ICNN
9.43e+01
9,30c+01
131.9
139.2
Notably, the coupled soft-plus function requires approximately 2x longer training time at d 20 and = 50 due to its (dd) pairwise interaction terms and longer gradient evaluation step. The relative RMSE remains below 2% for most cases, confirming that DLT provides accurate approximations suitable for downstream optimization tasks. These experiments establish that DLT with simple ResNet architectures can effectively handle diverse non-separable convex functions that would be intractable for classical decomposition methods.
Alternative NN approaches to computing Legendre transforms
DLT approximates the convex conjugate I* on D VS(C) by leveraging the exact Benchel-Young identity gradient-mapped points. The key insight is that for any x E with y = Vf(x), the identity (*) (I, 2) f(x) holds exactly.
22

---

Table B9: Test functions with sampling domains. C indicates the effective training region (primal domain), D = Vf(C') is the corresponding dual domain for gradient samples.
Function
Expression
C (Primal)
(Dual)
Quadratic SPD
f(z) = 32 T Qx
[-3,3]
--31Q112: 31Q
(random SPD, ||Q||
(std. normal 133)
(<2 [-333] when normalized)
Exp-minus-I.in
f(z) = e (a.2)) - (h, x)
3-3
[e [e-3vā - 1,, 1]a
(|a/|2 | = 1)
(std. normal +33)
ICNN (2-layer)
f(z) ICNN128, (2)
[-3,3]
Learnd, data-depenndet
(pretrained to quadratic)
(<s |3,3]*
Soft-plus pairs
f(z) = Eic, log(1+e"itt)) (pairwise coupling)
[-1,5 1.5)
~ d0
(uniform
Table B10: Performance of DLT with ResNet approximators for separable functions. Absolute RMSE, relative RMSE (normalized by function scale) and training times are reported as mean = standard deviation over 5 independent trials.
Function
RMSE
Relative RMSE
Time (s)
10
1.15e-2
3.7e-3
1.05e-2 + 3.9e-3
584 + 9
Quadratic SPD
20
5.13e-3 + 1.1e-3
3,98e-3 +: 6.3e-4
1016 + 10
50
2.01e-3 # 2.0e-3
1.05e-3 1.1e-3
2709 2
9.18e-2
3.66--
2.50e-2 1.0e-2
599 9
Exp-minus-Lin
20
5.78e-2 + 3.3e-2
1.60e-2 + 8.8e-3
1035 + 14
50
4.70--
1.41e-2 1.2e-2
2712 0.2
10
2.37e-2 + 7.1e-3
6.17e-3 1.4e-3
706 7
ICNN (2-layer)
20
1.92e
5.3e-3
6.06e-3 + 1.6e-3
1200 + 10
50
1.13e-2
+
3.69e-3 8.1e-4
2759 + 1
10
7.9e-3
5.21e-3 + 3.6e-3
965 + 9
Soft-plus pairs
20
2.63e-2 9.4e-3
4.14e-3 1.5e-3
2181 11
50
1.25e-1 + 4.3e-2
4.95e-3 + 1.7e-3
5498 + 13
We train neural network go D + R by minimizing the loss:
Iga V (=)) + f(x) (z, V(x)))
(DLT)
where x C is sampled to ensure = V f=)) adequately covers D. Crucially, this loss equals (g)(v) - (w)) for points on the manifold y = v (x)), providing an immediate accuracy certificate through the training residual itself.
In contrast,
consider a proxy method by using approximate target:
(*) (ho(v). f(hayy)), where ha (Vf)-
(Proxy)
This introduces potential bias since ha only an approximation of the true inverse gradient. More precisely, if a model is trained to the proxy target to (1) = (he y) f(())), then
go(v) (1) gv (v) (v) (v)
optim.
bias from
Even if optimization drives the first term down, the certificate still reflects the second term, which is zero only when ho (V 1 on D. While proxy methods can achieve low training loss on their approximate targets, evaluate all methods using DLT certificate (DLT), which measures the true approximation error regardless of the training approach used.
DLT provides several critical advantages over the Proxy method:
Guaranteed convexity: When using an ICNN architecture, DLT guarantees the convexity of the output ge since the training target is the true convex conjugate f*. The Proxy method
23

---

loses this guarantee because even with an ICNN, the approximator learns a biased target (ho(v). y) f(ho (v)) where he may not be the true inverse gradient.
2. No accuracy bottleneck: DLT's accuracy is limited only by the neural network's approxi- mation capacity, not by an intermediate operator approximation. The Proxy method's final accuracy is fundamentally bottlenecked by the precision of he z (Vf) regardless of the main network's 's capacity.
3. Efficient sampling: DLT uses the exact gradient mapping V for sampling, which is computationally efficient and exact. The Proxy method requires every new sample to pass through the learned inverse network h2, adding computational overhead and potential approximation errors at inference time.
4. Training overhead: While both methods may use inverse mapping for sampling control, the Proxy method critically depends on the quality of this mapping for its training targets, requiring more careful and time-consuming pretraining.
Table B11: DLT with approximate-inverss sampling vs. proxy. ha RMSE measures prelearned inverse quality; tx inverse pretraining time.
Function
d
Method
Train (s)
RMSE
ha RMSE
f/1 (8)
5
DLT iav.sampling
113.95
8.29e 03
8.23e-00
287.8
Quadratic
Proxy
104.60
1,49c-02
287.8
10
DLT inv.sampling
222.05
2.02e 02
7.4le 05
468.4
Proxy
209.22
2.10p-02
468.4
5
DLT inv.sampling
107.96
3.98c 02
5.79e 06
229.7
Neg. Log
Proxy
104.28
4.21e 02
229.7
10
DLT inv.sampling
209.41
1.28c-01
8.03e 05
390.2
Proxy
204.96
1.45e-01 01
390.2
5
DLT inv.sampling
110.70
3.27c-02
3.27e-05
326.5
Neg. Entropy
Proxy
106.79
3.25c-02
326.5
10 DLT + inv.sampling
221.21
3,53e 02
1,88e 04
475.0
Proxy
212.85
3.64c-02
475.0
Table B11 reports numerical resulty for the functions from Table B4 comparing DLT with approximate-inveee sampling with the inverse gradient proxy baseline for dimensions d' E (5, 10). Both methods share the same inverse pretraining time t, which ranges from approximately 4 to 8 minutes depending on the function and dimension. We learn the inverse mapping he to high precision (he RMSE between 10-6 and 10 10-4). As result, the certificate RMSE-which measures the true approximation error of becomes comparable in terms of training times between the two methods.
Handling out-of -atmaf- inverse mappings. When using approximate inverse mapping he D > Rdd for sampling, some points hy(y) may fall outside the constraint set C. Two approaches exist for handling these cases: (i) adding penalty term the loss function that penalizes violations of h2) & C, or (ii) discarding such samples from the training batch. While the penalty approach could
"Both methods employ ResNet architectures with blocks of 128 hidden units and GELJ activations. The key is two-stage procedure: first, we pretrain an approximate inverse mapping hump using an autoencoder-style cycle that maps - z g(y) + & = l + V (2)), ensuring & remains within the primal domain C through projection. This pretraining phase uses 200,000 steps for d = 5 and 400,000 steps for d = 20 with larger batch sizes (1.01x x the streaming batch). In the main training phase, both methods use streaming mini-batcher with fresh samples y Uniform( ()) drawn at each step, where the dual domain D carefully matched to the primal domain C' via the gradient mapping. To ensure numerical stability, we apply interior clipping with relative tolerance = 10- to keep gradients strictly within valid domains. The streaming batch size N is set explicitly (600 for d = 5. 1280 for d = 20) and scaled by 0.99x for both methods, with training steps compensated proportionally to maintain constant total sample complexity (30,000 base steps for d 5, 100,000 for = 20). Training employs AdamW optimization with learning rate 10-3 weight decay 10-6 and exponential decay schedule (halving every 20,000 steps). Out-of bounds gradient samples are dropped. Performance is evaluated on test set of 5,000 (z, 3) pairs measuring both mean squared error and relative La error || - f [[] the conjugate function approximation.
24

---

theoretically guide the network toward valid mappings, it introduces additional hyperparameters (penalty weight, barrier function choice) and may distort the training dynamics by mixing two abjectives-approiiaati the inverse and satisfying constraints. We adopt the discarding approach as it maintains the purity the training objective.
Hamilton-Jacohi equation and time-parameterized inverse gradient sampling
The Hamilton-Jacobi (HJ) equation a comeestone of classical mechanics, optimal control theory, and mathematical physies. It is given by the first-order nonlinear partial differential equation
0,u + ((V (V, ((x t)) 0,
where u(z, t) represents the action (value) function, and H denotes the Hamiltonian associated with a physical system or control problem.
Consider the Cauchy problem:
0eu + HVV.u(z,)) = =0, (x, 1) E Q, u(x, 0) = g(z),
E Rd
where Q Ed (0, o0).
The Hopf formula and function definitions
The Hopf formula provides an explicit solution to the Hamilton-Jacobi equation:
u(x, t) (g" (x) sup ((r (H)
(BI)
PERd
For our computational experiments, we consider classes of functions that exhibit different mathematical properties and computational challenges.
Quadratic functions.
These represent the classical case with explicit analytical solutions:
Initial condition: g(x) 1 Ei=
2
Conjugate: gp( -i=1 13
Hamiltonian: H(p) > p2
Σi # 2(1+0)
In this case, the analytical solution is u(x,t)
Exponential functions.
These provide more challenging case:
Initial condition: g(x) Ed t
Hamiltonian: H(p) Σd= ep
Conjugate: gp) Σi (p: log p:) (when p; > 0)
To approximate the HJ solution u(x,t) t) time-parameterized network uo(x, with input (x. t).
Time-parameterized inverse gradient sampling
Denote v(p,t) t) (p) + (p). Then the Hopf representation can be written
u(x, t) = o(- (x) sup ((x,7) p) - {)).
pER
Whenever u(; t) is differentiable (and, e.g., strictly convex so the maximizer is unique), the maximizer p p(x, 1) is characterized by the first-order optimality condition
t) = (p*(p) +tVH(p),
hence
p Vou(x, t) (Vpu(,t))
25

---

Time-DLT therefore requires samples of pairs (z, p) satisfying this relation. In the case where H and g are quadratic one has (p,, t) (1 t)p, so (1 - t)p (equivalently p 2/(1 + t)). which allows direct sampling. In non-quadratic cases the inverse map (Vpo(,
is not available in closed form: to generate samples with a prescribed distribution in z-space we learn a time- parameterized grrniiet network (o( (V,v(,t) (xx). Similar to Section 5 we train 1), (5.1)-Pli -ll2
loss funcitons, which only require evaluating Vp Vg* + IVH. In the implementation, this provides the inverse-gradient sampling step p R V (r, 1) (and in the quadratic case the explicit formula p 2/(1 + ()) used to form training triplets (z, p, 1).
Vo using N V. (2 (), 1) 2:12 and
Time-DLT implementation
Time-DLT directly approximates the Hopf formula through with neural network uo( 1) that learns:
uo(x,() u(z, () sup ((r,p) (H))
The training process minimizes the Legendre transform residual:
= E (m.p..)
[[uo(z t) - ((x, p) g" - (H()))
The sampling of triplets (x, p, t) depends on the function type: analytical formulas for quadratic functions and the trained inverse network for exponential functions. This approach preserves the structure of the Hopf formula while ensuring computational tractability.
Experimental framework
We evaluate both methods across increasing dimensions and function complexities. We sample (z,, t) Unif [Tmin Tmax]" dd x [0, tmau] with Tmin -2, Tmar 2, tmaz and report
in Tables B12 and BB the resulted RMSE and PDE residuals across different time moments t {0.1, 0.5, 1.0, 2.0).
Problem configurations. We evaluate both methods across the following configurations:
Dimensions: de (2,5, 5, 10) to assess scalability
Function combinations:
Quadratic-Quadratii:: baseline with known analytical solution
Quadratic-Expponnttaa Hamiltonian H is exponential
Exponential-Quadrati:: initial condition is exponential
Neural network architecture Both methods employ pamametramttri multilayer perceptrons (MLPs):
Input:
(z,t) E Rd R
Architecture: (d + 64 64
Activation: tanh
Optimizer: Adam with learning rate r -0-3
For Time-DLT, the inverse gradient network hS uses:
Input: (p,t) E pd x +
Architecture: (d + 1) 64 64 +
Two-phase training: 40% of epochs for hgg 60% for Fa
26

---

Implementation details
Ground truth computation. To assess accuracy when an explicit solution is unavailable, we compute reference solutions by evaluating for each point (æ:,) reference value by solving the inner optimization over in the Hopf representation (B1
Low dimensions (d S 2): Scipy's L-BFGS-B algorithm
Higher dimensions d < 5): PyTorch gradient-based optimization with Adam
Training specifications.
Batch size: 100 samples per iteration
Training epochs: 1500 for d 2, scaling 3000 for d 10
Sampling strategy: Adaptive mixing of initial condition and interior points
Regularization: 100 for initial condition enforcement
Gradient clipping: Maximum norm of 10.0 for stability
Evaluation metries
We assess performance using three complementary metrics:
1. 2 error against ground truth:
$
ua(xi () uer((xs f)
2. PDE residual:
EPDE V (2,1) 11::20 + H(Vzuo))
3. Initial condition error:
E / Ex [110(2, -)x()))
The results demonstrate several key findings:
1. Accuracy: Time-DLT consistently achieves significantly lower errors compared to DGM across all problem configurations.
2. Scaling with dimension: Both methods show increased errors dimension increases, but Time-DLT maintains superior performance even in higher dimensions.
The Time-DLT approach demonstrates that by explicitly modeling the time- panamarraattee convex conjugate (g" + f/1)", we can achieve more accurate and efficient solutions compared to direct PDE approximation methods. This approach is particularly effective for high-dimensional problems where g is known in closed form. If DLT is used to compute 9 as well, errors from such a double Legendre transform might accumulate, potentially degrading overall approximation accuracy.
27

---

Table B12: Comparison of DGM and Time-DLT on the Hamilton-Jacobt equation with quadratic initial condition and Hamiltonian
d
DGM
Time-DLT
Err. Ratio
2 Error
PDE Res.
Time (s)
L Error
PDE Res.
Time (s)
Mean
2
0.1
8.88e-3
9.63e-2
80.66
3.43e-3
6.81e-1
56.29
2.69
0.46
0.5
2.28e-2
4.98e-2
3,89e-3
3.67e-1
7.53
1.64
1.0
3.92e-2
3.36e-2
3.23e-3
2.07e-1
11.38
2.43
2.0
4.20e-2
1.93e-2
2.65e-3
9.22e-2
17.39
6.11
5
0.1
3.34e-2
2.47e-1
82.75
2.81e-2
1.10c+0
57.42
1.24
0.13
0.5
1.09c-1
1.27e-1
1.98e-2
5.87e-1
5.55
0.85
1.0
1.55e-1
8.27e-2
1.39c-2
3.31e-1
10.29
1.77
2.0
1.16e-1
3.97e-2
1.08e-2
1.49e-1
13.92
4.25
10
0.1
8.90e-2
4.82e-1
82.40
8,46e-2
1.60e+0
58.76
0.95
0.09
0.5
2.14e-1
3.04e-1
5.25e-2
9.17e-1
3.69
0.44
1.0
2.26e-1
1.85e-
3.68e-2
5.12e-1
7.32
1.30
2.0
1.64e-1
8.05e-2
2.72e-2
2.27e-1
8.64
2.57
20
0.1
2.53e-1
9.51e-1
87.39
2.43e-1
2.25e+0
88.85
1.03
0.09
0.5
3.78e-1
6.83e-1
1,49e-1
1.36e+0
2.96
0.34
1.0
3.87e-1
4.03e-1
1.13e-1
7.50e-1
4.22
0.97
2.0
4.39c-1
1.81e-1
7.36e-2
3.24c-1
6.14
0.96
30
0.1
3.45e-1
1.18++0
87.18
4.29e-1
2.69e+0
192.62
0.83
0.12
0.5
6.40e-1
8.85e-1
2.60e-1
1.82e+0
2.63
0.34
1.0
1.01e+0
4.68e-1
1.90e-1
9,69e-1
5.15
1.23
2.0
1.21e+0
2.56e-1
1.28e-1
3.98e-1
9.00
2.21
Notes: L Error represents the mean squared difference between predicted solutions and the analytical Σ x2
solution u(x, t) = 2(+1))
HJ Error measures residual of the Hamilton-Jacobi equation |0,u
H H(V.u)).
Table B13: Comparison of DGM and Time-DLT on the Hamilton-Jacoke equation with quadratic
initial condition and exponential Hamiltonian
DGM
Time-DLT
Error Ratio
L Error
PDE Residual
Time
L Error
PDE Residual
Time (s)
DGM / DLT
1.46c-01
2.25e-01
142.40
6.88e-02
2.24c-01
95.35
2.12
0.5
1.30e+00
1.27e-01
1.29c400
3.29e-02
1.01
2
1.0
3. 3,25c+00
8.21e-02
3.21c+00
7.19e-02
1.01
2.0
7.70c+00
6.48c-02
7.69c+00
6.87c-02
1.00
0.1
3.18e-01
5.37e-01
215.76
1.14++00
1.43e+00
145.49
0.28
0.5
5.53e-01
4.16e-01
4.94e-01
9.03e-01
1.12
1.32e+00
3.67e-01
4.69e-01
4.55e-01
2.80
2.0
2.64e+00
6.53e-01
4.90e-01
2.70e-01
5.39
Notes: 2 Error quantifies the mean squared deviation from the ground truth solution. PDE Residual the residual of the Hamilton-Jacobi equation 00 H.))
28

---

Commenivon at tsf.1
d+ H-exprr ental. IC-quadrakk
44
1.-
0.35
10-
0.19
8
20
1.
40
C6
C0
0.52
-10-
-3.
-30
0138
31
2.0
-2.
-20
0.03
Lie0.5
d:2 Hieeppmential, IQ-quadrat
Proper
TIme-OLT
2.
20
20
30
50
10
10
10
024
920
t0
10
60
036
08
-1.2
-16
-10
-25
000
-12
-2.
-20
-20
0.00
H-exxponenuaall..mmmm IC-quedrat
DGM Solution
Praper Time Solution
[DGM DUTI
2.
2.0
36
2.0
0.72
14
18
40
10
040
57
60
035
-16
00
10-
16-
50
0.0
-20
Figure B1: Comparison of DGM and Time-DLT at three different times = 0.11, 0.11,.5, for d = 2.
Approximate inverse sampling: additional experiments
In this section, consider two examples of approximate inverse sampling in the context of the Deep Legendre Transform.
The first example, when approximate inverse sampling might be particularly useful, is when the target distribution / is concentrated around certain points { E D of particular interest. Then, on the one hand, one has sample from distribution μ4, such that V μ has high density the points {y E D. On the other hand, if the points are distant enough, it might be inefficient to try to learn f* on the whole D vf(C) given that are interested in higher precision in the neighborhood of {yi ) E D. We illustrate that approximate inverse sampling can solve both problems, efficiently sampling (approximately) to the target distribution v.
all
In the second example, we explore the ability of KANs to find symbolic representations low- dimensional setups with and without approximate inverse sampling. The setup is similar to that of Section 6. of the main part. Surprisingly, KANs managed to find correct symbolic expression even with quite skewed sample.
29

---

The case of target distribution concentrated around specific points
We consider the case of the multidimensional negative entropy
d
f(x) Σ x; log(z;)
defined on the positive orthant C == {a E Rd xq > 0 for all o]. It is straightforward to check that
- -i= exp(x; 1)..
Consider two points in IR point with coordinates y: 1, for all i, and point with all coordinates equal to 4. The target distribution v modeled as fair mixture of normal distributions Ny,, 10-3 Ia) and (2, 10 10)3
Figures B3 and B3
Illustrate the process of learning the inverse gradient mapping in this setup.
Ipoch
Irverse sempling
50
60
120
Figure B2: Distributions μ and g (v)
Figure B3: Learning mixture of Gaussian distribution with the inverse sampling: leftmost train step 1, rightmost train step 3300.
Kolmogorov-Arnold neural networks and approximate inverse sampling
Inspired by the Kolmogorov- representation theorem, which shows that every continuous function 6: j0, 1 - IR can be written as
2d+1
Σ 2 Σ Een((es)
((x)
for one-dimensional functions d.j
[0, 1] + R and
0; IR + H. KAN is composition
4, Rdd IR.
of KAN-layers of the form
(x(x))
x. E Rdd-.
di 1 @ ( (Irr) -=
30

---

where the one-dimensional functions d/yi R + R are specified as
ΣaB c.
for given basis functions B; and learnable parameters G₁ IR.
The canonical choice of B; are polynomial splines with learnable grid-points. But, depending on the situation, some of the B can also be specified as special functions such as exp, log, sin, cos, tan.
31

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The abstract and introduction clearly and accurately state the main contribu- tions and scope of the paper, including the derivation of a posteriori error estimates, which are supported by the results in Sections 5 and the main theoretical result in Section 3.
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: The paper explicitly discusses limitations, such as the focus solely on differen- tiable convex functions, and does not address non-convex or non-diffcrentiable scenarios.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full of assumptions and complete (and correct) proof?
Answer: [Yes]
The paper clearly states the assumptions required for the main theoretical result in Section 3. including differentiability and convexity. The proof is clearly presented and fully justified by invoking standard results such as law large numbers.
Justification:
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Answer: [Yes]
Justification: Sections 5.1, 5.2, and the supplementary materials fully disclose model architectures, data distributions, training details, and evaluation metrics needed reproduce the experimental results.
5. Open access to data and code
Question: Does the paper provide open access the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: The paper explicitly provides access to the code, data, and sufficient instruc- tions the supplementary materials to reproduce all experiments.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Experimental settings, including network architectures, activation functions, optimizers, data splits, and sampling methods, are specified in detail in Section and supplementary materials.
Justification:
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: The paper includes appropriate error bars and statistical significance analysis for the main numerical results, detailed in the supplementary materials.
32

---

8. Experiments compute resources
Question: For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Answer: [Yes]
Justification: Section and supplementary materials provide comprehensive details regard- ing GPU and CPU resources utilized, memory requirements, and explicitly report training times for all numerical experiments.
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics http:////eeeeepp.///c/c////iiiaa
Answer: [Yes]
Justification: The research presented aligns fully with the NeurIPS Code of Ethics, involving no ethical conflicts or violations.
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: The paper primarily presents methodological advancements in numerical methods without direct societal implications.
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: The paper does not involve high-risk data models requiring safeguards.
12. Licenses for existing assets
Question:
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and the license and terms of use explicitly mentioned and properly respected?
Answer: [NA]
Justification: No external code, data, or assets are explicitly mentioned used in the paper.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: The paper does not introduce new datasets, models, or assets.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: The paper does not involve crowdsoureing or human subjects research.
15. Institutional review board (IRB) approvals equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
33

---

Answer: [NA]
Justification: The research does not involve human subjects or IRB-related matters.
16. Declaration LLM usage
Question: Does the paper describe the usage of LLMs if is an important, original, or non-standard component of the core methods in this research?
Answer: [NA]
Justification: The research does not use large language models similar tools in the core methodology.
34


Reviews:
[{"review_id": "ATa2Lko2TI", "reviewer": "Reviewer_oELF", "summary": "The paper introduces a new deep learning method for computing the Legendre transform of convex and differentiable functions. In particular, the method is based on an implicit form that allows computing the Legendre transform on the codomain of the gradient of the original function. This formulation allows an efficient machine learning formulation based only on a set of samples in the original domain of the function.\nAuthors also address the problem that the measure induced by the gradient of the original function $f$ might not be well distributed in its image, leading to biased datapoints when learning the Legendre Transform; they discuss how to learn an inverse function of the gradient map and use it to obtain well-distributed samples in the domain.\nFinally, authors present numerical experiments of their method with many different network architectures, such as KANs, ResNet, MLPs, and ICNN.", "questions": "- I do not understand the purpose of Section 3: it does not show and guarantee convergence rate, just exemplification of the law of large numbers and change of variable of an integral. Is something missing or in the Supplementary material?\n- Two questions related to the sampling from a distribution in the gradient space?\n\t- Could adding a penalty if $h_\\vartheta$ outputs outside of $C$, the domain of $f$, be beneficial to the learning instead of discarding those samples?\n\t- Most importantly, I do not get why, if you have learned a representation of the inverse of the gradient $h_\\vartheta$ , you can't just use Eq. (1.2) to directly compute the Legendre Transform, instead of computing it through DLT.\n- Could you give an example where $D$ is a proper subset of the domain of $f^*$ ? In that case, how does the DLT behave outside of its training samples?\n\n**Minors**\n- You use the acronymous ICNN before you specify what it is (line 64)\n- Line 77: why are ResNet convex?\n- Line 164: Why specifically two? Does not depend on the complexity of the function?", "limitations": "yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "The method is very straightforward and, most importantly works well, on the examples provided \n- I feel the paper would improve with additional motivation on why numerical Legendre transform is beneficial, maybe with a good numerical example where this method excel. The examples provided are not too impressive.\n- The numerical methods _Deep Legendre Transform_ is compared with are for sure less performant in high-dimension, but are more general: they do not require the function to be convex and differentiable. Moreover in the numerical experiment we have a pure time comparison that can be misleading: DLT requires the computation of the gradient that can hard to compute; in situation where the $f$ (and its gradient) are not available analytically, they could account to the most expensive part of the compute.", "was_revised": true, "final_justification": "The paper is not super well written, and it could be hard to understand its relevance at a first read.\nAuthors did a good rebuttal, you can clearly tell their work is solid from the answer they gave to me and others reviewers.\n\nStill, the work is not perfect, that's why my leaning for acceptance is mild. It could become more solid from the major revision they promised.", "rebuttal": "# Response to Reviewer oELF\n\nThanks for your suggestions. The Legendre transform is one of the fundamental methods of transforming one function into another one. It is often used to switch between different formulations of a problem that highlight different variables or properties. It is particularly useful in fields like thermodynamics, mechanics, and optimization. Our paper focuses on the formal introduction of the Deep Legendre Transform algorithm, the assumptions we need and an application to optimal-control (Hamilton–Jacobi) problems. Beyond that, the same machinery is useful in several areas—for instance, computing dual utilities and profit functions in economics; evaluating convex potentials in optimal transport; building Moreau envelopes; or plug-and-play methods, where one seeks proximal operators, which are also gradients of convex functions. We will extend the motivational part.\n\n## Response to Major Comments\n\n**2.** _\"The numerical methods DLT is compared with ... are more general: they do not require the function to be convex and differentiable.\"_\n\nThis is correct, DLT intentionally targets the convex-differentiable setting. This keeps the analysis clean and provides a formal stepping-stone for future extensions. We discuss all limitations explicitly in the paper and the conclusion. Extensions, based on this approach, will be discussed in future work.\n\n**3.** _\"I do not understand the purpose of Section 3.\"_\n\nSection 3 gives an a-posteriori certificate: the Monte-Carlo estimate (Eq. 3.1) converges to $ ||g-f^*||_{L^2(D,\\nu)} $. \n\nThis guarantee makes it possible to estimate errors of the approximation of the Legendre transform  $ f^* $ even in cases where $ f^* $ is unknown (e.g. for quadratic-over-linear functions $f$). It is true that $ ||g-f^*||_{L^2(D,\\nu)} $ cannot be computed exactly. But it can be estimated up to arbitrary precision with a Monte Carlo sum if it has enough terms. We will ...\n\n1. explain and motivate the definition of the certificate in Section 3 better\n2. add a short lemma showing that the estimator is unbiased and has variance $\\sigma^2/|\\mathcal{X}_{\\mathrm{te}}|$.\n\n**4a.**  *\"Could adding a penalty if $ h_\\theta (y) \\notin C$ be beneficial instead of discarding?\"*\n\nThanks for the suggestion. We tried both and observed similar convergence but higher variance with penalization because boundary violations dominate the loss early on. Discarding yields a cleaner signal and faster training without the need to tune any penalty strength/bias. We shall report this in the appendix.\n\n**4b.** _\"Why not plug the inverse into (1.2) and be done?\"_\n\nYes, one could do that in simple setups — specifically, when the inverse-gradient map can be learned to very high precision. But then ...\n\n- convexity of the outcome is no longer guaranteed, even if an ICNN is used to approximate $f^*$,\n- the final Legendre accuracy is bottlenecked by the precision of that operator-level approximation, and it takes extra time to train the mapping,\n- after training, every new sample passes through this learned network instead of using random sampling and the gradient mapping (see the comparisons in Table 3, where the approximate inverse operator is used for sampling).\n\nOur algorithm skips that burden. We do not need to learn a generalized inverse of the gradient mapping $\\nabla f$ unless it is heavily distorting or a specific $L^2$ sampling is required (see Section 5.2).\n\n- For sampling alone, a rough estimate of the inverse gradient mapping is sufficient.\n- Thanks to the result of Section 3, the test error of DLT is measured directly against the true $ f^* $ even if $ f^* $ is unknown. In particular, there are no hidden operator errors.\n- For a provably convex and accurate $ f^* $, it is not. DLT learns a globally consistent approximation, preserves convexity by construction, and avoids error propagation from an imperfect approximation of the inverse gradient map.\n\n**4c.** _\"Example where $D \\subsetneq \\text{dom } f^*$\"_\n\nYes, for example, the quadratic-over-linear function Eq. (5.1) in line 210. We do not expect (as is common in supervised learning) that DLT approximates well outside of the domain, where it was trained.\n\n## Minor Points & Clarifications\n\n| Reviewer Comment | Clarification/Change |\n|------------------|---------------------|\n| \"ICNN used before being defined\" | Thanks, we wil move the definition up. |\n| \"Why are ResNet convex?\" | Sorry, they aren't. We meant \"ICNN ensures convexity, whereas ResNet gives a better fit but no convexity guarantee\". We will rephrase line 77. |\n| \"Why specifically two hidden layers?\" | For all ICNN/MLP experiments, we used a depth of two as a fair baseline; deeper nets yield slightly better accuracy at the cost of longer training times. |\n\nIf any aspect remains unclear, we would be grateful for further guidance; otherwise, we thank you for reconsidering the manuscript in light of these revisions."}, {"review_id": "yKei7FEQxo", "reviewer": "Reviewer_rUGB", "summary": "Authors consider a problem of finding convex conjugate for differentiable convex function. To construct approximation they propose to use neural network as parametrisation for convex conjugate and loss function that is minimised by stochastic gradient descent. Loss function is chosen such that for some convex functions approximation of inverse Jacobian matrix is required to adjust distribution on domain of convex function of interest. This approximation is achieved by a second neural network. Authors evaluate their method on several problems and find it to be competitive to other approaches.", "questions": "1. **novelty and relevance**\n\n   Between lines 66 and 67 right before equation (1.4) authors formulate a natural loss that was already considered in https://arxiv.org/abs/2210.12153 for different sets of problems. Authors consider loss function (1.4) which is different from https://arxiv.org/abs/2210.12153 and apparently require to also approximate the inverse Jacobian matrix in practical cases. Essentially the primal proposition by the authors is to use a different loss function from https://arxiv.org/abs/2210.12153. In my view this significantly undermines the novelty of the contribution. It would be possible to consider the article from the point of view of applications, but only toy problems are chosen for evaluation (see below). Given that, I would like to ask the following questions:\n\n   a. Can the authors please defend the novelty of their approach in light of https://arxiv.org/abs/2210.12153?\n\n   b. Did the authors try to compare their loss function with the one obtained from first order optimality conditions (https://arxiv.org/abs/2210.12153)? Can the authors please provide results of such comparison? \n\n2. **a posteriori guarantees**\n\n   Authors used work \"a posteriori guarantees\" which usually means that we can ensure a certain level of error. For example, in a form of inequality that bound error from above based on the quantities that are computable without the use of exact solutions. On line 147 authors appeal to the law of large numbers to derive their error estimation.\n\n   a. Can the authors please elaborate how one can use the law of large numbers to obtain an upper bound on error?\n\n   b. Should not authors require some properties of probability distribution for the law of large numbers to hold?\n\n3. **weak benchmarks**\n\n   The first batch of examples is given in Table 1. In all cases functions are completely separable, i.e., both function and its conjugate are sums of one dimensional functions. This is also true for the chosen Hamilton–Jacobi equation considered in Section 5.3. The only non-separable function used to test the method is (5.1). It is clear that high-dimensional separable functions are easier to approximate (e.g., hyperbolic cross). Can the authors increase the number of non-trivial benchmarks?\n\n4. **misc**\n\n   a. How $L_2$ errors reported in Table 1 are computed?\n\n   b. How $L_2$ errors reported in Table 3 are computed? \n\n   c. The table on the right panel of Figure 2 contains the claim that authors computed $L_2$ which is not true since convex conjugate is not available.\n\n   d. How direct learning of high-dimensional functions is implemented?\n\n   e. Lines (204-205) \" The ResNet consistently provides the best approximation results despite not guaranteeing convexity ... \" Did the authors try to quantify convexity of learned functions?\n\n   f. Line 78 \"... resulting approximation of $f^{*}$ is automatically convex.\" Potentially misleading since convexity only holds for ICNN not the ResNet.\n\n   g. Equation between lines 66 and 67 right before equation (1.4): right bracket is missing.\n\n   h. Line 47, 48 \"A different approach, going back to a tweet by Peyré [2020]\" and in the footnote \"see e.g., Nesterov (2005), Beck and Teboulle (2012), Niculae and Blondel (2017) or Mensch and Blondel (2018)\" So is it from 2020 tweet or from 2005 Nesterov book (or article)? Why these references are given in plain text?", "limitations": "yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "**Strengths**\n1. Ideas proposed by authors are easy to understand\n2. Exhaustive context on the computation of convex conjugate is available which makes the article suitable for general audience\n3. Appendix contain thorough comparison with classical methods\n\n**Weaknesses**\n1. Novelty and relevance of the method are questionable\n2. Method lacks a posteriori approximation guarantees despite the claim of the authors\n3. Problems used for benchmarks are too simple\n\nSee below for details.", "was_revised": true, "final_justification": "Authors addressed my main concerns:\n1. The learning problem proposed in the papers seems to be original enough.\n2. Results for non-separable functions provide more challenging evaluation.\n\nThere are still some outstanding questions on guaranties and convexity, but in my view the quality of the paper is sufficient for publication.", "rebuttal": "# Response to Reviewer rUGB\n\n## 1. Novelty and Comparison with Existing Work\n\n**a.** _\"Essentially the proposition is just a different loss—please defend the novelty.\"_\n\nThanks, it's a good question. The paper \"On amortizing convex conjugates for optimal transport\" that you mention discusses different ways of solving the problem $x(y) = \\arg\\sup \\lbrace (x,y) - f(x) \\rbrace $, which, in our setup, amounts to approximating the inverse of the gradient mapping $\\nabla f$. This yields several interesting algorithms. But our approach is different. In particular, it avoids the approximation of the inverse of the gradient mapping $\\nabla f$ and the resulting approximation errors. In addition, if used with ICNNs as approximators, our approach ensures that the final approximation of the Legendre transorm is convex, as it should be. We summarize the differences between previous and our approaches in the following table:\n\n| Aspect | 2210.12153 (\"FOC-loss\") | Our work (DLT-loss) |\n|--------|-------------------------|---------------------|\n| Minimised identity | $\\|\\nabla f(h_\\theta(y)) - y\\|^2$ (first-order condition) | $[g_\\theta(\\nabla f(x)) + f(x) - \\langle x, \\nabla f(x) \\rangle]^2$ (Fenchel identity) |\n| What is learned | Inverse map $h_\\theta: D \\to C$ | Direct conjugate $g_\\theta: D \\to \\mathbb{R}$ |\n| Convexity guarantee | none (unless post-processed) | automatic when $g_\\theta$ is an ICNN |\n| A-posteriori error | not available (needs true $f^*$) | unbiased Monte-Carlo estimator (Sec. 3) |\n| Use as sampler | not covered | built-in (Sec. 4) |\n\nThus, DLT is not merely a \"different loss\"; it changes the target object (function conjugate instead of operator inverse), provides a certified error estimator, and preserves convexity, which the FOC-loss does not.\n\nPotentially, one can use $x_\\theta$ for and substitute it into Eq. (1.2), but this is suboptimal because \n- Convexity is no longer guaranteed, even if an ICNN is used as approximator.\n- The final Legendre accuracy is bottlenecked by the precision of that operator-level approximation, and it takes extra time to train the mapping.\n- After training, every new sample passes through this learned network instead of using random sampling and gradient mapping.\n\n**b.** _\"Did you compare losses?\"_\n\nYes, we did and can provide the results. In cases where the inverse gradient mapping exists and can be accurately approximated, the losses are comparable, otherwise the accuracy is bottlenecked by the precision of the approximation of the inverse gradient mapping. Moreover, convexity guarantees are lost and training times are longer (see Table 3) compared to DLT, which does not need to learn the inverse gradient mapping. \n\n## 2. A-posteriori Error Estimator\n\n_\"How can the law of large numbers provide an upper bound?\"_\n\nThanks for the question. Strictly speaking, it cannot compute the bound exactly. But it's very common to approximate complicated integrals with Monte Carlo sums. This provides unbiased estimators with confidence intervals that can be made arbitrality small if enough terms are used in the Monte Carlo sum. This provides accurate results at a relatively cheap computational cost.\n\n$ \\widehat{E} = \\frac{1}{|\\mathcal{X}_{\\text{te}}|}  \\sum_x   \\left[g(\\nabla f(x)) + f(x) - \\langle x, \\nabla f(x) \\rangle\\right]^2 \\quad \\Longrightarrow \\quad \\mathbb{E}[\\widehat{E}] = ||g - f^* ||^2_2 \\quad  \\text{  in } L^2(D,\\nu)$\n\nand\n\n$$\\text{Var}[\\widehat{E}] = \\sigma^2 / |\\mathcal{X}_{\\text{te}}|$$\n\nOf course, we need to assume that the $\\{x_i\\}$ are i.i.d. and integrable to make sure the Monte Carlo sum converges. For the confidence intervals we need square-integrability. This assumptions are satisfied in our examples. Thank you for pointing it out. We will explain this better in the paper.\n\n## 3. Benchmark Richness\n\nThanks for the suggestion. Our original paper over-emphasized separable functions. We mainly employed them since they have explicit Legendre transforms, which can be used as benchmarks. But our method also works for non-separable functions. As you mentioned, quadratic-over-linear functions are non-separable. We can e.g. add the following non-separable examples:\n\n1. **Quadratic form with random SPD matrix**: $f(x) = \\frac{1}{2} x^\\top Q x$.\n2. **Coupled soft-plus**: $f(x) = \\sum_{i<j} \\log\\left(1 + e^{x_i + x_j}\\right)$.\n3. **Exponential-minus-linear**: $f(x) = e^{\\langle a, x \\rangle} - \\langle b, x \\rangle$.\n\nIn all these examples, DLT has an $L^2$-error below $10^{-1}$ in up to 200 dimensions while tradiitional grid/FLT methods are intractable already in 10 dimensions.\n\n## 4. Miscellaneous Questions\n\n| Reviewer Question | Answer / Fix |\n|------------------|--------------|\n| **a./b.** How are $L^2$ errors (Tbl. 1 & 3) computed? | Monte-Carlo over 4096 i.i.d. test samples from the distribution $\\nu = \\text{Unif} \\circ (\\nabla f)^{-1}$; details will be added to the caption. |\n| **c.** Fig. 2 table claims \"computed $f^*$\" | The true $f^*$ is unknown there, but due to the results of section 3 the error is identical to the one obtained with the direct learning. |\n| **d.** \"Direct learning\" implementation | Same network architecture, trained on pairs $(y, f^*(y))$, where closed form is available. |\n| **e.** Did we quantify convexity loss in ResNet? | Not explicitly, but it is well known that deep NNs can give a wiggling approximation. |\n| **f.** Line 78 wording | Changed to \"if $g_\\theta$ is an ICNN the approximation is convex\". |\n| **g.** Missing bracket | Fixed. |\n| **h.** Tweet vs. Nesterov book | Sorry, the entropic log-sum-exp smoothing of a max can be traced back to the Gibbs variational principle; it was brought into modern convex-optimisation analysis by Nesterov (2005) and popularised for GPU-friendly Legendre transforms by Peyré's 2020 tweet. We will explain it more carefully and provide proper references. |"}, {"review_id": "eRe4L4ezCI", "reviewer": "Reviewer_gyT6", "summary": "This paper proposes a novel deep learning approach for approximating the convex conjugates of differentiable convex functions. The core algorithm leverages an efficient gradient-based framework to minimize approximation error, utilizing the implicit Fenchel formulation of convex conjugation. The numerical experiments show the effectiveness of the proposed method. Besides that, the implicit Fenchel formulation also facilitates a posteriori approximation guarantees.", "questions": "Please refer to **Strengths And Weaknesses** for my questions.", "limitations": "Please refer to **Strengths And Weaknesses** for the limitations.", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Summary of Claims and Evaluation**\n\n**Claim 1**: Training a machine learning model to learn the convex conjugate f* yields strong numerical performance compared to directly learning f* from its analytic form.\n\nThis claim is supported by Table 2, which shows that the L2-approximation error of the proposed DLT method is comparable to that of direct learning using the analytic expression of f*. This also serves as an informative ablation study for evaluating the suitability of different neural network architectures and validates the proposed learning formulation (1.5).\n\n**Claim 2**: Using ICNN or ResNet as the reparameterization network enables scaling to high-dimensional problems and ensures that the resulting approximation of f* is convex.\n\nI have two concerns regarding this claim:\n- The experiments only go up to 50 dimensions (e.g., as shown in Table 1), which may not be sufficiently high to convincingly support claims about scalability to high dimensions.\n\n- Line 94 states that convexity is not guaranteed when using ResNet, which appears to contradict Claim 2’s assertion that the approximation is automatically convex. This inconsistency should be clarified.\n\n\n**Claim 3**: KAN (Kernel Activation Network) used as the reparameterization network can recover exact solutions for certain convex functions.\n\nThis claim is supported by the results in Table 5. However, the experiments are limited to low-dimensional (e.g., 2D) problems, which may limit its practicality. \n\n**Claim 4**: The implicit Fenchel formulation enables a posteriori approximation error estimation.\n\nThis claim is supported by the main theoretical result in Section 3. The result appears sound, although it falls outside my current area of expertise, so I cannot fully assess its correctness.\n\n**Strengths**\n\n- The paper is clearly written and well organized.\n\n- The authors' claims are largely supported by numerical results.\n\n- The experimental evaluations are extensive and clearly presented.\n\n**Weaknesses**\n\n- The proposed method is introduced somewhat abruptly in both the abstract and introduction. I suggest expanding the first paragraph of the introduction to better motivate the significance and applications of computing convex conjugates.\n\n- The experiments only go up to 50-100 dimensions (e.g., as shown in Table 1), which may not be sufficiently high to convincingly support claims about scalability to high dimensions.\n\n- Line 94 states that convexity is not guaranteed when using ResNet, which appears to contradict Claim 2’s assertion that the approximation is automatically convex. This inconsistency should be clarified.\n\n- The KAN experiments in Table 5 are limited to low-dimensional (e.g., 2D) problems, which may limit the practicality. Please consider elaborating on how this approach might generalize to higher dimensions.\n\n- Minor comment: It would also be helpful to report uncertainties for metrics such as approximation errors and computation time, especially since the latter can vary significantly across runs due to hardware-level fluctuations.", "was_revised": true, "final_justification": "Some of my concerns have been addressed, and I remain overall positive about the paper. However, I still have concerns about the practicality of the proposed method. The experiments are limited to problems with up to 200 dimensions, which is significantly smaller than many real-world scenarios. Therefore, I cannot justify a substantial increase in score and retain my original rating of Borderline Accept.", "rebuttal": "# Response to Reviewer gyT6\n\n1. The Legendre transform is a classical and widely used tool for transforming functions and changing perspectives, which is widely used in thermodynamics, mechanics, and optimisation. Our paper formally introduces the **Deep Legendre Transform** (DLT), states the precise assumptions, and showcases its use in a Hamilton–Jacobi optimal-control problem. The same machinery is also used to derive dual utilities in economics, convex potentials in optimal transport, Moreau envelopes as well as plug-and-play proximal maps. We will expand this motivational part in the introduction.\n\n2. Scalability: our. approach scales well, we have included experiments in up to 200 dimensions in the supplements (Table 13).\n\n3. Convexity: sorry, we have to explain this more carefully. ICNNs lead to convex approximations; ResNets provide good approximations, but are in general not convex.\n\n4. KAN experiments: yes, our current examples are 2-dimensional. Our goal was to illustrate that in certain cases, the implicit Eq. (1.4) can yield exact results via symbolic regression, which we find interesting. In principle, there is no reason why this shouldn't work in higher dimensions. But it will depend on the implementation of KAN, which is a recent invention.\n\n5. Thanks, we report uncertainties (mean ± s.d.) throughout the supplement (see  Table 4 (main text) and Tables 9 and 10); it will be extended and explained more carefully."}, {"review_id": "gxJK0royV7", "reviewer": "Reviewer_5uuE", "summary": "The authors train a neural network to approximate the convex conjugate (Legendre–Fenchel transformation) of convex differentiable functions. Indeed, the convex conjugate $ f^* $ of $ f $ satisfies in this case $f^* (\\nabla f(x)) = \\langle x, \\nabla f(x) \\rangle - f(x)$, and therefore they minimise the distance between the output of the neural network and the expression $\\langle x, \\nabla f(x) \\rangle - f(x)$. However, this approximates the convex conjugate on $D = \\nabla f(C)$ only, where $C$ is the domain of $f$, and uniformly sampling training points in $C$ before transforming them through $\\nabla f$ sometimes leads to an uneven distribution over $D$. The authors address this second point by learning another neural network to approximate an inverse of $\\nabla f$, which allows to sample uniformly from $D$ and to then go back to $C$. However, $D$ itself is not always explicitly known as in the example in Figure 1 in the paper, which makes this solution impossible in this case, and makes it hard to know how much of the domain of $f^* $ we are covering.", "questions": "See point 2 in weaknesses.", "limitations": "Weakness 2 could be discussed further", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths:\n1. Well-written\n2. Theoretically sound\n3. Multiple architectures were tested for the approximating neural network\n\nWeaknesses:\n1. The method is not compared to any other method (except in the case of its application to the HJ equation, and to a naive learning of a neural network to directly approximate the conjugate when it is available is closed form). The authors mention (and test in the supplementary material) some existing non-deep learning methods for approximating convex conjugates, but do not compare them to their method. These numerical approximation methods do not scale to high dimensions, but they could have tested their method in lower dimensions to compare.\n2. Even though the fact that their method only learns an approximation of the convex conjugate on $D = \\nabla f(C)$ is acknowledged, I think that this point warranted more discussion. How big is this set compared to the domain of $ f^* $ in typical/useful cases? How does the learned neural network approximate $ f^* $ outside of $D$ anyway?\n3. A few more tests and comparisons on problems such as the HJ equation where convex conjugates appear could have been included.", "was_revised": true, "final_justification": "The rebuttals to the reviewers have been of a high quality therefore I have raised my score to 4. However, the proposed changes to the paper are so extensive that it is hard to evaluate the resulting quality of the paper and to raise the score to a strong accept.", "rebuttal": "# Response to Reviewer 5uuE \n\n## 1. Classical vs DLT in Low Dimensions\n\nThanks, the results in Tables 2 (main + supplementary) compare run times and errors. Classical grid/LLT are indeed faster in dimensions 1-4; in 6-7 dimensions, they become slower and less accurate owing to coarser grids. In more than 8 dimensions, they run out of memory while DLT still trains in seconds.\n\n## 2. Size of D\n\nFor essentially smooth convex functions (see Rockafellar 1970), one has\n\n$$\n\\text{ri}(\\text{dom} f^*  ) \\subseteq D  \\subseteq \\text{dom} f^*\n$$\n\nIn particular, any convex smooth function defined on $\\mathbb R^{d}$ is essentially smooth, so its entire effective dual domain is covered.  We have written this paper as a solid stepping‑stone for future extensions that will treat more general settings such as non-smooth functions and functions that are not defined on all of $\\mathbb R^{d}$ or explode at some boundary.\n\n## 3. Additional Benchmarks\n\nThanks for the suggestion, as mentioned to Reviewer rUGB, we probably put too much emphasis on separable functions. But our approach clearly also works for non-separable functions.\n\nWe will add more non-separable examples to the supplements to complement the existing separable and quadratic-over-linear examples."}]

Author rebuttal:

Meta-review: The paper proposes an algorithm for computing the Legendre-Fenchel conjugate of a function using deep learning. In contrast to previous approaches, it does not require the inverse gradient mapping and it guarantees that the resulting learned conjugate function is convex (as it should in this context) by using Input-Convex Neural Networks.

The reviewers were on the fence initially with many borderline scores. After the rebuttal period, many of the concerns of the reviewers were addressed by the authors, in particular explaining the originality compared to prior works and the improvement with respect to dimensionality of the problems. As a result, the scores were increased but still borderline - no strong accepts (5) were given.

In the post-rebuttal discussion, I explained to the reviewers my opinion that the paper was on the fence and likely to be rejected. To this, Reviewer oELF expressed that the paper should be accepted and is worthwhile. I also checked the final justification of Reviewer rUGB and they have also suggested that the paper is worthy of publication in their final assessment.

So, although none of the reviewers were willing to give a 5 score, it is now clear that multiple were ready and willing to defend the paper as worthy of acceptance; I therefore recommend to accept the paper.
Decision: Accept (poster)