Recently, there has been an increasing interest in exploring the application of multiobjective optimization (MOO) in machine learning (ML). The interest is driven by the numerous situations in real-life applications where multiple objectives need to be optimized simultaneously. A key aspect of MOO is the existence of a Pareto set, rather than a single optimal solution, which illustrates the inherent trade-offs between objectives. Despite its potential, there is a noticeable lack of satisfactory literature that could serve as an entry-level guide for ML practitioners who want to use MOO. Hence, our goal in this paper is to produce such a resource. We critically review previous studies, particularly those involving MOO in deep learning (using Physics-Informed Neural Networks (PINNs) as a guiding example), and identify misconceptions that highlight the need for a better grasp of MOO principles in ML. Using MOO of PINNs as a case study, we demonstrate the interplay between the data loss and the physics loss terms. We highlight the most common pitfalls one should avoid while using MOO techniques in ML. We begin by establishing the groundwork for MOO, focusing on well-known approaches such as the weighted sum (WS) method, alongside more complex techniques like the multiobjective gradient descent algorithm (MGDA). Additionally, we compare the results obtained from the WS and MGDA with one of the most common evolutionary algorithms, NSGA-II. We emphasize the importance of understanding the specific problem, the objective space, and the selected MOO method, while also noting that neglecting factors such as convergence can result in inaccurate outcomes and, consequently, a non-optimal solution. Our goal is to offer a clear and practical guide for ML practitioners to effectively apply MOO, particularly in the context of DL.

multiobjective optimization, machine learning, physics informed neural networks,
pareto front

68Q25, 68R10, 68U05

§ INTRODUCTION

Over the past decade, machine learning (ML) has seen remarkable growth and widespread adoption, both in research and applications. This surge is partly attributed to artificial neural networks (ANNs) because of their ability to model any continuous function [15, 39], given a sufficiently large network architecture. Another major catalyst for this growth owes to the availability of powerful computational resources.

By adding multiple layers to ANNs, we acquire a deep neural network (DNN) and hence transition into deep learning (DL), enabling the extraction of complex patterns from data. Generally speaking, the task in ML/DL is to train a model using a training data set, such that it is capable of making accurate predictions and generalizes to unseen data. This task translates to optimizing a performance metric, also referred to as loss function or an objective. There are numerous very efficient single objective optimization (SOO) techniques [53] for DNN training, the ADAM [33] algorithm being the most popular one.

For example, given the data about a product that is being produced in a factory, one could use DL with SOO to predict the parameters in the production line which can reduce the production costs of that product. This can be practically achieved by encoding the production cost as an objective to be minimized and using a gradient-based optimization to iteratively find a descent direction towards the minimum of the objective.

In classical optimization, the consideration of multiple, conflicting criteria is a long-known challenge. There are many examples where multiobjective optimization (MOO) is very important, for instance in production processes, where one wants to maximize the quality of a product while minimizing its cost. Hence, instead of a single optimal solution, there exists a set of optimal trade-offs called Pareto set [41]. The Pareto set represents the set of points for which no other feasible solutions exist that can improve one objective without degrading at least one of the other objectives. Areas where MOO plays a huge role include finance, economics, health systems, logistics, agriculture, and applied engineering.

As one can imagine, bringing MOO techniques and ML together is an exciting and promising research direction. Conflicting criteria naturally emerge in DL, such as the need to minimize loss, incorporate domain knowledge, reduce architectural size, or tackle distinct tasks in multi-task learning. Recently, this field has begun to gain momentum, supported by several studies [49, 6, 30, 54, 1]. However, there is a lack of accessible resources specifically tailored to introduce ML practitioners, to the application of MOO methods. On the contrary, we encountered studies that, upon close examination, reveal some common misconceptions in the application of MOO within DL. A notable example is the work [49], in which Pareto fronts for Physics-Informed Neural Networks (PINNs) [32, 34] have been produced using the weighted sum (WS) method. Our discussion will reveal that the conclusions drawn in this study may not fully align with the principles of MOO. Specifically, we observe potential gaps in the interpretation of the Pareto front and the methodology used to obtain it. In a related study [54] employing a multiobjective gradient descent algorithm (MGDA) for multi-task optimization, the network's sophisticated architecture effectively mitigated task conflicts such that the problem possesses a single solution that optimizes both tasks. These observations highlight the need for a more nuanced understanding of MOO in DL, as misinterpretations can lead to suboptimal decision-making and affect the validity of the identified trade-offs.

PINNs represent a compelling fusion of domain knowledge with data-driven approaches for tackling problems with complex dynamics. Such a fusion is achieved by defining two loss functions, the data loss representing the traditional supervised learning loss, and the physics loss, which integrates domain expertise through partial differential equations (PDEs). Despite the huge success in solving different systems like multi-physics systems [14], fluid dynamics [10], fractional PDEs [45], integro differential equations [64], parametrized PDEs [47], or stochastic PDEs [52], PINNs face several issues like hyper-parameter tuning [7] and slow convergence. The study [49] aims to address these issues by using MOO methods to optimize both the data and the physics losses. Our analysis suggests that there are aspects of their methodology and interpretation of results that demand a deeper investigation. The goal of this paper is to improve the understanding of how MOO can be used for DL problems in general, and PINN training in particular.

§.§ Current literature in the area of Multiobjective Deep Learning

The landscape of MOO in DL is marked by significant contributions across various methodologies and applications. A notable paper in the context of multi-task learning (MTL) is [54], where one neural network is trained to solve two learning tasks using MGDA. A MOO algorithm for DNNs, employing a modified Weighted Chebyshev scalarization to simplify the optimization of multiple tasks into sequential single-objective problems, is proposed in [30]. A model-agnostic method named Adamize [42] adapts the benefits of the Adam optimizer to multiobjective problems. Additionally, efforts have been made to develop approaches that generate the entire Pareto front in a single run [50]. For an in-depth analysis of MTL in computer vision, one can look into [58]. In addressing the balancing act between different losses in PINNs, several studies have explored various methodologies. While gradient-based approaches [49, 7, 63] offer one avenue for optimizing this trade-off, evolutionary algorithms present a gradient-free alternative strategy. These have found applications not just in PINNs but across a broader spectrum of deep learning applications [37, 31, 59, 12]. However, evolutionary algorithms tend to be very expensive, so the mentioned studies are mostly limited to smaller networks. Instead, they are popular for hyperoptimization tasks such as neural architecture search [38].

In contrast to multiobjective deep learning, people have also started using deep learning as an alternative method to accelerate the solution of MOPs [35, 44].

§.§ Organization

We begin in section <ref> with the fundamentals of MOO, giving a brief overview of various methods and focusing on two methods in more detail, the WS and the MGDA. In section <ref>, we delve into the basics of DL and explore the rationale behind employing MOO for DL challenges. Using PINNs as a case study, we frame the optimization of data and physics loss as an MOO problem. Section <ref> is dedicated to the key pitfalls we identified in applying MOO to DL. We present the results from our experiments of training PINNs using MOO for the logistic and the heat equation. The key features of the Pareto fronts obtained are discussed in detail. Apart from highlighting common pitfalls in MOO for DL, this section also brings to light additional challenges we encountered during our research. Finally, we conclude our study by summarizing our findings in section <ref>.

Our ultimate goal is for this document to serve as a beginner-friendly guide, assisting DL practitioners in effectively incorporating MOO into their research.

In most of the cases, ML practitioners are supposed to optimize a single criterion e.g., the cost of a product. By now there are numerous very efficient single objective optimization (SOO) techniques [53] to optimize an objective in a ML setting. However, often times, one encounters scenarios where it is required to optimize more than one objective. It is customary to simply add the different objective functions to convert the MOO problem into a SOO. While the straightforward addition of multiple objectives may seem appealing, it is crucial to acknowledge that this approach is only meaningful if the objectives share a common minimum, such that they have a common descent direction. If the objectives are conflicting, (i.e., do not have common minimum) then simply adding the objectives can land us in the wrong solution space <cit.>.

Consider the task of minimizing the cost of a product while maximizing its quality. It is clear that both of these objectives are conflicting since producing a better quality product would increase the cost. In other words, there is a trade off while optimizing the objectives simultaneously. Such a problem has to be dealt more carefully and one need to rely on the MO treatment. Hence, instead  a single optimal solution, there exist set of solutions called Pareto Set <cit.>. A Pareto set represents the set of solutions for which no other feasible solutions exist that can improve one objective without degrading the other. When plotted in the N dimensional objective space, we obtain the Pareto front which is a $N-1$ dimensional curve. Each point on the Pareto front represents a different trade-off among the objectives. As one can imagine, MOO has wide applications in those sectors where one has to optimize conflicting objectives. Some of the common sectors where MOO plays a huge role are finance, economics, health systems, logistics, agriculture, applied engineering etc.

The MOO techniques out there can be classified into two main broad categories, evolutionary algorithms <cit.> and gradient based approaches <cit.>. The Evolutionary algorithms iterate by initiating a population of  potential solutions and iteratively refining them through processes analogous to crossover, mutation, and selection found in natural genetics. Hence they are also refereed to as Genetic algorithms. Such a learning strategy enable evolutionary algorithms to converge towards the Pareto set, rather than a single point. A major drawback of using evolutionary algorithms lies in the increasing computational cost in higher dimensions. Furthermore, they are not mathematically well defined and hence do not provide any useful insight into the learning process and the underlying problem. Instead, most of the time they are used as a black box. More robust and mathematically well defined methods are the gradient based method, where one calculates the descent direction during the optimization process.

Gradient based MOO methods are the MO counterparts/versions of the single objective optimization methods e.g., steepest descent, the projected gradient, the Newton methods etc. In these methods, one calculates the descent direction at each iteration by solving a real-valued optimization problem. An obvious advantage of using such methods is that one can establish that the solution converges to a Pareto Optimal Solution. Within gradient bases approaches, a variety of methods exist which differ in performance and the nature of the problems at hand.

One of the popular strategies among gradient based methods is the scalarization approach <cit.> where one solves many instances of parameterized single-objective optimization problems. For example one can use the weighted sum approach where multiple objectives are combined into a single aggregated objective using weighted sums. The optimization is then done on this aggregated objective. By varying the weights, different points on the Pareto front can be traced. Due to various limitation of the traditional scalarization methods, one might opt for more advanced approaches like continuation methods <cit.> and multiobjective descent methods <cit.>.

As one could imagine, bringing MOO techniques to ML together is an exciting and promising adventure and to no surprise, recently this field has started gaining momentum <cit.>. However, to our surprise, we were unable to find any comprehensive and accessible resources specially tailored to introduce ML practitioners, especially those with background in DL, to the application of MOO methods. On the contrary, we encountered studies that, upon close examination, reveal some common misconceptions in the application of MOO within DL. A notable example is the work by Rohrhofer et al. [49], which claims to have identified the Pareto Front for Physics-Informed Neural Networks (PINNs) [32, 34]. Our analysis suggests that the conclusions drawn in this study may not fully align with the principles of MOO. Specifically, we observe potential gaps in the interpretation of the Pareto Front and the methodology used to obtain it. These observations highlight a critical need for a more nuanced understanding of MOO in DL, as misinterpretations can lead to sub-optimal decision-making and affect the validity of the trade-offs identified.

PINNs represent a compelling fusion of domain knowledge with data-driven approaches for tackling complex systems. Such a fusion is achieved by defining two loss functions, data loss representing the traditional supervised learning loss, and the physics loss, which integrates domain expertise through partial differential equations (PDEs). Despite the huge success in solving different systems like multi-physics systems[14], fluid dynamics[10], fractional PDE[45], integro differential equations[64], parametrized PDEs [47], stochastic PDEs [52] etc, PINNs face some issues like hyper-parameter tuning <cit.>  and convergence issues [49]. The study by Rohrhofer et al. [49] aims to address these issues by using MOO methods to optimize both the data and the physics losses. Our analysis suggests that while this approach is innovative, there are aspects of their methodology and interpretation of results that demand a deeper investigation. This paper aims to clarify these issues, providing insights that improve how MOO is understood and used with PINNs.

§ MULTIOBJECTIVE OPTIMIZATION

In this section, we give a brief introduction to MOO. In-depth discussions can be found in, e.g., [41, 18].

§.§ Foundations of multiobjective optimization

The goal of MOO is to minimize multiple objective functions. Opposing classical literature, we use the notation $\Loss_i$ instead of $f_i$ to denote objective functions in this paper, since we are interested in the minimization of losses in the context of machine learning. The decision variables which are the trainable weights of a deep neural network are denoted by $\theta$ since the variable $\mathbf{x}$ is reserved for input data. The general MOO problem reads

\begin{align}
\min_{\theta \in \R^n}\left[ \begin{array}{c}
\Loss_1(\theta)  \\
\vdots \\
\Loss_m(\theta)
\end{array}\right],
\label{eq:MOP}
\tag{MOP}
\end{align}

and formalizes the problem of jointly minimizing the objective functions $\Loss_i$ by choosing an appropriate decision vector $\theta \in \R^n$. Using the vector-valued function $\mathbfcal{L}: \R^n \to \R^m, \theta \mapsto \left( \Loss_1(\theta),\dots, \Loss_m(\theta)\right)^{\top}$ problem (<ref>) can also be written as

\begin{align*}
\min_{\theta \in \R^n} \mathbfcal{L}(\theta).
\end{align*}

The space of inputs $\R^n$ is called the decision space and the space $\R^m$ is called the image space. The vectors $\theta \in \R^n$ are decision vectors. We call the vectors $\mathbfcal{L}(\theta) \in \R^m$ the objective values and the set $\mathbfcal{L}(\R^n) \subseteq \R^m$ the attainable set. In many settings, it is desirable to constrain the decision space to a subset $\Theta \subseteq \R^n$. However, we only consider the unconstrained case here, as this is most common in deep learning.

Up to this point, it is not clear what is meant by solving (<ref>). In general, there is no decision vector $\theta^*$ that minimizes all objectives at the same time. Since there is no (natural) total order on the image space $\R^m$ we need a suitable definition of optimality for (<ref>).

(Pareto optimality and Pareto set)

i) A decision vector $\theta^* \in \R^n$ is (globally) Pareto optimal if there does not exist another decision vector $\theta \in \R^n$ such that $\Loss_i(\theta) \le \Loss_i(\theta^*)$ for all $i = 1,\dots,m$ and $\Loss_j(\theta) < \Loss_j(\theta^*)$ for some index $j \in \{1,\dots,m\}$.

ii) The set $\mathcal{P} \coloneqq \left\lbrace \theta \in \R^n \, : \, \theta \text{ is Pareto optimal } \right\rbrace \subseteq \R^n$ is called the Pareto set and the set $\mathcal{F} \coloneqq \mathbfcal{L}(\mathcal{P}) \subseteq \R^m$ is called the Pareto front.

iii) A decision vector $\theta^* \in \R^n$ is locally Pareto optimal if it satisfies i) on a neighborhood of $\theta^*$.

(Weak Pareto optimality and weak Pareto set)

i) A decision vector $\theta^* \in \R^n$ is (globally) weakly Pareto optimal if there does not exist another decision vector $\theta \in \R^n$ such that $\Loss_i(\theta) < \Loss_i(\theta^*)$ for all $i = 1,\dots,m$.

ii) The set $\mathcal{P}_w \coloneqq \left\lbrace \theta \in \R^n \, : \, \theta \text{ is weakly Pareto optimal } \right\rbrace \subseteq \R^n$ is called the weak Pareto set and the set $\mathcal{F}_w \coloneqq \mathbfcal{L}(\mathcal{P}_w) \subseteq \R^m$ is called the weak Pareto front.

iii) A decision vector $\theta^* \in \R^n$ is locally weakly Pareto optimal if it satisfies i) on a neighborhood of $\theta^*$.

Figure <ref> provides a visualization of the definition of Pareto optimal points. The Figure shows a schematic plot of the objective function values for an MOP with two objectives. The gray area is the attainable set. The green boundary on the lower left is the Pareto front. The elements $\theta_2$ and $\theta_3$ are Pareto optimal since there are no other objective function values to the lower left of $\mathbfcal{L}(\theta_2)$ and $\mathbfcal{L}(\theta_3)$ in the sense of Definition <ref>. On the other hand, $\theta_1$ is not Pareto optimal since $\mathbfcal{L}(\theta_2) < \mathbfcal{L}(\theta_1)$ which is clarified by the cone starting in $\mathbfcal{L}(\theta_1)$. In this setting, we say that decision vector $\theta_2$ dominates $\theta_1$. We use the same terminology for the objective vectors, i.e., $\mathbfcal{L}(\theta_2)$ dominates $\mathbfcal{L}(\theta_1)$.

For practical applications, it is impossible to verify whether a given point satisfies the conditions in Definitions <ref> and <ref> directly since one would have to compare infinitely many points. If the objective functions $\Loss_i$ are continuously differentiable we can use the following necessary condition for optimality which generalizes the Karush-Kuhn-Tucker conditions (KKT) from scalar optimization (i.e., $\nabla \mathcal{L}(\theta^*) = 0$) to MOO.

(Pareto critical points)

i) A decision vector $\theta^* \in \R^n$ is called Pareto critical if there exists a non-negative vector $\alpha \in \R_{\ge 0}^m$ with $\sum_{i=1}^m \alpha_i= 1$ satisfying

\begin{align*}
\sum_{i=1}^m \alpha_i \nabla \Loss_i(\theta^*) = 0.
\end{align*}

ii) The set $\mathcal{P}_c \coloneqq \left\lbrace \theta \in \R^n \, : \, \theta \text{ is Pareto critical } \right\rbrace \subseteq \R^n$ is called the Pareto critical set.

From a numerical point of view, it is not straightforward to work with Definitions <ref> and <ref>. In this setting, the definition of Pareto critical points is advantageous as it is easier to define algorithms which compute critical points. Computing only critical points is justified by the following observation. Any (weakly) Pareto optimal solution to (<ref>) is Pareto critical, i.e., Definition <ref> formulates a necessary condition for Pareto optimality, as do the standard KKT conditions for $m=1$.

For a general problem (<ref>) the Pareto set $\mathcal{P}$, the weak Pareto set $\mathcal{P}_w$ and the Pareto critical set $\mathcal{P}_c$ consist of infinitely many points. Computing and representing these sets numerically is not possible in a straightforward manner. One could argue that from a practical point of view, it is not necessary to compute the entire Pareto set $\mathcal{P}$, as an accurate finite approximation of $\mathcal{P}$ and $\mathcal{F}$ is sufficient.

§.§ Weighted sum method

Throughout this paper, we mainly focus on the WS method. This is a special instance of so-called scalarization methods, which follow the idea of turning (<ref>) into a family of singleobjective optimization problems. An optimal solution to the singleobjective optimization problem should give a Pareto optimal solution to (<ref>). Define the positive unit simplex $\Delta^m \coloneqq \left \lbrace \alpha \in \R_{\ge 0}^m \, : \, \sum_{i=1}^m \alpha_i = 1 \right \rbrace$. Given a weighting vector $\alpha \in \Delta^m$, we define the problem

\begin{align}
\min_{\theta \in \R^n} \sum_{i=1}^m \alpha_i \Loss_i(\theta).
\label{eq:MOP_weighted_sum}
\tag{WS}
\end{align}

The approach of transferring the problem (<ref>) into the scalar problem (<ref>) is popular. An optimal solution to (<ref>) will give a Pareto critical solution to (<ref>), since any critical point of (<ref>) will be Pareto critical. In fact an even stronger result holds. It can be shown that any (locally) optimal solution $\theta^*$ to (<ref>) is a (locally) weakly Pareto optimal solution for (<ref>) [20].

There is an intuitive relation between (<ref>) and its scalarization (<ref>) which is presented in Figure <ref>. In Figures <ref> - <ref> we visualize the vector $\alpha \in \Delta^m$ in image space $\R^m$ and show the contour of the attainable set $\mathbfcal{L}(\R^n)$ in red and black. Solving (<ref>) is equivalent to pushing the hyperplane orthogonal to $\alpha$ as far as possible in the direction of $-\alpha$. In the well-behaved convex setting in Figure <ref> we end up in a point where the hyperplane is tangential to the set $\mathbfcal{L}(\R^n)$. In the non-convex setting in Figure <ref> problems with this approach arise, as the $\alpha$ does no longer relate to a unique point on the front. As a consequence, we cannot find the non-convex part of the front, which is marked in red. Figure <ref> describes the worst case, where only the two black points can be attained using the WS method.

When we use the WS method to approximate the Pareto set, we first have to choose a fixed number of weights $\left\lbrace \alpha^1, \dots, \alpha^N \right\rbrace \subset \Delta^m$ and then solve the problem (<ref>) for all $\alpha \in \left\lbrace \alpha^1, \dots, \alpha^N \right\rbrace$ to compute optimal solutions $\theta^j \in \argmin_{\theta \in \R^n} \sum_{i=1}^m \alpha^j_i \Loss_i(\theta)$ for $j = 1,\dots, N$. Since any optimal solution to (<ref>) yields a weakly Pareto optimal solution to (<ref>), we know that $\tilde{\mathcal{P}} \coloneqq \left\lbrace \theta^1, \dots, \theta^N\right\rbrace \subset \mathcal{P}_w$. Therefore, $\tilde{\mathcal{P}}$ can be seen as an approximation of the weak Pareto set $\mathcal{P}_w$. A good approximation consists of decision vectors $\theta^j$ which are either evenly distributed in the decision space or in the image space. This is a challenge since for a general (<ref>) it is not clear how to choose the weights $\alpha^j$ a priori to find evenly distributed solutions $\tilde{\mathcal{P}}$. Since we have to solve (<ref>) for every $\alpha^j$, poorly chosen weights result in computational overhead. Additionally, one has to keep in mind how the WS method scales with the number of objective functions. The simplex $\Delta^m$ has $m-1$ dimensions. Therefore, the size of a uniform discretization of the space of weights $\Delta^m$ scales exponentially with the number of objective functions in (<ref>).

§.§ Multiobjective gradient descent algorithm

The multiobjective gradient descent algorithm (MGDA) introduced in [43, 22] is an iterative method that computes from an initial point $\theta_0$ a sequence

\begin{align}
\label{eq:general_iterative_method}
\theta_{k+1} = \theta_k + \eta_k v_k, \,\text{ for }\, k \ge 0,
\end{align}

given step size $\eta_k > 0$ and directions $v_k$. The step size and direction should be chosen in a way such that (<ref>) yields a descent for all objective functions, i.e., $\Loss_i(\theta_{k+1}) < \Loss_i(\theta_k)$ for all $k \ge 0$ and all $i = 1,\dots,m$. In singleobjective optimization, we can use the negative gradient as a step direction and compute a suitable step size by different strategies (e.g., constant step size, backtracking, diminishing step size). In MOO this is not possible since we have $m$ gradients $\nabla \Loss_i(\theta_k)$. Fortunately, there is a way to compute a common descent direction from these gradients using a subroutine. Defining

\begin{align}
\label{eq:mo_steepest_descent_direction}
v_k \coloneqq \argmin_{v \in \R^n} \max_{i=1,\dots,m} \langle \nabla \Loss_i(\theta_k), v \rangle + \frac{1}{2}\lVert v \rVert^2,
\end{align}

there exists a step size $\eta_k > 0$ such that $\Loss_i(\theta_{k+1}) < \Loss_i(\theta_k)$ for all $i = 1,\dots,m$, given the vector $\theta_k$ is not Pareto critical. It can be shown that the multiobjective steepest descent update (<ref>) with step direction (<ref>) and a suitable step size strategy converges to Pareto critical points of (<ref>). Note that problem (<ref>) is computationally not too demanding since it can be reduced to a quadratic optimization problem that is $m$ dimensional (i.e., the number of objectives).

In general, we do not apply plain gradient descent in deep learning, but more sophisticated methods using momentum, weight decay, and stochastic gradient information, most notably the ADAM Algorithm [33]. Similar adaptions as in the singleobjective setting are also possible in MOO. See, e.g., [36] for stochastic MGDA and [57, 56] for methods with momentum.

§.§ Brief overview of alternative methods

Despite the above-mentioned shortcomings, this paper focuses on the WS approach. However, there is a vast body of literature on alternative methods to find solutions to (<ref>) that all have their individual strengths and weaknesses.

In terms of scalarization, there are multiple methods trying to overcome the limitation of WS for non-convex problems. Examples are the $\varepsilon$-contraint method [28, 11], the elastic constraint methods [19], Benson's method [3] and the reference point method [61, 62]. There are further scalarization methods like Chebyshev scalarization [9, 30], Pascoletti-Serafini scalarization [46], the normal boundary intersection method [16] and hybrid methods [41, 60, 13].

Besides scalarization, there are multiple classes of algorithms that aim at directly solving (<ref>). A particularly popular class is evolutionary algorithms. However, they tend to scale poorly for very large parameter dimensions such that they are of little use in the context of deep learning. One of the most popular algorithms of this class is NSGA-II [17].

In addition to the gradient methods mentioned in subsection <ref>, there exist other first- and second-order methods. The resulting algorithms converge to single points on the Pareto set, and they are usually highly efficient. Other first-order methods are conjugate gradient methods [40], conditional gradient methods [2] or trust-region methods [55, 4, 5]. Typical higher-order methods are extensions of the Newton method to MOO, see, e.g.,  [21, 23]. While these methods do only compute a single element of the Pareto set it is possible to adapt them in order to compute multiple solutions [6, 1].

§ MULTIOBJECTIVE OPTIMIZATION AND DEEP LEARNING

At the heart of DL [27, 8] lies a neuron, a computational unit inspired by the human brain. A network of interconnected neurons in a layer-like fashion is called an artificial neural network (ANN). A fundamental property of ANNs is their capacity to approximate any continuous function [15] to arbitrary precision given a sufficiently large network architecture. A neuron processes a $d$-dimensional vector x as input and produces a scalar $a_i(\mathbf{x})$. The operation $a_i(\mathbf{x})$ consists of a linear transformation followed by a nonlinear one known as the activation function,

\begin{equation}
a_i(\mathbf{x}) = \sigma_i \big( \mathbf{w}^{(i)}(\mathbf{x})  + b^{(i)}\big),
\end{equation}

where $\mathbf{w}^{(i)} = \big( w_1^{(i)}, w_2^{(i)}, \cdots w_d^{(i)} \big)$ and $b^{(i)}$ correspond to the weights and the bias of the $i_{th}$ neuron, respectively. Together, we denote the parameters for the $i_{th}$ neuron as $\theta^{(i)} = \big( b^{(i)}, \mathbf{w}^{(i)} \big)$. Deep Neural Networks (DNNs) are constructed by stacking multiple layers of neurons together, which are typically divided into three types: the input layer, the hidden layers, and the output layer as shown in Figure <ref>. Indexing all neurons from all layers with $i$ we denote the model parameters as $\bm{\theta}$, consisting of the concatenation of all $\theta^{(i)}$.

Given a dataset  $D = \{(\mathbf{x}_1, u_1), (\mathbf{x}_2, u_2), \cdots , (\mathbf{x}_\mathrm{N}, u_\mathrm{N})\}$ where $\mathbf{x}_i \in \R^d$ are input vectors and $u_i \in \R $ are the corresponding labels, the goal in supervised learning is to estimate a predictor $\hat{u}(\cdot, \bm{\theta^*}): \R^d \longrightarrow \R$ with learned parameters $\bm{\theta^*}$ such that $\hat{u}(\mathbf{x}_i, \bm{\theta^*})\approx u_i$. The essence of learning entails finding $\bm{\theta^*}$ to match the input-output behavior on the training set $D$ as accurately as possible while ensuring that the learned pattern generalizes to unseen samples. The learning process in DL occurs iteratively, updating the parameters $\boldsymbol{\theta}$ in each step to minimize the discrepancy between the predicted outputs and the actual labels. Mathematically, this translates to minimizing a loss function $\Loss(\hat{u}(\mathbf{x}_i, \boldsymbol{\theta}), u_i)$ which encodes the difference between the model prediction $\hat{u}(\mathbf{x}_i, \bm{\theta})$ and the known label $u_i$. A very common objective function to minimize is the Mean Squared Error.

\begin{equation}
\text{MSE} = \frac{1}{N} \sum_{i=1}^{N} (\hat{u}_i - u_i)^2,
\label{eqn:mse}
\end{equation}

where $\hat{u}_i = \hat{u}(\mathbf{x}_i, \bm{\theta})$. The MSE in equation (<ref>) can be minimized using gradient-based methods, e.g., simple gradient descent [29], stochastic gradient descent [29], etc., where the backpropagation algorithm [51] is used to compute gradients. As we are interested in optimizing multiple objectives using various MOO methods, an interesting case of conflicting objectives arises in the optimization of PINNs, particularly when dealing with noisy data.

Because of the trade-off, instead of a single optimal point, one has to search for the Pareto front as discussed in the previous section. In such scenarios, one needs to resort to MOO techniques. In recent years, there has been an increasing interest in the DL community to employ gradient based MOO techniques  [49], [6], <cit.>, [54]. However, we realised that there are a many pitfalls that one is likely to fall into while using the MOO for DL. For example we found out that the study done in [49] is not only wrong at different levels but also very misleading especially to researchers who are just setting their foot in this field. In [49] weighted sum method has been used to find the Pareto front for PINNs. In this paper, we take a critical approach towards their work and highlight the various mistakes/inconsistencies that we have discovered. Although, we mainly this the study done in [49], our comments are valid for all instances where one uses MOO for DL. In the next section we give a brief introduction to PINNs and frame the MOO problem that we would like to solve.

§.§ Physics Informed Neural Networks

Recently, PINNs [48], [34] have caught a lot of attention in the ML community. This is because PINNs present a versatile and new approach to integrating domain knowledge and data-driven methods in solving complex scientific and engineering problems. PINNs leverage the universal approximation capabilities of ANNs [15] to approximate solutions to differential equations. What sets PINNs apart from traditional neural network applications is the incorporation of physical laws and governing equations in the loss function. This not only enables the neural network to learn from available data but also ensures that its predictions respect the underlying physics of the system that generates the data. With the integration of such physics-based constraints, PINNs can be effective even in data-scarce scenarios, which makes them particularly appealing for applications where data collection is expensive or challenging. Furthermore, the integration of the underlying physical laws acts as a regularization term, which is particularly useful when using noisy data. PINNs have been deployed to solve different kinds of PDEs, e.g., integer-order PDEs [48], fractional PDEs [45], integro-differential equations [64], parametric PDEs [47] or stochastic PDEs [52]. In the following, we give a brief introduction to the approach of PINNs. For a detailed overview, it is recommended to explore the work in [32]. In general, a PDE can be represented as:

\begin{equation}
\begin{split}
\mathcal{N}[u(\mathbf{x}, t)] + u_t &= 0, \hspace{1.2cm}\mathbf{x} \in \Omega, t \in [0, T], \\
u(\mathbf{x}, t_0) &= g_0(\mathbf{x}),  \hspace{.6cm} \mathbf{x} \in \Omega, \\
u(\mathbf{x}, t) &= g_{\partial \Omega}(\mathbf{x}, t), \hspace{.1cm} \mathbf{x} \in \partial\Omega,t \in [0, T],
\end{split}
\label{eqn:generalDifferentialEquation}
\end{equation}

where $\mathcal{N}[\cdot]$ is a nonlinear differential operator and $\Omega \subset \mathbb{R}^d$ is the spatial domain with boundary $\partial \Omega$. The functions $g_0(\mathbf{x})$ and $g_{\partial \Omega}(\mathbf{x})$ describe the initial condition (IC) and the boundary condition (BC), respectively.

Thus, \(L\) encapsulates all the spatial and temporal derivatives of the solution \(u(\textbf{x}, t)\). Specifically, if \(u\) is a function of spatial variables \(x_1, x_2, ..., x_d\) and time \(t\), then the derivatives contained in \(L\) can be formally defined as follows:

\[
L[u(\textbf{x}, t)] = \sum_{|\alpha| \leq m} a_{\alpha}(\textbf{x}, t)D^{\alpha}u(\textbf{x}, t) + \frac{\partial u(\textbf{x}, t)}{\partial t},
\]

where \(\alpha = (\alpha_1, ..., \alpha_n)\) is a multi-index, \(|\alpha| = \alpha_1 + ... + \alpha_n\) denotes the order of the derivative, \(a_{\alpha}(\textbf{x}, t)\) are coefficient functions, and \(D^{\alpha}u(\textbf{x}, t)\) represents the spatial derivatives of \(u\) with respect to \(x_1, x_2, ..., x_n\).

The solution to a PDE is a function (depending on space or, in the time-dependent case, both space and time), and it is thus a straightforward idea to use DNNs to express the PDE solution via a feed-forward DNN that takes as input the space and time coordinates $(\mathbf{x}, t)$ and produces the solution $\hat{u}(\mathbf{x}, t; \boldsymbol{\theta})$. Assuming that such a solution exists, one can define a residue function, $r(\hat{u})$ from equation (<ref>)

\begin{equation*}
r := \mathcal{N}[\hat{u}(\mathbf{x}, t; \boldsymbol{\theta})] + \hat{u}_t,
\label{eqn:residue}
\end{equation*}

which is ideally equal to zero, if we have the the right set of parameters $\boldsymbol{\theta}$ to represent the solution of the respective differential equation. This is achieved by encoding $r$ as a loss function called $\Loss_\mathrm{PDE}$ and using an optimization algorithm to minimize it. The ability to calculate the partial derivatives of the neural network output $\hat{u}$ with respect to the inputs (i.e., position $\mathbf{x}$ and time $t$) through backpropagation facilitates the evaluation of the partial differential operator $\mathcal{N}$ within the neural network setting.
Following the definition of MSE from equation (<ref>), one can write:

\begin{equation}
\Loss_\mathrm{PDE}(\mathbf{x}, t, \boldsymbol{\theta}) = \frac{1}{N_\mathrm{PDE}} \sum_{j = 1} ^{N_\mathrm{PDE}}\lvert r (x_j, t_j) \rvert^2,
\label{eqn:pdeLoss}
\end{equation}

where $N_\mathrm{PDE}$ is the number of data points ($x_j, t_j$) used to calculate $\Loss_\mathrm{PDE}$. The points are called collocation points and correspond to specific locations in the domain where $\Loss_\mathrm{PDE}$ is evaluated. Hence, by choosing the collocation points strategically and minimizing $\Loss_\mathrm{PDE}$, we are enforcing the neural network to satisfy the underlying PDEs at these locations. However, so far we have only encoded the residue $r(\hat{u})$ as a loss function. To obtain a particular solution, we also need to encode the IC and the BC from equation (<ref>) as individual losses to be minimized:

\begin{align*}
\Loss_\mathrm{IC}(\mathbf{x}^\mathrm{IC}, t = t_0, \boldsymbol{\theta}) &= \frac{1}{N_\mathrm{IC}} \sum_{j = 1} ^{N_\mathrm{IC}}\lvert \hat{u}(x_j^\mathrm{IC}, t_0) - g_0(x_j^\mathrm{IC}) \rvert^2 , \\
\Loss_\mathrm{BC}(\mathbf{x}^\mathrm{BC}, t, \boldsymbol{\theta}) &= \frac{1}{N_\mathrm{BC}} \sum_{j = 1} ^{N_\mathrm{BC}}\lvert \hat{u}(x_j^\mathrm{BC}, t_j) - g_{\partial \Omega}(x_j^\mathrm{BC}, t) \rvert^2 ,
\end{align*}

where $\mathbf{x}^\mathrm{IC}$ and $\mathbf{x}^\mathrm{BC} \in \partial \Omega$ are the collocation points for the IC and BC, respectively. $N_\mathrm{IC}$ and $N_\mathrm{BC}$ are the number of points sampled to calculate $\Loss_\mathrm{IC}$ and $\Loss_\mathrm{BC}$, respectively. Combining $\Loss_\mathrm{PDE}$, $\Loss_\mathrm{IC}$ and $\Loss_\mathrm{BC}$ yields the physics loss $\Loss_\mathrm{PHYSICS}$

\begin{align}
\Loss_\mathrm{PHYSICS} = \alpha_\mathrm{PDE}\Loss_\mathrm{PDE} + \alpha_\mathrm{IC}\Loss_\mathrm{IC} + \alpha_\mathrm{BC}\Loss_\mathrm{BC} ,
\label{eqn:physicsLoss}
\end{align}

where $\alpha_\mathrm{PDE}$, $\alpha_\mathrm{IC}$ and $\alpha_\mathrm{BC} \geq 0$ are the weights given to the respective losses.

It should be noted that so far, we have not used any labeled data in defining $\Loss_\mathrm{PDE}$ in equation (<ref>). If we additionally have access to a data set of measurements (e.g., from numerical simulations or experiments), then we can define the data loss, which simply is a standard loss function as we know it from supervised learning.

\begin{equation}
\Loss_\mathrm{DATA} = \frac{1}{N_\mathrm{DATA}} \sum_{j = 1}^{N_\mathrm{DATA}} \lvert \hat{u}(x_j, t_j) - \hat{u}_j \rvert^2.
\label{eqn:dataLoss}
\end{equation}

The total loss function $\Loss(\boldsymbol{\theta})$ is a linear combination of $\Loss_\mathrm{PDE}$ and $\Loss_\mathrm{DATA}$ as follows:

\begin{equation}
\Loss (\boldsymbol{\theta}) = \alpha_\mathrm{DATA}\Loss_\mathrm{DATA} + \alpha_\mathrm{PHYSICS}\Loss_\mathrm{PHYSICS} ,
\label{eqn:totalLoss}
\end{equation}

where $\alpha_\mathrm{DATA}, \alpha_\mathrm{PHYSICS} \geq 0$ are the weights given to the data loss and the physics loss respectively.

By tuning these weights, one can control the interplay between the data and the physics loss and hence affect the trainability and prediction of PINNs. However, choosing the weights is not straightforward and most of the time, either the weights are set to one or they are selected through trial and error to achieve the desired neural network predictions. In either case, it is unclear whether a particular weight is the best choice and why a specific weight might outperform others in terms of the predictive capabilities.

§.§ Multiobjective optimization of PINNs

To understand how the choice of the weights affects the optimization landscape and the independent losses, one needs to rely on the multiobjective treatment of minimizing $\Loss_\mathrm{DATA}$ and $\Loss_\mathrm{PHYSICS}$ as discussed in section <ref>. Hence, the problem of minimizing both losses jointly can be reformulated as an MOP

\begin{align*}
\min_{\boldsymbol{\theta} \in \R^n}\left[ \begin{array}{c}
\Loss_\mathrm{DATA}(\boldsymbol{\theta})  \\
\Loss_\mathrm{PHYSICS}(\boldsymbol{\theta})
\end{array}\right].
\label{eq:MOP_PINNS}
\end{align*}

We will use the WS method (subsection <ref>) and MGDA (subsection <ref>) to solve the MOP and thus find the Pareto front. For the WS method, we denote $\alpha_\mathrm{DATA}$ as $\alpha$ and $\alpha_\mathrm{PHYSICS}$ as $1 -\alpha$ in equation (<ref>). Additionally, we compare our results with NSGA-II [17], an evolutionary algorithm considered a benchmark in MOO.

This will require re-writing Equation <ref> as:

\begin{equation}
\Loss (\boldsymbol{\theta}) = \alpha\Loss_\mathrm{DATA} + (1-\alpha)\Loss_\mathrm{PHYSICS}
\label{eqn:weightedsumloss}
\end{equation}

and minimize $\Loss$ for each instance of $\alpha$ where $\alpha \in [0, 1]$ as discussed in Subsection <ref>.

§.§ Experimental Setup

In the following, we lay out the experimental settings to use the MOO of PINNs for two equations, the logistic equation and the heat equation. We thus start with an ordinary differential equation, as it allows us to visualize some of the key aspects, before moving on to a PDE.

§.§.§ Logistic equation

The logistic equation serves as a mathematical model for population growth over time, incorporating both the population's reproductive capacity and the growth limitations imposed by resource availability, competition, and other environmental factors. This equation finds widespread application in ecology, biology, and other disciplines to investigate population dynamics.

\begin{equation*}
\frac{du(t)}{dt} = r \, u(t)\left(1 - \frac{u(t)}{K}\right),
\label{eqn:logistic}
\end{equation*}

where $r$ is the intrinsic growth rate of the population and $K$ is the carrying capacity, which denotes the maximum population size that the environment can support without exceeding its resources. For simplicity, we will choose $K = 1$, and $u(0) = \frac{1}{2}$ as the IC. For this problem, there exists an exact analytic solution

\begin{align}
u_\mathrm{exact}(t) &= \frac{1}{1+ \exp(-rt)}.
\label{eqn:logisticexact}
\end{align}

§.§.§ Heat Equation

The heat equation serves as a mathematical model for the time-dependent distribution of heat in a given region. In one dimension (e.g., a rod of length $L$, where $\Omega = (0, L)$), the heat equation is written as:

\begin{equation*}
\frac{\partial u(x, t)}{\partial t} =  \kappa \frac{\partial^2 u(x, t)}{\partial x^2},
\label{eqn: heatequation}
\end{equation*}

where $\kappa$ is a positive constant called thermal diffusivity. We choose $u(x, 0) = \sin \left( \pi \frac{x}{L} \right)$ $ \forall$ $ x \in (0, L)$ as an IC, and $u(0, t) = u(L, t) = 0$ $\forall$ $t \in [0, T]$ as Dirichlet BC, respectively. With these conditions, the heat equation is well-posed with the following analytic solution:

\begin{equation}
u_\mathrm{exact} = \sin \left( \pi \frac{x}{L} \right) \cdot e^{\frac{-\kappa \pi^2}{^2}t}.
\label{eqn:heatsolutionexact}
\end{equation}

§.§.§ Experiments

We use equations (<ref>) and (<ref>) to calculate the data loss $\Loss_\mathrm{DATA}$ and the physics loss $\Loss_\mathrm{PHYSICS}$, respectively. The task of minimizing both losses is formulated as an MOP. For this optimization, we employ both the WS and MGDA techniques. The implementation of MGDA is an adaptation from [54] where the dual problem is solved instead of the primal[For a detailed discussion, please refer to [22].] problem. We use the Adam optimizer (with $ \eta=0.003$) to train feed-forward networks with layer sizes $(1,9,9,9,1)$ and $(2,50,50,50,50,1)$, respectively. $\Loss_\mathrm{DATA}$ and $\Loss_\mathrm{PHYSICS}$ are calculated on the same grid of equally spaced collocation points ($20$ points for the logistic equation, and $20\times20$ grid in the space-time domain of the heat equation). The collocation points that fall on the boundaries are used to calculate the loss for the boundary conditions. We corrupt the values taken from the exact solution $u_i$ at the collocation points by adding Gaussian noise with a mean of $0$ and standard deviation $\sigma$ (in most experiments with value $0.1$) to simulate real (noisy) data. In the subsequent chapter, we discuss the major pitfalls of using MOO for DL and demonstrate some aspects using our experiments.

§ PITFALLS

In this chapter, we delve into the outcomes of applying MOO for PINN training, pointing out significant challenges and pitfalls commonly encountered in the process. Our discussion aims to illuminate these challenges, providing a clearer path for those integrating MOO in DL. While our exploration focuses on the Pareto fronts of PINNs, the insights presented apply to any DL problem.

§.§ Identifying the Pareto front

Imagine a scenario where one uses MOO to optimize two objective functions, $\Loss_1$ and $\Loss_2$. Figure <ref> shows four different kinds of plots to highlight the existence of different kinds of Pareto fronts. In all four figures, the horizontal and the vertical axes represent the losses $\Loss_1$ and $\Loss_2$, respectively.  The gray area is the attainable set and the green boundary is the Pareto front. While in Figure <ref>, the Pareto front exhibits a convex shape, in Figure <ref>, the Pareto front has both convex and non-convex regions. Figure <ref> shows the existence of a discontinuous Pareto front. Finally, it is also possible to have fronts that are locally Pareto optimal but are dominated globally. Such a case is described by Figure <ref> where the blue parts are locally Pareto optimal but dominated globally. The globally Pareto optimal regions are again shown by the green parts at the boundary.

Hence, one should always keep in mind that any of the four cases could be present in one's specific application, such that it is likely possible to end up with incomplete Pareto fronts, clustered points, or only local optima. While these issues are easily visualized for low-dimensional problems, ignorance of these challenges can easily lead to poor solutions in high-dimensional problems, as is the case in DL. What adds even more to the challenge is the fact that DNN training is characterized by many local optima. In the worst case, one can even end up with a large fraction of non-optimal solutions (see, e.g., [49] Figure 3). A straightforward fix to mitigate these challenges is to employ a non-dominance test as a post-processing step. Therein, we eliminate all the points that are dominated by at least one other point in the identified Pareto set, meaning that it is superior in all objectives.

§.§ Check if the objectives are conflicting

In section <ref>, we discussed the importance of MOO. However, if the objectives are not conflicting, the Pareto set collapses to a single point, and applying MOO techniques has no additional benefit. To study this issue we have corrupted the solution data with different levels of Gaussian noise. For noiseless data, the two objectives should be non-conflicting, whereas the trade-off should grow for increasing levels of noise.

Logistic equation

Heat equation

Pareto fronts for logistic and heat equation with
varying noise level $\sigma \in \{0.0, 0.1, 0.2, 0.3\}$, indicated by different colors. The total loss is $\Loss(\boldsymbol{\theta}) = \alpha \Loss_\mathrm{DATA} + (1-\alpha)\Loss_\mathrm{PHYSICS}$ with $\alpha \in [0, 1]$.

Figures <ref> and <ref> show the Pareto fronts for the logistic and the heat equation, respectively. The conflict between the objectives is signified by different Pareto fronts corresponding to increasing standard deviations of the Gaussian noise. We observe that an increased noise leads to more trade-offs between the objectives. In the absence of noise, the Pareto front notably collapses to a single point, which is observable in Figure <ref> as the purple points. In [54], MOO is used to solve a multi-task problem. However, the chosen network architecture is so powerful that the conflict between the tasks is essentially resolved, which means that the usage of a method particularly tailored to MOO was not necessary. In a different study [49], where MOO is used to calculate the Pareto fronts for PINNs, no noise is added to the solution data which means that there is no conflict between the data and the physics loss. Hence, there was no need to use MOO. The conclusion of this section is thus that one should think carefully whether the considered objectives are truly conflicting.

§.§ Considering the right scaling

In Figure <ref>, we see that in general, a Pareto front can exhibit a convex or non-convex shape. Although most of the Pareto fronts that we have encountered in our experiments are convex, it is quite possible that one obtains a non-convex front. While running experiments to obtain the Pareto fronts, we noticed that the true shape of the Pareto front can change depending upon the scale that one uses in plotting. For example in Figure <ref>, we have obtained a convex Pareto front for the heat equation. However, the same front appears to be non-convex when plotted on a double logarithmic scale as shown in Figure <ref>. Such a behavior is expected since a non-linear transformation of the objective functions can change the convexity.

Pareto front on a linear scale

Pareto front on a double logarithmic scale

Pareto fronts obtained for heat equation for $L=5, k = 1$ plotted on different scales (a) normal scale, (b) double logarithmic scale

In the study conducted in [49], the Pareto fronts were illustrated using a logarithmic scale, and the analysis revealed both convex and non-convex characteristics. In addition, different interpretations of the convex and the non-convex parts of these fronts are given. In conclusion, the scaling can be highly deceptive for interpretation.

§.§ Know your optimization method

It is important to understand the theory behind the optimization method that one wants to employ. Especially, one should be well aware of what to expect of an optimization method. In Table <ref>, we summarize the strengths and weaknesses of some of the most popular MOO methods. Depending on the MOP at hand, it is important to choose a suitable method. Consider the WS method, for instance, which can only locate the convex regions of a Pareto front (subsection <ref>). This is consistent with the Pareto fronts that we have shown in Figure <ref> which are all convex.

A summary of the strengths and weaknesses of the weighted sum (WS), Weighted Chebyshev scalarization (WCS), multiobjective gradient descent algorithm (MGDA) and evolutionary algorithms (EAs).

Algorithm     Strength     Weakness

* One can use any existing gradient descent algorithm.

* Easy to implement.

* One can steer the direction to specific points by parameter variation.

* Only convex regions can be explored.

* Parametrization to obtain an equidistant distribution of points is very hard.

* Steering towards a specific desired tradeoff is difficult.

* Both convex and non-convex problems can be explored.

* One can steer the direction to specific points by parameter variation.

* Optimization problem with constraints $\Rightarrow$ more challenging to solve.

* Parametrization to obtain an equidistant distribution of points is very hard.

* Both convex and non-convex problems can be explored.

* Often numerically efficient.

* Strong theoretical foundation.

* Parameter free.

* Steering the direction to specific points is not possible.

* Multiple points on the front can be explored only via multi-start in combination with suitable initializations.

* Yields an approximation of the front by a population in a single run.

* Black box treatment of objectives $\Rightarrow$ easy implementation.

* Computationally expensive.

* Scale poorly w.r.t. the number of parameters.

* No convergence guarantees for a finite number of iterations.

In Figure <ref>, we present the evolution of the losses during training, which illustrates the formation of the Pareto front using the WS method. In this plot, for each $\alpha$, we start with the same point in the feasible space, i.e., with the same random initialization of the model, and record the losses during the training. We see that the weighted sum approach can obtain a convex Pareto front. To get a better view of the final Pareto front, we plot a zoomed-in version in Figure <ref> by omitting the first 50 epochs. In Figures <ref> and <ref>, we show the Pareto fronts for the heat equation for different domain sizes $L$ using the WS method. As expected, irrespective of the choice of $L$, the shape of the Pareto front is convex. In contrast, the Pareto fronts shown in [49] have both convex and non-convex fronts. As it is not possible to obtain non-convex fronts using the WS method, the consecutive analysis that the curvature of the Pareto front is related to the domain size appears to be questionable, which emphasizes again that a careful analysis of the optimization routine itself is extremely important.

No epochs omitted.

First 50 epochs omitted.

Loss evolution over training epochs demonstrating the formation of the Pareto front for the logistic equation using the weighted sum method.

$L = 5, k = 1$

$L = 1, k = 1$

Pareto fronts for heat equation with different settings using the weighted sum method.

Gradients not normalized.

Gradients normalized.

Pareto fronts for logistic equation using MGDA.

WS, MGDA, and NSGA-II.

Different network sizes for NSGA-II.

Comparing the Pareto fronts for the logistic equation obtained by (a) WS, MGDA and NSGA-II, (b) increasing size of network architectures for NSGA-II.

Even though we can compute the entire front by varying $\alpha$, the appropriate choice and spacing of the weights is extremely challenging. Through trial and error, we ended at a logarithmically scaled (with a base of 80) set for $\alpha$. Compared to the WS method, MGDA is a parameter-free MOO technique where we do not need to choose such a weight parameter (subsection <ref>). For each random initialization in the feasible space, the MGDA delivers a point on the Pareto front. On the downside, the lack of such a parameter means that we can not control where we land. Furthermore, to mitigate the effect of differences in magnitudes of gradients, it is important to normalize the gradients in each iteration before calculating the descent direction. Figure <ref> shows the Pareto fronts for the logistic equation obtained from MGDA, which are convex, hence validating the empirical results obtained from the WS method. We observe that in Figure <ref> where gradients are normalized, the points on the front are more evenly distributed compared to Figure <ref> where the gradients have not been normalized. We also notice that apart from the non-dominated points that constitute the Pareto front, we also have dominated points. Such points are a consequence of the non-convergence of solutions occurring due to unsuitable initialization.

Figure <ref> shows Pareto fronts obtained using the WS, MGDA, and NSGA-II. While both the WS and MGDA outperform NSGA-II, the WS method finds a large region of the front compared to MGDA. Even though EAs are global optimization methods, the algorithm appears to find a locally optimal solution only, in our case. Furthermore, we observed that compared to gradient-based methods, NSGA-II encountered overfitting more frequently. While we were able to define a stopping criterion (see subsection <ref>) for WS and MGDA, it is not straightforward to apply this to , as we can not keep track of individual points due to recombination. For more details about overfitting in evolutionary algorithms, refer to [25, 24, 26]. In Figure <ref>, we observe that as we increase the size of the DNN, the performance of NSGA-II decreases drastically. This finding supports what other studies have also found [1, 30], i.e., that evolutionary algorithms are less suited to train deep neural networks due to the large number of optimization variables.

Criterion     Scalarization     Gradient-Based     Evolutionary
Number of Parameters     Few to moderate     Many     Few

Number of Objectives     Few (Two)     Few     Many

Optimization Type     Gradient-Based     Gradient-Based     Gradient-Free

Pareto front landing     Controllable     Random     Controllable

Note     By default convex     Many adaptions     Hard to use in DL

Optimization Method Suitability

Criterion     WS     WC     MGDA     Evolutionary
Implementation     Easy     Moderate     Moderate     Easy

Problems suitable for     Convex     Convex/non-convex     Convex/non-convex     Convex/non-convex

Optimization Type     Gradient-Based     Gradient-Based     Gradient-Based     Gradient-Free

Pareto front landing     Controllable     Random     Random     Controllable

Note     By default convex     Many adaptions     Flexible, can be complex     Hard to use in DL

Computational cost     Controllable     Random     Random     Controllable

Note     By default convex     Many adaptions     Flexible, can be complex     Hard to use in DL

Comparison between different MOO techniques

§.§ Convergence

Neural networks require multiple training iterations, or epochs, to gradually refine their approximation of the desired function. Generally, we want to stop training when the model achieves optimal generalization, meaning that further training ceases to produce significant improvements without risking overfitting. To illustrate this concept, we consider the simpler logistic equation problem (<ref>) where we calculate the Pareto front using the WS method. A validation physics loss calculated at collocation points located right between the training points (for physics and data) serves as an indicator for overfitting behavior. Shown in Figure <ref> is the training loss evolution using a static learning rate. We observe that later in training, the losses start to jump unpredictably. We interpret this behavior as the learning rate being too high for that stage of training, indicating that we are close to a minimum. In this case, training for a fixed number of epochs can lead to stopping a model within such a spike, which would render the analysis of the results less meaningful. We also observed that this jumpy behavior was more pronounced and started earlier in training in the models with higher physics loss weight (low $\alpha$), which can be seen in the example given in Figure <ref>. This indicates that a model with high physics loss approaches the minimum faster. To mitigate the problem of the spikes in training loss, we deploy a learning rate scheduler that lowers the learning rate when the loss doesn't improve for a certain number of epochs which leads to the smooth loss evolution seen in Figure <ref>. For the second problem of different convergence rates with different $\alpha$ we deploy the validation loss mentioned above. We stop a model when the validation loss exceeds the (physics-)training loss by more than a fixed threshold ($0.1$) for the last time in training.

§.§.§ Problematic Plateaus

We observe in the smooth evolution (Figure <ref>) several small plateaus in the data loss (blue) and physics loss (orange) plots. Figure <ref> shows the model front after only 400 epochs of training, where the first plateau in training loss appears. Notably, several models in the upper left corner (red points) of the plot are not Pareto optimal (see Definition <ref> $i)$). This means, we cannot rely on a small plateau in the training evolution to indicate convergence.

Low physics loss weight

High physics loss weight

Loss histories over 50,000 epochs on logistic equation with validation loss tracking and constant learning rate

Loss history with learning rate scheduler

Solution front obtained after 400 epochs

Loss history and solution front

§.§.§ Learning Rate Scheduling

To address the jumpy behavior of the losses, we take a learning rate scheduler that lowers the learning rate when the loss doesn't improve for a certain number of epochs. Applying this in the training process results in a more smooth loss evolution, as shown in Figure
, those with higher weight assigned to the data loss,  Thus, even when the Pareto front appears well-formed, convergence may not have been achieved. Revisiting Figure <ref>, we observe that the losses continue to change significantly after epoch 400.

§.§ Overfitting

Overfitting is the phenomenon generally denoting that a model loses its generalization capability despite having a low training loss. This can happen usually when the model is too expressive and not regulated enough, such that it learns to fit individual data points. Generally, this is detected when using a validation set of data points that are not included in training the model. The performance on the validation set compared to the performance on the training set indicates the generalization of the model.

In our case of training with a physics loss and a data loss, overfitting happens when the model loses its generalization of the physics behavior in between training points, to fit the noisy training data as well as the physics directly at the training points. To address this overfitting issue [30]  we track a validation loss during training.

Physics Validation This loss is calculated by assessing the pde loss on validation points that are located right in between the original collocation points that are used for training on the physics and data loss.

In Figure <ref> we observe that the validation loss (green line) initially tracks closely with the physics loss (orange line) in the early stages of training, indicating that the physics behavior is generalized. However, at a certain point, the validation loss deviates abruptly from the physics loss and does not return to it. This divergence implies that the physics behavior of the model at the validation points deviates significantly from its behavior directly at the collocation points, suggesting overfitting. We leverage the difference between the validation loss and physics loss (for training) to detect overfitting while allowing the model to run sufficiently long to achieve convergence.
Figure <ref> showcases the Pareto front obtained after completing 50,000 training epochs. The Pareto front has a broken appearance in the middle and upper left part, showing how an overfitted Pareto front may look like.
Figure <ref> presents the Pareto front generated when applying a threshold rule.

§.§.§ Threshold Rule

This rule is a heuristic we devised to find the optimal stopping condition for training PINNs with noisy data.

Threshold Rule We train the models for the entire 50,000 epochs and then backtrack until the validation loss exceeds the physics loss, up to the defined threshold, for the last time in training.

Cases where the validation loss happens to be lower than the physics loss are ignored, as that does not indicate overfitting. The resulting Pareto front mostly corrects the broken front, but 4 models with higher data loss weight still appear to be incorrectly selected by this threshold rule. The following reasons could account for this.

* Sampling limitations: The points sampled for calculating the validation loss only provide an approximation of overfitting detection. Only if the overfitted portion of the solution contains a point for our validation loss calculation can overfitting be detected. This can be improved by using a denser validation set to increase the likelihood of capturing overfitting.

* Threshold rule selection: The arbitrary fixed-value threshold rule might not be the optimal choice. Alternative approaches could include using a percentage threshold, where the validation loss exceeds the physics loss by at least x%, instead of a fixed numerical value like in our experiment.

50,000 Epochs

Threshold Rule

Pareto fronts obtained for logistic equation using different rules

§.§ Understand model limitations

Although DNNs offer remarkable capabilities, it is crucial to acknowledge their limitations. Excessive expressiveness can lead to overfitting, while insufficient expressiveness can hinder performance. Suppose a model lacks the capacity to accurately model either the underlying physics or the given data. In that case, it may appear to trade-off between the two losses even if they are inherently compatible. We suspect that this could be one of the reasons for having obtained Pareto fronts in the work of Rohrhofer et al., [49] even though there is no real conflict between the objectives.

Another major limitation that DNNs suffer from is the spectral bias <cit.> which corresponds to the bias of ANNs towards learning low-frequency modes. Hence PINNs also struggle to solve high-frequency solutions, which has also been shown in the literature <cit.>. In Figure <ref> are shown Pareto fronts for two different settings which differ in the values of $L$ and $k$ respectively. As we can see in Figure <ref> where $L= 1$ and $k = 25$, we do not have a very well-defined Pareto front. To make sure whether the underlying physics is being accurately represented, we plot the prediction corresponding to the upper right corner blue point of Figure <ref> and the analytic solution of the heat equation with parameters $L=1, k=25$ in Figure <ref>. As can be seen from this plot, the model prediction and the analytic solution are not in agreement with each other. We suspect that the spectral bias of PINNs is the reason for not being able to make a good prediction since $L = 1$ corresponding to a higher frequency solution (check Equation <ref>).

$L = 5, k = 0.04$

$L = 1, k = 25$

Pareto fronts for Heat Equation with different settings

Finally, the performance of DNNs is influenced by several factors, such as the architecture and dimensions of the neural networks, the selection of activation functions, and the choice of optimizers. In the specific context of PINNs, accuracy is intrinsically linked <cit.> to the number of collocation points N, utilized for computing the physics loss defined in Equation <ref>. Consequently, a smaller N may deteriorate the physics loss, potentially leading to erroneous interpretations when represented on a Pareto front. Additionally, as documented in the literature <cit.>, PINNs are known to encounter convergence challenges. Therefore, ensuring the convergence of PINNs to a viable solution, as elaborated in Section <ref>, is of crucial importance in obtaining good solutions.

The accuracy of DNNs depends on many factors e.g., the type and size of the NNs used, activation function, optimizer etc. Furthermore, in case of PINNs, the accuracy is directly dependent on the number of collocation points N [reference] used to calculate the physics loss from Eqn. <ref>. Hence, a smaller N can result in higher physics loss which when plotted on a Pareto front can result to inaccurate conclusions. Last but not least, PINNs also suffer from convergence issues as discussed in the literature <cit.> <cit.> <cit.>. Hence, it is also very crucial to make sure that the PINNs that one is using is converging to a solution as was also discussed in section <ref>.

analytic solution

model prediction

Heat map of heat equation with $L=1, k=25$, of the blue dot (upper right corner)

* Conflict can also arise because of the following reasons:

* Model not being expressive enough (e.g., less parameters) either for the data loss or the physics loss or both.

* Not taking enough collocation points for the Physics loss.

* Spectral issue with PINNs.

* PINNs are not good at predicting higher frequency solutions.

* Hence, the Pareto front can also not be trusted. This is important to know to not make wrong decisions.

* show plots for low frequency along with other three plots.

* If possible, fix this issue/ show instances where we do get good solutions, otherwise provide references to fix the this issues.

§.§ Understand the problem at hand: Avoid jumping to conclusions.

* the goal should not to find the solution which is the closest to the true solution when there is conflict between the  objectives. Finding the solution which is closer to the exact solution only works when we have access to the exact/analytic solution. In the most realistic situation, we are going to have access to some noisy data and some/incomplete physics.

* just because you change the scale (M in heat equation) does not mean that the fronts should look different in nature. There is not reason to assume that and one should look for an alternative explaination in such cases.

Another major flaw which was noticed in the same paper is that a correlation is deducted between the characteristic length scale, M, parametrization coefficient $\kappa$ and the nature of the Pareto front i.e., whether the front is convex or non-convex. In Fig [reference], we show that irrespective of the choice of $M$ and $\kappa$, all the fronts are convex. We believe that that deriving such conclusions is the result of multiple pitfalls that we have discussed above. To be more precise, we summarize them as follows:

(i): Not having a concrete understanding of Pareto Optimality condition [Ref Pareto Optimal condition eqn.]. Hence, even the dominated points have been included in the Pareto front.

(ii): The objectives are not really conflicting as there is no noise added to the exact solution. So the only explanation to have obtained such fronts is because of the random noise and numerical errors.

(iii): Irrespective of how the Pareto Fronts have originated, plotting them on a log-log, scale turns out to be misleading. A Pareto front which is convex in reality may appear as a non-convex front on a log-log plot.

(iv): Finally, having used the weighted-sum method, it is not even possible to find a non-convex front.
Hence, it is very important to make sure that one has not fallen prey to the above mentioned pitfalls and hence drawn wrong conclusions. We believe that the study shown in [49] is wrong at multiple levels.

For example the nature of the fronts i.e., whether the fronts obtained are convex or non-convex has been shown to correlated to the choice of $\kappa$ and $M$. Given that we have explained above that weighted-sum is not supposed to find non-convex fronts, it is simply unreasonable to assume that these fronts arise due to the particular choices of $\kappa$ or $M$. Further, we know for sure that the diffusion equation which we have chosen is well behaved and has a unique solution. The reason to choose such a well posed problem is to make sure that our method is working properly, otherwise there is not need to re-solve equations where analytic solutions are well known. Hence, we recommend that one should not take the Pareto fronts obtained very seriously before investigating them thoroughly.

* Conflict can also arise because of the following reasons:

* Model not being expressive enough (e.g., less parameters) either for the data loss or the physics loss or both.

* Not taking enough collocation points for the Physics loss.

* Include plots demanded in the above text.

* mention the number of collocation points data points.

* sampling techniques if used.

* BC exclusion plot.

* Parameterized BC plot. Talk about that it is a better way to deine the physics loss as one avoids training.

§ CONCLUSION

We have highlighted key difficulties associated with the application of MOO in DL. Our exploration began with the fundamentals of MOO, focusing mostly on the WS method and MGDA. Using both approaches, we illustrate the approximation of the Pareto fronts, navigating the weighting between the data loss and physics loss within PINNs for the modeling of differential equations. We also compared the results obtained with NSGA-II which is an evolutionary algorithm. Although we used MOO for PINNs as an illustration, the challenges and solutions we've discussed are broadly relevant to any DL scenario involving MOO. The pitfalls identified and the strategies proposed for overcoming them can be universally applied to improve the effectiveness of MOO in various DL contexts, making our insights valuable across a wide range of applications.

Our analysis in subsection <ref> reveals how a basic misinterpretation of MOO can lead to erroneous identification of non-dominated points as parts of the Pareto front. In subsection <ref> the importance of understanding the objectives at hand is emphasized. We shine some light on possible misunderstandings in previous works and conclude that using MOO is only fruitful if the underlying objectives are conflicting in nature. In Subsection <ref> we have highlighted the importance of choosing appropriate scales for visualization, showing how a convex Pareto front might misleadingly appear non-convex when plotted on a logarithmic scale, potentially leading to incorrect interpretations. In subsection <ref>, we stress the need for an in-depth understanding of the chosen MOO methodologies. We construct the Pareto fronts via the WS, MGDA, and NSGA-II, illuminating the merits and drawbacks of these methods. Moreover, we offer a strategy for selecting an MOO approach based on the specific problem at hand. Guaranteeing the optimization method's adequate convergence is vital, as discussed in subsection <ref>.

Our research not only clarifies prevalent misunderstandings in the application of MOO within DL, but also lays the groundwork for future explorations, aiming to improve the efficacy and precision of MOO in sophisticated neural network designs. These include adapting more existing MOO techniques for DL and exploring Pareto fronts for more than two objectives, which could reveal more pitfalls yet unknown to our explorations in this paper. For example, one could add sparsity as an additional objective apart from the data loss and the physics loss in PINNs.

An interesting application would be to use the Pareto fronts to select the right hyper-parameter $\alpha$ to scale the data and the physics loss. Although producing a continuous appearing front is expensive especially with he weighted sum, method but given the simplicity of the method, one could use a few values of $\alpha$ to get a sense of the scale of contribution from different losses and hence make a more educated guess rather than choosing the coefficients of different losses by hit and trial.

Although we have used only two objectives in this work, can go beyond two objectives. For example we could add $L_2$ norm as the third objective. In this case, instead of a curve we the Pareto front will look like a surface. For more than three objectives, it visualizing the Pareto front will become more challenging.

§ ACKNOWLEDGMENTS

This research was conducted within the AI junior research group “Multicriteria Machine Learning” which is funded by the German Federal Ministry of Education and Research (BMBF).

[1]
A. C. Amakor, K. Sonntag, and S. Peitz, A multiobjective continuation method to compute the regularization path of deep neural networks, 2024, <https://arxiv.org/abs/2308.12044>.

[2]
P. Assunção, O. Ferreira, and L. Prudente, A generalized conditional gradient method for multiobjective composite optimization problems, Optimization,  (2023), pp. 1–31.

[3]
H. P. Benson, Existence of efficient solutions for vector maximization problems, Journal of Optimization Theory and Applications, 26 (1978), pp. 569–580, <https://api.semanticscholar.org/CorpusID:121336616>.

[4]
M. Berkemeier and S. Peitz, Derivative-free multiobjective trust region descent method using radial basis function surrogate models, Mathematical and Computational Applications, 26 (2021), p. 31.

[5]
M. Berkemeier and S. Peitz, Multi-Objective Trust-Region Filter Method for Nonlinear Constraints using Inexact Gradients,  (2023), <https://arxiv.org/abs/2208.12094>.

[6]
K. Bieker, B. Gebken, and S. Peitz, On the treatment of optimization problems with L1 penalty terms via multiobjective continuation, IEEE Transactions on Pattern Analysis and Machine Intelligence, 44 (2022), pp. 7797–7808, <https://doi.org/10.1109/TPAMI.2021.3114962>.

[7]
R. Bischof and M. Kraus, Multi-objective loss balancing for physics-informed deep learning,  (2021), <https://doi.org/10.13140/RG.2.2.20057.24169>, <http://rgdoi.net/10.13140/RG.2.2.20057.24169>.

[8]
C. M. Bishop and H. Bishop, Deep Learning - Foundations and Concepts, Springer, 2024.

[9]
V. J. Bowman, On the Relationship of the Tchebycheff Norm and the Efficient Frontier of Multiple-Criteria Objectives, in Multiple Criteria Decision Making, H. Thiriez and S. Zionts, eds., Berlin, Heidelberg, 1976, Springer Berlin Heidelberg, pp. 76–86.

[10]
S. Cai, Z. Mao, Z. Wang, M. Yin, and G. E. Karniadakis, Physics-informed neural networks (PINNs) for fluid mechanics: a review, Acta Mechanica Sinica, 37 (2021), pp. 1727–1738, <https://doi.org/10.1007/s10409-021-01148-1>.

[11]
V. Chankong and Y. Y. Haimes, Multiobjective Decision Making: Theory and Methodology, in Multiobjective Decision Making: Theory and Methodology, New York, NY, 1983, Elsevier Science Publishing Co.

[12]
C. A. C. Coello, G. B. Lamont, and D. A. V. Veldhuizen, Evolutionary Algorithms for Solving Multi-Objective Problems, Genetic and Evolutionary Computation, Springer New York, NY, 2 ed., 2007, <https://doi.org/10.1007/978-0-387-36797-2>, <https://doi.org/10.1007/978-0-387-36797-2>.

[13]
H. Corley, A new scalar equivalence for Pareto optimization, IEEE Transactions on Automatic Control, 25 (1980), pp. 829–830.

[14]
G. Coulaud and R. Duvigneau, Physics-Informed Neural Networks for Multiphysics Coupling: Application to Conjugate Heat Transfer, Tech. Report RR-9520, Université Côte d'Azur, Inria, CNRS, LJAD, Oct. 2023, <https://inria.hal.science/hal-04225990>.

[15]
B. C. Csáji et al., Approximation with Artificial Neural Networks, Faculty of Sciences, Etvs Lornd University, Hungary, 24 (2001), p. 7.

[16]
I. Das and J. E. Dennis, Normal-Boundary Intersection: A New Method for Generating the Pareto Surface in Nonlinear Multicriteria Optimization Problems, SIAM Journal on Optimization, 8 (1998), pp. 631–657, <https://doi.org/10.1137/S1052623496307510>.

[17]
K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, A fast and elitist multiobjective genetic algorithm: NSGA-II, IEEE Transactions on Evolutionary Computation, 6 (2002), pp. 182–197.

[18]
M. Ehrgott, Multicriteria Optimization, vol. 491, Springer Science & Business Media, 2005.

[19]
M. Ehrgott and D. M. Ryan, Constructing robust crew schedules with bicriteria optimization, Journal of Multi-Criteria Decision Analysis, 11 (2002), pp. 139–150.

[20]
G. Eichfelder, Adaptive Scalarization Methods in Multiobjective Optimization (Vector Optimization). 2008, Springer, Berlin, Heidelberg, 2008.

[21]
J. Fliege, L. G. Drummond, and B. F. Svaiter, Newton's method for multiobjective optimization, SIAM Journal on Optimization, 20 (2009), pp. 602–626.

[22]
J. Fliege and B. F. Svaiter, Steepest descent methods for multicriteria optimization, Mathematical Methods of Operations Research, 51 (2000), pp. 479–494.

[23]
E. H. Fukuda and L. M. G. Drummond, A survey on multiobjective descent methods, Pesquisa Operacional, 34 (2014), pp. 585–620.

[24]
I. Gonçalves and S. Silva, Balancing learning and overfitting in genetic programming with interleaved sampling of training data, in Genetic Programming, K. Krawiec, A. Moraglio, T. Hu, A. Ş. Etaner-Uyar, and B. Hu, eds., Berlin, Heidelberg, 2013, Springer Berlin Heidelberg, pp. 73–84.

[25]
I. Gonçalves, S. Silva, J. B. Melo, and J. M. B. Carreiras, Random sampling technique for overfitting control in genetic programming, in Genetic Programming, A. Moraglio, S. Silva, K. Krawiec, P. Machado, and C. Cotta, eds., Berlin, Heidelberg, 2012, Springer Berlin Heidelberg, pp. 218–229.

[26]
I. Gonçalves and S. Silva, Experiments on controlling overfitting in genetic programming, 10 2011.

[27]
I. Goodfellow, Y. Bengio, and A. Courville, Deep Learning, MIT Press, 2016.

[28]
Y. Haimes, On a Bicriterion Formulation of the Problems of Integrated System Identification and System Optimization, IEEE transactions on systems, man, and cybernetics,  (1971), pp. 296–297.

[29]
D. M. D. M. Himmelblau, Applied Nonlinear Programming, McGraw-Hill, New York, 1972.

[30]
S. S. Hotegni, M. Berkemeier, and S. Peitz, Multi-objective optimization for sparse deep multi-task learning, 2024, <https://arxiv.org/abs/2308.12243>.

[31]
Y. Jin and B. Sendhoff, Pareto-Based Multiobjective Machine Learning: An Overview and Case Studies, IEEE Transactions on Systems, Man, and Cybernetics, Part C (Applications and Reviews), 38 (2008), pp. 397–415, <https://doi.org/10.1109/TSMCC.2008.919172>.

[32]
G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, Physics-informed machine learning, Nature Reviews Physics, 3 (2021), pp. 422–440, <https://doi.org/10.1038/s42254-021-00314-5>, <http://www.nature.com/articles/s42254-021-00314-5> (accessed 2022-11-24).

[33]
D. P. Kingma and J. Ba, Adam: A Method for Stochastic Optimization, 2017, <https://arxiv.org/abs/1412.6980>.

[34]
I. Lagaris, A. Likas, and D. Fotiadis, Artificial Neural Networks for Solving Ordinary and Partial Differential Equations, IEEE Transactions on Neural Networks, 9 (1998), pp. 987–1000, <https://doi.org/10.1109/72.712178>.

[35]
K. Li, T. Zhang, and R. Wang, Deep Reinforcement Learning for Multiobjective Optimization, IEEE Transactions on Cybernetics, 51 (2021), p. 3103–3114, <https://doi.org/10.1109/tcyb.2020.2977661>, <http://dx.doi.org/10.1109/TCYB.2020.2977661>.

[36]
S. Liu and L. N. Vicente, The stochastic multi-gradient algorithm for multi-objective optimization and its application to supervised machine learning, Annals of Operations Research,  (2021), pp. 1–30.

[37]
B. Lu, C. B. Moya, and G. Lin, NSGA-PINN: A Multi-Objective Optimization Method for Physics-Informed Neural Network Training, Mar. 2023, <http://arxiv.org/abs/2303.02219> (accessed 2023-05-05).
arXiv:2303.02219 [cs].

[38]
Z. Lu, R. Cheng, Y. Jin, K. C. Tan, and K. Deb, Neural architecture search as multiobjective optimization benchmarks: Problem formulation and performance assessment, IEEE Transactions on Evolutionary Computation, 28 (2024), pp. 323–337, <https://doi.org/10.1109/TEVC.2022.3233364>.

[39]
Z. Lu, H. Pu, F. Wang, Z. Hu, and L. Wang, The Expressive Power of Neural Networks: A View from the Width, CoRR, abs/1709.02540 (2017), <http://arxiv.org/abs/1709.02540>, <https://arxiv.org/abs/1709.02540>.

[40]
L. R. Lucambio Pérez and L. F. Prudente, Nonlinear Conjugate Gradient Methods for Vector Optimization, SIAM Journal on Optimization, 28 (2018), pp. 2690–2720.

[41]
K. Miettinen, Nonlinear Multiobjective Optimization, Springer, New York, 1998.

[42]
B. Mitrevski, M. Filipovic, D. Antognini, E. L. Glaude, B. Faltings, and C. Musat, Momentum-based Gradient Methods in Multi-Objective Recommendation, 2021, <https://arxiv.org/abs/2009.04695>.

[43]
H. Mukai, Algorithms for Multicriterion Optimization, IEEE transactions on automatic control, 25 (1980), pp. 177–186.

[44]
A. Navon, A. Shamsian, G. Chechik, and E. Fetaya, Learning the Pareto Front with Hypernetworks, 2021, <https://arxiv.org/abs/2010.04104>.

[45]
G. Pang, L. Lu, and G. E. Karniadakis, fPINNs: Fractional Physics-Informed Neural Networks, SIAM Journal on Scientific Computing, 41 (2019), pp. A2603–A2626, <https://doi.org/10.1137/18M1229845>, <https://doi.org/10.1137/18M1229845>, <https://arxiv.org/abs/https://doi.org/10.1137/18M1229845>.

[46]
A. Pascoletti and P. Serafini, Scalarizing Vector Optimization Problems, Journal of Optimization Theory and Applications, 42 (1984), pp. 499–524.

[47]
M. Penwarden, S. Zhe, A. Narayan, and R. M. Kirby, A Metalearning Approach for Physics-Informed Neural Networks (PINNs): Application to Parameterized PDEs, Journal of Computational Physics, 477 (2023), p. 111912, <https://doi.org/https://doi.org/10.1016/j.jcp.2023.111912>, <https://www.sciencedirect.com/science/article/pii/S0021999123000074>.

[48]
M. Raissi, P. Perdikaris, and G. Karniadakis, Physics-Informed Neural Networks: A Learning Framework for Solving Forward and Inverse Problems Involving Nonlinear Partial Differential Equations, Journal of Computational Physics, 378 (2019), pp. 686–707, <https://doi.org/10.1016/j.jcp.2018.10.045>, <https://linkinghub.elsevier.com/retrieve/pii/S0021999118307125> (accessed 2023-04-30).

[49]
F. M. Rohrhofer, S. Posch, C. Gößnitzer, and B. C. Geiger, Data vs. Physics: The Apparent Pareto Front of Physics-Informed Neural Networks, IEEE Access, 11 (2023), pp. 86252–86261, <https://doi.org/10.1109/ACCESS.2023.3302892>.

[50]
M. Ruchte and J. Grabocka, Scalable Pareto Front Approximation for Deep Multi-Objective Learning, 2021, <https://arxiv.org/abs/2103.13392>.

[51]
D. E. Rumelhart, G. E. Hinton, and R. J. Williams, Learning representations by back-propagating errors, Nature, 323 (1986), pp. 533–536, <https://doi.org/10.1038/323533a0>, <https://doi.org/10.1038/323533a0>.

[52]
C. Salvi and M. Lemercier, Neural Stochastic Partial Differential Equations, CoRR, abs/2110.10249 (2021), <https://arxiv.org/abs/2110.10249>, <https://arxiv.org/abs/2110.10249>.

[53]
J. Sen, ed., Machine Learning - Algorithms, Models and Applications, IntechOpen, dec 2021, <https://doi.org/10.5772/intechopen.94615>, <https://doi.org/10.5772

[54]
O. Sener and V. Koltun, Multi-Task Learning as Multi-Objective Optimization, Advances in neural information processing systems,  (2018).

[55]
M. G. Shaojian Qu and B. Liang, Trust region methods for solving multiobjective optimisation, Optimization Methods and Software, 28 (2013), pp. 796–811.

[56]
K. Sonntag and S. Peitz, Fast Convergence of Inertial Multiobjective Gradient-like Systems with Asymptotic Vanishing Damping,  (2024), <https://arxiv.org/abs/2307.00975>.

[57]
K. Sonntag and S. Peitz, Fast Multiobjective Gradient Methods with Nesterov Acceleration via Inertial Gradient-Like Systems, Journal of Optimization Theory and Applications,  (2024), pp. 1–44.

[58]
S. Vandenhende, S. Georgoulis, W. Van Gansbeke, M. Proesmans, D. Dai, and L. Van Gool, Multi-Task Learning for Dense Prediction Tasks: A Survey, IEEE Transactions on Pattern Analysis and Machine Intelligence,  (2021), p. 1–1, <https://doi.org/10.1109/tpami.2021.3054719>, <http://dx.doi.org/10.1109/TPAMI.2021.3054719>.

[59]
L. Wang, A. H. C. Ng, and K. Deb, eds., Multi-objective Evolutionary Optimisation for Product Design and Manufacturing, Springer London, 1 ed., 2011, <https://doi.org/10.1007/978-0-85729-652-8>, <https://doi.org/10.1007/978-0-85729-652-8>.

[60]
R. E. Wendell and D. Lee, Efficiency in multiple objective optimization problems, Mathematical programming, 12 (1977), pp. 406–414.

[61]
A. P. Wierzbicki, A methodological guide to multiobjective optimization, in Optimization Techniques, K. Iracki, K. Malanowski, and S. Walukiewicz, eds., Berlin, Heidelberg, 1980, Springer Berlin Heidelberg, pp. 99–123.

[62]
A. P. Wierzbicki, The Use of Reference Objectives in Multiobjective Optimization, in Multiple Criteria Decision Making Theory and Application, G. Fandel and T. Gal, eds., Berlin, Heidelberg, 1980, Springer Berlin Heidelberg, pp. 468–486.

[63]
T. Yu, S. Kumar, A. Gupta, S. Levine, K. Hausman, and C. Finn, Gradient Surgery for Multi-Task Learning, CoRR, abs/2001.06782 (2020), <https://arxiv.org/abs/2001.06782>, <https://arxiv.org/abs/2001.06782>.

[64]
L. Yuan, Y.-Q. Ni, X.-Y. Deng, and S. Hao, A-PINN: Auxiliary physics informed neural networks for forward and inverse problems of nonlinear integro-differential equations, Journal of Computational Physics, 462 (2022), p. 111260, <https://doi.org/https://doi.org/10.1016/j.jcp.2022.111260>, <https://www.sciencedirect.com/science/article/pii/S0021999122003229>.