Newton Losses: Using Curvature Information for Learning with Differentiable Algorithms

Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 09:52:21

Abstract: When training neural networks with custom objectives, such as ranking losses and shortest-path losses, a common problem is that they are, per se, non-differentiable. A popular approach is to continuously relax the objectives to provide gradients, enabling learning. However, such differentiable relaxations are often non-convex and can exhibit vanishing and exploding gradients, making them (already in isolation) hard to optimize. Here, the loss function poses the bottleneck when training a deep neural network. We present Newton Losses, a method for improving the performance of existing hard to optimize losses by exploiting their second-order information via their empirical Fisher and Hessian matrices. Instead of training the neural network with second-order techniques, we only utilize the loss function's second-order information to replace it by a Newton Loss, while training the network with gradient descent.  This makes our method computationally efficient. We apply Newton Losses to eight differentiable algorithms for sorting and shortest-paths, achieving significant improvements for less-optimized differentiable algorithms, and consistent improvements, even for well-optimized differentiable algorithms.

Main paper: 
 introduction: 
 Traditionally, fully-supervised classification and regression learning relies on convex loss functions such as MSE or cross-entropy, which are easy-to-optimize by themselves. However, the need for large amounts of ground truth annotations is a limitation of fully-supervised learning; thus, weaklysupervised learning has gained popularity. Rather than using fully annotated data, weakly-supervised learning utilizes problem-specific algorithmic knowledge incorporated into the loss function. For example, instead of supervising ground truth values, supervision can be given in the form of ordering information. However, incorporating such knowledge into the loss can make it difficult to optimize, e.g., by making the loss non-convex in the model output, non-monotonic, introducing bad local minima, vanishing gradients, and/or exploding gradients, negatively impacting training [1] , [2] .Loss functions that integrate problem-specific knowledge can range from rather simple contrastive losses [3] to rather complex losses that require the integration of differentiable algorithms [4] . In this work, we primarily focus on the (harder) latter category, which allows for solving specialized tasks such as inverse rendering [5] - [7] , learning-to-rank [2] , [8] - [13] , self-supervised learning [14] , differentiation of optimizers [15] , [16] , and top-k supervision [10] , [12] , [17] . In this paper, we summarize these loss functions under the umbrella of algorithmic losses [18] as they introduce algorithmic knowledge into the training objective.While the success of neural network training is primarily due to the backpropagation algorithm and stochastic gradient descent (SGD), there has also been a strong line of work on second-order optimization for neural network training [19] - [28] . Compared to first-order methods like SGD, second-order optimization methods exhibit improved convergence rates and therefore require fewer training steps; however, they have two major limitations [26] , namely (i) computing the inverse of the curvature matrix for a large and deep neural network is computationally substantially more expensive than simply computing the gradient with backpropagation, which makes second-order methods practically inapplicable in many cases [29] ; (ii) networks trained with second-order information exhibit reduced generalization capabilities [30] .Inspired by ideas from second-order optimization, in this work, we propose a novel method for incorporating second-order information into training with non-convex and hard to optimize algorithmic losses. Loss functions are usually cheaper to evaluate than a neural network. Further, loss functions operate on lower dimensional spaces than the parameters of neural networks. If the loss function becomes the bottleneck in the optimization process because it is difficult to optimize, it suggests itself to use a stronger optimization method that requires fewer steps like second-order optimization. However, as applying second-order methods to neural networks is expensive and limits generalization, we want to train the neural network with first-order SGD. Therefore, we propose Newton Losses, a method for locally approximating loss functions with a quadratic via their second-order information. Thereby, Newton Losses provides a (locally) convex loss leading to better optimization behavior, while training the actual neural network with gradient descent. To clarify the notation, we note that, with "Newton Losses", we refer to the principle of converting an existing loss into a Newton loss.For the quadratic approximation of the algorithmic losses, we propose two variants of Newton Losses: (i) Hessian-based Newton Losses, which comprises a generally stronger method but requires an estimate of the Hessian. Depending on the choice of differentiable algorithm, choice of relaxation, or its implementation, the Hessian may, however, not be available. Thus, we also relax the method to (ii) empirical Fisher-based Newton Losses, which derive the curvature information from the empirical Fisher matrix, which only depends on the gradients. The Fisher variant can be easily implemented on top of existing algorithmic losses without consideration of computing second derivatives, while the Hessian variant requires computation of second derivatives and leads to greater improvements.We evaluate Newton Losses for an array of eight families of algorithmic losses on two popular algorithmic benchmarks: the four-digit MNIST sorting benchmark and the Warcraft shortest-path benchmark. We find that Newton Losses leads to performance improvements for each of the algorithms-for some of the algorithms (those which suffer the most from vanishing and exploding gradients) more than doubling the accuracy. In an ablation study, we demonstrate that Newton Losses does not harm training when applied to the convex cross-entropy loss in simple classification settings. 
background & related work: 
 The related work comprises algorithmic supervision losses and second-order optimization methods. To the best of our knowledge, this is the first work combining second-order optimization of loss functions with first-order optimization of neural networks, especially for algorithmic losses.Algorithmic Losses. Algorithmic losses, i.e., losses which contain some kind of algorithmic component, have become quite popular in recent machine learning research. In the domain of recommender systems, early learning-to-rank works already appeared in the 2000s [8] , [9] , [31] , but more recently Lee et al. [32] proposed differentiable ranking metrics, and Swezey et al. [13] proposed PiRank, which relies on differentiable sorting. For differentiable sorting, an array of methods has been proposed in recent years, which includes NeuralSort [33] , SoftSort [34] , Optimal Transport Sort [10] , differentiable sorting networks (DSN) [12] , and the relaxed Bubble Sort algorithm [18] . Other works explore differentiable sorting-based top-k for applications such as differentiable image patch selection [35] , differentiable k-nearest-neighbor [17] , [33] , top-k attention for machine translation [17] , differentiable beam search methods [17] , [36] , survival analysis [37] , and self-supervised learning [38] . But algorithmic losses are not limited to sorting: other works have considered learning shortest-paths [15] , [16] , [18] , learning 3D shapes from images and silhouettes [5] - [7] , [18] , [39] , [40] , learning with combinatorial solvers for NP-hard problems [16] , learning to classify handwritten characters based on editing distances between strings [18] , learning with differentiable physics simulations [41] , and learning protein structure with a differentiable simulator [42] , among many others.Second-Order Optimization. Second-order methods have gained popularity in machine learning due to their fast convergence properties when compared to first-order methods [19] . One alternative to the vanilla Newton's method are quasi-Newton methods, which, instead of computing an inverse Hessian in the Newton step (which is expensive), approximate this curvature from the change in gradients [26] , [43] , [44] . In addition, a number of new approximations to the pre-conditioning matrix have been proposed in the recent literature, i.a., [20] , [23] , [45] . While the vanilla Newton method relies on the Hessian, there are variants which use the empirical Fisher information matrix, which can coincide in specific cases with the Hessian, but generally exhibits somewhat different behavior. For an overview and discussion of Fisher-based methods (aka. natural gradient descent), see [46] , [47] . 
preliminaries: 
 We consider the training of a neural network f (x; θ), where x ∈ R n is the vector of inputs, θ ∈ R d is the vector of trainable parameters and y = f (x; θ) ∈ R m is the vector of outputs. As per vectorization, x = [x 1 , . . . , x N ] ⊤ ∈ R N ×n denotes a set of N input data points, and y = f (x; θ) ∈ R N ×m denotes the neural network outputs corresponding to the inputs. Further, let ℓ : R N ×m → R denote the loss function, and let the "label" information be implicitly encoded in ℓ. The reason for this choice of implicit notation is that, for many algorithmic losses, it is not just a label, e.g., it can be ordinal information or a set of encoded constraints. The training of a neural network may be expressed asEQUATIONwhere Θ ⊆ R d is the domain of the trainable parameters θ. We note that the formulation ( 1) is general and includes, e.g., optimization of non-decomposable loss functions (i.e., not composed of individual losses per training sample), which is relevant for some algorithmic losses like, e.g., ranking losses [2] .Typically, the optimization problem (1) is solved by using an iterative algorithm (e.g., SGD [48] , Adam [49] , or Newton's method [26] ) updating the weights θ by repeatedly applying a respective update step:EQUATIONHowever, in this work, we consider splitting this optimization update step into two alternating steps:EQUATIONEQUATIONThis split allows us later to use two different iterative optimization algorithms for (3a) and (3b), respectively. This is especially interesting for optimization problems where the loss function ℓ is non-convex and its minimization is a difficult optimization problem itself, and, in other words, those problems, where a more robust and stronger optimization method exhibits faster convergence.We note that a similar split (for the case of splitting between the layers of a neural network, and using gradient descent for both (3a) and (3b)) is also utilized in the fields of biologically plausible backpropagation [50] - [53] and proximal backpropagation [54] , leading to reparameterizations of backpropagation. In Supplementary Material D, we provide further discussion of the split for the special case of using the same optimizer for (3a) and (3b) from a theoretical perspective and show that: (1) for gradient descent the split does not change the parameter updates, (2) for Newton's method in a rank-1 setting the split does not change the parameter updates, and, (3) under simplifying assumptions, the sets of stationary points are invariant to the split for a broader class of deterministic strategies.Equipped with the two-step optimization method (3), we can now introduce Newton Losses:We propose optimizing (3a) with Newton's method, while optimizing (3b) with gradient descent.In the following, we formulate how we can solve optimizing (3a) with Newton's method, or, whenever we do not have access to the Hessian of ℓ, using a step pre-conditioned via the empirical Fisher matrix. This allows us to transform an original loss function ℓ into a Newton loss ℓ * , which allows optimizing ℓ * with gradient descent, making it suitable for common machine learning frameworks.We begin by considering the quadratic approximation of the loss function at the location y = f (x; θ)EQUATIONTo find the location z ⋆ of the stationary point of ly (z), we set its derivative to 0:EQUATIONHowever, when ℓ is non-convex or the smallest eigenvalues of ∇ 2 y ℓ(y) either become negative or zero, z ⋆ would not be the local minimum of ℓ but may instead be any other stationary point or lie far away from y, leading to exploding gradients downstream. To resolve this issue, we introduce Tikhonov regularization [55] with a strength of λ, which leads to a well-conditioned curvature matrix:EQUATIONUsing z ⋆ , we can plug the solution into (3b) to find the Newton loss ℓ * and compute its derivative asℓ * z ⋆ (y) = 1 2 (z ⋆ -y) ⊤ (z ⋆ -y) = 1 2 ∥z ⋆ -y∥ 2 2andEQUATIONVia this construction, we obtain the Newton loss ℓ * z ⋆ , a new convex loss, which itself has a gradient that corresponds to a Newton step of the original loss. Thus, when we optimizing the Newton loss with gradient descent, we obtain the proposed method. Definition 1 (Newton Losses (Hessian)). For a loss function ℓ and a given current parameter vector θ, we define the Hessian-based Newton loss via the empirical Hessian asℓ * z ⋆ (y) = 1 2 ∥z ⋆ -y∥ 2 2wherez ⋆ i = ȳi -1 N N j=1 ∇ 2 ȳj ℓ(ȳ) + λ • I -1 ∇ ȳi ℓ(ȳ)for all i ∈ {1, ..., N } and ȳ = f (x; θ).We remark that computing and inverting the Hessian of the loss function is usually computationally efficient. However, whenever the Hessian matrix of the loss function is not available, whether it may be due to limitations of a differentiable algorithm, large computational cost, lack of a respective implementation of the second derivative, etc., we may resort to using the empirical Fisher matrix as a source for curvature information. We remark that the empirical Fisher matrix is not the same as the Fisher information matrix [46] , and that the Fisher information matrix is generally not available for algorithmic losses. While the empirical Fisher matrix, as a source for curvature information, may be of lower quality than the Hessian matrix, it has the crucial advantage that it can be readily computed from the gradients, specifically,EQUATIONThis means that, assuming a moderate dimension of the prediction space m, computing the empirical Fisher comes at no significant overhead and may, conveniently, be performed in-place as we discuss later. Again, we regularize the matrix via Tikhonov regularization with strength λ and can, accordingly, define the empirical Fisher-based Newton loss as follows.Definition 2 (Newton Loss (Fisher)). For a loss function ℓ, and a given current parameter vector θ, we define the empirical Fisher-based Newton loss asℓ * z ⋆ (y) = 1 2 ∥z ⋆ -y∥ 2 2 where z ⋆ i = ȳi -1 N N j=1 ∇ ȳj ℓ(ȳ) ∇ ȳj ℓ(ȳ) ⊤ + λ • I -1 ∇ ȳi ℓ(ȳ)for all i ∈ {1, ..., N } and ȳ = f (x; θ).Before continuing with the implementation, integration, and further computational considerations, we can make an interesting observation. In the case of using the trivial MSE loss, i.e., ℓ(y) = 1 2 ∥y -y ⋆ ∥ 2 2 where y ⋆ denotes a ground truth, the Newton loss collapses to the original MSE loss. This demonstrates that Newton losses require non-trivial original losses. Another interesting aspect is the arising fixpoint-the Newton loss of a Newton loss is equivalent to a simple Newton loss.After introducing Newton Losses, in this section, we discuss aspects of implementation and integration into differentiable algorithm-based training architectures. In particular, we illustrate its integration in Algorithms 1 and 2. Whenever we have access to the Hessian matrix of the algorithmic loss function, it is generally favorable to utilize the Hessian-based approach (Algo. 1 / Def. 1), whereas we can utilize the empirical Fisher-based approach (Algo. 2 / Def. 2) in any case.In Algorithm 1, the differences to a regular training is that we use the original loss only for the computation of the gradient (grad) and the Hessian matrix (hess) of the original loss. Then we compute z ⋆ (z_star). Here, depending on the automatic differentiation framework, we need to ensure def backward(ctx, g): fisher = g.T @ g * g.shape[0] input_grad = g @ inverse(fisher + ctx.tik_l * eye(g.shape [1] )) return input_grad, None for data, label in data_loader:# apply a neural network model y = model(data) # inject the Fisher backward mod. y = InjectFisher.apply(y, tik_l) # compute the original loss l = loss(y, label) # backpropagate and optim. step l.backward() optimizer.step() not to backpropagate through the target z_star, which may be achieved, e.g., via ".detach()", ".stop_gradient()", or "with no_grad():" depending on the choice of library. Finally, the Newton loss l may be computed as a squared or MSE loss between the model output y and z_star and an optimization step on l may be performed.We note that, while we use a label from our data_loader, this label may be empty or an abstract piece of information for the differentiable algorithm; in our experiments, we use information about ordinal relationships between values as well as a shortest-path on a graph.In Algorithm 2, we sketch how to apply the empirical Fisher-based Newton Losses. In particular, due to the empirical Fisher matrix depending only on the gradient, we can compute it in-place during the backward pass / backpropagation, which makes this variant particularly simple and efficient to apply. This can be achieved via an injection of a custom gradient right before applying the original loss, which replaces the gradient in-place by a gradient that corresponds to Definition 2. The injection is performed by the InjectFisher function, which corresponds to an identity during the forward pass but replaces the gradient by the gradient of the respective empirical Fisher-based Newton loss.In both cases, the only additional hyperparameter to specify is the Tikonov regularization strength λ (tik_l). This hyperparameter heavily depends on the algorithmic loss function, particularly, on the magnitude of gradients provided by the algorithmic loss, which may vary drastically between different methods and implementations. Other factors may be the choice of Hessian / Fisher, the dimension of outputs m, the batch size N . Notably, for very large λ, the direction of the gradient becomes more similar to regular gradient descent, and for smaller λ, the effect of Newton Losses increases. We provide an ablation study for λ in Figure 4 for the ranking experiment. 
experiments: 
 For the experiments, we apply Newton Losses to eight methods for differentiable algorithms and evaluate them on two established benchmarks for algorithmic supervision, i.e., problems where an algorithm is applied to the predictions of a model and only the outputs of the algorithm are supervised. Specifically, we focus on the tasks of ranking supervision and shortest-path supervision because they each have a range of established methods for evaluating our approach. In ranking supervision, only the relative order of a set of samples is known, while their absolute values remain unsupervised. The established benchmark for differentiable sorting and ranking algorithms is the multi-digit MNIST sorting benchmark [2] , [10] , [12] , [33] , [34] . In shortest path supervision, only the shortest-path of a graph is supervised, while the underlying cost matrix remains unsupervised. The established benchmark for differentiable shortest-path algorithms is the Warcraft shortest-path benchmark [15] , [16] , [18] . As these tasks require backpropagating through conventionally nondifferentiable algorithms, the respective approaches make the ranking or shortest-path algorithms differentiable such that they can be used as part of the loss function. Finally, as an ablation study, we also apply Newton Losses to the trivial case of classification, where we do not expect a performance improvement, as the loss is not hard to optimize, but rather validate that the method still works. In this section, we explore ranking supervision [33] with an array of differentiable sorting-based losses. Here, we use the four-digit MNIST sorting benchmark [33] , where sets of n fourdigit MNIST images are given, and the supervision is the relative order of these images corresponding to the displayed value, while the absolute values remain unsupervised. The goal is to learn a CNN that maps each image to a scalar value in an order preserving fashion. As losses, we use sorting supervision losses based on the Neu-ralSort [33] , the SoftSort [34] , the logistic Differentiable Sorting Network [12] , and the monotonic Cauchy DSN [2] . NeuralSort and SoftSort work by mapping an input list (or vector) of values to a differentiable permutation matrix that is row-stochastic and indicates the order / ranking of the inputs. Differentiable Sorting Networks offer an alternative to NeuralSort and SoftSort. DSNs are based on sorting networks, a classic family of sorting algorithms that operate by conditionally swapping elements. By introducing perturbations, DSNs relax the conditional swap operator to a differentiable conditional swap and thereby continuously relax the sorting and ranking operators. We discuss the background of each of these differentiable sorting and ranking algorithms in greater detail in Supplementary Material A.Setups. The sorting supervision losses are cross-entropy losses defined between the differentiable permutation matrix produced by a respective differentiable sorting operator and the ground truth permutation matrix corresponding to a ground truth ranking. The Cauchy DSN may be an exception to the hard to optimize classification as it is quasi-convex [2] . We evaluate the sorting benchmark for numbers of elements to be ranked n ∈ {5, 10} and use the percentage of rankings correctly identified as well as percentage of individual element ranks correctly identified as evaluation metrics. For each of the four original baseline methods, we compare it to two variants of their Newton losses: the empirical Hessian and the empirical Fisher variant. For each setting, we train the CNN on 10 seeds using the Adam optimizer [49] at a learning rate of 10 -3 for 10 5 steps and a batch size of 100. Results. As displayed in Table 1 , we can see that-for each original loss-Newton Losses improve over their baselines. For NeuralSort, SoftSort, and Logistic DSNs, we find that using the Newton losses substantially improves performance. Here, the reason is that these methods suffer from vanishing and exploding gradients, especially for the more challenging case of n = 10. As expected, we find that the Hessian Newton Loss leads to better results than the Fisher variant, except for NeuralSort and SoftSort in the easy setting of n = 5, where the results are nevertheless quite close. Monotonic differentiable sorting networks, i.e., the Cauchy DSNs, provide an improved variant of DSNs, which have the property of quasi-convexity and have been shown to exhibit much better training behavior out-of-the-box, which makes it very hard to improve upon the existing results. Nevertheless, Hessian Newton Losses are on-par for the easy case of n = 5 and, notably, improve the performance by more than 1% on the more challenging case of n = 10. In summary, we obtain strong improvements on losses that are difficult-to-optimize, while on well-behaving losses the improvements are smaller. This aligns with our goal of improving performance on losses that are hard to optimize. In this section, we apply Newton Losses to the shortest-path supervision task of the 12 × 12 Warcraft shortest-path benchmark [15] , [16] , [18] . Here, 12 × 12 Warcraft terrain maps are given as 96 × 96 RGB images (e.g., Figure 2 left) and the supervision is the shortest path from the top left to the bottom right (Figure 2 right) according to a hidden cost embedding (Figure 2 center). The hidden cost embedding is not available for training. The goal is to predict 12 × 12 cost embeddings of the terrain maps such that the shortest path according to the predicted embedding corresponds to the ground truth shortest path. Vlastelica et al. [16] have shown that integrating an algorithm in the training pipeline substantially improves performance compared to only using a neural network with an easy-to-optimize loss function, which has been confirmed by subsequent work [15] , [18] . For this task, we explore a set of families of algorithmic supervision approaches: Relaxed Bellman-Ford [18] is a shortest-path algorithm relaxed via the AlgoVision framework, which continuously relaxes algorithms by perturbing all accessed variables with logistic distributions and approximating the expectation value in closed form. Stochastic Smoothing [56] is a sampling-based differentiation method that can be used to relax, e.g., a shortest-path algorithm by perturbing the input with probability distribution. Perturbed Optimizers with Fenchel-Young Losses [15] build on stochastic smoothing and Fenchel-Young losses [57] and identify the argmax to be the differential of max, which allows a simplification of stochastic smoothing, again applied, e.g., to shortest-path learning problems. We use the same hyperparameters as shared by previous works [15] , [18] . Relaxed Bellman-Ford.The relaxed Bellman-Ford algorithm [18] is a continuous relaxation of the Bellman-Ford algorithm via the Al-goVision library. To increase the number of settings considered, we explore four sub-variants of the algorithm: For+L 1 , For+L 2 2 , While+L 1 , and While+L 2 2 . Here, For / While refers to the distinction between using a While and For loop in Bellman-Ford, while L 1 vs. L 2 2 refer to the choice of metric between shortest paths. As computing the Hessian of the AlgoVision Bellman-Ford algorithm is too expensive with the PyTorch implementation, for this evaluation, we restrict it to the empirical Fisher-based Newton loss. The results displayed in Table 2 . While we can observe small improvements in three out out of four settings and the best result is achieved by Newton Losses, the differences are small. This can be attributed to (i) the high performance of the baseline algorithm on this benchmark, and (ii) only being able to use the empirical Fisher-based Newton loss which is expected to underperform compared to the Hessian Newton loss.Stochastic Smoothing. After discussing the analytical relaxation, we continue with stochastic smoothing approaches. First, we consider stochastic smoothing [56] , which allows perturbing the input of a function with an exponential family distribution to estimate the gradient of the smoothed function. For the baseline, we apply stochastic smoothing to a hard non-differentiable Dijkstra algorithm based loss function to relax it via Gaussian noise ("SS of loss"). We utilize variance reduction via the method of covariates. As we detail in Supplementary Material A.4, stochastic smoothing can also be used to estimate the Hessian of the smoothed function [56] . Based on this result, we can construct the Hessian Newton loss. As an extension to stochastic smoothing, we apply stochastic smoothing only to the non-differentiable Dijstra algorithm (thereby computing its Jacobian matrix) but use a differentiable loss to compare the predicted relaxed shortest-path to the ground truth shortest-path ("SS of algorithm"). In this case, the Hessian Newton loss is not applicable because the output of the smoothed algorithm is high dimensional and the Hessian of the loss becomes intractable. Nevertheless, we can apply the Fisher-based Newton loss. We evaluate both approaches for 3, 10, and 30 samples. In Table 3 , we can observe that Newton Losses improves the results for stochastic smoothing in each case with more than 3 samples. The reason for the poor performance on 3 samples is that the Hessian or empirical Fisher, respectively, is estimated using only 3 samples, which makes the estimate unstable. For 10 and 30 samples, the performance improves compared to the original method. In Figure 3 , we display a respective accuracy plot. When comparing "SS of loss" and "SS of algorithm", we can observe that the extension to smoothing only the algorithm improves performance for at least 10 samples. Here, the reason, again, is that smoothing the algorithm itself requires estimating the Jacobian instead of only the gradient; thus, a larger number of samples is necessary; however starting at 10 samples, smoothing the algorithm performs better, which means that the approach is better at utilizing a given sample budget.Perturbed Optimizers with Fenchel-Young Losses. Perturbed optimizers with a Fenchel-Young loss [15] is a formulation of solving the shortest path problem as an arg max problem, and differentiating this problem using stochastic smoothing-based perturbations and a Fenchel-Young loss. Berthet et al. [15] mentioned that their approach works well for small numbers of samples, which we can confirm as seen in Table 3 where the accuracy is similar for each number of samples. By extending their formulation to computing the Hessian of the Fenchel-Young loss, we can compute the Newton loss, and find that we can achieve improvements of more than 2%. However, for Fenchel-Young losses, which are defined via their derivative, the empirical Fisher is not particularly meaningful, leading to equivalent performance between the baseline and the Fisher Newton loss. An interesting observation is that perturbed optimizers with Fenchel-Young losses perform better than stochastic smoothing in the few-sample regime, whereas stochastic smoothing performs better with larger numbers of samples. Finally, as an ablation study, we want to ensure that Newton Losses do not impede training for the simple case of MNIST classification with a softmax cross-entropy loss. We remark that, as cross-entropy is convex and very well-behaved objective, we cannot expect improvements for this experiment. To facilitate a fair and extensive comparison, we benchmark training on 5 models and with 2 optimizers: We use 5-layer fully connected ReLU networks with 100 (M1), 400 (M2) and 1 600 (M3) neurons per layer, as well as the convolutional LeNet-5 with sigmoid activations (M4) and LeNet-5 with ReLU activations (M5). Further, we use SGD and Adam [49] as optimizers. To evaluate both early performance and full training performance, we test after 1 and 200 epochs. We run each experiment with 20 seeds, which allows us to perform significance tests (significance level 0.05). As displayed in Table 4 , we find that Newton Losses perform indistinguishable to regular training, specifically, in 9 out of 20 cases it is better and in 7 cases equal. We find that it is significantly better in a single case, which is to be expected from equal methods (on average 1/20 tests will be significant at a significance level of 0.05). In conclusion, we can observe that for the trivial case of classification with the convex cross-entropy loss, Newton Losses maintains the original training performance. We provide tables with runtimes for the experiments in Supplementary Material B. We can observe that the runtimes between the baseline and empirical Fisher-based Newton Losses are indistinguishable for all cases. For the analytical relaxations of differentiable sorting algorithms, where the computation of the Hessian can become expensive with automatic differentiation (i.e., without a custom derivation of the Hessian and without vectorized Hessian computation), we observed overheads between 10% and 2.6×. For all stochastic approaches, we observe indistinguishable runtimes for Hessian-based Newton Losses. In summary, applying the Fisher variant of Newton Losses has a minimal computational overhead, whereas, for the Hessian variant, any overhead depends on the implementation of the computation of the Hessian of the algorithmic loss function. For applications where the the output dimensionality m becomes very large such that inversion of the empirical Fisher becomes expensive, we refer to the Woodbury matrix identity [58] , which allows simplifying the computation via its low-rank decomposition. A corresponding deviation is included in Supplementary Material G. 
conclusion: 
 In this work, we focused on weakly-supervised learning problems that require integration of differentiable algorithmic procedures in the loss function. This leads to non-convex loss functions that exhibit vanishing and exploding gradients, making them hard to optimize. Thus, we proposed a novel approach for improving performance of algorithmic losses building upon the curvature information of the loss. For this, we split the optimization procedure into two steps: optimizing on the loss itself using Newton's method to mitigate vanishing and exploding gradients, and then optimizing the neural network with gradient descent. We simplified this procedure via a transformation of an original loss function into a Newton loss, which comes in two flavors: a Hessian variant for cases where the Hessian is available and an empirical Fisher variant as an alternative. We evaluated Newton Losses on a set of algorithmic supervision settings, demonstrating that the method can drastically improve performance for weakly-performing differentiable algorithms. We hope that the community adapts Newton Losses for learning with differentiable algorithms and see great potential for combining it with future differentiable algorithms in unexplored territories of the space of differentiable relaxations, algorithms, operators, and simulators. 
 
 Appendix: In this section, we extend the discussion of SoftSort, DiffSort, AlgoVision, and stochastic smoothing.A.1 SOFTSORT AND NEURALSORT SoftSort [34] and NeuralSort [33] are prominent yet simple examples of a differentiable algorithm. In the case of ranking supervision, they obtain an array or vector of scalars and return a row-stochastic matrix called the differentiable permutation matrix P , which is a relaxation of the argsort operator. Note that, in this case, a set of k inputs yields a scalar for each image and thereby y ∈ R k . As a ground truth label, a ground truth permutation matrix Q is given and the loss between P and Q is the binary cross entropy loss ℓ SS (y) = BCE (P (y), Q) . Minimizing the loss enforces the order of predictions y to correspond to the true order, which is the training objective. SoftSort is defined asEQUATION)where τ is a temperature parameter, "sort" sorts the entries of a vector in non-ascending order, ⊖ is the element-wise broadcasting subtraction, | • | is the element-wise absolute value, and "softmax" is the row-wise softmax operator, as also used in (32) (right). NeuralSort is defined similarly and omitted for the sake of brevity. In the limit of τ → 0, SoftSort and NeuralSort converge to the exact ranking permutation matrix [33] , [34] . A respective Newton loss can be implemented using automatic differentiation according to Definition 1 or via the Fisher information matrix using Definition 2.Differentiable sorting networks (DSN) [2] , [12] offer a strong alternative to SoftSort and NeuralSort. They are based on sorting networks, a classic family of sorting algorithms that operate by conditionally swapping elements [59] . As the locations of the conditional swaps are pre-defined, they are suitable for hardware implementations, which also makes them especially suited for continuous relaxation. By perturbing a conditional swap with a distribution and solving for the expectation under this perturbation in closed-form, we can differentially sort a set of values and obtain a differentiable doubly-stochastic permutation matrix P , which can be used via the BCE loss as in Section A.1. We can obtain the respective Newton loss either via the Hessian computed via automatic differentiation or via the Fisher information matrix.AlgoVision [18] is a framework for continuously relaxing arbitrary simple algorithms by perturbing all accessed variables with logistic distributions. The method approximates the expectation value of the output of the algorithm in closed-form and does not require sampling. For shortest-path supervision, we use a relaxation of the Bellman-Ford algorithm [60] , [61] and compare the predicted shortest path with the ground truth shortest path via an MSE loss. The input to the shortest path algorithm is a cost embedding matrix predicted by a neural network.Another differentiation method is stochastic smoothing [56] . This method regularizes a nondifferentiable and discontinuous loss function ℓ(y) by randomly perturbing its input with random noise ϵ (i.e., ℓ(y + ϵ)). The loss function is then approximated as ℓ(y)≈ ℓ ϵ (y) = E ϵ [ℓ(y + ϵ)].While ℓ may be non-differentiable, its smoothed stochastic counterpart ℓ ϵ is differentiable and the corresponding gradient and Hessian can be estimated via the following result. Lemma 1 (Exponential Family Smoothing, adapted from Lemma 1.5 in Abernethy et al. [56] ). Given a distribution over R m with a probability density function µ of the form µ(ϵ) = exp(-ν(ϵ)) for any twice-differentiable ν, thenEQUATIONEQUATIONA variance-reduced form of ( 11) and ( 12) isEQUATIONEQUATIONIn this work, we use this to estimate the gradient of the shortest path algorithm. By including the second derivative, we extend the perturbed optimizer losses to Newton losses. This also lends itself to full second-order optimization.A.5 PERTURBED OPTIMIZERS WITH FENCHEL-YOUNG LOSSES Berthet et al. [15] build on stochastic smoothing and Fenchel-Young losses [57] to propose perturbed optimizers with Fenchel-Young losses. For this, they use algorithms, like Dijkstra, to solve optimization problems of the type max w∈C ⟨y, w⟩, where C denotes the feasible set, e.g., the set of valid paths. Berthet et al. [15] identify the argmax to be the differential of max, which allows a simplification of stochastic smoothing. By identifying similarities to Fenchel-Young losses, they find that the gradient of their loss isEQUATIONwhere w ⋆ is the ground truth solution of the optimization problem (e.g., shortest path). This formulation allows optimizing the model without the need for computing the actual value of the loss function. Berthet et al. [15] find that the number of samples-surprisingly-only has a small impact on performance, such that 3 samples were sufficient in many experiments, and in some cases even a single sample was sufficient. In this work, we confirm this behavior and also compare it to plain stochastic smoothing. We find that for perturbed optimizers, the number of samples barely impacts performance, while for stochastic smoothing more samples always improve performance. If only few samples can be afforded (like 10 or less), perturbed optimizers are better as they are more sample efficient; however, when more samples are available, stochastic smoothing is superior as it can utilize more samples better.In this supplementary material, we provide and discuss runtimes for the experiments.In the differentiable sorting and ranking experiment, as shown in Table 5 , we observe that the runtime from regular training compared to the Newton loss with the Fisher is only marginally increased. This is because computing the Fisher and inverting it is very inexpensive. We observe that the Newton loss with the Hessian, however, is more expensive: due to the implementation of the differentiable sorting and ranking operators, we compute the Hessian by differentiating each element of the gradient, which makes this process fairly expensive. An improved implementation could make this process much faster. Nevertheless, there is always some overhead to computing the Hessian compared to the Fisher. In Table 6 , we show the runtimes for the shortest-path experiment with AlgoVision. Here, we observe that the runtime overhead is very small.In Table 7 , we show the runtimes for the shortest-path experiment with stochastic methods. Here, we observe that the runtime overhead is also very small. Here, the Hessian is also cheap to compute as it is not computed with automatic differentiation. For the experiments in the tables, select λ based one seed from the grid λ ∈ [0.001, 0.01, 0.1, 1, 10, 100, 100]. For the experiments in Tables 2 and 9 , we present the values in the Tables 8 and 9 , respectively. For the experiment in Table 2 , we use a Tikhonov regularization strength of λ = 1000, which was selected based on the first out of the four algorithms (For+L 1 ). Table 9 : Tikhonov regularization strengths λ for the experiment in Table 3 . We use the experimental setting from Table 1 for NeuralSort and SoftSort and n = 5 (top right, bottom left). Further, to ablate the hyperparameter of the base methods NeuralSort and SoftSort (τ ), we display the results for τ ∈ [1, 0.1]. Throughout the paper, we selected the optimal hyperparameter settings for each baseline, in particular, τ = 0.1 for NeuralSort and τ = 1 for SoftSort, and used the optimal baseline hyperparameter configurations also for Newton Losses.In the plots, we consider 13 values for λ from 0.001 to 1000. We can observe that Newton Losses are robust over multiple orders of magnitude for the hyperparameter λ.The plots are averaged over 2 runs and individual runs are indicated with reduced opacity [will be updated to 10 runs for the camera-ready].In this section, we consider the splitting of the optimization step from (2) into (3a) and (3b) in greater detail. First, we consider the equivalence between ( 2) and (3a) + (3b) for each case being gradient descent which has also been considered by [54] in a different context. Second, we consider the case of this equivalence for Newton's method and find that it is maintained in a rank-1 setting. Third, we reparameterize optimization steps to show that, under simplifying assumptions, the sets of stationary points are equivalent for a broader class of iterative optimization algorithms.Using gradient descent step according to (2) is equivalent to using two gradient steps of the alternating scheme (3), namely one step for (3a) and one step for (3b): Lemma 2 (Gradient Descent Step Equality between ( 2) and (3a)+(3b)). A gradient descent step according to (2) with arbitrary step size η coincides with two gradient descent steps, one according to (3a) and one according to (3b), where the optimization over θ has a step size of η and the optimization over z has a unit step size.Proof. Let θ ∈ Θ be the current parameter vector and let z = f (x; θ). Then the gradient descent steps according to (3a) and (3b) with step sizes 1 and η > 0 are expressed asEQUATIONEQUATIONCombining ( 16) and ( 17) eventually leads toEQUATIONwhich is exactly a gradient descent step of problem (1) starting at θ ∈ Θ with step size η.Here, we show that a corresponding equality also holds for a special case of the Newton step. Lemma 3 (Newton Step Equality between (2) and (3a)+(3b) for m = 1). In the case of m = 1 (i.e., a one-dimensional output), a Newton step according to (2) with arbitrary step size η coincides with two Newton steps, one according to (3a) and one according to (3b), where the optimization over θ has a step size of η and the optimization over z has a unit step size.Proof. Let θ ∈ Θ be the current parameter vector and let z = f (x; θ). Then applying Newton steps according to (3a) and (3b) leads toEQUATIONθ ← θ -η ∇ 2 θ 1 2 ∥z -f (x; θ)∥ 2 2 -1 ∇ θ 1 2 ∥z -f (x; θ)∥ 2 2 (20) = θ -η ∂ ∂θ ∂ f (x; θ) ∂ θ • (f (x; θ) -z) -1 ∂ f (x; θ) ∂ θ • (f (x; θ) -z) (21) = θ -η ∂ ∂θ ∂ f (x; θ) ∂ θ (f (x; θ)-z)+ ∂ f (x; θ) ∂ θ 2 -1 ∂ f (x; θ) ∂ θ • (f (x; θ) -z)Inserting (19) , we can rephrase the update above asEQUATIONBy applying the chain rule twice, we further obtain∇ 2 θ ℓ(f (x; θ)) = ∂ ∂θ ∂ f (x; θ) ∂ θ ∇ f ℓ(f (x; θ)) = ∂ ∂θ ∂ f (x; θ) ∂ θ ∇ f ℓ(f (x; θ)) + ∂ f (x; θ) ∂ θ ∂ ∂θ ∇ f ℓ(f (x; θ)) = ∂ ∂θ ∂ f (x; θ) ∂ θ ∇ f ℓ(f (x; θ)) + ∂ f (x; θ) ∂ θ ∇ f ∂ ∂θ ℓ(f (x; θ)) = ∂ ∂θ ∂ f (x; θ) ∂ θ ∇ f ℓ(f (x; θ)) + ∂ f (x; θ) ∂ θ 2 ∇ 2 f ℓ(f (x; θ)),which allows us to rewrite (22) asθ ′ = θ -(∇ 2 f ℓ(f (x; θ))) -1 ∇ 2 θ ℓ(f (x; θ)) -1 (∇ 2 f ℓ(f (x; θ))) -1 ∇ θ ℓ(f (x; θ)) = θ -(∇ 2 θ ℓ(f (x; θ))) -1 ∇ θ ℓ(f (x; θ)), which is exactly a single Newton step of problem (1) starting at θ ∈ Θ.In this section, we reparameterize update step of deterministic iterative optimization methods to show an equality of the sets of convergence points under model assumptions. Specifically, we can express optimization steps corr. to ( 2) and (3a) as θ ← arg minEQUATIONwhere Ω is a regularizer such that one step of a respective deterministic iterative optimization method (such as gradient descent) corresponds to the global optimum of the regularized optimization problems in (23) . More specifically, Ω 1 corresponds to the optimizer used in (2) and Ω 2 corresponds to the optimizer used in (3a). We note that the regularizer Ω is not a neural network regularizer like weight regularization or dropout; instead, the regularizers are reparameterizations of local optimization strategies. The reason for using the perspective of reparameterizing optimization steps as regularized cases of the global optimizations problems is exclusively for enabling a generality of our theory. We further note that a regularizer Ω is a valid choice if and only if it corresponds to a step of a respective deterministic iterative optimization method. As an example, the regularizer corresponding to a Newton step is Ω(z, f (x; θ)) = lf(x;θ) (z) -ℓ(z). Formally, given an input x we consider Y = {f (x; θ) , θ ∈ Θ} ⊆ R N ×m as the set of attainable outputs. Here, we assume Y = R N ×m , which generally holds for over-parameterized universal function approximators like neural networks, and we assume that such output remains attainable during training, which excludes a catastrophic collapse of the network, e.g., to a constant output. The regularizer Ω has the standard property that Ω(a, b) = 0 for any a = b. Note that the explicit form of the regularizer Ω does not need to be known nor computed. In SM E, we discuss explicit choices of Ω and provide further discussion.This allows us to express the set of points of convergence for the iterative optimization methods.Recall that an iterative optimization method has converged if it has arrived at a fixed point, i.e, the parameters do not change when applying an update. The set of points of convergence for (2) is as the method has converged if the update (3a) yields z ⋆ = z because the subsequent update (3b) will not change θ as z ⋆ = z = f (x; θ) already holds, and thus 1 2 ∥z ⋆ -f (x; θ)∥ 2 2 = 0. Now, we show that the iterative method (2) and the alternating method (3) lead to the same sets of convergence points, i.e., to the same sets possible solutions.A = θ | θ ∈ arg minFor illustrative purposes, in this section, we present examples of Newton losses of trivial losses, which we can solve in closed form.A popular loss function for classification is the softmax cross entropy (SMCE) loss, defined asEQUATIONwhere q i = exp(y i ) k j=1 exp(y j ).Example 1 (Softmax cross-entropy loss). For the SMCE loss, the induced Newton loss is given asEQUATION)where the element-wise Hessian variant isEQUATIONthe empirical Hessian variant is z ⋆ H = -E diag(q) -qq ⊤ -1 (q -p) + y ,and the empirical Fisher variant is z ⋆ F = -E (q -p)(q -p) ⊤ -1 (q -p) + y .A less trivial example is the binary cross-entropy (BCE) loss with ℓ BCE (y) = BCE(y, p) where p ∈ ∆ m is a probability vector encoding of the ground truth.Example 2 (Binary cross-entropy loss). For the BCE loss, the induced Newton loss is given asEQUATIONwhere the element-wise Hessian variant isEQUATIONthe empirical Hessian variant isEQUATIONand the empirical Fisher variant isEQUATIONand where ⊙ and ⊘ are element-wise operations.The BCE loss is often extended using the logistic sigmoid function to what is called the sigmoid binary cross-entropy loss (SBCE), defined as ℓ SBCE (y) = BCE(σ(y), p)where σ(x) = 1 1 + exp (-x) .Example 3 (Sigmoid Binary cross-entropy loss). For the SBCE loss, the induced Newton loss is given asEQUATIONwhere the element-wise Hessian variant is z ⋆ E = -diag σ(y) -σ(y) 2 -1 (σ(y) -p) + y ,the empirical Hessian variant is z ⋆ H = -diag E y σ(y) -σ(y) 2 -1 (σ(y) -p) + y ,and the empirical Fisher variant is z ⋆ F = -E y (σ(y) -p)(σ(y) -p) ⊤ -1 (σ(y) -p) + y ,and where ⊙ and ⊘ are element-wise operations.For the empirical Fisher method, where z ⋆ ∈ R N ×m is computed viaz ⋆ i = ȳi -1 N N j=1 ∇ ȳj ℓ(ȳ) ∇ ȳj ℓ(ȳ) ⊤ + λ • I -1 ∇ ȳi ℓ(ȳ)for all i ∈ {1, ..., N } and ȳ = f (x; θ), we can simplify the computation via the Woodbury matrix identity [58] . In particular, we can simplify the computation of the inverse to EQUATIONEQUATIONEQUATIONThis reduces the cost of matrix inversion from O(m 3 ) down to O(N 3 ). We remark that this variant is only helpful for cases where the batch size N is smaller than the output dimensionality m.Further, we remark that in a case where the output dimensionality m and the batch size N are both very large, i.e., the cost of the inverse is typically still small in comparison to the computational cost of the neural network. To illustrate this, the dimensionality of the last hidden space (i.e., before the last layer) h l is typically larger than the output dimensionality m. Accordingly, the cost of only the last layer itself is O(N • m • h l ), which is typically larger than the cost of the inversion in Newton Losses. In particular, assuming h l > m, O(N • m • h l ) > O(min(m, N ) 3 ).We remark that the Woodbury inverted variation can also be used in the InjectFisher implementation.For each step, z is initialized to the current model predictions f (x; θ). 


Decision: Given my own evaluation of this work and taking into account the reviewer's expertise, I am recommending to accept this interesting submission.

Reviews: review: paper_summary: paper_summary: The paper proposes a new method to optimize complex possibly non-smooth and algorithmic losses for neural networks. The approach is based on splitting the problem into two-step procedure, where in the first step we construct and optimize the so-called Newton loss and the second step is based on SGD-type procedure for MSE loss with the first step. The authors present a wide experimental comparison of the proposed Fisher and Newton approaches with existing methods.
strength_weakness: strength_weakness: The paper does not contain any proofs or convergence guarantees. The mathematical formulation of the main problem is also quite confusing for me.  For example, is vector $x$ fixed for all steps or is it a batch of data?  Is it a sum-type problem or an expectation problem? What are the properties for $l(\cdot)$? Is it differentiable, smooth? Because some parts of the text said that the loss is non-smooth and later we calculate the Hessian of such a function.  In Formulas 1 and 2, it is not clear what are the fixed parameters or data. Should $\theta$ in 2a be $\theta_{t-1}$? Also, I think the mention of Lemma 2 in the main text could be very helpful.

For the experimental section, personally, it feels that the most of space is taken by the description of the problems and the setup and not the actual comparison. As the paper is mostly experimental and empirical, one would expect a better comparison of the proposed methods with the multiple benchmarks. There are no convergence figures with the per-iteration or per-gradient performance. As the authors claim, the main issues in the existing approaches are vanishing and exploding gradients. However, I didn’t find any clipping method for the comparison, which are the possible solutions for exploding gradients.
questions: questions: See Weaknesses
limitations: limitations: n/a
paper_summary: paper_summary: The paper proposes second-order optimization with splitting for hard objectives that arise as smoothing of such hard problems as sorting and ranking to address the problem of vanishing/exploding gradients.
strength_weakness: strength_weakness: 1. Insufficient experiments. I'd appreciate adding a comparison here with the SFA technique from there, as it will rely only on first-order information: https://arxiv.org/pdf/2003.02122
questions: questions: 1. Considering Newton's method in the case of non-convex objectives is a mistake. No matter how much it is regularized, as long as regularization is L2, had authors considered cubic regularization? E.g., https://link.springer.com/article/10.1007/s10107-006-0706-8
2. Adding to the first question, for ranking objectives, Hessians are expected to converge to zero. Have you considered an increasing learning rate schedule? I am somewhat sure that this hessian/fisher type of method, due to vanishing gradients, also vanishes, resulting in effectively increasing the learning rate up to $\propto \lambda^{-1}$. I will appreciate experiments against the first-order method with the learning rate scheduler growing up to that value.
limitations: limitations: Yes
paper_summary: paper_summary: The paper presents an alternative backpropagation scheme for deep learning with algorithmic losses that combines a preconditioned step on the loss with a gradient step on a least square objective. Two preconditioning methods are investigated: using the Hessian, or the empirical Fisher. Experiments demonstrate that the proposed plugin consistently improve performance of algorithms across architectures and losses on two benchmarks. An ablation study of the potential additional tikhonov regularization is give, as well as a discussion of runtime comparisons.
strength_weakness: strength_weakness: - The soundness of the approach from a theoretical viewpoint is lacking. However, it is probably better to have clear experimental evaluations than wobbly theoretical explanations. And theoretical explanations can be given later.
questions: questions: - Why do the authors multiply an inverse of the Hessian with the gradients? Such operation amounts to solve a linear system and is always best handled with matrix-free linear algebra solvers that can make use of Hessian vector products.
- Can the approach be cast as appropriate block-coordinate minimization of the Lagrangian associated to the composition of the loss and the network? As the authors acknowledge if both steps in 2.a, 2.b are gradient steps, we retrieve the classical gradient back-propagation as demonstrated in [57]. Such a setting may help uncover the principles underlying the improved performance.
- I don't see how the learning rates were tuned in the experiments. Since it is generally a crucial hyperparameter for optimization methods, can the authors comment on this? 
- The authors keep on mentioning vanishing/exploding gradient phenomena. Could the authors actually demonstrate this issue numerically? I suppose the vanishing/exploding gradient phenomena depends on the number of steps taken by the algorithm defining the loss. Maybe one could correlate the number of steps taken by the underlying algorithm and the improved performance of the proposed approach.
limitations: limitations: Limitations (or rather scope) of the approach are discussed in Appendix A.

Comments: title: Acknowledging rebuttal
Thank you for your response and increasing your score.
comment: Dear Reviewer nJzH,

thank you for acknowledging our rebuttal and your positive remarks.

Thanks for the reference to [1], we will add a discussion to the camera-ready.
Also, thanks for the "second (uncentered) moments" suggestion; we will adjust the paper accordingly.

---
(Meta: Regarding your "Though I would not argue for a strong accept", we just wanted to make sure that you are aware that between the currently selected "weak accept" and a "strong accept", there is also the score "accept" this year.)
I thank the authors for their answers to my questions and comments.

- Surely, the paper would benefit from a better theoretical understanding.
- That said, the paper presents an extensive set of experiments and clear improvements with the proposed method. These experimental results may be a base for a better understanding of the approach from a mathematical viewpoint. Numerous experiments diagnose the approach to build upon it. Though I would not argue for a strong accept, I believe the value of this experimental work should not be completely dismissed and I maintain my score.
- Regarding the regularization comments of reviewer tok2, some recent works have analyzed the use of quadratic regularizations of Newton method with proven convergence guarantees (see [1] and further references of [1]). Most importantly, cubic regularizations can be expensive to compute and so not practical. Given that this paper is interested in practical improvements, analyzing the performance with either empirical Fisher* or Newton's method with quadratic regularization is valuable empirically.  

[1] Super-Universal Regularized Newton Method, https://arxiv.org/pdf/2208.05888v1
* To avoid further confusion about the Fisher matrix and its "empirical part", I would suggest the authors to rather speak about the second (uncentered) moments of the gradients.
Dear Reviewer tok2,

thank you very much for responding to our rebuttal, and for your clarification. 

> That is the point that you'll need to overregularize offsetting entirely your negative spectrum, which happens to be as large as positive one, as smoothed objectives are necessarily hyperbolic in nature. So while "you can" it's absolutely meaningless type of regularisation in this case.

We would like to kindly to point at the fact that the empirical Fisher is always PSD, which resolves this concern for the Newton Loss (Fisher).
For the Hessian case, you are right that it could require overregularizing, but the empirics show that using the Hessian still works very well, in most cases even better than the Fisher variant.

> [...] resulting in effectively increasing the learning rate up to $\propto \lambda^{-1}$. I will appreciate experiments against the first-order method with the learning rate scheduler growing up to that value.
> 
> [...] about the fact that it actually will have an effect if you make learning rate scheduled

We apologize as we might have misunderstood you earlier.
To address this concern with concrete experiments, we have now run a set of experiments with learning rate scheduling.
In particular, we schedule the learning rate, as suggested to grow the learning rate to initial learning rate $\cdot\lambda^{-1}$.
We use both linear as well as eponential growth for the learning rate. Further, we use NeuralSort (n=5), and ran experiments for both $\lambda=0.1$ and $\lambda=0.01$, covering the lambdas of both Hessian Newton Loss and Fisher Newton Loss.
We ran each experiment for 5 seeds.

| Method | Full Rankings Acc. | (Individual Ranks Acc.) |
|--|--|--|
| Baseline                                                   | 71.33±2.05 | 87.10±0.96 |
| Hessian Newton Loss (ours)                                 | 83.31±1.70 | 92.54±0.73 |
| Fisher Newton Loss (ours)                                  | 83.93±0.62 | 92.80±0.30 |
| GrowLR $\propto\lambda^{-1}$ ($\lambda=0.1$, linear)       | 69.99±3.17 | 86.39±1.58 |
| GrowLR $\propto\lambda^{-1}$ ($\lambda=0.1$, exponential)  | 70.66±1.21 | 86.73±0.58 |
| GrowLR $\propto\lambda^{-1}$ ($\lambda=0.01$, linear)      | 68.44±3.40 | 85.68±1.64 |
| GrowLR $\propto\lambda^{-1}$ ($\lambda=0.01$, exponential) | 68.19±1.73 | 85.48±0.86 |

This shows that introducing a learning rate scheduler as suggested does not improve performance, and that the method is not "resulting in effectively increasing the learning rate up to $\propto \lambda^{-1}$".

In case it is relevant to you, we also extend the table by loss factors (multiplying the loss by 0.1 or 10, and thus effectively multiplying the gradient by the same factor), which shows that loss factors do not have any significant effect.

| Method | Full Rankings Acc. | (Individual Ranks Acc.) |
|--|--|--|
| Loss factor 0.1                                            | 72.04±1.24 | 87.43±0.61 |
| Baseline (loss factor 1.0)                                 | 71.33±2.05 | 87.10±0.96 |
| Loss factor 10.0                                           | 71.43±1.25 | 87.17±0.60 |
>We thank the reviewer for highlighting Nesterov's work on cubic regularization. As noted in the mentioned paper, various regularization techniques can be applied when using the Newton method for non-convex or ill-conditioned problems. The most classical technique is Tikhonov regularization (or Levenberg-Marquardt regularization). Given that our Tikhonov-regularized version already demonstrates good empirical performance, we did not explore more sophisticated regularization techniques such as cubic regularization. However, this is indeed an interesting direction for future research. 

That is the point that you'll need to overregularize offsetting entirely your negative spectrum, which happens to be as large as positive one, as smoothed objectives are necessarily hyperbolic in nature. So while "you can" it's absolutely meaningless type of regularisation in this case.

>As the optimizer for all experiments in this paper is the Adam optimizer, such effects can be rather excluded, as multiplying the loss by a factor does not affect updates (apart from numerical precision effects).

It's about loss <-> lr scaling, not about the fact that it actually will have an effect if you make learning rate scheduled

I'll keep my score as is.
Dear Reviewer qoSU,

Thank you very much for the clarification! Accordingly, we will include both per-epoch and per-time plots in the appendix of the camera-ready version.

Please let us know in case any further questions or concerns come up.
Dear Authors,

Thank you! Regarding the “dependence on gradient/Hessian-vector products”, yes, I meant it as an alternative for the time plots, as some authors prefer this type of plots over the time one. I am not sure what is the best choice for your paper, that’s why, I listed both of them. 
Thank you!
Dear Reviewer qoSU,

Thank you very much for responding to our rebuttal and for increasing your score.

To address your clarification wrt. the plot request, we offer to include the same as Figure 3, but for all problems in the appendix in the camera-ready.
We can also make additional versions with time on the x-axis. Regarding plotting "dependence on gradient/Hessian-vector products", could you clarify what you mean such that we can also include it in the appendix? Was "dependence on gradient/Hessian-vector products" meant as an equivalent to the "time performance" or was it, e.g., number of floating-point operations?

Best regards,

The Authors
Dear Authors,

thank you for the detailed answers and comments.

*Could you please clarify what precisely you would like us to plot, considering that we show a per-epoch convergence plot in Figure 3?*

The same as Figure 3, but for all problems. They can be presented in the Appendix for the lack of space. Also, the dependence on gradient/Hessian-vector products or time performance could be interesting to have.

Regarding the rebuttals, I would like to increase my score from 4 to 5.