A provable control of sensitivity of neural networks through a direct parameterization of the overall bi-Lipschitzness

Authors: yuri  kinoshita
taro  toyoizumi
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-13 10:21:48

Abstract: While neural networks can enjoy an outstanding flexibility and exhibit unprecedented performance, the mechanism behind their behavior is still not well-understood. To tackle this fundamental challenge, researchers have tried to restrict and manipulate some of their properties in order to gain new insights and better control on them. Especially, throughout the past few years, the concept of *bi-Lipschitzness* has been proved as a beneficial inductive bias in many areas. However, due to its complexity, the design and control of bi-Lipschitz architectures are falling behind, and a model that is precisely designed for bi-Lipschitzness realizing a direct and simple control of the constants along with solid theoretical analysis is lacking. In this work, we investigate and propose a novel framework for bi-Lipschitzness that can achieve such a clear and tight control based on convex neural networks and the Legendre-Fenchel duality. Its desirable properties are illustrated with concrete experiments to illustrate its broad range of applications.

Main paper: 
 introduction 1.background: 
 Nowadays, neural networks have become an indispensable tool in the field of machine learning and artificial intelligence. While they can enjoy an outstanding flexibility and exhibit unprecedented performance, the mechanism behind their behavior is still not well-understood. To tackle this fundamental challenge, researchers have tried to restrict and manipulate some of their properties in order to gain new insights and better control on them. Especially, throughout the past few years, the concept of sensitivity has been proved as a beneficial inductive bias in many areas.Sensitivity can be translated into the concept of bi-Lipschitzness, which combines two different properties, namely, Lipschitzness and inverse Lipschitzness. The former describes the maximal, and the latter the minimal sensitivity of a function. Bi-Lipschitzness is attracting interests in various fields not only as it proposes a promising solution to avoid unexpected and irregular results caused by a too sensitive or too insensitive behavior of the trained function, but also as it achieves an approximate isometry that preserves geometries of the input dimension (Li et al., 2020) . It plays an essential role in generative models such as normalizing flows to guarantee invertibility (Behrmann et al., 2019) , in uncertainty estimation to prevent the recurrent problem of feature collapse (Liu et al., 2020a; Van Amersfoort et al., 2020) and in inverse problems (Kruse et al., 2021) to assure the stability of both the forward and inverse function (Behrmann et al., 2021) .Unfortunately, the appropriate design and effective control of bi-Lipschitz neural networks are far from simple, which hinders their application to prospective areas. First of all, the estimation of bi-Lipschitz constants is an NP-hard problem (Scaman and Virmaux, 2018) . Second, the design of bi-Lipschitz models is even harder as we cannot straightforwardly extend existing Lipschitz architectures that exploit some unique properties of the concept to inverse Lipschitzness and keep their advantages, and vice-versa. Finally, the control of bi-Lipschitz constants requires particular attention as it is a question of manipulating two distinct concepts in a preferably independent and simple manner.Currently existing bi-Lipschitz models still present some issues in terms of design or control. On the one hand, some lack theoretical guarantees because they impose soft constraints by adding regularization terms to the loss function (Van Amersfoort et al., 2020) , or because their expressive power is not well-understood and may be more limited than expected. On the other hand, others restrict bi-Lipschitzness on a layer-wise basis (Behrmann et al., 2019; Liu et al., 2020a) . Particularly, this means these approaches can only build in the essence a simple bi-Lipschitz function employed as a layer of a more complex neural network. In practice, this kind of parameterization impacts the generalization ability of the model or leads to loose control (Fazlyab et al., 2019) . They can also contain so many parameters affecting sensitivity to the same extent that controlling all of them is unrealistic and fixing some may affect the expressive power. See Section 2 for further details.Therefore, taking into account both the potential and complexity of this inductive bias, we first and foremost need a model that is precisely designed for bi-Lipschitzness realizing a direct and simple control of the constants of the overall function more complex than a single layer neural network equipped with solid theoretical analysis. In this work, we investigate an architecture that can achieve such a clear and tight control and apply it to several problem settings to illustrate its effectiveness.Our contributions can be summarized as follows. First, we construct a model bi-Lipschitz by design based on convex neural networks and the Legendre-Fenchel duality, as well as a partially bi-Lipschitz variant. This architecture provides a simple, direct and tight control of the Lipschitz and inverse Lipschitz constants through only two parameters, the ideal minimum, equipped with theoretical guarantees. These characteristic features are illustrated and supported by several experiments including comparison with prior models. Finally, we show the utility of our model in concrete machine learning applications, namely, uncertainty estimation and monotone problem settings and show that it can improve previous methods.Organization In Section 2, we will first explain in more detail existing architectures around bi-Lipschitzness. In Section 3, we will develop our model followed by theoretical analyses. The next Section 4 will be devoted to experiments and applications of our proposed method.Notation Throughout this paper, the Euclidean norm is denoted as ∥ • ∥ for vectors unless stated otherwise. Similarly, for matrices, ∥ • ∥ corresponds to the matrix norm induced by the Euclidean norm. For a real-valued function F : R m → R, ∇F is defined as (∂f (x)/∂x 1 , . . . , ∂f (x)/∂x n ) ⊤ , and ∇ ⊤ F as its transpose. When the function is a vector F : R m → R n , then its Jacobian is defined as ∇ ⊤ F = (∂f i /∂x j ) i,j . The subderivative is denoted as ∂ sub . 
preliminaries: 
 In this section, we explain mathematical backgrounds and existing bi-Lipschitz models to clarify the motivation of this work.Let us first start by the definition of bi-Lipschitzness.Definition 2.1 (bi-Lipschitzness). Let 0 < L 1 ≤ L 2 . f : R l → R t is (L 1 , L 2 )-bi-Lipschitz if L 1 ∥x -y∥ ≤ ∥f (x) -f (y)∥ ≤ L 2 ∥x -y∥ holds for all x, y ∈ R l . The right (left) inequality is the (inverse) Lipschitzness with constant L 2 (L 1 ). L 1 and L 2 are called together bi-Lispchitz constants. Several interpretations can be attributed to this definition. Bi-Lipschitzness appears when we want to guarantee a high sensitivity and a high smoothness for better generalization such as in uncertainty estimation. This can be also regarded as a quasi-isometry with a distortion of L 2 /L 1 , which means that the structure of the input is relatively inherited in the output as well. Moreover, since a bi-Lipschitz function is by definition invertible, the inverse Lipschitz can be interpreted as the Lipschitz constant of the inverse function. Therefore, controlling the bi-Lipschitz constants is beneficial for the stability of both the forward and inverse.In order to gain new insights of a problem and a predictable behavior of the bi-Lipschitz model, it is primordial that bi-Lipschitzness works as an inductive bias with high controllability and theoretical foundation. For this goal, the following points are of particular interest: (a) bi-Lipschitzness guaranteed by design, (b) theoretical bounds on the bi-Lipschitz constants, (c) tight bounds, (d) independent control of the Lipschitz and inverse Lipschitz constants, (e) a minimal number of hyperparameters to tune, (f) theoretical guarantee on the expressive power and (g) direct parameterization of the optimized variables (in the sense of Definition 2.3 of Wang and Manchester (2023) ).There are currently three methods that are mainly employed to achieve and control bi-Lipschitz neural networks. See Appendix A for a brief overview on Lipschitz and inverse Lipschitz architectures.Regularization The first one is based on a regularization method (Gulrajani et al., 2017; Van Amersfoort et al., 2020) . It adds to the loss a term that incites bi-Lipschitzness. While the method is simple, the resulting function may not be bi-Lipschitz at the end of the training, we have no theoretical bounds on the bi-Lipschitz constants and cannot achieve a real control of them.i-ResNet The invertible residual network (i-ResNet) of Behrmann et al. (2019) is based on the composition of several (1 -L gi , 1 + L gi )-bi-Lipschitz layers of the type f i (x) = x + g i (x), where g i is Lipschitz with constant L gi < 1. The Lipschitz g i can be constructed by layer-wise methods such as spectral normalization of weights (SN) (Miyato et al., 2018) . However, this kind of layer-wise control of (bi-)Lipschitzness is known to be loose in general (Fazlyab et al., 2019) . See Figure 1 for an illustration. The i-ResNet has limited expressive power as it cannot represent the function y = -x (Zhang et al., 2020) , and the construction of g i risks to introduce more restriction on the expressive power in practice (Anil et al., 2019) . A similar model, the Lipschitz monotone network (LMN) of Nolte et al. (2023) combines the GroupSort activation function with SN to create g i . These approaches requires to adjust the Lipschitzness of each weight during the training process, which may be sub-optimal (i.e., no direct parameterization).BiLipNet The concurrent work of Wang et al. (2024) also provides a bi-Lipschitz neural network (BiLipNet), which was mainly used for creating a Polyak-Łojasiewicz function useful in surrogate loss learning. It extends the Lipschitz architecture of Wang and Manchester (2023) . The BiLipNet is constructed as the composition of monotone Lipschitz layers and orthogonal layers. The former layer is realized based on a direct parameterization of the IQC theory of Megretski and Rantzer (1997) . While they can achieve a certified control of the bi-Lipschitz constant, theoretical guarantee on its expressive power is lacking. Moreover, BiLipNet composes many layers which may harm the tightness of the bounds.Therefore, current models fail to satisfy the desired features of Subsection 2.2. Notably, all these points cannot be satisfied as long as we rely on the simple layer-wise control which offers too conservative bounds and a number of tunable hyper-parameters proportional to that of the layers in general. In this work, we will establish a bi-Lipschitz architecture that satisfies almost all these expected attributes mentioned above. In terms of mathematical foundations, it takes inspiration from the inverse Lipschitz architecture of Kinoshita et al. (2023) . 
bi-lipschitz neural network: 
 We introduce our novel bi-Lipschitz model, followed by some theoretical analyses and discussion.We clarify notions of smoothness, convexity and the Lengendre-Fenchel transformation (LFT) which will be a core process in our construction. Definition 3.1. Let γ > 0. F : R m → R is γ-smooth if F is differentiable and ∇F is γ-Lipschitz. Definition 3.2. Let µ > 0. F : R m → R is µ-strongly convex if F (x) -µ 2 ∥x∥ 2 is convex. Definition 3.3. Let F : I → R a convex function over I ⊂ R m . Its Legendre-Fenchel transformation F * is defined as F * (x) := sup y∈I {⟨y, x⟩ -F (y)}.We extend the approach of Kinoshita et al. (2023) and propose a method that creates bi-Lipschitz functions, where it is sufficient to control two parameters to manipulate the overall bi-Lipschitzness without needing to dissect the neural network and control each layer, a feature that existing methods for (bi-)Lipschitzness cannot deliver. The first step is to notice that the gradient of a real-valued µ-strongly convex function becomes µ-inverse Lipschitz, and that of a γ-smooth function γ-Lipschitz by definition. Therefore, we aim to compose a function which is both strongly convex and smooth. Interestingly, the LFT of a 1/β-strongly convex function is β-smooth. This leads to our main theorem. Theorem 3.4. Let F be a closed 1/β-strongly convex function and α ≥ 0. Then the following function is α-strongly convex and α + β-smooth: sup y∈I {⟨y, x⟩ -F (y)} + α 2 ∥x∥ 2 . Thus, its derivative is (α, α + β)-bi-Lipschitz which equals f * (x) := argmax y {⟨y, x⟩ -F (y)} + αx.See Appendix B.2 for the proof. The term α/2∥x∥ 2 has the effect of turning a convex function into an α-strongly convex function by definition of strong convexity. This f * (x) constructed as described in the above theorem is precisely the bi-Lipschitz model we propose in this paper.Based on Theorem 3.4, we can create a bi-Lipschitz function parameterized by neural networks, and thus applicable to machine learning problems. Here, we clarify the implementation of a strongly convex function and the LFT. The overall explicit formulation of our bi-Lispchitz neural network (BLNN) is already shown in Algorithm 1 for convenience.Strongly Convex Neural Networks A µ-strongly convex neural network can be constructed by adding a regularization term to the output of the Input Convex Neural Network (ICNN) from Amos et al. (2017) . The resulting structure can be written asF θ (y) = z k + µ 2 ∥y∥ 2 , where z i+1 = g i (W (z) i z i + W (y) i y + b i ) (i = 0, . . . , k -1). {W (z) i } i are non-negative, W (z) 0= 0, and all functions g i are convex and non-decreasing. This architecture is a universal approximator of µ-strongly convex functions defined on a compact domain endowed with the sup norm (Chen et al., 2019) . Note that any other choice for the convex architecture is possible. This is only an example, also employed in Huang et al. (2021) and Kinoshita et al. (2023) . Indeed, we could opt for a convolutional version by using the convolutional ICNN proposed by Amos et al. (2017) .Algorithm 1: Forward pass of (α, β)-BLNN Input: input data x, an ICNN G θ with parameters θ and constants α, β ≥ 0 Output : image of x by an (α,α + β)-bi-Lipschitz function Step 0: Construct 1 β -strongly convex function F θ (y) = G θ (y) + 1 2β ∥y∥ 2 Legendre-Fenchel transformation Step: Find y * θ (x) = argmax y {⟨y, x⟩ -F θ (y)} Gradient Step: Compute f * θ (x) = ∇ x F * θ (x) + α 2 ∥x∥ 2 = y * θ (x) + αx return f * θ (x)Algorithms for LFT By Theorem 3.4, we only need to compute the optimal point of the LFT, i.e., y * θ (x) := argmin y {F θ (y) -⟨y, x⟩}, which is a strongly convex optimization. This computation is thus rather fast, and we can use various algorithms with well-known convergence guarantees such as the steepest gradient descent (Shamir and Zhang, 2013; Bansal and Gupta, 2017) . Such convex solvers will generate for a fixed x a sequence of points {y t (x)} t=0,...,T , and its last point will be an estimate of y * θ (x). However, this kind of discrete finite time approximation could compromise the bi-Lipschitzness of the whole algorithm. The bi-Lipschitz behavior of y t (x) can be explicitly described for the steepest gradient descent as follows. For other algorithms, see Appendix C.2.1. Theorem 3.5. Let the symbols defined as in Algorithm 1. Consider the steepest gradient descent of sup y {⟨y, x⟩ -F θ (y)} generating points {y t (x)} t at the t-th iterations and y * θ (x) is the global maximum. If F θ is µ-strongly convex and γ-smooth then the point lim t→∞ y t (x) is (1/γ, 1/µ)-bi-Lipschitz without any bias. Moreover, with η t = 1/(µ(t + 1)) as a step size and y 0 (x i ) = y 0 as initial point, then for all x i , x j , ∥y t+1 (x i ) -y t+1 (x j )∥ ≤ h(t)∥x i -x j ∥ where lim t→∞ h(t) = 1/µ. See Appendix C.2.2 for the concrete formulation of h(t) and the proof. The above theorem guarantees that the optimization scheme will ultimately provide a y t (x) that is bi-Lipschitz with the right constants. This was not trivial as a bias may have persisted due to the non-zero step size. More importantly, this statement offers a non-asymptotic bound for Lipschitzness. This is greatly useful when we theoretically need to precisely assure a certain degree of Lipschitzness such as in robustness certification against adversarial attacks (Szegedy et al., 2013) . The step size η t = 1/(µ(t + 1)) can be precisely calculated as µ = 1/β for an (α, β)-BLNN. Experimental results suggest a similar behavior for inverse-Lipschitzness. Note the strong convexity of F is always satisfied in our setting and smoothness can be fulfilled if we use an ICNN with softplus activation functions log(1 + e x ). Since this simple gradient descent has been thoroughly analysed and assured to properly behave in practice as well, we will employ it as the algorithm of LFT in the remainder of this paper.As we mentioned earlier, layer-wise approaches realize a bi-Lipschitz neural network by restricting the sensitivity of each layer. Nevertheless, this kind of construction may be sub-optimal as it does not take into account the interaction between layers, limiting more than expected the expressive power, which was often not considered in the original papers. Layer-wise bi-Lipschitz approaches cannot inherit the proofs and properties of the original network, which makes the understanding of their behavior more difficult. As for our model, the (α, β)-BLNN, we can guarantee the following universality theorem. Theorem 3.6. For any proper closed α-strongly convex and α + β-smooth function on a compact domain, there is a BLNN without taking the gradient at the end, i.e., sup y {⟨y, x⟩ -(G θ (y) + ∥y∥ 2 /(2β))} + α∥x∥ 2 /2 where G θ is an ICNN with ReLU or softplus-type activation function, that approximates it within ϵ in terms of the sup norm.Thus, after taking the gradient, we can create a sequence of functions that converges point-wise to any function which is α + β-Lipschitz, α-strongly monotone (Definition B.12) and the derivative of a real-valued function. See Appendix B.3 for the proofs and further discussion.The most straightforward way to compute the gradient of a BLNN is to track the whole computation of the forward pass including the optimization of the LFT and back-propagate through it. However, this engenders a crucial bottleneck since back-propagating over the for-loop of the optimization involves many computations of the Hessian and the memorization of the whole computational graph. Nevertheless, if the activation function of the ICNN is chosen as softplus, the convergence of the LFT is quite fast making this strategy scalable to moderately large data.Interestingly, when F θ := F (•; θ) is C 2 , the backward pass can be computed only with the information of the core ICNN, which means we do not need to track the whole forward pass involving many recurrent computations: Theorem 3.7. Suppose a loss function L : z → L(z), and the output of the BLNN is f* (x; θ) := ∇F * θ (x) + αx as defined in Algorithm 1. If F is C 2 and F * is differentiable, then the gradient ∇ ⊤ θ L(f * (x; θ)) can be expressed as -∇ ⊤ z L(z) ∇ 2 y F (y * (x; θ); θ) -1 ∂ ⊤ θ ∇ y F (y * (x; θ); θ). See Appendix B.4for the proof and formulation with more complex architectures. In practice, coding is simple since we can directly employ the backward algorithm provided in libraries such as Pytorch (Corollary B.19) . The advantage of this second method is threefold. First, it can reduce the computational and memory cost of the backward process as it does not depend on the iteration number and involves fewer matrix manipulations. Second, this leaves the freedom to choose the optimizer to calculate the LFT, while the brute force method requires the gradient to be tracked during the optimization process, which is a feature that many current solvers do not provide. Third, this method approximates the gradient by taking the derivative along the true curve at a point y(t) θ (x) close to y * θ (x), which is a fair approximation. In contrast, the brute force method approximates the gradient based on the recurrent computation of y (t) θ (x) but we do not know whether this is really a good estimate.Interestingly, our BLNN can be regarded as a bi-Lipschitz version of a broader family of models called deep equilibrium models (DEQs) from Bai et al. (2019) . A DEQ is an implicit model whose output z is defined as the fixed point of a neural network h θ , i.e., z = h θ (z, x), where x is the input. Similarly, the output of our algorithm can be re-formulated as finding the solution z of z = x -∇F θ (z) + z which corresponds to a DEQ with h θ (z, x) = x -∇F θ (z) + z. The general properties of a DEQ concerning the computational complexity (both time and space) are thus also inherited in our model. One major difference is that the iteration to find this fixed point z is in our case unique and guaranteed to converge, which addresses one of the main concerns of general DEQs. We believe this interpretation is promising for future work to increase the generality of our model as it enables to converge to the larger flow of work around DEQs and apply the various improvements for DEQs researched so far. However, we will not pursue this direction further since our primary goal is to develop a bi-Lipschitz model with direct and simple control.Despite all the theoretical guarantees mentioned above and the correspondence of the proposed method to DEQs, its main drawbacks persist in computational efficiency and expressive power. However, these weaknesses can be alleviated by imposing the bi-Lipschitz requirement on a limited number of variables. This can be realized by using the partially input convex neural network (PICNN) instead of the ICNN in our architecture (Amos et al., 2017) . A PICNN is convex with respect to a pre-defined limited set of variables. Based on this architecture, we can proceed similarly to Algorithm 1 and obtain a partially bi-Lipschitz neural network (PBLNN). See Appendix F.2 for further details. As a result, all heavy operations such as the LFT and the gradient are applied on this smaller set of variables, which makes the architecture much more scalable to higher dimensions. Moreover, a PICNN with k layers can represent any ICNN with k layers and any purely feedforward network with k layers, as shown in Proposition 2 of Amos et al. (2017) . This also enhances the expressive power of our model with larger liberty on the precise construction of the architecture.Regarding the time and space complexity of our model, it is largely equivalent to that of a DEQ (Bai et al., 2019) except that the core function is the derivative of a neural network, which adds some computational burden. In Figure 2 , we present a comparison of the computational complexity in floating point number operations (FLOPs) and memory for a single iteration between a traditional feedforward network and various BLNN variants: BLNN (with brute force backpropagation), BLNN with Theorem 3.7, PBLNN with only one variable constrained to be bi-Lipschitz. Comparing our model with Theorem 3.7 to a traditional feed-forward neural network, we can conclude that their difference is only a factor of order 10. Theorem 3.7 greatly contributes to reducing both computational and space requirements (improvement of order of 10 and 10 2 , respectively). This explains the scalability of our model to large datasets. Finally, the PBLNN clearly decreases the complexity of the model by limiting the number of variables we impose bi-Lipschitzness. See Appendix H.1 for details on the experimental setup.Our bi-Lipschitz model BLNN (and PBLNN) possesses interesting properties. First, it provides by design a theoretical bound of the bi-Lipschitz constants. Some methods cannot necessarily afford this. Furthermore, it enables a direct parameterization of the overall bi-Lipschitzness through the addition of two regularization terms with respective coefficients α and β at different outputs. This translation of bi-Lipschitzness into strong convexity plays an indispensable role as we do not disturb the formulation of the core function, do not track the transformation of the input down to the output and only have two hyper-parameters to monitor, the strict minimum. As a result, we can create bi-Lipschitz functions with known expressive power and more complex than a single layer which others had to compose many times to create a deep bi-Lipschitz neural network. We also avoid the risk of losing tightness in the bounds of the bi-Lipschitz constants caused by this kind of layer-wise parameterization. Therefore, our model satisfies the desired points of Subsection 2.2.While the PBLNN solves several issues of the BLNN, optimization inside the forward pass is still an expensive procedure. Nevertheless, there are plenty of approaches to accelerate the LFT and backpropagation through approximations as well. Since the objective function of the LFT is always strongly concave, its convergence speed does not depend on the dimension, and the only bottleneck is the computation of the gradient information of the objective function (i.e., the ICNN). Approximations by zeroth-order methods can improve the procedure, for example. Moreover, amortizing the LFT as Amos (2023) did may also be a promising simplification of the forward pass. On the other hand, backpropagation through the whole forward pass can also be simplified based on Theorem 3.7. For instance, the Hessian matrix can be replaced by its diagonal or upper triangular elements, making the inversion operation easier. Other improvements, including those for the backward pass, can be adapted from research on DEQs (Bai et al., 2019) .Another limitation, which is overcome by the PBLNN, is that the BLNN cannot represent a function that is not the gradient of a convex function. However, such a type of function is the core of machine learning problems related to optimal transport (Santambrogio, 2015) and some physical systems (Huang et al., 2021) . In the following chapter, we will show that our model can be applied to various problems and its performance is not overshadowed by these limitations. In short, it is a necessary price to pay to gain a bi-Lipschitz model with features such as known expressive power and high controllability so that it can outperform other methods with looser constraint but lower controllability and fewer guarantees. Indeed, if we define the unit of a bi-Lipschitz model as the basic architecture that requires the minimal number of hyperparameters, i.e., one for Lipschitzness and one for inverse-Lipschitzness, most existing models are limited to constructing simple (bi-)Lipschitz units with low expressive power (e.g., only linear) and they have to compose those units to achieve higher expressive power, which leads to looser bounds. In that sense, our model still has a higher expressive Table 1 : Tightness of Lipschitz bound when fitting f (x) = x (x < 0), x + 1 (x ≥ 0) . Mean over five trials. See Table 5 power and tighter bounds than other (bi-)Lipschitz units as ours can with only one unit produce complex functions and its parameterization is not layer wise. Now, if we can afford to sacrifice the tightness of the bounds to further increase expressive power, we can proceed like other methods by stacking multiple BLNNs or combining them with other architectures to suit the characteristics of the problem at hand (see Appendix F).Further substantial extensions of our method can be found in Appendix F, such as generalization to different input and output dimensions, to non-homogeneous functions and to other norms. 
experiments: 
 Experimental setups can be found in Appendix H and codes in https://github.com/yuri-k111/ Bi-Lipschitz-Neural-Network. Further detailed results are summarized in Appendix G.In this subsection, the goal is to empirically verify that (i) our model achieves a tight control of bi-Lipschitz constants, and (ii) it provides new beneficial behaviors different from other methods. We focus on simple problems as they effectively convey key ideas.Tightness of the Bounds Here, we verify the quality of the bi-Lipschitz bounds when the model undergoes training. Especially, we focus on the Lipschitz bound since it is the most affected by the approximation of LFT, and there exist many works to compare with it. Inspired by the experiment of Wang and Manchester (2023) , we aim to learn the function f (x) = x (x < 0), x + 1 (x ≥ 0) that has a discontinuity at x = 0. We take as comparison the LMN (Nolte et al., 2023) , the BiLipNet (Wang et al., 2024) and the i-ResNet network represented by its substructures: spectral normalization (SN) (Miyato et al., 2018) , AOL (Prach and Lampert, 2022) , Orthogonal (Trockman and Kolter, 2021) , SLL (Araujo et al., 2023) , Sandwich (Wang and Manchester, 2023) . The Lipschitzness of each model is constrained by a constant L. A model with a tight Lipschitz bound should achieve that upper bound around x = 0 in order to reproduce the behavior of f . The percentage between the empirical Lipschitz constant and the imposed upper bound L can be found in Table 1 . Interestingly, our method is the only one that achieves an almost perfect score for all settings, while for others the tightness is decreasing. This can be understood as the result of the direct control of bi-Lipschitzness without relying on the information of the individual layers and the construction which uses only direct parameterizations. See Figure 3 for a visualization and Appendix G.2 for more results. Flexibility of the Model If we can underestimate the Lipschitz constant of the target function as the previous experiment, we may also overestimate it. In this case, we observe that previous methods are influenced by the imposed Lipschitz bound presenting high fluctuations in the learned representation or slower learning speed. This is illustrated when learning the identity function y = x with a Lipschitz constraint of L = 1000 in Figure 4 . Our model can learn without any problem and the learning speed is hardly affected by L. This may be due to the way we control the Lipschitz constant of the model. In layer-wise models, the Lipschitz constant is adjusted by scaling the input, while in ours we add a regularization term at the end, resulting in different loss landscapes with seemingly better regularization performance for the latter strategy. See Appendix G.3 for more results.We now illustrate that our model can be efficiently used in areas where bi-Lipschitzness already plays an essential role as an inductive bias. In the estimation of uncertainty based on a single neural network, out-of-distributions and in-distributions points may overlap in the feature space making them indistinguishable. Imposing inverse Lipschitzness on the underlying neural network of the architecture is thus important as it can avoid this phenomenon called feature collapse. Lipschitzness is also beneficial to improve the generalization performance. In one of the state-of-the-art approaches called deterministic uncertainty quantification (DUQ) from (Van Amersfoort et al., 2020) , bi-Lipschitzness is constrained through a regularization method. Therefore, we can replace the gradient penalty with a hard restriction by changing the neural network to a BLNN, resulting in a method we call here DUQ+BLNN. See Appendix E for a precise formulation of this method and theme.The first experiment we lead is with the two moons dataset. Here, we compare DUQ+BLNN with the deep ensembles method (Lakshminarayanan et al., 2017) , DUQ and a DUQ with no bi-Lipschitz regularization (DUQ with no reg.). Results are plotted in Figure 5 . The ideal scenario is for the boundary of the yellow area to be as close as possible to the training data. As we can observe, by adding a bi-Lipschitz constraint, the area of certainty is decreased around the training data. DUQ+BLNN with (α, β) = (2, 4) achieves a tighter area than DUQ. We chose the value of α and β based on a grid search. We took the hyper-parameters with the highest accuracy and α since a higher α is expected to create a tighter yellow area as shown in Figure 5 . Note this strategy only relies on the training data, and such strategy is possible thanks to our direct and simple parameterization of the bi-Lipschitz constants. This clearly shows the advantage of the unique features of our model. See Appendix G.6.1 for further results and discussion on this experiment.Fashion-MNIST Next, we use real world data of FashionMNIST (Xiao et al., 2017) , MNIST (Le-Cun and Cortes, 2010) and NotMNIST (Bulatov, 2011) . The task is to learn to classify FashionMNIST, but at the end of training we verify whether the uncertainty of the model significantly increases when other types of data such as MNIST or NotMNIST is given to the model. This task to distinguish FashionMNIST from MNIST datasets is known to be a complicated task (Van Amersfoort et al., 2020) . We compute the AUROC for the detection performance. The result is shown in Table 2 . Our Sometimes, it happens that we have preliminary knowledge on some type of monotone behaviors of the dataset (Nolte et al., 2023) . For example, COMPAS (Angwin et al., 2016) , BlogFeedBack LoanDefaulter (Nolte et al., 2023) , HeartDisease (Janosi et al., 1988) and AutoMPG (Quinlan, 1993) are benchmark datasets that possess a monotone inductive bias on some variables. See Nolte et al. (2023) for further details on the dataset. As a result, it is more efficient to tune the architecture of the trained model so that it successfully reflects this inductive bias, and various models have been proposed to address this challenge. This is another field where we can apply our architecture that can create monotone (or inverse Lipschitz) functions. We can also control the Lipschitzness to improve the generalization performance. We compare our PBLNN with two state-of-the-art methods: LMN (Nolte et al., 2023) and SMNN (Kim and Lee, 2024) in Table 3 , and with other models in Table 9 as well. As we can observe, our model is competitive with the others.Generalization and Scalability Test Furthermore, we used the dataset provided by Nolte et al. (2023) , CIFAR101, which is a slight augmentation of the original CIFAR100 dataset and designed to exhibit a monotone behavior with respect to one variable. We adopted their training scheme, utilizing the entire dataset for training and intentionally overfitting the data in order to assess both the scalability and expressive power of the model. Successfully, we achieved a 0 loss and 100% accuracy for this experiment, and the convergence was faster than that of Nolte et al. (2023) . This illustrates the high scalability and expressive power of the PBLNN.The following holds for any x, y ∈ domF :(∇F (x) -∇F (y)) ⊤ (x -y) ≥ 1 γ ∥∇F (x) -∇F (y)∥ 2 .Convexity and strong convexity is defined as follows:Definition B.3. F : R m → R is convex if F (tx + (1 -t)y) ≤ tF (x) + (1 -t)F (y)holds for all t ∈ [0, 1] and x, y ∈ domF .Definition B.4. Let µ > 0. F : R m → R is µ-strongly convex if F (x) -µ 2 ∥x∥ 2 is convex.We also introduce the Lengendre-Fenchel transformation which will be a core process in our construction. The first step is to notice that the gradient of a real-valued L-strongly convex function F becomes L-inverse Lipschitz as pointed out by Kinoshita et al. (2023) .Proposition B.6. Let F be an α-strongly convex differentiable function. Then ∇F is α-inverse Lipschitz.Proof. Since F is α-strongly convex,F (y) ≥ F (x) + ∇F (x) ⊤ (y -x) + α 2 ∥y -x∥ 2 .Similarly,F (x) ≥ F (y) + ∇F (y) ⊤ (x -y) + α 2 ∥x -y∥ 2 .By summing both inequalities side by side, we obtainα∥x -y∥ 2 ≤ (∇F (x) -∇F (y)) ⊤ (x -y) ≤ ∥∇F (x) -∇F (y)∥∥y -x∥,where we used Cauchy-Schwarz inequality for the right inequality. As a result,∥∇F (x) -∇F (y)∥ ≥ α∥x -y∥.Therefore, an α-strongly convex neural network can be first built, and then its gradient calculated in order to construct a neural network which is guaranteed to be α-inverse Lipschitz. Now, since the derivative of a smooth function is Lipschitz by definition, we can similarly proceed to construct a function which is both Lipschitz and inverse Lipschitz. That is, we aim to compose a function which is both strongly convex and smooth.Interestingly, smoothness and strong convexity are closely related through the Legendre-Fenchel transformation.Proposition B.7. If F is a closed 1/β-strongly convex function. Then its Legendre-Fenchel transformation is β-smooth.See Zhou (2018) for a proof. Importantly, the smoothness of the Legendre-Fenchel transform does not depend on that of F as long as it is strongly convex. This means that this statement holds also for strongly convex neural networks with ReLU activation functions, and consequently non-differentiable. Proposition B.8. The resulting function F * of a Legendre-Fenchel transformation is also convex as long as its domain is convex.Proof. For all t ∈ [0, 1] and x 1 , x 2 ∈ domF * ,F * (tx 1 + (1 -t)x 2 ) = sup y∈I {⟨y, tx 1 + (1 -t)x 2 ⟩ -F (y)} = sup y∈I {t (⟨y, x 1 ⟩ -F (y)) + (1 -t) (⟨y, x 2 ⟩ -F (y))} ≤t sup y∈I {⟨y, x 1 ⟩ -F (y)} + (1 -t) sup y∈I {⟨y, x 2 ⟩ -F (y)} =tF * (x 1 ) + (1 -t)F * (x 2 ).This leads to the following statement. Theorem B.9. Let F be a closed 1/β-strongly convex function and α ≥ 0. Then the following function is α-strongly convex and α + β-smooth: Proof. Since F is a closed 1/β-strongly convex function, we know that F * is convex and β-smooth. Now clearly F * + α 2 ∥x∥ 2 is α-strongly convex by definition of strong convexity. Moreover, sinceF * (x) = sup y∈I {⟨y, x⟩ -F (y)} + α 2 ∥x∥ 2 .∇ F * (x) + α 2 ∥x∥ 2 -∇ F * (y) + α 2 ∥y∥ 2 = ∥∇F * (x) + αx -(∇F * (y) + αy)∥ ≤ ∥∇F * (x) -∇F * (y)∥ + ∥αx -αy∥ ≤ (β + α)∥x -y∥, F * + α 2 ∥x∥ 2 is α + β-smooth. Corollary B.10. Suppose F * is constructed as Theorem B.9, then ∇ F * is (α, α + β)-bi-Lipschitz, and f * := ∇ F * (x) = argmax y {⟨y, x⟩ -F (y)} + αx.See Zhou (2018) for a proof of ∇F * (x) = argmax y {⟨y, x⟩ -F (y)}. See Figure 6 for a summary of the construction flow. Remark B.11. This construction can only handle functions with the same input and output dimensions. However, this is a common problem in bi-Lipschitz architectures such as that of Behrmann et al. (2019) . It can be addressed by using a composition of bi-Lipschitz functions which will be explained in Appendix F. 
conclusion: 
 We built a model called BLNN based on convex neural networks and the LFT so that bi-Lipschitz functions more complex than a single layer can be constructed and its bi-Lipschitz constants manipulated through the coefficient of two regularization terms added at different outputs. That way, BLNN not only achieves such a tight, direct and simple control but also provides rigorous straightforward analysis on the expressive power and on approximations involved in practice. We illustrated with experiments its distinctive advantageous features compared to prior models. While the primary focus of this paper was to establish a framework suited for solid analyses and for the practical control of bi-Lipschitzness, it is, of course, essential to complement this effort with more varied machine learning applications in future work. We still believe this work on its own contributes to the further exploitation of bi-Lipschitzness in prospective fields and to the deeper understanding of neural networks by delivering a model with unique features for the control of this central inductive bias. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: See Appendix H.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Answer: [Yes] Justification:Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [NA]• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [Yes] Justification: See supplemental material.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
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 to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA] Justification: The paper does not involve crowdsourcing nor research with human subjects Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or 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?Answer: [NA] Justification: The paper does not involve crowdsourcing nor research with human subjects. 
 
 Appendix: In this appendix, we review existing architectures that aim to control the Lipschitzness and inverse Lipschitzness separately as existing bi-Lipschitz models borrow a lot of ideas from prior works on Lipschitzness and inverse Lipschitzness. Moreover, this appendix may constitute a reference for future bi-Lipschitz architectures.Let us start by the Lipschitz property which has been extensively researched over the past few years.A.1 LipschitznessA.1.1 Definition Definition A.1 (Lipschitzness). Let L > 0. f : R l → R t is L-Lipschitz (continuous) if ∥f (x) - f (y)∥ ≤ L∥x -y∥ holds for all x, y ∈ R l .L is called the Lipschitz constant and denoted as Lip(f ).As we can observe, this definition represents the concept of maximal sensitivity with respect to the input as changes in the input upper-bounds changes in the output. This is also referred to as smoothness since when the function is differentiable the Lipschitz constant upper bounds the L 2 -norm of its Jacobian.A convenient property which is often used in Lipschitz constrained neural networks is its closure under composition and addition. Thanks to these properties, the problem of Lipschitz regulation can be decomposed into smaller ones for neural networks constructed by the composition and addition of many simple layers.Due to its omnipresence in the field of machine learning, especially in the certification of robustness against adversarial attacks, Lipschitz constrained neural networks have been the focus of many prior works. Interestingly, their realization greatly varies, delivering different solutions for the same ultimate goal.One of the earliest works trying to control Lipschitzness in the context of machine learning proposes to clip the weights of the function to lie within a range of [-c, c] after back-propagation (Arjovsky et al., 2017) . Consequently, the set of possible parameters would be limited to a compact set, meaning that the overall function is guaranteed to be Lipschitz with a constant dependent on c. However, we cannot precisely know this value, and Miyato et al. (2018) pointed out this kind of brutal restriction favors low-rank weights, leading to a serious decrease of the expressive power of the neural networks.Regularization Another simple approach but quite effective is to add a regularization term to the loss function so that the behavior of the overall function is encouraged to become Lipschitz in a certain way. On the one hand, we can introduce a regularization term which penalizes the whole gradient of the network as Gulrajani et al. (2017) did in the context of Wasserstein GAN (Arjovsky et al., 2017) with the following term:E x [(∥∇F (x)∥ 2 -1) 2 ].(1) On the other hand, we can just focus on the weights. For example, Cisse et al. (2017) proposed a Lipschitz architecture called the Parseval network which uses the following regularization term:∥W ⊤ W -I∥ 2 . As a result, all the weight matrices are incited to become orthogonal, which means 1-Lipschitz.The main downside of these regularization methods is that they can only impose soft constraints and cannot provide any theoretical guarantees on the Lipschitzness of the neural network. This is crucial in some cases such as adversarial attacks where sensitivity has to be exactly controlled. While the implementation is not so complicated as we only need to add a term to the loss function, the computational cost heavily depends on the formulation of the regularization term. Moreover, penalty on the whole gradient like that of equation ( 1) may not be truly efficient since the expectation can only be computed for a limited number of points, which means that it is unclear whether or not the regularization effect will propagate throughout the whole function.The Lipschitz constant of a linear layer is equivalent to its largest singular value since ∥W x -W y∥ ≤ ∥W (x -y)∥ ≤ ∥W ∥∥x -y∥. As a result, the problem amounts to normalizing the spectral norm of the weights. This is one of the most famous approaches in this sector. The estimation of the largest singular value is thus crucial for this technique, and power iteration is often used in practice.Once we evaluated the spectral norm of each layer, we can for example normalize it (Miyato et al., 2018) : W = V /∥V ∥, or rescale only weights whose norm is above a certain threshold λ (Gouk et al., 2021) :W = 1 max(1, ∥V ∥ 2 /λ) V.This can be applied to most types of layers used so far in practice including convolutional layers which are only in essence linear operators (Farnia et al., 2019) .During training, since the weights are constantly changing, it is not realistic to exactly compute the spectral norm at each iteration. Nevertheless, it has been observed that weights evolve more slowly than the convergence of the power-iteration. As a consequence, we can only execute a few steps of the power-iteration at each step and at the end of training arrive at linear layers normalized as desired. This drastically reduces the computational cost.Spectral normalization enables us to simplify the problem of Lipschitz control of the overall function into that of the singular value of the linear units, as the overall Lipschitz constant becomes the product of the spectral norm of each layer. However, it has been repeatedly shown that this approach considerably overestimates the true Lipschitz constant of the whole function (Fazlyab et al., 2019) . For example, consider the composition of the following two matrices:A = 100 0 0 1/100 , B = 1/100 0 0 100 .The spectral norm is determined by the direction of the maximal expansion of the layer but this direction may not be aligned throughout all the layers, especially when there is a nonlinear activation function between each of them. As a result, it engenders a gap between the overall maximal sensitivity and the product of that of each layer, and an intentional control becomes really difficult.Similarly, since we are manipulating layers without taking into account the complex interaction between each other, the magnitude of the gradient may vanish (Anil et al., 2019) . Indeed, when the weights are restricted to be 1-Lipschitz, the norm of the gradient can only decrease through the layers during the backward pass, leading to potential vanishing gradients (Anil et al., 2019) . Notably, it was proved that this kind of 1-Lipschitz neural networks with spectral normalization and with ReLU activation functions cannot represent the absolute value |x| (Anil et al., 2019) which is 1-Lipschitz.In short, layer-wise spectral normalization risks overestimating the general Lipschitz constant, which is problematic if we want to control it, may result in the problem of vanishing gradient and introduce issues into the expressive power of the neural network.Orthogonalization The problem of vanishing gradient can be solved by restricting the matrices to be orthogonal since the eigenvalues of an orthogonal matrix always equal 1. The fact that all matrices are isotropic help to stabilize the training (Prach and Lampert, 2022) . This leads to a large body of work that investigates orthogonal linear layers.There exist several direct realizations of orthogonality. For instance, Xu et al. (2022) ; Huang et al. (2020) uses the form of Singla and Feizi (2021) provides a parameterization with a skew matrix V as follows:W = (V V ⊤ ) -1/2 V.W = exp V -V ⊤ .The latter method needs to approximate the exponential operator with a finite sum. Anil et al. (2019) also construct an orthogonal weight from matrix power series:A k+1 = A k I + 1 2 Q k + . . . + (-1) p -1/2 p Q p k ,whereQ k = I -A ⊤ k A k .Orthogonal convolutional layers are also studied (Li et al., 2019; Wang et al., 2020) , and other realizations were considered based on Cayley transform including convolutions (Trockman and Kolter, 2021; Yu et al., 2022) .The Almost-orthogonal layer (AOL) (Prach and Lampert, 2022) reduces the computational cost of creating orthogonal layers by approximating them by the following formulation:W = V diag   j |V ⊤ V | ij   -1/2 .Empirically, these weights were nearly orthogonal at the end of the training.In order to create L-Lipschitz functions, we can simply multiply the output or input by L. Nevertheless, while orthogonal layers are able to stabilize the training, this regularization destroys the information about the spectrum by setting all the singular values to one. Moreover, the problem of looseness of the Lipschitz bound persists since nonlinear activations are considered apart, and we do not know how they interact with the linear layers.Nonlinear Lipschitz Layer So far, we have reviewed methods that deal with linear layers. Recently, some that incorporate nonlinear functions have started to appear. Meunier et al. (2022) suggest to use residual type layers of the formf (x) = x -2/∥W ∥ 2 W ⊤ σ(W x + b), where W ⊤ σ(W x + b) is the derivative of s(W x + b) and s is a 1-smooth convex function so that s ′ = σ.This type of layer is guaranteed to be 1-Lipschitz. This has been generalized by Araujo et al. (2023) who built a nonlinear layer of the formf (x) = Hx + Gσ(W x + b).If there exists a diagonal matrix Λ with non-negative scalars that satisfiesγI -H ⊤ H + 2αβW ⊤ ΛW -H ⊤ G -(α + β)W ⊤ Λ -G ⊤ H -(α + β)ΛW 2Λ -G ⊤ G ⪰ 0, where σ is slope restricted to [α, β], then f (x) is √ γ-Lipschitz. Especially, when γ = 1, H = 1 and G = -2W T -1, where T is a diagonal with non-negative entries so that T ⪰ W ⊤ W , the corresponding layer is called SDP-based Lipschitz Layer (SLL). A choice of T isT ii = n j=1 |W ⊤ W | ij q j /q i ,where q i > 0. See Theorem 3 of their work for further detail. Fazlyab et al. (2019) proposes LipSDP which estimates the Lipschitz constant of a multi-layer neural network through a semi-definite programming (SDP). It is based on a SDP derived from an integral quadratic constraint. Suppose we reformulate the simple multi-layer neural network asBX = σ(AX + b), f (x) = CX + b L , where σ is [α, β] slope-restricted, X = (x ⊤ 0 , x ⊤ 1 , • • • , x ⊤ L ) ⊤ , A =     W 0 0 • • • 0 0 0 W 1 • • • 0 0 . . . . . . . . . . . . . . . 0 0 • • • W L-1 0     , B =     0 I n1 0 • • • 0 0 0 I n2 • • • 0 . . . . . . . . . . . . . . . 0 0 0 • • • I n L     ,andC = (0, . . . , 0, W L ), b = (b ⊤ 0 , . . . , b ⊤ L-1 ) ⊤ . Now, if there is a diagonal matrix Λ with non-negative entries such that A B ⊤ -2αβΛ (α + β)Λ (α + β)Λ -2Λ A B +     -γI n0 0 • • • 0 0 0 • • • 0 . . . . . . . . . . . . 0 0 • • • W ⊤ L W L     ⪯ 0 is satisfied, then f is √ γ-Lipschitz.Wang and Manchester (2023) provides a direct parameterization for this SDP. Its 1-layer version is called the Sandwich layer and formulated as follows:f (x) = √ 2A ⊤ Ψσ( √ 2Ψ -1 Bx + b),where A and B are produced from the Cayley transformation of an arbitrary matrix with correct dimensions, and Ψ is a diagonal matrix with positive entries.While these methods provide new interesting possibilities to parameterize Lipschitz functions, they are mainly layer-wise, still leading to overestimation of the overall Lipschitz constant, or use the typical structure of the Lipschitzness, meaning that we cannot extend them to bi-Lispchitzness for our purpose. Notably, the SDP-based approach that comes from Fazlyab et al. ( 2019) cannot handle more general structures such as skip connections.A.2.1 Definition Definition A.2 (inverse Lipschitzness). Let L ′ > 0. f : R l → R t is L ′ -inverse Lipschitz if ∥f (x) -f (y)∥ ≥ L ′ ∥x -y∥holds for all x, y ∈ R l . L ′ is called the inverse Lipschitz constant and denoted as invLip(f ).The above property implies that an inverse Lipschitz function is always injective, which means that it has an inverse which is 1/L ′ -Lipschitz.The inverse Lipschitzness is a mathematical description of the minimal sensitivity of the function. By increasing the inverse Lipschitz constant, we can dilute the function and make it more sensitive to small changes of the input. Unfortunately, the inverse Lipschitzness is not closed under addition as we can understand with the simple example of 0 = x -x.To the best of our knowledge, the only model that was built specially for the control of inverse Lipschitzness is that of Kinoshita et al. (2023) . The same model was used earlier by Huang et al. (2021) but for the construction of normalizing flows compatible with Brenier's theorem in optimal transport. Kinoshita et al. (2023) observed that the derivative of an α-strongly convex function is α-inverse Lipschitz. As a result, they propose a model of the typeEQUATIONwhere F is any convex function parameterized by a neural network. This approach has the interesting property that we do not need to know what happens between the layers as we only add a term to the output, which is possible thanks to the convexity of F . That is, a layer-wise control is not required, resulting in a simple and tight control of the inverse Lipschitz constant thanks to the fact that before taking the derivative F can also freely reproduce any convex function that is not necessarily strongly convex. In our work, we start from this model and extend it to bi-Lipschitzness.Finally, a function which is both Lipschitz and inverse-Lipschitz is called bi-Lipschitz. Note that by definition the inverse Lipschitz constant cannot exceed the Lipschitz constant.Definition A.3 (bi-Lipschitzness). Let 0 < L 1 ≤ L 2 . f : R l → R t is (L 1 , L 2 )-bi-Lipschitz if L 1 ∥x -y∥ ≤ ∥f (x) -f (y)∥ ≤ L 2 ∥x -y∥holds for all x, y ∈ R l . L 1 and L 2 are called together bi-Lispchitz constants.Several interpretations can be attributed to this definition. Here, we provide some examples with clearer explanations of the importance of bi-Lipschitzness which were omitted in the main paper. In these applications, both Lipschitzness and inverse Lipschitzness are useful, or sometimes indispensable, and Lipschitzness alone becomes insufficient.1. Injectivity and Out-of-Distribution Detection: Bi-Lipschitz functions are injective thanks to inverse Lipschitzness. This helps distinguish out-of-distribution and in-distribution data in the feature space, making uncertainty estimation possible (Van Amersfoort et al., 2020) . Without inverse Lipschitzness, the problem of feature collapse can occur and compromise the detection of outliers. The inverse Lipschitz constant can control the sensitivity to out-of-distribution points. Please see Appendix E for further details on this topic.2. Quasi-Isometry and Dimensionality Reduction: Bi-Lipschitz functions can be regarded as a quasi-isometry, meaning that the structure of the input is inherited in the output as well. This property is used for adequate dimensionality reduction or embeddings (Li et al., 2020) .3. Invertibility and Solving Inverse Problems: Bi-Lipschitz functions are invertible, and the inverse Lipschitz constant serves as the Lipschitz constant of the inverse function. In that sense, imposing inverse Lipschitzness helps guarantee good properties of the inverse function, just as Lipschitzness assures them for the forward function (Behrmann et al., 2021) . This aspect is used to create normalizing flows (Behrmann et al., 2019) and solve inverse problems (Kruse et al., 2021) .4. Balancing Sensitivity and Robustness: Bi-Lipschitz functions can avoid overly insensitive behavior with respect to their input by controlling the inverse Lipschitz constant. An overly invariant function is also vulnerable to adversarial attacks, as pointed out by Jacobsen et al. (2019) . To the best of our knowledge, the application of bi-Lipschitzness in this direction is underexplored, but this concept may provide an effective solution.In this section, we compare our model BLNN with other models in more details.Convex Potential Flow If we set α = 0, we obtain a BLNN whose output is ∇F * θ , the derivative of the LFT of a 1/β-strongly convex ICNN F θ in terms of Algorithm 1. Since the BLNN is still injective, it has an inverse which is only ∇F θ . As a consequence, in applications where we need to compute the inverse, this model can provide an interesting solution since we can create invertible functions with known inverse. This model was precisely used by Huang et al. (2021) in the context of normalizing flows. It is equivalent to that of Kinoshita et al. (2023) , but the motivation is different. Ours can thus be regarded as an extension of their model as well.Residual Network Interestingly, our model ultimately takes the form of αx + g(x), which can be compared with a residual network. The main differences are that the skip connection is scaled by α and that the formulation of g is restricted to derivatives of convex functions. These two features were crucial components for a direct parameterization of the bi-Lipschitz constants. Behrmann et al. (2019) composed many layers of the form of a residual network to guarantee high expressive power to their invertible residual network. This superposition of layers is only sub-optimal as it leads to a looseness in the bounds of the bi-Lipschitz constants of the overall function and it is not enough to represent some functions such as y = -x. In a sense, our work shows that by restricting ourselves to the derivative of convex functions, this kind of heuristics is not necessary at all and that this condition leads to tighter control.In this appendix, we provide further details and proofs of statements in the main paper.We first remind some definitions.Definition B.1. Let γ > 0. F : R m → R is γ-smooth if F is differentiable and ∥∇F (x) -∇F (y)∥ ≤ γ∥x -y∥ holds for all x, y ∈ R l . That is, ∇F is γ-Lipschitz.In the remainder of this chapter, we will concisely refer to a smooth function when we do not need to specify the smoothness constant. Furthermore, in this work, we will often deal with smooth convex functions. For this type of function, there exist four equivalent characterizations. See Appendix D for the proof.Theorem B.2. Let γ > 0 and F : R m → R a differentiable convex function on a convex domain. Then the following are equivalent:1. F is γ-smooth in the meaning of Definition B.1: ∥∇F (x) -∇F (y)∥ ≤ γ∥x -y∥ .2. The following holds for any x, y ∈ domF :(∇F (x) -∇F (y)) ⊤ (x -y) ≤ γ∥x -y∥ 2 .3. The following holds for any x, y ∈ domF :F (y) ≤ F (x) + ∇F (x) ⊤ (y -x) + γ 2 ∥y -x∥ 2 .In this subsection, we discuss the expressive power of BLNN. It is crucial for this type of construction to clarify this. As we mentioned earlier, layer-wise constraints realize a bi-Lipschitz neural network by restricting the sensitivity of each layer. Nevertheless, this kind of construction is only sub-optimal as it limits the expressive power of the function more than expected. For instance, neural networks with spectral normalization and ReLU activation functions cannot represent y = |x| (Anil et al., 2019) , and an invertible residual network cannot express the linear function y = -x (Zhang et al., 2020) , which were not mentioned in the original paper (Behrmann et al., 2019) . Layer-wise bi-Lipschitz approaches cannot directly inherit the proofs and properties of the original network, which makes the understanding of their behavior more difficult.As for our model, the (α, β)-BLNN, we can prove that 1. Before taking the gradient, our model can approximate any α-strongly convex and α + βsmooth functions on a compact domain endowed with the sup norm (Theorem 3.6). 2. After taking the gradient, we can create a sequence of functions that converges point-wise to any function which is α + β-Lipschitz, α-strongly monotone (see Definition B.12) and the derivative of a real-valued function.The proof is almost straightforward thanks to the fact that our model only makes modification at the output of the core neural network, and we can build on the proofs of previous works (Chen et al., 2019; Huggins et al., 2018) . Definition B.12. A monotone function f : R l → R t is defined as ⟨f Proof. From Proposition 3 of Huggins et al. (2018) , we already know that for any convex function on a compact domain, there is an ICNN with ReLU or softplus-type activation function that approximates it within ϵ in terms of the sup norm. 2 The key point in this proof is to show that the Legendre-Fenchel transformation does not deteriorate the approximation quality and that introducing strongly convexity and smoothness constraints does not cause any unexpected limitations on the expressive power.(x) -f (y), x -y⟩ ≥ 0. Furthermore, f is α-strongly monotone if ⟨f (x) -f (y), x -y⟩ ≥ α∥x -y∥ 2 . B.Let f be a proper closed α-strongly convex and α+ β-smooth function. Since f (x) is α-strongly convex and α + β-smooth, f (x) -α 2 ∥x∥ 2 becomes a convex β-smooth function. Indeed, by definition of strong convexity, f (x) -α 2 ∥x∥ 2 is convex. Moreover, ⟨∇f (x) -αx -(∇f (y) -αy), x -y⟩ ≤ β∥x -y∥ 2 .As a result, by definition of smoothness Rockafellar (1997) ). That is, the Legendre-Fenchel transforms is an involution. Taking f * , we obtain a 1/βstrongly convex function. Similarly, f * (y) -1 2β ∥y∥ 2 is convex. Therefore, there exists an ICNN f so that(Theorem B.2), f (x) -α 2 ∥x∥ 2 is β-smooth. Now, since f (x) := f (x) -α 2 ∥x∥ 2 is a proper closed convex function, f * * = f (Theorem 12.2 ofEQUATION)This precision is preserved after Legendre-Fenchel transformation. Indeed, if we define f * β as the Legendre-Fenchel transformation of f (y) + 1 2β ∥y∥ 2 , then supEQUATIONwhich is equivalent tosup x | f * β (x) -f (x)| < ϵ.This can be derived as follows:f * β (x) = sup y ⟨y, x⟩ -f (y) + 1 2β ∥y∥ 2 = sup y ⟨y, x⟩ -f * (x) + f * (x) -f (y) + 1 2β ∥y∥ 2 ≤ sup y ⟨y, x⟩ -f * (x) + sup y f * (x) -f (y) + 1 2β ∥y∥ 2 ≤ f * * (x) + ϵ,where we used the definition of f * * (x) and equation (4) for the last inequality. Since we can change the role of f and f * and the above inequality holds for all x, we obtain inequality (5). Finally, by adding α∥x∥ 2 /2 to f * β (x), we obtain a function that approximates f within ϵ supx | f * β (x) + α∥x∥ 2 /2 -f (x)| = sup x | f * β (x) -(f (x) -α∥x∥ 2 /2)| = sup x | f * β (x) -f (x)| ≤ϵ.In order words, f is the ICNN that creates a BLNN (before derivation) that approximates f within ϵ.As we can conclude, our model with parameters α and β (before taking the gradient) can approximate any α-strongly convex and α + β-smooth functions on a compact domain endowed with the sup norm. Importantly, the approximation quality of the BLNN equals that of the model we use to create the core convex function.The second point of this section is proved by the following theorem. Theorem B.14 (Huggins et al. (2018) Therefore, an (α, β)-BLNN can represent any α + β-Lipschitz α-strongly monotone (i.e., (α, α + β)bi-Lipschitz) function that is the derivative of another function almost everywhere.In this section, we briefly discuss the difference between the whole set of bi-Lipschitz functions and the expressive power of the BLNN and the possibility of designing an architecture that can approximate any bi-Lipschitz functions based on our model.There are three relevant classes of functions: (1) bi-Lipschitz functions, (2) monotone bi-Lipschitz functions, and (3) cyclically monotone bi-Lipschitz functions. Our model can represent any function of class (3) since cyclically monotone bi-Lipschitz functions are equivalent to the class of derivatives of (strongly) convex functions (Rockafellar, 1970) . In dimensions larger than 1, these three classes are different: f (x, y) = (-y, x) is in (1) but not in (2), and f (x, y) = (x + 2y, y) is in (2) but not in (3). Current bi-Lipschitz models are supposedly in class (1) or (2), like that of Behrmann et al. (2019) and Nolte et al. (2023) . Here, we will only focus on bi-Lipschitz functions with the same input and output dimensions, as most of the bi-Lipschitz models fall into this category. Let us now discuss how class (3) can be used to produce functions of class (1) (and (2)).First, if we suppose the function can be represented by the gradient of another real-valued function, class (1) and (3) are equivalent. Interestingly, under this condition, we can provide an even stronger statement based on Theorem 2.1 of Zlobec (2005) . It says that If f : R n → R is continuously differentiable and its derivative is Lipschitz on a convex set K with some Lipschitz constant L, then there are a convex function C(x) on K and a ≥ L so that f (x) = C(x) -ax ⊤ x/2. In other words, under the condition that a function can be written as the gradient of another real-valued function, class (3) is equivalent to all Lipschitz functions up to a quadratic concave term. Therefore, our model can still express a quite large family of functions.Next, the condition of being the gradient of a function can be interpreted as being rotation-free, i.e., ∇ × F = 0. It is known that any vector field can be decomposed into a rotation-free and divergence-free component (Hodge decomposition theorem) under some regularity conditions. In R 3 , this means that a function f in class (1) can be decomposed as follows:f = ∇ × A + ∇Bwhere A is a vector field, and B a scalar function. From the other direction, we can take∇B as a function of class (3), and A as an arbitrary function with bounded gradient. That way, we can reproduce a large variety of functions so that f is bi-Lipschitz without necessarily being representable as the gradient of another function. The choice of A is still ambiguous so that f is effectively bi-Lipschitz with explicit bi-Lipschitz parameter control, but we believe this approach is promising for the generalization of our method.In conclusion, while our model may have theoretically restrictive expressive power, it still leaves a large liberty to express more general bi-Lipschitz functions both in theory and practice.In this subsection, we provide the explicit formulation of the gradient of the BLNN with respect to the parameters. The first half considers the simple BLNN and the second a more complete case analysis with a complex architecture where BLNN is only one component of it.First of all, the gradient of the loss with respect to the parameters can be more concretely formulated as follows: Lemma B.16. Suppose a loss function L : z → L(z), and the output of the BLNN is f * (x; θ) := ∇ x F * θ (x) + αx as defined in Algorithm 1. Then the gradient can be expressed as follows:EQUATIONwhere z = f * (x; θ). This is only an application of the chain rule. Now, interestingly, ∇ θ F * (x; θ) can be written as a function of the core ICNN F as the following statements shows. Theorem B.17. Suppose F and F * are both differentiable, then the following holds: EQUATION= (∇ ⊤ θ y * (x; θ)) ⊤ x -∂ θ F (y * (x; θ); θ) -(∇ ⊤ θ y * (x; θ)) ⊤ ∇ y F (y * (x; θ); θ) = (∇ ⊤ θ y * (x; θ)) ⊤ x -∂ θ F (y * (x; θ); θ) -(∇ ⊤ θ y * (x; θ)) ⊤ x = -∂ θ F (y * (x; θ); θ),where we used equation (8) for the last equality. This results in the following representation of the gradient (Theorem 3.7). Theorem B.18. Suppose a loss function L : z → L(z), and the output of the BLNN is f * (x; θ) := ∇F * θ (x) + αx as defined in Algorithm 1. If F and F * are both differentiable, then the gradient can be expressed as follows:EQUATIONProof. From Lemma B.16 and Theorem B.17, we know that∇ ⊤ θ L(f * (x; θ)) =∇ ⊤ z L(z) ∇ ⊤ x ∇ θ F * (x; θ) ⊤ = -∇ ⊤ z L(z)∂ ⊤ θ ∇ x F (y * (x; θ); θ).Continuing the procedure, we obtain:∇ ⊤ θ L(f * (x; θ)) = -∇ ⊤ z L(z)∂ ⊤ θ ∇ x F (y * (x; θ); θ) = -∇ ⊤ z L(z)∇ x y * (x; θ)∂ ⊤ θ ∇ y F (y * (x; θ); θ) = -∇ ⊤ z L(z) ∇ 2 y F (y * (x; θ); θ) -1 ∂ ⊤ θ ∇ y F (y * (x; θ); θ), where in the last equality, we used the fact that y * θ = ∇ x F * θ and ∇ y F θ are inverse with respect to each other (Zhou, 2018) . -∇ ⊤ z L(f * (x; θ.requires_grad_(F))) ∇ 2 x F (y * θ (x); θ.requires_grad_(F)) -1 ∇ y F (y * θ (x); θ), where θ.requires_grad_(F) indicates that this θ is considered as a constant and not a parameter to be differentiated with respect to.In practice, a model may possess a more complex architecture where parameters to be optimized do not only come from the ICNN but also from other components. Consider the schema of Figure 7 . Mathematically, this can be formulated as follows. Suppose θ, ϕ and ψ are parameters, f (x; θ) is a BLNN, h(d; ϕ) is a parameterized function like a neural network, and L(z, w; ψ) is a realvalued function returning the loss. In this problem setting, the overall loss can be expressed as L(f (h(d; ϕ); θ), h(d; ϕ); ψ). This means, for a model that incorporates a BLNN, we can classify parameters into three groups:1. those that define the ICNN used in BLNN, 2. those that define the transformation whose output is transferred to the input of the BLNN, 3. and those that are unrelated to the BLNN.For instance, this formulation includes VAEs with a BLNN for the underlying neural network of the decoder. The derivative of the first group of parameters was already discussed above, and that of the last one is straightforward. As for the second family, we can proceed in a similar way as Theorem B.18.Theorem B.20. Suppose a model with loss included L(f * (h(d; ϕ); θ); ψ). The output of the BLNN is f * (x; θ) := ∇F * θ (x) + αx as defined in Algorithm 1. If F and F * are both differentiable, then the gradient of ϕ can be expressed as follows:∇ ϕ L(f * (h(d; ϕ); θ), h(d; ϕ); ψ) =∇ ⊤ z L(z, x; ψ) ∇ 2 y F (y * (x; θ); θ) -1 + αI ∇ ⊤ ϕ h(d; ϕ) + ∇ ⊤ x L(z, x; ψ)∇ ⊤ ϕ h(d; ϕ),where z = f * (h(d; ϕ); θ) and x = h(d; ϕ).Proof. We can proceed similarly as the previous theorems. By chain rule,∇ ϕ L(f * (h(d; ϕ); θ), h(d; ϕ); ψ) = ∇ ⊤ z L(z, w; ψ)∇ ⊤ ϕ f * (h(d; ϕ); θ) + ∇ ⊤ w L(z, w; ψ)∇ ⊤ ϕ h(d; ϕ),where z = f * (h(d; ϕ)) and w = h(d; ϕ). The second term does not require further computation. We will focus on the first term.∇ ⊤ ϕ f * (h(d; ϕ); θ) :=∇ ⊤ ϕ {∇ x F * θ (h(d; ϕ)) + αh(d; ϕ)} =∇ ⊤ ϕ y * (h(d; ϕ); θ) + ∇ ⊤ ϕ αh(d; ϕ) =∇ ⊤ x y * (x; θ)∇ ⊤ ϕ h(d; ϕ) + ∇ ⊤ ϕ αh(d; ϕ) = ∇ 2 y F (y * (x; θ); θ) -1 + αI ∇ ⊤ ϕ h(d; ϕ).This provides the desired result.To summarize, gradients used for the update of each type of parameters can be written as follows:Corollary B.21. Let L := L(f * (h(d; ϕ); θ), h(d; ϕ); ψ). Then under assumptions of Theorem B.20,∇ ⊤ θ L = -∇ ⊤ z L(z, x; ψ) ∇ 2 y F (y * (x; θ); θ) -1 ∂ ⊤ θ ∇ y F (y * (x; θ); θ), ∇ ⊤ ϕ L = ∇ ⊤ z L(z, x; ψ) ∇ 2 y F (y * (x; θ); θ) -1 + αI ∇ ⊤ ϕ h(d; ϕ) + ∇ ⊤ x L(z, x; ψ)∇ ⊤ ϕ h(d; ϕ), ∇ ⊤ ψ L = ∇ ⊤ ψ L(z, x; ψ),where x = h(d; ϕ) and z = f (h(d; ϕ); θ).In this appendix, we discuss the implementation of the LFT as an optimization problem and properties derived from the choice the optimization algorithm. As a reminder, we suppose that a convex solver generates for a fixed x a sequence of points {y t (x)} t=0,...,T based on the objective function ⟨y, x⟩ -F (y), converging to the true optimizer y * (x) := argmax y {⟨y, x⟩ -F (y)}.We first run experiments with the following common convex solvers: steepest gradient descent (GD), Nesterov's accelerated gradient (AGD) (Nesterov, 1983) , Adagrad (Duchi et al., 2011) , RM-Sprop (Hinton et al., 2012) , Adam (Kingma and Ba, 2015) and the Newton method. We calculated the bi-Lipschitz constants of the generated curve y t (x) for each optimization scheme at each iteration t. F was chosen as a two-dimensional ICNN with two hidden layers. The activation was set as the softplus function log(1 + exp(x)). The derivative of the softplus function is the sigmoid function, which means the overall convex neural network becomes smooth. At the output of the ICNN, we added a regularization term ∥x∥ 2 /(2 × 10) so that the overall function becomes 10-strongly convex as well. As a result, the Legendre-Fenchel transformation is also smooth and strongly convex with respective constants c ′ and c that we estimated beforehand. In other words, the Lipschitz and inverse Lipschitz constants of the function y t (x) with respect to x should converge to c ′ and c, respectively. As for the step size, we chose one so that the corresponding convex solver converged. For a µ-strongly convex objective function, it is known that GD converges with a decreasing step size of 1 µ(t+1) . Thus, we chose this for GD but also for Adam and RMSprop since it helped the algorithm to converge. For Adagrad, we set it as 1 µ , for the Newton method as 1 and for AGD as 1/c, the smoothness constant of the objective function. Results are shown in Figures 8 and 9 . See Appendix H for further details on the experimental setup. As we can observe, the simple gradient descent, AGD, Adagrad and the Newton method perform well in the perspective of conservation of bi-Lipschitzness while RMSprop and Adam are not able to provide such feature, especially for the inverse Lipschitz constant. This incites us to use GD for simplicity or Adagrad as convex solver of the Legendre-Fenchel transformation.Now that we have an idea of the behavior of well-known convex solvers, we proceed to the theoretical analysis of the approximation quality of optimization schemes concerning the bi-Lipschitz constants predefined in our model. The aim is to clarify when bi-Lipschitzness is preserved and to evaluate the effective bi-Lipschitz constants through the iterations.Here, we focus on deriving a general non-asymptotic bounds, i.e., the evaluation and evolution of the bi-Lipschitzness of y t (x) for a finite period of timer for general optimization algorithms. This is a crude bound. Theorem C.1. Let the symbols defined as in Algorithm 1. Consider an optimization scheme of sup y {⟨y, x⟩ -F θ (y)} generating points {y t (x)} t that achieves an error of ∥y t (x)-y * (x)∥ ≤ ϵ(t, x) after t iterations, where y * (x) is the global maximum. Then, for all x i , x j such that ∥x i -x j ∥ ≥ δ for δ > 0, α -(ϵ(t, x i ) + ϵ(t, x j ))/δ ≤ ∥f (t) θ (x i ) -f (t) θ (x j )∥ ∥x i -x j ∥ ≤ α + β + (ϵ(t, x i ) + ϵ(t, x j ))/δ,wheref (t) θ (x) := y t (x) + αx is the finite time approximation of f * θ (x).Proof. Concerning the upper bound,∥f (t) θ (x i ) -f (t) θ (x j )∥ =∥y t (x i ) + αx i -(y t (x j ) + αx j )∥ ≤α∥x i -x j ∥ + ∥y t (x i ) -y t (x j )∥ ≤α∥x i -x j ∥ + ∥y t (x i ) -y * (x i )∥ + ∥y * (x i ) -y * (x j )∥ + ∥y * (x j ) -y t (x j )∥ ≤α∥x i -x j ∥ + ϵ(t, x i ) + ∥y * (x i ) -y * (x j )∥ + ϵ(t, x j ) ≤α∥x i -x j ∥ + ϵ(t, x i ) + ∥∇F * θ (x i ) -∇F * θ (x j )∥ + ϵ(t, x j ) ≤α∥x i -x j ∥ + ϵ(t, x i ) + β∥x i -x j ∥ + ϵ(t, x j ) =α∥x i -x j ∥ + β∥x i -x j ∥ + ϵ(t, x i ) + ϵ(t, x j ),where F * θ is the Legendre-Fenchel transformation of F θ . As for the lower bound, we begin by evaluating ⟨y t (x i ) -y t (x j ), x i -x j ⟩. ⟨y t (x i ) -y t (x j ), x i -x j ⟩ =⟨y t (x i ) -y * (x i ), x i -x j ⟩ + ⟨y * (x i ) -y * (x j ), x i -x j ⟩ + ⟨y * (x j ) -y t (x j ),x i -x j ⟩ ≥ -∥y t (x i ) -y * (x i )∥∥x i -x j ∥ + ⟨∇F * θ (x i ) -∇F * θ (x j ), x i -x j ⟩ -∥y * (x j ) -y t (x j )∥∥x i -x j ∥ ≥ -ϵ(t, x i )∥x i -x j ∥ + ⟨∇F * θ (x i ) -∇F * θ (x j ), x i -x j ⟩ -ϵ(t, x j )∥x i -x j ∥ ≥ -ϵ(t, x i )∥x i -x j ∥ -ϵ(t, x j )∥x i -x j ∥,where for the last inequality we used that for any convex function ⟨∇F* θ (x i )-∇F * θ (x j ), x i -x j ⟩ ≥ 0. Now since f (t) θ (x) := y t (x) + αx ⟨f (t) θ (x i ) -f (t) θ (x j ), x i -x j ⟩ =⟨y t (x i ) + αx i -(y t (x j ) + αx j ), x i -x j ⟩ =⟨y t (x i ) -y t (x j ), x i -x j ⟩ + α∥x i -x j ∥ 2 ≥α∥x i -x j ∥ 2 -ϵ(t, x i )∥x i -x j ∥ -ϵ(t, x j )∥x i -x j ∥.By Cauchy-Schwarz inequality, ⟨f(t) θ (x i ) -f (t) θ (x j ), x i -x j ⟩ ≤ ∥f (t) θ (x i ) -f (t)θ (x j )∥∥x i -x j ∥. Thus, we finally obtain∥f (t) θ (x i ) -f (t) θ (x j )∥ ≥ α∥x i -x j ∥ -ϵ(t, x i ) -ϵ(t, x j ).In order to derive the inequalities of the statement, it suffices to divide all sides by ∥x i -x j ∥ and use ∥x i -x j ∥ ≥ δ on the two obtained bounds.Since ϵ(t, x) → 0 as t → ∞ for any optimization scheme that converges, we can recover an arbitrarily close approximation of α and α + β, the inverse Lipschitz and Lipschitz constants, with enough number of iterations. This is a general bound that can be applied to any algorithm with known behavior. Below, we provide two concrete evaluations applied to the gradient descent.Corollary C.2. Suppose F θ is µ-strongly convex and x -∥∂ sub F θ (y t )∥ ≤ G(x). If we employ gradient descent with η t = 1/(µ(t + 1)) as a step size and a fixed initial point y 0 = 0, ϵ(t, x) = 2G(x)/(µ √ t). Therefore,α -2(G(x i ) + G(x j ))/(δµ √ t) ≤ ∥f (t) θ (x i ) -f (t) θ (x j )∥ ∥x i -x j ∥ ≤ α + β + 2(G(x i ) + G(x j ))/(δµ √ t).Proof. This follows from Lemma 1 of Rakhlin et al. (2012) .A similar bound can be provided when x -∥∂ sub F θ (y t )∥ ≤ G(x) is satisfied almost surely. See Lemma 2 of Rakhlin et al. (2012) .Corollary C.3. Suppose F θ is µ-strongly convex and γ-smooth. If we employ gradient descent with η t = 1/γ as a step size and a fixed initial point y 0 , we can chooseϵ(t, x) = 1 -µ 2 γ 2 t/2 ∥y 0 (x) - y * (x)∥ 2 .Proof. Since F is γ-smooth, h(y) = ⟨y, x⟩ -F (y) is γ-smooth and µ-strongly concave. As a result,∥y t+1 -y * ∥ 2 =∥y t + η∇h(y t ) -y * ∥ 2 =∥y t -y * ∥ 2 -2η⟨∇h(y t ), y t -y * ⟩ + η 2 ∥∇h(y)∥ 2 =∥y t -y * ∥ 2 -2η⟨∇h(y t ) -∇h(y * ), y t -y * ⟩ + η 2 ∥∇h(y t ) -∇h(y * )∥ 2 ≤∥y t -y * ∥ 2 -2η 1 γ ∥∇h(y t ) -∇h(y * )∥ 2 + η 2 ∥∇h(y t ) -∇h(y * )∥ 2 =∥y t -y * ∥ 2 - 1 γ 2 ∥∇h(y t ) -∇h(y * )∥ 2 ≤∥y t -y * ∥ 2 - µ 2 γ 2 ∥y t -y * ∥ 2 = 1 - µ 2 γ 2 ∥y t -y * ∥ 2 = 1 - µ 2 γ 2 t+1 ∥y 0 -y * ∥ 2 ,where we used that ∇h(y * ) = 0 for the third equality, the co-coercivity of smooth convex functions for the first inequality, and inverse Lipschitzness for the second inequality. This is a crude estimation of the bi-Lipschitz constants. As long as we assume that all values can only take discrete ones with intervals bigger than δ, then this theorem provides a meaningful evaluation. For example, during training this does not pose a problem since the number of training data is usually finite. Technically speaking, we will not have any problem in practice either as the computer also deals with discrete values.The limitation of the above theorem appears once we authorize infinite precision as we could take two values arbitrarily close to each other, i.e., δ → 0, leading the lower and upper bound to explode. This makes the model prone to adversarial attacks. One way to remedy this is to evaluate the Legendre-Fenchel transformation for some discrete values and interpolate with some bi-Lipschitz functions (e.g., linear) so that the whole function satisfies the expected bounds.Ideally, we would like to derive a bound of the typeEQUATIONwhere h 1 and h 2 converge as ∥x i -x j ∥ → 0 and converge to 0 as t → ∞. We could not achieve this in Theorem C.1, which is predictable as we did not take into account the similarity of the optimization paths of x i and x j and derived the bounds by isolating each point. The next step is thus to take into account this closeness. Nevertheless, this kind of analysis is more involved and we do not know if it is applicable to all optimization schemes. This could explain the existence of some algorithms where the bi-Lipschitzness do not evolve as desired as we have pointed out in Figure 9 .The following theorem shows the first half of the statement of Theorem 3.5. That is, we show that in the limit there is no bias for GD in terms of bi-Lipschitz constants.Theorem C.4. Let the symbols defined as in Algorithm 1. Consider an optimization scheme of sup y {⟨y, x⟩ -F (y)} generating points {y t (x)} t at the t-th iterations. If F is µ-strongly convex, γ-smooth and twice differentiable and we employ gradient descent with a step size η t so that the discrete optimization converges to the global maximum y * (x), then for all x1 γ ⪯ ∇ ⊤ x y ∞ (x) ⪯ 1 µ .Proof. The recurrence relation of GD is as follows:y t+1 (x) = y t (x) + η t {x -∇ y F (y t (x))} .Since F is twice differentiable, we can take the Jacobian of both sides:∇ ⊤ x y t+1 (x) = ∇ ⊤ x y t (x) + η t I -∇ 2 y F (y t (x))∇ ⊤ x y t (x) .By replacing y t+1 (x) and y t (x) by y ∞ (x) as the sequence is converging by hypothesis, we obtain∇ ⊤ x y ∞ (x) =∇ ⊤ x y ∞ (x) + η t I -∇ 2 y F (y ∞ (x))∇ ⊤ x y ∞ (x) = I -η t ∇ 2 y F (y ∞ (x) ∇ ⊤ x y ∞ (x) + η t I Therefore, since∇ 2 y F (y) is invertible as µ ⪯ ∇ 2 y F (y), ∇ ⊤ x y ∞ (x) = ∇ 2 y F (y ∞ (x)) -1 .By using µ ⪯ ∇ 2 y F (y) ⪯ γ, we arrive at the desired result.Remark C.5. Note that y ∞ (x) and y * (x) are different in the sense that the former is defined by the limit (t → ∞) of the recurrence relation of the optimizer, while the latter is defined as the global maximum of sup y {⟨y, x⟩ -F (y)} which also satisfies ∇ y F (y * (x)) = x. In short, they differ in the formulation and explicit dependence on x.Since F is µ-strongly convex and γ-smooth, its Legendre-Fenchel transform is 1/γ-strongly convex and 1/µ-smooth, leading to a (1/γ, 1/µ)-bi-Lipschitz function by taking the derivative of the latter. This behavior is inherited by the optimization scheme in the limit of t → ∞ as the above theorem shows. The point is that there is no bias, which may have occurred if the influence of the step size persisted in the final result. A similar result can be obtained for AGD. For the others, the analysis is more complicated as the gradient is accumulated throughout the iterations. Since our main interest is non-asymptotic behavior, we will not further develop this.We will now prove the second half of the theorem and show that equation ( 10) can be actually established for some parameter settings at least for Lipschitzness.Theorem C.6. Let the symbols defined as in Algorithm 1. Consider an optimization scheme of sup y {⟨y, x⟩ -F (y)} generating points {y t (x)} t at the t-th iterations and y * (x) is the global maximum. If F is µ-strongly convex and γ-smooth and we employ gradient descent with η t = 1/µ(t + 1) as a step size and y 0 (x i ) = y 0 as initial point, then for all x i , x j∥y t+1 (x i ) -y t+1 (x j )∥ ≤ 1 -2η t µ + η 2 t γ 2 ∥y t (x i ) -y t (x j )∥ + η t ∥x i -x j ∥.As a result,∥y t+1 (x i ) -y t+1 (x j )∥ ≤ h(t)∥x i -x j ∥,wherelim t→∞ h(t) = 1 µ .Proof. Let us first prove the recurrence relation of the statement. Throughout the proof, we use the following notations:∆y t := y t (x i ) -y t (x j ), ∆x := x i -x j , ∆D t := ∇F (y t (x i )) -∇F (y t (x j )).Since ∥∆y t+1 ∥ = ∥∆y t -η t ∆D t + η t ∆x∥ ≤ ∥∆y t -η t ∆D t ∥ + η t ∥∆x∥, (11) it suffices to evaluate ∥∆y t -η t ∆D t ∥. Now,∥∆y t -η t ∆D t ∥ 2 =∥∆y t ∥ 2 -2η t ⟨∆y t , ∆D t ⟩ + η 2 t ∥∆D t ∥ 2 ≤∥∆y t ∥ 2 -2η t µ∥∆y t ∥ 2 + η 2 t ∥∆D t ∥ 2 ≤∥∆y t ∥ 2 -2η t µ∥∆y t ∥ 2 + η 2 t γ 2 ∥∆y t ∥ 2 = 1 -2η t µ + η 2 γ 2 ∥∆y t ∥ 2 ,where for the first inequality we used that F is µ-strongly convex, i.e., ⟨∇F (y) -∇F (y ′ ), y -y ′ ⟩ ≥ µ∥y -y ′ ∥ 2 and for the second inequality we used that F is γ-smooth.As µ ≤ γ always holds,1 -2η t µ + η 2 γ 2 = 1 -2η t µ + η 2 µ 2 + η 2 (γ 2 -µ 2 ) = (1 -η t µ) 2 + η 2 (γ 2 -µ 2 ) ≥ 0. As a result ∥∆y t -η t ∆D t ∥ ≤ 1 -2η t µ + η 2 γ 2 ∥∆y t ∥.By substituting this to equation ( 11), we obtain ∥∆y t+1 ∥ ≤ 1 -2η t µ + η 2 t γ 2 ∥∆y t ∥ + η t ∥∆x∥, which was the desired inequality.Next, we assume ∥∆x∥ > 0 and divide both sides by ∥∆x∥. Denoting α t := µ∥∆y t+1 ∥/∥∆x∥, we arrive at α t+1 ≤ 1 -2η t µ + η 2 t γ 2 α t + η t µ. We want to show that lim t→∞ α t = 1.By the above-mentioned argument when t ≥ 1, we can further develop this as follows:α t+1 ≤ 1 -2η t µ + η 2 t γ 2 α t + η t µ = (1 -η t µ) 2 + η 2 t (γ 2 -µ 2 )α t + η t µ ≤ (1 -η t µ) + 1 2 η 2 t (γ 2 -µ 2 ) 1 -η t µ 2 α t + η t µ = (1 -η t µ) + 1 2 η 2 t (γ 2 -µ 2 ) 1 -η t µ α t + η t µ.Since η t µ = 1/(t + 1), we getα t+1 = t t + 1 + γ 2 /µ 2 -1 2 1 t(t + 1) α t + 1 t + 1.By multiplying both sides by t + 1, defining A t := tα t and κ := γ 2 /µ 2 -1(≥ 0), we can conclude thatEQUATION)Let us prove that from equation 12, we have 13) for all t = 1, 2, . . .. We will proceed by mathematical induction.EQUATIONWhen t = 1, A 1 = α 1 . Since ∥∆y 1 ∥ = ∥∆y 0 -η 0 ∆D 0 + η 0 ∆x∥ and the initial point was a constant independent of x i , ∥∆y 1 ∥ = η 0 ∥∆x∥ = 1 µ ∥∆x∥. Thus, A 1 = α 1 = 1. This means, equation ( 13) is satisfied. Now suppose equation ( 13) holds for t = m. Then by equation 12,A m+1 ≤ 1 + κ 2 1 m 2 A m + 1 ≤ 1 + κ 2 1 m 2   m + κ 2 m-1 k=1 1 k k<i≤m-1 exp κ 2 1 i 2   + 1 ≤m + κ 2 1 m + 1 + κ 2 1 m 2   κ 2 m-1 k=1 1 k k<i≤m-1 exp κ 2 1 i 2   + 1 ≤m + 1 + κ 2 1 m + exp κ 2 1 m 2   κ 2 m-1 k=1 1 k k<i≤m-1 exp κ 2 1 i 2   ≤m + 1 + κ 2 1 m +   κ 2 m-1 k=1 1 k k<i≤m exp κ 2 1 i 2   ≤m + 1 +   κ 2 m k=1 1 k k<i≤m exp κ 2 1 i 2   .Therefore, by mathematical induction equation ( 13) holds for all t = 1, 2, . . ..j i=1 1 i 2 ≤ ∞ i=1 1 i 2 = π 2 6 ∀j = 1, 2, . . . and j i=1 1 i ≤ O(log j),we can conclude thatα t = A t t ≤ t + κ 2 t-1 k=1 1 k k<i≤t-1 exp κ 2 1 i 2 t = t + κ 2 t-1 k=1 1 k exp k<i≤t-1 κ 2 1 i 2 t ≤ t + κ 2 t-1 k=1 1 k exp κ 2 π 2 6 t ≤ t + O(log t) t →1 (t → ∞).This theorem proves that the generated curve y t (x) at a fixed time t is Lipschitz with a constant that converges to the true value with a convergence speed of O(log t/t). This theorem is interesting as it not only guarantees that y t (x) is bi-Lipschitz through the whole iteration and converges to the desired values but also is equipped with a non-asymptotic bound offering a concrete convergence speed and a value at each iteration. The advantage of this proof is that setting the step size as η t = 1/(µ(t + 1)) is realistic, since in our setting we know the strong convexity constant which corresponds to 1/β. In practice, the convergence of the bi-Lipschitz constants seems to be faster, which means that there may exist tighter bounds but we let this investigation for future work since our goal was to just assure a convergence with good properties.As for the lower bound, we could not prove a similar property for GD. Nevertheless, we can guarantee that it will not diverge to -∞ when ∥x i -x j ∥ → 0 as ⟨y t (x i ) -y t (x j ), x i -x j ⟩ ≤ -∥y t (x i ) -y t (x j )∥∥x i -x j ∥ and ∥y t (x i ) -y t (x j )∥ was just proved to be upper-bounded by a constant. We let further theoretical investigation for future work. Since GD shows and assures rather good performance, we will mainly use this algorithm in this paper.An analogue of Theorem 3.5 can be shown with a different step size, i.e., η = 1/γ, where γ is the smoothness constant of the objective function. This step size is optimal as it assures the fastest convergence. However, in practice, it is actually useless as estimating the smoothness constant γ for the step size is computationally demanding and during training this constant is constantly changing. Only after training, it suffices to estimate γ once, and all forward passes can be executed with this step size and convergence speed. The theorem is as follows: Theorem C.7. Let the symbols defined as in Algorithm 1. Consider an optimization scheme of sup y {⟨y, x⟩ -F (y)} generating points {y t (x)} t at the t-th iterations and y * (x) is the global maximum. If F is µ-strongly convex and γ-smooth and we employ gradient descent with η t = 1/γ as a step size and y 0 (x i ) = y 0 as initial point, then for all x i , x j∥y t+1 (x i ) -y t+1 (x j )∥ ≤ 1 - µ γ ∥y t (x i ) -y t (x j )∥ + η t ∥x i -x j ∥.As a result, ∥y t+1 (x i ) -y t+1 (x j )∥ ≤ h(t)∥x i -x j ∥, where lim t→∞ h(t) = 1 µ .Proof. Let us first prove the following inequality:EQUATION)If µ = γ, then by strong convexity of F ,⟨∇F (y) -∇F (y ′ ), y -y ′ ⟩ ≥ µ∥y -y ′ ∥ 2 = µγ µ + γ ∥y -y ′ ∥ 2 + µ 2 µ + γ ∥y -y ′ ∥ 2 .However, since µ = γ, we have µ∥y -y ′ ∥ = ∥∇F (y) -∇F (y ′ )∥. Thus, the above inequality leads to⟨∇F (y) -∇F (y ′ ), y -y ′ ⟩ ≥ µ∥y -y ′ ∥ 2 = µγ µ + γ ∥y -y ′ ∥ 2 + 1 µ + γ ∥∇F (y) -∇F (y ′ )∥ 2 .When γ > µ, since F is µ-strongly convex and γ-smooth, F (y) -µ 2 ∥y∥ 2 is convex and γ -µ smooth. This is straightforward from Definition 2 of Theorem B.2. Now, by applying Definition 4 of Theorem B.2 to F (y) -µ 2 ∥y∥ 2 , we obtain⟨∇F (y) -µy -(∇F (y ′ ) -µy ′ ) , y -y ′ ⟩ ≥ 1 γ -µ ∥∇F (y) -µy -(∇F (y ′ ) -µy ′ ) ∥ 2 .By developing both sides, we arrive at⟨∇F (y) -∇F (y ′ ), y -y ′ ⟩ -µ∥y -y ′ ∥ 2 ≥ 1 γ -µ ∥∇F (y) -∇F (y ′ )∥ 2 -2µ⟨∇F (y) -∇F (y ′ ), y -y ′ ⟩ + µ 2 ∥y -y ′ ∥ 2 .Rearranging the terms, we haveγ + µ γ -µ ⟨∇F (y) -∇F (y ′ ), y -y ′ ⟩ ≥ µγ γ -µ ∥y -y ′ ∥ 2 + 1 γ -µ ∥∇F (y) -∇F (y ′ )∥ 2 ,which leads to⟨∇F (y) -∇F (y ′ ), y -y ′ ⟩ ≥ µγ µ + γ ∥y -y ′ ∥ 2 + 1 µ + γ ∥∇F (y) -∇F (y ′ )∥ 2 .Let us now prove the recurrence relation of the statement. Throughout the proof, we use the following notations:∆y t := y t (x i ) -y t (x j ), ∆x := x i -x j , ∆D t := ∇F (y t (x i )) -∇F (y t (x j )).SinceEQUATIONit suffices to evaluate ∥∆y t -η t ∆D t ∥. Now,∥∆y t -η t ∆D t ∥ 2 =∥∆y t ∥ 2 -2η t ⟨∆y t , ∆D t ⟩ + η 2 t ∥∆D t ∥ 2 ≤∥∆y t ∥ 2 -2η t µγ µ + γ ∥∆y t ∥ 2 + 1 µ + γ ∥∆D t ∥ 2 + η 2 t ∥∆D t ∥ 2 = 1 - 2η t µγ µ + γ ∥∆y t ∥ 2 + η t η t - 2 µ + γ ∥∆D t ∥ 2 .where we used inequality ( 14) for the inequality.Since we set η t = 1/γ,∥∆y t -η t ∆D t ∥ 2 ≤ 1 - 2µ µ + γ ∥∆y t ∥ 2 + 1 γ 1 γ - 2 µ + γ ∥∆D t ∥ 2 ≤ 1 - 2µ µ + γ ∥∆y t ∥ 2 + 1 γ 1 γ - 2 µ + γ µ 2 ∥∆y t ∥ 2 = 1 - µ γ 2 ∥∆y t ∥ 2 ,where for the second inequality we used that 1 γ -2 µ+γ ≤ 0 (µ ≤ γ) and ∥∇F (y) -∇F (y ′ )∥ ≥ µ∥y -y ′ ∥ (inverse Lipschitzness). Substituting this to equation ( 15), we conclude that∥∆y t+1 ∥ ≤ 1 - µ γ ∥∆y t ∥ + η t ∥∆x∥ as µ/γ ≤ 1.Since ∥∆y 1 ∥ = ∥∆y 0 -η 0 ∆D 0 + η 0 ∆x∥ and the initial point was a constant independent of x i , ∥∆y 1 ∥ = η 0 ∥∆x∥ = 1 γ ∥∆x∥, and∥∆y t ∥ ≤ 1 γ t-1 i=0 1 - µ γ i ∥∆x∥.Finally, See Appendix H for further details on the experimental setup. We also leave this research of better approaches for future work.lim t→∞ ∥∆y t ∥ ≤ 1 γ ∞ i=0 1 - µ γ i ∥∆x∥ = 1/γ 1 -(1 -(µ/γ)) = 1 µ ∥∆x∥.In this Appendix, we prove Theorem B.2 by showing a more general theorem extended to dual norms.Dual norms are defined as follows: Definition D.1. The dual norm of a norm ∥ • ∥ is defined as∥y∥ * = sup {⟨x, y⟩ | ∥x∥ = 1} .Let us now prove the following theorem. When the norm in question is L 2 , its dual is also the same and we arrive at the initial statement. Theorem D.2. Let γ > 0, ∥ • ∥ and ∥ • ∥ * a pair of dual norms, and F : R l → R t a differentiable convex function on a convex domain. Then the following are equivalent:D1. The following holds for any x, y ∈ domF :EQUATIOND2. The following holds for any x, y ∈ domF :⟨∇F (x) -∇F (y), x -y⟩ ≤ γ∥x -y∥ 2 .D3. The following holds for any x, y ∈ domF :F (y) ≤ F (x) + ∇F (x) ⊤ (y -x) + γ 2 ∥y -x∥ 2 .The following holds for any x, y ∈ domF :(∇F (x) -∇F (y)) ⊤ (x -y) ≥ 1 γ ∥∇F (x) -∇F (y)∥ 2 * .Proof. D1.⇒ D2.When x = y, the inequality trivially holds. When x ̸ = y, this is straightforward by the Cauchy-Schwarz inequality, since⟨∇F (x) -∇F (y), x -y⟩ =∥x -y∥⟨∇F (x) -∇F (y), x -y ∥x -y∥ ⟩ ≤∥x -y∥ sup {⟨z, ∇F (x) -∇F (y)⟩ | ∥z∥ = 1} =∥x -y∥∥∇F (x) -∇F (y)∥ * ≤γ∥x -y∥ 2 ,where we used D1. in the last inequality. D2.⇒ D3.Consider the function G(t) := F (x + t(y -x)), which is well-defined as the domain of F is convex. Since F is differentiable andF (y) = G(1), F (y) =G(0) + 1 0 d dt G(t)dt =G(0) + 1 0 ⟨∇F (x + t(y -x)), y -x⟩dt =F (x) + ⟨∇F (x), y -x⟩ + 1 0 ⟨∇F (x + t(y x)) -∇F (x), y -x⟩dt ≤F (x) + ⟨∇F (x), y -x⟩ + 1 0 tγ∥x -y∥ 2 dt =F (x) + ⟨∇F (x), y -x⟩ + γ 2 ∥x -y∥ 2 ,where we used D2. for the inequality. D3.⇒ D4.Consider the following function:F x (z) := F (z) -⟨∇F (x), z⟩,which is a γ-smooth function. Moreover, since F is convex, x is the minimizer of F x . Consequently, the following PL inequality generalized to dual norms holds:F x (x) = inf w F x (w) ≤ inf w F x (z) + ⟨∇F x (z), w -z⟩ + γ 2 ∥w -z∥ 2 = inf ∥v∥=1 inf t F x (z) + ⟨∇F x (z), tv⟩ + γ 2 ∥tv∥ 2 = inf ∥v∥=1 inf t F x (z) + t⟨∇F x (z), v⟩ + γ 2 t 2 = inf ∥v∥=1 F x (z) - 1 2γ (⟨∇F x (z), v⟩) 2 =F x (z) - 1 2γ ∥∇F x (z)∥ 2 * ,where we used D3. for the inequality. Since this holds for all z, by substituting z = y, we obtain1 2γ ∥∇F (y) -∇F (x)∥ 2 * = 1 2γ ∥∇F x (y)∥ 2 * ≤ F x (y) -F x (x) = F (y) -F (x) -⟨∇F (x), y -x⟩.Since we can interchange x and y, we also have1 2γ ∥∇F (y) -∇F (x)∥ 2 * ≤ F (x) -F (y) -⟨∇F (y), x -y⟩.Summing side by side these two inequalities, we arrive at1 γ ∥∇F (y) -∇F (x)∥ 2 * ≤ ⟨∇F (x) -∇F (y), x -y⟩.D4.⇒ D1.When x = y, the inequality trivially holds. When x ̸ = y, by Cauchy-Schwarz inequality,1 γ ∥∇F (y) -∇F (x)∥ 2 * ≤ ⟨∇F (x) -∇F (y), x -y⟩ ≤ ∥∇F (x) -∇F (y)∥ * ∥x -y∥.Deep neural networks are nowadays used in many applications such as self-driving cars and large language models. However, in the real world, they are constantly subjected to a large amount of ambiguous scenarios, and assuring fail-safes is a top priority. One approach to realize this is to design the agent so that it quantifies the reliability of its own outputs and behavior, while achieving high performance. However, it is well-known that usual deep neural networks over-confidently extrapolate to unknown data or poorly perform in uncertainty tasks. As a result, the quantification of uncertainty has become the subject of many research, leading to Bayesian neural networks (Neal, 2012) , Monte Carlo dropouts (Gal and Ghahramani, 2016) and deep ensembles (Lakshminarayanan et al., 2017) . Unfortunately, these methods quickly become computationally expensive as they require to process multiple forward passes or to retain distribution samples over a large number of parameters. Furthermore, such ensemble methods cannot in the essence avoid that all constitutive members make the same mistake or place high confidence on the same out-of-distribution data.Therefore, a line of work that concentrates on the uncertainty estimation using a single neural network has recently started to draw attention (Liu et al., 2020a; Van Amersfoort et al., 2020) . In these works, they are interested in quantifying the uncertainty through a single forward pass and analyzing its behavior. On their own, they can not only provide new methods with low computational cost but also address the challenges of ensemble methods as successful individual agents can again be integrated in those models.Thanks to these works, some properties necessary for an accurate uncertainty quantification have been discovered. Notably, bi-Lipschitzness is an indispensable inductive bias. Intuitively, bi-Lipschitzness guarantees the neural network to be distance aware, i.e., distance is moderately preserved between the input and the feature space, resulting in correct detection of out-of-distribution points. Without this property, it is known that the problem of feature collapse occurs, which refers to the phenomenon that out-of-distribution and in-distribution points overlap in the feature space, making out-of-distribution detection impossible in principle (Van Amersfoort et al., 2020).We will focus on a recent model proposed by Van Amersfoort et al. (2020) as it is a model that has been shown to perform as well as deep ensemble methods that are the state-of-the-art in this area. It also performs well on the FashionMNIST dataset.Their model is called Deep Uncertainty Quantification (DUQ) and is mainly used for classification tasks. The idea is to create a class c represented by a centroid e c and to calculate the distance K c between a data point and all centroids. The data will be classified as the label of the closest centroid to that point. They draw inspiration from the radial basis function kernel. The mathematical formulation of K c is as follows:K c (f θ (x), e c ) = exp - 1 n ∥W c f θ (x) -e c ∥ 2 2 2σ 2 ,where f θ is a general neural network also called feature extractor, W c is a weight matrix, and σ is a hyper parameter called length scale.The loss function is the sum of the binary cross entropy over all classes. During the training, the centroids are updated by the following rule:N c,t =γN c,t-1 + (1 -γ)n c,t , m c,t =γm c,t-1 + (1 -γ)σ i W c f θ (x c,t,i ), e c,t =m c,t /N c,t ,where n c,t is the number of data points assigned to class c, x c,t,i is the i-th element of the minibatch corresponding to class c, and γ is a hyper parameter (van den Oord et al., 2017) .Bi-Lipschitzness was incorporated through a gradient penalty term as follows:λ ∥∇ x c K c ∥ 2 -1 .The discussion about this regularization can be found in Section 2. By deleting this regularization term, and replacing the neural network f θ by our bi-Lipschitz neural network, we obtain a model called here DUQ+BLNN.With only one BLNN, we can only provide an output that has the same dimension as the input. This is a common problem in bi-Lipschitz models. In this extension, we provide an architecture that extends to a more realistic case with different input and output dimensions. While we will lose many theoretical guarantees we have provided so far, some will be retained. Furthermore, throughout this paper, we will show that this method works well in practice.The idea follows Wang et al. (2021) and Kinoshita et al. (2023) . Let f * 1 (•; θ 1 ) and f * 2 (•; θ 2 ) be two BLNN with input (and output) dimension d 1 and d 2 , respectively. By composing them asEQUATIONwhere D is a d 2 × d 1 matrix with all diagonal components to 1, we obtain a function with input dimension d 1 and output dimension d 2 . Note that if d 1 ≤ d 2 then the whole bi-Lipschitz property is preserved.On the other hand, if d 1 ≥ d 2 then the Lipschitzness remains at least. Nevertheless, this parameterization ( 16) can represent any bi-Lipschitz function l : R d1 → R d2 with d 1 ≥ d 2 . It suffices to set f * 2 as the identity map and consider the extension l : R d1 → R d1 whose image coincides with that of l in the R d2 -dimensional subspace. This remains bi-Lipschitz like l. Therefore, this representation is reasonable to some extent.The gradient computation can also be explicitly formulated as follows: Theorem F.1. Suppose a model with loss includedL := L(f * 2 (Df * 1 (h(d; ϕ); θ 1 ); θ 2 ), h(d; ϕ); ψ), where f * 1 := ∇F * 1 (x 1 ; θ 1 ) + α 1 x 1 and f * 2 := ∇F * 2 (x 2 ; θ 2 ) + α 2 x 2 are two C 2 BLNN defined in Algorithm 1 with input dimensions d 1 and d 2 , respectively. D is a d 2 × d 1 diagonalmatrix with all diagonal elements set to 1. If F and F * are both differentiable, then the gradients with respect to each parameter can be expressed as follows:∇ ⊤ θ1 L = -∇ ⊤ z LΞ 2 D ∇ 2 y F 1 (y * 1 (x; θ 1 ); θ 1 ) -1 ∂ ⊤ θ1 ∇ y F 1 (y * 1 (x; θ 1 ); θ 1 ), ∇ ⊤ θ2 L = -∇ ⊤ z L(z, x; ψ) ∇ 2 y F 2 (y * 2 (w; θ 2 ); θ 2 ) -1 ∂ ⊤ θ2 ∇ y F 2 (y * 2 (w; θ 2 ); θ 2 ), ∇ ϕ L =∇ ⊤ z L(z, x; ψ)Ξ 2 DΞ 1 ∇ ϕ h(d; ϕ) + ∇ ⊤ x L(z, x; ψ)∇ ⊤ ϕ h(d; ϕ), ∇ ⊤ ψ L =∇ ⊤ ψ L(z, x; ψ),whereΞ 1 := ∇ 2 y F 1 (y * 1 (x; θ 1 ); θ 1 ) -1 + α 1 I, Ξ 2 := ∇ 2 y F 2 (y * 2 (w; θ 2 ); θ 2 ) -1 + α 2 I, z := f * 2 (Df * 1 (h(d; ϕ); θ 1 ); θ 2 ), w := Df * 1 (h(d; ϕ); θ 1 ) and x := h(d; ϕ).As pointed out by Jin and Lavaei (2020) , in some cases, we may have knowledge of specific physical constraints or requirements to design a non-homogeneous control of the sensitivity like in physics-informed neural networks (Raissi et al., 2019) . However, the current model imposes the same bi-Lipschitz constraint for each dimension. In the linear case, this corresponds to bounding all the eigenvalues of the matrix with the same constants from above and below respectively. For example, consider the following function:f (x 1 , x 2 ) = (2x 1 , 100x 2 ).If we can only choose 1 parameter for Lipschitzness, then it would be 100. However, this ignores the non-homogeneity of each dimension and their sparseness. In the above example, the Lipschitzness with respect to x 1 is only 2. Therefore, we would like to be able to control the bi-Lipschitzness of the function with respect to each dimension of the input.This can be executed by introducing a diagonal matrix A and B instead of α and β in our model. That is, instead of adding α∥x∥ 2 /2 and ∥y∥ 2 /(2β) (see Algorithm 1), we add x ⊤ A 2 x/2 and y ⊤ B -2 y/2, respectively, which makes possible a more flexible control of the sensitivity of the overall function. The drawback is that we obtain more parameters to tune, but this can be avoided by designing these matrices as parameters to learn throughout the training too.In the same idea, it may also be interesting to impose the convexity requirement on a limited number of variables. The Legendre-Fenchel transformation is then executed on those variables solely. As a result, we obtain partially bi-Lipschitz functions. This idea of incomplete convexity can be realized by the partially input convex neural network (PICNN) proposed by Amos et al. (2017) . A PICNN f (x, y; θ), convex only with respect to y with L layers, can be formulated as follows:u i+1 =g i Wi u i + bi , z i+1 =g i W (z) i z i ⊙ [W (zu) i u i + b (z) i ] + + W (y) i y ⊙ W (yu) i u i + b (y) i + W (u) i u i + b i (i = 0, . . . , L -1), f (x, y; θ) =z L , u 0 = x,where ⊙ denotes the Hadamar product, only W (z) i are non-negative, g i are convex non-decreasing and W (z) 0 = 0.In this paper, we mainly discussed notions of Lispchitzness, inverse Lispchitzness, strong convexity and smoothness in terms of the L 2 -norm. However, in some cases, it may be more interesting to work in other norms such as the L ∞ norms (Zhang et al., 2022) . In this extension, we briefly discuss how other norms can be introduced in our theoretical framework and its advantages. The previous extension can be regarded as such an example where the L 2 norm was changed to a weighted L 2 norm. Here, we treat the more general case. In this section, ∥ • ∥ will denote a general norm.The dual norm can be defined as follows. Definition F.2. Let ∥ • ∥ be a norm. Its dual ∥ • ∥ * is defined as follows:∥y∥ * := sup {⟨x, y⟩ | ∥x∥ = 1}For example, the dual of the L p -norm for any p ≥ 1 is the L q -norm such that 1/p + 1/q = 1. Definitions of strong convexity and smoothness are accordingly modified as follows:Definition F.3. Let µ > 0. F : R m → R is µ-strongly convex with respect to a norm ∥ • ∥ if for all t ∈ [0, 1] F (tx + (1 -t)y) ≤ tF (x) + (1 -t)F (y) - µ 2 t(1 -t)∥x -y∥ 2 . Definition F.4. Let γ > 0. F : R m → R is γ-smooth with respect to a norm ∥ • ∥ if F is differentiable and F (y) ≤ F (x) + ⟨∇F (x), y -x⟩ + γ 2 ∥x -y∥ 2 .Then, the following theorem holds: Theorem F.5 (Shalev-Shwartz and Singer (2010), Lemma 18). If F is a closed and µ-strongly convex function with respect to a norm ∥ • ∥, then its Legendre-Fenchel transformation F * is 1/µ-smooth with respect to the dual norm ∥ • ∥ * .Following Theorem D.2, the Legendre-Fenchel transformation of a closed and α-strongly convex function with respect to a norm ∥ • ∥ satisfies:EQUATIONThis means that by taking the derivative, we obtain a monotone Lipschitz function with respect to the dual norm.If ∥ • ∥ = ∥ • ∥ 1 , then ∥ • ∥ * = ∥ • ∥ ∞and we arrive at a Lipschitzness in terms of the L ∞ norm which is preferred in some situations (Zhang et al., 2022) . We can thus create convex smooth functions with a relatively wide choice on the norm.In many applications, since all norms are equivalent in a finite dimensional vector space, it may not be interesting to be able to deal with different norms. 3 In other words, controlling the of a function f in terms of a norm ∥ • ∥ can be executed simply through the calculation of the L 2 -norm asC 1 ∥x -y∥ ≤ ∥x -y∥ 2 ≤ C 2 ∥x -y∥.where C 1 and C 2 are constant independent of x. However, it turns out that usually C 1 and C 2 are dependent on the dimension. For example, between the L p and L q norms where p ≥ q, the following holds: ∥x∥ p ≤ ∥x∥ q ≤ n 1/q-1/p ∥x∥ p .As a result, in high dimensional spaces, the equivalence of norms may become impractical as we will need to deal with extremely large or small values when translating the requirements into L 2 norm, making the training unstable. The main advantage of our approach is that we can avoid this as we can directly characterize the function in any desired norm without any dependency on the dimension.On the other hand, the control of the lower bound (inverse Lipschitzness) becomes more difficult. Indeed, when p = 1 or p > 2, the L p norm is not strongly convex with respect to its own norm (Acu et al., 2023) . As a result, we cannot similarly proceed like our approach for the L 2 norm and introduce the generalized inverse Lipschitzness.Note that adding the squared norm α∥x∥ 2 2 /2 to a function F * satisfying equation ( 17), will lead for the lower bound to α∥x -y∥ 2 2 ≤ ⟨∇(F * (x) + α∥x∥ 2 /2) -∇(F * (y) + α∥y∥ 2 /2), x -y⟩ but for the upper bound to⟨∇(F * (x) + α∥x∥ 2 /2) -∇(F * (y) + α∥y∥ 2 /2), x -y⟩ ≤β∥x -y∥ 2 + α∥x -y∥ 2 2 ≤(C 2 α + β)∥x -y∥ 2 ,where the constant C 2 of the left hand side is highly dependent on the dimension. In high dimension, this may attenuate the importance of β, making control of the Lipschitzness difficult and unclear in practice.The above discussion is rather theoretical. We let concrete applications, as well as deeper investigations for future work. However, this may be a promising avenue since we can use other norms and the dependency on the dimension is completely dropped.Superposition As mentioned in the main paper, while the expressive power of our bi-Lipschitz unit is constrained to α + β-Lipschitz, α-strongly monotone functions which are themselves the derivative of a real-valued function, this can be alleviated by superposing several BLNNs. For example, we have tested to fit the sign function, and the composition of 5 BLNNs could achieve an accuracy around 0. Combination Another solution to improve the expressive power of our BLNN is to combine it with other architectures. While there may be various approaches, our model could be used as a pre-processing module for difficult tasks since the BLNN is by definition bi-Lipschitz which means that the geometric properties of the input data is relatively preserved in the output as well. Therefore, for some problems, we could implement a model where we first process the data through a BLNN and then transfer it to other networks. For example, we conducted an experiment using the convolutional version of BLNN (BLNNconv) on the problem of uncertainty estimation (as in Subsection 4.2) with the CIFAR-10 vs. SVHN dataset to illustrate the scalability of our method (Van Amersfoort et al., 2020) . For this problem, we implemented the model so that we first process the data through BLNNconv and then transfer it to the DUQ. The result compared to DUQ can be found in Table 4 . Our model is not only scalable to large-scale networks but also improves out-of-detection performance (the AUROC of SVHN). Using BLNNconv instead of the fully connected BLNN also improved the computation time (e.g., 2.5 times faster for the 5 first iterations).In this work, the Legendre-Fenchel transformation was used in order to provide a direct and single parameterization of the Lipschitzness in contrast to prior methods that controlled it on a layer-wise basis. At the expense of losing this benefit, it is also possible to characterize Lipschitzness at each layer. This may be more practical and faster to apply in some contexts. In this case, we will have to control the spectral norm of all the weights. Since we calculate the derivative, the relation between the overall Lipschitz constant and that of each layer is complex though. We provide it in the following proposition as reference.Proposition F.6. Define the ICNN G θ as followsz i+1 = g i (W (z) i z i + W (y) i y + b i ) (i = 0, . . . , k -1), G θ (y) = z k ,where g i is a 1-Lipschitz function, with bounded and Lipschitz derivative. Then the Lipschitz constant of ∇G θ can be obtained by solving the following set of inequalities, where α k is the Lipschitz constant:α 1 =Lip(∇g 1 )∥W (y) 1 ∥ 2 2 β 1 =∥W (y) 1 ∥ 2 α i+1 ≤∥∇g i ∥ ∞ ∥∥W (z) i ∥ 2 α i + Lip(∇g i )∥W (z) i ∥ 2 2 β i + 2Lip(∇g i )∥W (y) i ∥ 2 ∥W (z) i ∥ 2 β i + Lip(∇g i )∥W (y) i ∥ 2 2 β i+1 ≤∥W (z) i ∥ 2 β i + ∥W (y) i ∥ 2 (i = 1, . . . , k -1).Proof. Let us denote, z i (y) the i-th layer of the ICNN with input y. By construction,∥∇ y z i+1 (y 1 )-∇ y z i+1 (y 2 )∥ = ∇g i W (z) i z i (y 1 ) + W (y) i y 1 + b i W (z) i ∇ y z i (y 1 ) + W (y) i -∇g i W (z) i z i (y 2 ) + W (y) i y 2 + b i W (z) i ∇ y z i (y 2 ) + W (y) i . By defining h i (y) := ∇g i W (z) i z i (y) + W (y)i y + b i , we can further develop as follows:∥∇ y z i+1 (y 1 ) -∇ y z i+1 (y 2 )∥ ≤ h i (y 1 )W (z) i ∇ y z i (y 1 ) -h i (y 2 )W (z) i ∇ y z i (y 2 ) + h i (y 1 )W (y) i -h i (y 2 )W (y) i.The second term of the right hand side becomes:∥h i (y 1 ) W (y) i -h i (y 2 )W (y) i ≤Lip(∇g i ) W (z) i z i (y 1 ) + W (y) i y 1 + b i -W (z) i z i (y 2 ) + W (y) i y 2 + b i ∥W (y) i ∥ ≤Lip(∇g i )∥W (y) i ∥ ∥W (z) i ∥∥z i (y 1 ) -z i (y 2 )∥ + ∥W (y) i ∥∥y 1 -y 2 ∥ .As for the first term,∥h i (y 1 ) W (z) i ∇ y z i (y 1 ) -h i (y 2 )W (z) i ∇ y z i (y 2 ) = h i (y 1 )W (z) i ∇ y z i (y 1 ) -h i (y 2 )W (z) i ∇ y z i (y 1 ) +h i (y 2 )W (z) i ∇ y z i (y 1 ) -h i (y 2 )W (z) i ∇ y z i (y 2 ) ∥h i (y 1 ) -h i (y 2 )∥∥W (z) i ∥∥∇ y z i ∥ + ∥∇g i ∥ ∞ ∥W (z) i ∥∥∇ y z i (y 1 ) -∇ y z i (y 2 )∥ ≤Lip(∇g i ) ∥W (z) i ∥∥z i (y 1 ) -z i (y 2 )∥ + ∥W (y) i ∥∥y 1 -y 2 ∥ ∥W (z) i ∥∥∇ y z i ∥ + ∥∇g i ∥ ∞ ∥W (z) i ∥∥∇ y z i (y 1 ) -∇ y z i (y 2 )∥,where in the last inequality we used the definition of h i (y).On the other hand, we have∥∇ y z i ∥ ≤ ∥W (z) i-1 ∥∥∇ y z i ∥ + ∥W (y)i-1 ∥. By replacing ∥∇ y z i ∥ by β i , we obtainβ i ≤ ∥W (z) i-1 ∥β i-1 + ∥W (y) i-1 ∥. Since z i (y) is differentiable, ∥z i (y 1 ) -z i (y 2 )∥ ≤ β i .Finally, by setting α i := ∥∇ y z i+1 (y 1 ) -∇ y z i+1 (y 2 )∥/∥y 1 -y 2 ∥ and substituting the results into the first inequality, we obtain the desired result.Remark F.7. Note that in the above theorem g i has to be not only convex non-decreasing but also 1-Lipschitz with bounded and Lipschitz derivative. The layer-wise formulation requires thus further conditions on g i that were unnecessary in the BLNN (Algorithm 1). This means that our original model has more freedom and flexibility than the above simple approach.See Appendix H for further details on the setup of each experiment.First of all, we verify that BLNN behaves as expected with different values of α and β. We set all g i as ReLU or softplus functions. The Legendre-Fenchel transformation step is executed with GD (since it is simple and has good theoretical guarantees) until 1000 iterations are reached or the stopping condition, i.e., ∥∇ y (⟨y, x⟩ -F (y)) ∥ < 10 -3 , is satisfied. The BLNN has input dimension 2, and α is set to 4 as varying it does not severely impact the fundamental quality of the theoretical bounds. As we can observe, our model combined with gradient descent can generate functions with bi-Lipschitz constants that are inside the predefined bounds as expected for both softplus and ReLU activation functions. That is why, it seems reasonable to use the simple gradient descent in order to execute the Legendre-Fenchel transformation.Theoretically, since the initialization of weights are executed at random, the Lipschitz and inverse Lipschitz constants can be any value between α and α + β. Curiously though, the estimated Lipschitzness is close to the upper bound for almost every model. Note that we are generating the network and then directly estimating the bi-Lipschitz constants, which means this phenomenon is possibly closely related to the initialization of weights. To better understand this phenomenon, consider the linear network x → W x where W ∈ R n×n . In our architecture, this can be regarded as a simplification of ∇F θ , where F θ is the core ICNN. If we add x/β to this network and compute the inverse operation, corresponding to ∇F * θ , we obtain y → (W + I/β) -1 y. By letting σ min the smallest singular value of W , the Lipschitz constant becomes (σ min + 1/β) -1 . As a result, if the initialization provides a W with small σ min close to 0, we inevitably obtain a function with Lipschitz constant close to the upper bound β, independently of the distribution of the other eigenvalues. This seems to happen in Figures 11 and 12 . In other words, the construction employing the default initialization is biased, or more concretely, the derivative of an ICNN F θ with initialized weights is with high probability ϵ-inverse Lipschitz with ϵ << 1. This is supported by Figure 13 which calculated the bi-Lipschitz constants of a (4, 60)-BLNN initialized over 100 different trials. While the inverse Lipschitz constant is moderately distributed, the Lipschitz constant is concentrated on the theoretical maximum, i.e., 64.We also modified the initialization procedure of the non-negative weights W (z) i to confirm our hypothesis. The default approach was to draw each element from a uniform distribution proposed by Glorot and Bengio (2010) and clamp all negative elements to 0. Figure 14 As we can conclude, the choice of the initialization impacts the distribution of the starting bi-Lipschitz constants for a (α, β)-BLNN while theoretically any values are possible. In the default setting, we have a strong bias for the Lipschitz constant. Nevertheless, the fact that the bi-Lipschitz constants are close to the bounds can be regarded as the representation of a function with rich features. For example, if a matrix has the greatest eigenvalue and the smallest eigenvalue close to the limits, it implies that the other eigenvalues can assume any values in-between, allowing for considerable flexibility. This is a direct consequence of our meticulous construction of bi-Lipschitz neural networks respecting the desired constraints using convex neural networks and Legendre-Fenchel transformation.The adequate initialization of the weights of an ICNN in order to obtain BLNNs with a richer distribution of bi-Lipschitz constants is a whole new problem and outside the scope of our work. We will not investigate further but it seems to be an interesting research topic since an appropriate initialization may facilitate training. In this paper, we used the original setup. This subsection provides additional results of experiment 4.1 of the main paper. As a reminder, this experiment verifies the tightness of the proposed architecture by fitting the following function that has a discontinuity at the point x = 0:f (x) = x (x ≤ 0) x + 1 (x > 0)Therefore, when constraining the Lipschitzness of the model by a positive finite L, a model with tight Lipschitz bounds should achieve that upper bound around x = 0 in order to reproduce the behavior of f . We compare our model with other Lipschitz constrained architectures that were introduced in Appendix A. We take as comparison the spectral normalization (SN) (Miyato et al., 2018) , AOL (Prach and Lampert, 2022) , Orthogonal (Trockman and Kolter, 2021) , SLL (Araujo et al., 2023) and Sandwich (Wang and Manchester, 2023) . They can all regulate the Lipschitzness by an upper bound, but it is executed on a layer-wise basis.Models are built so that L is set to 2, 5, 10 and 50. For the compared prior works, this is realized by making all layers 1-Lipschitz, and scaling the input by L, which is the strategy often employed in practice. After learning was completed, the empirical Lipschitz constant was computed. The percentage between the empirical Lipschitz constant and the true upper bound L can be found in Table 5 . In order to be able to legitimately compare results, the number of parameters for each model were all of the same order. The difference in tightness is even clearer when setting the upper bound to 50 and plotting the result as shown in Figure 16 . Only our method can capture the presence of a discontinuous jump around 0. We can conclude that our method clearly achieves the tightest bound for Lipschitzness as well as the most coherent fit of the function.This subsection provides additional results of experiment 4.1 of the main paper. As a reminder, this experiment verifies the flexibility of our model. While we already provided a theoretical guarantee of the expressive power of our model, it is important to show that it can learn without problem functions with smaller Lipschitz and larger inverse Lipschitz constants than those we set as parameters. At the same time, this addresses one concern of the previous experiment since the tightness of our model may be interpreted as a result of the strong influence the regularization terms we add, effectively generating only α + β-Lipschitz functions. Therefore, we lead this experiment, where we overestimated the Lipschitz constant of the target function. The Lipschitz constant was set to 50 but the function we want to learn is a linear function with slope 1. Experimental setting was the same as the previous subsection. Results are shown in Figures 17 and 18 .Our model can without any problem fit the curve y = 50x and y = x even in the overestimation scenario. Interestingly, this is not the case for prior work. On the one hand, SN and AOL cannot learn the linear function with slope 50. On the other hand, others cannot learn the simple identity function even after convergence. While the overall trend is along y = x, we can observe small fluctuations. This may be due to the large scaling factor L we multiply at the input, leading to overfitting or interpolation which usually depends on the smoothness. That is, with high Lipschitzness, the function is likely to dynamically fluctuate around unseen points during training.In order to examine the influence of this scaling factor on the training we run the experiment again with even larger L, namely, 100 and 1000. The evolution of the loss and the fitted function for L = 1000 can be found in Figures 19 and 20 , respectively. As we can observe, only our method is hardly affected by the size of L and achieves a good generalization performance. While we do not have a rigorous mathematical explanation for this phenomenon, we believe our approach provides a better regularization than the previous methods. In other words, these layer-wise methods seem to provide a loss landscape that has many spurious minima where the function is driven. This phenomenon is aggravated with higher L. The difference in the smooth decrease of the loss function also supports this. The presence of these fluctuations due to the large scaling factor is thought to be related with the way we incorporated L into the model. For previous work, we had no choice but to scale the function by multiplying with L. However, our method is adding a regularization term without touching the core function. This is an interesting avenue for future research.In short, our model can perfectly learn both the function with slope 50 and 1 even though we greatly overestimated the Lipschitz constant.This experiment summarizes the two previous experiments. We propose to run a regression task of the function y = 50x with SN (representing the layer-wise methods) and our model, where the Lipschitz bound L is changed from 25 to 125. If the model provides tight bounds and perfect expressive power, then the loss should equal 0 for L ≥ 50 as theoretically y = 50x can be reconstructed, and should increase once L becomes smaller than 50 as the maximal slope we can reproduce is limited As we can observe, SN, a layer-wise model known to provide conservative bounds, only achieves a 0 loss from around L = 100, while ours exactly exhibits an increase of the loss from L = 50. This clearly shows that a layer-wise model has an expressive power that is largely lower than expected in the initial works. Notably, the convergence speed drastically decreases around L = 50 for our approach. This suggests choosing bi-Lipschitz constants as tight as possible tailored to the specific problem to ensure faster convergence.Interestingly, this difference of behavior around the true Lipschitz constant of our model could provide a new method to estimate the bi-Lipschitz constants. Indeed, by starting with a high Lipschitz constant α + β and decreasing it, we can find the largest point where the loss starts to increase. This value can become a good estimation of the Lipschitzness of the true function. We can similarly proceed for inverse Lipschitzness. We leave further investigation for future work. Since our model provide tight bounds and good regularization performance as shown in the precedent experiments, the annealing of its parameters could be effective in some suitable settings. Let us consider the case where we want to learn the exponential function. It is not Lipschitz in general since its slope is monotonically diverging to ∞ as x → ∞ and not inverse Lipschitz either since e x → 0 as x → -∞. We will start by a low Lipschitz parameter of the model and increase it by evaluating its effective Lipschitzness during training. The annealing will proceed in an elastic way. If the estimation is close to the upper bound then we will relax the bound since it means that our bound is too strong. That way, if we have a tight bound, the annealing process should progress faster and the overall function should rapidly converge to exp(x). The inverse Lipschitz constant is set to 0 and kept fixed. We compared our model with the Sandwich model. See Figure 22 . We can verify that the annealing process is faster than the layer-wise model. Ours constantly wants to increase the Lipschitz constant while the Sandwich layer abandons at some point and cannot keep up the increase of the Lipschitzness.This subsection provides additional results of experiment 4.2 of the main paper. As a reminder, we compare the uncertainty when learning the two moons dataset with not only DUQ and DUQ+BLNN but also the Deep Ensembles method and a DUQ with no bi-Lipschitz regularization. Results are plotted in Figure 23 . Uncertainty is calculated as the distance to the closest centroid (see Appendix E for the mathematical formulation). Blue indicates high uncertainty, and yellow low uncertainty.As mentioned in the main paper, DUQ+BLNN achieves a tighter bound which is possible because we are able to tune the bi-Lipschitz constant unlike DUQ. The experimental setup for DUQ was chosen like the original paper (Van Amersfoort et al., 2020) . Table 6 provides details on the accuracy of the learned model with different α and β, used for the grid search. As for the accuracy, the value of α does not matter as long as β is chosen large enough. Indeed, a BLNN with (α, β) = (0.0, 3.0) theoretically includes that with (α, β) = (1.0, 2.0). That is why they exhibit similar performance. Nevertheless, a too loose setting of these parameters, especially for β, still seems to affect the training, and a as tight as possible choice is beneficial as we can see for (α, β) = (2.0, 5.0). In a sense, decreasing β which is mainly in charge of the Lipschitz constant helps better generalization as pointed out in prior work as well (Bartlett et al., 2017) . In the same direction, a too large value of α leads to a nonsensical result as we can understand when α = 5. In this case, the model is randomly guessing which means that it is unsure everywhere. This is also translated in the uncertainty estimation as shown in Figure 24 . The area where the model is unsure increases as α is increased, meaning that it becomes more sensitive to unknown data (see Figure 25 ). Comparatively, setting α = 0 does not affect so much since the injectivity constraint is already beneficial in this type of task. Finally, a too low value of β means that the inverse Lipschitz and Lipschitz constants are close to each other and the function behaves more like a linear function. Consequently, learning becomes difficult as we can conclude from the table with β = 0.1. Similar trends could be also observed in more complex tasks such as FashionMNIST vs MNIST and FashionMNIST vs NotMNIST. This subsection provides additional results of experiment 4.2 of the main paper. In this experiment, we use real world data of FashionMNIST (Xiao et al., 2017) , MNIST (LeCun and Cortes, 2010) and NotMNIST (Bulatov, 2011) . The task is to learn to classify FashionMNIST but at the same time to detect out-of-distribution points from MNIST and NotMNIST. This task to distinguish FashionMNIST from MNIST datasets is known to be a complicated task (Van Amersfoort et al., 2020) . We compute the AUROC for the detection performance. Results for a downsampled dataset and full size dataset are shown in In this appendix, we provide further experimental details, including values of hyper-parameter, specific composition of used architectures and optimization schemes, which were omitted in the main text as a concern of clarity. All experiments were executed with Pytorch. We mainly used the GPU NVIDIA A100 with a memory of 80GB for the computation. Additional information for each experiments can be found in each subsection below and in the supplementary material. In this experiment, we were interested in the evolution of the bi-Lipschitz constants under different optimization schemes. Figures 8 and 9 .Data We generated 5000 2-dimensional points uniformly sampled at random in the interval [-1, 9]× [-1, 9].Architecture Our BLNN was created with β = 10. The value of α does not matter since we focused on the output of the optimization process which occurs before adding the regularization term α/2∥x∥ 2 . We used the softplus function as the activation function. The effective Lipschitz and inverse Lipschitz constants were estimated by a simple sampling with sample size 5000 chosen from the data.Optimization Schemes We used several optimization schemes, namely, steepest gradient descent (GD), Nesterov's accelerated gradient (AGD) (Nesterov, 1983) , Adagrad (Duchi et al., 2011) , RM-Sprop (Hinton et al., 2012) , Adam (Kingma and Ba, 2015) and the Newton method For Adagrad, RMSprop and Adam, all parameters were set as default except the learning rate.In this experiment, we were interested in the influence of the initial point on the convergence speed of the Legendre-Fenchel transformation. We used as a starting point either the final point of the previous epoch for each training data or the point (1, . . . , 1) ⊤ independently of the history. This experiment was run with the architecture explained in Subsection H.4.2. The codes were implemented in Python 3.11.4 with PyTorch 2.0.1+cu117. Data The FashionMNIST dataset was generated using the torchvision toolkit, with 60000 training data and 2000 test data. Batch size was set to 128. We used the original size of the dataset, namely 28 × 28.Optimization We used SGD with a learning rate of 0.05, momentum of 0.9 and weight decay of 10 -4 . We used the binary cross entropy as a loss function.Architecture We compared two models DUQ and DUQ with our method (DUQ+BLNN). The latter means that we used the DUQ framework but we deleted the gradient penalty of the loss function and directly replaced the vanilla neural network with our BLNN. We used the softplus function for DUQ+ours and ReLU for the other architectures. See Table 14 for further details on the architecture. As for our model, the best values of α and β were found through a grid search: α = 0.0 and β = 3.0.Particularly, we only used one BLNN in this experiment compared to the previous one, and projected it to a smaller sub-space of size 50 by a diagonal rectangular matrix with all components set to 1.In this experiment, we were interested in applying PBLNN to several tasks where data exhibited monotone behaviors with respect to some variables. The codes were imported and modified from (Nolte et al., 2023) .10 They were implemented in Python 3.10.14 with PyTorch 2.0.1+cu117.Corresponding Figures and Tables Table 3. Data COMPAS (Angwin et al., 2016 ) is a classification task with 13 features, where 4 of them have monotone inductive bias. It is important to highlight that this dataset is known to be racially biased (Angwin et al., 2016) . BlogFeedBack (Nolte et al., 2023 ) is a prediction task with 276 features, where 8 of them have monotone bias. LoanDefaulter (Nolte et al., 2023) a classification task with 28 features, where 5 of them have monotone inductive bias. HeartDisease (Janosi et al., 1988) is also a classification task with 13 features, including 2 with monotone inductive bias. AutoMPG (Quinlan, 1993) contains 7 features where 2 of them shows monotone inductive bias. CIFAR101 is an augmented dataset used by Nolte et al. (2023) .Optimization We used the same training scheme as Nolte et al. (2023) .Architecture We compared our model, the PBLNN, with that of Nolte et al. (2023) and that of Kim and Lee (2024) . As for the former, we constructed the PBLNN based on the inductive bias. For HeartDisease and AutoMPG, we took as the final output the average of the output of one PBLNN with bi-Lipschitzness imposed on the monotone variables. For COMPAS and LoanDefaulter, we created independent 1-dimensional PBLNNs whose bi-Lipschitzness is imposed on each variable that are monotone and then took the average. For CIFAR101, we used 101 PBLNNs with respect to the last variable. See Table 15 and the codes for further details.• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.38th Conference on Neural Information Processing Systems (NeurIPS 2024).We also tried to combine DUQ with BiLipNet but could not find competitive parameter settings.Here, a softplus-type activation function is a function s that satisfies the following conditions: s ≥ ReLU, s is convex and |s(x) -ReLU(x)| → 0 as |x| → ∞(Chen et al., 2019).Any two norms ∥ • ∥ (1) and ∥ • ∥ (2) on a finite dimensional space X over C are equivalent. That is, there exist constants C1 and C2 so that for all x ∈ X, C1∥x∥ (1) ≤ ∥x∥ (2) ≤ C2∥x∥ (1) .https://github.com/locuslab/deqhttps://github.com/MrYxJ/calculate-flops.pytorchhttps://github.com/y0ast/deterministic-uncertainty-quantificationhttps://github.com/niklasnolte/monotonic_tests 


Decision: This paper presents a novel bi-Lipschitz network approach. It is heavily theoretical with proof of concept ML experiments. Reviewers agreed that it is original and insightful. Productive rebuttal period yielded meaningful improvements and added essential information (e.g. time complexity). If accepted, please make sure to incorporate them.

Reviews: review: paper_summary: paper_summary: This paper proposes to control the bi-Lipschitzness of a neural-network by parameterizing the output by the Legendre-Fenchel-Dual. This involves parameterizing a strongly convex function and computing the minimum of that function in the forward pass. Several benchmarks are studied in simple regression tasks and uncertainty quantification.
strength_weakness: strength_weakness: -Computing the LFT seems to be quite expensive, hence why the experiments are only on simple 2d problems and fashion-MNIST. For this reason I'm doubtful that it will be used for any large-scale network training pipelines where tight Lipschitz control and estimation is challenging.

-The provable approximation class is limited to alpha-strongly monotone functions and is the derivative of some function almost everywhere. Lipschitz layers like AOL, SLL and Sandwich layer are all solutions to the LipSDP framework which only requires the activations themselves to be alpha-strongly monotone for alpha >= 0 (Fazlyab et al., 2019).
questions: questions: -Is it also necessary that BLNN is a strongly monotone function? It seems that many of the regression experiments involve monotone target functions (figure 2 and 3), but I'm not sure if that is because BLNN is not capable of representing monotone functions or just not a great representer due to your approximation theorem. If it can represent non-monotone functions, it would be interesting to see a simple regression comparison to SLL, AOL, etc. Answering this question will greatly help my evaluation.

-The Lipschitz parameterization of SLL, AOL, and Sandwich layers commonly uses compositions of 1-Lipschitz layers for the application of certified robustness. How would the BLNN parameterization compare to existing 1-Lipschitz layer networks in the certified robustness setting? I’d imagine BLNN might be much more expressive than composed 1-Lipschitz layers which could have a big impact.

-Have you considered amortizing the LFT computation as done in this paper? https://arxiv.org/abs/2210.12153 

-I'm curious if there is any possibility of extending BLNN to convolutional layers? These settings are interesting for larger image classification problems like CIFAR10 and Imagenet.
limitations: limitations: Limitation are adequately addressed.
paper_summary: paper_summary: This paper proposes a novel neural network architecture called BLNN (Bi-Lipschitz Neural Network) that allows direct control and parameterization of the overall bi-Lipschitzness of the network. The main contributions include: i) a framework that allows tight control of Lipschitz and inverse Lipschitz constants of networks via using convex neural networks and the Legendre-Fenchel transformation, ii) comprehensive theoretical analysis, iii) empirical evaluation showing the nice performance of BLNN on tasks like function fitting, out-of-distribution detection, and monotone regression.
strength_weakness: strength_weakness: 1. Computational Complexity: A detailed analysis of time and space complexity compared to traditional networks can be helpful.

2. Scalability and Practical Implications: There's insufficient exploration of how the method scales to very large networks or complex datasets (e.g. TinyImageNet).

3. Hyperparameter Sensitivity: More discussions on this issue will be beneficial.

4. The paper could be more explicit about scenarios where the theoretical guarantees might not hold, and could explore potential extensions to other network architectures beyond feedforward networks.
questions: questions: 1. How does the proposed method perform on larger, more complex datasets like TinyImageNet or ImageNet?

2. Can the authors clarify the computational complexity of their approach?

3. Can the authors provide a more comprehensive study on hyperparameter sensitivity?

4. Can the authors comment on other network structures (e.g. implicit models, DEQs, etc)?
limitations: limitations: It seems that an improved discussion on potential negative societal impacts or broader ethical considerations is still missing.
paper_summary: paper_summary: This paper investigates and proposes a novel bi-Lipschitz neural network architecture. This architecture provides a simple, direct and tight control of the Lipschitz and inverse Lipschitz constants through the use of two parameters, the ideal minimum, equipped with theoretical guarantees. To devise their architecture the authors exploit convex neural networks and the Legendre-Fenchel duality. The authors also propose a variant of their bi-Lipschitz architecture that is more scalable by exploiting partially input convex neural networks. Finally, the authors propose a set of experiments to showcase the utility of our model in concrete machine learning applications, namely, uncertainty estimation and monotone problem settings and show that it can improve previous methods.
strength_weakness: strength_weakness: - It would be interesting of the authors could provide experiments with both their architectures with respect to computational cost, and highlight time of training etc.
questions: questions: NA.
limitations: limitations: Lack of experiments wrt to computational cost and maybe an experiment with a more larger dataset than the current ones use in the paper.

Comments: comment: Thank you very much for considering our rebuttal and for your positive evaluation of our work! We greatly appreciate your insight regarding the amortization of LFT, as we also believe this could be a solution to improve the scalability of our method. We will mention this in the main paper and further analyze it in future work.
Thank you for clarifying the representation power of BLNN and its application to convolutional. I agree that the parameterization is quite convenient when considering compositions of BLNN blocks over other Lipschitz constrained layers. I think the framework is generally quite interesting so I will raise my score.

Regarding my question about amortization of LFT: I would certainly not expect you to produce these results during the rebuttal phase, but I thought it could be a helpful reference for alleviating the computational costs of BLNN.
Thank you once again to all the reviewers for evaluating our work. We extend additional thanks to those reviewers who have already reviewed our rebuttal, as well as to the area chair for managing the review process of our paper. We remain available until the end of the discussion period to discuss and provide further clarifications if any reviewers have remaining or additional concerns regarding our paper.
Thank you very much for considering our rebuttal! We greatly appreciate your positive position on the contributions of our paper.
Thank you for the rebuttal. I believe this paper to provide a good contribution and maintain my score.
Thank you very much for considering our rebuttal and for your positive position on the acceptance of our paper!
Thanks for the detailed response. I maintain my positive evaluation of this paper.