Harmonized Cone for Feasible and Non-conflict Directions in Training Physics-Informed Neural Networks
Authors: Dohyun Bu
Yujung Byun
Jong-Seok Lee
Venue: ICLR 2026 Poster (2026)
Keywords: Physics-Informed Neural Networks
Multi-loss Optimization
Gradient Conflict Resolution
Feasible Directions
Nonconvex Convergence
Primary area: optimization

Abstract: Physics-Informed Neural Networks (PINNs) have emerged as a powerful tool for solving PDEs, yet training is difficult due to a multi-objective loss that couples PDE residuals, initial/boundary conditions, and auxiliary physics terms. Existing remedies often yield infeasible scaling factors or conflicting update directions, resulting in degraded performance. In this paper, we show that training PINNs requires jointly considering feasible scaling factors and a non-conflict direction. Through a geometric analysis of per-loss gradients, we define the $\textit{harmonized cone}$ as the intersection of their primal and dual cones, which characterizes directions that are simultaneously feasible and non-conflicting. Building on this, we propose $HARMONIC$ (HARMONIzed Cone gradient descent), a training procedure that computes updates within the harmonized cone by leveraging the Double Description method to aggregate extreme rays. Theoretically, we establish convergence guarantees in nonconvex settings and prove the existence of a nontrivial harmonized cone. Across standard PDE benchmarks, $HARMONIC$ generally outperforms state-of-the-art methods while ensuring feasible and non-conflict updates.

Published as a conference paper at ICLR 2026
HARMONIZED CONE FOR FEASIBLE AND NON-CONFLICT DIRECTIONS IN TRAINING PHYSICS-INFORME NEURAL NETWORKS
Dohyun Bu' Yujung Byun' Jong-Seok Leef Korea Advanced Institute of Science and Technology fdohyun. bu, yujung4169, jongseok. lee) ac.kr
ABSTRACT
Physics-Infori Informed Neural Networks (PINNs) have emerged as a powerful tool for solving PDEs, yet training is difficult due to multi-objective loss that couples PDE residuals, initial/boundary conditions, and auxiliary physies terms. Exist- ing remedies often yield infeasible scaling factors or conflicting update directions, resulting in degraded performance. In this paper, we show that training PINNs requires jointly considering feasible scaling factors and non-conflict direction. Through a geometric analysis of per-loss gradients, we define the harmonized cone the intersection of their primal and dual cones, which characterizes di- rections that are simultaneously feasible and non-conflicting. Building on this, we propose HARMONIC (HARMONIzed Cone gradient descent), training proce- dure that computes updates within the harmonized cone by leveraging the Double Description method to aggregate extreme rays. Theoretically, we establish con- vergence guarantees in nonconvex settings and prove the existence of a nontriv- ial harmonized cone. Across standard PDE benchmarks, HARMONIC generally outperforms state-of-the oftth---rr- methods while ensuring feasible and non-conflict up- dates.
INTRODUCTION
Partial differential equations (PDEs) provide a foundational framework for modeling phenomena across science and engineering (Morton & Mayers 2005). However, for many PDEs of practical interest, analytic solutions seldom available. In practice. numerical methods such finite differences (Grossmann et al. 2007). finite elements (Bathe) 2006). and spectral methods (Boyd) 2001 approxiapte solutions to the governing equations by discretizing space and time. Such dis- cretizations are prone to error (Papež et al. 2014) and often become unstable on nonlinear or high- dimensional problems (Zhu et al. 2019). Neural networks, with their continuous representations well-suited for automatic differentiation (Baydin et al.) 2018 Yu et al. 2018). offer mesh-free al- ternative for solving PDEs. Building on this capability, Physics-Informed Neural Networks (PINNs)
(Raissi et al. 2019) incorporate the governing equations into the loss function by penalizing PDE residuals while sinultaneously enforcing initial and boundary conditions (Karniadakis al. 202
As PINNs have gained traction (Raissi et al.) 2020: Sahli Costabal et al.) 2020 Wang et al. 2020
Kochkov et al. 2021). research has focused on improving their accuracy through strategics such as objective reformulation (Eshaghi et al.) 2025), adaptive sampling (Wu et al. 2023). and robust initialization (Lee et al., 2025). In parallel, auxiliary terms terms-includde includdig integral constraints YYuan al. 2022) and output- conditione penalties YYu et al., (2022) have been incorporated, giving rise to additional loss components that more faithfully encode complex PDE-driven dynamics.
Even with this strong advantage, incorporating multiple losses inevitably complicates PINN train- ing. PDE residuals and initial/boundary conditions are enforeed simultaneously, and this multi-loss formulation often produces imbalanced gradients during backpropagation-ariss arising from numerical stiffness (Wang et al. 202T), mismatched convergence rates across different terms (Wang et al.) 2022), and incompatible optimal solutions among losses (Lin et al. (211)). To address these chal- lenges, two main streams research have emerged. The first focuses on reweighting, which seeks
"Equal contribution.
1 "oorrsponniig author.

---

Published as a conference paper at ICLR 2026
to minimize non-negative weighted sum of loss terms by adaptively adjusting their relative mag- nitudes (Xiang et al. 2022 Bischof & Kraus, 2025). While this alleviates scale-induced bias, fails to account for scenarios where gradients from different losses conflict in their directions. Since PINNs require all loss components to be driven simultaneously toward zero, this multi-loss formu- lation can also be understood in terms of multi objectives being optimized concurrentyy. From this perspective, the second pursues multi-objective optimization (MOO), which directly enforces non- conflict updates by ensuring nonnegative inner products with all loss-specific gradients (Hwang & Lim 2024 Liu et al. 2025). However, focusing solely on non-conflict often produces an infeasible direction that cannot be realized as the gradient of a reweighted objective with non-negative coef- ficients, as formally defined in Section 3.1] Although several methods have attempted to integrate both perspectives, there is still no consensus on an optimal strategy for combining reweighting and MOO.
To show the necessity of considering both infeasible and conflict, illustrate the bi-objective toy example in Fig. [1 defined by L1 and L over the parameter vector 0 (:1 02). a setting commonly used in the MOO literature (Liu et al. 2021a Senushkin et al. 2023 When updating solely along the conic combination -cc=A1 Ve L A2 Ve L2 with equal weights = the update is dom- inated by the larger-magnitude term and may converge
suboptimal point. Conversely, when updating only along a non- con-fffficct direction Jet satisfying VoL; Id=1 for j=1, the optimiza- tion avoids direct conflicts. Yet converges prematurely to a point outside the Pareto front, since ger can correspond to an infeusible direction. As shown in the right panel of Fig. [ combining the two requiiing the update to lie in (i) the feasible region and (ii) the non-conflict region defined by the loss gradients-guus guides the trajectory toward Pareto-optimal solution shared by both losses.
Given the necessity of avoiding both infeasible and conftict, we propose a gradient-based method using the HARMONIZed Cone gradient descent (HARMONIC). Concretely, the harmo- nized cone is defined as intersection of the primal and dual cones of the per-loss gradients, such that any update direction within this region can be expressed as a conic combination of the loss gradients and simultaneously has nonnegative inner products with all of them.
13-
Parato HARMONIC
Bky 8h.
B
V1
Fly
9
thn
-16
Figure I: Visualization of optimization trajectories with respect to the objective values and the update directions, where each arrow represents the update direction un- der specific loss-gradient conditions, and --- indicates the line orthogonal to the loss gradients. Initialization points are indicated by o and the Pareto front is shown corresponds to gee> which converges suboptimal point. - denotes 9ef. where one trajectory reaches the Pareto front but becomes over-emphasized toward L1, while the other diverges in the opposite direction away from the Pareto front. represents HARMONIC (gx), with all trajectories successfully con- verging to the Pareto front.
This work is the first to provide geomet- ric analysis that explicitly characterizes both the scaling factors of the losses and the directions of the corresponding gradi- ents, offering principled understanding of infeasible and conflict in general multi- loss settings. Furthermore, we establish theoretical convergence of HARMONIC, showing that the proposed intersection re- tains the desirable Pareto-stationary con- vergence behavior. HARMONIC unifies reweighting and non-conflict optimization without extra asymptotic cost. Experiments on the PIN- Nacle benchmark (Zhongkai et al. 2024) show that it generally outperforms Mate-of-the-ar methods within this unified framework.
in
2 RELATED WORKS
Reweighting Strategies. A central line of research aims to alleviate the issue of imbalanced gra- dients through adaptive loss weighting. LRA [Wang et al. 202T rescales individual losses based on gradient statistics, while GradNorm (Chen et al.) 2018) introduces trainable scaling factors that bal- ance the relative rate of loss reduction with respect to their initial magnitudes. SoftAdapt (Heydari et al. 2019) dynamically adjusts weights by measuring the per-loss change between consecutive
2

---

Published as a conference paper at ICLR 2026
Table Comparison of representative approaches against four criteria: applicability to Multiple loss terms (at least three), guarantee of Feasible update gradients with respect to all loss terms, guarantee of Non-Co Con-Conffict update directions with respect to all loss terms, and guarantee of Non- non-Degenerate scaling factors for all loss terms (V = satisfied, X = not satisfied).
Method Class
Multiple
Feasible
Non-Conflict
Non-Degenerate
LRA. NTK, ReLoBRaLo
MGDA
PCGrad, CAGrad, IMTI.-G Aligned-MTT. ConFIG
X
X
DCGD (Center)
X
HARMONIC
iterations. Recently, ReLoBRaLo (Bischof & Kraus 2025 integrates these strategies into a unified framework, achieving robust performance across PDEs. In parallel, leveraging eigenvalue of Neural Tangent Kernel (NTK) (Wang et al. 2022). they dynamically update the loss weights to balance the average convergence rates among losses.
MOO Perspectives. Recent advances in multi-task learning (MTL) recast gradient-based up- dates through the lens of multi-objective optimization. MGDA (Sener & Koltun 2018) seeks the minimum-norn vector inside the convex hull of task gradients. PCGrad (Yu et al. (2020] reduces conflicts by projecting pairwise gradients before aggregation. CAGrad (Liu et al. 202 la) maximizes the minimum loss decrease across tasks, bounded within Euclidean ball around the mean gradi- ent. IMTL-G (Liu et al. 2021b) aims to balance task gradients with respect to both their directions and magnitudes. Aligned-MTL (Senushkin et al. (023) enforces spectral lower bound gra- dient alignment. These perspectives have also been extended to PINNs (Zhou et al. 2023), which is made possible by the fact that each loss term in PINNs is required to approach zero simultane- ously. Consequently, this motivates the use of MOO-base based gradient strategies the PINN setting. ConFIG (Liu al. 2025) constructs a pseudo-inverse-be basee update that guarantees strictly positive inner products with all gradients and equal projection lengths across losses. DCGD (Hwang & Lim 2024 introduces a geometric analysis of the mean gradient, aiming to balance not only directional conflicts but also magnitude disparities
Limitations. Table summarizes the limitations of existing approaches. Reweighting strategies do not explicitly account for gradient conflicts, and thus may result in conflicting update directions. MGDA guarantees confoonfiic updates, but it does not enforce non non-degenerate scaling factors for all loss terms, which can result in some gradients being entirely suppressed. PCGrad, CAGrad, and IMTL.-G attempt to mitigate such sacrifices by heuristically modifying gradient directions, yet they provide no guarantees. Aligned-MTTL and ConFIG further ensure non-conflict for all loss terms, but their updates may fall outside the feasible region. DCGD satisfies the last three criteria but is dif- ficult to extend beyond two losses. To the best of our knowledge, prior work has not simultaneously addressed all four desiderata.
PRELIMINARIES
3.1 PRIMAL AND DUAL GRADIENT CONE
Notation. Let m, denote the number of loss functions and the dimension of the parameter vector 0, respectively, with d > m. For each j=l,..., define the loss gradient VoL E Rdd Collecting these, we form the gradient matrix G:-(g...... .m.] E dx which we assume to be of full rank. Define matrix D:=lg-1.-.. ...0............. g-m] E 12dxm as the Moore-Penrose pseudo-inverse of G and let A(G) E Rd denote the update direction. Additionally, O, and 1m denote the n-diiensional vectors whose entries are all and all 1, respectively.
Primal Gradient Cone Region. A primal gradient cone IK is generated by G (Vertex repre- sentation, V-representation) and simultaneously represented by D (Half-space representation, H- representation).
Definition 1 (Primal Gradient Cone). We define the primal gradient cone
K : (G> X A (ET))
3

---

Published as a conference paper at ICLR 2026
equivalently
k: (a E R D' mm
where RR denotes the m-dimensional nonnegative orthant.
Ej (L, )يم,)) for weights satisfying A 2 for all J. which defines the feasible region. Conversely, A(G) lies outside K, some weights become negative, leading infeasible region. The detailed proof is provided in the Appendix A.I
When A(G) lies within K, it can be written as A(G) Vo
Dual Gradient Cone Region. dual gradient cone defined as the set of vectors whose inner products with all 9j are nonnegative.
Definition 2 (Dual Gradient Cone). The dual gradient cone of lK is defined
K E Rd
G y 0m)
or equivalently
K (Dw RT).
A(G) inside K', A(G)) ensures that none of the loss terms increase, defining the non-conflict region. Conversely, when A(G) falls outside K* some loss inevitably increases, leading to conflict region. The detailed proof is provided the Appendix A.2
3.2 EMPIRICAL OBSERVATIONS AND ISSUES IN TRAINING PINNS
Limitations of Single-Aspect Strategies. As noted in Section 2 prior approaches have mainly relied on either the reweighting or the non-conflict strategies. Using the one-dimensional Burg- ers equation, Fig. - shows training curves and test MSE, highlighting the drawbacks of con- sidering only one of these aspects. Fig. 2a shows that when reweighting update direction (G): 1 プ g falls outside the dual cone K", the test MSE starts to increase. This indicates that conftict arises, causing one of the training losses to increase as the update is aligned opposite to its corresponding loss gradient. In contrast, Fig. 2b demonstrates that although all training losses decrease, the test MSE rises whenever non-conflict update direction (G), with GT (G) leaves the primal cone K. This indicates that even when conflict is mitigated, infeasible directions lead to degraded test performance.
3,3
5
580
1200
200
(a) Reweighting Method
(b) Non-Conflict Method
Figure 2: Training curves for the initial loss Li, boundary loss Lan PDE residual loss Ly. and
test MSE for one-dimensional Burgers equation. trajectory, -- denotes the Ly trajectory, indicates the test MSE when the update follows a feasible and non-conflict direction, indicates the test MSE when the update follows a conflict direction, and indicates the test MSE when the update follows an infeasible direction.
denotes the L; trajectory,
denotes the L
Infeasible and Conflict Directions Exposed by Cosine Values. To further investigate these chal- lenges in training PINNs, we record the cosine values between A(G') and each g_, (Fig. 3a
3b
as well as with 9, (Fig. 3e 3d). for all j throughout training. As shown in Fig. 3al and Fig 30
under the reweighting method, the cosine values between the PDE loss gradient and A.o frequently close to 1, indicating that the updates are often biased toward the PDE loss. In particular, Fig, 3e shows that more than half of the cosine values between the initial loss gradient and Ace (G) are negative, suggesting that, from the perspective of the initial loss, the updates are often made in deteriorating direction. For the non-conflict method, Fig. 3d shows that the cosine values between
4

---

Published as a conference paper at ICLR 2026
(G) and both the PDE loss gradient and the initial loss gradient mostly close to indicating relatively less improvement compared to the boundary loss. In Fig. 3b the cosine values between A, (G) and the extreme rays of IK' corresponding the initial and PDE loss gradients are negative, indicating that (Gf(G) moves in an infeaxible direction.
METHODOLOGY
In this section, define harmonized cone, the region that is simultaneously feasible and non- conflict, and present a geometric analysis considering both the scaling factors and directions of all loss gradients. Then, we introduce algorithm, which ensures that A(G)) remains inside this region.
4.1 HARMONIZED CONE REGION
A harmonized cone IHI is defined as the set of vectors that have nonnegative inner product values with both 9j and 9 for all j. Definition (Harmonized Cone). The set
H {= E IRd
G 2 0m, DJ= 2 0m)
or equivalently
: Rox a E
is called the harmonized cone of G, where R E IR denotes the generating matrix l.. H...
When A(G) lies within IHI, it is both feasible and non-conflict. Conversely, if A(G) & HI, then it either infeasible, which may lead to degraded test performance, or one the losses inevitably increases even with sufficiently small step sizes due to the conflict. The detailed proof is provided in the Appendix B.I
Theorem ] establishes the necessary and sufficient conditions for A(G) to remain inside H., ex- pressed terms of the relative magnitudes between the gradients of all loss terms. addition, Theorem provides the same characterization from dual perspective.
Theorem 1. Suppose that the gradient matrix at iteration is denoted by G(O, and let (H()) be the harmonized cone induced by G Define = ..... TT E ET as scaling vector, where each component X, is the factor corresponding to the loss Lj. Then, a feasible update direction G (t) E Hi(() if and only if
G (1))
0m
Theorem Suppose that the iteration is denoted by D(0,, and let H() be the harmonized cone induced by D Define (=()) (t) WT E Tem scaling vector, where each component wj is the scaling fac- tor corresponding to g-g. Then, a non-conftict update direction (( e H()) if and only if
pseudo-inverse of transposed gradient matrix at
D (t)) W (t) 2 0 0m-
Theorems I and establish clear criteria for when reweighting or non-conflict can lead to adverse training in PINNs, with detailed proofs provided in the Appendix B.2 B.3 Together, these results highlight the harmonized cone as a unifying concept that integrates reweighting and non-conflict approaches.
10
(a) Reweighting (in K))
(b) Non-Conflict (in K)
(c) Reweighting (in K*)
(d) Non-Conflict (in K*)
Figure 3: Histograms of cosine values (o,) between A(G) and 95,, and between A(G) and gj-
5

---

Published as a conference paper at ICLR 2026
Algorithm HARMONIzed Cone gradient descent (HARMONIC)
Require: G E pdxm .7()) (learning rate), ((t) (network parameters).
Im
1: Form H-representation (Theorem 1: A G
2: Double Description (H-V): compute the Verepresentation of {
Initialize (0) fei...... .C....... Cm }. where e; is the j-th standard basis vector.
For each constraint row af of G G /||/
Partition R(t-1) into (R R2 R
Update R(() + R U {(a n
π - (a π ) x RØ, क R
Remove duplicates to obtain II [
T (1) Rmxp.
3: Map rays to parameter space: for each
j= ................... p. set (1) t m/ E ER...
4: Normalize directions: Set d(t) ンy ; (1) / //lr) | and d d /11///1///
5: Aggregate rays: set An (G( + (1 (O) )
6: Update: ((++1) + (()) 11(t) (G(!)))
4.2 HARMONIZED CONE GRADIENT DESCENT
With the theoretical foundation in Section we now propose HARMONTZed Cone gradient de- scent (HARMONIC), which ensures that A(G)) remains inside H. A representation matrix A (t) with respect to (²) is derived from 20m and GO ( x(0 Through the Double Description Method (Motzkin et al. 1953), the H-representation (2) is converted into the V-representation, given by the generating matrix [I()= (4) row at of G ((9) ()) 122 the current ray set is partitioned into R { a 0), R2 = {TT a 0}, and R, ={T a. π 0}. and new candidate rays are generated from all (m,,
...., .....,Tp " " Specifically, for each constraint
R, via (ag n x) π (a] --) .-. The resulting candidate is then pruned to obtain the final generator matrix IT(() Each G n an extreme ray of H())) and corresponds to column vector of in Definition Finally, the update Ar (G!)) obtained by averaging and normalizing the extreme rays of ((( to form a unified direction, and the direction then scaled by the sum of the inner products with all loss gradients. The proposed method is summarized in Algorithm
4.3 COMPARISON WITH RESPECT TO EXISTING APPROACHES
In this section, we illustrate how HARMONIC updates its search direction in comparison with the other approaches. We then examine a toy example in which we visualize the trajectory converging toward the Pareto front and compare it with the trajectories generated by the other methods. Finally, we demonstrate the advantages of HARMONIC through additional illustrative examples, building on the same experimental setting applied to the other methods in Section 3.2
CanRi
HORNORE
HARMOWK
2.10
-10
(a) NTK
(b) MGDA
(c) IMTL-G
(d) ConFIG
Figure 4: Unsatisfied cases of the existing methods. - denotes the gradient of each loss. -->
represents A(G) of the existing methods. represents A(G) of HARMONIC.

---

Published as a conference paper at ICLR 2026
Fig. visualizes, in dirres-imeniional parameter space with three loss gradients, how A()) is selected under each of the four methods in comparison with HARMONIC. We compare HARMONIC with NTK (Wang al., 2022), MGDA (Sener & Koltun 2018), IMTL-G (Liu et al. 2021b), and ConFIG (Liu et al. 2025). From Fig. 4a NTT can lead to updates that lie within the primal gradient cone K, f.e., the feasible region constructed by the loss gradients. However, such updates may still incur conflict, as the update direction can form an obtuse angle with some of the individual gradients. On the other hand, shown in Fig. 4b) MGDA produces updates that are not only feasible but also non-conffict. Yet, when the inner product between any pair of loss gradients is nonnegative, the update direction is biased toward the gradient with the smallest magnitude. In such cases, some loss terms may fail to decrease, since the update direction can align almost exclusively with the gradient of the smallest magnitude. IMTL--, illustrated in Fig- 4c may also result in updates that lie in an infeasible region, potentially inducing conflict. Finally, as shown in Fig. 4dl ConFIG ensures non- conflict updates; however, the resulting direction may fall outside the feasible region. Further details and discussions are provided in Appendix E.
(a) HARMONIC
(h) CAGrad
(c) Aligned-MTTI
(d) ConFIG
Figure 5: Visualization of optimization trajectories for the toy example. The starting points of the trajectories are denoted by (0.8, 0.8, 0.5, 0.5), (0.3, 0.9, 0.2, 0.8), and o (0.1, 0.8, 0.1, 0.7).
And the trajectories are shown fading from (start) to (end). The Pareto front is given by
We construct a toy example in a four-dimensional parameter space to evaluate whether recent MOO methods can successfully converge to the Pareto front for a three-objective minimization task. Comparing against recent MOO algorithms such as CAGrad (Liu et al., (2021a), Aligned- MTL (Senushkin et al., 2023), ConFIG, we evaluate HARMONIC under the same conditions. All methods are initialized from the same three starting points, and the optimization trajectories are visualized over 20,000 steps with learning rate of 0.001. Details of the experimental setup are provided in Appendix E.2
As shown in Fig. 5 HARMONIC achieves convergence to the Pareto front from all initial points. In contrast, since CAGrad does not guarantee non-conflict updates, its iterates are pulled toward the mean-loss gradient even after reaching the Pareto front, eventually converging regions dominated by single objective. On the other hand, the trajectories of Aligned-MTLL and ConFIG, both of which enforee non-conflicting updates, converge stably once they reach the Pareto front. Yet their updates can occasionally step outside the feasible region, causing premature convergence before the Pareto front for certain initializations.
Taken together, these comparisons highlight that HARMONIC resolves the challenges faced by exist- ing methods across all loss gradient scenarios. In Fig. [aa| we empirically demonstrate that updates with HARMONIC decrease all training losses without degrading test performance. Furthermore, Fig. 6b 6c show that, for all j. both 9j and 9-, have strictly positive cosine values with A(G), which only a feasible but also an non-conflict direction, indicating that A(G) always lies in the interior of both IK and IK*.
4.4 THEORETICAL ANALYSES
We show that HARMONIC can converge to a Pareto stationary point even in possibly nonconvex settings, and theoretically establish that the intersection of K. and KK always exists, excluding the trivial solution.
Theorem 3. Let L1, Izy..., .., Iem be differentiable and possibly nonconvex objectives, and assume that the total gradient ΣØ=1 VoLi is Lipschitz continuous with constant > O. Consider updates along HARMONIC, denoted by An (G), with step size 1) p. where M > 0 denotes the minimum co-
2
7

---

Published as a conference paper at ICLR 2026
3.
13
300
(a) HARMONIC
(b) HARMONIC (in K)
(c) HARMONIC (in IK")
Figure 6: (left) Training curve (a) and histograms of cosine values cos( ()) between A(G') and 9-g (b), and between A(G) and (c).
sine value between As GG) and each 9,. Then the sequence either converges a Pareto-stationary point, or, after T iterations, satisfies the bound
1
m D
T
VOL
11
-1,
T
nMT j=1
j=1
where T denotes the total number of iterations.
of ((1vVT) in the nonconvex setting, The detailed proof of the theorem is provided in the Appendix B.4
Theorem 3 states that HARMONIC converges to stationary point
Theorem 4. Let the convex gradient hull be defined as U = {G Rd X 2 O..
1 X = 1},
which is contained the primal gradient cone K. Then there always exists vector
U
arg min |
UEU
such that | > O. Consequently, u belongs to HI, and hence H. containing necessarily exists.
This guarantees that HARMONIC can always operate within a nontrivial harmonized cone, ensuring its applicability across all training scenarios. Additional details and proofs deferred to Ap- pendix B.5
EXPERIMENTS
In this section, we compare HARMONIC with various methods, including MultiAdam Yao al.
2023
LRA (Wang et al. 2021). ReLoBRaLo (Bischof & Kraus 2025), MGDA (Sener & Koltun
2018
PCGrad (Yu al. 2020). CAGrad (Liu al. 2021a). IMTL-G (Liu al. 2021b), Aligned-
MTL
Senushkin et al. 2023). and ConFIG (Liu et al. 20225.
We evaluate performance the PINNacle benchmark (Zhongkai et al. 2024) for PDEs and the A-PINN (Yuan et al. 2022) benchmark for the integro-differential equation (IDE). these exper- iments, consider five benchmark equations: popular PDEs (Waveld-C, Poisson2d-C), high- dimensional PDE (HNd), an IDE (Volterrald), and inverse problem (HInv). All experimental set- tings follow PINNacle, except that trained for 50,000 iterations. All methods are evaluated under two learning rate settings (0.001, and 0.0001), and report the best performance for each method. Details of the experimental settings and additional results on more equations are provided in Appendix
network architecture with 50 neurons and 3 layers,
and H.I
Comparision on Benchmark Equations. In Table we highlight the best and the second best methods based on the mean over five random seeds. The numbers of objectives for each bench- mark are: 4 (Waveld-C, Volterrald), and (Poisson2d-C. HNd, HInv). All objectives are treated as independent losses. Experimental results demonstrate that HARMONIC generally outperforms all benchmarks. Consistent with prior empirical findings namely, that DCGD (Hwang & Lim) 2024) (reporting results with two losses) and ConFIG (Liu al. (025)) (reporting results with three losses) achieve performance gains through non-conflict strategies- our experiments confirm that mon-conffict-batt methods outperform reweighting approaches not only with three but also with four loss terms. In particular, for Poisson2d-C. all methods except HARMONIC exhibit poor perfor- mance, whereas HARMONIC achieves substantially superior performance.
8

---

Published as a conference paper at ICLR 2026
Table 2: Average relative I,2 errors (with standard deviations) estimated across random seeds for each algorithm.
Method
Waveld-C
Polsson2d-C
HNd
Volterrald
Hinv
MultiAdam
1.0958 (0. 1004)
0.6733 (0.0820)
0.0040 (0.0028)
0.0004 (0.0002)
0.1446 (0.0245)
LRA
0.2608 (0.0542)
0.5712 (0.15
0.0020 (0.0011)
0,0070 (0.0065)
0.0546 (0.0131)
ReLoBRaLo
0.3341 (0.0089)
0.6602 (0.0221)
0.0004 (0. [0.0000)
0.0002 (0.0000)
1.5999 (0.0576)
MGDA
0.3995 (0.0135)
0.6777 (0.0239)
0.0009 (0.0001)
0.0011 (0.0019)
2.4224 (1.6765)
PCGrad
0.2256 (0.0701)
0.6918 (0.0190)
0.0005 (0.0002)
0.0004 (0.0004)
0.0660 (0.0152)
CAGrad
0.1186 (0.0071)
0.6969 (0.0197)
0.0005 (0.0001)
0.0014 (0.0024)
0.6435 (1.0862)
IMTL-G
0.2724 (0. (373)
0.5039 (0.0684)
0.0006 (0.0001)
0.0735 (0.0989)
0.0498 (0.0075)
Aligned-MTI.
0.2302 (0.0514)
0.2847 (0. 2363)
0.0006 (0.0002)
0.0037 (9.0044)
0.0804 (0.0457)
CoaFiG
0,0668 (0.0279)
0.6856 (0.0274)
0.0006 (0.0001)
0.0018 (0.0024)
0.0166 (0.0068)
HARMONIC
0.0655 (0.0293)
0.0214 (0.0179)
0.0005 (0.0000)
0.0003 (0.0001)
0.0461 (0.0098)
Mitigating Unsatisfied Cases. We conducted additional experiments to assess the effect of con- straining A(G)) lie inside H.. In Fig. 77 the left panels of Fig. [7a] and Fig. 7b) visualize, for ReLoBRaLo and ConFIG respectively on Poisson2d-C, whether A(G) enters IH during training. In both cases, the test relative L error decreases at first but then remains nearly constant beyond a cer- tain point, indicating that the training has settled into an undesirable suboptimal point. By contrast, the right panels of Fig. 7a and Fig. 7b]ssatt from the same initialization but, whenever the baseline update would drive A(G outside H., we instead apply HARMONIC to keep the update inside H. In- voking Theorem 2
2 we mathematically certify whether A(G) lies inside or outside H and, based on this certification, choose the update rule used to form A(G). Under this intervention, the test relative L error converges close to zero, indicating that optimization does not stall at a suboptimal point but instead progresses toward more desirable minimum.
.1611
#61 W
1900
1480 10:00
(a) (left) ReLoBRaLo only (right) with HARMONIC
(b) (left) ConFIG only (right) with HARMONIC
Figure 7: Training curve for Poisson2d-C, where denotes A(G) E Hi and denotes A(G)) & H.
We then extend the visualization experiment to the same settings as Table 2 For each baseline method (ReLoBRaloo, ConFIG, and CAGrad), whenever HH) we enforce the update to remain inside HI by applying HARMONIC. We denote the resulting variants HH-ReLBBRaLo, H--CArrad, and H-ConFIG, and their performance is reported in Table 3 The results in Table shoow that these variants consistently outperform their original counterparts. In particular, all methods demonstrate performance improvements on the benchmark equations, with especially substantial gains observed for Poisson2d-C.
Table 3: Average relative La errors (with standard deviations) estimated across 5 random seeds for each algorithm. Red-colored entries indicate results that are better than their counterparts.
Method
Waveld-C
Poisson2d-C
HNd
Volterrald
Hin
H-RoLoBRaLo
0.3274 (0.0392)
0.09418 (0.1763)
0.0030 (0.0056)
0.0179 (0.0351)
0.0167 (0.0126)
H-CAGrad
0.1265 (0.0463)
0.0312 (0.0122)
0.0013 (0.0024)
0.0007 (0.0004)
0.0518 (0.0070)
H. ConFiG
0.0407 (0.0149) 0.0094 (0.0022)
0.0006 (0.0002)
0.0064 (0.0096)
0.0558 (0.0124)
Computational Cost. Tepoot the average seconds per 100 epochs for each method across the benchmark datasets. Across all benchmarks. MGDA and CAGrad consistently required longer runtimes, as each epoch involves an iterative optimization procedure. By comparison, HARMONIC, despite enforcing both the feasible and non-confticr conditions, is generally faster, or at least main- tains a computational cost comparable to ConFIG. In particular, computing the V-representation of the 2m x m matrix in HARMONIC is practical in terms of wall-clock time. Appendix H.3 H.3 further confirms that HARMONIC scales similarly even with eight loss components.

---

Published as a conference paper at ICLR 2026
Table Average computational costs (in seconds) with standard deviations from 100 epochs across 100 repeats. The worst methods are highlighted.
Method
Waveld-C
Poisson2d-C
HNd
Volterrald
Hllv
MultiAdam
1.53 (0.06)
1.08 (0.07)
5.29 (0.01)
1.03 (0.05)
2.36 (0.01)
LRA
1.39 (0.05)
1.03 (0.04)
5.46 (0.03)
1.10 (0.05)
2.61 (0.05)
ReLoBRaLo
0.51 (0.01)
0.46 (0.02)
2.24 (0.01)
0.41 (0.02)
1.09 (0.02)
MGDA
3.19 (0.52)
3.00 (0.53)
6.40 (0.32)
4.01 (1.38)
3.13 (0.13)
PCGrad
1.58 (0.07)
1.16 (0.08)
5.40 (0.02)
1.15 (0.05)
2.48 (0.03)
CAGrad
2.48 (0.24)
1.90 (0.14)
6.41
1.91 (0.21)
3.45 (0.13)
IMTL-G
1.56 (0.06)
1.09 (0.05)
5.34 (0.01)
1.02 (0.05)
2.41 (0.02)
Aligned-MTI.
1.49 (0.07)
0.97 (0.03)
5.35 (0.01)
1.01 (0.05)
2,42 (0.02)
ConFIG
1.54 (0.05)
1.12 (0.06)
5.39 (0.03)
1.06 (0.04)
2.44 (0.03)
HARMONIC
1.56 (0.07)
(0.07)
5.35 (0.01)
(0.03)
2.42 (0.02)
CONCLUSIONS
this work, we proposed the concept of the harmonized cone to address challenges that arise from multiple losses in the training of PINNs. Defined the intersection of the primal and dual cones of per-loss gradients, the harmonized cone characterizes a region that simultancously guarantees fea- sible and non-conflict. Building on this geometric analysis, introduced HARMONIC, gradient- based method that leverages the Double Description method to aggregate extreme rays, therehy enforcing that the update direction remains within the interior of the harmonized cone. From the- oretical perspective, we established the existence of a nontrivial harmonized cone and proved that HARMONIC converges Pareto-stationary solutions in nonconvex settings. Empirical results on the PDEs and IDE benchmarks demonstrate that HARMONIC consistently achieves either superior performance compared the state-of-the-art baselines. Furthermore, showed that HARMONIC incurs no substantial computational overhead compared to the existing methods.
For future work, compelling direction is to incorporate objective preferences, as studied in MOO (Chen et al. (2024). to guide update directions within the harmonized cone. Another promising direc- tion is extend the framework beyond first-order optimization and explore second-order methods for PINNs (Wang et al., 2025)
7 REPRODUCTIBIITY STATEMENT
The proofs of the main results are provided in Appendix A and B Details on how the figures in the main text were generated are described in Appendix C D and Appendix [F) provides complete descriptions of the benchmarks used in our experiments as well as the experimental settings, while Appendix GerenndxxGppreseee presents the implementation details of the baseline methods considered for comparison. In addition, we provide the benchmark data files and fully reproducible code the supplementary material.
8 ACKNOWLEDGEMENN
This research was supported in part by the National Research Foundation of Korea grant funded by the Korea government (MSIT) (Grant number RS-2020-NR049544, RS-2022-NR06877 and RS- 2024-00361377) and in part by the Korea Evaluation Institute of Industrial Technology grant funded by the Korea government (MOTIR) (Grant number RS-2002525458052),
10

---

Published as a conference paper at ICLR 2026
REFERENCES
Klaus-Jürgen Bathe. Finite element procedures. Klaus-Jurgen Bathe, 2006.
Atilim Gunes Baydin, Barak A Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind. Automatic differentiation in machine learning: a survey. Journal of machine learning research, 18(153):1-43, 2018.
Rafael Bischof and Michael A Kraus. Multi-objective loss balancing for physics-informed deep learning. Computer Methods in Applied Mechanics and Engineering, 439: 199111 2025.
John Boyd. Chebyshev and Fourier spectral methods. Courier Corporation, 2001.
Lisha Chen, A Saif, Yanning Shen, and Tianyi Chen. Ferero: A flexible framework for preference- guided oulti-objective learning. Advances in Neural Information Processing Systems, 37:18758- 18805, 2024.
Zhao Chen, Vijay Badrinarayanan, Chen- Yu Lee, and Andrew Rabinovich. Gradnorm: Gradient normalization for adaptive loss balancing in deep multitask networks. In International conference on machine learning, pp. 703. PMLR, 2018.
Mohammad Sadegh Eshaghi, Cosmin Anitescu, Manish Thombre, Yizheng Wang. Xiaoying Zhuang, and Timon Rabczuk. Variational physics-informed neural operator (vino) for solving partial differential equations. Computer Methods in Applied Mechanics and Engineering, 437: 117785, 2025.
Xavier Glorot and Yoshua Bengio. Understanding the difficulty of training deep feedforward neural networks. In Proceedings of the thirteenth international conference on artificial intelligence and statistics, pp. 249-256. JMLR Workshop and Conference Proceedings, 2010.
Christian Grossmann, Hans-Gorg Roos, and Martin Stynes. Numerical Treatment of Partial Differ- ential Equations: Translated and revised by Martin Stynes. Springer, 2007.
Ali Heydari, Craig Thompson, and Asif Mehmood. Softadapt: Techniques for adaptive loss weighting of neural networks with multi part loss functions. arXiv preprint arXiv:1912.123555 2019.
Youngsik Hwang and Dongyoung Lim. Dual cone gradient descent for training physics-informed neural networks. Advances in Neural Information Processing Systems, 37:98563-995 998166539999,----- 2024.
George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, 3(6):422-440, 2021.
Dmitrii Kochkov, Jamie A Smith, Ayya Alieva, Qing Wang. Michael P Brenner, and Stephan Hoyer. Machine caaming-accelennan computational fluid dynamics. Proceedings of the National 2021.
Academy of Sciences, 118(21):c210178444
Hyunwoo Lee, Hayoung Choi, and Hyunju Kim. Robust weight initialization for tanh neural net- works with fixed point analysis. In The Thirteenth International Conference Learning Repre- sentations, 2025.
Xi Lin, Hui-Ling Zhen, Zhenhua Li, Qing-Fu Zhang, and Sam Kwong. Pareto multi-task learning. Advances in neural information processing systems, 32, 2019.
Bo Liu, Xingchao Liu, Xiaojie Jin, Peter Stone, and Qiang Liu. Conflict-averse gradient descent for multi-task learning. Advances in Neural Information Processing Systems, 34:18878-18890, 2021a.
Liyang Liu, Li, Zhanghui Kuang, Jing- Hao Xue, Yimin Chen, Wenming Yang, Qingmin Liao, and Wayne Zhang. Towards impartial multi-task learning. In International Conference on Learning Representations, 2021b.
Qiang Liu, Mengyu Chu, and Nils Thuerey. ConFIG: Towards conflict-free training of physics informed neural networks. In The Thirteenth International Conference on Learning Representa- fions, 2025.
11

---

Published as a conference paper at ICLR 2026
Keith W Morton and David Francis Mayers. Numerical solution of partial differential equations: an introduction. Cambridge university press, 2005.
T. S. Motzkin, H. Raiffa, G. L. Thompson, and R. M. Thrall. The double description method. In H. W. Kuhn and A. W. Tucker (eds.), Contributionn to the Theory of Games, Volume II, number 28 in Annals of Mathematies Studies, pp. 51-73. Princeton University Press, Princeton, NJ, 1953.
Yurii Nesterov. Lectures on Convex Optimization. Springer Publishing Company, Incorporated, 2nd edition, 2018. ISBN 3319915770.
J Papež, Liesen, and Z Strakoš. Distribution of the discretization and algebraic error in numerical solution of partial differential equations. Linear Algebra and its Applications, 449:89-114, 2014.
Maziar Raissi, Paris Perdikaris, and George Karniadakis. Physies-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-707, 2019.
Maziar Raissi, Alireza Yazdani, and George Em Karniadakis. Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations. Science, 3cience,367(6481 (666 2020.
Francisco Sahli Costabal, Yibo Yang, Paris Perdikaris, Daniel E. Hurtado, and Ellen Kuhl. Physics- informed neural networks for cardiac activation mapping. Frontiers in Physics, 8:42, 2020.
Ozan Sener and Vladlen Koltun. Multi-task learning as multi-ohjectivv optimization. Advances in neural information processing systems, 31, 2018.
Dmitry Senushkin, Nikolay Patakin, Arseny Kuznetsov, and Anton Konushin. Independent compo- nent alignment for multi-task learning. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 20083-20093, 2023.
Rui Wang, Karthik Kashinath, Mustafa Mustafa, Adrian Albert, and Rose Yu. Towards physics- informed deep learning for turbulent flow prediction. In Proceedings of the 26th ACM SIGKDD international conference knowledge discovery & data mining, pp. 1457-1466, 2020.
Sifan Wang, Yujun Teng, and Paris Perdikaris. Understanding and mitigating gradient flow patholo- gies in physics-informed neural networks. SIAM Journal on Scientific Computing, A3355 A3081, 2021.
Sifan Wang, Xinling Yu, and Paris Perdikaris. When and why pinns fail train: A neural tangent kernel perspective. Journal of Computational Physics, 449:110768 2022.
Sifan Wang, Ananyae Kumar Bhartari, Bowen Li, and Paris Perdikaris. Gradient alignment in physics-informed neural networks: A second-order optimization perspective. arXiv preprint arXiv:2502.00666 00604, 2025.
Chenxi Wu, Min Zhu, Qinyang Tan, Yadhu Kartha, and Lu Lu. A comprehensive study of non- adaptive and residual- based adaptive sampling for physics-informed neural networks. Computer Methods in Applied Mechanics and Engineering, 403: 15671, 2023. ISSN 0045-7825.
Zixue Xiang, Wei Peng, Xu Liu, and Wen Yao. Self-adaptive loss balanced physics-informe neural networks. Neurocomputing, 496: 1 1-34, 2022.
Jiachen Yao, Chang Su, Zhongkai Hao, Songming Liu, Hang Su, and Jun Zhu. Multiadam: Parameter-wise scale-invariant optimizer for multiscale training of physics-informed neural net- works. In International Conference on Machine Learning, pp. 39702-3977 PMLR, 2023.
Bing Yu et al. The deep ritz method: deep learning- ceaseing-based numerical algorithm for solving varia- tional problems. Communications in Mathematics and Statistics, 6(1): 2018
Jeremy Yu, Lu Lu, Xuhui Meng, and George Em Karniadakis. Gradient-enhanced physics-informt neural networks for forward and inverse pde problems. Computer Methods in Applied Mechanics and Engineering, 393:1 14823, 2022.
12

---

Published as a conference paper at ICLR 2026
Tianhe Yu, Saurabh Kumar, Abhishek Gupta, Sergey Levine, Karol Hausman, and Chelsea Finn. Gradient surgery for multi-task learning. Advances in neural information processing systems, 33: 5824-5836, 2020.
Lei Yuan, Yi-Qing Ni, Xiang-Yun Yun Deng, and Shuo Hao. A-pinn: Auxiliary physies informed neural networks for forward and inverse problems of nonlinear integro-differential equations. Journal of Computational Physics, 462:1 462:111260, 2022.
Hao Zhongkai, Jiachen Yao, Chang Su, Hang Su, Ziao Wang, Fanzhi Lu, Zeyu Xia, Yichi Zhang, Songming Liu, Lu Lu, et al. Pinnacle: A comprehensive benchmark of physics-informed neu- ral networks for solving pdes. Advances in Neural Information Processing Systems, 37:76721- 76774, 2024.
Taotao Zhou, Xiaoge Zhang, Enrique Lopez Droguett, and Ali Mosleh. A generic physics-informed neural network-basee framework for reliability assessment of multi-state systems. Reliability Engineering System Safety, 229: 108835, 2023.
Yinhao Zhu, Nicholas Zabaras, Phaedon- Stelioos Koutsourelakis, and Paris Perdikaris. Physics- constrained deep learning for high- high-dimensional surrogate modeling and uncertainty quantification without labeled data, Journal of Computational Physics, 394:56-81, 2019.
13

---

Published as a conference paper at ICLR 2026
A PROOFS FOR SECTION 3
A.1 DETAIL PROOF FOR DEFINITION
Definition (Primal Gradient Cone). We define the primal gradient cone as
K (G> X RT).
or equivalently
K : (o R D mm}
(2)
where RV denotes the m-dimensional nonnegative orthant.
Proof. A cone K. is generated by conic combinations of finitely many column vectors of the V- representation G.
Theorem (Weyl's Theorem), Any nonempty finitely generated cone is finitely constrained.
By Theorem S the cone IK can also be expressed in terms of an H-representation H with finitely many columns.
Corollary (Weyl Duality). Given matrix G d x consider K {GX X R fxr Rd G 0. Then K" = L, L K, where denotes dual relationship.
H corresponds the V-representation of the dual cone K Consider GT and its pseudo-inverse D. Since G has full rank, we have
By Corollary [
G D I
(3)
This implies
9k 9- 1,
k J k: *
(4)
g) 9-3 O,
(5)
This is equivalent to the process of activating all constraints, except the one associated with the j-th column of G, to yield the extreme ray. Hence, each g is extreme ray of the dual cone K*
Consequently, can be regarded as one of the valid H-representations of K, and we may equiva- lently write
K = (x E IRd
D x 2 P.
A.2 DETAIL PROOF FOR DEFINITION
Definition (Dual Gradient Cone), The dual gradient cone of K. is defined as
{y Rdd G y 2 0mts
or equivalently
K" (Dw E R
Proof. By definition, the dual cone K can represented in the H-representation by finitely many column vectors of G.
Theorem (Minkowski's Theorem), Any polyhedral cone is nonempty and finitely generated.
By Theorem [ the inequality representation
{y Rd G G 0m)
admits an equivalent V-representation V with finitely many columns.
By Corollary [ corresponds the H-representation of the primal cone KK. Furthermore, from the result of Proof A.I since D is H-representation of K, it can equivalently be regarded valid V-representation of K.
14

---

Published as a conference paper at ICLR 2026
B PROOFS FOR SECTION
B.I DETAIL PROOF FOR DEFINITION Definition 3
(Harmonized Cone). The set
H = (= E Rdd
=² 0 D'= z 2 m
(10)
or equivalently
H == Rox a E RR)
(11)
called the harmonized cone of G, where R E pdxp denotes the generating matrix of H.
Proof. By definition, the harmonized cone Hi is the intersection of KK and IK.
Since H must simultaneously satisfy the constraints given by the H-representations of both IK and K it can written as
= Rdd G
Om: D Om)
(12)
By Theorem [ the inequality representation admits an equivalent V-representation R with finitely many columns.
B.2 DETAIL PROOF FOR THEOREM
Theorem 1. Suppose that the gradient matrix at iteration denoted by GO and let H()) be the harmonized cone induced by ((t), Define = ( E Ry as scaling vector, where each component X is the scaling factor corresponding to the loss Lş- Then, aeeasiile update direction (t)) e H[(t) if and only if
G (
2 Om-
(13)
Proof. A feasible update direction (L)) (1), which lies in the primal cone K, must satisfy Defini- tion 3 in order to belong to the harmonized cone.
we denote
G (t)), then must hold that
G(OT (() 0m
(14)
(()T(())) 2 Om-
(15)
Since (()) G(t)=Im, for (+) defined in RT, the condition reduces to requiring G x())(()))
2
Om Thus, x belongs to the harmonized cone if and only if this condition is satisfied.
B.3 DETAIL PROOF FOR THEOREM 2
Theorem 2. Suppose that the
pseudo-inverse of transposed gradient matrix at iteration t is denoted by D(1), and let (()) be the harmonized cone induced by D(t), Define w(()= [1] W E RT as a scaling vector, where each component wy is the scaling fac- tor corresponding to g Then, a non-conftict update direction ( Hi( if and only if
D D w > 0m-
(16)
Proof. A non-conflict update direction D (1) which lies in the dual cone K*, must satisfy Defi nition in order to belong to the harmonized cone.
we denote ()() then must hold that
( (t)
> 0 0ms
(17)
( (t) (1) > 0m-
(18)
Since (()) (+)) for (1) defined in RT, the condition reduces to requiring ( D (1)) (t)
>
Om Thus, (1), (1) belongs the harmonized cone if and only if this condition is satisfied.
15

---

Published as a conference paper at ICLR 2026
B.4 DETAIL PROOF FOR THEOREM
Theorem 3. Let L1, Lg..... Li be differentiable and possibly nonconvex objectives, and assume that the total gradient Ej=1 VoL, is Lipschitz continuous with constant > 0. Consider updates along HARMONIC, denoted by (G), with step size rj S /. where M denotes the minimum cosine value between A, (G) and each 95. Then the sequence either converges to a Pareto-stationary point, after T iterations, satisfies the bound
T Σ
10 S nMT
Σ - - j=1
(19)
where T denotes the total number of iterations.
Proof. By the descent lemma
(Nsttrovv. 2018 an update along An (G(()) with a step size
S 22
gives us
(++1)
< Σ 1 An (G
(20)
j=1 j=1
where the superscript indicates the optimization iteration index.
Recall from Algorithm that d(t) == ( /d)) /||d0/ satisfies || = and d H()) C K*()), so
that 9j (t)T (t)) > 0 for all j. Also, Algorithm defines
(G())) (1 GLO d(0) d Σ
9
(O)
d d
(21)
j=1
Therefore
A (G
= |1 || (Σ
m
(t)T d(t)
d
9
j=1
(1T d(t) )d
>
(0) d(t)
9j
(U)T g00
||
1=1
(V, 9
(T) g(t) 0)
1=1
A (G
An (G
j=1
9 (OT (,
(22)
1=1
9
(c)
A [G
16

---

Published as a conference paper at ICLR 2026
And the triangle inequality where
||
9/l|, we have
m
(++1)
(1)
!
(G
Σ
m
L
"
Σ
m
19
(0 2
[9 ( )²
9, An GG
2
j=1
1=1
11
m
9, (n(U0)))
$
Σ
9
2
|
AA (G())]
j=1
j=1
9j (G(
(9)
VoL (1) 12-
g An (G )²
(23)
2
j=1
Summing this inequality over t = 0,1,... T
results in
T-1
T-1
1
|1 9) An (G
VoL 10 12
Σ
( ).
(24)
t=
0= 1=1
By defining minogisT-1
9. (o((0) (1g
)²
M, we have
T-1
2 m
VoL 112 v
Σ
(
- Lm7)
(25)
nM
t=
j=1
j=1
Finally, averaging both sides leads to
T-1 m
2
Σ j=1
m
(L)
(0)
L ).
(26)
nMT
As co, either the right-hand side in the inequality goes to 0, implying that the minimal gradient norm converges to zero, which means a stationary point, or the inequality above holds.
B.5 DETAIL PROOF FOR THEOREM
Theorem 4. Let the convex gradient hull be defined as U= { GA E Rd
A Om: I A
1},
which is contained in the primal gradient cone K. Then there always exists a vector
u arg min |ze|||
(27)
ueU
such that ||u*|| || > 0. Consequently, a belongs to H. and hence HI containing u² necessarily exists.
Proof. For each i, define
U. = -"
(28)
Since U is convex, it holds that
u" + co, E U, VEE [0. 1].
(29)
By optimality of ut we have Expanding the left-hand side gives
|u" - =vuil? > 1|
(30)
20(u* v,) + + 2
(31)
For sufficiently small > 0, this implies
u v 2 0.
(32)
Substituting v; = 9: yields
| 2 0.
(33)
Because G is full rank, u cannot be trivial solution. Moreover, since 2 is nontrivial, the set HI containing u² cannot be empty. Equivalently, the cone generated by G and its dual necessarily intersect.
17

---

Published as a conference paper at ICLR 2026
C DETAILED DISCUSSION FOR SECTION 1
In this section, we provide detailed account of the bi-objective toy example introduced in Fig.
We consider 0 = (01,02) 02) E g² and define two objectives:
L (0) (ci(0s) Sia 0) c2(02) [1b(0)) 8.
(34)
(2 ci (02) f2a (2(2) (0),
(35)
where
(0) = log max (0.500000000000 - - tanh -02), 5x 10 10)00006))))))))) +
(36)
fza (0) (max (0.5(- (00 3) tanh( (anh(-82))))) + 5x10 + 6,
(37)
f10(0) (-01 7) + 0.1 8)²
20,
(38)
10
(2) (-01 7)² + (. 02 8)2
20,
(39)
10
and the switching coefficients
ca (02) =
0),
C2 (2( = -aa(tant(-0-055 0).
(40)
The composite loss is then
Lo(0) = (0) + (0).
(41)
Optimization is carried out using Adam with learning rate 0.001 for 80,000 epochs. We initialize (-8.5, (.5)) and (9.0,94 9.0), and trajectories are plotted every 1,000 epochs. The hyperparameter is varied investigate how rescaling L affects convergence behavior.
Pareto Opt HARMONIC 9or Osc
13-
Parwto
Paretn HARMONICC
HARMONIIC
B
11
-4
4
(a) 0.1
0.5
(c) 8 1.5
Pareto Opt HARMONIC
Pareto Opt HARMORIC
24
-33
46
-33
77
-5
91
61
(d)
(e) 8.0
Figure 8: Optimization trajectories for the toy example under different values of s. Initialization points are indicated by o and the Pareto front shown in indicates gee- denotes gef> and represents HARMONIC.
The observations across different settings are summarized as follows:
0.1: From the initial point (-8.5, 7.5). all methods converged to the Pareto set. From
(9.0,9.0 9.0), only HARMONIC succeeded in reaching the Pareto set.
8 = 0.5: From (-8.5, 7.5), both HARMONIC and the non-conflict update 9.f converged to
the Pareto set. From (9.0, 9.0), only HARMONIC reached the Pareto set.
18

---

Published as a conference paper at ICLR 2026
1.5: From (-8.5, 7.5), HARMONIC alone converged to the Pareto set. From (9.0, 9.0), both HARMONIC and 9ef converged.
8 6.0: From (-8.5, 7.5), only HARMONIC converged. From (9.0, 9.0), all methods
successfully converged to the Pareto set.
8 8.0: From (-8.5, only HARMONIC converged. From (9.0, 9.0), both HAR- MONIC and the conic-combination update Jec converged.
Overall, in all tested cases HARMONIC consistently converged to the Pareto set, whereas alterna- tive update rules exhibited sensitivity to both the scaling factor and the initialization point. This demonstrates the necessity of jointly enforcing feasible and non-conflict, as argued in Section] [
19

---

Published as a conference paper at ICLR 2026
D EXPERIMENTAL SETUP FOR FIG. 2 FIG. 3 AND FIG.
6
To produce the training curves and cosine-value histograms in Fig. 2 Fig. 3 and Fig. 6 trained a physics-informed neural network (PINN) on the one-dimensional Burgers equation under the fol- lowing controlled setting. The network architecture consisted of four fully connected layers with 32 neurons per layer using hyperbolic tangent activation function, and weights were initialized by Xavier normal initialization with zero bias. Training was performed for total of 2,000 epochs using stochastic gradient descent (SGD) with learning rate 0.001. We sampled 10,000 points each for the PDE residual, initial condition, and boundary condition terms.
20

---

Published as a conference paper at ICLR 2026
E DETAILED DISCUSSION FOR SECTION 4
E. UNSATISFIED CASES OF THE EXISTING METHODS (SECTION 4.3
In this section, we provide detailed numerical evidence complementing the visualizations. For each method, we present the gradient matrix G' = [g 92, gs| We then report the update direction of HARMONIC A₁ (G) and that of the comparison method A. Next, we compute inner products with two types of test vectors:
K* -test: the inner products (g). A(G)) between each gradient 9j and A(G). If any of them is negative, update conflicting.
K-test: the inner products A A())), where 9-3 denotes the defining rays of the dual gradient cone K" If any of them is negative, the update lies in an infeasible region. A degenerate case where only one inner product equals while the others are 0 indicates that the update ignores all but one gradient.
NTK HARMONIC
0.0
-0.5
1.0
0.5
0.0
-1.0
0.5 8
-1.0
0.0
0.5
Figure 9: Unsatisfied cases of NTK. - denotes the gradient of each loss. represents (G) of the existing methods. represents A, (G) of HARMONIC.
NTK (Fig. 9
0.4278
-0.8925
0.1431
G²
0.5179
0.7754
0.3613
-0.0979
0.9941
-0.0472
(G = [-0.4067, 0.1670, 0.6891]
A.(G) = |-0.0026, 0.1524]
testet
Ax) = (0.1236, 2.1667,0017338 > 0, (g),A.(G) = (-0. 03800.00.000. 0.2804, 0.2837).
NTK induces conflict with 91.
MGDA (Fig, 10).
T0.0137
0.1432
0.3131
0.3571
0.5491
0.1414
00.8233
0.9361
0.0552
(G) (0.4320, 0.6414, 0.3868J
A.(G) [0.0137, 0.143,
K-test:
An (G)) (0.9674, 0.4974, 0.2455) > 0 Ag (1.0000, 0.000,,00000 >
MGDA produces degenerate update, aligning only with 91-
21

---

Published as a conference paper at ICLR 2026
MGDA HARMONIC
0.0
-0.5
1.0
0.5
0.0
0.0 -1.0
0.5 6
-1.0
0.5
Figure 10: Unsatisfied cases of MGDA. - denotes the gradient of each loss. represents (G))) of the existing methods. -- represents A, (G) of HARMONIC.
IMTL-G HARMONIC
-05
1.0
0.5
-1.0
0.5
-0.5
0.0
0.5
Figure 11: Unsatisfied cases of IMTL-G. denotes the gradient of each loss. represents (.) of the existing methods. represents A, (G) of HARMONIC.
IMTL-G (Fig. 11
0.1323
1.0037
-0.20901
G
0.3367
0.4206
0.3984
0.9204
0.5455
0.8364
(G) = 1-044360, 0.6668, 0.3361
A.(G) = [0.8040, 00.127,
K-test:
(g AG)) = (0.3225, 0.4974, 0.495,,02445) > O, A.G)) = (- 4.6542, -2.3954).
-testt
(g) = (0.6566, 0.5612, 1.0462) > 0, Ag) (-0.7451, 4830, - 00.9788) < 0.
IMTL-G yields infeasible and conflicting updates.
22

---

Published as a conference paper at ICLR 2026
ConFIG
HARMONIC
0.0
-0.5
1.0
0.5
0.0
-1.0
0.5 6
-1.0
0.0
0.5
Figure 12: Unsatisfied cases of ConFIG. - denotes the gradient each loss. -- represents (G)) of the existing methods. represents A, (G) of HARMONIC.
ConFIG Fiig.
-
0.1807
0.2457
0.9523
G²
=
0.4623
0.6682
-0.5829
-..055
0.5583
0.8230
(G) (0.2071, 0.w.69.7.70.... 0.3001
A.(G)
[5.5958,
00.4348 2.22411
K ttet::
An = (0.4187, 0.6966, 0.3735) > 0, Ac()) = (75.9019, 26. 7388, 166.1925).
ConFIG violates feasibility.
Overall, these numerical checks confirm the unsatisfied cases of existing approaches, and demon- strate that HARMONIC uniquely guarantees both feasible and non-conflict simultaneously.
E.2 TOY EXAMPLE FOR COMPARING MULTIIGRADIENT METHODS (SECTION 4.3
To provide a controlled environment for analyzing gradient interactions among objectives, struct a three-objective optimization problem
f : R R x
(z) T2. Fz, .A) + f(x) = ((()
with the search domain restricted [0, 1Ja
Radial terms in (4) We begin by defining three quadratic distances
g( = (x3 0.5) + (x4 - 0.5)², 91 (æ) = (x3 0.45)² + (x4 0.55)² (2(2) = (x3 0.55)² + (x4 0.45)
(42)
(43)
(44)
The terms and ge create two nearby basins of attraction that determine which objective, fi or f2. dominates the local behavior.
Switching mechanism in (2) We introduce smooth logarithmic transform
(ba) log (nax((a),e)) + 6,
x 10-6
(45)
and two gating functions
(x) max(tanhi max(tann(20(12 0.5)), 0),
(2(2)
0.5)), 0). (46)
23

---

Published as a conference paper at ICLR 2026
These gates satisfy 61 for 22 > 0.5 and cg for 22 0.5, switching the curvature of the landscape across the line 22 = 0.5.
Using these components, we define four shape functions:
21
Joy (x) (0..3( 0.5) 7) - tanh -20002 0.5))).
(47)
22
log (z) = @(0.5 0.5) + 3) + tanh -20 - 0.5)) + +2,
(48)
(x) = (-20(x) - 0.5) 0.5) +)² + 0.1( 10.1-200(2:)) 0.5) 20)
20,
(49)
10
(x) (-20(x1 - 0.5) 7)2
(.10-20((2))0 0.5) 20)2
20.
(50)
10
The two composite shape functions are
21 = 21 01 + 2
22
22 c1
log
+ 22 C2-
(51)
Definition of the objectives ments), we define
With a global scaling parameter > (we use s = 6 in our experi-
fx()) = (x)
10 911(a);
(52)
f2(x) = 22(x²
(2(x),
(53)
10
fs(x)) +1+ g(=)) sin nh 21
sin 2
(54)
The optimization problem is therefore
min (Ji(z), f2(x), fs(x)) =
(55)
Sign structure and role reversal of gi and 92 For typical initial points used in experiments, the values of fi and f2 are positive, so reducing Si and la decreases both objectives. Along the Pareto front, however, f1 and f2 become negative, which reverses the effect of the multiplicative factors: and then act as maximization terms. This sign change results in a characteristic "role switch" of the gradients in (4),, creating sharp transitions in the interaction between the objectives.
Visualization To illustrate the structural properties of the toy problem, we additionally provide contour plots of the individual objectives under three representative slices of the four-dimensional domain. Each slice fixes two coordinates of æ and visualizes f1. f2, and fa as functions of the remaining two coordinates:
(1) Slice in (I1, x2) (Fig. 13ah : Fix (I3, =1) = (0.5, 0.5), and visualize fi, f2, and fa over (I1, [0, This slice highlights how the log-to-quadratic transition across 22 0.5 shapes the curvature of the first two objectives.
(2) Initial-condition slice in (I3, 24) (Fig, 13b): Fix (z1, x2) = (0.3, 0.8) and visualize fis f2. and fs over (I3, x.) 1J2 For this configuration, f1 and fa are positive, so the factors 91 and 92 act as minimization terms, creating distinct attraction regions centered near (0.45, 0.55) and (0.55, 0.45).
(3) Pareto-condition slice in (Iz, .1) (Fig. 13c: Fix (I1, .2) = (0.3, 0.1) and visualize fis fe, and fa over x1) [0, 1]2. In this region, fi and f2 become negative, reversing the effect of 9i and 92 so that they now behave as maximization terms. This results a qualitatively different landscape and induces sharp gradient-direction changes relevant to multi-gradient optimization.
These three complementary slices allow us to clearly visualize (i) the curvature switching in (2), (ii) the initial positive-objective regime in (æz, 2.a), and (iii) the lignflipped regime along the Pareto front.
24

---

Published as a conference paper at ICLR 2026
(a) Slice (I1, x2)
(b) Initial-condition slice in (I3, 274)
(c) Pareto-condition slice in (23, Ji4)
Figure 13: Three-dimensional surface plots of the objective functions under different slices of the
decision space. (a) Slice with respect to (I X2) while fixing T3 = za 0.5. (b) Slice with
respect to (xz, 24) while fixing D1 0.3. 22 0.8. (c) Slice with respect to (æs, 244) while fixing
Ii = 0.3, 22 = 0.1.
25

---

Published as a conference paper at ICLR 2026
F DETAILS OF PDES
F.1 BURGERS' EQUATIONS
ID Burgers Equation (Burgersid)
The PDE is given by
+ u Ug = Uzz
(56)
The domain is defined as The initial and boundary conditions are
(a, = 1-1, x (0, 1].
u(x, 0) = sin(
(57)
u( 11,1)t) t(1,t) O.
(58)
The parameter is
0.01
(59)
π
0.75
0.50
0.25
-0.25
- 0.50
0.75
0.0
11.0
-0.75
0.00
0.75
1.00
x
Figure 14: Reference solution of Burgersid using PINNacle reference data.
2D Coupled Burgers Equation (Burgers2d)
The governing equations are given by
+ uu, - uuly = u(uzz + Uyy);
(60)
1t + ut + voly v(Uzz + Uyg).
(61)
The domain is defined as The initial conditions are given by
(x, y. t) E x 11 [0.
(62)
w v) Σ Ljj=-
Σ
ais sin (2r(ix + jy)) + bij (2s jui)
(63)
[u(x, v. 0), v(z, y, 0)] = 2 w( y) +
(64)
where aij. bij, (N. 1) E R². The boundary conditions are given by
u(0, y. t) = u(L, y, t),
(65)
u( 0,t) u(x L,t):
(66)
26

---

Published as a conference paper at ICLR 2026
v00,.tt) vLL, y, t), 0,t = v(z, L).
(67)
(68)
The parameters are
L= 4, - 0.001.
(69)
30
30
15
25
15
15
10
0.
0.5
05
10
15
20
25
35
40
3.0
05
1.6
15
20
25
30
25
30
3.
14
255
125
15
400
15-
-835
10
0.6
95
10
15
20
25
24
35
1.0
05
1.6
15
2.
25
30
35
Figure 15: Reference solution of Burgers2d at timesteps 0, 1, where the left panels show u and right panels show from the PINNacle reference data.
Experiment details The model architecture consisted of 3 layers with 50 neurons, and training was performed for 50,000 epochs over 5 independent runs with a learning rate of 10-3 We em- ployed the tanh activation function, Glorot normal initialization (Glorot Bengio 2010), and the Adam optimizer.
Table 5: Sampling sizes used in our experiments, following the default configuration of PINNacle.
BurgersId
Burgers2d
Collocation points
8192
32768
Initial points
2048
8192
Boundary points
2048
8192
Test points
8192
32768
With this experimental setting, the results are presented as follows.
F.2 POISSON EQUATIONS
Poisson 2D Classic (Poisson2d-C)
The PDE is given by
Muz + Myy 0.
(70)
27

---

Published as a conference paper at ICLR 2026
Table Average of relative I errors (with standard deviation) estimated across 5 random seeds for each algorithm. We highlight the best and the second methods for Burgers' Equation.
Method
Burgersid
Burgers2d
MultiAdam
0.8526 (0. (9.1900)
0.5491 (0.0000)
LRA
0.0427 (0.0267)
0.5170 (0.0031)
ReLoBRaLo
0.0159 (0.0017)
0.5282 (0.0045)
MGDA
0.0505 (0.0737)
0.5231 (0.0054)
PCGrad
0.0218 (0.0101)
0.5228 (0.0009)
CAGrad
0.0280 (0.0097)
0.5244 (0.0017)
IMTL-G
0.3221 (0.2031)
1.2998 (1.4134)
Aligned-MTT.
0.0280 (0.0097)
0.5248 (0.0018)
ConFIG
0.0264 (0.0108)
0.5245
(0.0174)
HARMONIC
0.0164 (0.0023)
0.5439 (0.0126)
The domain is rectangle minus four circles (ar 3) E Ω is the rectangle and R; denotes four circle areas:
Qren Ri, where Pree = [-0.5, 0.5]²
R1 {x : (r 0.3) + (y 0.3) s0.1")
(71)
R2 {<,,y) (r) 0.3) (y (y-0.3) 0.3)² <..1").
(72)
R3 {(x, v): (x (y+ s0.1).
(73)
R4 = ((,,)) (r) (y + 0.3)* 5..1").
(74)
The boundary condition is
u(x,,) = 0, (r. v) E RR, u(r,9) 1, (x. y) E Nn
(75)
(76)
1.0
0.8
0.6
0.4
0.0
-02
-0.4
-02
-0.4
-0.2
0.2
0.4
Figure 16: Reference solution of Poisson2d-C using PINNacle reference data.
Poisson-Boltzmann (Helmholtz) 2D Irregular Geometry (Poisson2d-CG)
The PDE is
Vu + 22 f(x,8)
(77)
where f(x) is defined as
f(x) A. + sin()rm()) (
(78)
28

---

Published as a conference paper at ICLR 2026
The domain is (z, E Ω [-1, 1]2 with several circles removed. The circles Qzec 1 are
R1 = {(z y) (z 0.5) + (y =.
(79)
R₂ ((t,)) (r 0.4) (y + 0 =.4"),
(80)
R3 ((z,)) (r 0.2) + (y7 < 0.1%)
(81)
R4 ((x, v): (r+ 0.6) (v 0.5)² < 0.3").
(82)
The boundary conditions
u(r,,) 0.2,
(r, 02wes u(r, v) = 1, (r, E eirrel--
(83)
(84)
The parameters are
P1 L, μ1 4, 8, A 10.
(85)
1.25
0.75
1.00
0.50
0.75
ubb.y)
0.00
0.25
0025
0.00
-025
-0.75
-0.50
-0.75
-0.25
Figure 17: Reference solution of Poisson2d-CG using PINNacle reference data.
Experiment details The model architecture consisted of 3 layers with 50 neurons, and training was performed for 50,000 epochs over 5 independent runs with a learning rate of 10-3, We em- ployed the tanh activation function, Glorot normal initialization, and the Adam optimizer.
Table 7: Sampling sizes used in our experiments, following the default configuration of PINNacle.
Poisson2d-C
Poisson2d-CG
Collocation points
8192
8192
Initial points
2048
2048
Boundary points
2048
2048
Test points
8192
8192
With this experimental setting, the results are presented as follows.
F.3 HEAT EQUATIONS
2D Heat with Varying Coefficients (Heat2d-VC)
The governing equation is
U1
(a((xr, f(.,1,t).
(86)
29

---

Published as a conference paper at ICLR 2026
Table Average of relative I errors (with standard deviation) estimated across random seeds for each algorithm. We highlight the best and the second methods for Poisson equations.
Method
Poisson2d-C
Poisson2d-CC
MultiAdam
0.6733 (0.0820)
0.8104 (0.0120)
LRA
0.5712 (0.1595)
0.3154 (0.0448)
ReLoBRaLo
0.6602 (0.0221)
0.6338 (0.0259)
MGDA
0.6777 (0.0239)
0.7486 (0.0281)
PCGrad
0.6918 (0.0190)
0.1914 (0.0450)
CAGrad
0.7806 (0. 1553)
0.3615 (0.0621)
IMTL-G
0.5039 (0.0684)
1.0175 (1.1243)
Aligned Aligned-MTL
0.2847 (0.2363)
0.7024 (0.0416)
ConFIG
0.6954 (0. (2.6)
0.1252 (0. (9.1941)
HARMONIC
0.0214 (0.0179)
0.1460 (0.1810)
The domain is
(z, y) x = [0, 1]2 x [0,5].
(87)
The function a(x, v) is chosen similarly to Darcy flow but with an exponential GRF.
The function fx, f) is defined as
f(x y, A sin(minx) sin(m2ry) sin(m))t).
(88)
The initial and boundary conditions are
u(x, y, 0) =0, 0, (x, y) Ω u( = 0, (x, v) an.
(89)
(90)
The parameters are
A 200, m. 1, m2 5, m3 1.
(91)
0.75
6.155
0.050
8.0100
8.0500
-0.025
-0.000
0.0000
-0.073
2A
0.6
0.8
0.0000
0.0300
0.0000
0.2225
1005
0.2300
0.0000
0.2075
-0 6075
0.2000
0.4
0.2025
02-
0.-
0.2000
0.159
0.0003
0.
0.0
Figure 18: Reference solution of Heat2d-VC at timesteps 0,0.3, 0.7, 1 ussin PINNacle reference data.
2D Heat Multi-Scale (Heat2d-MS)
30

---

Published as a conference paper at ICLR 2026
The PDE equation is given by
W1 (500z)² Maa
x²2 Kyy
(92)
The domain is (z, E Ω T [0, 1|2 [0.
The initial and boundary conditions are
u(x, y, 0) = sin(20Tx) sin( (2x, y. t) O,
(2, E Ω
(93)
(x, v) E in, te T.
(94)
0.35
0.75
0.50
04
-.25
0.4
0.5
02
0.2
0.0
0.0
0.50
0.55
0.25
-050
0.50
0.
0.6
0.5
1.0
0.8
0.2
9.4
0.
0.8
1.0
Figure 19: Reference solution of Hear2d-MS at timesteps ( = 0, 1,25 using PINNacle reference data.
2D Heat Complex Geometry (Heat2d-CG)
The PDE is
Vu=0.
(95)
The domain defined as
( (1-8, 8, 8) x [-12, 12] U:R:) x [0.3].
(96)
The initial condition
u(I, 0) O.
(97)
The boundary conditions are
-n- (- Vu) 5 --1
(98)
(Vu) u,
(99)
n- (a-(-vww Vu) = 0.1 u.
(100)
Dimensional Heat Equation (HNd)
The governing equation is
4 kAu + f(x,t),
(101)
where
f(x,t) = | exp 1
(102)
31

---

Published as a conference paper at ICLR 2026
18-
23
1.5
-5
-10
-18
-10
-12
Figure 20: Reference solution of Heat2d-CG at timesteps = 0, 2,33 using PINNacle reference data.
The domain is (æx, 1) E Q x T. The geometric domain is a unit sphere in d-dimensional space and T = [0, 1j. We choose d as 5.
The initial and boundary conditions are
n Vu g(x, (), (x, 1) E an x [0, 1], u( ( 0) g(x, 0), x E S
(103)
(104)
where
g g(x, t) exp( ||
(105)
We choose dimension d = 5 and the parameter is
(106)
d'
Experiment details The model architecture consisted of 3 layers with 50 neurons for Heat2d-VC, Heal2d-CG and HNd, while Hear2d-MS was configured with layers and 100 neurons. Training was performed for 50,000 epochs over 5 independent runs with a learning rate of 10-3, We employed the tanh activation function, Glorot normal initialization, and the Adam optimizer.
Table 9: Sampling sizes used in our experiments, following the default configuration of PINNacle.
Heat2d-VC
Hcat2d-MS
Heat2d-CG
HNd
Collocation points
32768
32768
32768
32768
Initial points
8192
8192
8192
8192
Boundary points
8192
8192
8192
8192
Test points
32768
32768
32768
32768
With this experimental setting, the results are presented as follows.
Table 10: Average of relative L errors (with standard deviation) estimated across 5 random seeds for each algorithm. We highlight the best and the second methods for Heat.
Method
Hcat2d-VC
Hcat2d-MS
Heat2d-CG
HNd
MultiAdam
1.1038 (0.0708)
0.9987 (0.0021)
0.7049 (0.0107)
0.0322 (0.0338)
IRA
0.2984 (0.0474)
0.0278 (0.0070)
0.2930 (0.1242)
0.0020 (0.0011)
ReLoBRaLo
1.8305 (0.1208)
0.0355 (0.0212)
0.1057 (0.0034)
0.0034 (0.0067)
MGDA
1.2157 (0.1851)
0.0824 (0.1416)
0.1205 (0.0063)
0.0048 (0.0088)
PCGrad
1.1198 (0.1314)
0.0324 (0.0161)
0.0991 (0.0063)
0.0028 (0.0019)
CAGrad
0.3355 (0.0437)
0.0442 (0.0176)
0.1173 (0.0013)
0,0005 (0.0001)
IMTL-G
1.2430 (0.4677)
0.9933 (0.0027)
0.0989 (0.0037)
0.0006 (0.0002)
Aligned-MTL.
0.9080 (0.6703)
0.9849 (0.4311)
0.1189 (0.0070)
0.0006 (0.0002)
ConFIG
0.2812 (0.1329)
0.8121 (0.4211)
0.1195 (0.0107)
0.0006 (0.0001)
HARMONIC
0.1927 (0.0071)
0.0317 (0.0331)
0.0918 (0.0062)
0,0005 (0.0000)
32

---

Published as a conference paper at ICLR 2026
F.4 WAVE EQUATION
ID Wave Equation (Waveld-C)
The governing equation is
Utt Unn 0.
(107)
The domain is (a, 1) 7 [0 x [0.
The initial condition
u(2, 0) = sin( + sin a
(108)
2
(e( 0) = 0.
(109)
The boundary conditions are
u(0,t) t) = u(1,t) 0.
(110)
The analytical solution of this problem
u(x, t) = sin (TT) cos() + 2
sin(ATx)) cos( 8Tt).
(111)
1.0
0.8
0.5
0.6
005
-1.0
x
Figure 21: Reference solution of Waveld-C using PINNacle reference data.
Experiment details The model architecture consisted of 3 layers with 50 neurons. Training was performed for 50,000 epochs over independent runs with learning rate of 10-3 We employed the tanh activation function, Glorot normal initialization, and the Adam optimizer.
Table 11: Sampling sizes used in experiments, following the default configuration of PINNacle.
Waveld-C
Collocation points
8192
Initial points
2048
Boundary points
2048
Test points
8192
With this experimental setting, the results are presented as follows.
F.5 NAVIER- SHIERRSTOKKS EQUATIONS
2D NS lid-driven flow (NS2d-C)
33

---

Published as a conference paper at ICLR 2026
Table 12: Average of relative I errors (with standard deviation) estimated across 5 random seeds for each algorithm. We highlight the best and the second methods for Wave equation.
Method
Waveld-C
MultiAdam
1.1540 (0.0630)
LRA
0.2608 (0.0542)
ReLoBRaLo
0.4102 (0.0476)
MGDA
0.4627 (0.1359)
PCGrad
0.2256 (0.0701)
CAGrad
0.1186 (0.0071)
IMTL-G
0.7664 (0.4757)
Aligned-MTL
0.3801 (0.1741)
ConFIG
0.0668 (0.0279)
HARMONIC
0.0655 (0.0293)
The PDE is given by
1
u2 + olly + Px
(Uzx + Uyy),
(112)
Re
1
ut + voy + Py
(Uew + Uyu),
(113)
Re
Ur + Uyy
0.
(114)
The domain is (x, E Ω [0, 1Ja. the top boundary is IT, the left, right and bottom boundary is T. The boundary conditions are
u(x, v) 4x(1 - x),,
E Tı,
(115)
u(x, 4) 0,
y lT,
(116)
"(r, y) O, v(x, v) = O.
E..
(117)
p(0,0) 0)
(118)
The Reynolds number Re 100.
Figure 22: Reference solution of NS2d-C, where the left panel shows u, the middle pannel shows v, and the right panel shows p from the PINNacle reference data.
2D Back Step Flow (NS-CG)
The PDE is given by
1
Wu, + ruy + Ps.
(uza + Ugg),
(119)
Re
1
ut wDy + Py
(Uaw + Uyy))
(120)
Re
ur + U1
(121)
The domain is defined as (x, v) E Ω = j0, 4] x [0.2] ((0, 2] x [1,2] R,) (excluding the top-left quarter). The inlet boundary is T and the outlet boundary Tout-
34

---

Published as a conference paper at ICLR 2026
The boundary conditions
w(0, y) 4y(1 -),
(122)
v(0, y) O,
(123)
p(4, y) 0,
(124)
u(x; y) O, v(x, v) = 0,
I. E Fin UT out-
(125)
The Reynolds number of Re 100.
15
12
15
Figure 23: Reference solution of NS-CG, where the left panel shows u, the middle pannel shows v, and the right panel shows p from the PINNacle reference data.
Experiment details The model architecture consisted of layers with 100 neurons. Training was performed for 50,000 epochs over independent runs with a learning rate 10-3, We employed the lanh activation function, Glorot normal initialization, and the Adam optimizer.
Table 13:
Sampling sizes used in experiments, following the default configuration of PINNacle.
NS2d-C
NS-CG
Collocation points
8192
8192
Initial points
2048
2048
Boundary points
2048
2048
Test points
8192
8192
With this experimental setting, the results are presented as follows.
Table 14: Average of relative L errors (with standard deviation) estimated across 5 random seeds for each algorithm. We highlight the best and the second methods for Navier-Stokes equations.
Method
NS2d-C
NS-CG
MultiAdam
1.3164 (0.3269)
0.8175 (0.0260)
LRA
0.5027 (0.1197)
0.2760 (0.0129)
ReLoBRaLo
0.0954 (0.0366)
0.1212 (0.0093)
MGDA
0.4316 (0.2851)
0.1999 (0.0605)
PCGrad
0.0593 (0.0283)
0.29 (0.0171)
CAGrad
0.1372 (0.0945)
0.2027 (0.0147)
IMTL-G
0.6960 (0.2076)
0.3463 (0.0984)
Aligned-MTT
0.1980 (516)1)
0.197 (0.0424)
ConFIG
0.9274 (0. 4688)
0.2889 (0.0669)
HARMONIC
0.0767 (0.0639)
0.1547 (0.0144)
F.6 INVERSE PROBLEMS
Poisson inverse problem (PInv)
The governing PDE is
-V.(aV)) f(x,))
(126)
35

---

Published as a conference paper at ICLR 2026
The source term f is
2x sin T.. sin(Ty)
f(z 1 x² 22² 1)2 + (y-14 1)
2n (2r 1) cos (Tw)) sin(Tr) (2y 1) ( cos(Tt))
(127)
(1 + x 3² (x 1)2 - (y 1)2)2
The geometric domain is (z. E Qt = [0, 11 and
u(x, y) = sin ( sin(
(128)
To ensure the uniqueness of the solution, we impose a boundary condition of a(z, y),
a(x,y) v) 1 + (x
(xz, v) E an
(129)
We sample data of uz y) with 2500 uniformly distributed 50 x 50 points and add Gaussian noise (0 0.1) to it. The ground truth of a(x, is
a(x,y/ y) 1 + 1)2 1)-
(x, v) E Q.
(130)
Heat (Diffusion) inverse problem (HInv)
The PDE is given by
V-(a f(x
(131)
where the source function f(x, v, is
f(x,,,t ((4T 1) sin (тι)) sTy(
+ (xx) (Tσ) sin (my) cos? (πx) sin din2(n (xy) sin sin)( cos (Ty)))
(132)
The geometric domain is (x, y/, 1) E Ω x T |-1,
1 1J2 x [0, 1}. and u = sin (Tz) sin(Tr).
(133)
We impose boundary condition for the diffusion coefficient field
a(x, 9) Oz E Ω.
(134)
We sample data of u(r,, y,() I) randomly with 2500 points from the temporal domain Ω x T and add Gaussian noise (0, 0.1) to it. The ground truth is
a(x; y) + (in( sin()
(a. E Ω
(135)
Experiment details The model architecture consisted of layer with 100 neurons for Plnv and layers with 50 neurons for HInv. Training was performed for 50,000 epochs over 5 independent runs with learning rate or 10-3 We employed the tanh activation function, Glorot normal initialization, and the Adam optimizer.
Table 15: Sampling sizes used in our experiments, following the default configuration of PINNacle.
Plnv
Hlnv
Collocation points
8192
32768
Initial points
2048
8192
Boundary points
2048
8192
Test points
8192
32768
With this experimental setting, the results presented as follows.
36

---

Published as a conference paper at ICLR 2026
Table 16: Average of relative I errors (with standard deviation) estimated across 5 random seeds for each algorithm. We highlight the best and the second methods for Inverse problems.
Method
Pfnv
Hmny
MultiAdam
0.2047 (0.0724)
0.1857 (0.0652)
LRA
0.1324 (0.0273)
0.0546 (0.0131)
ReLoBRaLo
0.5176 (0.9037)
1.5999 (0.0576)
MGDA
0.2434 (0.0528)
2.4224 (1.6765)
PCGrad
0.0909 (0.0216)
0.1269 (0.0267)
CAGrad
0.2080 (0.0484)
0.6435 (1.0862)
IMTL-G
0.3775 (0.2857)
0.0515 (0.0105)
Aligned-MTL.
0.0915 (0.0194)
0.0804 (0.0457)
ConFIG
0.1344 (0.0320)
0.0525 (0.0132)
HARMONIC
0.1183 (0.0271)
0.0506 (0.0121)
F.7 VOLTERRA PROBLEMS
1D nonlinear Volterra
The IDE expressed as
du(x;)
dx + u(x)) c² u(t) dt.
feettalt
(136)
If we define
v(x) /-" c²~² u(t)
(137)
then the IDE can be written as
du(xx) dx
+ u(x) \v(x).
(138)
Also, v(z) should satisfy
dv(x)) u(x) / e e²==x(tt u(t) dt = u(r)- v(x).
(139)
dr
The domain is x E Ω = 10, 5].
The equation has an exact solution
u(x) = cosh(a).
(140)
The initial condition is
(0)) = 1, (0) =
(141)
(142)
The parameter is À 1.
Volterra system
The IDEs system is
deu(x) dx2
x3
1 do(())
+ >1 I (x²2t)++
=1
+ + dt,
(143)
3
dr
dv(xx)
-1 + +x² x² - x xu(x)) + A2 1" (12(1) (()))~2(())
=2) dt.
(144)
dx²
(145)
As we define
w(x) = 1 (c²(t) (v2(t)+ dt.
(146)
p(x) /
dt.
(147)
37

---

Published as a conference paper at ICLR 2026
0.9
0.8
0.6
Figure 24: Reference solution of Volterrald.
we can rewrite IDEs as
d du(x))
x² dv((x)
+ \Iw(x),
(148)
du²
2 dr
d=v(xx)
-1 + x² x ru(x) )2p(1))
(149)
dx
(150)
As follows, auxiliary outputs w(x) and p(x:) are governed by
dw(x) u² (x) + v² (x)) dr
(151)
dp(a:) u² (x) u² (z)
(152)
The domain is Ω = =0,, 1].
The exact solutions are
u(x) + c²
(153)
m(x) x 6"
(154)
The initial conditions are
du(
dx 2, du(0)
dr 0.
(155)
and
u(0) 1, v(0) 1, w(0) 0, p(0) 0.
(156)
The parameters are
1 2
A1
A2
(157)
Experiment details The model architecture consisted of 3 layers with 50 neurons for Volterra ld and 5 layers with 100 neurons for Volterra system. Training was performed for 50,000 epochs over independent runs with a learning rate of 10 3, We employed the tanh activation function, Glorot normal initialization, and the Adam optimizer.
With this experimental setting, the results presented as follows.
F.8
OTHER PDES
We conducted experiments the PDEs included in the PINNacle benchmark, however, did not evaluate certain PDEs that do not align with the purpose our study. First, in the two-loss setting, the concept is identical to DCGDDcenter), and even the computational procedure coincides with that
38

---

Published as a conference paper at ICLR 2026
Table 17: Sampling sizes used in our experiments.
1D Volterra
Volterra system
Collocation points
8192
8192
Initial points
2048
2048
Boundary points
2048
2048
Test points
8192
8192
Table 18: Average of relative L errors (with standard deviation) estimated across 5 random seeds for each algorithm. We highlight the best and the second methods for Volterra.
Method
Volterrald
Volterra system
MultiAdam
0.0007 (0.0004)
0.9998 (0.0000)
LRA
0.0070 (0.0065)
0.3120 (0.1856)
ReLoBRaLoo
0.0036 (0.0044)
0.0009 (0.0005)
MGDA
0,2806 (0. (2.3)
0.0119 (0.0148)
PCGrad
0,0012 (0.0010)
0.0039 (0.0084)
CAGrad
0.0014 (0.0024)
0.0059 (0.0092)
IMTL-G
1.3180 (1.3065)
0.6734 (0.1375)
Aligned-MTII
0.0071 (0.0073)
0.0268 (0.0575)
ConFIG
0.0018 (0.0024)
0.1082 (0.1274)
HARMONIC
0.0003 (0.0001)
0.0007 (0.0008
of ConFIG. Therefore, excluded Poisson 3D Complex Geometry with Two Domains (Poisson3d- CG), 2D Poisson Equation with many subdomains (Poisson2d-MS) N-Dimensional Poisson Equa- tion (PNd), and the Kuramoto- Kuramoto-Sivashinsky Sivashinsky Equation (KS), since even when accounting for the gov- erning equation and initial/boundary conditions, they involve only two losses. Second, as reported in PINNacle, the relative 7,2 error for the 2D Heat Long Time (Heat2d-LT),, 2D Davier-Stokes Long Time (NS2d-LT), 2D Wave Equation in Heterogencous Medium (Wave2d-CG), 2D Multi- Scale Long Time Wave Equation (Wave2d-MS), and the 2D Diffusion-Reaction Gray-Scott Model (GS) typically exceeds or remains close to it across most methods. For these PDEs, we con- ducted experiments with bot our method and the comparison methods, but the overall performance was unsatisfactory. Since under such circumstances it is not meaningful to claim that one method outperforms another, we chose not to include these results in the report.
G DETAILS OF COMPARISON METHODS
We aim introduce the considerations taken into account for reporting the hyperparameters of each method used in the experiments and for comparing their performance. In addition, we did perform hyperparameter tuming, since in the multi-loss setting relative superiority cannot be captured by single quantitative validation metric, making PDE-specific hyperparameter selection unavailable. Accordingly, we adopted the hyperparameter settings that were either commonly used or explicitly recommended in the original papers.
MultiAdam As Adam applies Adam to each loss in a multi-loss setting, it takes the expo- nential decay rates for the first-moment and second-moment estimates as hyperparameters. In our experiments, we followed the recommendation of MultiAdam (Yao et al. (023) by setting both hyperparameters to 0.99.
ReLoBRaLLo ReLoBRaLoo estimates the scaling factors using the following procedure.
39

---

Published as a conference paper at ICLR 2026
L,(t)
exp
TL,(t)
A hal (t, = m
i f1 f1....,,m m}
(158)
((((1)
2j=1 m
exp
(Ly(tT)
Ahin (t) = pd;(t 1) + (1 -1) Ax(t) hist (1 a)) (t.t 1)
(159)
(160)
In the above equations, m denotes the number of losses, and T. and a are hyperparameters. The hyperparameters were
0.1, P 0.9999, and 0.999, following the configuration used in the Burgers Id experiments of ReLoBRaLo (Bischof & Kraus 2025).
MGDA In MGDA, the gradient vectors are first normalized, and the update step is taken in the direction of the minimum-norm point within the convex hull formed by their normalized convex combination. We employed the l normalization, proposed MGDA (Sener & Koltun 2018). which has the advantage of providing stable definition of the update direction.
CAGrad In CAGrad (Liu et al., 2021a), the hyperparameter E [0, 1) constrains the update di- rection remain within certain distance from the average of the loss gradients, and set 0.5 in our experiments.
ConFIG In ConFIG, the update direction originally obtained by solving a least-squares problem to find a vector whose inner product with each loss gradient equals 1. However, this procedure produced NaN values all five independent runs on NS2d-C, Volterrald, and Volterra system. Therefore, in our implementation, instead of using the unstable least-squares method, we computed the pseudo-inverse the loss gradient matrix and multiplied it by 1m to obtain the update direction instead.
NTK At each iteration, the scaling factors are estimated using the singular values of the kernel matrix, where the kernel matrix is defined
K JJ
(161)
and the Jacobian matrix J has elements corresponding to the derivatives of each loss with respect to the network parameters 8. In practice, our experiments employed networks ranging from 3 layers with 50 neurons to 5 layers with 100 neurons. Thus, computing the Jacobian matrix and subse- quently forming kernel matrix at every iteration would incur a prohibitive computational cost. Moreover, since the number of losses in our experiments can be as large as 10, we considered such computations impractical. fact, even for the simplest case of the one dimensional Burgers equa- tion, evaluating the Jacobian and constructing the kernel matrix result in average computational speed of only 0.2 iterations per second, excluding the cost of training itself. By contrast, all other methods except NTK achieved training speeds exceeding 50 iterations per second average, indi- cating that the NTK approach is prohibitively inefficient under our experimental setting.
40

---

Published as a conference paper at ICLR 2026
H ADDITIONAL EXPERIMENTS
H.1 COMPARISON ACROSS LEARNING RATES
We conducted experiments on the WavelD-C, Poisson2D-C. HNd, VolterralD, and HInv problems using two learning-rate configurations: 10 10-3, and 10-4, The results for each setting are reported accordingly.
Table Average of relative T,a errors (with standard deviations) estimated across 5 random seeds for each algorithm. For every learning rate, we highlight the best and the overall best methods for Waveld-C.
Method
10 10-3
10
MultiAdam
1.1540 (0.0630)
1.0958 (0.1004)
LRA
0.2608 (0.0542)
0.3385 (0.0277)
ReLoBRaLo
0.4102 (0.0476)
0.3542 (0.0089)
MGDA
0.4627 (3.1359)
0.3995 (0.0135)
PCGrad
0.2256 (0.0701)
0.3124 (0.0198)
CAGrad
0.1186 (0.0071)
0.3313 (0.0428)
IMTL-G
0.7664 (0.4757)
0.2724 (0.3372)
Aligned-MTI.
0.3801 (0.1741)
0.2302 (0.0514)
ConFIG
0.0668 (0.0279)
0.1690 (0.0423)
HARMONIC
0.0655 (0.0293)
0.2499 (0.0456)
Table 20: Average of relative _ errors (with standard deviations) estimated across 5 random seeds for each algorithm. For every learning rate, we highlight the best and the overall best methods for Poisson2d-C.
Method
10 10-3
10
MultiAdam
0.6733 (0.0820)
1.2376 (0.1976)
LRA
0,5712 (0.1595)
0.6908 (0.0212)
ReLoBRaLo
0.6602 (0.0221)
0.6629 (0.0208)
MGDA
0.6777 (0.0239)
0.6840 (0.0192)
PCGrad
0.6918 (0.0190)
0.9286 (0.2068)
CAGrad
0.7806 (0.1553)
0.6969 (0.0197)
IMTL-G
0.5039 (0.0684)
0.9743 (0.6459)
Aligned-MTL
0.2847 (0.2363)
0.7137 (0.0324)
ConFIG
0.6954 (0.4296)
0.6856 (0.0274)
HARMONIC
0.0214 (0.0179)
0.3579 (0. (6.4681)
Table 21: Average of relative L errors (with standard deviations) estimated across 5 random seeds for each algorithm. For every learning rate, we highlight the best and the overall best methods for HNd.
Method
10 10-3
10
MultiAdam
0.0322 (0.0338)
0.0040 (0.0028)
LRA
0.0020 (0.0011)
0.0036 (0.0057)
ReLoBRaLo
0.0034 (0.0067)
0.0004 (0.0000)
MGDA
0.0048 (0.0088)
0.0009 (0.0001)
PCGrad
0.0028 (0.0019)
0.0005 (0.0002)
CAGrad
0.0005 (0.0001)
0.0006 (0.0000)
IMTL-G
0.0006 (0.0002)
0.0006 (0.0001)
MLT.
0.0006 (0.0002)
0.0016 (0.0014)
ConFIG
0.0006 (0.0001)
0.0008 (0.0008)
HARMONIC
0.0005 (0.0000)
0.0038 (0.0076)
41

---

Published as a conference paper at ICLR 2026
Table 22: Average of relative I,2 errors (with standard deviations) estimated across 5 random seeds for each algorithm. For every learning rate, we highlight the best and the overall best methods for Volterrald.
Method
10-3
10
MultiAdam
0.0007 (0.0004)
0.0004 (0.0002)
LRA
0.0070 (0.0065)
0.0080 (0.0039)
ReLoBRaLoo
0.0036 (0.0044)
0.0002 (0.0000)
MGDA
0.2806 (0. 1263)
0.0011 (0.0019)
PCGrad
0.0012 (0.0010)
0.0004 (0.0004)
CAGrad
0.0014 (0.0024)
0.0017 (0.0027)
IMTL-G
1.3180 (1.3065)
0.0735 (0.0989)
Aligned-MTL
0.0071 (0.0073)
0.0037 (0.0044)
ConFIG
0.0018 (0.0024)
0.0042 (0.0069)
HARMONIC
0.0003 (0.0001) 0.0004 (0.0001)
Table 23: Average relative L errors (with standard deviations) estimated across 5 random seeds for each algorithm. For every learning rate, we highlight the best and the overall best methods for HInv.
Method
10
10
MultiAdam
0.1857 (0.0652)
0.1446 (0.0245)
LRA
0.0546 (0.0131)
0.0626 (0.0061)
ReLoBRaLo
1.5999 (0.0576)
1.7514 (0.0447)
MGDA
2,4224 (1.6765)
2.8915 (1.6877)
PCGrad
0.1269 (0.0267)
0.0660 (0.0152)
CAGrad
0.6435 (1.0862)
2.2489 (0.4788)
IMTL-G
0.0515 (0.0105)
0.0498 (0.0075)
Aligned-MTL.
0.0804 (0.0457)
0.1902 (0.0532)
ConFIG
0.0525 (0.0132)
0.0466 (0.0068)
HARMONIC
0.0506 (0.0121)
0.0461 00.0098)
H.2 COMPARISON UNDER TWO LOSS COMPONENTS
We conducted experiments using loss configurations: the average of PDE losses, and the average of initial and boundary condition losses.
Table 24: Average relative 2 errors (with standard deviations) estimated across 5 random seeds for each algorithm.
Method
Waveld-C
Poisson2d-C
HNd
Volterrald
Hin
MultiAdam
1.0020 (0.0035)
1.3607 (0.0003)
14.0535 (1.7338)
0.8254 (0.0009)
0.1221 (0.0359)
I.R
0.3537 (0.0214)
0.4903 (5.1)
0.0010 (0. 0005)
0.0010 (0.0010)
0.0687 (0.0094)
ReLoBRaLo
0.1566 (0.0385)
0.6948 (0.0173)
0.0010 0016)
0.0025 (0.0037)
1.4010 (0.0415)
MGDA
0.4282 (0.00
0.6927 (0.0223)
0.0045 (0 (0.0030)
0.0037 (0.0013)
1.4780 (0.0726)
PCGrad
0.4065 (0.0120)
0.5727 (3.93
0.0012 (0.00 (0.0016)
0.0003 (0.0002)
0.1587 (0.0259)
CAGrad
0.3692 (0.0236)
0.6967 (0.0174)
0.0013 (0. (0.0019)
0.0008 (0.0010)
2.0905 (0.2565)
IMTL-G
0.3593 (0.0312)
0.7119 (0.0168)
0.0004 (0. 0.000)
0.0026 (0.0036)
0.0514 (0.0135)
Aligned- MTT.
014865 1849)
0.6572 (0.0664)
0.0005 (0.0000)
0.0028 (0.0036)
0.2769 (0.4351)
ConFIG
0.3979 (0.0374)
0.3991 (0.3600)
0.0005 (0. 0001))
0.0028 (0.0057)
0.0521 (0.0120)
DCGD
0.4022 (0.0309)
(0.3755)
0.0005 (0.0001)
(0.0015)
0.0515 (0.0122)
HARMONIC
0.3882 (0.0215)
0.4242 (0.3757)
0.0005 0000)
0000 (0.0001)
0.0521 (0.0119)
As can be seen in Table 24 under the two-loss setting, DCGD, ConFIG, and HARMONIC all lie within the harmonized cone, resulting in highly similar performance across the methods.
42

---

Published as a conference paper at ICLR 2026
H.3 SCALABILITY WITH MORE LOSS COMPONENTS
In addition to the problems presented in the main text( Wa Poisson2D-C, HNd, Volterra and HInv) we also evaluate the computational cost of our method NS2D-C, which involves eight loss components.
Table 25: Average computational costs (in seconds) with standard deviations from 100 epochs across 100 repeats. The worst methods are highlighted.
Method
NS2d-C
MultiAdam
5.6655 (0.1410)
LRA
6.1194 (0.1649)
ReLoBRaLo
1.2052 (0.0556)
MGDA
9.8877 (1.6511)
PCGrad
6,1200 (0.3040)
CAGrad
6.6949 (0.1894)
IMTL-G
5.5782 (0.3040)
Aligned-MTL
5.6768 (0.1596)
ConFIG
5.6790 (0.1842)
HARMONIC
5.6746 (0.1412)
As shown in Table 25 even with large number of losses, the computational cost remains compa- rable to or even lower than that other methods.
43

---

Published as a conference paper at ICLR 2026
I
LLM USAGE
Large Language Models (LLMs) were used in limited capacity aid in the preparation of this manuscript. Specifically, LLMs were employed solely for writing support tasks such as grammar checking, polishing sentence structure, and improving readability of author drafts. They were not involved in research ideation, methodological design, data analysis, or result interpretation. The authors take full responsibility for the content of this paper.
44


Reviews:
[{"review_id": "vt9aZs4uw5", "reviewer": "Reviewer_wpQu", "summary": "This paper proposes a geometric framework to address two key challenges in the multi-loss training of Physics-Informed Neural Networks: infeasible update directions caused by loss scaling and conflicting gradients across objectives. The authors introduce the harmonized cone, defined as the intersection of feasible and non-conflicting directions, and develop the HARMONIC algorithm based on this idea. Using the Double Description method, they compute valid update directions that satisfy both conditions. Theoretical results guarantee convergence to Pareto stationary points, and experiments on the PINNacle benchmark, particularly on the challenging Poisson2d-C problem, demonstrate superior robustness and stability compared to existing methods.", "questions": "Please refer to the Weaknesses section for the main points. In addition, I have the following specific questions:\n\n- Would it be possible to conduct an ablation study to isolate the effect of the infeasibility component?\n\n- The proposed method performs remarkably well on the Poisson2d-C problem, suggesting that avoiding infeasibility plays a crucial role in this particular case. Could the authors provide further explanation on why infeasibility is especially critical for this benchmark and what characteristics of Poisson2d-C make this effect more pronounced?", "limitations": "", "rating": 6, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- (Clear motivation and good presentation) The paper convincingly explains the difficulties of PINN training into two factors, infeasibility and conflict, supported by well-designed motivating examples and illustrative figures. The theoretical framework is presented in a structured manner, and the connection between geometric intuition and optimization behavior is clearly articulated. Overall, the paper is well written.\n\n- (New insight into the training of PINNs) This paper identifies a new and important issue, infeasibility, which has not been explicitly addressed by existing conflict-free methods. The paper provides a principled framework to resolve this issue by combining the idea of dual cone gradient descent methods. \n\n- (Acceptable computational overhead) Although the Double Description method used to compute the harmonized cone can be computationally expensive in theory, the number of loss terms in PINN is typically very small. Given this practical context, the computational cost remains acceptable.", "weaknesses": "- (Limited novelty): The overall novelty of this paper is somewhat limited. The proposed method heavily relies on the dual cone gradient descent (DCGD) framework and differentiates itself mainly by introducing the additional infeasibility restriction compared to methods such as ConFIG. \n\n- (Unclear effect of infeasibility restriction): The effect of the infeasibility restriction remains unclear. While the proposed harmonized cone integrates both feasibility and conflict-avoidance constraints, the paper does not provide sufficient empirical or theoretical analysis isolating the specific contribution of the infeasibility component.\n\n- (Computational cost) The experiments in this paper appear to focus primarily on cases with three or more loss terms. However, in practical PINN applications such as inverse problems, the number of losses can easily increase to four or more. In such cases, the computational cost associated with the Double Description method may grow significantly and could become a clear disadvantage (acceptable though) compared to other existing methods.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "r954Q7Nb6G", "reviewer": "Reviewer_mqgc", "summary": "This paper proposes Harmonized Cone Gradient Descent (HARMONIC) for multi-loss optimization. HARMONIC constructs the update vector by aggregating the rays of the Harmonized Cone, which is defined as the intersection of the primal gradient cone (formed by conic combinations of individual loss gradients) and the dual gradient cone. The rays of the Harmonized Cone are computed using the Double Description method. HARMONIC outperforms existing loss-balancing and gradient-manipulation methods on multiple PDE benchmarks while maintaining comparable computational cost. Furthermore, applying the harmonized-cone constraint to existing methods (including ConFIG, ReLoBRaLo) improves their performance.", "questions": "**Questions**\n1. HARMONIC appears to be applicable not only to PINNs but also to other domains with multiple losses (for example, multi-task learning). Have you applied HARMONIC to such fields or conducted preliminary experiments beyond PINNs? (This question will not affect my rating since I understand that the rebuttal period is short.)\n\n2. According to the experimental results provided in the Appendix, there are several benchmarks (e.g., Burgers, PInv, Poisson2d-CG) where HARMONIC does not outperform other competitors. In these cases, algorithms that do not satisfy the feasibility or non-conflict properties sometimes achieve better performance. Could you explain the limitations of HARMONIC and under what conditions such results are likely to occur?\n\n3. In my understanding, HARMONIC appears to be a generalized version of DCGD, since DCGD also satisfies the characteristics listed in Table 1 but only for the two-loss setting. Therefore, as in the experiments of ConFIG [1], it would be valuable to include two-loss scenarios to clarify this relationship. More importantly, demonstrating the advantages of HARMONIC when extending to three or more losses, beyond the simple separation of the PINN loss into PDE residual and other term, would provide stronger evidence for the algorithm’s contribution.\n\n[1] Liu, Q., Chu, M., & Thuerey, N. ConFIG: Towards Conflict-free Training of Physics Informed Neural Networks. In The Thirteenth International Conference on Learning Representations.", "limitations": "", "rating": 4, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "**Strengths**\n1. The method of generating and aggregating the rays of the Harmonized Cone to determine the update direction is interesting and appears novel. Moreover, by showing performance improvements when combined with existing loss-balancing methods (e.g., ReLoBRaLO) and gradient-manipulation methods (e.g., ConFIG), the proposed approach empirically validates the motivation and effectiveness of enforcing the harmonized cone.\n\n2. Although the proposed HARMONIC method is demonstrated on Physics-Informed Neural Networks (PINNs), it is widely applicable to other deep learning tasks involving multiple losses, suggesting potential contributions to various domains.\n\n3. The proposed HARMONIC method shows superior performance compared to existing approaches in the presented experimental results.", "weaknesses": "**Weakness**\n\n1. The explanation of the proposed algorithm is insufficient. The paper does not clearly describe how the Double Description method computes the rays of the Harmonized Cone, and it lacks references or related works that would help readers understand this process. It also remains unclear whether the Double Description method can always compute the rays exactly or if there are cases where approximation or numerical instability may occur. Furthermore, it would be helpful to explicitly include the Double Description step in the pseudocode to clarify how it is integrated into the overall algorithm.\n\n\n2. The experiments presented in the paper were conducted on relatively simple benchmark problems. It would strengthen the work to evaluate HARMONIC on more challenging PDE benchmarks where PINNs are known to struggle, such as the Navier–Stokes or Kuramoto–Sivashinsky equations (see [1]). Demonstrating that HARMONIC can improve PINN performance in such difficult cases would provide stronger evidence of the algorithm’s effectiveness. Furthermore, it is important to analyze whether HARMONIC remains effective when combined with various PINN variants (for example, [2], [3]), as this would highlight its robustness and general applicability.\n\n3. The example in Figure 1 is not sufficiently convincing in demonstrating the advantage of HARMONIC, as several compared algorithms such as CAGrad and Aligned-MTL also show convergence to the Pareto front. To better highlight the unique characteristics and advantages of HARMONIC, it would be more appropriate to include an example where other algorithms fail to converge while HARMONIC succeeds. Such a case would more clearly illustrate the distinctive behavior and strength of the proposed method.\n\n4. The paper proposes a new optimizer and compares it with several existing methods, but it seems that the learning rate was fixed at 1e-3 for all experiments. Since different optimizers often require different optimal learning rates, it would be fairer to compare them after tuning the learning rate individually for each optimizer.\n\n[1] Wang, S., Sankaran, S., & Perdikaris, P. (2024). Respecting causality for training physics-informed neural networks. Computer Methods in Applied Mechanics and Engineering, 421, 116813.\n\n[2] Zhao, Z., Ding, X., & Prakash, B. A. PINNsFormer: A Transformer-Based Framework For Physics-Informed Neural Networks. In The Twelfth International Conference on Learning Representations.\n\n[3] Cho, J., Nam, S., Yang, H., Yun, S. B., Hong, Y., & Park, E. (2023). Separable physics-informed neural networks. Advances in Neural Information Processing Systems, 36, 23761-23788.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Ly2boguzvE", "reviewer": "Reviewer_CSQr", "summary": "This paper focuses on improving training physics-informed neural networks (PINN). PINN is known to be challenging for training due to different losses, e.g., PDE loss, boundary loss and initial value loss. Prior works can be categorized into two main streams: 1. Reweighting the losses; 2. Pursuing multi-objective optimization (MOO). It turns out that these two directions can be combined together, which is the method proposed in this paper.\n\n1. The problem can be reformulated as seeking parameter updating direction that is both *feasible* and *non-conflict*, which is called *harmonic* direction. Feasible means nonnegative linear combination of gradients wrt each loss. Non-conflict means the direction in which none of losses increases. \n2. Algorithm 1 *HARMONIC* is proposed to find a harmonic direction, which is built upon Double Description Method (DDM).\n3. Theoretic analysis of algorithm 1 (Theorem 3) shows that convergence to Pareto front is guaranteed.\n4. Experiment on benchmark PINNacle shows that the method is highly competitive without extra computation resource.", "questions": "1. What is \"the descent lemma\" in the proof of Theorem 3 in Appendix B.4? I understand that equation (20) holds for gradient update rules. Please clarify this.\n2. Also in the proof of Theorem 3 in line 776, why does the norm of $\\mathcal{A}_h(G)$ equal to the sum? If $g_j$'s are not orthogonal, why is this true?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "**Originality** is high. This paper summarizes existing works and sharply points out the limitation of either \"feasible\" or \"non-conflict\" methods. The proposed method is built on a principled framework, see summarization.\n\n**Quality** is high. The method seems to be theoretically sound and the comprehensive experiment shows strong evidence of effectiveness.\n\n**Clarity** is high. I found the paper mostly clear to follow.\n\n**Significance** is high. The paper provides a unified and effective method that sets new SOTA to PINN methods.", "weaknesses": "See questions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)