mode=buildnew

# M-HOF-Opt: Multi-Objective Hierarchical Output Feedback Optimization via
Multiplier Induced Loss Landscape Scheduling

Xudong Sun<EMAIL_ADDRESS>
code: https://github.com/marrlab/DomainLab/tree/mhof Institute of AI for
Health, Computational Health Center, Helmholtz Munich, Munich, Germany Nutan
Chen  Machine Learning Research Lab, Volkswagen Group, Munich, Germany  Alexej
Gossmann  U.S. FDA Center for Devices and Radiological Health
Silver Spring, MD, USA  Yu Xing  KTH Royal Instiute of Technology
Stockholm
Sweden  Carla Feistner  Institute of AI for Health, Computational Health
Center, Helmholtz Munich, Munich, Germany Emilio Dorigatti  Institute of AI
for Health, Computational Health Center, Helmholtz Munich, Munich, Germany
Felix Drost  Institute of AI for Health, Computational Health Center,
Helmholtz Munich, Munich, Germany Daniele Scarcella  Institute of AI for
Health, Computational Health Center, Helmholtz Munich, Munich, Germany Lisa
Beer  Institute of AI for Health, Computational Health Center, Helmholtz
Munich, Munich, Germany Carsten Marr  Institute of AI for Health,
Computational Health Center, Helmholtz Munich, Munich, Germany

###### Abstract

We address the online combinatorial choice of weight multipliers for multi-
objective optimization of many loss terms parameterized by neural works via a
probabilistic graphical model (PGM) for the joint model parameter and
multiplier evolution process, with a hypervolume based likelihood promoting
multi-objective descent. The corresponding parameter and multiplier estimation
as a sequential decision process is then cast into an optimal control problem,
where the multi-objective descent goal is dispatched hierarchically into a
series of constraint optimization sub-problems. The subproblem constraint
automatically adapts itself according to Pareto dominance and serves as the
setpoint for the low level multiplier controller to schedule loss landscapes
via output feedback of each loss term. Our method is multiplier-free and
operates at the timescale of epochs, thus saves tremendous computational
resources compared to full training cycle multiplier tuning. It also
circumvents the excessive memory requirements and heavy computational burden
of existing multi-objective deep learning methods. We applied it to domain
invariant variational auto-encoding with 6 loss terms on the PACS domain
generalization task, and observed robust performance across a range of
controller hyperparameters, as well as different multiplier initial
conditions, outperforming other multiplier scheduling methods. We offered
modular implementation of our method, admitting extension to custom definition
of many loss terms.

##  1 Introduction

In many deep learning fields including domain generalization (Gulrajani and
Lopez-Paz, 2020), the loss function for neural networks amounts to a summation
of multiple terms, with a multiplier weighting each term. As the number of
loss terms increases, the space of choice for those multipliers grows
dramatically. In addition, when the validation set is drawn from a different
distribution than the target domain, the validation performance does not
faithfully reflect the target domain performance (Chen et al., 2022b) and can
easily reach and maintain a saturated value during training. All these factors
complicate and add difficulties to the process of model selection using
conventional hyperparameter tuning (Probst et al., 2019) methods.

Instead of optimizing a fixed weighted summation of many neural network loss
terms using a single objective optimization algorithm which relies on heavy
hyperparameter tuning of the weight multipliers, multi-objective optimization
on these loss terms has raised attention in the deep learning community
recently (Sener and Koltun, 2018; Lin et al., 2019; Mahapatra and Rajan, 2021;
Chen et al., 2022b). Previous works in this direction (Mahapatra and Rajan,
2021; Chen et al., 2022b) addresses the problem via careful choice of
multipliers at each optimization iteration instead of fixed for the whole
training cycle, but suffers from excessive memory requirement and heavy extra
computation burden, which hinders their applicability to large dataset and
large scale neural networks Chen et al. (2022b). Furthermore, it is restricted
to Exact Pareto optimization (Mahapatra and Rajan, 2021) where the solution
should be aligned to a user input preference vector, while how to set that is
not clear in the general case. Their performance further relies on careful
tuning of learning rate of single objective low level optimization algorithms
Chen et al. (2022b). Single objective optimization algorithms like ADAM
(Kingma and Ba, 2015) were analyzed predominantly with convex optimization
theory (Kingma and Ba, 2015; Reddi et al., 2019), which does not hold globally
(Milne, 2019). From a control theory point of view, these low level
optimization algorithms drive the model parameters iteratively, which forms a
dynamic system. Due to lack of comprehensive theoretical underpinnings of the
loss function property with respect to model parameters, this dynamic system
has unreliable open-loop behavior.

Inspired by feedback control theory (Doyle et al., 2013) and feedback
optimization in power systems (Picallo et al., 2020; He et al., 2023;
Hauswirth et al., 2024), we treat the dynamic system of model parameters
driven by the low level optimization algorithm as an uncontrolled plant, with
model parameter as state and the many loss terms as output. We use state
dependent multiplier induced loss landscapes as control input and designed a
hierarchical control structure (Dietterich et al., 1998; Kulkarni et al.,
2016; Sun et al., 2020) to drive the system towards a multi-objective descent
goal, without modifying the internal dynamic of low level neural network
optimization algorithms.

This hierarchical control structure allows for a more refined and flexible
approach to high-level optimization of neural networks, by automatically
balancing the trade-off between a task-specific empirical loss and various
regularization losses, which we detail in Equation 1. We depicted our idea in
the control diagram in Figure 1 where the motivation of choosing multiplier
induced landscape scheduling as control input is illustrated in Figure 4.

We demonstrated that our method removes the combination curse of
dimensionality of multi-dimensional multipliers, via training domain invariant
variational auto-encoding (Ilse et al., 2020) on the PACS dataset. Our method
only has a small number of hyperparameters for our controllers while achieving
rubust performance when changing them. The baselines include different warmup
multiplier scheduling strategies, which still need explicit specification of
multiplier ultimate values. Their performance varies drastically with changes
in those ultimate values. We attribute the improved performance of our method
to the automatic trade-off of different loss terms. Our major contributions
are:

  * •

We proposed a probabilistic graphical model (PGM) for depicting the joint
multiplier and model parameter adaptation process with a likelihood function
imposed, promoting multi-objective descent of each loss term, which can be
used to described general multi-objective deep learning. We proposed to
transform the problem of parameter inference of this PGM into a sequential
decision process through an optimal control formulation.

  * •

We proposed to solve the sequential decision process for the multi-objective
descent goal via breaking the goal into a series of constraint optimization
sub-goals which forms a hierarchical controller that automatically adjusts the
value of both the multiplier and the model parameters during the training
process. The constraint bound updates itself via Pareto dominance. The multi-
objective descent of the empirical risk and multidimensional penalty
(regularization) loss in Equation 1 offers additional model selection criteria
along with validation set performance. Compared to Bayesian optimization, our
controller operates at the timescale of model parameter dynamic epoch instead
of a complete training, which saves tremendous computational resources.
Compared to existing multi-objective deep learning methods, it circumvent
excessive memory requirements and heavy extra computational burden.

  * •

When treating the neural network model parameter dynamic system driven by low
level optimization algorithm as uncontrolled plant, we proposed to use state
dependent multiplier induced loss landscapes as control input (see Figure 4)
and offered initial theoretical analysis of the multi-objective descent
behavior of our closed loop system . We offered modular implementation of our
method 111https://github.com/marrlab/DomainLab/tree/mhof and showed robustness
of our method with respect to different multiplier initial conditions and
controller hyperparameters for domain generalization on the PACS dataset with
domain invariant auto-encoding.

##  2 Preliminaries

###  2.1 Structural Risk Minimization

Structural Risk Minimization (SRM) is a principle in learning algorithms that
balances training performance and model complexity via optimizing the
penalized loss in Equation 1.

$\displaystyle\begin{split}L(\theta,\mu,\mathcal{D}_{tr})&=\ell(\theta,\mathcal{D}_{tr})+\mu^{T}R(\theta,\mathcal{D}_{tr})\\\
\mu,R(\cdot)&\in\mathbb{R}_{+}^{d}\end{split}$ (1)

In Equation 1, we consider loss $L$ with training data $\mathcal{D}_{tr}$ and
model parameters $\theta$. $\ell(\theta,\mathcal{D}_{tr})$ represents the
empirical risk and $R(\theta,\mathcal{D}_{tr})$ represents the regularization
term of the loss function, also referred to as the penalty function and thus
$L$ the penalized loss. $\mathbb{R}_{+}^{d}$ is the $d$ dimensional positive
octant. We refer to the multiplier $\mu$ as the penalty multiplier, as a
special case of general hyperparameters (e.g. learning rate) in deep learning.

##### Notation

In the following sections, for brevity, we omit the implicit dependence of
$\ell(\cdot)$ and $R(\cdot)$ on $\mathcal{D}_{tr}$ and $\theta$. We also write
$R(\theta,\cdot)$ where we need an explicit dependence on $\theta$ but use
$\cdot$ only to represent data $\mathcal{D}_{tr}$ and/or other arguments. We
use subscript of $R(\cdot)$ to indicate the component of $R(\cdot)$
corresponding to that subscript, e.g., if $\mu=[\beta,\gamma]$, then
$R_{\beta}(\cdot)$ corresponds to the component of $R(\cdot)$ associated with
$\beta$. We use $X,Y\sim\mathcal{D}$ to represent a random variable $X$ and
its supervision signal $Y$ from distribution $\mathcal{D}$ and we use the same
symbol $\mathcal{D}$ to denote a corresponding dataset as well. Let $x$ be an
instance of $X$, so is $y$ for $Y$.

We use superscript $k$ in bracket to index the optimization iteration, which
by default corresponds to one epoch for our method (it can also correspond to
one training cycle in Equation 2 to iteratively achieve Equation 3 for
conventional hyperparameter optimization, see Remark 3.1).

######  Definition 2.1.

When optimizing Equation 1 with multiplier $\mu$ iteratively (e.g. with
gradient based methods like ADAM), the model parameter evolution dynamic
$\theta^{(k+1)}=f_{\theta}(\mu,\theta^{(k)},\mathcal{D}_{tr})$ brings
$\theta^{(k)}$ to its next value $\theta^{(k+1)}$. $\theta^{(k+1)}$ then
becomes the initial value for further iterations. For notational simplicity,
we replace the explicit dependence of this dynamic on training data
$\mathcal{D}_{tr}$, multiplier $\mu$ (potentially varying value for each
iteration) and/or even the initial model parameter $\theta$ with dots. This
convention results in notations $f_{\theta}(\theta,\cdot)$,
$f_{\theta}(\mu,\cdot)$ and $f_{\theta}(\cdot)$. We also write
$\theta^{+}=f_{\theta}(\mu,\theta,\cdot)$ and
$\theta^{+}=f_{\theta}(\theta,\cdot)$ (potentially varying $\mu$ for each
iteration) to represent the model parameter dynamic system without explicitly
using $k$.

###  2.2 Multi-objective optimization

######  Definition 2.2.

We use $R(\theta_{1},\cdot)\prec R(\theta_{2},\cdot)$ (see Equation 1) to
indicate each component of the $d$ dimensional vector $R(\theta_{1},\cdot)$ is
$\leq$ the corresponding component of $R(\theta_{2},\cdot)$. We use $\nprec$
as the negation of $\prec$. The clause
$\ell(\theta_{1},\cdot),R(\theta_{1},\cdot)\nprec\ell(\theta_{2},\cdot),R(\theta_{2},\cdot)$
AND
$\ell(\theta_{2},\cdot),R(\theta_{2},\cdot)\nprec\ell(\theta_{1},\cdot),R(\theta_{1},\cdot)$
defines an equivalence relation which we denote as
$\theta_{1}\sim\equiv\theta_{2}$.

######  Definition 2.3.

Starting from $\theta^{(0)}$, under dynamic system
$\theta^{+}=f_{\theta}(\theta,\cdot)$, we use
$s\mathcal{R}\left(\theta^{(0)},f_{\theta}(\cdot)\right)$ to represent the set
of model parameters ($\theta$ points) that can be reached. Accordingly, we use
$v\mathcal{R}_{\ell(\cdot),R(\cdot)}\left(\theta^{(0)},f_{\theta}(\cdot)\right)$
to represent the corresponding set of multi-objective function values.

######  Definition 2.4.

If $\theta_{1}\sim\equiv\theta_{2}$ where $\theta_{1}\in
s\mathcal{R}(\theta^{(0)},f_{\theta}(\cdot))$ and $\theta_{2}\in
s\mathcal{R}(\theta^{(0)},f_{\theta}(\cdot))$ as in Definition 2.3, we define
$\theta_{2}\in\mathcal{C}_{\ell(\cdot),R(\cdot)}\left(\theta_{1},f_{\theta}(\cdot),\theta^{(0)}\right)$
the map from a representing element $\theta_{1}$ to its equivalence class,
which contains $\theta_{2}$. We use
$\mathcal{E}_{\ell(\cdot),R(\cdot)}\left(\theta_{1},f_{\theta}(\cdot),\theta^{(0)}\right)$
to represent the set of $\ell(\cdot),R(\cdot)$ values for each element of
$\mathcal{C}_{\ell(\cdot),R(\cdot)}\left(\theta_{1},f_{\theta}(\cdot),\theta^{(0)}\right)$
.

###  2.3  Hyperparameter tuning

The multiplier $\mu$ in Equation 1 is a hyperparameter that affects the
evolution of the model parameters $\theta$ in the optimization process for $L$
defined in Equation 1, which leads to optimized model parameter
$\theta_{\mu,\mathcal{D}_{tr}}$ below, where we omit the potential dependence
of $\theta_{\mu,\mathcal{D}_{tr}}$ on its initial condition $\theta^{(0)}$ for
notational simplicity.

$\displaystyle\theta_{\mu,\mathcal{D}_{tr}}$
$\displaystyle=\arg\min_{\theta}L(\theta,\mu,\mathcal{D}_{tr}).$ (2)

Hyperparameter optimization aims to adjust the multiplier $\mu$ in alignment
with specific performance metrics ${O}$ (e.g. validation set accuracy for a
classification task) in Equation 3 evaluated on the validation set
$\mathcal{D}_{val}$, which lead to the selected hyperparameter
$\mu_{\mathcal{D}_{tr},\mathcal{D}_{val}}$.

$\displaystyle\mu_{\mathcal{D}_{tr},\mathcal{D}_{val}}=\arg\min_{\mu}O(\theta_{\mu,\mathcal{D}_{tr}},\mathcal{D}_{val}).$
(3)

The process of hyperparameter optimization in Equation 3 leads to selected
model $\theta_{\mu_{\mathcal{D}_{tr},\mathcal{D}_{val}},\mathcal{D}_{tr}}$,
which can be conceptualized as a problem of algorithm configuration (Hutter et
al., 2007, 2014; López-Ibáñez et al., 2016; Sun et al., 2020) where the
hyperparameter configures a machine learning algorithm.

To solve Equation 3 iteratively, each iteration relies on a complete training
cycle to optimize Equation 2. The selection of $\mu$ depends on validation set
$\mathcal{D}_{val}$, which can only be drawn from a different distribution
compared to the target domain for the domain generalization problem (Section
2.4). Additionally, the performance metric $O$ on the in-domain validation set
can reach and maintain a saturated value during training with neural networks
of high expressive power.

###  2.4 Domain generalization

Domain generalization aims at enabling the trained neural network to achieve
robust generalization to unseen domains exhibiting distribution shifts
(Gulrajani and Lopez-Paz, 2020; Sun et al., 2019b). Many domain generalization
methods (Ganin et al., 2016; Li et al., 2018; Levi et al., 2021; Carlucci et
al., 2019; Ilse et al., 2020; Sun and Buettner, 2021; Rame et al., 2022)
promotes this goal via domain invariant representation by adding domain
invariant regularization losses $R(\cdot)$ upon task-specific losses
$\ell(\cdot)$ (e.g. classification loss) on the training data in Equation 1.

As an example of $\ell(\cdot)+\mu^{T}R(\cdot)$ loss for domain generalization,
we label each component of the loss of Domain Invariant Variational
Autoencoding (DIVA) (Ilse et al., 2020) in the following equations, where we
use $x$ to denote the input instance (e.g. input image), $y$ to denote the
corresponding supervision class label, and $d$ the domain label. DIVA (Ilse et
al., 2020) aims to disentangle domain specific information with learned
features $z_{d}$ and class specific information with $z_{y}$, and has two
corresponding classification losses $\mathbb{E}_{q_{\phi_{y}}(z_{y}|x)}[\log
q_{\omega_{y}}(y|z_{y})]$ (classifying the correct class label based on
$z_{y}$) and $\mathbb{E}_{q_{\phi_{d}}(z_{d}|x)}[\log
q_{\omega_{d}}(d|z_{d})]$ (classifying the correct domain label based on
$z_{d}$).

$\displaystyle L(\theta,\mu,x,y,d)=\ell(\cdot)+\mu^{T}R(\cdot)=$ (4)
$\displaystyle\gamma_{y}\mathbb{E}_{q_{\phi_{y}}(z_{y}|x)}[\log
q_{\omega_{y}}(y|z_{y})]+$ (5)
$\displaystyle\mu_{recon}\mathbb{E}_{q_{\phi_{d}}(z_{d}|x),q(z_{x}|x),q(z_{y}|x)}\log
p_{\theta_{r}}(x|z_{d},z_{x},z_{y})+$ (6)
$\displaystyle[-\beta_{x}KL(q_{\phi_{x}}(z_{x}|x)||p_{\theta_{x}}(z_{x}))]+$
(7)
$\displaystyle[-\beta_{y}KL(q_{\phi_{y}}(z_{y}|x)||p_{\theta_{y}}(z_{y}|y))]+$
(8)
$\displaystyle[-\beta_{d}KL(q_{\phi_{d}}(z_{d}|x)||p_{\theta_{d}}(z_{d}|d))]+$
(9) $\displaystyle\gamma_{d}\mathbb{E}_{q_{\phi_{d}}(z_{d}|x)}[\log
q_{\omega_{d}}(d|z_{d})]$ (10)

Here, $KL$ denotes Kullback-Leibler divergence, $\mathbb{E}$ denotes
expectation, $p(\cdot)$ stands for the prior distribution or the distribution
of the generative model, and $q(\cdot)$ stands for the approximate posterior
distribution. For more details, refer to (Ilse et al., 2020). In the above
situation, we have
$\mu=[\mu_{recon},\beta_{x},\beta_{y},\beta_{d},\gamma_{d}]$ and
$\theta=[\phi_{y},\phi_{d},\theta_{r},\theta_{d},\theta_{x},\theta_{y},\phi_{x},\phi_{y},\phi_{d},\omega_{d},\omega_{y}]$
with each entry representing the weight (with bias) of a neural network.

In (Ilse et al., 2020), $\mu_{recon}=1.0$ in Equation 6. Additionally,
$\gamma_{y}$ and $\gamma_{d}$ are maintained as constants. $\gamma_{y}$ is
corresponding to
$\ell(\cdot)=\gamma_{y}\mathbb{E}_{q_{\phi_{y}}(z_{y}|x)}[\log
q_{\omega_{y}}(y|z_{y})])$ in Equation 5, and $\gamma_{d}$ in Equation 10 is
associated with one component of $R(\cdot)$. However, the combined choice of
$\gamma_{y}$ and $\gamma_{d}$ significantly influences the generalization
performance, as shown by our experimental findings in Section 5.

For the multiplier $\beta_{x},\beta_{y},\beta_{d}$ corresponding to other
components of $R(\cdot)$, Ilse et al. (2020) used a warm-up strategy to
increase them gradually from a small value to a pre-defined value. This,
however, still requires a choice of the ultimate values for
$\beta_{x},\beta_{y},\beta_{d}$. In (Ilse et al., 2020), these ultimate values
are simply set to one. We coin this kind of multiplier scheduling strategy a
_feedforward_ scheme. Note that this strategy does not reduce the number of
hyperparameters to be selected, since ultimate values for the multipliers
still have to be specified.

##  3 Methods

In this section, we elaborate on the hierarchical output feedback control
strategy for the multidimensional multiplier adaptation in Equation 1 (see
Figure 1) towards multi-objective optimization. We propose a Probabilistic
Graphical Model (PGM) for the joint model parameter and multiplier adaptation
process in Section 3.1, which leads to an optimal control formulation of
parameter estimation, extending the classical hyperparameter tuning-based
approach to parameter-multiplier co-evolution in Section 3.2.

Figure 1: A control diagram illustrating the hierarchical output feedback
optimization process we proposed with multi-objective setpoint adaptation. The
uncontrolled plant corresponds to an open loop dynamic
$\theta^{+}=f_{\theta}(\mu,\theta,\cdot)$ via optimizing Equation 1 with a low
level optimization algorithm (e.g. ADAM) lacking feedback. The $\mu$
controller adjusts the multiplier $\mu$, thus schedules loss landscape (see
Illustration in Figure 4), based on the difference $e$ between the setpoint
$b$ (defined in Equation 23) and the measured output components $R(\cdot)$,
guiding the optimization of model parameters $\theta$ through a feedback loop.
The setpoint $b$ is adjusted via feedback from $\ell(\cdot)$ and $R(\cdot)$
(Equation 24a) which forms a higher hierarchy.

###  3.1  Probabilistic graphical model for parameter and multiplier adaptive
optimization process

Figure 2: A probabilistic graphical model for the sequential decision process
of joint multiplier and model parameter adaptation in multiple-domain (data-
site specific distribution) supervised learning.

Gradient based optimization algorithms for neural network defines a dynamic
system for model parameters where the weight multiplier in Equation 1
configures such a dynamic system. Instead of using fixed multipliers (or fixed
multipliers’ ultimate values), we propose to choose the multiplier values also
dynamically. Due to the stochastic nature of the gradient based optimization
process in deep learning with different subsample and mini-batch distributions
for each iteration, it is appropriate to use a Probabilistic Graphical Model
(PGM) to describe such a stochastic process. PGM has been used to describe
sequential decision processes like reinforcement learning (Levine, 2018; Sun
and Bischl, 2019). We propose Figure 2 to describe the sequential decision
process of model parameter and multiplier adaptation in the form of a
probabilistic graphical model (PGM) (Koller and Friedman, 2009). In the upper
part of Figure 2, we model the generative process of $N$ observations of the
instance $X$ (e.g. an image) and the corresponding supervision signal $Y$
(e.g. a class label). The path $\beta\rightarrow D_{ob}\rightarrow X$, where,
$D_{ob}$ is the data-site Sun et al. (2019a) specific prior parameter (domain
variable, characterizing distribution) with hyper-prior $\beta$, models the
information which is not captured by the supervision signal $Y$ of $X$.

In the lower part of the graphical model, with superscript $k$ representing
the iteration for the decision process, we use $\mu^{(k)}$ to represent the
multiplier of the loss in Equation 1 at the $k$th iteration. The multiplier
$\mu^{(k)}$ serves as a configuration parameter for the optimization process
of Equation 1 that affects the value of $\theta_{\mu}^{(k)}$, evolved from its
previous value $\theta^{(k-1)}$. $\theta_{\mu}^{(k)}$ then becomes the initial
value for the next iteration. This can be described in Equation 14 with
explicit dependence of $\theta^{(k)}$ on $\mu^{(k)}$ omitted. The model
parameter $\theta^{(k)}$ and observed data $X$, supervision signal $Y$ co-
parent the performance indicator as output ${O}^{(k)}$.

The plate replicates $B$ in the lower part represents the number of
optimization iterations (each iteration correspond to $N$ observations). The
adaptive generation of the next multiplier $\mu^{(k+1)}$ depends on the
previous value $\mu^{(k)}$ (the self-loop dashed arrow in Figure 2) and the
feedback information of the output $O^{(k)}$ (long dashed arrow), as in
Equation 11.

$\mu^{(k+1)}=f_{\mu}(\mu^{(k)},O^{(k)},\cdot)$ (11)

######  Remark 3.1.

An example of $f_{\mu}$ can be based on the acquisition function in Bayesian
optimization (BO) (Garnett, 2023) for approximating Equation 3, which leads to
a BO realization of this PGM in (Sun et al., 2019a), where $O$ in Equation 11
is usually chosen to be the validation set prediction performance and
iteration index $k$ corresponds to a whole training cycle of Equation 2 to
iteratively achieve Equation 3.

Different from the BO realization of this PGM, which generates multipliers at
a much slower timescale (usually once at the end of training) than that for
the evolution of model parameters in neural network training, the hierarchical
output feedback optimization we proposed (Figure 1 and elaborated in following
sections) operates at the timescale of epoch, thus forms another realization
of this PGM with iteration index $k$ corresponding to epochs. In addition, we
have $O(\cdot)=\ell(\cdot),R(\cdot)$ for the PGM in Figure 2 in contrast to
validation set performance used in BO.

We introduce an objective function for estimating $\theta_{\mu^{(B)}}$ based
on the profile likelihood in eq. 12 in accordance with the PGM in Figure 2.
This likelihood promotes a high $e\mathcal{HV}(\cdot)$ in Definition 3.1 with
respect to $\theta_{\mu^{(B)}}$. $\Xi$ is the normalization factor.

$\displaystyle\mathcal{P}(\theta_{\mu^{(B)}})$
$\displaystyle=\frac{1}{\Xi}\exp\left(e\mathcal{HV}_{O(\cdot)}(\theta_{\mu^{(B)}},\theta^{(0)},f_{\theta},\cdot)\right)$
(12)

######  Definition 3.1.

$e\mathcal{HV}_{O(\cdot)=\ell(\cdot),R(\cdot)}(\theta,\theta^{(0)},f_{\theta},\cdot)$
maps $\theta$ to the dominated hypervolume (Zitzler and Künzli, 2004; Zitzler
et al., 2007; Guerreiro et al., 2021) of
$\mathcal{E}_{\ell(\cdot),R(\cdot)}\left(\theta,f_{\theta}(\cdot),\theta^{(0)}\right)$
in Definition 2.4 with respect to reference point
$\ell(\theta^{(0)},\cdot)\in\mathbb{R},R(\theta^{(0)},\cdot)\in\mathbb{R}^{d}$.
As illustrated in Figure 3.

The PGM we defined here can also be used to describe other multi-objective
deep learning methods Mahapatra and Rajan (2021) via online combinatorial
choice of multipliers.

######  Remark 3.2.

The PGM also provides a potential possibility to impose a Bayesian
interpretation to the joint multiplier tuning and model training process for
frequentist machine learning models, the deeper investigation of which we
leave for future work.

Figure 3: Illustration to $e\mathcal{HV}$ in Definition 3.1: We take
$[\ell(\theta^{(0)},\cdot)\in\mathbb{R},R(\theta^{(0)})\in\mathbb{R}^{d=1}]$
as reference point. Suppose $\\{A,B,C,D,E,F\\}$ in the illustration constitute
$v\mathcal{R}\left(\theta^{(0)},f_{\theta}(\cdot)\right)$ in Definition 2.3.
For any $\theta$ corresponding to the point A in the illustration with
coordinate $[\ell(\theta,\cdot),R(\theta,\cdot)]$, $e\mathcal{HV}$ maps
$\theta$ to
$\mathcal{C}_{\ell(\cdot),R(\cdot)}\left(\theta,f_{\theta}(\cdot),\theta^{(0)}\right)$,
then to
$\mathcal{E}_{\ell(\cdot),R(\cdot)}\left(\theta,f_{\theta}(\cdot),\theta^{(0)}\right)$
(see Definition 2.4), which is the point set {A,B,C,D} , then it calculates
the dominated hypervolume with respect to the reference point (union of the
shaded rectangles).

###  3.2 Hierarchical control formulation of joint parameter multiplier
optimization

####  3.2.1 Optimal control formulation

Analogous to how optimal control Lewis et al. (2012) can be formulated as
maximum likelihood estimation in a probabilistic graphical model (Levine,
2018; Sun and Bischl, 2019), the estimation of $\theta^{(B)}$ based on the
profile likelihood like objective in Equation 12 for Figure 2 is also a
sequential decision process. From a statistical inference point of view, we
are conducting a sequential inference (i.e., finding a sequence of $\mu$ in
Equation 13), such that the last step $B$ has the best possible hyper-volume
with respect to output $O$ (in a stochastic sense, see Figure 3), where $O$
can be a performance metric like accuracy on the validation set in Bayesian
optimization, or defined in Equation 15 for our method.

To reformulate the estimation problem as an optimal control problem, we
present the sequential decision on choosing a sequence of multipliers
$\mu^{(1)},\cdots,\mu^{(B)}$ defined in Equation 1 with respect to the
objective function in Equation 13 at iteration $B$ which promotes Equation 12.

$\displaystyle\min_{\mu^{(1)},\cdots,\mu^{(k)},\cdots,\mu^{(B)}}$
$\displaystyle
e\mathcal{HV}_{O(\cdot)}(\theta_{\mu^{(B)}},\theta^{(0)},f_{\theta},\cdot)$
(13) $\displaystyle\mathrm{s.t.}\leavevmode\nobreak\ $
$\displaystyle\theta^{(k)}=f_{\theta}(\mu^{(k)},\theta^{(k-1)},\cdot)$ (14)
$\displaystyle
O^{(k)}=\left[\ell(\theta^{(k)}|\cdot),R(\theta^{(k)}|\cdot)\right]$ (15)
$\displaystyle k=0,1,\cdots,B$ (16)

Here, we use $\theta$ to denote the neural network weights, which can be
regarded as the state of a dynamic system, where the state transition is
governed by a low level model parameter optimization dynamic $f_{\theta}$ in
Equation 14 (e.g. ADAM (Kingma and Ba, 2015)). From a control theory point of
view, $f_{\theta}$ depicts the uncertain dynamic of a plant to be controlled,
while the task here is to define a controller to generate a $\mu$ sequence to
guide the uncertain dynamic of $\theta$ with respect to plant output $O$, as
depicted in Figure 1.We treat the loss terms $R(\cdot)$ and $\ell(\cdot)$ as
the output of the system in Equation 15. From a statistical point of view,
$f_{\theta}$ corresponds to the parent structure of $\theta^{(k)}$ in Figure
2, which drives $\theta$ to the next value.

Although we have no complete knowledge of $f_{\theta}$ from Equation 14
describing low level optimization algorithms like ADAM, in practice, the
penalized loss from Equation 1 often descent after some number of iterations
(possibly with oscillations in between), which leads to Definition 3.2.

######  Definition 3.2.

A penalized loss descent operator $\mathcal{G}$ satisfies that if
$\theta^{+}\in\mathcal{G}_{\mu}(\theta;\ell(\cdot),R(\cdot))$, then
$\ell(\theta^{+})+\mu R(\theta^{+})<\ell(\theta)+\mu R(\theta)$

Therefore, we take the frequent occurrences of operators in Definition 3.2
during neural network training as a mild assumption in our following
discussion. For instance, we could run $f_{\theta}$ in Equation 14 in one step
to have
$\ell(\theta^{(k+1)})+\mu^{(k+1)}R(\theta^{(k+1)})<\ell(\theta^{(k)})+\mu
R(\theta^{(k)})$, or in more than one steps to achieve Definition 3.2.

Our final goal, however, is to ensure a joint descent of $\ell(\cdot)$ and
$R(\cdot)$, which leads to Definition 3.3:

######  Definition 3.3.

A Pareto-descent operator takes $\theta$ to $\theta^{(+)}$ (via one or several
iterations) which descents $R(\cdot)$ and $\ell(\cdot)$ simultaneously, i.e.,

$\displaystyle\ell(\theta^{+})$ $\displaystyle<\ell(\theta)$ (17)
$\displaystyle R(\theta^{+})$ $\displaystyle\preceq R(\theta)$ (18)

####  3.2.2 Multi-objective optimization via shrinking the reference signal in
constrained optimization

A concrete form of Equation 11 will detail how $\mu$ is updated by a
controller $f_{\mu}$. Combining Equations 11 and 14, we have

$\displaystyle\theta^{(k+1)}$
$\displaystyle=f_{\theta}(f_{\mu}(\mu^{(k)},O^{(k)}),\theta^{(k)},\cdot)$ (19)
$\displaystyle=f_{\theta,\mu}^{k}(\mu^{(0)},\theta^{(0)},\cdot)$ (20)

where we use $f_{\theta,\mu}^{k}$ to represent the $k$-times compound
recursive evaluations of $f_{\theta}$ and $f_{\mu}$ in Equation 19 until the
dependence is only on $\theta^{(0)},\mu^{(0)}$. We call Equation 20 the
_closed-loop_ dynamics for $\theta$.

The question then arises: how can we design such a closed-loop dynamic system
to ensure a Pareto descent of $O(\cdot)$ (both $R(\cdot)$ and $\ell(\cdot)$),
leading to multi-objective optimization in Definition 3.3 or Equation 13?

To resolve this, we reduce it to a simpler problem, i.e., when bounding the
regularization term $R(\cdot)$ from above by a time-varying reference
$b^{(k)}$ (a.k.a. setpoint) in eq. 26, what value can we achieve for
$\ell(\cdot)$ at Equation 25. If we could ensure the bound changes
monotonically, multi-objective optimization can be achieved.

First, we define the time-varying reference signal being governed by the
following:

$\displaystyle b^{(0)}$ $\displaystyle=\rho R^{(0)}$ (21) $\displaystyle 0<$
$\displaystyle\rho\in\mathbb{R}<1$ (22) $\displaystyle b^{(k)}$
$\displaystyle=g_{b}(\ell^{(0:k)}(\cdot),R^{(0:k)}(\cdot))$ (23)

where $g_{b}({\cdot})$ represents the mapping from the previous value at the
training step $k-1$ to the new value at step $k$, defined as

$\displaystyle b^{(k)}=$ $\displaystyle R^{(k)},\leavevmode\nobreak\ \text{if
}R^{(k)}\prec b^{(k-1)}\land\ell^{(k)}<\min_{j}\ell^{(j)}$ (24a)
$\displaystyle b^{(k)}=$ $\displaystyle j=0,\cdots,k-1$ (24b) $\displaystyle
b^{(k)}=$ $\displaystyle b^{(k-1)},\leavevmode\nobreak\ \text{otherwise}$
(24c)

######  Remark 3.3.

Note that although $b^{(k)}$ is only defined for $R^{(k)}(\cdot)$, the
shrinkage of $b^{(k)}$ also depends on $\ell(\cdot)$ decrease, thus indicates
multi-objective Pareto-descent in Definition 3.3.

With $b^{(k)}$ defined above, under $s_{k}$ iterations (number of iterations
$s_{k}$ depends on $\theta^{(k)}$), with $b^{(k)}$ fixed, i.e.
$b^{(k)}=b^{(k+1)}=\cdots=b^{(k+s_{k}-1)}$, suppose the following constrained
optimization Bertsekas (2014) problem starting with $\theta^{(k)}$ has a
solution and could be approximated:

$\displaystyle\ell^{(k+s_{k})}(\cdot)=\min_{\mu^{(k+1)},\cdots,\mu^{(k+s_{k})}}\leavevmode\nobreak\
\ell(\cdot)$ (25) $\displaystyle\quad\mathrm{s.t.}\leavevmode\nobreak\ $
$\displaystyle R^{(k+s_{k})}(\cdot)\prec b^{(k)}$ (26) $\displaystyle
R^{(j)}(\cdot)\nprec b^{(k)}$ (27)
$\displaystyle\theta^{(j+1)}=f_{\theta}(\mu^{(j+1)};\theta^{(j)})\leavevmode\nobreak\
(\lx@cref{creftypecap~refnum}{eq:controller_theta_abs})$ (28) $\displaystyle
j=k,k+1,\cdots,k+s_{k}-1$ (29)

Then, regardless if the minimization in Equation 25 is attainable or not, as
long as we have

$\displaystyle\ell^{(k+m)}(\cdot)$ $\displaystyle<\ell^{(k)}(\cdot),$ (30)
$\displaystyle\exists\leavevmode\nobreak\ m$ $\displaystyle\geq s_{k}$ (31)

, due to Equation 26, Equation 24a, we still achieve multi-objective decrease
in Definition 3.3 at step $k+m$ compared to step $k$.

The evolution process of the model parameters in this case can be regarded as
a plant with uncertain dynamic $f_{\theta}$ in Equations 14 and 28, where the
only control we have at hand is the sequence
$\mu^{(k+1)},\cdots,\mu^{(k+s_{k})}$, which provides different loss landscapes
at different iterations with different $\theta$ values as depicted in Figure
4. How can we design such $\mu$ or loss landscape sequence to promote Equation
26 to happen?

To simplify the discussion, first consider
$\mu,R(\cdot)\in\mathbb{R}_{+}^{d=1}$. Before the constraint in
LABEL:{eq:feasibility_r} is satisfied, design
$\mu^{(k+1)},\cdots,\mu^{(k+s_{k})}$ to be an increasing sequence Bazaraa et
al. (2013), which gives more weight to the corresponding regularization term
in $R(\cdot)$, to approximate the constraint optimization in Equation 25 by
optimizing Equation 1 for each $\mu$ until Equation 26 is satisfied. However,
will the increased weight on the regularization term in $R(\cdot)$ result in
deteriorated $\ell(\cdot)$?

Figure 4:  Illustration of the concept of multiplier induced loss landscape
scheduling we proposed, used as control signal in Figure 1: The lifted
landscape $\ell(\cdot)+\mu^{(k+1)}R(\cdot)$ (dotted curve) scheduled at
iteration $k+1$, enables the model parameter dynamic to overcome the local
minimum $\bar{\theta^{(k)}_{\mu^{(k)}}}$ of the old loss landscape
$\ell(\cdot)+\mu^{(k)}R(\cdot)$ (solid curve) scheduled at iteration $k$, as
can be imagined via the two balls with different colors rolling along the
corresponding loss landscapes. It showcases a scenario corresponding to
Definition 3.4: In comparison to $\theta^{(k)}$, $\theta^{(k+1)}$ corresponds
to a decreased $\ell(\cdot)+\mu^{(k+1)}R(\cdot)$ value but increased
$\ell(\cdot)+\mu^{(k)}R(\cdot)$ value.

One possibility can be, although $f_{\theta}(\mu^{(k+1)},\cdot)$ brings
$\theta^{(k)}$ to the next value $\theta^{(k+1)}$ which decreases the
penalized loss in Equation 1 with the new $\mu^{(k+1)}$, but the new
$\theta^{(k+1)}$ evaluates to a deteriorated penalized loss with respect to
the old $\mu^{(k)}$ as shown in Figure 4. To describe this behavior during the
increasing process of $\mu$, we have Equation 33 in Definition 3.4. This
situation is of particular interest due to the possibility of overcoming local
minima of a fixed loss landscape, as indicated in Figure 4.

######  Definition 3.4.

A reg-Pareto slider $\mathcal{A}(\cdot)$ with respect to function
$\ell(\cdot)$ and $R(\cdot)$ is defined to be a set map, s.t. if
$\\{\theta^{(k+1)},\mu^{(k+1)}\\}\in\mathcal{A}(\mathcal{\ell}(\cdot),R(\cdot);\theta^{(k)},\mu^{(k)})$
then

$\displaystyle\ell(\theta^{(k+1)})+\mu^{(k+1)}R(\theta^{(k+1)})$
$\displaystyle\leq\ell({\theta}^{(k)})+\mu^{(k+1)}R({\theta}^{(k)})$ (32)
$\displaystyle\ell(\theta^{(k+1)})+\mu^{(k)}R(\theta^{(k+1)})$
$\displaystyle\geq\ell({\theta}^{(k)})+\mu^{(k)}R({\theta}^{(k)})$ (33)
$\displaystyle\mu^{(k+1)}$ $\displaystyle\geq\mu^{(k)}$ (34)

In this situation, take $R(\cdot)\in\mathbb{R}^{d=1}$, we have the following
conclusion in Corollary 3.1 where we show the increase of $\ell(\cdot)$ as an
expense to decrease of $R(\cdot)$ is bounded.

######  Corollary 3.1.

Let

$\displaystyle\\{\theta^{(k+1)},\mu^{(k+1)}\\}\in\mathcal{A}(\mathcal{\ell}(\cdot),R(\cdot);\theta^{(k)},\mu^{(k)})$
(35)

from definition 3.4 then,

$\displaystyle R(\theta^{(k+1)})$ $\displaystyle\leq R({\theta}^{(k)})$ (36)
$\displaystyle\ell(\theta^{(k+1)})$ $\displaystyle\geq\ell({\theta}^{(k)})$
(37) $\displaystyle\ell(\theta^{(k+1)})-\ell({\theta}^{(k)})$
$\displaystyle\leq\mu^{(k+1)}(R({\theta}^{(k)})-R(\theta^{(k+1)}))$ (38)

###### Proof.

Multiply $-1$ to both sides of Equation 33, add the results to Equation 32, we
have

$(\mu^{(k+1)}-\mu^{(k)})R(\theta^{(k+1)})\leq(\mu^{(k+1)}-\mu^{(k)})R({\theta}^{(k)})$
(39)

which implies

$(\mu^{(k+1)}-\mu^{(k)})\left(R(\theta^{(k+1)})-R({\theta}^{(k)})\right)\leq
0$ (40)

Since $\mu^{(k+1)}>\mu^{(k)}$

$R(\theta^{(k+1)})-R({\theta}^{(k)})\leq 0$ (41)

Equation 33 and 32 give

$\displaystyle 0\leq\mu^{(k)}(R({\theta}^{(k)})-R(\theta^{(k+1)}))\leq$ (42)
$\displaystyle\ell(\theta^{(k+1)})-\ell({\theta}^{(k)})\leq\mu^{(k+1)}(R({\theta}^{(k)})-R(\theta^{(k+1)}))$
(43)

∎

######  Remark 3.4.

In the case of $R(\cdot)\in\mathbb{R}^{d>1}$, when all components of $\mu$
increases, Equation 32 and Equation 33 imply at least one component of
$R(\cdot)$ will decrease. We defer to Remark 3.6 and Section 3.2.3 for the
discussion on when some components of $\mu$ increases while other components
decreases.

######  Remark 3.5.

Definition 3.4 and Corollary 3.1 gives a condition of decreasing $R(\cdot)$ at
the expense of increasing $\ell(\cdot)$ (increase bounded though) at iteration
$k$. When the gradient component from $\ell(\cdot)$ and $R(\cdot)$ agrees with
each other at another iteration, we could expect both objectives to decrease
as in Definition 3.3. It can also be the case at a particular iteration,
$\ell(\cdot)$ decreases at the expense of $R(\cdot)$ increases.

What is left unclear, however, is after the accumulative effect of several
iterations, if $\exists m$ to ensure Equation 30. To answer this, first,
Corollary 3.2 gives a bound on the accumulative change of $\ell(\cdot)$ with
$m=B$.

######  Corollary 3.2.

Suppose that for all $k=0,\dots,B-1$

$\displaystyle\ell(\theta^{(k+1)})+\mu^{(k+1)}R(\theta^{(k+1)})$
$\displaystyle\leq\ell(\theta^{(k)})+\mu^{(k+1)}R(\theta^{(k)})$

and $\mu^{(k+1)}>\mu^{(k)}$. Then

$\displaystyle\ell((\theta^{(B)})\leq\ell(\theta^{(0)})+S_{>}+S_{<},$ (44)

where

$\displaystyle S_{>}$
$\displaystyle=\sum_{k\in\mathcal{K}_{>}}\mu^{(k+1)}(R(\theta^{(k)})-R(\theta^{(k+1)})),$
$\displaystyle S_{<}$
$\displaystyle=\sum_{k\in\mathcal{K}_{<}}\mu^{(k)}(R(\theta^{(k)})-R(\theta^{(k+1)})),$
$\displaystyle\mathcal{K}_{>}$
$\displaystyle=\\{k\in\\{0,\dots,B-1\\}:\lx@cref{creftypecap~refnum}{eq:deteriorate_old_lambda}\text{
holds.}\\},$ $\displaystyle\mathcal{K}_{<}$
$\displaystyle=\\{0,\dots,B-1\\}\setminus\mathcal{K}_{>}.$

###### Proof.

For $k\in\mathcal{K}_{>}$ (Equation 33 holds), from Corollary 3.1 we have that

$\displaystyle\ell(\theta^{(k+1)})-\ell(\theta^{(k)})$
$\displaystyle\leq\mu^{(k+1)}(R(\theta^{(k)})-R(\theta^{(k+1)})).$

For $k\in\mathcal{K}_{<}$, it holds that

$\displaystyle\ell(\theta^{(k+1)})+\mu^{(k)}R(\theta^{(k+1)})<\ell(\theta^{(k)})+\mu^{(k)}R(\theta^{(k)}),$

implying

$\displaystyle\ell(\theta^{(k+1)})-\ell(\theta^{(k)})<\mu^{(k)}(R(\theta^{(k)})-R(\theta^{(k+1)})).$

Decomposing

$\displaystyle\ell((\theta^{(B)})-\ell(\theta^{(0)})$
$\displaystyle=\sum_{k=0}^{B-1}[\ell((\theta^{(k+1)})-\ell(\theta^{(k)})]$
$\displaystyle=\Big{(}\sum_{k\in\mathcal{K}_{>}}+\sum_{k\in\mathcal{K}_{<}}\Big{)}[\ell((\theta^{(k+1)})-\ell(\theta^{(k)})]$

yields the conclusion. ∎

Based on Corollary 3.2, we have the following conjecture to ensure Equation 30

######  Conjecture 3.1.

Further decompose $S_{<}$ from Corollary 3.2 into $S_{<}^{-}$ and $S_{<}^{+}$
via decomposing $\mathcal{K}_{<}$ in Corollary 3.2 into $\mathcal{K}_{<}^{-}$
and $\mathcal{K}_{<}^{+}$. Here $\mathcal{K}_{<}^{-}$ corresponds to the
situations of $\ell(\cdot)$ descent with $S_{<}^{-}<0$. And
$\mathcal{K}_{<}^{+}$ corresponds to the $\ell(\cdot)$ ascent cases but
upperbounded by
$S_{<}^{+}=\sum_{k\in\mathcal{K}_{<}^{+}}\mu^{(k)}(R(\theta^{(k)})-R(\theta^{(k+1)}))$.
Then $\exists$ sequence $\\{\mu^{(k)}\\}$ s.t.

$\displaystyle\ell((\theta^{(B)})\leq\ell(\theta^{(0)})+S_{>}+S_{<}^{+}+S_{<}^{-}<0$
(45)

######  Remark 3.6.

In the case of $d>1$, as long as a loss component of $R(\cdot)^{(k)}$, say
$R_{1}^{(k)}(\cdot)$ is not bounded by its corresponding setpoint component
$b_{1}^{(k)}(\cdot)$, we can still increase the value of the $\mu$ component
$\mu^{(k)}_{1}$ to give more weight to the corresponding loss component, until
the loss component in question decreases below the setpoint. However, an
increased $\mu$ component $\mu^{(k)}_{1}$ corresponds to more weight on the
gradient component corresponding to $R_{1}^{(k)}(\cdot)$, which makes it
harder for other components of $R(\cdot)$ and $\ell(\cdot)$ to decrease.
Therefore, instead of always increasing each component of $\mu$, if one
component of $R$ has overshoot over the setpoint (i.e. constraint switched
from being not satisfied to satisfied with respect to a particular loss
component), we decrease the $\mu$ value to give the other components of $\mu$
and $R$ more feasible space to adjust themselves.

######  Remark 3.7.

We could extend the conjecture above to $\exists\mu$ sequence such that after
$B$ iterations, Pareto descent in Definition 3.3 could be achieved
(equivalently, setpoint should shrink see Remark 3.3). We did observe such a
Pareto descent behavior, defined in Definition 3.3 in our experiment for the
$d>1$ case of Equation 4, e.g. see Figure 5 and Figure 6.

####  3.2.3 Output feedback PI-like controller

How exactly should the $\mu$ sequence discussed above look like? With the
adaptive law from Equation 23 for setpoint $b$, we design the following
controller for $\mu$:

$\displaystyle e^{(k)}$ $\displaystyle=R(\cdot)-b^{(k)}$ (46)
$\displaystyle\delta_{I}^{(k+1)}$
$\displaystyle=(1-\xi_{d})\delta_{I}^{(k)}+\xi_{d}e^{(k)}$ (47)
$\displaystyle\mu^{(k+1)}$
$\displaystyle=\max(\mu^{(k)}\exp^{\max(K_{I}\delta_{I}^{(k+1)},v_{sat})},\mu_{clip})$
(48) $\displaystyle K_{I}$ $\displaystyle>0$ (49) $\displaystyle K_{I}$
$\displaystyle\in\mathbb{R}^{d}$ (50) $\displaystyle k$
$\displaystyle=0,1,\cdots,B$ (51)

In Equation 46, we calculate how far the current output is away from the
setpoint $b$, which gets passed through a moving average in Equation 47 with
$\xi_{d}$ being the coefficient of moving average in Equation 47 (Rezende and
Viola, 2018). $\mu$ is the output of the controller in Equation 48, $K_{I}$ is
the control gain for PI (proportional-integration) Johnson and Moradi (2005)
like control, $K_{I}\delta_{I}^{(k+1)}$ is component wise multiplication,
$\exp$ and $\max$ are computed component wise, $v_{sat}$ is the exponential
shoulder saturation, $\mu_{clip}$ is the upper bound for $\mu$. Note that
$\mu_{clip}$ determines the dynamic range of $\mu$ while $v_{sat}$ determines
an upper bound about how fast $\mu$ changes.

######  Remark 3.8.

To avoid arbitrary choice of $K_{I}\in\mathbb{R}^{d}$, we use
$\delta^{(0)}=R^{(0)}-b^{(0)}$ in Equation 47 where $b^{(0)}=\rho R^{(0)}$ (a
percentage $\rho$). We divide by $\eta v_{sat}$ with $0<\eta<1$ to get the
value for $K_{I}$ so that the hyperparameter for our algorithm is
$\rho,\eta\in\mathbb{R}$ instead of $K_{I}\in\mathbb{R}^{d}$.

The whole process is summarized in Algorithm 1 and Figure 1.

Algorithm 1 Multi-objective Hierarchical Output Feedback Optimization

1:procedure M-Hof-Opt($\mu^{(0)},\theta^{(0)},\rho,\cdot$)

2: Initialize $\mu^{(0)}$, calculate $b^{(0)}=\rho R^{(0)}$

3: Compute $K_{I}$ based on Remark 3.8

4: while budget $B$ not reached do

5: Update $\mu$ from controller in Equation 48

6: Update $\theta$ via running Equation 14 with the updated $\mu$ value

7: Adapt setpoint according to Equation 23.

8: end while

9: return $\theta^{(B)}$ $\triangleright$

10:end procedure

####  3.2.4 Multi-objective setpoint based model selection

In general, the feasibility in Equation 26 becomes more difficult to be met
each time the setpoint $b$ shrinks. After several setpoint shrinkages, the
difficulty of attaining the feasibility can lead to an oscillating behavior of
$R(\cdot)$ around the setpoint with uncontrolled amplitude, as shown in Figure
5. To circumvent this behavior, the up-to-event best model is selected at the
last setpoint shrinkage.

##  4 Related work and Discussion

Multipliers warmup (Sønderby et al., 2016; Ilse et al., 2020) or phase-in
(Sicilia et al., 2023) requires ultimate multiplier values for each loss term,
thus necessitates hyperparameter (ultimate multiplier values) optimization.
Other works, such as (Rezende and Viola, 2018; Klushyn et al., 2019), has
introduced constraint optimization by representing the weights of loss terms
using Lagrange multipliers, designed to prevent over-regularisation of a
single term within the loss function, which is extended to multiple loss terms
in (Chen et al., 2022a). In addition, Shao et al. (2020) (only deals with 2
loss terms),Chen et al. (2022a) proposed similar multiplier adaptation schemes
as in Section 3.2.3 where the multiplier changes with a rate proportional to
how the corresponding loss component distances from a constraint bound.
However, their methods relied on a fixed constraint bound for each loss term
and fixed proportionate gain for each loss term (defined to be $K_{I}$ in
Equation 48 in our work), necessitating a hyperparameter search, which is
challenging when the number of terms in the loss increases. Moreover, for some
regularizers, such as the ones based on KL divergence, it is challenging to
estimate the constraint bound range. Specifically in the field of domain
generalization, it is not obvious how to choose the constraint bound value and
other hyperparameters due to the lack of observations from the target domain,
where the validation set from the training domain can easily reach a saturated
value due to the high expressive power of modern neural networks.

In comparison, we deal with many loss terms with multidimensional multipliers
and introduced Equation 23 to progressively adapt the constraint bound and
Remark 3.8 to choose the multidimensional controller gain value, thus avoid
exhaustive hyperparameter search. In addition to the control theory
formulation for multiplier adjustment, we also provided a Probabilistic
Graphical Model interpretation in Section 3.1 and a landscape scheduling
interpretation in Figure 4. Furthermore, we discussed our method under the
multi-objective optimization scheme compared to single-objective consideration
from other works.

Earlier multi-objective deep learning (Sener and Koltun, 2018; Lin et al.,
2019) suffers from increased computational cost with the number of losses,
explicit computation of the gradient with respect to each loss term and
gradient norms, making it practically infeasible Mahapatra and Rajan (2021)
for large neural networks and large datasets. For instance, Pareto Invariance
Risk Minimization (PAIR) (Chen et al., 2022b) studies domain generalization
with the use of the multi-objective optimization technique Exact Pareto
Optimal Search (EPOS) (Mahapatra and Rajan, 2021), which requires explicit
full gradient information for each loss term and computation of the $C$ matrix
with complexity $m^{2}n$ (Mahapatra and Rajan, 2021). With $m$ being the
number of loss terms (equivalent to $d+1$ in our paper), and $n$ the dimension
of gradient, which is super high in modern neural networks (e.g. 200 million
for ResNet50), thus results in heavy memory requirements and excessive
computation burden, and hinders the deployment to large-scale networks Chen et
al. (2022b). In comparison, our multi-objective optimization method (M-HOF-
Opt) need only $d$ scalar multiplications, thus the extra computation is
negligible and can be used in large-scale neural networks. EPOS (Mahapatra and
Rajan, 2021) also requires the user to provide a preference vector which can
be difficult to set, say with 6 objectives for our case. In addition, Chen et
al. (2022b) reported hyperparameter tuning for the preference vector and step
length is needed for the method to work properly. In comparison, our method,
M-HOF-Opt is preference vector free. It ensures the Pareto descent due to the
update rule of our setpoint. There is no requirement to search for step
length.

Note that the EPOS method Mahapatra and Rajan (2021) is searching for Exact
Pareto optimal (Mahapatra and Rajan, 2021), which might not exist. Chen et al.
(2022b) divides the training into two phases where the first phase uses
single-objective optimization favoring ERM loss while the second phase is for
balancing different loss terms. However, how many computation resources should
one set for the first phase is not clear. Furthermore, the first stage (phase)
can affect the exact Pareto optimality.

##  5 Experiments

Figure 5:  Our method drives the different components of $R(\cdot)$ loss terms
down at different scales and rates towards the setpoint, which further
promotes the setpoint shrinkage. In the Top row, we show the multiplier
dynamic as controller output signal in Equation 48 with respect to training
epochs. In the Bottom row, we present the tracking behavior of the
corresponding regularization loss $R(\cdot)$ in Equation 4 with respect to
setpoint $b$ defined in Equation 23. This figure is from training DIVA (Ilse
et al., 2020) with our multi-objective hierarchical feedback optimization
scheme on the PACS dataset (Li et al., 2017). Left column: For the
regularization loss of negative KL-divergence corresponding to multiplier
$\beta_{y}$ in Equation 8 with
$R_{\beta_{y}}=-KL(q_{\phi_{y}}(z_{y}|x)||p_{\theta_{y}}(z_{y}|y))$ (blue
curve), we see that with the setpoint (red curve) initially set to $\rho=0.99$
in Equation 21 of the initial $R_{\beta_{y}}$ value, in the first few epochs
$R_{\beta_{y}}$ overshoots to an even worse value, but our controller pulls it
back towards the setpoint in later epochs via increase in $\beta_{y}$ (top
figure in the same column) thus adding more weight on $R_{\beta_{y}}(\cdot)$
in the penalized loss. When $R_{\beta_{y}}$ is controlled to be bounded by its
setpoint, the corresponding multiplier value $\beta_{y}$ also decays to give
space for other objectives to decrease (top figure in the same column). With
coordinated $\mu$ dynamic from each component of the controller output, the
setpoint shrinks when all $R(\cdot)$ components Pareto dominates the setpoint.
After the setpoint shrinks several times, feasibility in Equation 26 becomes
harder to maintain and our controller still keeps $R_{\beta_{y}}$ oscillating
around the setpoint. Our model selection criteria in Section 3.2.4 selects the
model at the last occurrence of the setpoint shrinkage, thus the oscillating
behavior at the later epochs can be ignored. Middle column: For the domain
classification loss $R_{\gamma_{d}}$ corresponding to multiplier $\gamma_{d}$
in Equation 10 (blue curve), we observe a similar behavior, but the setpoint
shrinkage is more visible. Right column: Although the output component for the
loss $R_{\mu_{recon}}$ in Equation 6 oscillates with relatively large
amplitude after the last setpoint shrinkage, our model selection criteria at
Section 3.2.4 will exclude this oscillating behavior.  Figure 6: Our training
scheme adapts the multi-dimensional multiplier $\mu$, steering the
optimization process towards a configuration that balances the trade-off
between minimizing output $\ell(\cdot)=\mathbb{E}_{q_{\phi_{y}}(z_{y}|x)}[\log
q_{\omega_{y}}(y|z_{y})]$ in Equation 5 and output
$R_{\gamma_{d}}(\cdot)=\mathbb{E}_{q_{\phi_{d}}(z_{d}|x)}[\log
q_{\omega_{d}}(d|z_{d})]$ in Equation 10. In this plot, the solid large point
corresponds to the initial output. We use gradually changing colors to
indicate the training iterations (epochs), indicated by the color bar. For
improved visualization, we exclusively utilized the initial 120 epochs and
plotted data points solely for every 10 epochs. This figure corresponds to
Figure 5.

###  5.1 Experimental Setting

Figure 7: Our automatic multiplier adjustment scheme ensures robust out-of-
domain test accuracy when training DIVA (Ilse et al., 2020) with the modified
ResNet50 from (Gulrajani and Lopez-Paz, 2020) on PACS dataset (testing on
domain sketch and training on domain photo, art-painting and cartoon). Each
panel (i.e., subplot) corresponds to a training scheme. Each box plot inside a
panel corresponds to a specific hyperparameter (e.g. multiplier value or
controller hyperparameters) combination. Inside each box plot, we repeat the
experiment with different random seeds corresponding to the dots scatter. In
the first panel, we warm up only $\beta_{y}$, $\beta_{d}$ in Equations 9 and 8
to their ultimate value $1$ while sampling fixed multiplier $\gamma_{d}$ and
$\gamma_{y}$ while fixing $\mu_{recon}=1.0$. In the second panel,
additionally, we also warm up $\mu_{recon}$ in Equation 6 to $1.0$ and warm up
$\gamma_{d}$ in Equation 10 to the sampled value, and keep $\gamma_{y}$ fixed
to be the sampled value. In the third panel, we use fixed constant multipliers
without warm-up while keep $\beta=1$ and $\mu_{recon}=1$ and only sample
$\gamma_{d},\gamma_{y}$ combinations. The last two panels correspond to
results utilizing our multi-objective hierarchical feedback optimization
training method, and we let $\mu$ include all multipliers except
$\gamma_{y}=1$ in Equation 5. For simplicity, when we write $\mu^{(0)}=1.0$
where $\mu^{(0)}\in\mathbb{R}_{+}^{d}$, we mean each component of $\mu^{(0)}$
equals $1.0$.

This section presents our experiments training the domain generalization model
DIVA (Ilse et al., 2020) with $6$ loss terms defined in Equation 4 on the
widely used domain generalization benchmark dataset PACS (Li et al., 2017)
with different training schemes. Due to the excessive memory requirements and
heavy computation of earlier multi-objective deep learning methods discussed
in Section 4, which hinders their application, we can only compare our
training scheme with other multiplier scheduling techniques. Given that the
sketch domain within the PACS dataset is considered the most inherently
challenging, it was chosen as the sole leave-one-out domain to test the
performance of domain generalization. We run experiments to answer the
following questions:

  * •

Is our multiplier controller capable of adjusting multidimensional multipliers
automatically to drive each component of $R(\cdot)$ at different scales to the
setpoint? Does the setpoint shrink, which implies multi-objective descent? We
answer these questions with Figure 5 and Figure 6.

  * •

Will wrong multiplier choices for warmup and fixed multiplier training schemes
result in catastrophic effects? Although our method in contrast is multiplier-
free, it does need multiplier initial condition $\mu^{(0)}$ and controller
hyperparameters. Does varying those have detrimental effects on the
predicative performance as well? We answer this question with the benchmark
setting in Section 5.1.1 and the results presented in Figure 7.

####  5.1.1 Benchmark setting

In the benchmark, we sample different combinations of hyperparameters,
(including multiplier for baselines and controller hyperparameter for our
method) for each method to be compared. For baselines, we sample
$\gamma_{d}\in\\{1,1001,100001\\}$ and $\gamma_{y}\in\\{1,1001,100001\\}$
value combinations, and set the ultimate values of $\beta_{y},\beta_{d}$ to be
$1.0$, following (Ilse et al., 2020). Our method does not need to set ultimate
values for $\mu$, so we sample controller hyperparameters, as well as
different initial conditions $\mu^{(0)}$ (Each component of the $\mu^{(0)}$
vector set to be the same value).

###  5.2 Results and Discussion

Regarding the dynamics of how the setpoint $b$ decreases: Initially, Equation
27 holds, in which case we do not update the setpoint, even if some but not
all the blue curves (corresponding to $R(\cdot)$ components) are below the red
curves (components of the setpoint). After all constraints for each component
of $R(\cdot)$ are satisfied as described in Equation 26, the setpoint (red
curves) also gets adapted to a new value (which we term shrinkage), as a new
goal to be reached by the $R(\cdot)$ loss in Equation 10. Note that the
shrinkage of setpoint depends on Pareto dominance in Equation 24a, while non-
dominance defined an equivalent relation in Definition 2.2. Our hierarchical
output feedback optimization training scheme produced the dynamics of
$R(\cdot)$ shown in Figure 5 together with the corresponding setpoints $b$ and
multipliers. Note that the multiplier $\beta_{y},\gamma_{d},\mu_{recon}$ and
their corresponding loss components operate at different numerical ranges, but
our controller still manages to drive the different $R(\cdot)$ loss terms down
at different scales and rates.

Figure 6 shows the output portrait of $\ell(\cdot)$ versus a component of
$R(\cdot)$ of our training scheme in correspondence to Figure 5. With the
multi-objective setpoint based model selection criteria in Section 3.2.4, our
method admits a multi-objective descent of the selected model compared to the
initial output, such that the uncontrolled behavior at the end of iterations
can be safely ignored.

From Figure 7, we can observe that the generalization performance of baselines
is highly sensitive to hyperparameter selection. In contrast, our method, as
supported by the adaptive behaviors visualized in Figure 5, automatically
adjusts the multipliers during the training process, leading to a more robust
performance with respect to changes of controller hyperparameters and
multiplier initial condition $\mu^{(0)}$. Our method does not need
hyperparameter searching of the multiplier, thus saves a fair amount of
computational resources.

##  6 Conclusion

This work addresses the issue of combinatorial choice for multidimensional
multipliers weighting many terms in neural network loss function. We proposed
a probabilistic graphical model for joint multiplier and model parameter
optimization with respect to a multi-objective descent of all loss terms. The
estimation of model parameter leads to an automatic adjustment scheme for the
multipliers, adopting a hierarchical control scheme, which breaks the multi-
objective optimization problem into a series of constraint optimization sub-
problems. Each sub-problem is configured with a self-adaptive multi-objective
setpoint updated via Pareto dominance. A PI-like multiplier controller drives
the loss term to satisfy the setpoint constraint for each sub-problem. Our
method operates at the timescale of epoch level during training, thus
circumvents the need for exhaustive multiplier search and saves tremendous
computational resources compared to methods like Bayesian optimization. It
also circumvents the excessive memory requirements and heavy computational
burden of existing multi-objective deep learning methods. Our method
demonstrates robust generalization performance in the PACS domain
generalization task with Domain Invariant Variational Autoencoding Ilse et al.
(2020), compared to other multiplier choice or scheduling schemes which
produce very unstable behavior in the training due to the need for a
combinatorial choice of multipliers. Our method is multiplier-free since the
multiplier is the output of our controller, and it has robust performance with
respect to different multiplier initial values. For future work, we are
interested in testing the performance of our hierarchical feedback
optimization on other domain generalization methods, as well as on other deep
learning fields in the form of Equation 1. Theoretically, it is interesting to
validate Conjecture 3.1, derive sufficient conditions, incorporate more
mathematical analysis into multidimensional case for $R$, and investigate the
potential $\mu$ convergence behavior.

##  7 Contributions And Acknowledgement

XS proposed the idea, developed the theories and algorithms, implemented the
code, designed and carried out the experiments, processed experimental data
and generated the figures, wrote the manuscript. NC proposed the PI-like
multiplier adaptation in Section 3.2.3, helped XS with PI-like controller
implementation, held various discussions, proofread and improved the
manuscript. AG initiated the visualization code for Figure 5 and Figure 6,
refined by XS, proofread and improved the manuscript. YX developed Corollary
3.2 and corresponding remarks with XS, proofread and improved the constraint
optimization part of the theory, discussed on control theory part. CF
initiated the benchmark code with refinements from XS, proofread Corollary
3.1. ED and FD tested and improved the code for benchmark, proofread the
manuscript. FD further improved the tensorboard visualization of the
algorithm. LB proofread Corollary 3.1, tested the benchmark code. DS tested
the benchmark code, assisted XS in collecting experimental results. CM
supervised the project, offered intensive and in depth proofreads of the
manuscript with major quality improvements.

We sincerely thank Dr. Changxin Liu and Dr. Prof. Jann Rolfes for spending
hours on proofreading the manuscript and offering many helpful suggestions to
the manuscript. We thank Tomas Chobola for Inkscape support. We thank the
anonymous reviewers for offering suggestions to enhance the readability of our
paper.

## References

  * Bazaraa et al. [2013] Mokhtar S Bazaraa, Hanif D Sherali, and Chitharanjan M Shetty.  _Nonlinear programming: theory and algorithms_.  John wiley & sons, 2013.
  * Bertsekas [2014] Dimitri P Bertsekas.  _Constrained optimization and Lagrange multiplier methods_.  Academic press, 2014.
  * Carlucci et al. [2019] Fabio M Carlucci, Antonio D’Innocente, Silvia Bucci, Barbara Caputo, and Tatiana Tommasi.  Domain generalization by solving jigsaw puzzles.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 2229–2238, 2019.
  * Chen et al. [2022a] Nutan Chen, Patrick van der Smagt, and Botond Cseke.  Local distance preserving auto-encoders using continuous knn graphs.  In _Topological, Algebraic and Geometric Learning Workshops 2022_ , pages 55–66, 2022a.
  * Chen et al. [2016] Xi Chen, Yan Duan, Rein Houthooft, John Schulman, Ilya Sutskever, and Pieter Abbeel.  Infogan: Interpretable representation learning by information maximizing generative adversarial nets.  In _Advances in neural information processing systems_ , pages 2172–2180, 2016.
  * Chen et al. [2022b] Yongqiang Chen, Kaiwen Zhou, Yatao Bian, Binghui Xie, Bingzhe Wu, Yonggang Zhang, Kaili Ma, Han Yang, Peilin Zhao, Bo Han, et al.  Pareto invariant risk minimization: Towards mitigating the optimization dilemma in out-of-distribution generalization.  _arXiv preprint arXiv:2206.07766_ , 2022b.
  * Dietterich et al. [1998] Thomas G Dietterich et al.  The maxq method for hierarchical reinforcement learning.  In _ICML_ , volume 98, pages 118–126, 1998.
  * Doyle et al. [2013] John C Doyle, Bruce A Francis, and Allen R Tannenbaum.  _Feedback control theory_.  Courier Corporation, 2013.
  * Ganin et al. [2016] Yaroslav Ganin, Evgeniya Ustinova, Hana Ajakan, Pascal Germain, Hugo Larochelle, François Laviolette, Mario Marchand, and Victor Lempitsky.  Domain-adversarial training of neural networks.  _The journal of machine learning research_ , 17(1):2096–2030, 2016.
  * Garnett [2023] Roman Garnett.  _Bayesian optimization_.  Cambridge University Press, 2023.
  * Guerreiro et al. [2021] Andreia P Guerreiro, Carlos M Fonseca, and Luís Paquete.  The hypervolume indicator: Computational problems and algorithms.  _ACM Computing Surveys (CSUR)_ , 54(6):1–42, 2021\.
  * Gulrajani and Lopez-Paz [2020] Ishaan Gulrajani and David Lopez-Paz.  In search of lost domain generalization.  _arXiv preprint arXiv:2007.01434_ , 2020.
  * Hauswirth et al. [2024] Adrian Hauswirth, Zhiyu He, Saverio Bolognani, Gabriela Hug, and Florian Dörfler.  Optimization algorithms as robust feedback controllers.  _Annual Reviews in Control_ , 57:100941, 2024.
  * He et al. [2023] Zhiyu He, Saverio Bolognani, Jianping He, Florian Dörfler, and Xinping Guan.  Model-free nonlinear feedback optimization.  _IEEE Transactions on Automatic Control_ , 2023.
  * Hutter et al. [2007] Frank Hutter, Holger H Hoos, and Thomas Stützle.  Automatic algorithm configuration based on local search.  In _Aaai_ , volume 7, pages 1152–1157, 2007.
  * Hutter et al. [2014] Frank Hutter, Manuel López-Ibánez, Chris Fawcett, Marius Lindauer, Holger H Hoos, Kevin Leyton-Brown, and Thomas Stützle.  Aclib: A benchmark library for algorithm configuration.  In _Learning and Intelligent Optimization: 8th International Conference, Lion 8, Gainesville, FL, USA, February 16-21, 2014. Revised Selected Papers 8_ , pages 36–40. Springer, 2014.
  * Ilse et al. [2020] Maximilian Ilse, Jakub M Tomczak, Christos Louizos, and Max Welling.  Diva: Domain invariant variational autoencoders.  In _Medical Imaging with Deep Learning_ , pages 322–348. PMLR, 2020\.
  * Johnson and Moradi [2005] Michael A Johnson and Mohammad H Moradi.  _PID control_.  Springer, 2005.
  * Kingma and Ba [2015] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _International Conference for Learning Representations_ , 2015.
  * Klushyn et al. [2019] A. Klushyn, N. Chen, R. Kurle, B. Cseke, and P. van der Smagt.  Learning hierarchical priors in VAEs.  _Advances in Neural Information processing Systems_ , 32, 2019.
  * Koller and Friedman [2009] Daphne Koller and Nir Friedman.  _Probabilistic graphical models: principles and techniques_.  MIT press, 2009.
  * Kulkarni et al. [2016] Tejas D Kulkarni, Karthik Narasimhan, Ardavan Saeedi, and Josh Tenenbaum.  Hierarchical deep reinforcement learning: Integrating temporal abstraction and intrinsic motivation.  _Advances in neural information processing systems_ , 29, 2016.
  * Levi et al. [2021] Matan Levi, Idan Attias, and Aryeh Kontorovich.  Domain invariant adversarial learning.  _arXiv preprint arXiv:2104.00322_ , 2021.
  * Levine [2018] Sergey Levine.  Reinforcement learning and control as probabilistic inference: Tutorial and review.  _arXiv preprint arXiv:1805.00909_ , 2018.
  * Lewis et al. [2012] Frank L Lewis, Draguna Vrabie, and Vassilis L Syrmos.  _Optimal control_.  John Wiley & Sons, 2012.
  * Li et al. [2017] Da Li, Yongxin Yang, Yi-Zhe Song, and Timothy M Hospedales.  Deeper, broader and artier domain generalization.  In _Proceedings of the IEEE international conference on computer vision_ , pages 5542–5550, 2017.
  * Li et al. [2018] Da Li, Yongxin Yang, Yi-Zhe Song, and Timothy Hospedales.  Learning to generalize: Meta-learning for domain generalization.  In _Proceedings of the AAAI conference on artificial intelligence_ , volume 32, 2018.
  * Lin et al. [2019] Xi Lin, Hui-Ling Zhen, Zhenhua Li, Qing-Fu Zhang, and Sam Kwong.  Pareto multi-task learning.  _Advances in neural information processing systems_ , 32, 2019.
  * López-Ibáñez et al. [2016] Manuel López-Ibáñez, Jérémie Dubois-Lacoste, Leslie Pérez Cáceres, Mauro Birattari, and Thomas Stützle.  The irace package: Iterated racing for automatic algorithm configuration.  _Operations Research Perspectives_ , 3:43–58, 2016.
  * Mahapatra and Rajan [2021] Debabrata Mahapatra and Vaibhav Rajan.  Exact pareto optimal search for multi-task learning: Touring the pareto front.  _arXiv e-prints_ , pages arXiv–2108, 2021.
  * Milne [2019] Tristan Milne.  Piecewise strong convexity of neural networks.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Picallo et al. [2020] Miguel Picallo, Saverio Bolognani, and Florian Dörfler.  Closing the loop: Dynamic state estimation and feedback optimization of power grids.  _Electric Power Systems Research_ , 189:106753, 2020.
  * Probst et al. [2019] Philipp Probst, Anne-Laure Boulesteix, and Bernd Bischl.  Tunability: Importance of hyperparameters of machine learning algorithms.  _The Journal of Machine Learning Research_ , 20(1):1934–1965, 2019.
  * Rame et al. [2022] Alexandre Rame, Corentin Dancette, and Matthieu Cord.  Fishr: Invariant gradient variances for out-of-distribution generalization.  In _International Conference on Machine Learning_ , pages 18347–18377. PMLR, 2022.
  * Reddi et al. [2019] Sashank J Reddi, Satyen Kale, and Sanjiv Kumar.  On the convergence of adam and beyond.  _arXiv preprint arXiv:1904.09237_ , 2019.
  * Rezende and Viola [2018] D. J. Rezende and F. Viola.  Taming VAEs.  _CoRR_ , 2018.
  * Sener and Koltun [2018] Ozan Sener and Vladlen Koltun.  Multi-task learning as multi-objective optimization.  _Advances in neural information processing systems_ , 31, 2018.
  * Shao et al. [2020] Huajie Shao, Shuochao Yao, Dachun Sun, Aston Zhang, Shengzhong Liu, Dongxin Liu, Jun Wang, and Tarek Abdelzaher.  Controlvae: Controllable variational autoencoder.  In _International Conference on Machine Learning_ , pages 8655–8664. PMLR, 2020.
  * Sicilia et al. [2023] Anthony Sicilia, Xingchen Zhao, and Seong Jae Hwang.  Domain adversarial neural networks for domain generalization: When it works and how to improve.  _Machine Learning_ , pages 1–37, 2023.
  * Sønderby et al. [2016] Casper Kaae Sønderby, Tapani Raiko, Lars Maaløe, Søren Kaae Sønderby, and Ole Winther.  Ladder variational autoencoders.  _Advances in neural information processing systems_ , 29, 2016.
  * Sun and Bischl [2019] Xudong Sun and Bernd Bischl.  Tutorial and survey on probabilistic graphical model and variational inference in deep reinforcement learning.  In _2019 IEEE Symposium Series on Computational Intelligence (SSCI)_ , number 1908.09381, 2019.
  * Sun and Buettner [2021] Xudong Sun and Florian Buettner.  Hierarchical variational auto-encoding for unsupervised domain generalization.  _ICLR 2021 RobustML, arXiv preprint arXiv:2101.09436_ , 2021.
  * Sun et al. [2019a] Xudong Sun, Andrea Bommert, Florian Pfisterer, Jörg Rahnenführer, Michel Lang, and Bernd Bischl.  High dimensional restrictive federated model selection with multi-objective bayesian optimization over shifted distributions.  _arXiv preprint arXiv:1902.08999_ , 2019a.
  * Sun et al. [2019b] Xudong Sun, Alexej Gossmann, Yu Wang, and Bernd Bischl.  Variational resampling based assessment of deep neural networks under distribution shift.  In _2019 IEEE Symposium Series on Computational Intelligence (SSCI)_ , number 1906.02972, 2019b.
  * Sun et al. [2020] Xudong Sun, Jiali Lin, and Bernd Bischl.  Reinbo: Machine learning pipeline conditional hierarchy search and configuration with bayesian optimization embedded reinforcement learning.  In _Machine Learning and Knowledge Discovery in Databases: International Workshops of ECML PKDD 2019, Würzburg, Germany, September 16–20, 2019, Proceedings, Part I_ , pages 68–84. Springer International Publishing, 2020.
  * Zitzler and Künzli [2004] Eckart Zitzler and Simon Künzli.  Indicator-based selection in multiobjective search.  In _International conference on parallel problem solving from nature_ , pages 832–842. Springer, 2004.
  * Zitzler et al. [2007] Eckart Zitzler, Dimo Brockhoff, and Lothar Thiele.  The hypervolume indicator revisited: On the design of pareto-compliant indicators via weighted integration.  In _Evolutionary Multi-Criterion Optimization: 4th International Conference, EMO 2007, Matsushima, Japan, March 5-8, 2007. Proceedings 4_ , pages 862–876. Springer, 2007.