Glocal Hypergradient Estimation with Koopman Operator

## Abstract

Abstract Gradient-based hyperparameter optimization methods update hyperparameters using hypergradients, gradients of a meta criterion with respect to hyperparameters.
Previous research used two distinct update strategies: optimizing hyperparameters using global hypergradients obtained after completing model training or local hypergradients derived after every few model updates.
While global hypergradients offer reliability, their computational cost is significant; conversely, local hypergradients provide speed but are often suboptimal.
In this paper, we propose glocal hypergradient estimation, blending “global” quality with “local” efficiency.
To this end, we use the Koopman operator theory to linearize the dynamics of hypergradients so that the global hypergradients can be efficiently approximated only by using a trajectory of local hypergradients.
Consequently, we can optimize hyperparameters greedily using estimated global hypergradients, achieving both reliability and efficiency simultaneously.
Through numerical experiments of hyperparameter optimization, including optimization of optimizers, we demonstrate the effectiveness of the glocal hypergradient estimation.

## 1 Introduction

A bi-level optimization problem is a nested problem consisting of two problems for the model parameters ${\bm{\theta}}\in\mathbb{R}^{p}$ and the meta-level parameters called hyperparameters ${\bm{\phi}}\in\mathbb{R}^{q}$ as

$$ $\displaystyle{\bm{\phi}}^{\ast}\in\operatorname*{argmin}_{{\bm{\phi}}}\tilde{ \ell}({\bm{\theta}}^{\ast}({\bm{\phi}});\tilde{\mathcal{D}})$ (1) $\displaystyle\text{such that}\quad{\bm{\theta}}^{\ast}({\bm{\phi}})\in \operatorname*{argmin}_{{\bm{\theta}}}\ell({\bm{\theta}},{\bm{\phi}};\mathcal{ D}).$ (2) $$

Its inner-level problem (Equation 2) is to minimize an inner objective $\ell({\bm{\theta}},{\bm{\phi}};\mathcal{D})$ with respect to ${\bm{\theta}}$ on data $\mathcal{D}$.
The outer-level or meta-level problem (Equation 1) aims to minimize a meta objective $\tilde{\ell}({\bm{\theta}},{\bm{\phi}};\tilde{\mathcal{D}})$ with respect to ${\bm{\phi}}$ on data $\tilde{\mathcal{D}}$, and usually $\mathcal{D}\cap\tilde{\mathcal{D}}=\emptyset$.
A typical problem is hyperparameter optimization (Hutter et al., ), where $\ell$ and $\tilde{\ell}$ are training and validation objectives, and $\mathcal{D}$ and $\tilde{\mathcal{D}}$ are training and validation datasets.
Another example is meta learning (Hospedales et al., ), where $\ell$ and $\tilde{\ell}$ correspond to meta-training and meta-testing objectives.
In the deep learning context, which is the main focus of this paper, ${\bm{\theta}}$ corresponds to neural network parameters, which are optimized by gradient-based optimizers, such as SGD or Adam (Kingma & Ba, ), using gradient $\nabla_{\bm{\theta}}\ell({\bm{\theta}},{\bm{\phi}})$.

Similarly, gradient-based bi-level optimization aims to optimize the hyperparameters with gradient-based optimizers by using hypergradient $\nabla_{\bm{\phi}}\tilde{\ell}({\bm{\theta}}({\bm{\phi}}))$ (Bengio, ; Larsen et al., ).
Although this hypergradient is not always available, when it can be obtained or estimated, gradient-based optimization is efficient and scalable to even millions of hyperparameters (Lorraine et al., ), surpassing the black-box counterparts scaling up to few hundreds (Bergstra et al., , ).
As a result, gradient-based approaches are applied in practical problems that demand efficiency (Choe et al., ), such as neural architecture search (Liu et al., ; Zhang et al., ; Sakamoto et al., ), optimization of data augmentation (Hataya et al., ), and balancing several loss terms (Shu et al., ; Li et al., ).

Hypergradients in previous works can be grouped into two categories: global hypergradients and local hypergradients.
Gradient-based bi-level optimization with global hypergradients uses the gradient of the meta criterion value after the entire model training completes with respect to hyperparameters (Maclaurin et al., ; Micaelli & Storkey, ; Domke, ).
This optimization can find the best hyperparameters to minimize the final meta criterion, but it is computationally expensive because the entire training loop needs to be repeatedly executed.
Contrarily, gradient-based hyperparameter optimization with local hypergradient leverages hypergradients of the loss values after every few iterations of training and updates hyperparameters on-the-fly (Franceschi et al., ; Luketina et al., ).
This strategy can optimize model parameters and hyperparameters alternately and achieve efficient hyperparameter optimization, although the obtained hypergradients are often degenerated because of “short-horizon bias” (Wu et al., ; Micaelli & Storkey, ).

In this paper, we propose *glocal hypergradient* that leverages the advantages and avoids shortcomings of global and local hypergradients.
This method estimates global hypergradients using a trajectory of local hypergradients, which are used to update hyperparameters greedily (see Figure 1).
This leap can be achieved by the Koopman operator theory (Koopman, ; Mezić, ; Brunton et al., ), which linearizes a nonlinear dynamical system, to compute desired global hypergradients as the steady state from local information.
As a result, gradient-based bi-level optimization with glocal hypergradients can greedily optimize hyperparameters using approximated global hypergradients.
We verify its effectiveness in numerical experiments, namely, hyperparameter optimization of optimizers and data reweighting.

Our specific contributions are as follows:

- 1.
We propose the *glocal* hypergradient estimation that combines the advantages of global and local hypergradient techniques.
Specifically, it leverages the Koopman operator theory to predict the global hypergradients from local information, thereby retaining the efficiency of local methods while capturing global information.
- 2.
We provide a comprehensive theoretical analysis of the computational complexities involved in the glocal approach, revealing its computational efficiency compared to the global method.
Additionally, we present an error bound quantifying the accuracy of the proposed glocal estimation relative to the actual global hypergradient.
- 3.
We numerically demonstrate that the proposed glocal hypergradient estimation achieves performance comparable to global approach while maintaining the efficiency of the local method.

The remaining text is organized as follows: Section 2 explains gradient-based bi-level optimization and the Koopman operator theory with related work, Section 3 introduces the proposed glocal hypergradient estimation, Section 4 demonstrates its empirical validity with analysis,
and Section 5 concludes this work.

Figure: Figure 1: The schematic view of hypergradients. We want to use global hypergradient ${\bm{h}}_{T}$ to update hyperparameters, but it needs to wait for the completion of the entire training process. Updating hyperparameters ${\bm{\phi}}_{s}$ with a local hypergradient ${\bm{h}}_{s\tau}$ is efficient but may lead to suboptimal solutions. To leverage both advantages, we propose *glocal* hypergradient estimation that approximates global hypergradient ${\bm{h}}_{T}$ using a local hypergradient trajectory ${\bm{h}}_{t}$ for $t\in[(s-1)\tau+1,\dots,s\tau]$, enabling to update ${\bm{\phi}}_{s}$ using ${\bm{h}}_{T}$ efficiently.

## 2 Background

### 2.1 Gradient-based Bi-level Optimization

Gradient-based bi-level optimization solves the outer-level problem, Equation 1, using gradient-based optimization methods by obtaining hypergradient $\nabla_{\bm{\phi}}\tilde{\ell}({\bm{\theta}}^{\ast}({\bm{\phi}}))$.

When the inner-level problem (Equation 2) involves the training of neural networks, which is our main focus in this paper, computing its minima is infeasible.
Thus, we instead truncate the original inner problem to the following $T$-step optimization process:

$$ $\displaystyle{\bm{\phi}}^{\ast}\in\operatorname*{argmin}_{{\bm{\phi}}}\tilde{ \ell}({\bm{\theta}}_{T}({\bm{\phi}});\tilde{\mathcal{D}})$ (3) $\displaystyle\text{such that}~{}{\bm{\theta}}_{t+1}({\bm{\phi}})=\Theta({\bm{ \theta}}_{t},{\bm{\phi}};\mathcal{D}),$ (4) $$

for $t=0,1,\dots,T-1$.
$\Theta$ is a gradient-based optimization algorithm, such as $\Theta({\bm{\theta}}_{t},{\bm{\phi}};\mathcal{D})={\bm{\theta}}_{t}-\eta\nabla
_{\bm{\theta}}\ell({\bm{\theta}}_{t},{\bm{\phi}};\mathcal{D})$ in the case of vanilla gradient descent, where $\eta$ is a learning rate and can be an element of ${\bm{\phi}}$.

Similarly, we focus on the case that Equation 3 is also optimized by iterative gradient-based optimization

$$ $\displaystyle{\bm{\phi}}_{s+1}=\tilde{\Theta}({\bm{\theta}}_{T},{\bm{\phi}}_{s };\tilde{\mathcal{D}})\quad\text{such that }{\bm{\theta}}_{t+1}({\bm{\phi}})= \Theta({\bm{\theta}}_{t},{\bm{\phi}}_{s};\mathcal{D}).$ (5) $$

where $s=0,1,\dots,S-1$.
The optimization algorithm $\tilde{\Theta}$ in Equation 5 adopts hypergradient $\nabla_{\bm{\phi}}\tilde{\ell}({\bm{\theta}}_{T}({\bm{\phi}}))$, which is referred to as *global* hypergardient (Maclaurin et al., ; Micaelli & Storkey, ; Domke, ).
This design requires $T$-iteration model training $S$ times, which is called non-greedy (Micaelli & Storkey, ).
A greedy approach that alternately updates model parameters and hyperparameters is also possible: hyperparameters are updated every $\tau~{}(<T)$ iteration by hypergradients obtained by a playout until the $T$-th model update $\nabla_{{\bm{\phi}}_{s}}\tilde{\ell}({\bm{\theta}}_{T}({\bm{\phi}}_{s}))$.
In both cases, gradient-based bi-level optimization using global hypergradients requires a computational cost of $O(ST)$, which is computationally challenging for a large $T$.

Instead of waiting for the completion of model training to compute global hypergradient, *local* hypergradient $\nabla_{{\bm{\phi}}_{s-1}}\tilde{\ell}({\bm{\theta}}_{s\tau}({\bm{\phi}}_{s-1}))$ obtained every $\tau$ iteration can also be used to greedily update ${\bm{\phi}}$ (Franceschi et al., ; Luketina et al., ; Wu et al., ).
This relaxation replaces Equations 3 and 4 as

$$ $\displaystyle{\bm{\phi}}_{s+1}=\tilde{\Theta}({\bm{\theta}}_{(s+1)\tau},{\bm{ \phi}}_{s};\tilde{\mathcal{D}})\quad\text{such that }{\bm{\theta}}_{t+1}({\bm{ \phi}})=\Theta({\bm{\theta}}_{t},{\bm{\phi}}_{s};\mathcal{D}),$ (6) $$

for $t\in\mathcal{I}_{s}=[(s-1)\tau,(s-1)\tau+1,\dots,s\tau-1]$.
By setting $\tau\approx T/S$, that is, $S\approx T/\tau$, this approach approximately optimizes ${\bm{\phi}}$ in an $O(T)$ computational cost.
A downside of this approach is that the local hypergradients may be biased, especially when the inner optimization involves stochastic gradient descent (Wu et al., ; Micaelli & Storkey, ).

### 2.2 Computation of Hypergradients

To compute such global or local hypergradients, a straightforward approach is to differentiate through the $T$-step optimization process (Finn et al., ; Grefenstette et al., ; Domke, ).
This unrolling approach is applicable to any differentiable hyperparameters.
Yet, it suffers from large memory requirements, $O(Tp)$, when reverse-mode automatic differentiation is used.
This challenge may be alleviated by using forward-mode automatic differentiation (Micaelli & Storkey, ; Franceschi et al., ; Deleu et al., ).
Moreover, differentiating through long unrolled computational graphs suffers from gradient vanishing/explosion, limiting its applications.

Alternatively, implicit differentiation can also be used with an assumption that ${\bm{\theta}}_{T}$ reaches close enough to a local optimum.
The main bottleneck of this approach is the computation of inverse Hessian with respect to model parameters (Bengio, ), which can be bypassed by iterative linear system solvers (Pedregosa, ; Rajeswaran et al., ; Blondel et al., ), the Neumann series approximation (Lorraine et al., ), and the Nyström method (Hataya & Yamada, ) along with matrix-vector products.
Although this approach is efficient and used in large-scale problems (Choe et al., ; Hataya et al., ; Zhang et al., ), its application is limited to hyperparameters that directly change inner-level loss functions.
In other words, the implicit-differentiation approach cannot be used for other hyperparameters, such as the learning rate of the inner-level optimizer $\Theta$.

The proposed glocal hypergradient estimation can rely on the unrolling approach but differentiating through only $\tau~{}(\ll T)$ iterations; thus, this approach is applicable to various hyperparameters as the unrolling approach while achieving the efficiency like the implicit differentiation approach in terms of being independent of $T$ to compute hypergradients.

### 2.3 Koopman Operator Theory

Here, we roughly introduce the Koopman operator theory.
For a complete introduction, refer to, for example, Brunton et al. (); Kutz et al. ().

Consider a discrete-time dynamical system on $\mathbb{R}^{m}$ represented by $f:\mathbb{R}^{m}\to\mathbb{R}^{m}$ such that ${\bm{x}}_{t+1}=f({\bm{x}}_{t})$ for $t\in\mathbb{Z}$
Then, given a measurement function $g$ in some function space $\mathcal{G}$, the Koopman operator is defined as the following infinite-dimensional linear operator $\mathcal{K}$ such that $\mathcal{K}g=g\circ f$.
In other words, the Koopman operator advances via an observable $g$ the dynamical system one step forward:

$$ $\mathcal{K}g({\bm{x}}_{t})=g(f({\bm{x}}_{t}))=g({\bm{x}}_{t+1}).$ (7) $$

We suppose that $\mathcal{K}$ is invariant in $\mathcal{G}$, i.e., $\mathcal{K}g\in\mathcal{G}$ for any $g\in\mathcal{G}$.
Then, an observation at any time $t$ can be represented with pairs of eigenfunctions $\varphi_{i}:\mathbb{R}^{m}\to\mathbb{C}$ and eigenvalues $\lambda_{i}\in\mathbb{C}$ of the Koopman operator $\mathcal{K}$ as

$$ $\displaystyle g({\bm{x}}_{t})=\mathcal{K}^{t-1}g({\bm{x}}_{0})=\mathcal{K}^{t- 1}\sum_{j}\varphi_{j}({\bm{x}}_{0})v_{j}=\sum_{j}\lambda_{j}^{t-1}\varphi_{j}( {\bm{x}}_{0})v_{j},$ (8) $$

where $v_{j}=\left<\varphi_{j},g\right>(\in\mathbb{C}$), often referred to as a Koopman mode.
For sufficiently large $t$, terms with $|\lambda_{j}|\neq 1$ diverge or disappear.
Consequently, if the dynamics involves no diverging modes, the steady state of $g({\bm{x}}_{t})$ can be written as

$$ $g({\bm{x}}_{t})\xrightarrow[t\to\infty]{}\sum_{j:|\lambda_{j}|=1}\lambda_{j}^{ t-1}\varphi_{j}({\bm{x}}_{0})v_{j}.$ (9) $$

Terms with $|\lambda_{j}|=1$ but $\lambda_{j}\neq 1$ will oscillate in the state space, and those with $\lambda_{j}=1$ will converge to a fixed point.

Numerically, given a sequence ${\bm{g}}({\bm{x}}_{0}),{\bm{g}}({\bm{x}}_{1}),\ldots,{\bm{g}}({\bm{x}}_{t})$ for a set of measurement functions ${\bm{g}}:=[g_{1},\ldots,g_{n}]^{T}$ ($g_{i}\in\mathcal{G})$, Koopman operator $\mathcal{K}$ can be approximated with a finite-dimensional matrix ${\bm{K}}$ by using dynamic mode decomposition (DMD, ) that solves

$$ ${\bm{K}}:=\operatorname*{argmin}_{{\bm{A}}\in\mathbb{C}^{n\times n}}\sum_{t=0} ^{t-1}\left\|{\bm{A}}{\bm{g}}({\bm{x}}_{t})-{\bm{g}}({\bm{x}}_{t+1})\right\|_{ 2}^{2}$ (10) $$

so that $\mathcal{K}g\approx{\bm{c}}^{H}{\bm{K}}{\bm{g}}$ for any $g={\bm{c}}^{H}{\bm{g}}$ (${\bm{c}}\in\mathbb{C}^{n}$).

The Koopman operator theory has demonstrated its effectiveness in the deep learning literature (Hashimoto et al., ; Konishi & Kawahara, ; Naiman & Azencot, ).
In particular, some works used it in the optimization of neural networks (Dogra & Redman, ; Manojlović et al., ) and network pruning (Redman et al., ).
However, these methods require DMD on the high-dimensional neural network parameter space, making its applications to large-scale problems difficult.
Our method also relies on the Koopman operator theory, but we used DMD for the lower-dimensional hypergradients to indirectly advance the dynamics of neural network training, allowing more scalability.

## 3 Glocal Hypergradient Estimation

As explained in Section 2.1, both global and local hypergradient have pros and cons.
Specifically, global hypergradient can optimize the desired objective (Equation 3), but it is computationally demanding.
On the other hand, local gradient can be obtained efficiently, but it may diverge from the final objective.

This paper proposes *glocal* hypergradient that leverages the virtues of these contrastive approaches.
Namely, glocal hypergradient approximates the global hypergradient from a trajectory local hypergradients using the Koopman operator theory to achieve

$$ $\displaystyle{\bm{\phi}}_{s+1}=\tilde{\Theta}(\hat{{\bm{\theta}}}_{\infty}^{(s )},{\bm{\phi}}_{s};\tilde{\mathcal{D}})\quad\text{such that }{\bm{\theta}}_{t} ({\bm{\phi}})=\Theta({\bm{\theta}}_{t-1},{\bm{\phi}}_{s};\mathcal{D}).$ (11) $$

for $t\in\mathcal{I}_{s}$.
$\hat{{\bm{\theta}}}_{\infty}^{(s)}$ is the estimate of the model parameter after the model training $\theta_{T}$ so that $\nabla_{\bm{\phi}}\tilde{\ell}(\hat{{\bm{\theta}}}_{\infty}^{(s)}({\bm{\phi}})
)\approx\nabla_{\bm{\phi}}\tilde{\ell}({\bm{\theta}}_{T}({\bm{\phi}}))$ only by using local hypergradients obtained in $\mathcal{I}_{s}$.
In the remaining text, the superscript indicating the outer time step (^(s)) is sometimes omitted if not confusing for brevity.
Figure 1 schematically illustrates the proposed glocal hypergradient estimation with global and local hypergradients.

To approximate global hypergradient from a trajectory of local hypergradients, we use the Koopman operator theory.
We regard the transition of local hypergradients during $\mathcal{I}_{s}$ as

$$ $\displaystyle{\bm{h}}_{t}\coloneqq\nabla_{{\bm{\phi}}_{s}}\tilde{\ell}({\bm{ \theta}}_{t})=\nabla_{\bm{\theta}}\tilde{\ell}({\bm{\theta}}_{t})\frac{\mathrm {d}{\bm{\theta}}_{t}}{\mathrm{d}{\bm{\phi}}_{s}},$ (12) $$

as a nonlinear dynamical system in terms of $t$.
When using forward mode automatic differentiation that computes Jacobian-matrix product $J_{f}:({\bm{x}},{\bm{V}})\mapsto\partial f({\bm{x}}){\bm{V}}$, where $f:\mathbb{R}^{n\times m}$, ${\bm{x}}\in\mathbb{R}^{n}$, and ${\bm{V}}\in\mathbb{R}^{n\times k}$ for some $n,m,k\in\mathbb{N}$, the right-hand side can be computed iteratively as

$$ $\displaystyle\underbrace{\frac{\mathrm{d}{\bm{\theta}}_{t}}{\mathrm{d}{\bm{ \phi}}_{s}}}_{\eqqcolon{\bm{Z}}_{t}}=\frac{\partial\Theta({\bm{\theta}}_{t-1}, {\bm{\phi}}_{s})}{\partial{\bm{\theta}}_{t-1}}\underbrace{\frac{\mathrm{d}{\bm {\theta}}_{t-1}}{\mathrm{d}{\bm{\phi}}_{s}}}_{={\bm{Z}}_{t-1}}+\frac{\partial \Theta({\bm{\theta}}_{t-1},{\bm{\phi}}_{s})}{\partial{\bm{\phi}}_{s}}=J_{ \Theta(\cdot,{\bm{\phi}}_{s})}({\bm{\theta}}_{t-1},{\bm{Z}}_{t-1})+J_{\Theta( \theta_{t-1},\cdot)}({\bm{\phi}}_{s},{\bm{I}}_{q})$ $$

where ${\bm{Z}}_{1}\in\mathbb{R}^{p\times q}$ is zero.
We assume that there exist a Koopman operator $\mathcal{K}$ and an observable ${\bm{g}}$ that forwards this dynamical system one step as Equation 7, *i.e*., $\mathcal{K}{\bm{g}}(\nabla_{\bm{\phi}}\tilde{\ell}({\bm{\theta}}_{t}))=
\mathcal{K}{\bm{g}}({\bm{h}}_{t})={\bm{g}}({\bm{h}}_{t+1})={\bm{g}}(\nabla_{
\bm{\phi}}\tilde{\ell}({\bm{\theta}}_{t+1}))$, and the elements of ${\bm{g}}$ can be represented by $\mathcal{K}$’s eigenfunctions as in Equation 8.
We also assume that $\mathcal{K}$ can be approximated by a finite-dimensional matrix ${\bm{K}}$ using DMD.

Then, we can estimate the global hypergradient from local hypergradients:

$$ $\displaystyle{\bm{g}}({\bm{h}}_{T})=\mathcal{K}^{T-t}{\bm{g}}({\bm{h}}_{t})= \mathcal{K}^{T-t}\sum_{j}\varphi_{j}({\bm{h}}_{t}){\bm{v}}_{j}\approx{\bm{K}}^ {T-{t}}\sum_{j}b_{j}{\bm{u}}_{j}=\sum_{j}b_{j}{\lambda_{j}}^{T-t}{\bm{u}}_{j},$ (13) $$

where $t\in\mathcal{I}_{s}$, $b_{j}\in\mathbb{C}$ and ${\bm{u}}_{j}\in\mathbb{C}^{n}$ is the $j$-th eigenvector with respect to the eigenvalue $\lambda_{j}$ of ${\bm{K}}$, *c.f*., Equation 8.

If the spectral radius, the maximum norm of eigenvalues, is larger than $1$, the global hypergradient will diverge, suggesting the current hyperparameters are in inappropriate ranges.
Also, if there exists $k$ such that $\lvert\lambda_{k}\rvert=1$ but $\lambda_{k}\neq 1$, the global hypergradient oscillates, suggesting the instability of the current hyperparameter choices.
Thus, we suppose that the spectral radius is not greater than $1$, and if it is $1$, then $\lambda_{k}=1$ for all $k$ and ignore terms with $\lambda_{j}<1$ as they will disappear for sufficiently large $\tau$.
Indeed, this assumption holds in practical cases as shown in Section 4.3.
Subsequently, we approximately obtain glocal hypergradient

$$ $\hat{{\bm{h}}}_{\infty}\coloneqq\nabla_{\bm{\phi}}\tilde{\ell}(\hat{{\bm{ \theta}}}_{\infty})={\bm{g}}^{-1}(\sum_{j:\lambda_{j}=1}b_{j}{\bm{u}}_{j}),$ (14) $$

as Equation 9 and updating hyperparameters by a gradient-based optimizer, such as vanilla gradient descent ${\bm{\phi}}_{s+1}=\tilde{\Theta}(\hat{{\bm{\theta}}}_{\infty}^{(s)},{\bm{\phi}
}_{s};\tilde{\mathcal{D}})={\bm{\phi}}_{s}-\tilde{\eta}\hat{{\bm{h}}}^{(s)}_{\infty}$, where $\tilde{\eta}$ is a learning rate.

In summary, Algorithm 1 shows the pseudocode of the glocal hypergradient estimation with the setting that vanilla gradient descent is used for both inner and outer optimization.

Figure: Algorithm 1 Pseudocode of the glocal hypergradient estimation

Training neural networks for $\tau$ steps, as the function Training in Algorithm 1, requires time $O(\tau p)$ and space $O(p)$ complexities using the standard reverse-mode automatic differentiation.
To compute the hypergradients, we have two options, *i.e*., forward-mode automatic differentiation and reverse-mode automatic differentiation (Baydin et al., ).
The first approach needs the evaluation of the Jacobian-matrix product in 6, resulting in time $O(q\cdot\tau p)$ and space $O(p+\tau q)$ complexities by repeating Jacobian-vector product.
On the other hand, the reverse-mode approach involves $\tau$ evaluations of Training, computed in time $O(\tau\cdot(\tau p+q))$ and space $O(\tau(p+q))$ (Franceschi et al., ).
In the experiments, we adopt forward-mode automatic differentiation to avoid reverse-mode’s quadratic complexities with respect to $\tau$.
However, when the number of hyperparameters $q$ is large, and $\tau$ is limited, reverse-mode automatic differentiation would be preferred.

Additionally, the DMD algorithm at 11 requires space $O(q^{2})$ and time $O(\min(q^{2}\tau,q\tau^{2})+\min(q,\tau)^{3})$ complexities, where each term comes from singular value decomposition (SVD) to solve Equation 10 and eigendecomposition of ${\bm{K}}$.
Compared with the cost to evaluate the Jacobian, these computational costs are negligible.

Table 1 compares the complexities of gradient-based hyperparameter optimization with global, local, and glocal hypergradients, when forward-mode automatic differentiation is adopted for the outer-level problem.
Because the computation of global hypergradients corresponds to the case where $\tau=T$, and it is used $S$ times, its time complexity is $O(STpq)$.
As can be seen, the proposed estimation is as efficient as computing a local hypergradient.

**Table 1: Time and space complexities of gradient-based hyperparameter optimization with global, local, and glocal hypergradients using forward-mode automatic differentiation to compute hypergradients.**
|  | Global | Local | Glocal |
| --- | --- | --- | --- |
| Time | $O(STpq)$ | $O(S\tau pq)$ | $O(S\tau pq)$ |
| Space | $O(p+q)$ | $O(p+q)$ | $O(p+\tau q)$ |

The error of the proposed glocal hypergradient from the actual global hypergradient is bounded as follows.

###### Theorem 3.1 .

Assume that there exists a finite-dimensional Koopman operator $\bar{{\bm{K}}}$ governing the trajectory of hypergradients in each $\mathcal{I}_{s}$, which is approximated with ${\bm{K}}$ using DMD, and the spectral radii of $\bar{{\bm{K}}}$ and ${\bm{K}}$ are 1. Then,

$$ $\displaystyle\lVert{\bm{h}}_{T}-\hat{{\bm{h}}}_{\infty}\rVert_{2}$ $\displaystyle\leq\lVert{\bm{U}}\rVert_{2}\lVert{\bm{U}}^{-1}\rVert_{2}\{\lVert {\bm{e}}_{\tau}\rVert_{2}+(T-s\tau)\varepsilon_{\tau}\}+\sum_{j:\lvert\lambda_ {j}\rvert<1}\lvert b_{j}\rvert\lvert\lambda_{j}\rvert^{\tau}\lVert{\bm{u}}_{j} \rVert_{2},$ (15) $$

where ${\bm{U}}=[{\bm{u}}_{1},\dots,{\bm{u}}_{q}]$, ${\bm{e}}_{\tau}={\bm{h}}_{s\tau}-{\bm{g}}^{-1}\left({\bm{K}}{\bm{g}}({\bm{h}}_
{s\tau-1})\right)$, and $\varepsilon_{\tau}$ is a constant only depends on the number of local hypergradients $\tau$.

###### Proof.

The left hand side of Equation 15 can be decomposed as
$\lVert{\bm{h}}_{T}-\hat{{\bm{h}}}_{\infty}\rVert_{2}\leq\lVert{\bm{h}}_{T}-
\hat{{\bm{h}}}_{T}\rVert_{2}+\lVert\hat{{\bm{h}}}_{T}-\hat{{\bm{h}}}_{\infty}
\rVert_{2}$,
where $\hat{{\bm{h}}}_{T}$ is obtained from the DMD algorithm as
$\hat{{\bm{h}}}_{T}=\sum_{j}b_{j}\lambda_{j}^{T-s\tau}{\bm{u}}_{j}$.
Then, we get

$$ $\lVert{\bm{h}}_{T}-\hat{{\bm{h}}}_{T}\rVert_{2}\leq\lVert{\bm{U}}\rVert_{2} \lVert{\bm{U}}^{-1}\rVert_{2}\{\lVert{\bm{e}}_{\tau}\rVert_{2}+(T-s\tau) \varepsilon_{\tau}\}$ (16) $$

by using Theorem 3.6 of Lu & Tartakovsky (), and

$$ $\displaystyle\lVert\hat{{\bm{h}}}_{T}-\hat{{\bm{h}}}_{\infty}\rVert_{2}=\lVert (\sum_{j}-\sum_{j:\lvert\lambda_{j}\rvert=1})b_{j}{\lambda_{j}}^{T-s\tau}{\bm{ u}}_{j}\rVert_{2}=\lVert\sum_{j:\lvert\lambda_{j}\rvert<1}b_{j}{\lambda_{j}}^{ T-s\tau}{\bm{u}}_{j}\rVert_{2}\leq\sum_{j:\lvert\lambda_{j}\rvert<1}\lvert b_{ j}\rvert\lvert\lambda_{j}\rvert^{\tau}\lVert{\bm{u}}_{j}\rVert_{2},$ (17) $$

since $T-s\tau\geq\tau$.
∎

###### Remark 3.2 .

The first term of Equation 15 decreases as the outer optimization step $s$ proceeds.
In this term, $\lVert{\bm{e}}_{\tau}\rVert$ is minimized by the DMD algorithm as Equation 10.
$\varepsilon_{\tau}$ and the second term of Equation 15 decrease as we use more local hypergradradients $\tau$ for estimation, because ${\bm{K}}$ converges to $\bar{{\bm{K}}}$ as $\tau$ increases (Korda & Mezić, ).
When $T$ is constant, increasing $\tau$ improves the quality of the glocal hypergradient estimation while decreasing the number of the outer optimization step, and we need to take their trade-off in practice.

###### Theorem 3.1 .

###### Proof.

###### Remark 3.2 .

## 4 Experiments

This section empirically demonstrates the effectiveness of the glocal hypergradient estimation.

We implemented neural network models and algorithms, including DMD, using JAX (v0.4.28) (Bradbury et al., ), optax (DeepMind et al., ), and equinox (Kidger & Garcia, ).
Reverse-mode automatic differentiation is used for model gradient computation.
Forward-mode automatic differentiation is adopted for hypergradient computation.

In the following experiments, instead of using the vanilla DMD shown in Equation 10, we adopt Hankel DMD (Arbabi & Mezic, ) that augments the observable ${\bm{g}}$ with time-delayed coordinates as ${\bm{g}}({\bm{h}}_{t})=[{\bm{h}}_{t};{\bm{h}}_{t+1};\dots;{\bm{h}}_{t+m-1}]\in
\mathbb{R}^{mq}$
for some positive integer $m$.

For DMD, we use the last $\sigma$ hypergradients out of $\tau$ hypergradients in each $\mathcal{I}_{s}$ because the dynamics immediately after updating hyperparameters may be unstable.
Throughout the experiments, we optimize model parameters 10k times and hyperparameters every 100 model updates, *i.e*., $\tau=100$, using Adam optimizer with a learning rate of $0.1$ (Kingma & Ba, ).
We set $m=10$ and $\sigma=80$, except for the analysis in Section 4.3.

When computing hypergradients, we use hold-out validation data $\tilde{\mathcal{D}}$ from the original training data, following Micaelli & Storkey ().
To minimize the effect of stochasticity, we use as large a minibatch size as possible for the validation.
Average performance on test data over three different random seeds is reported.
Further experimental details can be found in Appendix B.

We adopt local and global baselines that greedily update hyperparameters using global and local hypergradients, respectively.
As the global baseline is prone to diverge and computationally expensive, we can only present its results on MNIST-sized datasets.

The source code is included in the supplemental and will be released upon publication.

### 4.1 Optimizing Optimizer Hyperparameters

Appropriately selecting and scheduling hyperparameters of optimizers, in particular, learning rates, is essential to the success of training machine learning models (Bergstra et al., ).
Here, we demonstrate the validity of the glocal hypergradient estimation in optimizing such hyperparameters.
Unlike (Micaelli & Storkey, ), we adopt forward-mode automatic differentiation to differentiate through the optimization steps, which allows to use any differentiable optimizers, including SGD and Adam, for the inner optimization without requiring any manual reimplementation.

First, we train LeNet (LeCun et al., ) using SGD with learnable learning rate, momentum, and weight-decay rate hyperparameters, and Adam (Kingma & Ba, ) with learnable learning rate and betas.
The logistic sigmoid function is applied to these hyperparameters to limit their ranges in $(0,1)$.
MNIST variants, namely, MNIST (Le Cun et al., ), Kuzushiji-MNIST (KMNIST, Clanuwat et al. ()) and Fashion-MNIST (FMNIST, Xiao et al. ()) are used as datasets.
The LeNet has 15k parameters and is trained for 10k iterations.

Figures 2 and C.3 shows learning curves with the transition of hyperparameters of optimizers on the FMNIST dataset.
We can observe that the development of hyperparameters of the glocal approach exhibits similar trends with the global baseline, and the glocal and global approaches show nearly identical performance.
Contrarily, the local baseline changes hyperparameters aggressively, resulting in suboptimal performance.
Other results can be found in Appendix C.

Next, we train WideResNet 28-2 (Zagoruyko & Komodakis, ) on CIFAR-10/100 Krizhevsky et al. () using SGD with learnable learning rate and weight-decay rate hyperparameters.
The used model has 1.5M parameters and is trained for 10k iterations.
The results are presented in Figures 3 and C.4
Contrarily to the local baseline drastically changing hyperparameters, the glocal estimation adjusts the hyperparameters gradually and yields better performance, partially because it can predict the future state.

Figure: Figure 2: Test accuracy and the transition of the hyperparameters of SGD and Adam. The proposed local approach shows similar hyperparameter development to the global baseline.
Refer to caption: x1.png

Figure: Figure 3: The transition of the learning rate hyperparameter and test accuracy curves of WideResNet 28-2 on CIFAR-10.
Refer to caption: x2.png

### 4.2 Data Reweighting

Data reweighting task trains a meta module $\mu_{\bm{\phi}}$ to reweight a loss value to each example to alleviate the effect of class imbalanceLi et al. (); Shu et al. ().
$\mu_{\bm{\phi}}:\mathbb{R}\to(0,1)$ is an MLP, and the inner loss function of the task is $\ell({\bm{\theta}};\mathcal{D})=\sum_{({\bm{x}},y)\in\mathcal{D}}L({\bm{x}},y;
{\bm{\theta}})\cdot\mu_{\bm{\phi}}(L({\bm{x}},y;{\bm{\theta}}))$, where $L:\mathcal{D}\to\mathbb{R}$ is cross-entropy loss.
${\bm{\phi}}$ is trained on balanced validation data.

We train WideResNet 28-2 on imbalanced CIFAR-10 and CIFAR-100 Cui et al. (), which simulate class imbalance.
Specifically, the imbalanced data with an imbalance factor of $f$ reduces the number of data in the $c$-th class to $\lfloor\frac{1}{f^{c/C}}\rfloor$, where $C$ is the number of categories, *e.g*., 10 for CIFAR-10.
As the meta module, we adopt a two-layer MLP with a hidden size of 128, consisting of 385 parameters.
Table 2 demonstrates that the glocal approach surpasses the local baseline, indicating its effectiveness to a wide range of problems.

**Table 2: Test accuracy of WideResNet 28-2 trained on imbalanced datasets.**
| Dataset / Imbalance Factor | Local | *Glocal* |
| --- | --- | --- |
| CIFAR-10 / 50 | 0.520 | 0.715 |
| CIFAR-10 / 100 | 0.445 | 0.595 |
| CIFAR-100 / 50 | 0.315 | 0.350 |

### 4.3 Analysis

Below, we analyze the factors of the glocal hypergradient estimation using the task in Section 4.1 on the FMNIST dataset.

Figure: Figure 4: Left The eigenvalues obtained by the Hankel DMD for a hypergradient trajectory. Eigenvalues nearly close to $1$ are highlighted in orange. Middle The magnitude of modes of the estimated hypergradient corresponding to the learning rate $b_{j}\lambda_{j}^{t}u_{j}$ for $j:\lambda_{j}\neq 1$ over 100 iterations of $t$. Right The comparison of validation performance of the proposed estimation with different configurations. LeNet is trained on FMNIST with an initial learning rate of $0.1$.
Refer to caption: x3.png

**Table 3: Runtime comparison of optimizing optimizer hyperparameters with LeNet in second. “w/o HPO” indicates training without hyperparameter optimization process.**
| w/o HPO | Local | Global | *Glocal* |
| --- | --- | --- | --- |
| 11.4 | 58.2 | 2891.2 | 58.3 |

In Section 3, the existence of eigenpairs whose eigenvalues equal to one was assumed for efficient approximation of a global hypergradient.
To see whether it holds in practice, the left panel of Figure 4 shows the eigenvalues obtained by the Hankel DMD at $s=1$.
We can observe two eigenvalues nearly close to $1$ (highlighted in orange).
Additionally, the middle panel of Figure 4 illustrates that modes with other eigenvalues decay rapidly.
Because the magnitudes of the modes associated with the eigenvalue of 1 is an order of $10^{-3}$, other modes can be ignored, numerically supporting the validity of our assumption.

The right panel of Figure 4 compares the proposed method with different configurations to see how the configurations of the Hankel DMD algorithm, specifically, the number of hypergradients ($\sigma$) and the number of stacks per column ($m$), affect the estimation.
We can observe that it performs better when DMD uses the last 80 hypergradients out of 100.
Although the estimation is expected to improve as more hypergradients are used, they may be unstable after the change of hyperparameters, which degenerates the quality of the prediction; discarding the first 20 hypergradients may avoid such an issue.

Table 3 compares runtime on a machine equipped with an AMD EPYC 7543P 32-Core CPU and an NVIDIA RTX 6000 Ada GPU with CUDA 12.3.
Note that because of the limitation of JAX that eigen decomposition for asymmetric matrices on GPU is not yet supported, the computation of the adopted DMD algorithm is suboptimal.
Nevertheless, the proposed method achieves significant speedup compared with the global baseline, revealing its empirical efficiency.

## 5 Discussion and Conclusion

This paper introduced *glocal* hypergradient estimation, which leverages the virtues of global and local hypergradients simultaneously.
To this end, we adopted the Koopman operator theory to approximate a global hypergradient from a trajectory of local hypergradients.
The numerical experiments demonstrated the validity of bi-level optimization using glocal hypergradients.

In this work, we have implicitly assumed that the meta criterion converges after enough iterations, and so do hypergradients.
This assumption may be strong when a certain minibatch drastically changes the value of the meta criterion, such as adversarial learning.
Luckily, we did not suffer from such a problem, probably because we focused our experiments only on supervised learning with sufficiently sized datasets and used full-batch validation data.
Introducing stochasticity by the Perron-Frobenius operator (Korda & Mezić, ; Hashimoto et al., ), an adjoint of the Koopman operator, would alleviate this limitation, which we leave for future research.

## Appendix A Additional Discussion on Theorem 3.1

The first half of the proof of Theorem 3.1 depends on the theorem 3.6 of Lu & Tartakovsky (), which shows Equation 16 without equality holds if $\bar{{\bm{K}}}$’s spectral radius $\rho(\bar{{\bm{K}}})<1$.
This condition can be relaxed to $\rho(\bar{{\bm{K}}})\leq 1$, and then Equation 16 holds.
$\varepsilon_{\tau}$ in Theorem 3.1 can be given as

$$ $\varepsilon_{\tau}=\left\{c_{\tau}\left(\lVert{\bm{h}}_{1}\rVert_{2}^{2}+\sum_ {t=1}^{\tau}\lVert{\bm{f}}_{t}\rVert_{2}^{2}\right)\right\}^{1/2},$ (18) $$

where ${\bm{f}}_{t}=\bar{{\bm{K}}}{\bm{g}}({\bm{h}}_{t})-{\bm{g}}({\bm{h}}_{t+1})$ and $c_{\tau}\geq\lVert\bar{{\bm{K}}}-{\bm{K}}\rVert_{2}^{2}$ depends on the number of hypergradients.

## Appendix B Detailed Experimental Settings

Throughout the experiments, the batch size of training data was set to 128.

### B.1 Optimizing Optimizer Hyperparameters

We set the number of filters in each convolutional layer to 16 and the dimension of the following linear layers to 32. The leaky ReLU is used as its activation.

We modified the original WideResNet 28-2 in Zagoruyko & Komodakis () as follows: replacing the batch normalization with group normalization Wu & He () and adopting the leaky ReLU as the activation function.

We separated 10% of the original training data as validation data.

### B.2 Data reweighting

We modified the original WideResNet 28-2 in Zagoruyko & Komodakis () as follows: replacing the batch normalization with group normalization and adopting the leaky ReLU as the activation function.
The model was trained with an inner optimizer of SGD with a learning rate of 0.01, momentum of 0.9, and weight decay rate of 0.

We separated 1000 data points from the original training dataset to construct a validation set.

## Appendix C Additional Experimental Results

Figures C.2 and C.1 present the test accuracy curves and the transition of SGD’s hyperparameters of KMNIST and MNIST.
As the results presented in the main text in Section 4.1, the proposed glocal approach yields similar behaviors as the global one, while maintaining the efficiency of the local approach.
As can be seen from Figures C.2 and C.3, the global method often fails because of loss explosion, as discussed in the literature (Micaelli & Storkey, ).
Similarly, Figure C.4 shows the results on the CIFAR-100 dataset.

Figure: Figure C.1: The transition of the SGD’s hyperparameters and test accuracy curves of LeNet on MNIST with an initial learning rate of 0.01.
Refer to caption: x5.png

Figure: Figure C.2: The transition of the SGD’s hyperparameters and test accuracy curves of LeNet on KMNIST with an initial learning rate of 0.01.
Refer to caption: x6.png

Figure: Figure C.3: The transition of the SGD’s hyperparameters and test accuracy curves of LeNet on FMNIST with an initial learning rate of 0.1.
Refer to caption: x7.png

Figure: Figure C.4: The transition of the SGD’s hyperparameters and test accuracy curves of WideResNet 28-2 on CIFAR-100.
Refer to caption: x8.png

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: 1. Can the proposed method handle changing dynamics, given that the training dynamics change over time (initially the loss drops significantly and then less so)?
2. For experiments, how do the results change over different initial learning rates? Can Glocal HGD always converge to a good learning rate no matter what the initial learning rate is?
3. Why does Glocal outperform global HGD towards the end of training?

Also see Weaknesses above.
soundness: 2
strengths: The paper studies an important problem, since hyperparameter optimization is a common challenge in training neural nets. The paper's approach using Koopman operator theory is clever and connects hyperparameter optimization to nonlinear dynamical systems. The algorithm and the computational complexities are clearly written, and I appreciate the diagnostic plots in the experimental section.
summary: This paper proposes a hypergradient descent method that uses local hypergradients to estimate the global hypergradient. The method is based on Koopman operator theory which uses a finite-dimensional linearization of the nonlinear dynamics of the hypergradient. Authors characterize the estimation error and provides detailed comparisons on the computational complexities of Glocal HGD, local HGD, and global HGD. Finally, they present experiments that show the practical efficacy of the method.
weaknesses: 1. Important design choices are not given: (1) how should we select the dimension of the Koopman operator $n$? Intuitively, $n$ should depend on properties of the underlying dynamical system, and it would be helpful to have some guidelines. (2) how should we select $\textbf{g}$? Authors use Hankel DMD in the experiments, and it would be great to provide some justification. 

2. Experiments: (1) the experiments are relatively small scale. Authors mention that global hypergradients are difficult to obtain for larger models, but still it would be good to compare Glocal to other baselines, for example the best tuned SGD. Without larger scale experiments, it is difficult to know how the proposed method scales with model size. (2) it would be useful to compare Glocal with the best tuned SGD/Adam as a baseline (3) it would be great to verify empirically that Glocal can indeed estimate the global hypergradient. If I understand correctly, one can compare the estimated gradient with the true global hypergradient. 

3. Theoretical analysis: Theorem 3.1 assumes a Koopman operator of dimension $n$ exists. Without specifications on $n$, this result is less meaningful as $n$ can be very large. However, a large $n$ can make the Glocal method impractical. It would be helpful to provide justification and guidance on $n$.

4. Section 2.3 is not clear: (1) Line 172, what is the relationship between $\phi$ and $\varphi$? The decomposition of $g$ is written in $\phi$ on line 172 but in $\varphi$ in eq. (7). (2) Eq. (8) is not precise, because large $\lambda$'s will diverge, but the arrow notation usually means convergence. (3) The notation for the set of measurement functions $\textbf{g}$ and the individual functions $g_i$ is confusing because they can be very different from the original function $g$ and yet use the same letter.

---

confidence: 3
contribution: 3
initial_rating: 5
presentation: 2
questions: 1. In equation (9), could you provide an analysis of the computational complexity of solving the optimization problem w.r.t $A\in \mathbb{C}^{n\times n}$? Furthermore, are there any potential challenges in its implementation?

2. In equation (13), the computation of $g^\dagger$ involves solving a system of linear equations $g^\dagger g(x) = x$. How does the computation affect the overall complexity of Algorithm 1?

3. In Section 4.1, why is the Global curve missing in Figure 3? Could you discuss any implications this might have for the comparison between methods?

4. How does the proposed approach compare with specific state-of-the-art hyperparameter optimization methods, such as Bayesian optimization or evolutionary algorithms, in terms of performance metrics like accuracy or convergence speed, as well as computational overhead?

There are some typos:
1. In equation (1), the meta-level function is expressed as $\tilde{l}(\theta^*(\Phi);\tilde{D})$. However, the meta objective in line 37 is written as  $\tilde{l}(\theta,\Phi;\tilde{D})$

2. In equation (9), the expression $\sum\limits_{t=0}^{t-1}$ has a conflict because the variable $t$ is used both as the index and the limit of summation.

3. In line 266, "requires time $O(\tau p)$ and space $O(p)$ complexities" should be "requires time complexity of $O(\tau p)$ and space complexity of $O(p)$".

4. The sentence in line 325 seems to be awkwardly structured.
soundness: 2
strengths: 1. The integration of Koopman operator theory to enhance hypergradient estimation is a novel approach, offering a fresh perspective on hyperparameter optimization.

2. The method significantly reduces computational costs compared to traditional global hypergradient approaches.

3. The approach is scalable to large-scale problems, making it applicable to real-world deep learning tasks. Furthermore, the paper provides numerical experiments demonstrating the method's effectiveness in various scenarios.
summary: This paper proposes a novel method called glocal hypergradient estimation for hyperparameter optimization. This method combines the computational efficiency of local hypergradients with the reliability of global hypergradients. It utilizes Koopman operator theory to approximate global hypergradients from the trajectory of local hypergradients, enabling efficient and effective hyperparameter updates. Finally, the paper validates the effectiveness of the proposed method in hyperparameter optimization through experiments.
weaknesses: 1. Algorithm 1 and Theorem 3.1 rely on assumptions about the spectral radius and stability, which may not hold in all cases.

2. The theoretical foundation involving Koopman operators may be complex for practitioners unfamiliar with the concept.

3. The experiments are somewhat limited. Could additional datasets be included, or could comparative experiments be conducted on other models as well?

4. The presentation of the experimental results is somewhat unclear. For example, all the experimental results in Section 4.2 are only summarized in Table 2, with no accompanying figures. Could Table 3 offer a more detailed explanation?

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: Comments & questions
1. I suggest mentioning, at least in passing, the initialization of $\theta_0$ beneath equation (4).
2. Minor notational note: in your problem formulation, you are adopting the convention of outer-level symbol carrying an additional tilde. The parameters are the exception ($\theta$ vs $\phi$); however, for the optimization algorithm, you are using capital $\Theta$ and $\tilde \Theta$. I like having separate symbols for the parameters, so for consistency I would recommend using a capital $\Phi$ instead of $\tilde \Theta$ for the outer-level gradient step. Or find a non-theta symbol to use to denote an optimization step. Just a suggestion, I'll leave it up to the authors.
3. In equation (8), why do terms which diverge ($|\lambda_j| > 1$) no longer appear in the expression? I would think that they would dominate all other modes.
4. For what trajectory $x_t$ is (9) being solved? Wouldn't it make sense to look at many trajectories from different initial conditions?
5. I'm confused about the implications of the bound in Theorem 3.1. The bound aggregates terms corresponding to non-unit eigenvalues, presumably becasue these were discarded previously (see Weakness 1). Wouldn't your theorem suggest that including these terms would result in less estimation error?
6. I don't understand the measurement defined in line 345. How can the measurement function depend on future hypergradients, some of which potentially haven't been computed yet? For example, how would you compute $g(h_{t^*})$ where $t^*$ is the biggest index in $\mathcal{I}_s$?

Minor notes:
* Line 123: "global hypergardient"
* Line 167: the inputs and outputs of $g$ should be clarified; I presume $g: R^m \to R$
* Line 172: mixing up $\phi$ and $\varphi$
* Line 234: "as in Section 2.3"
* Line 334: "outer steps"
* Figure 4: periods after left/middle/right
soundness: 3
strengths: I'm on the fence on this submission -- the method is well-motivated and the derivation is for the most part clear, but I felt that the experimental results are somewhat of a let down.

* The paper is overall well-written, with a few minor typos. The authors do an admirable job of making their theory tractable and easy to read.
* Meta optimization is an important problem and the research setting is well-motivated.
* The authors provide a thorough runtime comparison and associated discussion, which shows that their approach is as computationally efficient as using local hypergradient.
summary: This work considers the problem of hyperparameter optimization for bilevel optimization problems. By leveraging Koopman operator theory, the authors provide a method to estimate global hypergradients (usually only available after fully solving the inner problem) with local hypergradient trajectories (partial solution trajectories of the inner prolbem).
weaknesses: 1. I'm somewhat suspicious of the handwaving around non-unit eigenvalues. Specifically, consider the section from line 246 - 252, where basically all eigenvalues besides those which are equal to one are discarded. Is there any theoretically grounded explanation from doing so? If the koopman operator says that the global hypergradients should oscillate, when intervene and artificially eliminate those modes? Similarly, when solving the DMD for $K$ as in (9), I don't see why you would get modes with an eigenvalue exactly equal to one as there is random noise in the collected hypergradient data.
2. The experiment settings are minimal, and only involve relatively simple image classification tasks. 
3. The experimental results somewhat conflict with the aims of the paper. Namely, the performance of the global hyperparameter estimation degrades severely in the bottom-left panel of Figure 2. In Figures C.2 and C.3 in the appendix are even worse. The appendix briefly mentions "loss explosion," but a further discussion is warranted, since my impression is that the whole point of this paper is to provide a fast way of approximating the "gold-standard" global hypergradients!
4. Additional baselines of estimating the global hypergradient should be considered. At least: simply using the hypergradient at iteration $\tau$ as a substitute for $h_T$ (just once, not continuing to train as in local hypergradient estimation).

Meta-review: This paper proposes a novel method called glocal hypergradient estimation for hyperparameter optimization. It uses Koopman operator theory to approximate global hypergradients from the trajectory of local hypergradients, enabling efficient and effective hyperparameter updates.  

The reviewers generally agreed that the paper was well-written and presented a novel approach to an important problem.  They also found the theoretical contributions to be sound and the experimental results to be promising.  However, there were some concerns about the clarity of the presentation and the limited scope of the experiments.  

Overall, there was consensus among the reviewers that the empirical results are not quite compelling enough to justify publication yet.

Decision: Reject