Documents Live, a web authoring and publishing system
If you see this, something is wrong
Table of contents
First published on Tuesday, Jun 23, 2026 and last modified on Tuesday, Jun 23, 2026 by François Chaplais.
Department of Industrial & Systems Engineering, KAIST Email
Department of Industrial & Systems Engineering, KAIST Email
Discipline of Business Analytics, The University of Sydney Email
Department of Mathematical Sciences, Seoul National University Email
Abstract
Physics-informed neural networks (PINNs) are a promising approach for solving partial differential equations (PDEs). Their training, however, is often difficult because multiple loss terms induced by PDE residuals and boundary or initial conditions must be optimized simultaneously. To address this difficulty, existing approaches often construct update directions by explicitly enforcing particular desirable properties, such as scale robustness and simultaneous descent. While effective in many cases, such property-by-property designs can make it unclear which conditions are essential, what geometric principle determines the selected update direction, and how different methods are structurally related. In this work, we formulate update-direction selection for PINN training as a Chebyshev-center problem in the dual cone. The proposed formulation selects a normalized direction that maximizes the minimum distance to the cone facets. The resulting formulation admits an efficient dual problem in a much lower-dimensional space and yields a convergence guarantee in the nonconvex setting. It also recovers the key desirable properties targeted by existing approaches without imposing them separately; rather, they follow from the single geometric criterion underlying the formulation. This makes the selected direction interpretable through a single geometric rule and provides a unified basis for systematically comparing related direction-selection methods. Experiments on several PINN benchmarks further demonstrate strong empirical performance of the proposed method.
1 Introduction
Partial differential equations (PDEs) are fundamental tools for modeling a wide range of phenomena in science and engineering, including fluid dynamics, heat transfer, and wave propagation. Since analytically solving complex PDEs is often intractable, a variety of numerical methods have been developed to approximate their solutions. Among recent learning-based approaches, physics-informed neural networks (PINNs) have emerged as a promising framework for solving PDEs [1, 2]. In PINNs, a neural network is trained to approximate the PDE solution using objectives constructed from the governing equation together with boundary and initial conditions. This flexibility has enabled PINNs to be used for a broad class of forward and inverse problems. At the same time, however, their practical success depends critically on how effectively they can be trained [3, 4].
Training PINNs is challenging in practice, as their performance can be highly sensitive to data sampling, model design, and training dynamics [5, 6, 7]. These challenges have motivated a broad literature on sampling strategies [8, 9, 10], network design [11, 12, 13], and optimization [14, 15]. In this work, we focus on the optimization aspect of PINN training, where multiple loss terms induced by PDE residuals, boundary conditions, and initial conditions must be optimized simultaneously. This multi-loss structure can produce imbalanced gradient magnitudes across objectives, while also creating conflicts among their preferred update directions [16, 17]. As a result, the key optimization challenge is how to construct an update direction from multiple loss terms so that different physical constraints are reflected in a balanced manner.
Existing approaches to multi-loss optimization in PINNs can be broadly categorized into adaptive weighting methods and direction-selection methods. Adaptive weighting methods primarily aim to adjust the relative importance or the scales of different loss terms [7, 14, 15], whereas direction-selection methods seek to construct an update direction from multiple loss-specific gradients [16, 17, 13]. Many such direction-selection methods are motivated by desirable properties of the update direction, such as equalized progress across objectives, simultaneous descent, or robustness to gradient-scale differences. While these constructions are meaningful and often effective, designing directions by imposing such properties one by one can make it less clear which conditions are more essential, what geometric principle determines the selected direction, and how these methods are structurally related. These observations motivate a direction-selection formulation based on a single geometric criterion, rather than on separately imposed properties.
To this end, we formulate the direction-selection problem as a Chebyshev-center problem over the dual cone of the set of loss gradients. The dual cone provides natural candidates for update directions as it collects vectors that are descent directions for every loss function. However, not all directions in the dual cone are equally desirable: a direction close to a facet may give only limited improvement to the corresponding objective. Motivated by this, based on the idea of Chebyshev centers, we select a normalized direction that serves as the center of the largest ball inscribed in the dual cone. Rather surprisingly, the single geometric formulation achieves several desirable properties all at once, which have been the individual targets of previous methods: scale robustness via normalized gradient comparisons, balanced treatment across objectives, and simultaneous descent. Moreover, the geometric interpretation, together with the corresponding primal and dual formulations, allows us to systematically compare related direction-selection methods.
In particular, we derive the dual of the Chebyshev center formulation, defined over a low-dimensional simplex, and show how to recover the corresponding primal update direction from the dual solution. This formulation also leads to a natural connection with Pareto stationarity, which enables a stopping criterion and a convergence analysis for the proposed algorithm. Then we compare our method with some representative methods, including DCGD [16], HARMONIC [13], IMTL-G [18], ConFIG [17], MGDA [19] and GAPO [20], and the comparison shows that the proposed formulation can explain prior cone-based choices, recover desirable properties without imposing them separately, and relax restrictive equality-based formulations. Finally, we validate the proposed method through numerical experiments on a range of PINN problems.
2 Preliminaries
2.1 Physics-informed neural networks
Physics-informed neural networks (PINNs) [1] approximate the solution \(u\) of a partial differential equation by a neural network \(u_\theta\) with parameters \(\theta\). Consider a PDE defined on a domain \(\Omega\), written as \(\mathcal{N}[u](x,t)=0\) for \((x,t)\in\Omega\), together with boundary and initial conditions \(\mathcal{B}[u](x,t)=0\) on \(\partial\Omega\) and \(\mathcal{I}[u](x,0)=0\) on \(\Omega\). The goal is to find \(\theta\) such that \(u_\theta\) approximately satisfies these physical constraints. To this end, PINNs sample collocation points from the interior of the domain, the boundary, and the initial slice, and define separate loss terms for each constraint, e.g., \( \mathcal{L}_{\mathrm{r}}(\theta) = \frac{1}{N_r}\sum_{j=1}^{N_r} \|\mathcal{N}[u_\theta](x_j^r,t_j^r)\|_2^2, \) \( \mathcal{L}_{\mathrm{b}}(\theta) = \frac{1}{N_b}\sum_{j=1}^{N_b} \|\mathcal{B}[u_\theta](x_j^b,t_j^b)\|_2^2, \) and \( \mathcal{L}_{\mathrm{i}}(\theta) = \frac{1}{N_i}\sum_{j=1}^{N_i} \|\mathcal{I}[u_\theta](x_j^i,0)\|_2^2. \) Thus, PINN training naturally involves multiple loss terms associated with different physical constraints, and balancing them is often critical for stable and effective training. Related studies have been conducted to address this issue, which we review next.
2.2 Multi-objective perspectives on training PINNs
To address the difficulty of balancing multiple PINN losses, prior work has explored several directions. One line focuses on adaptive weighting methods for PINNs, such as LRA [7], NTK-based weighting [14], and ReLoBRaLo [15]. Another line builds on classical multicriteria optimization, where a dual cone-based view of simultaneous descent was implicit in the characterization of descent directions [21]. Related to this classical viewpoint, gradient-based multi-objective and multi-task learning methods, including MGDA [19], PCGrad [22], CAGrad [23], and IMTL-G [18], combine or modify gradients to reduce conflict among objectives. More recently, gradient-based multi-objective methods have also been developed for constrained reinforcement learning, preference-based MOO/MORL, and multi-objective LLM alignment, including GCPO [24], PMGDA [25], GAPO [20], and RACO [26]. In the PINN setting, DCGD [16], ConFIG [17], and HARMONIC [13] are closely related methods that construct update directions from multiple loss gradients. Other PINN-specific optimization strategies include the parameter-wise adaptive optimizer MultiAdam [27] and the minimax-based method BGDA [13]. In this paper, we focus on gradient-based methods for determining update directions and compare methods that are closely related to ours in Section 4.
2.3 Chebyshev centers of polyhedral sets
The Chebyshev center is a classical notion for defining the center of a constrained geometric set (see [28]). Given a norm \(\|\cdot\|\), the depth of \(x \in C\) for a bounded set \(C \subset \mathbb{R}^n\) with nonempty interior is \(\inf_{y\in \mathbb{R}^n\setminus C}\|x-y\|\), namely the radius of the largest norm ball centered at \(x\) and contained in \(C\). A Chebyshev center is a point of maximum depth, or equivalently the center of a largest inscribed norm ball. For a polyhedron \(C=\{x\mid a_i^\top x\le b_i, i=1,…,m\}\), requiring a norm ball of radius \(R\) centered at \(x\) to lie in \(C\) is equivalent to \(a_i^\top x+R\|a_i\|_*\le b_i\) for all \(i\), where \(\|\cdot\|_*\) denotes the dual norm of \(\|\cdot\|\). Thus, the Chebyshev center gives a notion of centrality for polytopes. At the same time, it naturally extends to polyhedral cones, explained in Section 3.1.
3 Methodology
3.1 Chebyshev center of the dual cone for direction selection
We consider a multi-objective optimization problem with \(m\) loss functions \(\{\mathcal{L}_i\}_{i=1}^m\). At iteration \(t\), let \(g_i := \nabla \mathcal{L}_i(\theta_t)\) denote the gradient of \( \mathcal{L}_i\) . Under standard gradient descent, the rate of change in \(\mathcal{L}_i\) along a direction \(d\) is determined by \(g_i^\top d\); in particular, if \(g_i^\top d < 0\), \( d\) is not a descent direction for \(\mathcal{L}_i\). Therefore, depending on the choice of \(d\), some losses may decrease while others increase. To avoid such conflicts, we consider directions satisfying \(g_i^\top d \ge 0\) for all \(i=1,\dots,m\). This simultaneous descent condition can be naturally expressed in terms of the dual cone.
Definition 1 (Dual cone)
Let \(\mathbf{K} \subseteq \mathbb{R}^n\) be a cone. Then, the dual cone of \(\mathbf{K}\) is defined as
Let \(\mathbf{K}_t := \operatorname{cone}(g_1,\dots,g_m)\) be the cone generated by the gradients at iteration \(t\), and let \(\mathbf{K}_t^*\) denote its dual cone. If \(d \in \mathbf{K}_t^*\), then \(x^\top d \ge 0\) for all \(x \in \mathbf{K}_t\). In particular, \(g_i^\top d \ge 0\) for all \(i=1,\dots,m\). Thus, choosing the update direction \(d\) from \(\mathbf{K}_t^*\) guarantees a descent direction for all losses. As the dual cone may contain infinitely many feasible directions, an additional criterion is needed to choose an update direction. Inspired by this, in this work, we develop a centrality notion within the dual cone as the selection criterion, favoring directions that lie away from the cone boundaries.
Note that hyperplane \(g_i^\top d = 0\) defines a face of the cone. Since \(g_i^\top d\) determines the rate of change in \( \mathcal{L}_i\) , choosing a direction near this face yields a small rate of decrease for the corresponding loss. This motivates a centrality criterion that favors directions away from the cone boundary. To formalize this, we use the Chebyshev center, a geometric notion that defines a center of a polyhedral set. To elaborate, we consider a pair of conjugate \(\ell_p/\ell_q\) norms with \(1\leq p\leq\infty\) and \(q=p/(p-1)\). Let \(\hat g_i^{(p)} := {g_i}/{\|g_i\|_p}\) denote the gradient normalized with respect to the \(\ell_p\) norm. Throughout the paper, we assume that \(\|g_i\|_p>0\) for all \(i\), so that each normalized gradient is well defined. Since \(\hat g_i^{(p)}\) is a positive rescaling of \(g_i\), the hyperplane \(g_i^\top v = 0\) is equivalently written as \((\hat g_i^{(p)})^\top v=0\). Then (the absolute value of) the distance from \(v\) to this hyperplane is \({g_i^\top v}/{\|g_i\|_p} = (\hat g_i^{(p)})^\top v\). Here, the Chebyshev-center criterion seeks a center \(v\) that maximizes the smallest of the distances to the hyperplanes. Because the dual cone is scale-invariant, this criterion requires normalization of \(v\); otherwise, the distances can grow trivially along the same ray. Since our interest is in the update direction rather than its magnitude, we fix the scale by imposing \(\|v\|_q \le 1\). Note that finding a Chebyshev center boils down to computing the maximum admissible radius of an \( \ell_q\) ball contained in the dual cone. Hence, we deduce the following formulation.
(1)
Note that \((v,r)=(0,0)\) is always feasible for Problem (1), and hence the optimal value satisfies \(r^\star \ge 0\). Therefore, any optimal solution \(v^\star\) satisfies \((\hat g_i^{(p)})^\top v^\star \ge 0\) for all \(i=1,\dots,m\). Since \(\hat g_i^{(p)}\) is a positive rescaling of \(g_i\), it follows that \(g_i^\top v^\star \ge 0\) for all \(i\), and hence \(v^\star \in \mathbf{K}_t^*\).
The radius \(r\) is determined by the closest facets to the center \(v\); hence, only the closest facets are active at the optimum. Other ones may remain strictly farther away, as illustrated in Figure 1. Therefore, the formulation does not enforce that the distances to the hyperplanes are equal, but only requires them to be at least \(r\). This captures the largest admissible inscribed \( \ell_q\) ball centered at the selected direction.
3.2 Dual formulation to get optimal solution efficiently
We next derive the dual formulation of (1) for efficient computation of the update direction.
Proposition 1
The dual of (1) is given by
(2)
Moreover, strong duality holds between (P) and (D).
By Proposition 1, (1) can be solved equivalently through its dual Problem (2). This is advantageous in our setting because the primal variable \(v\) lies in the parameter space \(\mathbb{R}^n\), whose dimension \(n\) is typically very large in neural network training. In contrast, the dual variable \(\alpha\) lies in \(\mathbb{R}^m\), where \(m\) is the number of loss terms and is usually much smaller than \(n\). In such case, solving the dual is computationally more favorable than solving the primal problem directly, since it reduces the optimization to a simplex-constrained convex problem in a much lower-dimensional space.
Although (1) can be solved efficiently through its dual formulation, our main interest is in the primal optimal direction \(v^\star\), which is used as the update direction. Solving the dual problem yields the optimal dual variable \(\alpha^\star\), and strong duality ensures agreement of the optimal values, but it does not directly provide the primal direction itself. The following proposition shows that, when \(r^\star>0\), the primal optimal direction can be recovered explicitly from the dual optimal solution.
Proposition 2
Assume that \(1<p<\infty\). Let \( (v^\star, r^\star) \) be an optimal solution to the primal (1), and let \( \alpha^\star \) be an optimal solution to the dual (2). Define \( w^\star := \sum_{i=1}^m \alpha_i^\star \hat{g}_i^{(p)}\) . If \( r^\star > 0 \), then \(w^\star\neq 0\) and the primal optimal direction is recovered as \( v^\star ={ \operatorname{sgn}(w^\star)\odot |w^\star|^{p-1} }/{ \|w^\star\|_p^{p-1} }\) , where \(\operatorname{sgn}(\cdot)\), \(|\cdot|\), and the power are applied componentwise, and \(\odot\) denotes componentwise multiplication.
Note that when \(p=2\), the recovery formula in Proposition 2 reduces to the simple Euclidean normalization \(v^\star=\frac{w^\star}{\|w^\star\|_2}\). For \(1<p<\infty\), Proposition 2 shows that the primal update direction can be obtained directly from the dual optimal solution. Therefore, in practice, it is sufficient to solve (2), form \(w^\star = \sum_{i=1}^m \alpha_i^\star \hat g_i^{(p)}\), and recover \(v^\star\) using the formula in Proposition 2. This provides an efficient way to compute the update direction without solving the primal problem explicitly.
Remark 1
The cases \(p=1\) and \(p=\infty\) are also well defined for both the primal and dual formulations. However, the closed-form map in Proposition 2 no longer directly recovers the primal direction from the dual solution, as the corresponding primal direction may not be uniquely determined. Related \(\ell_\infty\)-norm constraints also arise in classical multicriteria steepest-descent methods; in that setting, the resulting direction-selection problem can be formulated as a linear program [21]. In Appendix A.3, we spell out the corresponding recovery step explicitly for our formulation.
The same primal–dual framework also provides a natural stopping criterion for the method. To formalize this criterion, we use Pareto stationarity as the first-order notion of optimality in multi-objective optimization.
Definition 2 (Pareto stationarity)
For the multi-objective problem \( \min_{\theta \in \mathbb{R}^n} \; \bigl(\mathcal{L}_1(\theta), \dots, \mathcal{L}_m(\theta)\bigr)\) where each \(\mathcal{L}_i\) is differentiable, a point \(\theta^\star \in \mathbb{R}^n\) is Pareto-stationary if there exists \(\alpha \in \mathbb{R}_+^m\) such that \( \sum_{i=1}^m \alpha_i = 1\) and \( \sum_{i=1}^m \alpha_i \nabla \mathcal{L}_i(\theta^\star) = 0\) .
Proposition 3
Let \( (v^\star, r^\star) \) be an optimal solution to the primal (1) at \(\theta\). If \( r^\star = 0 \), then \(\theta\) is Pareto-stationary.
Proposition 3 shows that the primal optimal value provides a natural stopping criterion for our method. In practice, however, requiring \(r^\star = 0\) exactly is too restrictive, so we use an approximate criterion instead. Specifically, we refer to a point \(\theta_t\) as \(\epsilon\)-Pareto-stationary if the corresponding primal optimal value satisfies \(r^\star_t \le \epsilon\) for a prescribed tolerance \(\epsilon > 0\). By strong duality, this can be checked through \(\|\sum_{i=1}^m \alpha_{t,i}\hat g_i^{(p)}(\theta_t)\|_p,\) which is computed in the algorithm. With this stopping criterion, our method can be summarized in Algorithm 1. Since the dual is a convex problem over a simplex, one may use the Frank–Wolfe algorithm, as for MGDA [19]. At the same time, in the \(\ell_2\) case with \(m\) is small, there exists a faster dual solver; see Appendix C. In the final step, the recovered \(\ell_q\)-unit direction is scaled by the sum of its inner products with all loss gradients to form the final direction.
Algorithm 1 Chebyshev Center-Based Direction Search via the Dual Formulation (DualChebyshev)
Remark 2
For the cases \(p=1\) and \(p=\infty\), the implementation differs from Algorithm 1. The dual objective becomes nonsmooth, so the dual problem is better handled through an LP reformulation rather than the Frank–Wolfe step used above. Moreover, as discussed in Remark 1, these cases require an additional recovery step because the closed-form map in Proposition 2 does not directly apply. Appendix A.5 describes how to solve the dual problem in the \(p=1\) and \(p=\infty\) cases, and Appendix A.3 describes the corresponding recovery step.
3.3 Convergence analysis
The following theorem provides a convergence guarantee for Algorithm 1 in the nonconvex setting under a smoothness assumption on the total loss that is compatible with the chosen \(\ell_q\) geometry. Here, \(\theta^\star\) denotes a reference point satisfying \(\mathcal{L}(\theta^\star) \le \mathcal{L}(\theta_t)\) for all \(t\).
Theorem 1
Assume that the loss functions \(\mathcal{L}_1, \mathcal{L}_2, \dots, \mathcal{L}_m\) are differentiable and the total loss \( \mathcal{L}(\theta) := \sum_{i=1}^m \mathcal{L}_i(\theta)\) is \(\beta\)-smooth with respect to the \( \ell_q\) norm, i.e., \(\mathcal{L}(\theta') \le \mathcal{L}(\theta) + \nabla \mathcal{L}(\theta)^\top(\theta'-\theta) + \frac{\beta}{2}\|\theta'-\theta\|_q^2 \) for all \(\theta,\theta'\). Let \(\{\theta_t\}_{t=0}^{T}\) be the sequence generated by Algorithm 1. Then, for a step size \(\eta = {1}/{\beta}\), the sequence either terminates at an \(\epsilon\)-Pareto-stationary point, or satisfies
This implies that, unless Algorithm 1 terminates earlier at an \(\epsilon\)-Pareto-stationary point, the best iterate among \(\{\theta_t\}_{t=0}^{T-1}\), in terms of the total-loss gradient norm, satisfies an \(O(1/\sqrt{T})\) bound.
4 Comparative analysis
In this section, we compare our method with several closely related direction-selection methods. Recall that ours selects a direction through the Chebyshev-center criterion over the dual cone, rather than by imposing individual properties. We will see that based on the primal and dual formulations for the Chebyshev center, we may systematically analyze and compare the related methods.
Since the related methods considered below are formulated in the Euclidean geometry, we focus in this section on the default case \(p=q=2\). Accordingly, we write \(\hat g_i := \hat g_i^{(2)}={g_i}/{\|g_i\|_2}\) throughout this section. For ease of comparison, we express all methods in this Euclidean notation. Figure 2 gives a geometric overview of the four direction-selection methods compared in this section. The construction details for this illustration are provided in Appendix B.1.
4.1 Geometric comparison with cone-based methods
Cone-based methods use conic geometry to construct update directions from multiple loss gradients. In this subsection, we compare our method with DCGD [16] and HARMONIC [13], focusing on how their constructions relate to our Chebyshev-center formulation.
4.1.1 DCGD
Within the dual cone, three direction-selection variants are considered: Projection, Average, and Center [16]. Among them, the Center variant uses the angle bisector and was reported to achieve the best empirical performance. In the two-objective nondegenerate case, the angle bisector coincides with the direction selected by our Chebyshev-center formulation. Thus, for \(m=2\), the reported advantage of DCGD (Center) is consistent with the centrality criterion used in our method.
For \(m \ge 3\), however, the notion of an angle bisector is no longer canonical, whereas our Chebyshev-center formulation may take any number of objectives.
4.1.2 HARMONIC
HARMONIC [13] constructs the update direction from the harmonized cone \(\mathbf{H}_t:=\mathbf{K}_t\cap \mathbf{K}_t^*\), where \(\mathbf{K}_t\) is the primal gradient cone and \(\mathbf{K}_t^*\) is its dual cone. Candidates are parameterized as \(d=G_t\lambda\), where \(G_t=[g_1,\dots,g_m]\) and \(\lambda\in\mathbb{R}^m_+\), so that \(d\in\mathbf{K}_t\). Then the condition \(d\in\mathbf{K}_t^*\) is imposed as \(G_t^\top d=G_t^\top G_t\lambda\ge0_m\). It then enumerates the extreme rays of the resulting coefficient-space cone via the Double Description method, maps them back through \(G_t\), and averages the normalized rays. Thus, HARMONIC enforces membership in both \(\mathbf{K}_t\) and \(\mathbf{K}_t^*\) by construction.
In contrast, our method neither constructs \(\mathbf{H}_t\) nor enumerates extreme rays. The Chebyshev-center problem (1) directly selects a direction in \(\mathbf{K}_t^*\). Moreover, Proposition 2 gives \(v^\star = w^\star/\|w^\star\|_2\), where \(w^\star\) is a conic combination of the normalized gradients and hence lies in \(\mathbf{K}_t\). Therefore, the \(\mathbf{H}_t\) membership, explicitly enforced by HARMONIC, automatically holds in our method as an optimality consequence of the Chebyshev-center formulation.
4.2 Comparison with equality-based methods based on the primal formulation
We next compare our formulation with methods that determine the update direction through gradient–direction inner products. In particular, IMTL-G [18] and ConFIG [17] both aim to equalize these inner products across objectives, but they differ in whether the resulting direction is guaranteed to lie in the dual cone. Through (1), we view these methods as equality-based counterparts of our inequality-based Chebyshev-center formulation.
4.2.1 IMTL-G
IMTL-G [18] seeks a direction \(d=\sum_{i=1}^m \alpha_i g_i\), with \(\sum_{i=1}^m \alpha_i=1\), whose normalized-gradient inner products are equal across objectives: \(\hat g_1^\top d=\hat g_2^\top d=…=\hat g_m^\top d.\) Under this condition, IMTL-G yields a closed-form expression for the coefficients \(\alpha\). Specifically, letting \(\alpha=[\alpha_2,…,\alpha_m]^\top\), \(U^\top=[\hat g_1-\hat g_2,…,\hat g_1-\hat g_m]\), and \(G^\top=[g_1-g_2,…,g_1-g_m]\), the resulting coefficients are given by \(\alpha=g_1^\top U(G^\top U)^{-1}\).
Thus, IMTL-G explicitly normalizes gradients and imposes equality of their inner products with the update direction. However, as noted in ConFIG [17], the common value can be negative; if \(\hat g_1^\top d=…=\hat g_m^\top d=c<0\), then the resulting direction does not lie in \(\mathbf{K}_t^*\).
In contrast, in our formulation, normalization is not introduced as a separate designed property; it follows from measuring distances to the hyperplane of the dual cone. Moreover, the update direction is selected within \(\mathbf{K}_t^*\), avoiding the negative inner product value described above.
4.2.2 ConFIG
ConFIG [17] normalizes the gradients and sets their inner products with the update direction to the same positive value. In our notation, it solves \(\hat G^\top w=\mathbf{1}_m\) by the pseudo-inverse, \(w=(\hat G^\top)^+\mathbf{1}_m\), where \(\hat G=[\hat g_1,…,\hat g_m]\), and uses the normalized direction \(v=w/\|w\|_2\). A positive scalar rescaling is then applied to form the final update. In fact, one may equivalently formulate ConFIG as the following equality-constrained counterpart of (1):
(3)
By a change of variables, Problem (3) reduces to the pseudo-inverse minimum-norm system used by ConFIG. The detailed derivation is provided in Appendix B.2.
Moreover, the following theorem establishes another relationship between (1) and (3).
Theorem 2
Let \((v^\star,r^\star)\) be an optimal solution of (1), and let \((\tilde v,\tilde r)\) be an optimal solution of its equality-constrained counterpart (3). Assume that (P\( _{=}\) ) is feasible. Then \( r^\star \ge \tilde r\) . Consequently, for every \(i=1,\dots,m\), we have \( \hat g_i^\top v^\star \ge \hat g_i^\top \tilde v\) .
Theorem 2 shows that, for every objective, our method dominates ConFIG with respect to the rate of change along the computed direction. Moreover, the improvement can be strictly larger when some inequality constraints in (1) are inactive at the optimum, as illustrated in the following example.
Example 1
Let \(g_1 = c_1(1,0,0)\), \(g_2 = c_2(0,1,0)\), and \(g_3 = c_3(2,2,1)/3\), where \(c_1,c_2,c_3>0\). Then the normalized gradients are \(\hat g_1 = (1,0,0)\), \(\hat g_2 = (0,1,0)\), and \(\hat g_3 = (2,2,1)/3\). For (3), the optimal direction is \( \tilde v = (1,1,-1)/\sqrt{3}\) , which gives rise to the ConFIG direction. In contrast, for (1), the optimal direction is \( v^\star = (1,1,0)/\sqrt{2}\) . Moreover, the corresponding inner products satisfy \(\hat g_1^\top \tilde v = \hat g_2^\top \tilde v = \hat g_3^\top \tilde v = 1/\sqrt{3}\), whereas \(\hat g_1^\top v^\star = \hat g_2^\top v^\star = 1/\sqrt{2}\) and \(\hat g_3^\top v^\star = 2\sqrt{2}/3\). Therefore, our method yields strictly larger rates of change for all three objectives in this example. This illustrates that enforcing equalities for all inner products, as in ConFIG, can be more restrictive than the inequality-based formulation underlying our method.
We also note that in the two-objective case, the optimum of (1) necessarily activates both constraints, so our method and ConFIG coincide. Overall, ConFIG’s property-driven construction imposes stronger conditions than necessary. In contrast, the Chebyshev-center geometry naturally recovers the intended properties in a relaxed form, where equality is required only for active constraints.
4.3 Comparison with minimum-norm methods based on the dual formulation
We next compare our formulation with minimum-norm convex-combination methods. In particular, we compare with MGDA and GAPO, which construct update directions by minimizing the norm of a convex combination of gradients or normalized gradients.
4.3.1 MGDA
MGDA [19] determines the update direction as the minimum-norm convex combination of the original gradients. In our notation, its formulation can be written as \( \min\{\| \sum_{i=1}^m \alpha_i g_i \|_2: \sum_{i=1}^m \alpha_i = 1, \alpha_i \ge 0, i=1,\dots,m.\}\) . This formulation is closely related to our dual formulation (2). The key difference is that MGDA uses the original gradients \(g_i\), whereas our dual uses the normalized gradients \(\hat g_i = g_i / \|g_i\|_2\). This distinction is particularly important for training PINNs. In PINNs, different loss terms often have substantially different scales, and such imbalance is known to interfere with training [7]. Because MGDA is based on the original gradients, loss-scale imbalance can directly affect the selected direction. In contrast, our dual formulation determines the direction using normalized gradients, thereby removing the influence of gradient magnitudes. Such magnitude-independent direction selection is well suited to PINN training.
4.3.2 GAPO
GAPO [20] normalizes the gradients in the MGDA objective to reduce influence of loss-scale imbalance. In our notation, the formulation can be written as \( \min\{\| \sum_{i=1}^m \alpha_i g_i/\|g_i\|_2^\rho \|_2: \sum_{i=1}^m \alpha_i = 1, \alpha_i \ge 0, i=1,\dots,m.\}\) where \(\rho\) controls the degree of gradient normalization. When \(\rho=1\), the normalized gradients become \(\hat g_i=g_i/\|g_i\|_2\), so GAPO coincides with the dual problem (2) of our method. The difference lies in the derivation and the resulting hyperparameter dependence: GAPO introduces \(\rho\) as a tunable normalization exponent in the MGDA objective, whereas in our formulation the corresponding normalization is determined by the Chebyshev-center geometry and therefore does not depend on such a hyperparameter.
5 Numerical experiments
We conduct numerical experiments to evaluate the proposed method in both forward and inverse PINN settings. We first test its performance on five forward PDE benchmark problems and compare it with a broad set of optimization baselines, including Adam [29], LRA [7], ReLoBRaLo [15], MGDA [19], PCGrad [22], IMTL-G [18], DCGD [16], ConFIG [17], and HARMONIC [13].
We further consider an inverse problem focused on coefficient recovery under different data conditions, comparing the proposed method with the same set of baselines where applicable. We run each experiment over three independent trials with different random seeds and report the mean and standard deviation of the final relative \(L^2\) error. For forward problems, the error is computed on the solution, whereas for inverse problems, it is computed on the recovered coefficient.
5.1 Comparison on forward benchmark equations
We evaluate the \(\ell_2\)-norm formulation of the proposed method on five forward PDE benchmarks solved with vanilla PINNs: two steady-state boundary-value problems (Helmholtz and 2D Kovasznay flow) and three time-dependent initial-boundary-value problems (Burgers’ equation, the 2D varying-coefficient heat equation, and the Klein–Gordon equation). We compare it against the baselines listed above; implementation details are provided in Appendix D.2. Table 1 reports the mean and standard deviation of the final relative \(L^2\) errors. Lower values indicate better accuracy. In each column, we highlight the best and second-best results. A dash denotes an omitted entry: DCGD is omitted for three-term objectives, and ConFIG is omitted for two-term objectives where it coincides with ours.
| Method | Helmholtz | Kovasznay | Burgers’ | Heat 2D VC | Klein–Gordon |
| Adam | 0.0307 (0.0162) | 0.0044 (0.0010) | 0.1659 (0.0358) | 2.9813 (0.6427) | 0.0267 (0.0111) |
| LRA | 1.0000 (0.0000) | 0.0087 (0.0014) | 0.4032 (0.0190) | 0.9999 (0.0000) | 1.0174 (0.0413) |
| ReLoBRaLo | 0.0570 (0.0295) | 0.0048 (0.0004) | 0.1561 (0.0465) | 2.6150 (0.4534) | 0.0389 (0.0230) |
| MGDA | 0.7542 (0.4058) | 0.0095 (0.0014) | 0.1258 (0.0114) | 0.9370 (0.0219) | 0.1494 (0.0478) |
| PCGrad | 0.0075 (0.0014) | 0.0057 (0.0012) | 0.0706 (0.0046) | 0.4185 (0.0829) | 0.0466 (0.0194) |
| IMTL-G | 0.0086 (0.0007) | 0.0039 (0.0007) | 0.2591 (0.0235) | 1.0099 (0.0128) | 0.0871 (0.0950) |
| DCGD | 0.0097 (0.0061) | 0.0052 (0.0025) | — | — | — |
| ConFIG | — | — | 0.0708 (0.0191) | 0.2401 (0.0359) | 0.0339 (0.0068) |
| HARMONIC | 0.0116 (0.0060) | 0.0061 (0.0045) | 0.0677 (0.0065) | 0.2361 (0.0124) | 0.0143 (0.0031) |
| Ours | 0.0072 (0.0033) | 0.0017 (0.0005) | 0.0599 (0.0042) | 0.2316 (0.0190) | 0.0218 (0.0074) |
Overall, these results demonstrate the superior empirical performance of the proposed method on both steady-state and time-dependent forward PDE problems.
5.2 Coefficient recovery in an inverse heat problem
We further evaluate the proposed method on an inverse heat problem, where the goal is to recover the unknown diffusion coefficient \(a(x,y)\) from noisy observations of \(u(x,y,t)\). We consider three data settings: Standard, Higher noise, and Fewer data. We compare against the same baselines listed above except DCGD, with further implementation details provided in Appendix D.3.
Figure 3 compares the coefficient reconstructions obtained by Adam and our method. Adam produces a visibly distorted coefficient field with large structured errors, whereas our method recovers the ground-truth coefficient field with uniformly small reconstruction errors. Table 2 reports the mean and standard deviation of the final relative \(L^2\) errors for the recovered coefficient \(a\). Lower values indicate better coefficient recovery. In each column, we highlight the best and second-best results. Overall, our method shows strong coefficient-recovery performance across the tested data settings.
| Method | Standard | Higher noise | Fewer data |
| Adam | 6.1510 (3.7224) | 6.1409 (3.6977) | 6.0555 (3.6322) |
| LRA | 0.2390 (0.0000) | 0.2390 (0.0000) | 0.2390 (0.0000) |
| ReLoBRaLo | 3.7730 (2.1329) | 3.8501 (2.0031) | 3.0611 (1.0022) |
| MGDA | 0.2389 (0.0001) | 0.2389 (0.0000) | 0.2389 (0.0000) |
| PCGrad | 0.1371 (0.0529) | 0.1272 (0.0494) | 0.1247 (0.0633) |
| IMTL-G | 0.0599 (0.0005) | 0.1069 (0.0085) | 0.0686 (0.0026) |
| ConFIG | 0.0612 (0.0024) | 0.1038 (0.0042) | 0.0655 (0.0034) |
| HARMONIC | 0.0596 (0.0014) | 0.0988 (0.0110) | 0.0694 (0.0008) |
| Ours | 0.0593 (0.0002) | 0.1040 (0.0049) | 0.0631 (0.0026) |
6 Conclusions
We proposed a Chebyshev center-based direction selection framework for multi-loss PINN training. We formulated update direction selection as finding a central direction that maximizes the minimum distance to the dual cone facets under a general \( \ell_p/\ell_q\) geometry. We derived its low-dimensional dual problem, showed how to recover the corresponding primal update direction, and obtained an efficient algorithm with a convergence guarantee in the nonconvex setting. The geometric interpretation and primal–dual formulations also enabled a systematic comparison with representative direction-selection methods, clarifying how our approach relates to existing constructions. Experiments on forward and inverse PINN benchmarks demonstrated the effectiveness of the proposed method.
This work also leaves several directions for future study. Although the formulation is developed for general \( \ell_p/\ell_q\) geometries, our experiments focused on the Euclidean case. A deeper theoretical and empirical understanding of how the choice of \( p\) interacts with PDE structure, loss landscapes, and training dynamics would be valuable. Beyond PINNs, the same dual-cone and Chebyshev-center perspectives may also be useful for other multi-objective learning problems, including multi-task learning and recent applications such as multi-objective LLM alignment.
A Additional detail and proofs for section 3
A.1 Proof of proposition 1
Recall the primal problem
We first derive its dual. Rewriting the constraints as
the Lagrangian is
where \(\alpha_i\ge 0\) for all \(i\) and \(\beta\ge 0\). Rearranging terms gives
Hence the dual function is
For the supremum over \(r\) to be finite, the coefficient of \(r\) must vanish. Therefore,
Otherwise, if \(1-\sum_i\alpha_i>0\), then letting \(r\to +\infty\) yields \(\phi(\alpha,\beta)=+\infty\), and if \(1-\sum_i\alpha_i<0\), then letting \(r\to -\infty\) yields \(\phi(\alpha,\beta)=+\infty\).
Next, we examine when the supremum over \(v\) is finite. Define
Then the part of the Lagrangian involving \(v\) is
By Hölder’s inequality,
Thus,
Therefore, if \(\|w\|_p>\beta\), then there exists \(s\) with \(\|s\|_q=1\) and \(w^\top s>\beta\). Hence, along \(v=\tau s\),
so the supremum over \(v\) is \(+\infty\). Thus, for the dual function to be finite, it is necessary that
Under this condition, we now compute the value of the supremum. If \(\|w\|_p\le \beta\), then
for all \(v\), and equality is attained at \(v=0\). Hence
Therefore, the dual function is
subject to
Minimizing over \(\beta\) yields
and thus the dual problem reduces to
It remains to show strong duality. The primal (1) can be written as the minimization of \(-r\) subject to
and
The objective \(-r\) is linear, each inequality \(r-(\hat g_i^{(p)})^\top v\le 0\) is affine in \((v,r)\), and the constraint \(\|v\|_q\le 1\) defines a convex set. Hence, (1) is a convex optimization problem in the equivalent minimization form. Moreover, Slater’s condition holds, because \((v,r)=(0,-1)\) is strictly feasible:
Therefore, by Slater’s theorem, strong duality holds between (P) and (D).
A.2 Proof of proposition 2
Let \((v^\star,r^\star)\) be an optimal solution of the primal (1), and let \(\alpha^\star\) be an optimal solution of the dual (2). Define
and let \(\beta^\star := \|w^\star\|_p\).
Assume that \(r^\star>0\). Then \(v^\star\neq 0\), since if \(v^\star=0\), feasibility of \((v^\star,r^\star)\) would imply
for all \(i\), and hence \(r^\star\le 0\), a contradiction.
From the proof of Proposition 1, the Lagrangian is
Since \(1<q<\infty\), the \(\ell_q\) norm is differentiable at \(v^\star\neq 0\). Therefore, the KKT stationarity condition with respect to \(v\) gives
Equivalently,
We next show that \(\|v^\star\|_q=1\). By strong duality, the primal and dual optimal values coincide. From Proposition 1, the dual optimal value is \(\beta^\star\), while the primal optimal value is \(r^\star\). Hence
Under the assumption \(r^\star>0\), we therefore have \(\beta^\star>0\) and \(w^\star\neq0\).
Now consider the complementary slackness condition associated with the constraint \(\|v\|_q\le 1\):
Since \(\beta^\star>0\), it follows that
Since the gradient of the \(\ell_q\) norm at \(v^\star\) is
where the operations are applied componentwise, and \(\|v^\star\|_q=1\) by \((2)\), equation \((1)\) becomes
Therefore, componentwise,
Since \(q=p/(p-1)\), we have \(1/(q-1)=p-1\), and hence
Since \(\beta^\star=\|w^\star\|_p\), this gives
This completes the proof.
A.3 The \(p=1\) and \(p=\infty\) cases for the recovery step
The main text focuses on \(1<p<\infty\), where the primal direction is uniquely recovered from the dual solution in closed form, as shown in Proposition 2. The cases \(p=1\) and \(p=\infty\) are also well defined for both the primal and the dual problem. The difference is that the recovery step is no longer given by a closed-form map from the dual solution to a unique primal direction. Indeed, when \(p=1\) or \(p=\infty\), the support direction associated with a dual vector can be non-unique, because the corresponding unit \( \ell_q\) ball has flat faces.
Given an optimal dual solution \(\alpha^\star\), define
By strong duality, \(r^\star\) is the optimal value of the primal problem. Hence, a valid primal optimal direction can be recovered by finding \(v^\star\) such that
and
Any \(v^\star\) satisfying these conditions is primal optimal, since it attains the optimal value \(r^\star\). Thus, in the \(p=1\) and \(p=\infty\) cases, the dual solution provides the optimal radius \(r^\star\), while the corresponding primal direction can be selected from the set of directions that attain this radius.
For \(p=1\), we have \(q=\infty\). In this case, the norm constraint becomes
Therefore, the recovery step can be written as the linear feasibility problem
Equivalently, \(v^\star\) can be obtained directly by solving the primal problem (1) as the following linear program:
For \(p=\infty\), we have \(q=1\). The recovery step is
Introducing auxiliary variables \(z_j\ge 0\), the \(\ell_1\)-constraint can be written as
Thus, the recovery step becomes the linear feasibility problem
Similarly, \(v^\star\) can be obtained directly by solving the primal problem (1) as the following linear program:
These formulations show that the \(p=1\) and \(p=\infty\) cases can be handled through LP-based direction selection. However, unlike the \(1<p<\infty\) case, the primal direction is not recovered from \(w^\star\) by a single closed-form map, and the LP formulation may still have multiple optimal directions. When this happens, the Chebyshev-center criterion alone does not select a unique direction, so an additional selection criterion may be required.
A.4 Proof of proposition 3
Let \((v^\star,r^\star)\) be an optimal solution of the primal (1) at \(\theta\), and assume that
Let \(\alpha^\star\) be an optimal solution of the dual (2). By strong duality from Proposition 1, the primal and dual optimal values coincide. Hence,
which implies
Since \(\hat g_i^{(p)} = g_i/\|g_i\|_p\), where \(g_i=\nabla \mathcal{L}_i(\theta)\), we obtain
Now define
Since \(\alpha_i^\star\ge 0\) and \(\sum_{i=1}^m \alpha_i^\star=1\), we have \(\bar\alpha_i\ge 0\) for all \(i\) and
Moreover,
Therefore, there exists \(\bar\alpha\in\mathbb{R}_+^m\) such that
By Definition 2, \(\theta\) is Pareto-stationary.
A.5 Solving the dual problems for \(p=1\) and \(p=\infty\)
The dual problem (2) remains convex at the cases \(p=1\) and \(p=\infty\), but the objective is nonsmooth. Instead of applying the Frank–Wolfe step in Algorithm 1 directly, these dual problems can be written as linear programs.
For \(p=1\), the dual problem is
Using auxiliary variables \(z_j\ge 0\), this is equivalent to
The optimal value is then
For \(p=\infty\), the dual problem is
Introducing a scalar variable \(\tau\), this is equivalent to
The optimal value is
Thus, in both \(p=1\) and \(p=\infty\) cases, an optimal dual solution \(\alpha^\star\) and the corresponding optimal value \(r^\star\) can be obtained by solving a linear program.
A.6 Proof of theorem 1
For each iteration \(t\), let \(\alpha_t\) be an optimal solution of the dual (2) at \(\theta_t\), and define
Also let
and let \(r_t^\star\) denote the optimal value of the primal (1) at \(\theta_t\).
If Algorithm 1 terminates at some \(t\in\{0,\dots,T-1\}\), then it returns an \(\epsilon\)-Pareto-stationary point by construction. We therefore consider the remaining case where the algorithm does not terminate for any \(t=0,\dots,T-1\). Then
By strong duality, \(\|w_t\|_p=r_t^\star\). Hence
We first establish two auxiliary relations used below.
Step 1: Lower bound on \(\|d_t\|_q\).
Since \(\alpha_t\) is optimal for the dual (2) and \(r_t^\star>0\) by (1), Proposition 2 implies that \(v_t\) is the corresponding optimal primal direction for (P). Therefore,
Multiplying both sides by \(\|g_i(\theta_t)\|_p\) yields
Summing over \(i\), we obtain
where the second inequality follows from the triangle inequality. Since \(\|v_t\|_q=1\) and \(\sum_{i=1}^m g_i(\theta_t)^\top v_t \ge 0\), it follows that
Combining this with \((1)\) gives
Step 2: Relation between \(\nabla \mathcal{L}(\theta_t)\) and \(d_t\).
By definition,
On the other hand, since \(\|v_t\|_q=1\),
Therefore,
Step 3: Descent estimate and final bound via telescoping.
Now use the \(\beta\)-smoothness of \(\mathcal{L}\) with respect to \(\|\cdot\|_q\). Since
we have
Applying \((3)\), this becomes
Now choose \(\eta=1/\beta\). Then
Using \((2)\), we further obtain
Summing \((4)\) over \(t=0,\dots,T-1\) yields
Rearranging,
Since \(\mathcal{L}(\theta^\star)\le \mathcal{L}(\theta_t)\) for all \(t\), in particular \(\mathcal{L}(\theta^\star)\le \mathcal{L}(\theta_T)\). Hence
Dividing by \(T\) and using
we conclude that
Therefore, unless Algorithm 1 terminates earlier at an \(\epsilon\)-Pareto-stationary point, the best iterate among \(\{\theta_t\}_{t=0}^{T-1}\) satisfies the stated bound.
B Additional detail and proof for section 4
B.1 Details of the comparative toy example
We consider a three-objective toy example in \(\mathbb{R}^3\). Let
Then \(\|g_1\|_2=5\), \(\|g_2\|_2=3\), and \(\|g_3\|_2=1\). The normalized gradients are
MGDA solves
For the gradients in this example, the solution is
and hence
Although \(v_{\rm {MGDA}}\in \mathbf{K}^*\), it lies close to the facet of the dual cone. This shows that MGDA is affected by gradient-scale imbalance, causing the selected direction to be biased toward \(g_3\).
For IMTL-G, using its closed-form coefficient rule yields
The resulting aggregated direction before normalization is
After normalization, this gives
Therefore, the normalized inner products are
Thus \(v_{\rm {IMTL}\text{-}G}\notin \mathbf{K}^*\).
For ConFIG, the direction enforces equal normalized inner products,
Solving these equalities and normalizing the direction gives
Indeed,
Therefore,
Our method solves
In this example, the solution is
The normalized inner products are
Hence
The third constraint is inactive because
Therefore,
This illustrates that enforcing equality for all objectives can be more restrictive than the max–min formulation.
B.2 Derivation of the equivalence between ConFIG and the equality-constrained problem
Suppose \(r>0\), and define a new variable \(w = v/r\). Then the equality constraints become
or equivalently,
Moreover, since \(v = rw\) and \(\|v\|_2 \le 1\), we have
which implies
Therefore, maximizing \(r\) in (P\( _{=}\) ) is equivalent to minimizing \(\|w\|_2\) subject to \(\hat G^\top w = \mathbf{1}_m\), that is,
The minimum-norm solution of this system is given by
Hence, the optimal direction of (P\( _{=}\) ) is obtained by normalizing \(w^\star\):
Therefore, the normalized direction used by ConFIG coincides with the optimal direction of the equality-constrained counterpart of our primal (1).
B.3 Proof of theorem 2
Let \(\mathcal{F}\) and \(\mathcal{F}_{=}\) denote the feasible sets of Problems (P) and (P\( _{=}\) ), respectively.
Since every equality implies the corresponding inequality, we have
Therefore, (1) is a relaxation of Problem (P\( _{=}\) ). Since both problems maximize the same objective \(r\), the optimal value over the larger feasible set cannot be smaller than that over the smaller one. Hence,
By feasibility of \((v^\star,r^\star)\) for (P), we have
Combining this with \(r^\star \ge \tilde r\), we obtain
Since \((\tilde v,\tilde r)\) is feasible for Problem (P\( _{=}\) ), its equality constraints imply
Therefore,
This completes the proof.
C Pseudocode for Algorithms
Algorithm 2 Exact Dual Solver for Two Loss Terms
Algorithm 2 follows from the closed-form structure of the two-loss dual problem. In this case, the dual objective reduces to minimizing
Since both gradients are \( \ell_2\) -normalized, this objective is minimized at \( \alpha=1/2\) whenever \( \hat g_1\ne \hat g_2\) . If \( \hat g_1=\hat g_2\) , every \( \alpha\in[0,1]\) gives the same direction, so the symmetric choice \( \alpha=1/2\) remains valid.
Algorithm 3 Exact Dual Solver for Three Loss Terms
Algorithm 3 uses the simple geometric structure of the \( \ell_2\) dual problem over \( \Delta^3\) . Its solution is either a zero-norm convex combination, an interior simplex point, or an edge solution. The first case corresponds to \( 0\in\operatorname{conv}\{\hat g_1,\hat g_2,\hat g_3\}\) . In the second case, the KKT condition implies equal inner products with all three normalized gradients, which yields the candidate \( \tilde\alpha=G^\dagger\mathbf{1}/(\mathbf{1}^\top G^\dagger\mathbf{1})\) . If this candidate is infeasible, the minimizer lies on an edge. Because the gradients are \( \ell_2\) -normalized, the exact edge solution is the midpoint of the pair with the smallest Gram entry \( G_{ij}\) .
Algorithm 4 DualChebyshev-Exact-Adam
Algorithm 4 combines the exact dual solvers above with an Adam-style parameter update. When \( m=2\) or \( m=3\) , the dual weights are computed by Algorithms 2 and 3, respectively. For larger numbers of loss terms, the dual subproblem is solved using a generic Frank–Wolfe procedure. The resulting aggregate direction is then used to recover the primal direction and the final scaled update direction, exactly as in Algorithm 1, before applying Adam moments to the update.
D Experimental details
D.1 Software and hardware environments
We conduct all experiments in a conda environment with Python 3.11, PyTorch 2.7.0+cu128, and CUDA 12.8. Our experiments are run on two servers: a Rocky Linux 9.7 server with 2 Intel Xeon Silver 4210R CPUs and 4 NVIDIA GeForce RTX 3090 24 GB GPUs, and an Ubuntu 24.04.4 LTS server with 2 Intel Xeon Silver 4210R CPUs and 4 NVIDIA GeForce RTX 2080 Ti 11 GB GPUs. Each complete experiment is run entirely on a single server to ensure fair comparison across methods.
D.2 Details for section 5.1
D.2.1 Benchmark equations
The Helmholtz problem is defined by
We use the exact solution
which induces the forcing term
In our experiments, we set \(k=1\), \(a_1=1\), and \(a_2=4\).
The 2D Kovasznay benchmark is governed by the steady incompressible Navier–Stokes equations
where
Dirichlet boundary conditions are imposed from the exact solution
with
We set \(\nu = 1/40\).
The viscous Burgers benchmark is defined by
In our experiments, we use \(\nu = 0.01/\pi\). Since this benchmark does not admit a closed-form exact solution under the above initial and boundary conditions, we evaluate it using a high-resolution numerical reference solution generated offline.
The 2D varying-coefficient heat benchmark is given by
where
Here \(a(x,y)\) denotes the spatially varying diffusion coefficient taken from the PINNacle benchmark dataset [30]. For evaluation, we use the corresponding reference solution provided in the benchmark as a COMSOL-exported solution field.
The Klein–Gordon benchmark is defined by
We use the manufactured exact solution
which yields
D.2.2 Network architecture and training setup
For all forward benchmark problems, we use a fully connected multilayer perceptron with three linear layers, hidden width 50, and hyperbolic tangent activations. All weights are initialized with Xavier normal initialization [31], and all biases are initialized to zero. Each model is trained for 50,000 optimization steps on a single GPU. Following the implementation used throughout experiments, all methods adopt Adam [29] as the base optimizer with fixed parameters \( \beta_1 = 0.9\) , \( \beta_2 = 0.999\) , and \( \epsilon = 10^{-8}\) , and no learning-rate decay is used.
D.2.3 Sampling strategy
For all forward benchmark problems, training points are resampled at every optimization step. Interior collocation points are drawn by uniform random sampling over the corresponding spatial or spatio-temporal domain. For boundary constraints, points are sampled uniformly on each boundary component; for the rectangular 2D domains, the boundary samples are distributed approximately evenly across the four edges, while for the 1D time-dependent problems they are sampled on the left and right spatial boundaries with uniformly sampled time coordinates. For initial conditions, the time coordinate is fixed at \(t=0\) and the spatial coordinates are sampled uniformly over the spatial domain. In the experiments, we use \(N_r=5120\) and \(N_b=256\) for the two steady-state problems, and \(N_r=2560\), \(N_b=256\), and \(N_i=256\) for the three time-dependent problems.
D.2.4 Optimization hyperparameters
To account for differences in optimization behavior across methods, we select the learning rate separately for each method and benchmark from the set \( \{10^{-3},\,3\times 10^{-4}\}\) . Among the two candidates, we use the one that yields the better final relative \( L^2\) error. The selected learning rates are summarized in Table 3. For LRA, we set \( \alpha=0.1\) and update the loss weights every 10 steps. For ReLoBRaLo, we use \( \alpha=0.95\) , \( \beta=0.99\) , \( \tau=1.0\) , and \( \epsilon=10^{-8}\) . For MGDA, we use a Frank–Wolfe subsolver with maximum iteration 20 and tolerance \( 10^{-6}\) . For our method, we set the stopping threshold to \( 10^{-6}\) .
| Method | Helmholtz | Kovasznay | Burgers | Heat 2D VC | Klein–Gordon |
| Adam | \(10^{-3}\) | \(10^{-3}\) | \(10^{-3}\) | \(10^{-3}\) | \(10^{-3}\) |
| LRA | \(3\times 10^{-4}\) | \(10^{-3}\) | \(10^{-3}\) | \(3\times 10^{-4}\) | \(10^{-3}\) |
| ReLoBRaLo | \(10^{-3}\) | \(10^{-3}\) | \(10^{-3}\) | \(10^{-3}\) | \(3\times 10^{-4}\) |
| MGDA | \(10^{-3}\) | \(10^{-3}\) | \(10^{-3}\) | \(10^{-3}\) | \(10^{-3}\) |
| PCGrad | \(10^{-3}\) | \(10^{-3}\) | \(10^{-3}\) | \(3\times 10^{-4}\) | \(3\times 10^{-4}\) |
| IMTL-G | \(10^{-3}\) | \(3\times 10^{-4}\) | \(3\times 10^{-4}\) | \(10^{-3}\) | \(3\times 10^{-4}\) |
| DCGD | \(3\times 10^{-4}\) | \(3\times 10^{-4}\) | – | – | – |
| ConFIG | – | – | \(10^{-3}\) | \(10^{-3}\) | \(3\times 10^{-4}\) |
| HARMONIC | \(3\times 10^{-4}\) | \(3\times 10^{-4}\) | \(10^{-3}\) | \(10^{-3}\) | \(3\times 10^{-4}\) |
| Ours | \(3\times 10^{-4}\) | \(10^{-3}\) | \(10^{-3}\) | \(10^{-3}\) | \(3\times 10^{-4}\) |
D.2.5 Implementation detail for our method
The main text presents our method using a generic dual subproblem solver. In the forward experiments in Section 5.1, which use the \(\ell_2\) setting, we use the exact variant of our method for the two-loss and three-loss cases. In particular, for all reported forward runs of our method, we use the implementation described in Appendix C, which solves the dual problem exactly when \(m=2\) or \(m=3\).
D.2.6 Evaluation metric
We report the relative \( L^2\) error as the evaluation metric. Given predictions \( \hat{u}(x_i)\) and reference values \( u_{\mathrm{ref}}(x_i)\) evaluated at \( N\) points, we compute
For scalar-valued problems, the inner \( \ell_2\) norm reduces to the absolute value. For the Helmholtz equation, 2D Kovasznay flow, and Klein–Gordon equation, \( u_{\mathrm{ref}}\) is given by the corresponding analytic or manufactured solution. For Burgers’ equation, \( u_{\mathrm{ref}}\) is the offline high-resolution numerical reference described above. For the 2D varying-coefficient heat equation, \( u_{\mathrm{ref}}\) is the benchmark reference solution field taken from the PINNacle dataset.
D.3 Details for section 5.2
D.3.1 Inverse problem setup
We consider an inverse diffusion problem in which the solution \(u(x,y,t)\) and the spatially varying diffusion coefficient \(a(x,y)\) are learned simultaneously. The governing equation is
We construct the benchmark from the reference fields
and choose \(f(x,y,t)\) accordingly so that \((u^\star, a^\star)\) satisfies the PDE.
D.3.2 Observation data and data settings
To define the inverse problem, we use noisy pointwise observations of \(u\) together with boundary supervision for \(a\). The observation locations are drawn from a fixed pool of candidate points and then subsampled without replacement. More precisely, given observation points \(\{(x_j,y_j,t_j)\}_{j=1}^{N_{\mathrm{obs}}}\), we construct the observed values as
The resulting noisy observation set is fixed throughout training. In addition, we impose coefficient boundary supervision by matching \(a(x,y)\) to \(a^\star(x,y)\) on \(\partial\Omega \times [0,1]\).
We consider three data settings in Section 5.2: Standard with \(N_{\mathrm{obs}}=2500\) and \(\sigma=0.1\), Higher noise with \(N_{\mathrm{obs}}=2500\) and \(\sigma=0.2\), and Fewer data with \(N_{\mathrm{obs}}=1250\) and \(\sigma=0.1\).
D.3.3 Network architecture and training setup
For the inverse heat problem, we use a parallel fully connected network with two separate branches, one for the solution field \(u(x,y,t)\) and the other for the coefficient field \(a(x,y)\). Each branch is implemented as a multilayer perceptron with five linear layers, hidden width \(50\), and hyperbolic tangent activations. To reflect the fact that the coefficient is time-independent, the \(a\)-branch takes only the spatial coordinates as effective inputs, while the time coordinate is masked out. All weights are initialized with Xavier normal initialization, and all biases are initialized to zero. Each model is trained for 30,000 optimization steps.
D.3.4 Sampling strategy
At each optimization step, the collocation points and the coefficient-boundary points are resampled. The collocation points are drawn by uniform random sampling over \(\Omega \times [0,1]\), while the boundary points for \(a\) are sampled uniformly on the four spatial boundary components of \(\partial\Omega \times [0,1]\), with the samples distributed approximately evenly across the four edges. In contrast, the noisy observation set for \(u\) is generated once before training and then kept fixed throughout optimization. In the final experiments, we use \(N_r=5120\) collocation points and \(N_{b,a}=1024\) coefficient-boundary points per step. The number of observation points is \(N_{\mathrm{obs}}=2500\) for the Standard and Higher-noise settings, and \(N_{\mathrm{obs}}=1250\) for the Fewer-data setting.
D.3.5 Optimization hyperparameters
The optimization hyperparameters are chosen in the same manner as in the forward experiments. In particular, for each method and data setting, we select the learning rate from \(\{10^{-3},\,3\times 10^{-4}\}\) and use the one that yields the better final relative \(L^2\) error. The selected learning rates are summarized in Table 4. Method-specific hyperparameters are the same as those described above for the forward experiments.
| Method | Standard | Higher noise | Fewer data |
| Adam | \(3\times 10^{-4}\) | \(3\times 10^{-4}\) | \(3\times 10^{-4}\) |
| LRA | \(3\times 10^{-4}\) | \(3\times 10^{-4}\) | \(10^{-3}\) |
| ReLoBRaLo | \(10^{-3}\) | \(10^{-3}\) | \(10^{-3}\) |
| MGDA | \(3\times 10^{-4}\) | \(3\times 10^{-4}\) | \(3\times 10^{-4}\) |
| PCGrad | \(3\times 10^{-4}\) | \(3\times 10^{-4}\) | \(3\times 10^{-4}\) |
| IMTL-G | \(10^{-3}\) | \(10^{-3}\) | \(10^{-3}\) |
| ConFIG | \(10^{-3}\) | \(10^{-3}\) | \(10^{-3}\) |
| HARMONIC | \(10^{-3}\) | \(10^{-3}\) | \(10^{-3}\) |
| Ours | \(10^{-3}\) | \(10^{-3}\) | \(10^{-3}\) |
D.3.6 Evaluation metric
As in the forward experiments, we use the relative \(L^2\) error as the evaluation metric. For the inverse problem, however, the reported error is computed on the recovered coefficient field \(a(x,y)\) rather than on the solution field \(u(x,y,t)\).
E Additional experimental results
E.1 Ablation Study on the Norm Parameter
We conduct an ablation study on the norm parameter \(p\) in the proposed \(\ell_p\)-based formulation. The experiment is performed on the Burgers forward benchmark using the same setting as in the forward experiments, with \(2560\) residual points, \(256\) boundary points, and \(256\) initial points. We compare \(p\in\{1.5,2,3,4\}\), and each setting is repeated over three independent runs. Table 5 reports the mean and standard deviation.
Overall, the choices of \(p\) other than \(p=2\) do not improve performance on this benchmark and generally lead to worse results.
| Norm parameter \(p\) | Relative \(L_2\) error |
| \(1.5\) | 0.1317 (0.0522) |
| \(2\) | 0.0599 (0.0042) |
| \(3\) | 0.0966 (0.0238) |
| \(4\) | 0.1628 (0.0069) |
E.2 Computational Cost
To compare the practical computational cost of the competing methods, we record the mean wall-clock time required to complete the full training run for each setting. The forward benchmark problems are run for \(50{,}000\) training steps, whereas the inverse benchmark settings are run for \(30{,}000\) training steps. The measured runtimes, reported in seconds, are summarized in Tables 6 and 7, respectively.
| Method | Helmholtz | Kovasznay | Burgers | Heat 2D VC | Klein–Gordon |
| Adam | 268.07 | 508.71 | 293.77 | 351.97 | 442.28 |
| LRA | 283.39 | 525.18 | 320.22 | 334.14 | 429.72 |
| ReLoBRaLo | 323.73 | 526.12 | 361.85 | 372.15 | 473.33 |
| MGDA | 288.43 | 530.61 | 616.87 | 408.62 | 574.87 |
| PCGrad | 290.90 | 513.03 | 361.78 | 395.85 | 461.49 |
| IMTL-G | 302.55 | 525.45 | 338.28 | 371.99 | 477.12 |
| DCGD | 282.54 | 567.91 | – | – | – |
| ConFIG | – | – | 342.06 | 387.18 | 463.90 |
| HARMONIC | 295.28 | 578.07 | 331.61 | 363.10 | 534.58 |
| Ours | 311.12 | 540.58 | 366.63 | 405.21 | 501.02 |
| Method | Standard | Higher noise | Fewer data |
| Adam | 642.93 | 574.93 | 587.18 |
| LRA | 632.66 | 578.48 | 597.58 |
| ReLoBRaLo | 614.49 | 607.36 | 632.04 |
| MGDA | 649.07 | 676.26 | 669.19 |
| PCGrad | 647.41 | 610.27 | 626.52 |
| IMTL-G | 612.95 | 601.61 | 596.25 |
| ConFIG | 633.47 | 613.84 | 655.04 |
| HARMONIC | 637.30 | 605.30 | 609.60 |
| Ours | 676.37 | 648.24 | 624.63 |
References
[1] Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations Journal of Computational physics 2019 378 686–707
[2] Physics-informed machine learning Nature Reviews Physics 2021 3 6 422–440
[3] Characterizing possible failure modes in physics-informed neural networks Advances in neural information processing systems 2021 34 26548–26560
[4] Scientific machine learning through physics–informed neural networks: Where we are and what’s next Journal of Scientific Computing 2022 92 3 88
[5] A comprehensive study of non-adaptive and residual-based adaptive sampling for physics-informed neural networks Computer Methods in Applied Mechanics and Engineering 2023 403 115671
[6] On the eigenvector bias of Fourier feature networks: From regression to solving multi-scale PDEs with physics-informed neural networks Computer Methods in Applied Mechanics and Engineering 2021 384 113938
[7] Understanding and mitigating gradient flow pathologies in physics-informed neural networks SIAM Journal on Scientific Computing 2021 43 5 A3055–A3081
[8] Efficient training of physics-informed neural networks via importance sampling Computer-Aided Civil and Infrastructure Engineering 2021 36 8 962–977
[9] DAS-PINNs: A deep adaptive sampling method for solving high-dimensional partial differential equations Journal of Computational Physics 2023 476 111868
[10] Failure-informed adaptive sampling for PINNs SIAM Journal on Scientific Computing 2023 45 4 A1971–A1994
[11] Extended physics-informed neural networks (XPINNs): A generalized space-time domain decomposition based deep learning framework for nonlinear partial differential equations Communications in Computational Physics 2020 28 5
[12] Separable physics-informed neural networks Advances in Neural Information Processing Systems 2023 36 23761–23788
[13]
[14] When and why PINNs fail to train: A neural tangent kernel perspective Journal of Computational Physics 2022 449 110768
[15] Multi-objective loss balancing for physics-informed deep learning Computer Methods in Applied Mechanics and Engineering 2025 439 117914
[16] Dual cone gradient descent for training physics-informed neural networks Advances in Neural Information Processing Systems 2024 37 98563–98595
[17] ConFIG: Towards Conflict-free Training of Physics Informed Neural Networks The Thirteenth International Conference on Learning Representations 2025
[18] Towards impartial multi-task learning International conference on learning representations 2021
[19] Multi-task learning as multi-objective optimization Advances in neural information processing systems 2018 31
[20] Gradient-adaptive policy optimization: Towards multi-objective alignment of large language models Proceedings of the 63rd Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers) 2025 11214–11232
[21] Steepest descent methods for multicriteria optimization Mathematical methods of operations research 2000 51 3 479–494
[22] Gradient surgery for multi-task learning Advances in neural information processing systems 2020 33 5824–5836
[23] Conflict-averse gradient descent for multi-task learning Advances in neural information processing systems 2021 34 18878–18890
[24] Gradient-adaptive pareto optimization for constrained reinforcement learning Proceedings of the AAAI Conference on Artificial Intelligence 2023 37 9 11443–11451
[25] Pmgda: A preference-based multiple gradient descent algorithm IEEE Transactions on Emerging Topics in Computational Intelligence 2025
[26] Reward-free Alignment for Conflicting Objectives arXiv preprint arXiv:2602.02495 2026
[27] Multiadam: Parameter-wise scale-invariant optimizer for multiscale training of physics-informed neural networks International conference on machine learning 2023 39702–39721 PMLR
[28] Convex optimization Cambridge university press 2004
[29] Adam: A method for stochastic optimization arXiv preprint arXiv:1412.6980 2014
[30] Pinnacle: A comprehensive benchmark of physics-informed neural networks for solving pdes Advances in Neural Information Processing Systems 2024 37 76721–76774
[31] Understanding the difficulty of training deep feedforward neural networks Proceedings of the thirteenth international conference on artificial intelligence and statistics 2010 249–256 JMLR Workshop and Conference Proceedings