Second-Order Forward-Mode Automatic Differentiation for Optimization

## Abstract

Abstract This paper introduces a second-order hyperplane search, a novel optimization step that generalizes a second-order line search from a line to a k 𝑘 k italic_k -dimensional hyperplane.
This, combined with the forward-mode stochastic gradient method,
yields a second-order optimization algorithm that consists of forward passes only, completely avoiding the storage overhead of backpropagation.
Unlike recent work that relies on directional derivatives (or Jacobian–Vector Products, JVPs), we use hyper-dual numbers to jointly evaluate both directional derivatives and their second-order quadratic terms. As a result, we introduce forward-mode weight perturbation with Hessian information (FoMoH). We then use FoMoH to develop a novel generalization of line search by extending it to a hyperplane search.
We illustrate the utility of this extension and how it might be used to overcome some of the recent challenges of optimizing machine learning models without backpropagation. Our code is open-sourced at https://github.com/SRI-CSL/fomoh .

## 1 Introduction

There is a growing interest in investigating the practical plausibility of forward-mode automatic differentiation (AD) as an alternative to
reverse-mode AD (aka backpropagation) for neural network optimization.
Forward gradient descent (FGD) [] relies on sampling tangent vectors to update function parameters. These parameters are updated by subtracting the tangents that are scaled by their directional derivatives. As a result, their approach only requires forward passes, and avoids the computation and memory costs associated with implementing the backwards pass of reverse-mode AD. While there is still an accuracy gap between forward-mode and reverse-mode optimization approaches [], there are several recent efforts that have explored alternative ways of leveraging forward-mode AD, with a focus of reducing the variance of the gradient estimator with the increase in dimensions [, ].
Thus, the large performance gap between the forward-mode approaches and backpropagation (BP) has been shrinking.

An interesting direction to explore, which builds on the existing work in forward-mode AD for optimization, is the introduction of Hessian information.
Second-order derivative information provides optimization routines with information about the local curvature. These methods often require access to the Hessian, whose storage is not feasible for high-dimensional problems.
As an example, a function, $f:\mathbb{R}^{D}\rightarrow\mathbb{R}$, has a quadratic storage cost of $\mathcal{O}(D^{2})$ and a linear compute cost of $\mathcal{O}(D)$ reverse-mode gradient evaluations to build the Hessian [, ].
Instead of building the full Hessian through the use of reverse-mode AD, we introduce an approach incorporating Hessian information, known as FoMoH, which significantly outperforms the first-order gradient descent method, FGD. Furthermore, we demonstrate how our approach effectively bridges the gap between a line search and a full Newton step. By incorporating Hessian information, FoMoH leverages second-order derivatives to provide more accurate curvature approximations of the objective function. This allows the method to adaptively adjust the step size and direction with greater precision compared to a simple line search, which only uses gradient information. Simultaneously, FoMoH avoids the computational complexity of a full Newton step, which requires calculation of $(\nabla^{2}f)^{-1}\mathbf{v}$ where $\nabla^{2}f$ is the Hessian. As a result, our approach balances the efficiency of a line search with the accuracy of a full Newton method, offering a robust and versatile optimization technique.

The central contributions of this paper are as follows:

- •
We introduce three new optimization approaches that use forward-mode second-order information.
- •
We show how one of these approaches, Forward-Mode Hyperplane Search, generalizes from a line search all the way to Newton’s method, without any need for backpropagation.
- •
We demonstrate the performance of the proposed new approaches on optimization problems of different parameter sizes and difficulty to show the advantage of second-order information.
- •
We release an AD backend in PyTorch that implements nested forward AD and interfaces with PyTorch models: [https://github.com/SRI-CSL/fomoh](https://github.com/SRI-CSL/fomoh).

The rest of the paper is organized as follows. §2 and §3 summarize the related work and relevant preliminary information regarding AD, which is then followed by §4 that outlines what is required to extend forward-mode AD to higher order derivatives. §5 introduces FoMoH and the generalization to our forward-mode hyperplane search. Finally, §6 provides experimental results that explore the behavior of FoMoH. We then conclude in §7.

## 2 Related Work

There has been considerable interest in developing approaches that avoid reverse-mode AD and its backwards pass. In moving away from backpropagation, it might be possible to build optimization algorithms that more closely align with biological systems [, ], or enable neural networks to run on emerging hardware architectures, such as analog optical systems []. Baydin et al. [] introduced FGD as a possible replacement for backpropagation in neural network training by relying on weight perturbations. This removed the truncation error of previous weight perturbation approaches [] by moving to forward-mode AD. While FGD is a promising backpropagation-free approach,
the scaling of FGD to high-dimensional models is challenging due the variance of the gradient estimator.
This challenge has led to multiple efforts that have focused on reducing this variance [, , ]. In particular, a common approach has been to rely on local backpropagation steps to reduce the variance and/or to provide a better guess for the perturbation direction. In our work, we focus on second-order forward-mode optimization approaches, but highlight that these other approaches on variance reduction are orthogonal to our approach and could also be combined together and generalized to FoMoH.

Becker et al. [] were one of the first to the use the second-order information to optimize neural networks. This approach leverages a local backpropagation step [] to capture the diagonal Hessian terms resulting in what they call a “pseudo-Newton step”. Both LeCun et al. [, §6–§9] and Goodfellow et al. [, §8.6] provide discussions of the use of Hessian information for neural network training and cover Newton’s method, as well as the Levenberg–Marquardt algorithm [, ], conjugate gradients, and BFGS. A key objective of these approaches is to investigate whether second-order information can be leveraged for neural network training without the prohibitively high cost of a full Hessian evaluation. Additionally, the research question then explored is whether the approximated Hessian is still good enough to give the desired advantage over first-order methods. Examples of effective approaches often rely on diagonal preconditioners, with some approximating or directly using the diagonal of the Hessian [, , ], and others leveraging momentum and a variant of the diagonal of the empirical Fisher [, , ]. While extending a gradient preconditioner to the full inverse Hessian is generally infeasible, there are approaches that use: block-diagonal approximations, low-rank approximations, and Krylov-subspace based approximations [, ]. Additional references for these preconditioners can be found in Martens []. Finally, although it is challenging to scale second-order approaches to large neural network models, a recent work, Sophia [], has managed to show success in using such an approach. Like previous works, they use a diagonal approximation of the Hessian. Importantly, they show that Sophia requires half the number of update steps as Adam [] to train a language model (GPT-2 []). It is worth noting that none of the described second-order optimization methods rely solely on forward-mode AD like the one being proposed in this paper.

## 3 Automatic Differentiation

In this section we summarize the common definitions of forward-mode and reverse-mode automatic differentiation. For a complete introduction, please refer to Griewank and Walther []. We use $\mathbf{x}\in\mathbb{R}^{D}$ to denote a column vector.

Forward-Mode AD [] applies the chain rule in the forward direction. The forward-mode evaluation, $F(\bm{\theta},\mathbf{v})$, requires an additional tangent vector, $\mathbf{v}\in\mathbb{R}^{D}$, along with the parameter vector $\bm{\theta}\in\mathbb{R}^{D}$ for a function $\boldsymbol{f}:\mathbb{R}^{D}\rightarrow\mathbb{R}^{O}$. The result of the evaluation, $[\boldsymbol{f}(\bm{\theta}),\nabla\boldsymbol{f}(\bm{\theta})\mathbf{v}]$, is a function evaluation and the corresponding Jacobian vector product (JVP), where $\nabla\boldsymbol{f}(\bm{\theta})\in\mathbb{R}^{O\times D}$. For a unidimensional output function, the JVP is the directional derivative, $\nabla\boldsymbol{f}(\bm{\theta})\cdot\mathbf{v}$. The time and space (memory cost) complexity
are linear, both approximately twice that of a single function evaluation.(^1^11More precisely, the basic time complexity of a forward evaluation is constant $\in[2,2.5]$ times that of the function call [].) A common implementation of forward-mode AD is to use dual numbers. A dual number $a+b\epsilon\in\mathbb{D}(\mathbb{R})$ contains a real (primal) component, $a\in\mathbb{R}$, and a dual component, $b\in\mathbb{R}$. We can think of this as representing a truncated Taylor series, $a+b\epsilon+\mathcal{O}(\epsilon^{2})$, notationally simplified by the rule $\epsilon^{2}=0$. Using this, $f(a+b\epsilon)=f(a)+\nabla f(a)b\epsilon$. A simple example can be shown for the function, $f(a_{1},a_{2})=a_{1}\times a_{2}$. Using dual numbers, $(a_{1}+b_{1}\epsilon)(a_{2}+b_{2}\epsilon)$, we retrieve the function evaluation and the corresponding well-known product rule: $a_{1}a_{2}+(a_{1}b_{2}+b_{1}a_{2})\epsilon$. This can be extended to multiple dimensions, and is the basis of forward-mode AD: lift all real numbers $\mathbb{R}$ to dual numbers $\mathbb{D}(\mathbb{R})$.

Reverse-Mode AD requires both a forward pass and a reverse pass. The reverse-mode evaluation $R(\bm{\theta},\boldsymbol{u})$, also requires the additional adjoint vector, $\boldsymbol{u}\in\mathbb{R}^{O}$, which is often set to $1$ for scalar-valued functions. Using the same notation as for forward-mode, an evaluation of reverse mode results in the vector-Jacobian product, $\boldsymbol{u}^{\top}\nabla\boldsymbol{f}(\bm{\theta})$, as well as the function evaluation. When $\boldsymbol{u}=1$, this results in the gradient $\nabla\boldsymbol{f}(\bm{\theta})$. Reverse-mode is required to store intermediate values on the forward pass that are used during the reverse pass. This results in a higher time and space complexity, that is higher computational cost and memory footprint per call of $R(\bm{\theta},\boldsymbol{u})$. However, for the scalar-valued functions ($O=1$) that are common for ML optimization problems, only a single call of $R(\bm{\theta},\boldsymbol{u})$ is needed to collect all the gradients compared to $D$ (dimension of inputs) calls of $F(\bm{\theta},\mathbf{v})$. This is one of the key reasons for the widespread adoption of the reverse-mode AD in current gradient-based machine learning methods, despite the higher memory footprint.

## 4 Higher-Order Forward Mode Automatic Differentiation

As described in the previous section, forward-mode automatic differentiation implementations can
use dual numbers, $\bm{\theta}+\mathbf{v}\epsilon$. Dual numbers can be extended to truncate at a higher order, or to not truncate at all; and to allow nesting by supporting multiple distinct formal $\epsilon$ variables []. Specifically focusing on second-order terms is often referred to as hyper-dual numbers (for example in the Aeronautics and Astronautics community []). A hyper-dual number is made up from four components, which is written as $\bm{\theta}+\mathbf{v}_{1}\epsilon_{1}+\mathbf{v}_{2}\epsilon_{2}+\mathbf{v}_{
12}\epsilon_{1}\epsilon_{2}$. In the same manner that we look at imaginary and real parts of complex numbers, we can look at the first derivative parts of a hyper-dual number by inspecting the $\epsilon_{1}$ and $\epsilon_{2}$ components, and we can look at the second derivative by inspecting the $\epsilon_{1}\epsilon_{2}$ component. To understand how this formulation arises, we can introduce the definitions $\epsilon_{1}^{2}=\epsilon_{2}^{2}=(\epsilon_{1}\epsilon_{2})^{2}=0$ and replace the Taylor series expansion of a function, $f:\mathbb{R}^{D}\rightarrow\mathbb{R}$ around $\bm{\theta}$ with perturbation, $\mathbf{d}\in\mathbb{R}^{D}$, with an evaluation of a hyper-dual number:(^2^22Note, we have left our function definition as being a scalar output for pedagogical reasons but nothing precludes a vector or matrix output, which is required for the composition of functions in most ML architectures.)

$$ $\displaystyle f(\bm{\theta}+\mathbf{d})$ $\displaystyle=f(\bm{\theta})+\nabla f(\bm{\theta})\mathbf{d}+\frac{1}{2} \mathbf{d}^{\top}\nabla^{2}f(\bm{\theta})\mathbf{d}+\cdots$ $\displaystyle f(\bm{\theta}+\mathbf{v}_{1}\epsilon_{1}+\mathbf{v}_{2}\epsilon_ {2}+\mathbf{v}_{12}\epsilon_{1}\epsilon_{2})$ $\displaystyle=f(\bm{\theta})+\nabla f(\bm{\theta})\mathbf{v}_{1}\epsilon_{1}+ \nabla f(\bm{\theta})\mathbf{v}_{2}\epsilon_{2}$ $\displaystyle\quad+\nabla f(\bm{\theta})\mathbf{v}_{12}\epsilon_{1}\epsilon_{2 }+\mathbf{v}_{1}^{\top}\nabla^{2}f(\bm{\theta})\mathbf{v}_{2}\epsilon_{1} \epsilon_{2}+\cdots.$ $$

An alternative but isomorphic view is to regard $a+b\epsilon_{1}+c\epsilon_{2}+d\epsilon_{1}\epsilon_{2}=(a+b\epsilon_{1})+(c+d
\epsilon_{1})\epsilon_{2}$ as an element of $\mathbb{D}(\mathbb{D}(\mathbb{R}))$, with subscripts to distinguish the inner vs outer $\mathbb{D}$s; from an implementation perspective, hyperduals can be regarded as inlining the nested structures into a single flat structure.

### 4.1 Implications of Hyper-Dual Numbers for Machine Learning

A function evaluation with hyper-dual numbers takes an input vector, $\bm{\theta}+\mathbf{v}_{1}\epsilon_{1}+\mathbf{v}_{2}\epsilon_{2}+\mathbf{0}
\epsilon_{1}\epsilon_{2}$, with the $\epsilon_{1}\epsilon_{2}$ part set to zero. A typical setting to get exact gradient and Hessian elements is to set $\mathbf{v}_{1}=\mathbf{e}_{i}$ and $\mathbf{v}_{2}=\mathbf{e}_{j}$, where $\mathbf{e}_{i}$ and $\mathbf{e}_{j}$ are each a basis of one-hot unit vectors. Therefore, these basis vectors select the corresponding elements of the gradient and Hessian:

$$ $f(\bm{\theta})+\nabla f(\bm{\theta})_{i}\epsilon_{1}+\nabla f(\bm{\theta})_{j} \epsilon_{2}+\nabla^{2}f(\bm{\theta})_{ij}\epsilon_{1}\epsilon_{2}=f(\bm{ \theta}+\mathbf{e}_{i}\epsilon_{1}+\mathbf{e}_{j}\epsilon_{2}+\mathbf{0} \epsilon_{1}\epsilon_{2}).$ $$

A single loop over the input dimension provides the exact gradient, whereas a nested loop provides the full Hessian. As a side note, a single loop also can give the Hessian vector product. This is done by setting one of the tangent vectors to the chosen vector, and looping through the basis for the other tangent vector.(^3^33While interesting, these results might not immediately seem attractive to the ML community. The Hessian calculation of a model with parameters, $\bm{\theta}\in\mathbb{R}^{D}$, would require $D\cdot(D+1)/2$ function evaluations, whereas Hessian vector products are widely available in many AD libraries leveraging tricks such as a forward over reverse routine [, ].) However, the key advantage is that with a single forward pass, we get both first order and second order information in the form of directional derivatives and the local curvature information respectively. We have already seen from Baydin et al. [] and follow up works [, ] that despite the scalar form of the directional derivative, it can still be used to build optimization routines that appear competitive with backpropagation.
In this paper, we investigate whether the additional access to local curvature through forward-mode AD can enable improvements over the current FGD algorithm.

#### Local Curvature: 𝐯 1 ⊤ ⁢ ∇ 2 f ⁢ ( 𝜽 ) ⁢ 𝐯 2 superscript subscript 𝐯 1 top superscript ∇ 2 𝑓 𝜽 subscript 𝐯 2 \mathbf{v}_{1}^{\top}\nabla^{2}f(\bm{\theta})\mathbf{v}_{2} bold_v start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ⊤ end_POSTSUPERSCRIPT ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_f ( bold_italic_θ ) bold_v start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT .

The Hessian contains the curvature information at a point, $\bm{\theta}$, in the form of the second order partial derivatives. When we evaluate a function over hyper-dual numbers we get the bilinear form, $\mathbf{v}_{1}^{\top}\nabla^{2}f(\bm{\theta})\mathbf{v}_{2}$. This is a function that operates over two vectors, such that it is linear in each vector separately. The bilinear form
is common in optimization routines, such as for conjugate gradient descent to assess whether two vectors are conjugate with respect to the Hessian. The value of $\mathbf{v}_{1}^{\top}\nabla^{2}f(\bm{\theta})\mathbf{v}_{2}$ tells us about how curvature covaries along the two vectors. In the case where $\mathbf{v}_{1}=\mathbf{v}_{2}=\mathbf{v}$, we arrive at the quadratic form that describes the curvature, which indicates the rate of change of the slope as you move in the direction of $\mathbf{v}$. The quadratic form also provides information on the convexity of the function at that point. If $\mathbf{v}^{\top}\nabla^{2}f(\bm{\theta})\mathbf{v}>0,\forall\mathbf{v}\in
\mathbb{R}^{D}_{>0}$, then the function is convex at that point and the Hessian is positive definite. The curvature also indicates the sensitivity of moving in certain directions. For example, when using gradients to optimize a function, taking a gradient step in a region of low curvature is likely to increase the value of the function. However the same sized step in a region of large curvature could significantly change the value of the function (for better or worse).

#### Computational Cost.

The cost of a single forward pass with hyper-dual numbers is of the same order of time and space complexity
as the original function call. This reduces the memory cost compared to reverse-mode AD. The forward pass with hyper-dual numbers scales in the same way as the original function call to the number of parameters. However, there is a constant overhead price (no change in computational complexity) to be paid in the form of evaluating the second order terms throughout a forward pass. An example is that the addition of two scalars now requires $4$ additions, and multiplication now requires $9$ products and $5$ additions. However, unlike for reverse-mode, these intermediate values can be overwritten once they have been used to propagate gradient information.

## 5 Forward-Mode Optimization with Second Order Information

We now introduce the three new optimization routines that incorporate forward-mode second order information:

- 1.
Forward-Mode Line Search (FoMoH)
- 2.
Forward-Mode Line Search with Backpropagation (FoMoH-BP)
- 3.
Forward-Mode Hyperplane Search (FoMoH-$K$D)

### 5.1 Forward-Mode Line Search: FoMoH

We introduce a new gradient-based optimization routine leveraging Forward-Mode automatic differentiation with Hessian information (FoMoH). We begin with the one dimensional case, which uses a second-order line search [, ]. For a given update direction $\mathbf{v}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$ and learning rate $\eta$, FoMoH normalizes all gradient steps by the curvature, giving:

$$ $\bm{\theta}^{\prime}=\bm{\theta}+\eta\frac{\nabla f(\bm{\theta})\cdot\mathbf{v }}{|\mathbf{v}^{\top}\nabla^{2}f(\bm{\theta})\mathbf{v}|}\mathbf{v}.$ (1) $$

The directional derivative in the numerator and the curvature (the quadratic form) in the denominator are both provided with one forward-pass of the function $f(\bm{\theta}+\mathbf{v}\epsilon_{1}+\mathbf{v}\epsilon_{2}+\mathbf{0}\epsilon
_{1}\epsilon_{2})$. The normalization via this quadratic form results in accounting for the unit distance at the location $\mathbf{x}$ in the direction $\mathbf{v}$. Therefore, this update step takes into account the distance metric defined at $\bm{\theta}$. In regions of high curvature the step size will be smaller, which is a desirable behavior. Like Newton’s method, setting $\eta\approx 1.0$ seems to work in many cases, but by analogy with trust regions we suggest the inclusion of a learning rate, although we found the method to be reasonably robust to its value.

### 5.2 Forward-Mode Line Search with Backpropagation: FoMoH-BP

A line search starts with identifying a descent direction, followed by determining the value of the step size to move in that direction. Therefore, the second approach that we propose is to perform a line search on the gradient. Thus, this combines forward-mode and reverse-mode to build an optimization routine that provides the step-size for the ground truth gradient obtained from backpropagation. This additional step sets $\mathbf{v}=\nabla f(\bm{\theta})$ in Equation (1). The result is an update step:

$$ $\bm{\theta}^{\prime}=\bm{\theta}+\eta\frac{\nabla f(\bm{\theta})\cdot\nabla f( \bm{\theta})}{|\nabla f(\bm{\theta})^{\top}\nabla^{2}f(\bm{\theta})\nabla f( \bm{\theta})|}\nabla f(\bm{\theta}).$ (2) $$

Unlike FoMoH (and FoMoH-$K$D in the next section), FoMoH-BP includes a single reverse-mode AD evaluation. Specifically, this requires a single backpropagation step followed by forward-mode step that sets the tangent vectors to the gradient.

### 5.3 Forward-Mode Hyperplane Search: FoMoH- K 𝐾 K italic_K D

The final algorithm that we propose is the Forward-Mode $K$-Dimensional Hyperplane Search, FoMoH-$K$D. Rather than performing a second order line search along direction $\mathbf{v}$, we perform a $K$-dimensional hyperplane search. Starting with the $K=2$ example, if we take two search directions, $\mathbf{v}_{1}$ and $\mathbf{v}_{2}$, we can build a $2\times 2$ matrix to form a Hessian in the plane defined by $\bm{\theta}+\kappa_{1}\mathbf{v}_{1}+\kappa_{2}\mathbf{v}_{2}$. We evaluate a function, $f(\cdot)$, with a hyper-dual number using the pairs $\{\mathbf{v}_{1},\mathbf{v}_{1}\}$, $\{\mathbf{v}_{1},\mathbf{v}_{2}\}$, and $\{\mathbf{v}_{2},\mathbf{v}_{2}\}$ for the $\epsilon_{1},\epsilon_{2}$ coefficients. The result is the Hessian, $\mathbf{\tilde{H}}_{2\times 2}$, in the $2\times 2$ plane, and corresponding step sizes, $\kappa_{1}$ and $\kappa_{2}$, to take in each search direction:

$$ $\displaystyle\mathbf{\tilde{H}}_{2\times 2}$ $\displaystyle=\left[\begin{array}[]{cc}\mathbf{v}_{1}^{\top}\nabla^{2}f(\bm{ \theta})\mathbf{v}_{1}&\mathbf{v}_{1}^{\top}\nabla^{2}f(\bm{\theta})\mathbf{v} _{2}\\ \mathbf{v}_{2}^{\top}\nabla^{2}f(\bm{\theta})\mathbf{v}_{1}&\mathbf{v}_{2}^{ \top}\nabla^{2}f(\bm{\theta})\mathbf{v}_{2}\end{array}\right],$ $\displaystyle\left[\begin{array}[]{c}\kappa_{1}\\ \kappa_{2}\end{array}\right]$ $\displaystyle=\mathbf{\tilde{H}}_{2\times 2}^{-1}\left[\begin{array}[]{c} \mathbf{v}_{1}^{\top}\nabla f(\bm{\theta})\\ \mathbf{v}_{2}^{\top}\nabla f(\bm{\theta})\end{array}\right]$ (9) $$

As a result, we formulate a new update step,

$$ $\bm{\theta}^{\prime}=\bm{\theta}+\kappa_{1}\mathbf{v}_{1}+\kappa_{2}\mathbf{v} _{2}.$ $$

We then extend the above result to any $K$-dimensional hyperplane by sampling $K$ search directions and evaluating the corresponding update step:

$$ $\bm{\theta}^{\prime}=\bm{\theta}+\eta\sum_{k=1}^{K}\kappa_{k}\mathbf{v}_{k}.$ (10) $$

This resulting generalized hyperplane update step allows one to trade-off computational cost with the size of the search space. For example, the cost of evaluating a $K$-dimensional Hessian and then invert it is $\mathcal{O}(K^{3})$, which is feasible for small enough $K$.(^4^44For instances where $\mathbf{\tilde{H}}_{K\times K}$ is not invertible, we add jitter to the diagonal. This seems to work well.) Our new forward-mode hyperplane search, FoMoH-$K$D, opens up the possibility of transitioning between a line search, when $K=1$, all the way to a full Newton step, when $K=D$, which we demonstrate in §6.1. The pseudo-code for a single update step is given in Algorithm 1. The overall FoMoH-$K$D routine is given in Algorithm 2.

Figure: Algorithm 1 FoMoH-$K$D, hyperplane update step for function, $f$, and parameters, $\bm{\theta}\in\mathbb{R}^{D}$.

## 6 Experiments

### 6.1 Rosenbrock Function

In this section, we test FoMoH and its variants on the Rosenbrock function [] which has a global minimum at $f(\mathbf{1})=0$. The Rosenbrock function, $f(\bm{\theta})=\sum_{i=1}^{D-1}(100(\theta_{i+1}-\theta_{i}^{2})^{2}+(1-\theta
_{i})^{2})$, is designed to be a challenging test case for non-convex optimization where the solution falls inside a narrow valley. The learning rate for the FoMoH variants is set to $1.0$ for the Rosenbrock experiments.

As an initial illustration of the behavior of each forward-mode approach, we show how a single step looks from a randomly chosen starting point in Figure 1(a) for FGD, FoMoH, and FoMoH-BP. We plot the expected (average) step across $10{,}000$ samples for all approaches in the 2D Rosenbrock function. These steps are shown with solid lines. For each approach, we also plot the sampled steps by superimposing a scatter plot in the corresponding approach’s color. All methods are compared to the Newton step, $(\nabla^{2}f(\bm{\theta}))^{-1}\,\nabla f(\bm{\theta})$ shown in red. For FGD, we see that the expected descent direction is the same as the gradient at that point, hence the alignment with FoMoH-BP that directly calculates the gradient. This plot highlights the reliance on a well-chosen learning rate for FGD, whereas FoMoH-BP’s step size is automatically normalized by the local curvature along the gradient. FoMoH (blue), on the overhand, has a descent direction that differs from the gradient and is governed by the distribution of samples that fall on the ellipse defined by the Hessian, $\nabla^{2}f(\bm{\theta})$. For this point, the Newton step is the descent direction that falls on this ellipse and corresponds to the local minimum of the quadratic approximation. Another insight gained from this figure is that the variance of FGD’s descent direction is less constrained than FoMoH’s descent direction (blue), where the sample direction is controlled by the local Hessian. As a final note, we do not plot FoMoH-2D in Figure 1(a) as it directly aligns with the Newton step, with very little variance (see §A, Figure 5).

In Figure 1(b), we now compare the performance of the competing optimization routines for the 2D Rosenbrock function. These results are shown for the same $10$ randomly sampled starting locations for all approaches initialized with different random seeds. We highlight with the thicker line the median performing run. Both axes are on the log scale and show the significant advantage of FoMoH-$K$D, with $K=2$. We also note the advantage (at least for this 2D example) of all FoMoH approaches that use second-order information. Additionally, the two forward-mode only approaches actually outperform the optimization routines that include backpropagation.

Figure: (a) Expected step taken by the stochastic approaches of FoMoH and FGD. Included is a single Newton step for reference. The samples show how the curvature constrains the step size, compared to the sensitivity of FGD to the learning rate.
Refer to caption: x1.png

#### 10D Rosenbrock Function.

We now focus on the performance of FoMoH-$K$D as we increase $K$ from $2$ to the input dimension of the function. Figure 2 shows this comparison for the 10D Rosenbrock function, where we use 10 random initializations for the different $K$. The median performance for each $K$ is then shown, where we see a perfect ordering of performance that aligns with the dimension of the hyperplane. The best performing FoMoH-$K$D is for $K=D$, with the worst corresponding to the lowest dimension implemented, $K=2$. Overall this figure highlights how FoMoH-$K$D trends towards Newton’s method as $K$ tends to $D$, where we actually see the median performances of FoMoH-10D and Newton’s method aligned.

Figure: Figure 2: Performance of FoMoH-$K$D for $K=2\dots 10$ on the 10D Rosenbrock function. Solid lines represent the median, with transparent lines corresponding to the each of the 10 random seeds. There is a clear pattern of higher dimensions performing better, with the performance of $K=10$ coinciding with Newton’s Method (black dotted line).
Refer to caption: x3.png

### 6.2 Logistic Regression

We now compare the performance of FoMoH for logistic regression applied to the MNIST dataset []. Table 1 displays mean and standard deviation performance for each approach, where the forward-mode-only approaches are highlighted separately from the methods that include reverse-mode steps. Additional details on hyperparameter optimization are included in Appendix A.3. Figure 3 displays the training and validation curves for both accuracy and negative log-likelihood (NLL) for the forward-mode approaches (Figure 6 in appendix includes additional reverse-mode approaches). We see an improvement in speed of convergence as $K$ increases. However, we also see that FoMoH and FoMoH-$K$D, with a fixed learning rate, degrade in performance after reaching their local minimum (NLL) or maximum (accuracy). We therefore introduce a learning rate scheduler to improve on this behavior. For this task, FGD is competitive with the FoMoH variants but is slower to converge. In §6.3 we show how FGD degrades in performance for a larger parameter space.

**Table 1: Logistic regression results for MNIST. When comparing the forward-mode-only approaches in the upper section of the table, we see improvement in performance with increasing hyperplane dimension for FoMoH-$K$D. For this logistic regression example, FoMoH-3D and FoMoH-2D with learning rate schedulers, are competitive with FGD. However we will see this result change with a larger dimensional problem in §6.3. Both reverse-mode approaches in the lower section of the table have similar performance, and are included for reference.**
| Approach | Training Loss | Validation Loss | Training Accuracy | Validation Accuracy |
| --- | --- | --- | --- | --- |
| FGD | $0.2976\pm 0.0007$ | $\mathbf{0.2949\pm 0.0017}$ | $0.9154\pm 0.0003$ | $0.9163\pm 0.0019$ |
| FoMoH | $0.3223\pm 0.0013$ | $0.3186\pm 0.0019$ | $0.9073\pm 0.0011$ | $0.9110\pm 0.0021$ |
| FoMoH (LR-Sch.) | $0.3192\pm 0.0012$ | $0.3160\pm 0.0025$ | $0.9085\pm 0.0011$ | $0.9118\pm 0.0021$ |
| FoMoH-2D | $0.3010\pm 0.0018$ | $0.3015\pm 0.0031$ | $0.9144\pm 0.0009$ | $0.9149\pm 0.0015$ |
| FoMoH-2D (LR-Sch.) | $0.2921\pm 0.0014$ | $0.2951\pm 0.0027$ | $0.9174\pm 0.0005$ | $\mathbf{0.9170\pm 0.0010}$ |
| FoMoH-3D | $0.3449\pm 0.0017$ | $0.3343\pm 0.0034$ | $0.8999\pm 0.0008$ | $0.9036\pm 0.0023$ |
| FoMoH-3D (LR-Sch.) | $\mathbf{0.2893\pm 0.0017}$ | $0.3054\pm 0.0034$ | $\mathbf{0.9195\pm 0.0007}$ | $0.9153\pm 0.0010$ |
| FoMoH-BP | $0.2312\pm 0.0001$ | $0.2679\pm 0.0003$ | $0.9366\pm 0.0001$ | $0.9267\pm 0.0002$ |
| Backpropagation | $0.2265\pm 0.0000$ | $0.2710\pm 0.0001$ | $0.9381\pm 0.0001$ | $0.9267\pm 0.0003$ |

Figure: Figure 3: Forward-mode training and validation curves for the logistic regression model on the MNIST dataset. Average and standard deviation is shown for five random initializations.
Refer to caption: x4.png

### 6.3 Convolutional Neural Network

We now move from a model with $7{,}850$ parameters to a convolutional neural network (CNN) with $431{,}080$ parameters. As before we use the MNIST dataset and leave the details on hyperparameter optimization using Bayesian optimization to the Appendix A.4.
Both Table 2 and Figure 4 highlight the advantage of FoMoH-$K$D over FGD. For the larger parameter space FGD requires more epochs to converge compared to all FoMoH variants. The learning rate scheduler further improves FoMoH and FoMoH-$K$D by helping to avoid getting stuck in low performance regions. Here, we see the clear trend that the best performing forward-mode-only approach comes from the largest $K$, which was $K=3$ for this experiment. As expected, both optimizers FoMoH-BP and Backpropagation outperform the forward-mode-only approaches. Overall, these results highlight that second-order information helps scale the performance of forward-mode optimization to larger dimensions.

**Table 2: CNN results for MNIST. The forward-mode-only approaches in the upper section of the table show that FoMoH’s performance improves with the dimension of $K$, especially when used with the learning rate scheduler. FoMoH-3D outperforms FoMoH-2D, FoMoH, and FGD. The reverse-mode approaches in the lower section outperform forward-mode, with BP slightly better than FoMoH-BP.**
| Approach | Training Loss | Validation Loss | Training Accuracy | Validation Accuracy |
| --- | --- | --- | --- | --- |
| FGD | $0.1211\pm 0.0097$ | $0.1104\pm 0.0086$ | $0.9641\pm 0.0038$ | $0.9677\pm 0.0030$ |
| FoMoH | $0.1663\pm 0.0069$ | $0.1571\pm 0.0092$ | $0.9518\pm 0.0011$ | $0.9550\pm 0.0002$ |
| FoMoH (LR-Sch.) | $0.1617\pm 0.0119$ | $0.1575\pm 0.0158$ | $0.9515\pm 0.0035$ | $0.9539\pm 0.0034$ |
| FoMoH-2D | $0.1015\pm 0.0041$ | $0.1016\pm 0.0066$ | $0.9691\pm 0.0011$ | $0.9693\pm 0.0020$ |
| FoMoH-2D (LR-Sch.) | $0.0900\pm 0.0070$ | $0.0913\pm 0.0041$ | $0.9731\pm 0.0022$ | $0.9718\pm 0.0014$ |
| FoMoH-3D | $0.1085\pm 0.0105$ | $0.1073\pm 0.0133$ | $0.9674\pm 0.0022$ | $0.9686\pm 0.0028$ |
| FoMoH-3D (LR-Sch.) | $\mathbf{0.0809\pm 0.0061}$ | $\mathbf{0.0923\pm 0.0132}$ | $\mathbf{0.9759\pm 0.0014}$ | $\mathbf{0.9734\pm 0.0022}$ |
| FoMoH-BP | $0.0093\pm 0.0016$ | $0.0310\pm 0.0006$ | $0.9981\pm 0.0005$ | $0.9903\pm 0.0003$ |
| Backpropagation | $0.0053\pm 0.0034$ | $0.0329\pm 0.0032$ | $0.9990\pm 0.0009$ | $0.9909\pm 0.0004$ |

Figure: Figure 4: Forward-mode training and validation curves for the CNN on the MNIST dataset. Average and standard deviation is shown for three random initializations. Note how FGD (blue) is much slower to converge, with FoMoH-$K$D improving in performance with increasing $K$.
Refer to caption: x5.png

## 7 Conclusion, Broader Impact, and Limitation Discussion

The results in §6 highlight the potential of the use of second-order forward-mode AD for optimization tasks. For the Rosenbrock function, we illustrated the behavior of our three new optimization routines: FoMoH, FoMoH-BP, and FoMoH-$K$D, and we compared them to Newton’s method. In particular we were able to show that as one increases the hyperplane dimension of FoMoH-$K$D the method tends to Newton’s method, without the need for any backpropagation. This significant result is shown in Figure 2. For the learning tasks of logistic regression and CNN classification, we see how the first-order optimization approach of FGD degrades with increasing dimension of the parameter space. We do not see this degradation for FoMoH-$K$D, and we also observe that the second-order information means fewer epochs are needed to reach a better performance. This has the broader impact of improving efficiency, reducing cost, and increasing accuracy in ML optimization routines.

In conclusion, we introduced a novel approach that uses second-order forward-mode AD for optimization. We have introduced: forward-mode line search (FoMoH); forward-mode line search with Backpropagation (FoMoH-BP); and forward-mode hyperplane search (FoMoH-$K$D). We have shown how these approaches compare to the previous first-order forward-mode approach of FGD, as well stochastic gradient descent for multiple optimization problems over a wide range of dimensions. Furthermore, FoMoH-$K$D behaves closer to the performance of Newton’s method as $K$ increases. In addition to contributing the new second-order forward-mode optimization routines, we provide a Python package that implements the AD backend and interfaces with PyTorch.
Our work is a further step in the direction
of showing the value of cheap second-order information in optimization with the need to scale to even larger dimensions.
We expect that future work will be able to mix the advantages gained from first-order approaches with that of second-order approaches.

## Acknowledgments and Disclosure of Funding

This material is based upon work supported by the United
States Air Force and DARPA under Contract No. FA8750-23-C-0519. Any opinions, findings
and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect
the views of the United States Air Force and DARPA.

## Appendix A Additional Results

### A.1 FoMoH- K 𝐾 K italic_K D Algorithm

Figure: Algorithm 2 FoMoH-$K$D

### A.2 Rosenbrock Function

Figure: Figure 5: Histogram over expected step taken by the stochastic approaches of FoMoH, FGD, and FoMoH-2D corresponding to Figure 1(a). Noteworthy is that the variance of the 2D hyperplane search step is significantly smaller and expectation is close to Newton step.
Refer to caption: x6.png

### A.3 Logistic Regression

Table 3 includes the final hyperparameter selection for the experimental results in §6.2. We used [] to perform a grid search with 100 iterations, where the batch size choice was between [128, 512, 1024, 2048]. For the learning rate scheduler, we reduced the learning rate at the epoch where the NLL starts to increase. For FoMoH-3D we multiplied the learning rate by $0.8$, whereas for the other approaches we multiplied the learning rate by $0.1$. There is likely room for improvement on the parameters of the learning rate scheduler, but that would only improve the current results.

Figure 6 includes the reverse-mode training and validation curves for Backpropagation and FoMoH-BP in addition to the curves shown in Figure 3.

**Table 3: Hyperparameter Optimization for Logistic Regression.**
| Approach | Learning Rate | Learning Rate Bounds | Batch Size |
| --- | --- | --- | --- |
| FGD | 0.00006497 | [0.00001, 0.1] | 128 |
| FoMoH | 0.1362 | [0.001, 1.0] | 1024 |
| FoMoH (LR-Sch.) | 0.1362 | [0.001, 1.0] | 1024 |
| FoMoH-2D | 0.04221 | [0.001, 1.0] | 512 |
| FoMoH-2D (LR-Sch.) | 0.04221 | [0.001, 1.0] | 512 |
| FoMoH-3D | 0.1 | [0.001, 1.0] | 512 |
| FoMoH-3D (LR-Sch.) | 0.1 | [0.001, 1.0] | 512 |
| FoMoH-BP | 0.04688 | [0.01, 1.0] | 2048 |
| Backpropagation | 0.03561 | [0.01, 0.5] | 2048 |

Figure: Figure 6: Training and validation curves for logistic regression model on the MNIST dataset. Average and standard deviation is shown for five random initializations.
Refer to caption: x7.png

### A.4 CNN

Table 4 includes the final hyperparameter selection for the experimental results in §6.3. We used [] to perform Bayesian optimization with 100 iterations, where the batch size choice was fixed to 2048. For FoMoH-3D, we used the same hyperparameters as for FoMoH-2D as this gave sufficient performance (and still outperformed the other forward-mode approaches). All learning rate schedulers reduced the learning rate by 10 every 1000 epochs.

Figure 7 includes the reverse-mode training and validation curves for Backpropagation and FoMoH-BP in addition to the curves shown in Figure 4.

**Table 4: Hyperparameter Optimization for Logistic Regression.**
| Approach | Learning Rate | Learning Rate Bounds | Batch Size |
| --- | --- | --- | --- |
| FGD | 0.0001376 | [0.00001, 0.1] | 2048 |
| FoMoH | 0.542 | [0.001, 1.0] | 2048 |
| FoMoH (LR-Sch.) | 0.542 | [0.001, 1.0] | 2048 |
| FoMoH-2D | 0.3032 | [0.001, 1.0] | 2048 |
| FoMoH-2D (LR-Sch.) | 0.3032 | [0.001, 1.0] | 2048 |
| FoMoH-3D | 0.3032 | [0.001, 1.0] | 512 |
| FoMoH-3D (LR-Sch.) | 0.3032 | [0.001, 1.0] | 2048 |
| FoMoH-BP | 0.04688 | [0.01, 1.0] | 2048 |
| Backpropagation | 0.03561 | [0.005, 0.2] | 2048 |

Figure: Figure 7: Training and validation curves for CNN on the MNIST dataset. Average and standard deviation is shown for three random initializations.
Refer to caption: x8.png

### A.5 Computational Resources

All experiments are run on a NVIDIA RTX 6000 GPU. The main compute cost came from both the grid search and the Bayesian optimization that we ran over the six different optimization routines for the different experiments.

## Peer Reviews

confidence: 5
contribution: 3
initial_rating: 3
presentation: 3
questions: - First it would be great to revise the proofs to make them rigorous.
- Could you add a full mathematical definition of FoMoH-BP?
- What is the logistic regression model? I suppose it is not a regression but a classification problem first? Then if you use 7850 parameters it's probably not a simple linear model but some form of Multi-Layer Perceptron?
- Detail the CNN architecture used in the experiment.
- In the abstract, you mention alternative (orthogonal to be exact) methods for forward gradients. They are never compared in the experiments. It would be great to have them.
- By curiosity how can analog optical systems compute derivatives of intermediate functions to implement forward mode automatic differentiation? The reference provided by the authors does not mention AD, at least from its abstract.
- You mention "linesearch" in line 243 but there is no linesearch at all in the algorithm. What do you mean by linesearch?
soundness: 2
strengths: - The idea of exploiting second order directional derivatives is original and could be further explored.
- The proposed method shows clearly superior performance than a simple forward gradient descent.
summary: The authors present a new optimization method relying on forward mode automatic differentiation (AD). Namely, the authors propose to use second order directional derivatives computed along random directions to precondition stochastic estimates of the gradients obtained by forward mode automatic differentiation along these same random directions. The proposed Forward Mode Second-order Hyperplane Search (FoMoH) interpolates between using an approximate Cauchy stepsize (that is a Cauchy stepsize computed with a quadratic approximation of the objective) along a random direction and an approximate Newton step. The proposed method is shown to outperform a standard forward gradient descent on the Rosenbrock function, a logistic regression problem, and the MNIST image classification task with a CNN.
weaknesses: Unfortunately the paper does not give justice to the potential of the main idea.
Improvements are necessary and possible:
- The method does not require introducing dual numbers. Computing second order directional derivatives can easily be done with nested forward mode autodiff:
```
import jax

def hqp(fun, w, v1, v2):
  def dir_der_v1(w):
    return jax.jvp(fun, (w,), (v1,))[1]
  return jax.jvp(dir_der_v1, (w,), (v2,))[1]
```
The above may be slightly slower than the implementation with dual numbers (the difference is probably minimal) but it is also much simpler to implement than having to adopt a new kind of automatic differentiation library. If one want the best possible implementation, then Taylor-mode automatic differentiation can also be used but again the benefits are really minor. Algorithm 1 is then unnecessarily complicated when it could be written in just K calls to the above function to define the approximate Hessian. Presenting the algorithm with a simple implementation that any user could recode in at most 50 lines of code in jax or pytorch would greatly improve the potential adoption of the method.
- Unfortunately, the proofs are not rigorous, nor are the claims.
  - Theorem 1 is a corollary of Theorem 2 so no need to present it. Also results like $\lim_{t \rightarrow +\infty} \theta_t = \theta^*$ are meaningless: we don't want to wait the time of the universe to see convergence. Rates like the ones provided in Theorem 2 are relevant.
  - In all claims, detail the setting: what algorithm is used, what is theta_t, what is the expectation taken against etc... You may not do that in the main text by lack of space but at least make sure that the appendix contains a result that details all assumptions.
  - The proof of theorem 2 is unfortunately not well detailed:
    - Please give a detailed proof that $\tilde \theta_{t+1} = \tilde \theta_t + P(\tilde \theta^* - \tilde \theta_t)$. We are dealing with quadratics so the proof should boil down to simple linear algebra. Avoid intuitive arguments, just write the equations one by one showing the results. I personally will refuse this paper to be accepted without detailed proofs.
    - Give a proper reference for the fact that the expectation of a projection matrix defined from Gaussian variables is the scaled identity
    - First line of last set of equations of the proof of theorem 2 should read $\tilde \theta_{t+1} - \tilde \theta^* = (1- K/D)(\tilde \theta_t - \tilde \theta^*)$
    - Please be rigorous when you write expectations. You need to detail each time with respect to which randomness you are taking the expectation. For example at one point you write $\tilde \theta_{t+1} = \mathbb{E}[\tilde \theta_t + P(\tilde \theta^* - \tilde \theta_t)]$ but so then $\tilde \theta_{t+1}$ is not random. Unless you meant $E[\tilde \theta_{t+1}] = $. Such lack of rigor is detrimental to the potential of the idea.
    - Section B.4 contains multiple errors:
      - Reaching a critical point does not imply that you reach a minimum unless you make additional assumptions like convexity.
      - Again be rigorous in the use of expectations, one usually use conditional expectations conditioned on the previous iterate for example.
      - The provided rate is clearly not linear. Consider rereading in details the reference for example.
- Second order methods may generally be very sensitive to the batch-size. It would be great to plot a sensitivity analysis of the method with the batch-size for e.g. a given learning rate.
- Consider another dataset than MNIST. MNIST is well known to be particularly easy and may not reflect potential challenges that the method can have.

---

confidence: 2
contribution: 3
initial_rating: 5
presentation: 3
questions: 1. Is it correct to say that the method performs Newton steps in the sampled subspaces (assuming learning rate being 1)? If not, what's the relationship between the two?
2. the $\kappa_i$'s need not be positive, correct?
3. in the learning tasks, are also examples being smapled (i.e. mini-batch updates?). If so, is the proposed method sensitive to the mini-batch size?
4. Do you think the proposed method could be relevant also for forward-mode gradient-based hyperparameter optimization?
soundness: 2
strengths: - The significance of developing a preforming optimization method that does not rely on reverse-mode differentiation is very high, and the proposed method seems to make a concrete step in this direction 
- The paper is mostly well written and easy to follow, the notation is clear (although a few passages could be made clearer, see below)
- The work does a good job introducing the concept of dual and hyper-dual numbers, which I expect the community to benefit from 
-  The theoretical results clearly show the advantages of the proposed method over FGD, as well as the larger scale experiment.
summary: The paper presents a forward-mode-only optimization method that uses second order information on randomly sampled hyperplanes. 
The method makes use of (K^2 + K)/2 calls to a "double" forward-mode AD, implemented with hyper-dual numbers, computing along the way the Hessian projected onto the drawn plane. The paper presents theoretical results on convex and quadratic functions that also relates it to the Newton method, alongside some empirical validation on a test function and two learning problems.
weaknesses: - One main promise of forward mode differentiation is to vastly decrease the memory requirement for large models, however the paper does not empirically quantifies the advantage of FoMoH in this regard. It would be nice to include memory footprint comparisons (as well as perhaps a table summarizing the runtime and memory complexity of key algorithms)
- I would have appreciated some larger scale experiment with transformer architectures, e.g. for fine-tuning LLMs, which could be an interesting application of the proposed method 
- Some details could be better specified in the main paper:
    -  It is not entirely clear to me by reading the paper if the (hyper-)dual numbers allow for the computation of JVPs and bilinear hessian products by just "tracking the epsilons" while computing the function, or if one has actually to implement the above-mentioned operations (I think it's the first). One "implementation" example would help focus ideas.
    - I'm a bit confused about the origin of Equation (1) right. Do the expressions for the $\kappa$ come from manipulation of the resulting Taylor expansion or is it "set by design" (to mimic Newton updates)? In general, I would have appreciated some more details around lines 227 to 240 
   - the learning rate scheduler seems like an important addition for empirical performance, some more details (e.g. which scheduler, how did you choose its hyperparamters) in the main paper would be welcomed.
   - some more discussion on relations between FoMoH-1d, FoMoH-BP and FGD would have been nice. For me it is not immediately clear why FoMoH-1d consistently underperforms FGD on the learning tasks and why FoMoH-BP consistently performs on par with BP. 
  - what does BP stand for in the experiments? Is it plain (stochastic?) gradient descent or some other adaptive method?

---

confidence: 2
contribution: 2
initial_rating: 5
presentation: 3
questions: 1. Why is this paper suitable for publication as machine learning research? What is the difficult point of the optimization method in machine learning, and how does the paper address it? If the approximation of a second-order method using dual numbers is new, why has no literature in optimization theory or operations research discussed it?

1. Does not the overhead of the proposed method increase the runtime in one iteration? If it does, the proposed method is still faster than baselines in terms of runtime until convergence?

1. How is the scalability of the proposed method?
soundness: 3
strengths: 1. This paper is well-written and easy to follow. Assumptions and mathematical expansions are explained in detail.
1. Experimental results demonstrate the proposed method is more efficient than baseline first-order methods in terms of iterations of parameter updates.
1. Convergence properties are proved in theory, and the proposed method is guaranteed to be consist with Newton's method.
summary: This paper investigates the effectiveness of the second-order differentiation in forward automated differentiation.
This paper and regards the forward propagation with dual number as the objective function applied second-order Taylor-series expansion to.
By using dual numbers, the proposed method estimates the Hessian matrix and applies the approximated Newton's method for the optimization.
weaknesses: 1. I am not very sure that this paper addresses a new topic in machine learning.
The explanation of the proposed method seems to address general optimization problems but does not seem to focus on the problem in machine learning.
For examples, this paper does not explicitly address the difficulty of accessing the objective function due to large datasets like SGD, and not address the difficulty of non-convexity caused by nonlinear models.
Since I do not have much expertise in optimization theory or operations research, I cannot evaluate the novelty of this paper well in the context of optimization theory. 
Even so, I doubt the proposed method is very novel because the used mathematical tools are fundamental and the addressed objective function does not seem very difficult. 
Is the research problem specialized for machine learning problems? And, is the paper new even in the context of the optimization problem?


1. While this paper evaluates the convergence of the proposed method in terms of iterations, it does not evaluate the runtime of the proposed method.
The proposed method requires an inverse of Hessian matrices, and I think its computational cost can be high.
Does the overhead of the proposed method not increase the runtime in one iteration? If it does, the proposed method is still faster than baselines in terms of runtime until convergence?
To emphasize the practical usefulness of the proposed method, this paper needs the evaluation of runtime.

1. It is not clear that the proposed method is scalable for recent deep neural network model architectures.
(Fournier et al., 2023) seems to show that the first-order forward gradient method is applicable to ResNet18. Is the proposed method applicable to such modern architectures?

Meta-review: The paper introduces a new second-order forward-mode automatic differentiation method for optimization. It generalizes a second-order line search to a K-dimensional hyperplane search using hyper-dual numbers. The method is evaluated on the Rosenbrock function, logistic regression, and learning a CNN classifier. 

According to the reviewers, the paper is well-written and easy to follow. The proposed method is shown to be more efficient than baseline first-order forward-mode methods in terms of iterations. 

However, reviewers pointed out the lack of coverage of related works, that the method is not scalable for recent deep neural network model architectures, and that the paper does not evaluate the runtime of the proposed method. It was also questioned whether dual numbers are really necessary to present the proposed method. More importantly, a reviewer pointed out holes in the proof of Theorem 2 that the authors were unable to fix.

As a result, we believe the paper is not ready for publication yet.

Decision: Reject