Nonlinear Splitting for Gradient-Based Unconstrained and Adjoint Optimization

TL;DR

This paper introduces nonlinear splitting techniques for gradient-based optimization, enhancing accuracy and efficiency in high-dimensional and nonconvex problems across various fields, including PDE-constrained optimization.

Contribution

The paper proposes a novel nonlinear splitting framework for gradient and constraint evaluations, improving convergence and reducing computational costs in complex optimization tasks.

Findings

Outperforms existing methods in accuracy and runtime.

Reduces the number of high-dimensional solves by a factor of three.

Demonstrates effectiveness on diverse optimization problems.

Abstract

High dimensional and/or nonconvex optimization remains a challenging and important problem across a wide range of fields, such as machine learning, data assimilation, and partial differential equation (PDE) constrained optimization. Here we consider gradient-based methods for solving unconstrained and constrained optimization problems, and introduce the concept of nonlinear splitting to improve accuracy and efficiency. For unconstrained optimization, we consider splittings of the gradient to depend on two arguments, leading to semi-implicit gradient optimization algorithms. In the context of adjoint-based constrained optimization, we propose a splitting of the constraint $F(\mathbf{x},\theta)$, effectively expanding the space on which we can evaluate the ``gradient''. In both cases, the formalism further allows natural coupling of nonlinearly split optimization methods with acceleration techniques, such as Nesterov or Anderson acceleration. The framework is demonstrated to outperform existing methods in terms of accuracy and/or runtime on a handful of diverse optimization problems. This includes low-dimensional analytic nonconvex functions, high-dimensional nonlinear least squares in quantum tomography, and PDE-constrained optimization of kinetic equations, where the total number of high-dimensional kinetic solves is reduced by a factor of three compared with standard adjoint optimization.

Tables

Table 1. Table 1. List of nonlinear splitting methods for unconstrained and equality-constrained optimization

Name	Optimization type	Method	Standard counterpart
NS-GD	Unconstrained	(3.2)	GD (2.2)
NS-NAG	"	(3.10)	NAG (3.9)
NS-GD-Newton$\left(1\right)$	"	(3.12)	GD
NS-NAG-Newton$\left(1\right)$	"	(3.13)	NAG
NS-AA-Newton$\left(1\right)$	"	(3.14)	AA applied to GD
NS-Adj-GD	Equality-constrained	Algorithm˜2	Adj-GD (see Section˜4.2.1)
NS-Adj-NAG	"	Algorithm˜3	Adj-NAG (see Section˜4.2.1)
NS-Adj-AA	"	Algorithm˜4	Adj-AA (see Section˜4.2.1)

Table 2. Table 2. Comparison of runtime and cost for the three benchmark functions, averaged over 100 randomized initial states, for NS-GD, NS-NAG, and MATLAB optimizers.

	Rastrigin		Rosenbrock		Beale
Method	Runtime (s)	Cost	Runtime (s)	Cost	Runtime (s)	Cost
NS-GD	0.0064	4.3e-12	1.42	2.09e-9	0.12	2.27e-9
NS-NAG	0.0089	3.49e-11	1.45	2.46e-5	0.11	1.07e-10
FM	0.0012	4.97	0.0013	9.85	0.0013	1.02e-9
PS	0.018	3.41e-11	0.026	8.01e-10	0.017	1.03
GA	0.010	1.06	0.012	1.24	0.0085	0.18
PartS	0.0085	0.22	0.011	0.18	0.0055	0.0093
Surr	1.36	0	2.61	0.014	1.59	0.013
GS	0.059	0	0.31	0	0.11	3.93e-14

Equations

$$\underset{\mathbf{x}}{\min }C(\mathbf{x}),$$

$$\underset{\theta }{\min }J(\mathbf{x};\theta ),\text{subject to }F(\mathbf{x};\theta )=0.$$

$$\dot{\mathbf{x}}=-\nabla C(\mathbf{x}),$$

$${\mathbf{x}}^{k+1}$$

$${\mathbf{x}}^{k+1}={\mathbf{x}}^k-\frac{1}{L}{A}^{-1}\nabla C({\mathbf{x}}^k),$$

$$\ddot{\mathbf{x}}+{\alpha }_1\dot{\mathbf{x}}+{\alpha }_2\nabla C(\mathbf{x})=0,$$

$${\text{prox}}_{\gamma C}(\mathbf{x})=\underset{\mathbf{v}}{\mathrm{argmin}}\left(C(\mathbf{v})+\frac{1}{2\gamma }\Vert \mathbf{x}-\mathbf{v}{\Vert }^2\right).$$

$${\mathbf{x}}^{k+1}={\text{prox}}_{{\gamma }_{k}C}({\mathbf{x}}^k).$$

$${\mathbf{x}}^{k+1}+{\gamma }_k\nabla C({\mathbf{x}}^{k+1})={\mathbf{x}}^k.$$

$${\mathbf{x}}^{k+1}={\text{prox}}_{{\gamma }_k{C}_2}({\mathbf{x}}^k-{\gamma }_k\nabla C_1({\mathbf{x}}^k))=\underset{\mathbf{v}}{\mathrm{argmin}}\left(C_2(\mathbf{v})+\frac{1}{2{\gamma }_k}\Vert {\mathbf{x}}^k-{\gamma }_k\nabla C_1({\mathbf{x}}^k)-\mathbf{v}{\Vert }^2\right).$$

$${\mathbf{x}}^{k+1}+{\gamma }_k\nabla C_2({\mathbf{x}}^{k+1})={\mathbf{x}}^k-{\gamma }_k\nabla C_1({\mathbf{x}}^k).$$

$${\mathbf{y}}^{k+1}={\mathbf{y}}^k+\Delta t\tilde{F}({\mathbf{y}}^k,{\mathbf{y}}^{k+1}).$$

$${\mathbf{y}}^{k+1}={\mathbf{y}}^k+\Delta t{J}_2[{\mathbf{y}}_k{]}^{-1}\tilde{F}({\mathbf{y}}^k,{\mathbf{y}}^k),$$

$$\tilde{\nabla }C(\mathbf{x},\mathbf{x})=\nabla C(\mathbf{x})$$

$${\mathbf{x}}^{k+1}={\mathbf{x}}^k-{\gamma }_k\tilde{\nabla }C({\mathbf{x}}^k,{\mathbf{x}}^{k+1}).$$

$$\nabla C(\mathbf{x})=L(\mathbf{x})\mathbf{x}+g(\mathbf{x}),$$

$$\tilde{\nabla }C(\mathbf{x},\mathbf{y}):=L(\mathbf{x})\mathbf{y}+g(\mathbf{x}).$$

$${\mathbf{x}}^{k+1}={\mathbf{x}}^k-{\gamma }_k(L({\mathbf{x}}^k){\mathbf{x}}^{k+1}+g({\mathbf{x}}^k)),$$

$${\mathbf{x}}^{k+1}=(I+{\gamma }_{k}L({\mathbf{x}}^k){)}^{-1}({\mathbf{x}}^k-{\gamma }_{k}g({\mathbf{x}}^k)).$$

$$\tilde{\nabla }C(\mathbf{x},\mathbf{y})={\nabla }_{\mathbf{y}}{D}_{\mathbf{x}}(\mathbf{y})\text{for}{D}_{\mathbf{x}}(\mathbf{y})={\mathbf{y}}^T\mathbf{A}(\mathbf{x})\mathbf{y}+\mathbf{B}(\mathbf{x})\mathbf{y}+\mathbf{b}.$$

$$\tilde{\nabla }C(\mathbf{x},\mathbf{y})$$

$$C(\mathbf{x},\mathbf{y})=\frac{1}{2}a(\mathbf{x}{)}^T(\mathbf{y}\odot \mathbf{y})$$

$$⟨\tilde{\nabla }C(\mathbf{x},\mathbf{y}),\mathbf{y}-\mathbf{x}⟩=C(\mathbf{y})-C(\mathbf{x}),$$

$$C({\mathbf{x}}^{k+1})\le C({\mathbf{x}}^k),$$

$$\tilde{\nabla }C(\mathbf{x},\mathbf{y}):={\int }_0^1\nabla C(\xi \mathbf{x}+(1-\xi )\mathbf{y})d\xi .$$

$${\mathbf{x}}^{k+1}$$

$${\mathbf{v}}^{k+1}$$

$${\mathbf{x}}^{k+1}$$

$${\mathbf{v}}^{k+1}$$

$$G(\mathbf{x}):=\mathbf{x}-\gamma \nabla C(\mathbf{x}).$$

$$G_A(\mathbf{x}):=\mathbf{x}-\gamma A(\mathbf{x}{)}^{-1}\nabla C(\mathbf{x}),$$

Full text

1. 1 Introduction

2. 2 Gradient flow and nonlinear splitting

3. 3 Gradient-based optimization with nonlinear splitting

4. 3.1 Nonlinear splitting for unconstrained optimization

5. 3.1.1 Acceleration

6. 3.2 Nonlinear splitting with equality constraints

7. 3.2.1 Acceleration

8. 4 Numerical examples

9. 4.1 Unconstrained examples

10. 4.1.1 Benchmark nonconvex optimization problems

11. 4.1.2 Nonlinear least squares: pure state quantum tomography

12. 4.2 Examples with equality constraints

13. 4.2.1 Nonlinear neutron transport equation

14. 4.2.2 Splitting of the initial condition for an initial value problem

15. 5 Conclusion

16. A The adjoint of the neutron transport equation

1. Introduction

High dimensional and/or nonconvex optimization remains a challenging and important problem across a wide range of fields, such as machine learning, data assimilation, and partial differential equation (PDE) constrained optimization. This paper is concerned with gradient-based methods for solving two classical optimization problems, namely the unconstrained optimization problem

[TABLE]

and the constrained optimization problem

[TABLE]

For both problems, we will employ nonlinear splittings to improve the efficiency and accuracy of gradient-based optimization. Specifically, given a function $f(\mathbf{x})$, we define a nonlinear splitting of $f$ to be a function $g(\mathbf{x},\mathbf{y})$ that satisfies $g(\mathbf{x},\mathbf{x})=f(\mathbf{x})$ and is differentiable in both arguments. In the context of unconstrained optimization, we will consider splittings of the gradient $\nabla C$ to depend on two arguments, leading to semi-implicit gradient optimization algorithms. For the constrained optimization (1.2), the constrained objective is often sufficiently complex that treating components of the gradient implicitly may not be tractable. Instead, we propose a splitting of the constraint $F(\mathbf{x};\theta)$ within an adjoint-based optimization. This effectively expands the space on which we can evaluate the “gradient,” and induces a sort of splitting in the gradient as well. In both cases, the formalism further allows natural coupling of nonlinearly split optimization methods with acceleration techniques, such as Nesterov or Anderson acceleration.

This paper is organized as follows. We begin by motivating nonlinear splitting in the context of gradient flow and classical algorithms like gradient descent, proximal point method, and proximal gradient method in Section˜2. In Section˜3, we develop the nonlinear splitting framework for gradient-based optimization. Section˜3.1 focuses on the unconstrained setting, where we introduce the notion of a nonlinear splitting of a gradient and use this to define the nonlinear splitting gradient descent method. In Section˜3.1.1, we discuss momentum-based and fixed-point-based acceleration of nonlinear splitting gradient descent. Throughout, we discuss connections to existing methods and generalizations thereof. In Section˜3.2, we discuss nonlinear splitting in the setting of equality-constrained optimization. Here, we introduce the notion of a nonlinear splitting of the equality constraint and a nonlinear splitting of the gradient. Using a modified adjoint approach, we show how to systematically construct a nonlinear splitting of the gradient, given a nonlinear splitting of the equality constraint. A general framework is provided where any nonlinear approximation or solver written in fixed-point form can naturally induce such a splitting. In Section˜3.2.1, we provide momentum-based and fixed-point-based accelerated variants of the equality-constrained nonlinear splitting gradient method.

In Section˜4, we demonstrate these proposed methods on a variety of numerical examples. For unconstrained optimization, we consider two examples in Section˜4.1 that highlight the advantage of the linear stability gained by using nonlinear splitting gradient descent. Specifically, in Section˜4.1.1, we consider several two-dimensional non-convex benchmark unconstrained optimization problems, and demonstrate how the linear stability from these nonlinear splitting methods allows one to take large steps in order to escape local minima and find the global minima of these non-convex functions. In terms of objective value versus runtime, these nonlinear splitting methods outperform most of MATLAB’s built-in optimizers, with MATLAB’s global search algorithm being the only one that is comparable. That being said, the nonlinear splitting methods are implemented via naive looping in MATLAB scripts and thus we expect some performance degradation compared with the optimized methods provided in MATLAB. In Section˜4.1.2, we consider a nonlinear least squares problem to reconstruct a pure quantum state arising in quantum tomography. This example demonstrates linearization of the nonlinear splitting methods, using a single Newton iteration, as well as Anderson acceleration of the corresponding method. For equality-constrained optimization, we consider two examples in Section˜4.2. In Section˜4.2.1, we consider a source optimization problem for the nonlinear neutron transport equation, where we demonstrate that our proposed methods significantly reduce the number of linear transport solves in order to achieve the same objective value and equality constraint residual, compared to their standard counterparts. In Section˜4.2.2, we show how the nonlinear splitting framework can be applied in the case when the equality constraint is given by the dynamics of an ordinary differential equation (ODE), which leads us to a momentum-based scheme for problems with time dependence.

2. Gradient flow and nonlinear splitting

The idea of using ordinary differential equation (ODE) theory and techniques to solve nonlinear root-finding and optimization problems for cost function $C(\mathbf{x})$ dates back many decades, and many optimization methods can be seen as certain discrete approximations to continuous gradient flow equations,

[TABLE]

or continuous modifications thereof. It is also common for optimization algorithms to be based on solving local optimization problems during each step/iteration towards the global solution. For example, suppose $C(\mathbf{x})$ has a Lipschitz continuous differential $DC(\mathbf{x})$ with constant $L$, then [26, Lemma 54] $C(\mathbf{x}+\mathbf{h})\leq C(\mathbf{x})+DC(\mathbf{x})\mathbf{h}+\frac{L}{2}\|\mathbf{h}\|^{2}$. Taking a greedy minimization with respect to $\mathbf{h}$ yields the gradient descent algorithm,

[TABLE]

which is exactly a forward Euler discretization of the continuous gradient flow (2.1). If we have a matrix-induced $A$-norm for symmetric positive definite $A$ (or natural gradients with a locally defined inner product [4]), $\|\mathbf{h}\|_{A}^{2}=\mathbf{h}^{T}A\mathbf{h}$, and minimization with respect to $\mathbf{h}$ yields the correction

[TABLE]

which is a linearly preconditioned gradient descent [34; 32; 18]. This remains an explicit approximation to the continuous gradient flow, except here we take the gradient with respect to the $A$ inner product, where $\nabla_{A}C(\mathbf{x})=A^{-1}\nabla C(\mathbf{x})$. Newton’s method can be derived from continuous gradient flow with respect to a Hessian-induced inner product (or duality pairing for non-SPD Hessians [50]), and Levenberg–Marquardt follows from explicit discrete approximations to the modified continuous gradient flow $\dot{\mathbf{x}}=-H(\mathbf{x})^{-1}\nabla B(\mathbf{x})$, e.g. [1; 12; 19].

Damping or momentum can also be introduced via second-order formulations. The heavy-ball with friction method, which has had a resurgence of interest due to the development of large-scale optimization and machine learning problems, was developed in 1964 [41]. There, they consider multistep iterative methods to solve linear and nonlinear systems, and relate the multistep methods to the continuous second-order ODE

[TABLE]

for objective function $C(\mathbf{x})$. More general second order equations for optimization are typically of the form $\mu(t)\ddot{\mathbf{x}}+\beta(t)\dot{\mathbf{x}}+\nabla C(\mathbf{x})/\mu(t)=0$, e.g., [25; 2; 13; 3; 7]. The continuous limit of Nesterov’s accelerated gradient (NAG) descent, for example, is given by $\ddot{\mathbf{x}}+\frac{3}{t}\dot{\mathbf{x}}+\nabla C(\mathbf{x})=0,$ with $\mathbf{x}(0)=\mathbf{x}_{0},\dot{\mathbf{x}}(0)=0$ and a change of variable in time to decay the learning rate [48].

Above are all explicit approximations to continuous gradient flows. Although each step can be posed as a certain local optimization problem, the solution is either analytic or linear in nature. In contrast, proximal methods are based around optimizing or solving some “proximal,” i.e. nearby problem. The typical $\ell^{2}$ proximal operator $\textnormal{prox}_{F}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ with constant $\gamma_{k}>0$ takes the form

[TABLE]

Much of the literature on proximal methods is specifically for convex functions, where $\mathbf{x}^{*}$ minimizes $C(\mathbf{x})$ if and only if $\mathbf{x}^{*}=\textnormal{prox}_{C}(\mathbf{x}^{*})$, suggesting the fixed-point algorithm [43]

[TABLE]

This is known as the proximal point method. If we assume that $C(\mathbf{x})$ is a closed convex function with unique minimum such that $\nabla C(\mathbf{x}^{*})=\mathbf{0}$ (or we can simply relax the objective to finding a local rather than global minimum in (2.5)), note that $\textnormal{prox}_{\gamma_{k}C}(\mathbf{x})=\left\{\mathbf{v}\textnormal{ : }\mathbf{v}+\gamma_{k}\nabla C(\mathbf{v})=\mathbf{x}\right\}$. Then the proximal point algorithm (2.6) corresponds to an implicit backward Euler discretization of the continuous gradient flow (2.1) [33]:

[TABLE]

Combining these two approaches, the proximal gradient method considers a bi-objective minimization $\min_{\mathbf{x}}C_{1}(\mathbf{x})+C_{2}(\mathbf{x})$, where $C_{1}$ and $C_{2}$ are typically assumed to be convex, and $C_{2}(\mathbf{x})$ is usually a regularization in the unconstrained minimization of $C_{1}(\mathbf{x})$. The proximal gradient method then defines iterations [39]

[TABLE]

Expressing the $\operatorname*{arg\,min}_{\mathbf{v}}$ in (2.8) as finding a zero gradient (again either by assumption of convexity or relaxing the minimization to finding a local minimum) yields $\mathbf{x}^{k+1}=\{\mathbf{v}\textnormal{ : }\mathbf{v}+\gamma_{k}\nabla C_{2}(\mathbf{v})=\mathbf{x}^{k}-\gamma_{k}\nabla C_{1}(\mathbf{x}^{k})\}.$ This exactly corresponds to a first-order implicit-explicit formulation (also referred to as a forward-backward splitting) applied to the continuous gradient flow (2.1) for $C(\mathbf{x})=C_{1}(\mathbf{x})+C_{2}(\mathbf{x})$:

[TABLE]

Note that, when $C_{2}\equiv 0$, this reduces to gradient descent for $C_{1}$ (2.2), and when $C_{1}\equiv 0$, this reduces to the proximal point algorithm for $C_{2}$ (2.6).

Implicit-explicit and splitting-based time integration schemes have been widely successful for simulation of partial differential equations, particularly when there are multiple coupled scales or physics (e.g., [6]). Broadly, explicit schemes (like gradient descent) can be intractably slow due to stability restrictions on the stepsize, while implicit schemes (like the proximal point method) are too expensive due to the required nonlinear implicit solve at each step. Implicit-explicit (IMEX) methods attempt to find a balance between the two for efficient (or tractable) simulation. The art of IMEX integration in practice is determining an appropriate splitting (or “partitioning” as it is often referred to in the time integration literature) of the problem that ensures stable and accurate integration. The overwhelming majority of literature on splitting in numerical PDEs assumes a linear additive or component splitting of the problem, as in the proximal gradient method (2.9). However, not all problems have a natural linear splitting, something discussed in detail for thermal radiation transport (TRT) and radiation hydrodynamics in [45; 46]. This motivated the development of a new class of time integration methods for nonlinear splitting of operators [15; 51], which follows from the semi-implicit integrators of [11] and a certain approach used for decades in the TRT literature [28]. For the ODE $\dot{\mathbf{y}}=F(\mathbf{y})$, one first defines a nonlinear splitting of $F$, which is a function ${\widetilde{F}}$ of two arguments such that $\widetilde{F}(\mathbf{y},\mathbf{y})=F(\mathbf{y})$. The nonlinearly split Euler method for approximately evolving the the system in time is then

[TABLE]

The timestep (2.10) is explicit in the first argument of $\widetilde{F}$ and implicit in the second argument. Unlike the fully implicit method (2.7), nonlinear splitting allows one to distribute the implicitness in a more computationally efficient fashion. For example, if $F(y)=y^{2}$, we can define ${\widetilde{F}}(u,v)=uv\implies{\widetilde{F}}(y,y)=y^{2}=F(y)$. Then each step of (2.10) requires a linear solve in $v$, but not a nonlinear solve.

If $\widetilde{F}$ is nonlinear in its second argument, then the timestep (2.10) requires a nonlinear solve. By performing only a single Newton iteration as an inexact solve, we obtain the Rosenbrock method

[TABLE]

where $J_{2}[\mathbf{y}_{k}]=D_{2}F(\mathbf{y}_{k},\mathbf{y}_{k})$ denotes the Fréchet derivative of ${\widetilde{F}}$ with respect to the second argument, evaluated at $\mathbf{y}_{k}$ in both arguments. Replacing the fully nonlinear solve with a single linear solve does not affect the order of accuracy and both (2.10) and (2.11) are first-order integrators for approximating the exact flow. Moreover, in (2.11), the nonlinear splitting takes a more subtle role in defining the linear preconditioner $J_{2}$ that is used to stabilize the otherwise explicit update. In optimization, where it is critical to minimize cost per iteration in method design, only requiring a single linear implicit solve is appealing, and later we demonstrate this to be effective on unconstrained nonlinear least squares problems. The unconstrained optimization methods proposed in Section˜3.1 can largely be seen as generalizations or modifications of (2.10).

Lastly, we remark that higher order nonlinearly split methods have been developed for the purposes of numerical PDEs [15; 51], however, it is not clear the benefit of higher order methods for gradient flow interpretations of optimization problems.

3. Gradient-based optimization with nonlinear splitting

3.1. Nonlinear splitting for unconstrained optimization

In this section, we consider unconstrained optimization, i.e., to minimize a function $C:\mathbb{R}^{N}\rightarrow\mathbb{R}$, which we assume to be $C^{1}$. The gradient descent (GD) method is given by (2.2) with stepsize $1/L\mapsto\gamma_{k}>0$, and the method is specified with an initial iterate $\mathbf{x}^{0}\in\mathbb{R}^{N}$. We now introduce the notion of a nonlinear splitting of the gradient.

Definition 3.1**.**

A nonlinear splitting of the gradient of $C:\mathbb{R}^{N}\rightarrow\mathbb{R}$ (or simply, a nonlinearly split gradient) is a mapping $\widetilde{\nabla}C:\mathbb{R}^{N}\times\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ such that

[TABLE]

for all $\mathbf{x}\in\mathbb{R}^{N}$. Furthermore, we assume that ${\widetilde{\nabla}}C$ is continuously differentiable in both arguments.

Given a nonlinearly split gradient $\widetilde{\nabla}C$ of $C$, we then define the nonlinearly split gradient descent (NS-GD method) as

[TABLE]

As a first example of a nonlinear splitting of a gradient, we define a linearly-implicit nonlinear splitting. The basic idea is to choose a splitting such that the implicit equation (3.2) for $\mathbf{x}^{k+1}$ for NS-GD only requires a linear solve. This is inspired from the time integration literature, particularly, linearly-implicit Runge–Kutta schemes [16; 10].

Example 3.1** (Linearly-implicit nonlinear splitting of a gradient).**

Suppose that the gradient of $C$ can be expressed in the following form

[TABLE]

where $L:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N\times N}$ is a state-dependent linear operator and $g:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ is a generally nonlinear map. The linearly-implicit nonlinear splitting of the gradient is given by

[TABLE]

We refer to this as a linearly-implicit nonlinear splitting since, in the NS-GD (3.2) method, a linear solve is needed to obtain the updated iterate. Specifically, with this choice of splitting, the NS-GD method (3.2) becomes

[TABLE]

and thus, the updated iterate is obtained by a single linear solve

[TABLE]

For fixed $\mathbf{x}$, viewing $\widetilde{\nabla}C(\mathbf{x},\mathbf{y})$ as a one-form with respect to $\mathbf{y}$, if its exterior derivative is zero, then there is an associated cost function $D_{\mathbf{x}}(\mathbf{y})$ such that ${\nabla_{\mathbf{y}}D_{\mathbf{x}}(\mathbf{y})=\widetilde{\nabla}C(\mathbf{x},\mathbf{y})}$. In the case of a linearly-implicit splitting, this implies that there exist matrices $\mathbf{A}(\mathbf{x})$, $\mathbf{B}(\mathbf{x})$ and vector $\mathbf{b}$ such that

[TABLE]

In this scenario, each step of NS-GD applied to $\widetilde{\nabla}C(\mathbf{x},\mathbf{y})$ is equivalent to solving a new quadratic optimization problem using implicit gradient descent. If the nonlinear splitting is selected so that the minimum of each quadratic function $D(\mathbf{y};\mathbf{x})$ is also the global minimum of $C(\mathbf{x})$ (in general highly non-trivial to do), then NP-GD will escape all local minimum and move towards the global min.

For example, consider the simplified case where a nonlinear split gradient has the form

[TABLE]

where $\odot$ denotes the Hadamard product. Each iteration of NS-GD is equivalent to implicit gradient descent applied to the function

[TABLE]

For all $\mathbf{x}$ for which $a(\mathbf{x})_{i}>0$, the iteration will move the solution towards zero. An example of this phenomena can be seen in our numerical experiments on the Rastrigin function in Section˜4.1.1 where NS-GD obtains the global minimum in less than 10 iterations.

As we previously alluded to in Section˜1, the nonlinear splitting framework defines a broad class of gradient-based descent methods that include existing methods, for example linear additive splitting such as in the proximal gradient method as discussed in Section˜2. As another example, the discrete gradient method [20] can be interpreted as a nonlinear splitting gradient descent method.

Example 3.2** (Discrete gradient methods).**

Following [20], a discrete gradient is precisely a nonlinear splitting of the gradient ${\widetilde{\nabla}}C(\mathbf{x},\mathbf{y})$ with an additional requirement that the splitting satisfies the mean value property

[TABLE]

where $\langle\cdot,\cdot\rangle$ denotes the standard inner product on $\mathbb{R}^{N}$. The additional requirement (3.6) ensures that the nonlinear splitting gradient descent step (3.2) is monotone decreasing in the cost, i.e.,

[TABLE]

for any stepsize $\gamma_{k}\geq 0$.

An example of a discrete gradient is the average vector field of $\nabla C$ defined by

[TABLE]

$\widetilde{\nabla}C(\mathbf{x},\mathbf{y})$* is interpreted as the average of the gradient over the line segment connecting $\mathbf{x}$ and $\mathbf{y}$.*

3.1.1. Acceleration

We discuss two types of acceleration of the NS-GD method (3.2), by using momentum and by viewing the NS-GD method as a fixed point iteration.

Momentum. The Nesterov accelerated gradient (NAG) method [37; 38] is given by

[TABLE]

where $\mu>0$ is the momentum factor and $\gamma_{k}>0$ is again the stepsize; the method is initialized with an initial iterate $\mathbf{x}^{0}\in\mathbb{R}^{N}$ and an initial step direction $\mathbf{v}^{0}=0\in\mathbb{R}^{N}$. The NAG method (3.9) can be interpreted as first doing a “look-ahead" step $\mathbf{x}^{k}\mapsto\mathbf{x}^{k}+\mu\mathbf{x}^{k}$, followed by a gradient step, with the gradient evaluated at the look-ahead point, $\nabla C(\mathbf{x}^{k}+\mu\mathbf{v}^{k})$. We express NAG in the form (3.9) since the generalization to a nonlinearly splitting form is straightforward, although there are other equivalent formulations.

By replacing the gradient appearing in NAG by the nonlinearly split gradient, the nonlinearly split NAG (NS-NAG) method is given by

[TABLE]

Fixed-points and linear preconditioning from a nonlinear splitting. Another approach for accelerating gradient-based descent is to view GD as a fixed-point iteration $\mathbf{x}^{k+1}=G(\mathbf{x}^{k})$, with fixed-point operator

[TABLE]

A fixed-point of $G$ then corresponds to a stationary point of $C$. One can also generalize and consider a fixed-point operator with linear preconditioning of the gradient

[TABLE]

where we assume $A(\mathbf{x})$ is invertible. Fixed-point iterations of $G_{A}$ correspond to the linearly preconditioned gradient descent method (2.3). Note that the fixed-points of $G_{A}$ are precisely the same as the fixed-points of $G$; however, the convergence properties of the fixed-point iterations are generally different, and $A(\mathbf{x})$ is usually chosen to improve convergence and/or stability [34; 18].

As discussed in Section˜2, one can consider a linearized solve of a nonlinear splitting method resulting in a linearly-implicit Rosenbrock method (2.11). Here, we will show that such an approach for a NS-GD step can be interpreted as a linear preconditioning of gradient descent. Performing a single Newton iteration, linearized about $\mathbf{x}^{k}$, for the nonlinear solve required in the NS-GD step (3.2) yields precisely a linearly preconditioned gradient descent method,

[TABLE]

where $\widetilde{G}$ denotes the corresponding fixed-point operator. Note that the preconditioner $A(\mathbf{x}):=I+\gamma D_{2}{\widetilde{\nabla}}C(\mathbf{x},\mathbf{x})$ is defined from the choice of nonlinear splitting. We refer to (3.12) as NS-GD-Newton(1).

Naturally, one can consider accelerated variants of (3.12). For example, using Nesterov momentum, an accelerated method can be constructed by instead linearizing the Newton iteration about the look-ahead point $\overline{\mathbf{x}}^{k}\coloneqq\mathbf{x}^{k}+\mu\mathbf{v}^{k}$, which yields

[TABLE]

We refer to (3.13) as NS-NAG-Newton$(1)$. As before, the fixed-points of the iteration (3.13) are the same as those of $G$ and $\widetilde{G}$. Alternatively, one can consider acceleration methods applied to the nonlinear splitting fixed-point operator $\widetilde{G}$ defined in (3.12). For example, a common method for accelerating fixed-point iterations is Anderson Acceleration (AA) [40; 54; 44]. Let $g(\mathbf{x}):=\widetilde{G}(\mathbf{x})-\mathbf{x}$, $\mathbf{x}^{1}=\widetilde{G}(\mathbf{x}^{0})$, and $m\geq 1$ be an integer. AA applied to the fixed-point operator $\widetilde{G}$, for $k=1,2,\dots$, is given by

[TABLE]

We refer to (3.14) as NS-AA-Newton(1). Compared with NS-GD-Newton(1) (3.12), NS-AA-Newton(1) has an additional least squares solve per iteration; the least squares solve in (3.14c) is assumed to be very low dimensional compared to the parameter space, $m\ll N$, and thus the total cost of the AA update remains $\mathcal{O}(N)$ for $\mathbf{x}\in\mathbb{R}^{N}$, a cost that becomes marginal with more complex preconditioners or gradient computations.

We conclude this section with another example, where we show that, for nonlinear least squares with a particular choice of nonlinear splitting, NS-GD-Newton$(1)$ is equivalent to the Levenberg–Marquardt (LM) algorithm [30; 35].

Example 3.3** (Levenberg–Marquardt algorithm as NS-GD-Newton$(1)$).**

Consider the nonlinear least-squares problem, given by

[TABLE]

where $\mathbf{v}\in\mathbb{R}^{N}$ is fixed, $\mathbf{f}:\mathbb{R}^{M}\rightarrow\mathbb{R}^{N}$ is a differentiable and generally nonlinear mapping, and the residual is defined as $\mathbf{r}(\mathbf{x})\coloneqq\mathbf{f}(\mathbf{x})-\mathbf{v}$. Denoting the Jacobian of $\mathbf{r}$ (equivalently, of $\mathbf{f}$) as

[TABLE]

the Levenberg–Marquardt algorithm is given by

[TABLE]

where $\mathbb{J}_{k}\coloneqq\mathbb{J}(\mathbf{x}^{k})$ denotes the Jacobian of the residual at the $k^{th}$ iterate and $\sigma_{k}>0$ is a damping parameter.

Now, we show that the LM algorithm is equivalent to NS-GD-Newton$(1)$ with the following choice of nonlinear splitting, after an appropriate rescaling. The gradient of the nonlinear least squares cost is

[TABLE]

We define a nonlinear splitting of this gradient by

[TABLE]

Then, $D_{2}{\widetilde{\nabla}}C(\mathbf{x},\mathbf{y})=\mathbb{J}(\mathbf{x})^{T}\mathbb{J}(\mathbf{y})$. In this case, (3.12) becomes

[TABLE]

Dividing both sides by $\gamma_{k}>0$ and defining $\sigma_{k}=1/\gamma_{k}$ yields precisely the LM algorithm (3.16).

If we consider instead NS-NAG-Newton$(1)$, we have with this choice of nonlinear splitting

[TABLE]

where $\overline{\mathbf{x}}^{k}=\mathbf{x}_{k}+\mu\mathbf{v}_{k}$ again denotes the look-ahead point and $\overline{\mathbb{J}}_{k}=\mathbb{J}(\overline{\mathbf{x}}^{k})$ denotes the Jacobian linearized about the look-ahead point. This reproduces the “inertial Levenberg–Marquardt" algorithm recently proposed in [29]. Similarly, we can view LM as a fixed-point algorithm and accelerate it with AA. We will compare NS-GD-Newton$(1)$, NS-NAG-Newton$(1)$, and NS-AA-Newton$(1)$ for a nonlinear least squares problem in the second numerical example, Section˜4.1.2.

3.2. Nonlinear splitting with equality constraints

Throughout, $X$ and $\Theta$ will denote vector spaces. Additionally, we assume that $X$ is a Hilbert space with inner product $\langle\cdot,\cdot\rangle$.

Consider the following equality-constrained optimization problem

[TABLE]

where the optimization parameter $\theta\in\Theta$, the state vector $\mathbf{x}\in X$, the objective function $J:X\times\Theta\rightarrow\mathbb{R}$, and the equality constraint $F:X\times\Theta\rightarrow X$.

Analogous to the unconstrained case discussed in Section˜3.1, we are interested in a nonlinear splitting of the gradient for this problem, We denote the gradient of this equality-constrained problem by $\nabla_{\theta}J(\mathbf{x};\theta)$ where $\mathbf{x}$ satisfies $F(\mathbf{x};\theta)=0$. More precisely, suppose that there exists a parameter-to-state (also known as control-to-state) mapping $\theta\mapsto\mathbf{x}[\theta]$ such that $\mathbf{x}[\theta]$ is a solution to $F(\mathbf{x}[\theta];\theta)=0$. Then, the gradient is

[TABLE]

A standard way to compute this gradient, without needing such a parameter-to-state mapping, is through the adjoint method. Let us briefly recall the standard adjoint method [53; 22; 27; 31; 21; 9] for the equality-constrained optimization problem (3.19). The adjoint approach computes the gradient by solving the following system

[TABLE]

where $\partial/\partial\theta$ denotes the derivative with respect to explicit dependence on $\theta$, and the adjoint variable $\boldsymbol{\lambda}\in X$111The adjoint variable is more naturally interpreted as an element of the dual space $X^{*}$; we use the inner product on $X$ throughout to identify the dual space with $X$. For a more detailed discussion of duality in the context of adjoint equations, see [50; 52]. satisfies the adjoint equation (3.20b), $(\cdot)^{*}$ denotes the adjoint with respect to the inner product.

We will now introduce the notions of nonlinear splittings of the equality constraint and of the gradient in the equality-constrained setting, and provide a systematic construction of the nonlinear splitting using a modified adjoint approach.

Definition 3.2** (Nonlinear splitting of an equality constraint).**

A nonlinear splitting of the equality constraint is a mapping

[TABLE]

such that ${\widetilde{F}}(\mathbf{x},\mathbf{x};\theta)=F(\mathbf{x};\theta)$ for all $\mathbf{x}\in X$ and $\theta\in\Theta$. In particular, this implies

[TABLE]

For each $\theta\in\Theta$, we define $K^{\widetilde{F}}(\theta)$ as the zero level set of ${\widetilde{F}}(\cdot,\cdot;\theta)$, i.e.,

[TABLE]

We now define the notion of a nonlinear splitting of the gradient for the equality-constrained optimization problem (3.19).

Definition 3.3** (Nonlinear splitting of the gradient for problem (3.19)).**

Given a nonlinear splitting of the equality constraint ${\widetilde{F}}$, a corresponding nonlinear splitting of the gradient for the problem (3.19) is given by a mapping

[TABLE]

such that, if $\mathbf{x}=\mathbf{y}$, then

[TABLE]

Again, $\nabla_{\theta}J(\mathbf{x};\theta)$ denotes the gradient of problem (3.19) where $\mathbf{x}$ satisfies $F(\mathbf{x};\theta)=0$.

The intuition behind the above definitions is as follows. For standard gradient-based equality-constrained optimization, for a given $\theta$, the derivative corresponding to the problem (3.19) is evaluated at a solution to the equality constraint $\mathbf{x}\in X$ to $F(\mathbf{x};\theta)=0$. By splitting the equality constraint, we have enlarged the space on which we can evaluate the “gradient”. Hence, given a nonlinear splitting of the gradient, we can consider descent algorithms which update the state vector $\mathbf{x}$ and the parameters $\theta$, without requiring that the state vector solves the equality constraint at each iteration. Note that when $\mathbf{x}=\mathbf{y}$, the equality constraint $F(\mathbf{x};\theta)=0$ is satisfied if and only if $(\mathbf{x},\mathbf{x})\in K^{\widetilde{F}}(\theta)$. We require (3.24) to further ensure consistency of the split gradient, reducing to the true gradient if $\mathbf{x}=\mathbf{y}$.

Nonlinear splitting from a nonlinear solver. It is often the case for high-dimensional constrained optimization (e.g., PDE-constrained optimization) that solving the equality constraint at each optimization iteration is a major computational bottleneck. A broad class of numerical methods for solving such equality constraints can be interpreted as nonlinear solvers, as discussed in [14], which consist of applying an iterative fixed-point method to find the solution of the equality constraint. To mitigate this computational bottleneck, we discuss how a single iteration of such a nonlinear solver gives rise to a nonlinear splitting and hence, yields an efficient modified adjoint method to compute the gradient for the constrained objective.

Consider the solution of the equality constraint $F(\mathbf{x};\theta)=0$ via a nonlinear solver. The nonlinear solver framework of [14] is a general framework which includes many common iterative methods for solving nonlinear equations, including, for example, fixed-point iterations, line searches, nonlinear multigrid, and Newton–Krylov methods. For a given $\theta$, we express one iteration of the nonlinear solver as $\mathbf{x}^{k+1}=G(\mathbf{x}^{k};\theta)$. We assume that a solution of the original equality constraint corresponds to a fixed-point of the nonlinear solver $G$, which covers most of the examples of nonlinear solvers discussed in [14]. We emphasize that $G$ is abstract and very general as an operator – it can be include complex linear and nonlinear relations and solves within its action, and simply represents the output of a single iteration of the nonlinear fixed-point solver $G$. Note that this also includes the limiting case of $G$ being an exact solve of the equality constraint, wherein the nonlinearly split method reduces to a classical adjoint. Thus, we can redefine the equality constraint to be $F(x;\theta):=G(x;\theta)-x$. Then, define a nonlinear splitting of the equality constraint ${\widetilde{F}}$ of $F$ by

[TABLE]

We modify the adjoint system (3.20) by using a single iteration of the nonlinear solver corresponding to $G$ as the state variable data for the adjoint equation. It is straightforward to verify that this defines a nonlinear splitting of the equality constraint. In particular, this gives the following nonlinear splitting adjoint gradient (NS-Adj-Gradient):

[TABLE]

where at each iteration, the equations are computed from top to bottom. Note that the nonlinear splitting of the gradient is subtle, with the term $\partial F/\partial\theta$ depending on both $\mathbf{x}^{k}$ and $\mathbf{x}^{k+1}$ because $\mathbf{x}^{k+1}$ does not satisfy the equality constraint exactly. This algorithm is summarized in Algorithm˜1.

Remark 3.1**.**

Algorithm˜1* returns the updated state, updated iterate and the nonlinearly split gradient; we utilize all outputs when they are wrapped into accelerated methods, discussed in Section˜3.2.1.*

More generally, one can consider any nonlinear splitting (not necessarily one defined from a nonlinear solver). Given any nonlinear splitting of the equality constraint, ${\widetilde{F}}$, we define a nonlinear splitting of the gradient by

[TABLE]

It is straightforward to verify that this construction satisfies the definition of a nonlinear splitting of the gradient (Definition˜3.2). To see this, if $\mathbf{y}=\mathbf{x}$, we have $0={\widetilde{F}}(\mathbf{x},\mathbf{x};\theta)=F(\mathbf{x};\theta)$, and the modified adjoint equation for the gradient (3.27c) reduces to the standard adjoint equation (3.20b), since ${\widetilde{F}}(\mathbf{x},\mathbf{x};\theta)=F(\mathbf{x};\theta)$ implies

[TABLE]

Thus, given a nonlinear splitting of the equality constraint, we have shown how to obtain a nonlinear splitting of the gradient using a modified adjoint method.

3.2.1. Acceleration

We consider accelerated gradient descent methods to wrap around the nonlinearly split adjoint gradient Algorithm˜1, using Nesterov accelerated gradient [37; 38] and Anderson acceleration (AA) [40; 54; 44].

For comparison, we first state the adjoint gradient descent algorithm, modified to use the nonlinear split gradient of Algorithm˜1. We now must have two stopping criteria. The first is that the objective is sufficiently small, analogous to standard GD. The second stopping criterion enforces that the equality constraint is satisfied to sufficient tolerance by checking the residual of the equality constraint is sufficiently small. This is done by noting at the $k^{th}$ iterate of NS-GD, we have

[TABLE]

The NAG variant of the nonlinear splitting adjoint gradient descent is given in Algorithm˜3.

The AA variant is given by viewing the gradient descent step

[TABLE]

as a fixed-point iteration. A fixed-point then corresponds to a vanishing (nonlinear splitting) gradient. The AA variant is given in Algorithm˜4.

Remark 3.2**.**

In Anderson acceleration, it is common to use a (economy) QR decomposition to solve $G\xi=g^{k}$, interpreted in the least squares sense. Note that the QR decomposition does not need to be recomputed at each iteration, but rather updated when a column of $G$ is added or removed. For more details, see [40; 54; 44].

In the third numerical example, Section˜4.2.1, we will compare these methods to their standard counterparts, for a source optimization problem where the equality constraint is given by a nonlinear neutron transport equation.

4. Numerical examples

We consider several numerical examples that demonstrate the nonlinear splitting framework. In Section˜4.1, we demonstrate nonlinear splitting methods on examples in unconstrained optimization; particularly, in Section˜4.1.1, we consider several benchmark nonconvex functions and in Section˜4.1.2, we consider a nonlinear least squares problem arising in pure state quantum tomography. Subsequently, in Section˜4.2, we demonstrate the nonlinear splitting methods on examples in equality-constrained optimization; particularly, we consider the stationary neutron transport equation as the equality constraint in Section˜4.2.1 and an initial value problem as the equality constraint in Section˜4.2.2.

For convenience, we summarize the methods introduced thus far in Table˜1.

4.1. Unconstrained examples

4.1.1. Benchmark nonconvex optimization problems

As examples of challenging unconstrained non-convex optimization using nonlinear splitting, we consider several two-dimensional non-convex benchmark problems, namely, the Rastrigin, Rosenbrock, and Beale functions, which all have a global minimum value of zero, and complex nonconvex loss landscapes. The functions and their global minima are shown in Figure˜1. For each of these functions, we define a linearly-implicit nonlinear splitting (Example˜3.1), and we compare NS-GD and NS-NAG based on these splittings to GD, NAG, discrete gradient NS-GD, discrete gradient NS-NAG, and to the built-in MATLAB [24] optimizers fmin (FM), Pattern search (PS), Genetic algorithm (GA), Particle swarm (PartS), Surrogate (Surr), Global search (GS).

The Rastrigin function, its gradient, and a linearly-implicit nonlinear splitting of the gradient are given by

[TABLE]

The Rosenbrock function, its gradient, and a linearly-implicit nonlinear splitting of the gradient are given by

[TABLE]

The Beale function, its gradient, and a linearly-implicit nonlinear splitting of the gradient are given by

[TABLE]

In the following, NS-GD and NS-NAG refer to the nonlinear splitting methods with the above linearly-implicit splittings, whereas NS-GD(AVF) and NS-NAG(AVF) refer to the nonlinear splitting methods using the average vector field (Example˜3.2) as the splitting.

Figures˜2, 3 and 4 show the cost as a function of stepsize and iteration for these three different problems. For all three cases, classical GD and NAG are unable to achieve $<\mathcal{O}(1)$ loss in 10,000 iterations for any stepsizes tested across $\gamma\in[10^{-6},0.1]$ and furthermore, are unstable for larger stepsizes. In contrast, the linearly-implicit methods NS-GD and NS-NAG are stable for larger step sizes for all three problems, and are able to converge to the true minimum. The AVF variants are only successful for the Rosenbrock problem Figure˜3, although there they do converge notably faster than the linearly-implicit methods.

Table˜2 shows a comparison of average runtime and average cost between NS-GD, NS-NAG and built-in MATLAB optimizers. The runs are averaged over 100 random initial conditions, using the same random initial conditions for all methods. For the Rastrigin function, $\mathbf{x}^{0}\in[0,3]\times[0,3]$; for Rosenbrock, $\mathbf{x}^{0}\in[-2,2]\times[-1,3]$; for Beale, $\mathbf{x}^{0}\in[-4,4]\times[-4,4]$. NS-GD and NS-NAG use $\gamma=0.1$ and are run until either 4000 iterations or the cost goes below $10^{-10}$. The most important point from Table˜2 is that the only methods to obtain a small average cost across all three problems (i.e. find the true minimum more or less independent of initial condition) are NS-GD, NS-NAG, and the MATLAB Global Search (GS). In terms of runtime, the NS methods are $6-9\times$ faster for Rastrigin, approximately the same for Beale, and $4-5\times$ slower for Rosenbrock. That being said, note that the NS methods are implemented via naive looping in MATLAB scripts, and we expect a relatively significant performance degradation compared with the optimized methods provided in MATLAB.

4.1.2. Nonlinear least squares: pure state quantum tomography

In this example, we consider a problem in pure state quantum tomography [55], which is to construct a state $|\Psi\rangle$ such that the expectation value of a set of operators $\{A_{i}\}_{i=1}^{N}$ on this state, $\langle\Psi|A_{i}|\Psi\rangle$, agree with a set of observed values $\{a_{i}\}_{i=1}^{N}$. We approach this problem as an unconstrained state reconstruction. This can be formulated as a nonlinear least squares problem of the form

[TABLE]

For this example, we consider $n$-qubit states

[TABLE]

where $|0\rangle$ and $|1\rangle$ denote the standard single qubit ground and excited states, and $c_{i_{1}\dots i_{n}}\in\mathbb{C}$. The dimension of this vector space is $2^{n}$ as a complex vector space and $2^{n+1}$ as a real vector space. A basis for the space of linear operators on this space can be constructed by considering the extended Pauli basis

[TABLE]

A basis for operators on $n$-qubit states is then given by $n$-fold tensor products of the extended Pauli basis, i.e.,

[TABLE]

For this example, we consider $n=6$ which corresponds to a real vector space of dimension $2^{6+1}=128$ with $4^{6}=4098$ operators defined by (4.3). To generate the values $a_{i}$ for the problem (4.1), we constructed a random quantum state $|\Psi_{*}\rangle$ and computed its expectation value with all 4098 operators. In Figure˜5, we compare NS-GD-Newton(1), NS-NAG-Newton(1), and NS-AA-Newton(1), using the nonlinear least squares splitting defined in Example˜3.3, to their standard counterparts GD, NAG, and AA, all using stepsize $\gamma=10$. In Figure˜5(A), we compare the cost versus iteration of these methods with 100 randomly chosen operators. Similarly, in Figure˜5(B), we compare these methods for 500 randomly chosen operators. In Figure˜5, we compare the fidelity of the optimized state after 100 iterations versus the number of operators chosen, where the fidelity of the optimized state $|\Psi\rangle$ versus the prepared state $|\Psi_{*}\rangle$ is defined as

[TABLE]

with a fidelity of $1$ corresponding to equivalent states. In Figure˜5, standard GD and NAG are immediately unstable. On the other hand, all of the nonlinear splitting methods are stable. Furthermore, both NS-GD-Newton(1) and NS-AA-Newton(1) converge to the true minimum of the cost. For sufficiently many operators, the fidelity of the reconstructed state is close to 1 for both NS-GD-Newton(1) and NS-AA-Newton(1). For a lower number of operators, despite both NS-GD-Newton(1) and NS-AA-Newton(1) reaching the true minimum of the cost, the fidelity is relatively compromised due to random selection of operator subset without any systematic selection scheme.

4.2. Examples with equality constraints

4.2.1. Nonlinear neutron transport equation

As an example of equality-constrained optimization using nonlinear splitting, we consider a source optimization problem for the nonlinear neutron transport equation. This source optimization problem arises in optical molecular imaging [47] and bioluminescence tomography [8]; well-posedness of the optimization problem with linear transport is investigated in [8; 47; 23].

The equality constraint for this optimization problem is given by the following nonlinear neutron transport equation

[TABLE]

where ${\boldsymbol{\Omega}}\in\mathbb{S}^{2}$ is the transport angle, $d\Omega$ is the standard measure on $\mathbb{S}^{2}$, $D\ni x$ is the spatial domain, $q(x,{\boldsymbol{\Omega}})$ is the source, $\psi_{\textup{in}}$ is the inflow boundary data, $\sigma_{t}(x)$ is the total cross section, and $\sigma_{s}(x,\phi(x))$ is the nonlinear scattering cross section which depends on $x$ and the scalar flux $\phi$, defined as the zeroth angular moment of $\psi$, i.e.,

[TABLE]

Due to the dependence of the scattering cross section on $\phi$ (and hence, $\psi$), (4.4) is a nonlinear equality constraint.

We consider the following equality-constrained optimization problem

[TABLE]

where the norm is induced from the inner product

[TABLE]

i.e., to find a source $q(x,{\boldsymbol{\Omega}})$ such that the solution to the nonlinear transport equation (4.4) matches a target state $\psi_{\text{tar}}$. Additional details on the adjoint equation of the neutron transport equation for computing the gradient are provided in Appendix˜A (see also [42]).

For this example, we take $D=[0,0.5]\times[0,1]$ using a uniform $24$-by-$48$ grid, $\sigma_{t}(x)=2$ for all $x$, $\sigma_{s}(x,\phi(x))=1+0.02\phi(x)^{2}$, and inflow boundary data

[TABLE]

where $x=(x_{1},x_{2})$ denote the spatial coordinates and ${\boldsymbol{\Omega}}=(\Omega_{1},\Omega_{2},\Omega_{3})$ denote the angular coordinates using the embedding $\mathbb{S}^{2}\hookrightarrow\mathbb{R}^{3}$. The target state $\psi_{\text{tar}}$ was generated by solving the nonlinear neutron transport equation (4.4) with the above boundary data, cross sections, and a source given by

[TABLE]

The spatial discretization is constructed with the finite element library MFEM [5; 36] using linear $L^{2}$ discontinuous finite elements, with $24\times 48\times 4=4608$ spatial degrees of freedom. The angular dependence is discretized using the $S_{N}$ method (see, e.g., [17; 49]) with the Level Symmetric $S_{20}$ angular quadrature rule which has 220 discrete directions, yielding $4608\times 220=1013760$ degrees of freedom for functions of space and angle, particularly for the state vector $\psi$ and the optimization parameter (source) $q$.

The nonlinear splitting of the equality constraint is given by linearizing the scattering cross section about the previous optimization iterate, i.e.,

[TABLE]

We compare NS-Adj-GD (Algorithm˜2), NS-Adj-NAG (Algorithm˜3), NS-Adj-AA (Algorithm˜4) with their standard counterparts Adj-GD, Adj-NAG, and Adj-AA, where the standard methods solve the nonlinear equality constraint to a tolerance of $10^{-6}$ at each optimization iteration. A constant $\gamma=5\times 10^{-9}$ is used for all optimization methods. NS-Adj-NAG and Adj-NAG both use a momentum factor $\mu=0.9$. NS-Adj-AA and Adj-AA both use $m=25$. All runs are initialized with zero initial state, $\psi^{0}=0$, and zero initial source, $q^{0}=0$. Figure˜6(A) shows the relative objective $J/J^{0}$ versus iteration, where $J^{0}$ is the objective value with the initial state and source $\psi^{0}=0=q^{0}$. An efficiency plot is shown in Figure˜6(B), which compares the relative objective versus number of linear transport solves in the forward pass; in particular, the nonlinear splitting variants are able to achieve the same cost as the standard methods with $2-3\times$ less linear transport solves in the forward pass. The equality constraint residual $\|F\|$ versus iteration is shown in Figure˜6(C).

A comparison of the scalar flux $\phi$ of the target state versus the states obtained by NS-Adj-GD, NS-Adj-NAG, and NS-Adj-AA (using the lowest values of the objective function in Figure˜6) is shown in Figure˜7. Similarly, a comparison for the magnitude of the current, $\mathbf{j}$, defined as the first angular moment of $\psi$, i.e.,

[TABLE]

is shown in Figure˜8.

4.2.2. Splitting of the initial condition for an initial value problem

For our final example, we will consider an equality-constrained optimization problem, where the equality constraint is given by a simple linear ODE (more precisely, a time integration method applied to the ODE; for simplicity, we will use backward Euler). We will show how the nonlinear splitting framework can be used to introduce momentum into the initial condition for an optimization problem constrained by an initial value problem, and demonstrate this for a simple linear ODE.

Consider the linear initial value problem

[TABLE]

on a finite-dimensional vector space $X\ni\mathbf{x}(t)$, $A\in L(X,X)$ is a linear map from $X$ to $X$, $\mathbf{f}:[0,T]\rightarrow X$ is a time-dependent external force, and the initial condition $\mathbf{x}_{0}\in X$ is given. As before, we assume that $X$ is equipped with an inner product $\langle\cdot,\cdot\rangle$ and denote the adjoint of the operator $A$ with respect to this inner product by $A^{*}$. Consider the backward Euler method applied to (4.7), using a stepsize $\Delta t=T/N$,

[TABLE]

The backward Euler method (4.8), over all unknowns $\vec{\mathbf{x}}:=\{\mathbf{x}_{1},\dots,\mathbf{x}_{N}\}$, forms the equality constraint under consideration. We will consider an external force optimization problem

[TABLE]

subject to the equality constraint given by the backward Euler method (4.8). The corresponding adjoint equation is

[TABLE]

with unknowns $\vec{\boldsymbol{\lambda}}:=\{\boldsymbol{\lambda}_{0},\dots,\boldsymbol{\lambda}_{N-1}\}$. Now, consider the following family of nonlinear splittings of the equality constraint

[TABLE]

where the family is parametrized by $\alpha\in\mathbb{R}$. Of course, when $\alpha=0$, this reduces to the original equality constraint, and thus, the corresponding nonlinear splitting method (Equation˜4.12 below) reduces to standard adjoint gradient descent. On the other hand, for $\alpha\neq 0$, this can be interpreted as using momentum from the previous optimization iteration to modify the initial condition for each time step. The connection to momentum can be seen by rewriting $(1-\alpha)\mathbf{x}^{k+1}_{n}+\alpha\mathbf{x}^{k}_{n}=\mathbf{x}^{k+1}_{n}-\alpha(\mathbf{x}^{k+1}_{n}-\mathbf{x}^{k}_{n})$, which is the standard initial condition at each timestep $\mathbf{x}^{k+1}_{n}$ minus a momentum correction proportional to $\alpha$. As we will see in the simple example below, the momentum can stabilize the gradient descent. Interestingly, for $\alpha=1$, (4.11) can be computed completely parallel across all timesteps.

For simplicity, we will consider here just adjoint gradient descent for the problem (4.9). The nonlinear splitting adjoint gradient descent for a given $\alpha$, NS-Adj-GD($\alpha$), with this splitting takes the form

[TABLE]

where $\vec{\mathbf{f}}:=\{\mathbf{f}_{1},\dots,\mathbf{f}_{N}\}$ denotes all of the optimization parameters.

Consider the one-dimensional diffusion equation, where $A$ is the standard second-order central difference of the Laplacian on the domain $[0,1]$ with 50 interior nodes, with zero Dirichlet boundary conditions; the problem has an initial condition $\sin(6\pi x)$ and is run to a final time $T=1$ with $N=100$ timesteps. We optimize over the time-dependent external force $\vec{\mathbf{f}}$ initialized at zero, with target terminal state $\sin(2\pi x)$.

In Figure˜9, we compare the objective versus iteration of NS-Adj-GD($\alpha$) for choices of $\alpha=0,0.5,1$ (noting that $\alpha=0$ corresponds to the standard adjoint method). For the non-zero values of $\alpha$, we have included both a “true" and a “proxy" objective value. The proxy value at the $k^{th}$ iterate is simply $J(\mathbf{x}_{N}^{k})$; however, since we are optimizing over the external force $\vec{\mathbf{f}}$, the true value would be given by solving the equality constraint (4.11) with the current external force and then evaluating $J$. In practice, one would use the proxy as opposed to the true cost (since it does not require re-solving the equality constraint exactly), but we have included both, to show that they are in close agreement. As can be seen in Figure˜9, for a smaller stepsize $\gamma=10$, the three methods $\alpha=0,0.5,1$ behave similarly. On the other hand, for a larger stepsize $\gamma=50$, the standard adjoint method $\alpha=0$ becomes unstable, whereas the non-zero $\alpha$ values are able to reach the true minimum. This simple example suggests that using nonlinear splitting, e.g., in the initial conditions as considered here, can be used to stabilize an optimization method. We plan to explore this idea further in future work.

5. Conclusion

In this paper, we introduced a nonlinear splitting framework for gradient-based unconstrained and constrained optimization. In the unconstrained setting, we introduced a nonlinear splitting of the gradient, yielding methods with increased stability, allowing the optimization to explore the loss landscape more globally, as demonstrated by the unconstrained examples. In the constrained setting, we instead applied the nonlinear splitting to the equality constraint, which expands the space on which one can evaluate the gradient of the constrained objective, yielding a more efficient algorithm to obtain an adjoint gradient. We demonstrated the efficacy of this approach on a source optimization problem constrained by a high-dimensional kinetic equation. Interestingly, nonlinear splitting of the equality constraint can also be used to add momentum into the equality constraint; for example, we showed how the nonlinear splitting framework can be used to add momentum in the initial conditions when the constraint is an initial value problem, and provided an example of the stabilizing effect of this added momentum on a linear diffusion equation.

For future research directions, we plan to investigate the convergence and stability of nonlinear splitting optimization methods, noting that, since the nonlinear splitting framework defines a broad family of optimization methods, such analysis would require assumptions to narrow down the class of methods under consideration. It would also be interesting to investigate other applications, both unconstrained and constrained, and develop nonlinear splittings for those applications.

Acknowledgements

BKT was supported by the Marc Kac Postdoctoral Fellowship at the Center for Nonlinear Studies at Los Alamos National Laboratory. BSS was supported by the DOE Office of Advanced Scientific Computing Research Applied Mathematics program through Contract No. 89233218CNA000001, and the National Nuclear Security Administration Interlab Laboratory Directed Research and Development program under project number 20250861ER. DBC was supported by the Center for Nonlinear Studies through the Los Alamos Summer Student Program 2025. SAS was supported by Center for Nonlinear Studies. TB was funded by the National Science Foundation under grant NSF-OIA-2327484. Los Alamos National Laboratory Report LA-UR-25-28812.

Appendix A The adjoint of the neutron transport equation

We provide some additional details on the adjoint of the nonlinear neutron transport equation considered in Section˜4.2.1. The adjoint equation corresponding to the nonlinear neutron transport equation (4.4) is given by

[TABLE]

where $p_{\psi}(x,{\boldsymbol{\Omega}})$ is the adjoint variable, and the scattering cross section is modified for the adjoint equation to be

[TABLE]

Note that each optimization iteration involves an extra linear transport solve for the adjoint equation; these additional adjoint linear transport solves are not included in the efficiency plot, Figure˜6(B), since they are shared across all methods.

Adjoint of inflow boundary conditions. Adjoint boundary conditions are determined by finding necessary and sufficient conditions for boundary integral terms to vanish. The adjoint is chosen to satisfy zero outflow boundary conditions (1.1b) so no integration by parts terms arise. To see this, one linearizes the equality constraint (4.4) with zero inflow boundary conditions $\delta\psi(x,{\boldsymbol{\Omega}})=0$ for all $x\in\partial D$ and ${\boldsymbol{\Omega}}\cdot\mathbf{n}(x)<0$. Given $x\in\partial D$, denote by $\mathbb{S}^{2}_{in}(x)\coloneqq\{\Omega\in\mathbb{S}^{2}:{\boldsymbol{\Omega}}\cdot\mathbf{n}(x)<0\}$; we define $\mathbb{S}^{2}_{out}(x)$ analogously. Then, the boundary term between the adjoint and the linearized equation can be expressed

[TABLE]

where $dS$ denotes the surface measure on $\partial D$. The first term vanishes by the boundary conditions of the variational equation. The second term vanishes for arbitrary $\delta\psi$ on the outflow boundary if and only if the intensity adjoint $p_{\psi}$ vanishes on the outflow boundary.

Adjoint source. The source for the adjoint equation is defined in terms of the objective, the state, and optimization parameters $\theta$. For an optimization problem

[TABLE]

with equality constraint given by (4.4), by setting $g[\psi,\theta]=D_{\psi}J(\psi,\theta)$, the gradient of the objective with respect to $\nu$ can be computed via the adjoint method by

[TABLE]

where the equality constraint is expressed as $F=0$ with

[TABLE]

For the source optimization problem considered in Section˜4.2.1,

[TABLE]

the optimization parameter is the source $\theta=q$. The gradient of the equality-constrained objective with respect to the source is simply (minus) the adjoint variable, which can be seen by

[TABLE]

since $\partial F/\partial q(x,{\boldsymbol{\Omega}})=-\delta(x-x^{\prime})\delta({\boldsymbol{\Omega}}-{\boldsymbol{\Omega}}^{\prime})$.