Incorporating Sum Constraints into Multitask Gaussian Processes
Authors: Philipp Pilar
Carl Jidling
Thomas B. Schön
Niklas Wahlström
Venue: Accepted by TMLR ()

Abstract: Machine learning models can be improved by adapting them to respect existing background knowledge. In this paper we consider multitask Gaussian processes, with background knowledge in the form of constraints that require a specific sum of the outputs to be constant. This is achieved by conditioning the prior distribution on the constraint fulfillment. The approach allows for both linear and nonlinear constraints. We demonstrate that the constraints are fulfilled with high precision and that the construction can improve the overall prediction accuracy as compared to the standard Gaussian process.

Published in Transactions on Machine Learning Research (11/2022)
Incorporating Sum Constraints into Multitask Gaussian Processes
Philipp Pilar
philipp-pilardee pilar su.se
Deportment of Information Technology
Uppsala University
Carl Jidling
carl.jidling uu.se
Department of Information Technology
Uppsala University
Thomas B. Schön
Woomasschon@ii..c wu.se
Department of Information Technology
Uppsala University
Niklas Wahlström
niklas. wahlstrom Oit su.se
Deportment of Information Technology
Uppsala University
Reviewed on OpenReview: https: // openrevie.. net/ bet/ffornn forum? id= g2u42bBYY7
Abstract
Machine learning models can be improved by adapting them to respect cxisting background knowledge. this paper we consider multitask Gaussian processes, with background knowledge in the form of constraints that require a specific sum of the outputs to be constant. This achieved by conditioning the prior distribution on the constraint fulfillment. The approach allows for both linear and nonlinear constraints. We demonstrate that the constraints are fulfilled with high precision and that the construction can improve the overall prediction accuracy compared to the standard Gaussian process.
1 Introduction
Many real world problems come with background knowledge known a priori, for instance that the outputs must be positive at all times fulfill certain differential equation. The constraints are often known to near perfect precision. Any model would certainly benefit from having such knowledge hardcoded in advance instead of having to rediscover it, as the additional information would allow for the exclusion of the majority of possible outputs.
In this work we consider the Gaussian process (GP) Rasmussen Williams 2006 which is a popular and powerful machine lcarning model. Some assumptions about the underlying function, c.g. regarding its smoothness can be encoded in relatively straightforward way into the kernel of the GP. However, is usually trickier to include more specific prior knowledge and constrained GPs (or, for that matter, constrained machine learning methods) constitute relevant and active area of research (Willard et al. 2021 Swiler et al. 2020
this work, we focus on constraints that take the form of a sum over the outputs of a multitask GP. Constraints of this form arisc, for example, when considering conservod quantities in physics such as energy and momentum, where the sum over the energies or momenta of all subcomponents of closed system must remain constant. As toy example, we consider the harmonic oscillator, which is ubiquitous in physics; the

---

Published in Transactions on Machine Learning Research (11/2022)
expression for the energy takes the form
E = (t) + (xin()) /2+ + m
(1)
where Epo and Ekin denote potential and kinetic energy, respectively. We assume that the displacement from the rest position 2 and the velocity v are the outputs of multitask GP, whereas the time t serves as input. While the input this example is one-dimensional, the results we derive in this paper also apply higher dimensional inputs.
We have developed a method that allows nonlinear constraints like [ to be incorporated into the GP. First, we show how nonlinear constraints can be reduced to linear ones via suitable transformation of the outputs of the GP. Then we proceed to condition the joint prior of the GP on the constraints, which in turn results in a coistrained predictive distribution In the next section, we start by providing a formal definition of the problem.
2 Problem Formulation
2.1 Background on the GP
A GP is formally defined "a collection of random variables, any finite number of which have joint Gaussian distribution" (Rasmussen & Williams 2006). Formally, we write f(x) ~ GP(mmxx), x')), where m(x) = [x] and x(x,x')= [[((() - m m(x))] are the mean and the covariance function of the GP, respectively. The dataset available for training the GP consists of inputs X {xx J and noisy
outputs y = f(xx) + ck: where we assume Gaussian noise E ~ N((0, o)).
We use y to denote a vector storing
all N outputs.
In the following we consider the multitask setting (Bonilla al. 2008 Skolidis & Sanguinetti 2011 where vector f(x) of NN outputs is learned. The overall GP framework remains unchanged but the output vector f (and observation vector y) has to be interpreted as an extended vector consisting of the concatenated multitask outputs fa fxxx) that is f = [ [ for which it holds that ~ (X), Kr (X,X') X'))
When constructing the mean and covariance function, the different tasks need to be taken into account [Alarrz al. 2012). We write the mcan as
mF () = [ma((i) (mx)))))))))) mm(x) (xxx)))
where ma(-) is the data mean and m-(-) is the task mean. The task mean returns column vector of length Ny. The covariance matrix becomes
Kain (x1,
kaiekt (x11 kazake (x2, X2)
Kr.r (X.
kazıke (x2, x1)
(3)
where kaij = ka(X,,,)), and where ka(ss)) :) and k.(-.)) denote the data and task kernels, respectively. Note that the task kernel returns a matrix of size (N,, Ny).
The task mean and kernel often assumed to be position independent (although this assumption is not necessary for our method to work); then my and Ke can be written Kronecker products
mx(X) = ma X) Ø m
(4a)
Kr (X,X') X') = ka X, X') Σ..
(4b)
Given the expressions for the mean and the kernel, the predictive distribution is formed through the standard
procodure; see Section [.I in the supplementary material for details. See also Section [32] for dotails on how
to deal with the case of incomplete measurements, i.c. when there are data points yk for which only some of the output tasks have been measured.
2

---

Published in Transactions on Machine Learning Research (11/2022)
2.2 Sum Constraint
The main concern of this work is to show how constraints on the sum of some (nonlinear) transformations h,(-) of the outputs fi can be incorporated into the GP. Formally, we define this class of constraints as
F(f(x)] = E (x)xx(())(( C(x).
where the functions (x) serve prefactors to the various terms in the sum, i indexes the outputs of the GP, and C(x) specifies what value the sum over the outputs should equal at position x in the input space. In the following we refer to constraints of this form sum constraint.
In the general case 5 we consider input-dependent constraints C(x) and (x)).. This requires knowledge of the functions ((x) and a,(x), which could be practically infeasible. Hence, an important special case of 5 is the constant sum constraint
F[f(x)] Σ axha(fi(x)) -,
with constant prefactors a, and constant sum C.
One example constant sum constraint is the previously mentioned energy conservation for the harmonic oscillator ① There we have 01 k/2, 02 m/2, (2) = 2² h2((v) v² and C E. Other situations where sum constraints arise include learning of probabilities that must sum one, and the of mechanical equilibrium where the sum of acting forces must be zero each point.
3 Method
Let now develop the methodology required to incorporate constraints defined in Section 2.2 into the GP. In Section 3.1.1 consider the case where all the outputs of the GP enter the sum constraint via a monotonic (invertible) nonlinearity and show how to reduce it to linear sum constraint. In Section 3. 1.2 we extend the procedure to sum constraints with non-monotonie nonlincarities. Finally, we show in Section 3.2 how to include linear sum constraints into the GP and hence, via the aforementioned reductions, also nonlinear sum constraints.
3.1
Reduction to Linear Constraint
3.1.1
Monotonically Increasing Nonlinearity
while the constraint is nonlinear in terms of the outputs, it is linear in terms of the transformed outputs hu(fi); defining Si ha(fs) and substituting it into (5) yields
Consider the sum constraint 5
F[f((x)) = Σ (x) = C(x),
which is linear in the transformed outputs H. Hence, we can train GP to predict the transformed outputs obeying the linear constraint 7 and backtransform to the original outputs via f. = -f() Note that this GP needs to be trained on transformed data y', where y/ = hu(y:). This approach requires that the nonlinear functions h.).) are invertible, otherwise it is not possible to unambiguously recover the ft- See also Snelson et al. 22044..
However, is not necessary for h; to be invertible on its entire domain. Consider the case where it is known that the output fo is restricted to invertible subregion of the domain of h2 then we solve the problem by choosing the backtransformation h-1 corresponding to this subregion. For example, in case of the square function, we can consider the case where f; is known to be always positive (or always negative). Then we can just restrict the domain of the nonlinearity h; to the positive (negative) half-axis where the function is fact. invertible.
When employing the transformation 7 it is important keep in mind that the GP prior now has to be chosen in way suitable for the transformed outputs instead of f; depending the transformations h(.)
3

---

Published in Transactions on Machine Learning Research (11/2022)
Algorithm 1 The Constrained GP: High-level Procedure
Step 1: train an unconstrained GP on the data y obtain the auxiliary outputs faux (optional) use the posterior mean of faux to create virtual mcasurements Step 2: train the constrained GP on the transformed data y' (for details, see Algorithm 2 to obtain f'
(ooptiona) (re)learn the auxiliary outputs together with the constrained outputs
Step 3: backtransform the transformed outputs f' using the posterior mean of faux from Step 1
involved, this could prove to be more challenging. We recover credible intervals for f in the same way as we recover f, by backtransforming them; for more details, see Section B.7 in the Supplementary material.
Furthermore, the noise corresponding to the transformed data y' will in general not be normally distributed anymore, which means that GP regression loses its analytical tractability due to the resulting non- con-Gaussian likelihood. Methods to deal with non-Gaussian likelihoods include the Laplace approximation (Williams & Barber [998] Vanhatalo et al. 2009), variational inference (Blei et al. 2017: Tran et al., 2016). and expectation propagation (Minka) (0011). Due to its simplicity, in this work we use the Laplace approximation to deal with this issue, where applicable. It enables to approximate non- non-Gaussian distributions with Gaussian; Appendix B.3 for details.
3.1.2 Non-monotonically Increasing Nonlinearity
In the previous section we showed how to reduce nonlinear sum constraints to linear ones, as long as the nonlinearities are monotonic. However, this a rather limiting assumption as it would exclude e.g. the square function h(f) f² from the admissible transformations. Here we describe a way of circumventing this problem.
The idea underlying our solution is to introduce one (or multiple) auxiliary variables that allow for a unique backtransformatton. Typically, the auxiliary variables will keep track of where in the domain of k(-) it is that I' lies, such that the correct local inverse can be chosen when backtransforming. In case of the square function, we can add the auxiliary output faux = f and retrieve the initial output f via f = () = sign(faux) While the initial output f is a practical choice here, this is in general not necessary and faux can be chosen arbitrarily.
There is no guarantee that learned values f' will always fall within the domain of the backtransformation. If it happens that a predicted value lies outside, a pragmatic solution is to approximate f' with the closest valid value; for example zero in case of negative valued predictions for square values.
Somctimes more information can be extracted from faux and used to ameliorate the transformed data y', for instance when the backtransformation switches from one local inverse to another; then we can add virtual measurements for f' at those points and force the constrained GP towards values consistent. with fauxs which can significantly reduce artefacts in the backtransformed outputs f. Note that this can come at the cost of overconfident credible intervals in the vicinity of the virtual mcasurements.
In Algorithm 1 we summarize this procedure. In most cases, it is advantageous to learn the auxiliary outputs in a separate GP in Step 1, independently of the constrained outputs; when virtual measurements are to be created, this is required. Optionally, auxiliary outputs can be (re)learned in Step 2; for some examples, this can stabilize the hyperparameter learning of the constrained GP. However, when virtual measurements are involved, the prediction faux from Step should also be used for the backtransformattin..
We illustrate the approach by returning to the harmonic oscillator [ with the transformed outputs fi = 2² and f v2 (see also the last paragraph in Section 2.2 We choose the auxiliary outputs faux and f2 = v, which we use to extract the sign when backtransforming fl and fí; furthermore, we use the auxiliary outputs to create virtual measurements for Si and S at the zero crossings of the posterior mean of f and fax. In order to fit the transformed outputs of the GP, the observations yą = [x]" are transformed analogously obtain y's = [22 vZ. Ex, vz] T; 2 and U/ are part of y' since we chose relearn them together with the constrained outputs to improve the performance. The virtual measurements are also included in the transformed data y'. In terms of the transformed outputs the const can be written compactly
4

---

Published in Transactions on Machine Learning Research (11/2022)
Algorithm Constraining the GP (Section [.11 refers to the Supplementary material)
Input: mean mr mr(-); kernel (--;( constraints (F,, S); (transformed) data y,y"; points of prediction X.
Output: constrained predictive distributiin f, ,', X.
Note: During hyperparameter optimization X, = [ and hence f. {
Step 1:
Construct the joint prior distribution for [f', f.]" (No, Ev) according
B.1
-omit noise term of Step 2: Construct Foot Stot according to 9b
Step 3: Use Frot, Stot to calculate constrained μ E' according to 8b
Step 4: Remove entries in μt, E' corresponding to incomplete measurements as detailed in Section B.2 if Hyperparameter optimization then
Step 5: Calculate the log marginal likelihood according to B.7cl
Step 6: Perform optimization step
else if Prediction then
Step 5: Calculate the predictive distribution f' X, according to B.7a
end if
Ff' C, where F = [01, 02 0, 0, For more details on the harmonic oscillator dataset, see Section C.1 in the supplementary material.
3.2 Solving with Linear Constraints
Having shown how to reduce nonlinear sum constraints to linear ones, proceed to describe how to incorporate linear sum constraints into the GP. The idea is to make use of the fact that sampling from a GP is equivalent to sampling from a multivariate Gaussian distribution, where the mean and covariance are obtained by evaluating the mean and the kernel of the GP at the points of interest.
Let the random vector f' (n, Σ); we are interested in the conditional distribution F' F' =ff/= More generally, include multiple sum constraints, we want find the distribution l' S, where the rows of the matrix F contain the coefficients for each of the NF constraints to be included, and the elements of the vector S contain the corresponding sums; compare equation (9a) below.
The required conditional distribution can be calculated analytically (Majumdar & Majumdar 2019 as
(f' = S) ~ NNW, E'),
(8a)
where
n' = A E' = ATA, D DFEE" A DFF.
(8b)
Of course, need to enforce the constraint at all Noo data points that end, construct the blockdiagonal matrix Fiot and the vector Stot according
a1 (x)
(x(x)
Cx(x)
F(x)= =
(x)
b2(x)
S(x)
Cb(x)
(9a)
F Foot = diag FF (x11), F(x2).....
Siot
(9b)
We use Nost in two different contexts: during the hyperparameter optimization, Niat denotes the number of data points; whereas during prediction, Nos denotes the number of both data and predictive points.
Algorithm 2 summarizes the practical procedure of constructing the covariance and the mean, both during hyperparameter optimization and when forming the constrained predictive distribution of the GP. In case of position dependent constraint, important to note that the values of the functions C(x) and (x(x) must be known at all points for which the constraint should be enforced; in our case, this means all No points.
5

---

Published in Transactions on Machine Learning Research (11/2022)
Note that in Step 1 of the algorithm we first omit the noise term, since the constraints only hold exactly for noiseless data; the noise then enters in Step 4, after the constraints have been taken into account.
Mathematically, the constraint is enforced by conditioning the Gaussian distribution on it. While the method is not strictly global in the sense of providing a constrained kernel for the GP, is global for practical purposes as the constraint is enforced at all points of prediction of the GP.
Due to the matrix inversion in (8b), the computational complexity of the algorithm cubie with leading order term O(Ng.Nio)) during both hyperparameter optimization and prediction.
3.2.1 Special Case of Constant Constraints
In the special case of constant constraints and constant inter-task dependencies of the GP mean and kernel, the constraints can be incorporated more efficiently. Here, the kernel of the GP factorizes into data and task kernel as in (4) and the procedure above simplifies: it now suffices enforce the constraints (F, ()) on the task mean and covariance matrix and to subsequently perform the Kronecker product with the data mean and covariance matrix to obtain the constrained distribution.
Formally, this can be written as follows: let be and Ee be the task mean and covariance matrix, respectively; then the constrained quantities 14 and Σ are calculated via (8b, using F and (since the task mean and covariance matrix are constrained directly, it is not necessary to construct F Foo and Stor). Finally, the full constrained mean and covariance matrix are constructed via μ m @ 14 and E' K @ En, where m and K are the data mean and covariance matrix, respectively. Due to the constant constraint, the data mean is also required to be constant. Without loss of generality, we choose it as m 1N.... (compare B.6.1 This procedure is summarized in Algorithm 3 in the Supplementary material. Furthermore, we provide proof that this approach is indeed equivalent to the more general approach from Section 3.2 in Appendix B.6.1
Now, the complexity of the matrix inversion involved in 8b is reduced to ~ O(N)); since 44 and Σ, are constrained directly no longer depends on No (compare also 4 This constitutes a significant improvement over the general algorithm as usually Np I N₁ < Neot where Nr the number of constraints and N, the number of tasks. Whenever applicable, is preferable to use this way of incorporating the constraint, since it is more efficient and numerically more stable than the general procedure given in Algorithm
4 Experimental Results
In this section, we demonstrate our method at the hand of two simulation experiments and one real data experiment | They have in common that the constraints involved are constant (see Section 3.2.1 ) for examples of the non-constant case, see Sections A.2 and A.3 in the Supplementary material.
4.1 Toy Problem Revisited
We gave a formulation of the auxiliary variables approach for the harmonic oscillator in Section 3.1.2 and detailed information on the datasot can be found in Section C.I the Supplementary matcrial. Figure illustrates this approach. The constrained GP achieves higher overall accuracy around extremal points, where the prediction is more robust with regard to the influence of random noise. In addition, the constrained GP manages to mitigate the negative effect of incomplete measurements, i.e. data points where only one of the two output dimensions has been measured, better than the unconstrained one (compare left part of Vaux in the figure). This is natural, since the constrained GP has implicitly added a correlation between the two outputs, which the unconstrained GP is lacking.
The credible intervals in Figure I clarify another advantage of the constrained GP: when multiple outputs are learned a different degree of certainty, information can be transferred from high- to cow-credibility outputs, thereby narrowing the credible intervals also for the latter. This is elcarly visible in areas with incomplete measurements. On the other hand, credible intervals lend to be overconfident in the vicinity of
The code used for the experiments is available at httee////ehhhaacco/////////a/nnn

---

Published in Transactions on Machine Learning Research (11/2022)
Zaux
V
2
V2
2E
V
1.0
1.
155
3.-
-05-
-45-
-10
-10
-15
-15
Figure 1: Demonstration of the auxiliary variables approach for the harmonic oscillator. The quantities Raux, and Caux; v refer to the position and velocity of the harmonic oscillator, respectively. We distinguish between Zaux: and Vaux. V to emphasize that, while they aim to approximate the same curve, they are learned by different GPs. The posterior means of the GPs are depicted, together with the 2a credible intervals. The dotted lines represent the true curves and the big dots/crosses correspond to the data available to the GPs. Left: Results for the unconstrained GP are shown. For this example, these outputs coincide with the auxiliary outputs required for the constrained GP. Middle: The transformed outputs learned by the constrained GP are depicted, together with the constraint 2E kz mu² (where k = = 1). The results for the auxiliary outputs have been employed to create virtual measurements zero crossings (differently colored squares) in order to force the quadratic functions towards zero. Right: The backtransformed outputs of the constrained GP are shown, where the auxiliary outputs Zaux and Vaux have been used to recover the signs.
0.05
On = 0.1
0.3
GP-c
GP-u
GP-c
GP-u
GP-c
GP-1
RMSE
2.3+0.6
3.2
4.4
6.411.1
13.7+3.7
18.5#3.3
(o-2)
fa = 0
ACI
0.0+0.0
3.140.77
0.0+0.0
6.541..4
0.1
18.744.4
(e-2)
RMSE
3.443.0
4.822.8
5.4+1.9
8.0.2.5
17.0
23.0
(e-2)
fa = 0.2
ACI
0.0+0.1
4.3+1.5
0.1+0.2
7.8+222
0.2+0.4
22.6+6.33
(c-2)
Table 1: Comparison of the performance of the constrained GP (GP-c) and the unconstrained GP (GP-u) for the harmonic oscillator. Shown the rool mean squared error (RMSE) of the prediction as well as the mean absolute violation of the constraint, |AC]. The standard deviation of the noise given by o, whereas fa is the probability with which output components have been omitted at random from the data. The values have been obtained by averaging over 50 datasets and are given plus-or-me or-minus one standard deviation. Bold font highlights best performance.
virtual measurements. Due to the nonlinear, piecewise backtrransfrmmatio, some discontinuities have been introduced in the credible intervals of the constrained GP near the zero crossings.
In Table values for both the root mean squared error (RMSE) and the average absolute violation of the constraint ACI are given for various noise levels on, both with complete and incomplete measurements; case of incomplete measurements, the output components have been omitted at random with probability 0.2. The values have been obtained by averaging over 50 datasets. We observe that the constrained GP fulfills the constraint with to two orders of magnitude higher accuracy and also performs slightly better in terms of RMSE.
The reason why the constraint is not fulfilled with yet higher accuracy for the constrained GP is that around zero crossings it can occur that invalid values are predicted by the constrained GP (that is, negative values for z² and v²), which we pragmatically put to zero. This is also the origin of the small artefacts visible in that region of the mean curves in the right plot of Fig.

---

Published in Transactions on Machine Learning Research (11/2022)
On
GP-c
GP-u
GP-tr
RMSE
3.3+:.2
4.8+0..3
14430
(o-3)
12
le-4
|AC|
0.3
1.940.1
1.842.7
(e-3)
10-
a =
RMSE
5.5+1.0
5.010.3
9.2416
(o-3)
1e-3
Q =4
|AC]
0.8
2.0
1.611.3
(e-3)
y
RMSE
4.2:0.88
1.640.2
3.9
(c-2)
le-2
Ct 0
|AC
6.241..0
6.411.44
8.6:1.4
(e-3)
-10
10
Table 2: Results for the length constraint applied to the triangle in the plane. We compare results for the constrained GP (GP-c), the unconstrained GP (GP-u) and the unconstrained GP trained on the transformed outputs (GP-tr) (Salznann & Urtasun 2010a For small values of noise On, the sum constraint improves the performance of the GP. The values have been ob- tained by averaging over 50 datasets and are given plus-or-minus one standard deviation.
Figure 2: Visualization of the triangle in the plane. The task for the GP is give the location of the corners of the triangle (blue dots) when given the parameter a, which parameterizes different poses of the triangle.
4.2 Pose Estimation
Here we demonstrate how our approach can incorporate length constraints (Perriollat et al. 2011), inspired by applications such as pose estimation. In essence, the length constraint states that the distance Lim between two adjacent points (indexed by I and m) in a rigid body constant, irrespective position and orientation of the body. When the position is given terms of Cartesian coordinates 2;, the length constraint takes the following form
2 21
2mn + 2m Lim
(10)
This constraint is no longer an instance of the sum constraint defined in (5), since the middle term depends multiple outputs. However, with more elaborate transformation procedure, the sum constraint can still be applied.
make this more concrete, we cousider the example of triangle in the plane. Here, the outputs of interest. are the coordinates of the triangle corners, f [Ziz; Fly; R2r; 22ys R3n; 23y]- The input is continuous parametrization of different poses of the triangle in the plane. Although a is one-dimensional in this example, the approach generalizes to higher dimensional inputs. In our choice of transformed outputs, we follow the approach by Salzmann & Urtasun (2010a), where pairwise products of the original outputs are learned and subsequently transformed back via a singular value decomposition (SVD); for more details on the technicalities we refer to Section C.5 in the Supplementary material.
A visualization of the problem provided in Figure 2 where different poses a of the triangle are depicted; the blue points represent the corners of the triangle, the positions of which are learned by the GP. As can be seen from the data in Table 2 our approach here performs best for low noise levels. When the noise is very small, 2 1e-3, the constrained approach achieves about the same overall accuracy in terms of RMSE as the unconstrained GP, whereas the error in the constraint is reduced by factors of 2-6.
This reduction is not simply a result of the particular parameterization of the problem, which enforces the constraint implicitly for noiseless observations, as shown by Salzmann Urtasun (2010a). To see that, we included the results for a GP that is trained on the transformed outputs, but where the constraint is not enforced explicitly. Table shows that the result is improved when enforcing the constraint in addition using the transformed outputs.
8

---

Published in Transactions on Machine Learning Research (11/2022)
7.5
Up
3.-
Ekin
E
25-
-2
Epat
00-
-50-
x
0.00
1.00
2.00
t
Figure 3: Left: Trajectory of the double pendulum. Note that the trajectory shown here is longer than the sequences of motion considered in the plots to the right. Middle: Kinctic energy Ekin: potential energy Epot and total energy E of the double pendulum are shown. It is apparent that for the considered segment of the motion the energy is constant for practical purposes, except for fluctuations in the contribution of the kinctic energy due to measurement error. An estimate É of the energy is obtained by averaging over E. Right: Positions 20, and velocities vo, V of the masses (four components each) as learned by the unconstrained (left inset) and the constrained GP (right inset), respectively; the posterior means of the GPs are depicted together with the 20 credible intervals. The dotted lines represent the available data, where the subset of big dots has been used for training.
4.3 Real Data Experiment: Double Pendulum
In this section we consider the 'Double Pendulum Chaotic' dataset (Asseman et al. 2018); this dataset. consists of 21 different two dimensional trajectories double pendulum and contains annotated positions of the masses attached the ends of the two pendula. Each trajectory consists of about 17000 measurements, taken a frequency of 500 Hz. For more information the parameters of the double pendulum, see Section D.I in the Supplementary material. We attempt to construct GP that models both positions Az; 2y and velocities Uar By of the two masses (i.e. 8 outputs), while at the same time respecting the law of energy conservation; the time t serves as input. As friction is present, we consider a limited section of the trajectory during the second half of the motion where we can assume constant energy (compare Figure 3 energy conservation here takes the form
E mbg2by + gyg2gy + (vī + eByy)
mg 2
(u)
(11)
where denotes the gravitational acceleration on earth, and where the indices b and g refer to the blue and the green pendulum, respectively. The constraint is incorporated into the GP in analogy to the harmonic oscillator. In terms of (6). we identify 01 0, 02 m,g; a3 a4 mgg, as m2/2, ag m az m,/2 /2, ag 2, h2(2by) Tby, ha(2gy) Zgy: hs ha(()e) uh.. ho(Uby) vby. h7((gz) Ugun hs(Ugu) 2 and C = F;; note that the coefficients a1, correspond to the outputs 2hr, Egz: which are not part of the constraint (11..
We pick a sequence of 200 data points (which are fairly close together) from one of the trajectories; 15 of these points are used during hyperparameter optimization, and to receive an estimate É of the energy. The remaining 185 points are used test data to compare the performance of constrained and unconstrained GP, both in terms of constraint fulfilmmen and in terms of RMSE with respect to the data.
Results for one individual sequence are shown in the rightmost plot of Figure 3 We observe that the constrained GP is better at learning the precise shapes of the extrema of the velocity curves, although some artefacts arise close to zero crossings due to inaccurately learned square values. For values close to zero, the credible intervals of the unconstrained GP are often smoother and thinner than those of the constrained GP.
Averaging the results over 50 sequences chosen random from the second half of the trajectories (with less friction), the RMSE for the constrained GP is 0.31 + 0.14, whereas for the unconstrained GP is 0.33 + 0.16.

---

Published in Transactions on Machine Learning Research (11/2022)
In terms of constraint fulfillmen,, the constrained GP clearly performs better with |AC| = 0.17 0.14 compared to |AC| 0.91 + 0.61 for the unconstrained GP. The values here are given plus-or-minus one standard deviation.
5 Related Work
Several research projects have considered incorporating constraints into the GP: examples include boundary
conditions (Solin & Kok 2019), inequality constraints (Veiga & A.Marrel 2012 Maatouk & Bay 2017) and
differential equation constraints (Jidling et al. 2017; Raissi et al.) 2017 2018). The recent review by Swiler jt al. 2020 provides good overview of the existing literature on constrained GPs. So far, most of the efforts have been concentrated on the single-task GP. The sum constraint, however, is qualitatively very different from constraints on single-lask GPs, in that it explicilly enforces relationship between different outputs instead of acting on individual outputs. Hence, in this section, we focus on works that consider constraints on the outputs of multitask GPs.
Prior knowledge about vector fields have been imposed into GPs through special divergence-free and curl-free kernels [Wahlström et al., 2013). Jidling et al. (2017) developed a more general method to include linear operator constraints into the kernel of the GP: this possible by using the property that GPs are closed under linear transformations (Papoulis & Pillai 2001) and relating the GP to a suitable latent GP, resembling the of potential functions in physics. See also Lange-Hegermanni (2018) for discussion of this approach from more mathematical perspective. Practical applications include modelling of electromagnetic fields
(Solin et al., 2018) and reconstruction of strain fields (Jidling et al. 2018 Hendriks et al.) 2019b 2020b
Geist & Trimpe (2020) consider affine constraints on the dynamics of mechanical systems and construct a GP satisfying Gauss' principle of least constraint.
There is connection between our method and the method by Jidling et al. (211)): while they do not consider affine constraints, their approach can be extended to include those in the context of the constant linear sum constraint (compare also Hendriks et al. (2020a), where the same idea is applicd to neural networks). The two works attack the problem from different angles: whereas Jidling et al. (2017) start by directly constructing a covariance matrix out of vectors spanning the nullspace of the constraining operator, we start with the covariance matrix and subsequently constrain it. More details on these parallels are given in Appendix F An advantage of our approach that is straightforward to inchude additional structure in the task kernel, such as in B.12). Furthermore, we consider the general case of nonconstant, nonlinear sum constraints.
Constructing kernels that are invariant with respect certain symmetry transformations has proven fruitful in the fields of atomic and molecular physics. Glielmo al. 2017 consider GPs to model interatomic force fields; they construct 'covariant kernel' by including symmetries of the force, such rotation and reflection. Methods for constructing invariant kernels are given by Haasdonk & Burkhardt (2007). : whereas Chmiela ct al. (2200) use a similar approach to construct a kernel that allows for simultancous prediction of energies and forces in molecules.
Pose estimation constitutes another area where constrained multitask GPs are of importance; in the case of rigid pose estimation, the lengths are required to be constant. A method to explicitly enforce the ronstraints during inference is given by Salzmann & Urtasun (2010b), whereas Salzmann & Urtasun (2010a) propose method to implicitly enforce the fixed-length constraint by learning transformed outputs in which the constraint is linear. We followed this latter approach in the example with the rotated triangle in Section 4.2 in addition to using the transformed outputs we also imposed the length constraint explicitly, which (at least in principle) should allow for training points that do not fulfill the constraint exactly.
6 Conclusions and Future Work
We have derived a way of incorporating both linear and nonlinear sum constraints into multitask GPs. This is achieved by learning transformed outputs and by conditioning the prior distribution of the GP on the constraint. The toy problem of the harmonic oscillator demonstrated the potential of the method; it showed that the constraint. is fulfilled with high accuracy and that the constrained GP can mitigate detrimental
10

---

Published in Transactions on Machine Learning Research (11/2022)
effects of noise or of incomplete measurements. Our experiment with the triangle in the plane showed that the sum constraint improved the method by Salamann & Urtasun (2010a) of including the length constraint into pose estimation problems; so far, these results are particularly promising in the low-noise setting. The results for the double pendulum datasct showed that our method also works well in case of real-world, noisy data, given a way of estimating the constraint with sufficient accuracy.
In light of the results received for the triangle in the plane in Section 1.2 it appears if it would be worth investigating the applicability of this approach to pose estimation problems further; especially, in cases where the approach by (Salzaan & Urtasun (2010a) gives good results, our constrained GP could potentially improve the performance. To incrcase the suitability of the approach for big datasets, combining the sum constraint framework with methods such as sparse variational inference (Hensman et al.) 2013) appears to be a fruitful direction of inquiry. Finding general methods to incorporate constraints similar to the length constraint 10) into the GP, where nonlincarities may depend on more than one of the outputs at once, constitutes another interesting avenue of future research and would widen the range of possible applications.
Acknowledgements
The work is financially supported by the Swedish Research Council (VR) via the project Plysics-inffomme machine learning (registration number: 2021-04321) and by the Kjell och Märto Beijer Foundation.
References
Mauricio A. Alvarez Lorenzo Rosasco, and Neil D. Lawrence. Kernels for vector-valued functions: a review. arXiv: 106.6251, 2012.
A. Asseman, T. Kormuta, and A. Ozcan. Learning beyond simulated physics. In Newral Information Processing Systems, Modeling and Decision-making in the Spatiotemporal Domain Workshop, 2018.
David M. Blei, Alp Kucukelbir, and Jon D. McAuliffe. Variational inference: A review for statisticians. Journal of the American Statistical Association, 112(518):859 877, 2017.
E. V. Bonilla, K. Ming, A. Chai, and C. Williams. Multi-task Gaussian process prediction. Tn Advances in Neural Information Processing Systems 20, pp. 153 160, 2008.
S. Chmiela, H. F. Sauceda, A. Tkatchenko, and K. R. Müller. Accurate Molecular Dynamics Enabled by Efficient Physically Constrained Machine Learning Approaches. Springer International Publishing, 2020.
J. R. Gardner, G. Pleiss, D. Bindel, K. Q. Weinberger, and A. G. Wilson, Gpytorch: Blackbox matrix-matrix Gaussian process inference with GPU acceleration. In Advances in Newral Information Processing Systems 31, pp. 7576 7586, 2018.
A. R. Geist and S. Trimpe. Learning constrained dynamics with Gauss principle adhering Gaussian processes. CoRR, abs/2004.11238 2020.
A. Glielmo, P. Sollich, and A. De Vita. Accurate interatomic force fields via machine learning with covariant kernels. Physical Review B, 95 (21), 2017.
B. Haasdonk and H. Burkhardt. Invariant kernel functions for pattern analysis and machine learning. Machine Learning, 68:35-61, 2007.
J. Hendriks, Alex Gregg, Chris Wensrich, and Adrian Wills. Implementation of traction constraints Bragg-edge neutron transmission strain tomography. Strain, 55 (5), 2019a.
J. N. Hendriks, C. M. Wensrich, A. Wills, V. Luzin, and A. W. T Gregg. Robust inference of two-dimensional strain fields from diffraction-basee measurements. Nuclear Instrumends and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, 444:80-99 90, 2019b.
J. N. Hendriks, C. Jidling, A. Wills, and T. B. Schön. Linearly constrained neural networks. ar Xr::v:2002.0000 2020a.
11

---

Published in Transactions on Machine Learning Research (11/2022)
J. N. Hendriks, C. M. Wensrich, and A. Wills. A Bayesian approach to triaxial strain tomography from high-energy x-ray diffraction. Strain, 56 (3), 2020b.
James Hensman, Nicolò Fusi, and Neil D. Lawrence. Gaussian processes for big data. In Proceedings of the Twenty-Ninth Conference on Uncertainty Artificial Intelligence, pp. 282-290, 2013.
C. Jidling, N. Wahlström, A. Wills, and T. B. Schön. Linearly constrained Gaussian processes. In Advances in Neural Information Processing Systems 31, pp. 1216 1224, 2017.
C. Jidling, Hendriks, N. Wahlström, A. Gregg, T.B. Schön, C. Wensrich, and A. Wills. Probabilistic modelling and reconstruction of strain. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, 436:141-156 2018
M. Hegermaan. Algorithmic linearly constrained Gaussian processes. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
A. Lindholm, N. Wahlström, F. Lindsten, and T. B. Schön. Machine Learning A First Course for Engineers and Scientists. Cambridge University Press, 2021.
Maatouk and X. Bay. Gaussian process emuulators for computer experiments with inequality constraints. Mathematical Geosciences, 49:557-558 552, 2017.
D. J. C. MacKay, Introduction to Gaussian processes. In C. M. Bishop (ed.), Neural Networks and Machine Learning, pp. 133 165. Springer, 1998.
R. Majumdar and S. Majumdar. On the conditional distribution of a multivariate normal given a transformation the linear case. Heliyon, 5(c01136), 2019.
A. G. G. Matthews, J. Hron, R. E. Turner, and Z. Ghahramani. Sample-then-opimize posterior sampling for Bayesian linear models. In NewrIPS Workshop on Advances in Approzimate Bayesion Inference, 2017.
Thomas Minka. A family of algorilloss for approcimate Bayesian inference. PhD thesis, Masachusetis Institute of Technology, 2001.
A. Papoulis and S. U. Pillai. Probability, Random Variables and Stochastic Processes. McGraw-Hill Education, New York, 2001.
M. Perriollat, R. Hartley, and A. Bartoli. Monocular template-based based reconstruction of inextensible surfaces. International Journal of Computer Vision2, 95:124 137, 2011.
M. Raissi, P. Perdikaris, and G. E. Karniadakis. Machine learning of linear differential equations using Gaussian processes. Journal of Computational Physics, 348:683-699 693,, 2017.
M. Raissi, P. Perdikaris, and G. E. Karniadakis. Numerical Gaussian processes for Lime-dependent and nonlinear partial differential equations. SIAM Journal on Scientific Computing, A0(1):A172-A1S:) 2018.
C. E. Rasmussen and C. K. T. Williams. Gaussian Processes for Machine Learning. MIT Press, 2006.
M. Salzmann and R. Urtasun. Implicitly constrained Gaussian process regression for monocular non-rigid pose estimation. In Advances in Neural Information Processing Systems 23, pp. 2063-2073 2010a.
M. Salzmann and R. Urtasun. Combining discriminative and generative methods for 3D deformable surface and articulated pose reconstruction. IEEE Conference on Computer Vision and Pottern Recognition, 2010b.
G. Skolidis and G. Sanguinetti. Bayesian multitask classification with Gaussian process priors. IEEE Transactions on Neural Nelworks, 22(12), 2011.
E. Snelson, C.E. Rasmussen, and Z. Ghahramani. Warped Gaussian processes. In Advances in Neural Information Processing Systems 16, pp. 337-344, 2004.
12

---

Published in Transactions on Machine Learning Research (11/2022)
A. Solin and M. Kok. Know your boundaries: Constraining Gaussian processes by variational harmonic features. In Proceedings of Machine Learning Research2, volume 89, pp- 2193 2193-2202, 2019.
A. Solin, M. Kok, N. Wahlström, T. B. Schön, and S. Särkkä. Modeling and interpolation of the ambient magnetic field by Gaussian processes. IEEE Transactions Robolics, 34 (4):1112-1127, 2018.
Laura P. Swiler, Mamikon Gulian, Ari Frankel, Cosmin Safta, and John D. Jakeman. survey of constrained Gaussian process regression: Approaches and implementation challenges. Journal of Machine Learning for Modeling and Computing, 1:119-156, 2020.
Michalis Titsias and Neil D. Lawrence. Bayesian Gaussian process latent variable model. In Yee Whye Teh and Mike Titteringlon (eds.), Proceedings of the Thirleenth International Conference on Arlificial Intelligence and Statistics, volume 9 of Proceedings of Machine Learning Research, pp. 844-854 04-151, 2010.
Dustin Tran, Rajesh Ranganath, and David M. Blei. The variational Gaussian process. In 4th International Conference on Learning Representations, 2016.
Jarno Vanhatalo, Pasi Jylänki, and Aki Vehtari. Gaussian process regression with Student-t likelihood. In Advances in Neural Information Processing Systems, volume 22, 2009.
S. Da Veiga and M.rerel Gaussian process modeling with inequality constraints, In In Annales de la Faculté des Sciences de Toulouse, volume 21, pp. 529 555, 2012.
N. Wahlström, M. Kok, T. B. Schön, and F. Gustafsson. Modeling magnetic fields using Gaussian processes. IEEE International Conference on Acoustics, Speech and Signal Processin, pp. 3522 3526, 2013.
Jared Willard, Xiaowei Jia, Shaoming Xu, Michael Steinbach, and Vipin Kumar. Integrating scientific knowledge with machine learning for engineering and environmental systems. arXtv:203.049199, 2021.
C.K.I. Williams and D. Barber. Bayesian classification with Gaussian processes. IEEE Trunsactions on Pottern Analysis and Machine Intelligence, 20(12):1342 1351, 1998. doi: 10.1109/34.735507.
13

---

Published in Transactions on Machine Learning Research (11/2022)
Z
Vaux
Figure 4: Demonstration of the auxiliary variables approach for the free fall. The quantities 2o, and Vauxs refer to the position and velocity of the mass, respectively. We distinguish between 20, and Vauxs V to emphasize that, while they aim to approximate the same curve, they are learned by different GPs. The posterior means of the GPs are depicted, together with the 2-sigma credible intervals. The dotted lines represent the true curves and the big dots/crosses correspond to the data available to the GPs. Left: Results for the unconstrained GP are shown. For this example, the u-curve coincides with the auxiliary variable Vaux required for the constrained GP. Middle: The transformed outputs learned by the constrained GP are depicted. The results for the auxiliary variable have been employed to create a virtual measurement for v at the zero crossing of Vaux (differently colored square) in order to force the quadratic function towards zero. Right: The backtransformed outputs of the constrained GP are shown, where the auxiliary output Vaux has been used to recover the sign of v.
0.05
On = 0.1
= 0.3
GP-c
GP-u
GP-c
GP-u
GP-c
GP-u
RMSE
1.9::..5
2.340.4
3.1
4.7+1.22
10.113.6
13.243.3
(e-2)
fa 0
|A|
0.7+2.6
35.3+13.6
0.1
79.0=31.8
0.0+0.1
216.1+84.0
(c-2)
RMSE
2.5+0.9
3.411.2
3.941.6
6.0+2.0
15.317.9
20.0
(e-2)
fa = 0.3
AC|
0.4+0.8
47.5++66.8
0.140.1
93.7+35.1
0.7+2.1
286.9+123.8
(o-2)
Table 3: Comparison of the performance of the constrained GP (GP-c) and the unconstrained GP (GP-u) for the free fall. Shown the root mean squared error (RMSE) of the prediction well the mean absolute violation of the constraint, |A||. The standard deviation of the noise is given by o, whereas fa is the probability with which output components have been omitted at random from the data. The values have been obtained by averaging over 50 datascts and are given mins-or-minus one standard deviation. Bold font highlights best performance.
A Additional Examples
In this Section, take look at additional examples where the constraint be applied. The free fall dataset in Section A.I is another example from physics, where one of the outputs enters linearly into the constraint, instead of quadratically. The damped harmonic oscillator in Section A.2 constitutes a variation of the harmonic oscillator toy example and demonstrates the case of a non-constant constraint. In we invostigate an example where the constraint includes differont nonlincaritics.
Scction A.3
A.1 Free Fall
In addition the harmonic oscillator (see Section 4.1 we investigated the simple example of a mass free fall as second toy problem. Here, the output of the GP consists in position and velocity of the mass,
14

---

Published in Transactions on Machine Learning Research (11/2022)
f [], whereas the time serves as imput. Then the constraint takes the following form
F(f(t)] mgz(t) + m v(t)² = Epot (t) + Ekin (t) = E.
(A.1)
In terms of (G), we identify 01 mg, 02 - m h1 (2) = z; 2, h2(v)= - x² and C- E. Hence, we receive for the transformed outputs fi = fi = 2 and f/2 =². We choose the auxiliary output flux = u, which we use to extract the sign when backtransforming the posterior mcan of faux- In order to fit the transformed outputs of the GP, the observations Yą = [ are transformed analogously to obtain y/2 = [ze, v ux]", tg is part of the constrained outpuls, as this improves the performance for this example. The virtual measurements are also included in the transformed data y'. In terms of the transformed outputs the constraint can be written compactly as Ff' = C, where F = [01, a2 0]. For more details on the free fall dataset, see Scction
f/2 and to create virtual measurements for $2 at the sero crossings of
In Figure 4 results for both constrained and unconstrained GP, applied to the free fall dataset, are depicted. When comparing the left and the right plot, it is apparent, that the constrained GP manages to mitigate detrimental effects of both noise and incomplete measurements, where some of the observed output components have been omitted at random, better than the unconstrained GP (compare area around peak of æn in the figure). When looking at the 2o credible intervals, we get the same picture as before for the harmonic oscillator: the constrained GP can utilize higher cortainty in one output and transfer it to the other one, resulting in overall slimmer intervals. Close to the zero crossing of v, however, some artefacts are present due to the piecewise, nonlinear backtransformation, which are absent for the unconstrained GP; furthermore, confidence intervals are stretched bit due to the backtransformation via the square root.
In Table 3 values for both the root mean squared error (RMSE) and the average absolute violation of the constraint ACI are given for various noise levels on, both with complete and incomplete measurements; in case of incomplete measurements, the output components have been omitted at random with probability fa = 0.3. The values have been obtained by averaging over 50 datasets. We observe that the constrained GP fulfills the constraint with up to two orders of magnitude higher accuracy and also performs better in terms of RMSE.
A.2 Damped Harmonic Oscillator
Next, we investigate slight variation of the harmonic oscillator, the damped harmonic oscillator. The formal treatment remains mostly unchanged and details can be found in Section 3.1.2 in the main paper; the main difference is that damping has been added to the model of the oscillator. As a consequence, the energy is no longer constant and the amplitude of the oscillation decays over time; see Section C.2 for more details. Hence, this example constitutes an instance of the non-constant sum constraint F(f(t)) = E(I), where Algorithm applies.
In Figure 5 results for both constrained and unconstrained GP are depicted. The findings are similar the undamped harmonic oscillator, and it is apparent that the constr GP can mitigate the detrimental effects of noisy or incomplete measurements better than the unconstrained GP. In Table 4 the performance 50 random datasets is evaluated. The outputs of the constrained GP fulfill the constraint with up to two orders of magnitude higher accuracy than the unconstrained one, and also perform slightly better in terms of RMSE This example demonstrates that, given similar datasets, the performance of our method is very similar, both in case of constant and non-constant constraints (compare Section 4.1
A.3 Non-square Nonlinearity
Finally, we investigate an example where nonlinearities other than the square-nonlinaariyy are involved in the constraint. We consider the outputs f [f1, f2] on which we want to enforce the constraint
F[[(x)] = log (f (x)) (2(=))(()))=))))))))))))())))))) C(x).
(A.2)
In terms of 6 we identify a1 1, a2 1, h1 (1(1) = log(f1), he(fe) sin(f2) and C C(x). Here, we assume that the true value C(x) is known. Hence, receive for the transformed outputs f log(fi) and f/ sin(f2)-
15

---

Published in Transactions on Machine Learning Research (11/2022)
Zaux
Vaux
2E
1.50
05
0.75
0.
0.00
Figure 5: Demonstration of the auxiliary variables approach for the damped harmonic oscillator. The quantities Saux, and Vaux: V refer to the position and velocity of the damped harmonic oscillator, respectively. We distinguish between Raux, 2 and Vaux, to emphasize that, while they aim to approximate the same curve, they are learned by different GPs. The posterior means of the GPs are depicted, together with the 2o credible intervals. The dotted lines represent the true curves and the big dots/crosses crosses correspond to the data available to the GPs. Left: Results for the unconstrained GP are shown. For this example, these outputs coincide with the auxiliary outputs required for the constrained GP. Middle: The transformed outputs learned by the constrained GP are depicted, together with the constraint 2E(t) k2² + mw² (where and m = 1). The results for the auxiliary outputs have been employed to create virtual measurements at zero crossings (differently colored squares) in order to force the quadratic functions towards zero. Right: The backtransformed outputs the constrained GP are shown, where the auxiliary outputs 2aux and Vaux have been used to recover the signs.
0.05
On = 0.1
0.3
GP-c
GP-u
GP-c
GP-u
GP-c
GP-u
RMSE
3.141.3
3.2+0.6
5.64
6.5+1.2
13.4
17.74+3.8
(e-2)
fa
ACI
0.0+0.0
2.440.6
0.0.0.0
5.111.1
0.110.1
13.343.3
(e-2)
RMSE
4.1.2.33
4.8.2.22
6.2.1..6
7.81211
20.7+11.4
23.5
(e-2)
fa = 0.2
ACI
0.110.1
3.3
0.1
5.8=1.4
0.2+0.3
16.344.8
(e-2)
Table 4: Comparison of the performance of the constrained GP (GP-c) and the unconstrained GP (GP-u) for the damped harmonic oscillator. Shown are the root mean squared error (RMSE) of the prediction as well as the mean absolute violation of the constraint, |AC]. The standard deviation of the noise is given by o, whereas fa is the probability with which output components have been omitted at random from the data. The values have been obtained by averaging over 50 datasets and are given plus-or-minus standard deviation. Bold font highlights best performance.
We choose the auxiliary output faux - S2, which we to disambiguate the backtranssormattion via the aresine, that is keep track of how many multiples of = /2 the output fax has crossed. We also use the auxiliary output to create virtual measurements for f at points where the posterior mean of fl crosses multiples of f in order to reduce artefacts caused by the discontinuity in the backtransformation
In Figure 6 results for both constrained and unconstrained GP are depicted. It apparent that, while not perfect, the constrained GP outperforms the unconstrained one. In Table 5 the results averagod over datasets are given. We see, that the constrained GP outperforms the unconstrained one in terms of RMSE, and it fulfills the constraint with up to 30 times higher accuracy.
16

---

Published in Transactions on Machine Learning Research (11/2022)
log(f1)
Sin(f2)
2C
-10
-0. 2.0
2.0 -2.0 -05 0.0
2.3 -10 -0.5
8.0
6.5
20
Figure Demonstration of the auxiliary variables approach for the example with non-square nonlincarity. The quantities SP, and f2 fa refer the same respective outpuls of the GPs. We distinguish between f fi and f2, fo to emphasize that, while they aim to approximate the same curve, they are learned by different GPs. The posterior means of the GPs are depicted, together with the 2-sigma credible intervals. The dotted lines represent the true curves and the big dots/ dots/crosses crosses correspond to the data available to the GPs. Left: Results for the unconstrained GP are shown. For this example, the fl--urve coincides with the auxiliary output faux required for the constrained GP. Middle: The transformed outputs learned by the constrained GP are depicted, together with the constraint C log(fi) + sin(f2). The result for the auxiliary output has been employed to create virtual measurement at the point where Saux crosses -x/2 (differently colored square). Right: The backtransformed outputs of the constrained GP are shown, where the auxiliary output. faux has been used to disambiguate the backtransformation via the arcsine.
0.05
0.1
On 0.15
GP-c
GP-u
GP-c
GP-u
GP-c
GP-u
RMSE
3.0
3.5
4.6+L..
7.311.1
7.0+1.0
10.842.2
(e-2)
fa =
ACI
1.342.3
5.941.4
0.9
13.9+4.8
0.5+0.33
18.0+6.3
(e-2)
RMSE
3.141.0
8.6+99.9
5.641.88
11.7+
8.4+2.1
16.2+7.3
(e-2)
fa 0.2
|AC|
1.140.8
12.2+14.2
0.740.6
17.2411.4
0.4
23.5+100.8
(e-2)
Table Comparison of the performance of the constrained GP (GP-c) and the unconstrained GP (GP-u) for the example with non-square nonlinearity. Shown are the root mean squared error (RMSE) of the prediction as well as the mean absolute violation of the constraint, |AC]. The standard deviation of the noise given by on whereas fa the probability with which output components have been omitted at random from the data. The values have been obtained by averaging over 50 datasets and are given plus-or-minus one standard deviation. Bold font highlights best performance.
A.4 Comparison of approximation methods
In this section, we give a brief comparison of different methods of approximate inference at the example of the harmonic oscillator. The approximation methods under consideration are the Laplace approximation B.3 and variational inference B.4
Fig. shows the predictive performance of the the unconstrained GP, variational inference and the Laplace approximation. While it is apparent that both approximate GPs fulfill the constraint with high precision, the variational approach tends to overfit to the data. On the other hand, overconfident credible intervals seem to be less of an issue for the variational approach than for the Laplace approximation.
In Table [ results obtained when averaging over 20 runs are given for different noise settings. is apparent that the constrained GP with Laplace approximation performs best. While the constrained GP utilizing variational inference performs worst in terms of root-mean-qquar error, the constraint still fulfilled with high precision. In case of the variational approach, it might be possible to improve upon these results by
17

---

Published in Transactions on Machine Learning Research (11/2022)
Zaux Vaux
15
M
65-
-05
-05
-05
-10
-30
Figure 7: Comparison of the performance of the unconstrained GP (Left), the variational approach (Middle), and the Laplace approximation (Right) for the example of the harmonic oscillator. The quantities Raux; a and aaxx refer the position and velocity of the harmonic oscillator, respectively. The posterior means of the GPs depicted, together with the 2o credible intervals. The dotted lines represent the true curves and the big dots/ crosses correspond to the data available to the GPs.
= 0.05
on = 0.1
GP-c L
GP-c
GP-u
GP-c L
GP-c
GP-u
RMSE
2.4
3.840.6
3.240.4
4.211.11
7.0+1.6
6.340.9
(e-2)
La = 0
|AC|
0.0
0.0+0.0
3.340.7
0.0
0.1
6.411.4
(e-2)
RMSE
3.0=1.2
5.543.0
4.541.6
6.744.0
9.1
8.322.
(e-2)
fa = 0.2
AC||
0.0:0.00
0.2+0..3
4.240.9
0.1+0.2
0.4+0.6
8.1+1.7
(c-2)
Table 6: Comparison of the performance of the constrained GP with Laplace approximation (GP-c L), the constrained GP using variational inference (GP-c var), and the unconstrained GP (GP-u) for the harmonic oscillator. Shown the root violation of the constraint, AA]]. The standard deviation of the noise is given by on whereas fa is the probability with which output components have been omitted at random from the data. The values have been obtained by averaging over 20 datasets and are given plus-or-minus one standard deviation. Bold font. highlights best performance.
squared error (RMSE) of the prediction well the
absolute
trying different, parameterizations of the variational distribution, by finding better suited optimization scheme.
B Technicalities
B.1 Background on GP Regression
In this section give very brief overview of some important GP regression formulas. For more detailed accounts see c.g. Rasmusson & Williams (2006); Lindholm et al. (2021). Given the mean function m(-) and kernel ((---) of the GP, the predictive distribution of the GP can be calculated by first constructing the joint distribution between observations y and function values at test locations f.,
y F.
K + 2I K, K.
). K..
W
(B.1)
m.
where mX),, m. mX.), k= K(X,)) X), K. K(X,X.) and K., K(X,,X,) X.).
18

---

Published in Transactions on Machine Learning Research (11/2022)
Then, the conditional distribution f. X. y, X. is constructed as follows:
f. X y. X, N (f., cov(f.)) where
(B.2a)
f, - Eff. X y. X.]
= m, + K! K + - - m).
(B.2b)
cov(f.) = K., K K + o -K, K,
(B.2c)
The log-marginal likelihood, which is used for hyperparameter optimization, is given by logp(y/
1 j - m) "() I1 m)
2 log/ K + I
1
" log2x.
(B.3)
2
B.2 Accommodating incomplete measurements
Throughout the paper, we often consider the case of incomplete measurements, iee. data points where measurements are available only for a subset of the tasks. This can be taken into account by considering equation (B.1 and removing the the rows and columns on the right-hand side corresponding to missing entries in y.
To make this more concrete, let us assumee that the j-th entry of is missing on the left-hand side of B.1 Then we also delete the j-th row of m, (K I and K., well as the j-th column of (K + o and K before explicitly constructing the joint distribution. We proceed analogously when training the constrained GP on the transformed data y'.
B.3 Laplace approximation
The Laplace approximation can be employed when the noise distribution corresponding to the (transformed) observations y' is non-Gaussian in order to obtain analytical expressions for the predictive equations and for the log-marginal likelihood. Following again Rasmussen & Williams (2006), we approximate the posterior
y(f'ly') x p(y' (f)(()))pf)) via yff) qf' l') N f f. (K + W) where
f' K (vr pog (f)) -=
(B.4a)
W - Vr Vr ply f') (")
(B.4b)
Newton's method is employed to iteratively determine f from B.4al via the update rule
p new
f' - r (Vr V('))
(B.5)
7 + (1 7)f' + (KK WV- logply + w(f' -m)),
(B.6)
where ((')
log ((y log (f and where y is the step size. In terms of these quantities, the B.2 from the previous section become
expressions
F. m m K.K 1f -m).
(B.7a)
cov (f.) = K., K. K + w-
(B.7b)
logp(y' X) (f' m) K 1(p' 1K-1(pp m) + logp(y' P') [ | W)).
(B.7c)
For details on the derivation of these formulas, see Section 3.4 in Rasmussen & Williams 2006
For the Laplace approximation, the likelihood Py (y [f') of the transformed data y' h(y) has to be known.
We assume the original data y to be contaminated by Gaussian noise. In case of the square nonlinearity
19

---

Published in Transactions on Machine Learning Research (11/2022)
where y' y'= y², the likelihood is then given by the pdf of noncentral chi-squared distribution. In the case of arbitray nonlinearities hj, the likelihood can be obtained via
dh dn
Py (y 1f") Pv, Pvt, (1ylffl) Puts (171)))))))
(B.8)
dy
where the indices and j denote data points and tasks, respectively.
There is no guarantee that Newton's method will determine the correct maximum f' in case of multimodal distributions, that the resulting Gaussian distribution will constitute good approximation of the true posterior. For these reasons, it has to be decided a case by case basis whether the Laplace approximation should be employed or not. Visual inspection of the GP predictions often gives good idea on whether the Laplace approximation performs well or not. In cases where it does not perform well, standard GP regression might still produce reasonable results. Alternatively, different methods such as variational inference (?) or expectation propagation (Minka 200T) could be employed; the equations (B.7) will then need to be replaced by expressions corresponding these techniques. A brief discussion on variational inference is given in Appendix B.4 well as comparison with the Laplace approximation in Appendix A.4
Throughout the paper, we used the Laplace approximation for the harmonic oscillator the damped harmonic oscillator A.2 and the example with non-square nonlinearity
the free fall A.1
A.3 In case of the
double pendulum 4.3 and the triangle in the plane 4.2 we chose standard GP regression over the Laplace approximation.
B.4 Variational Inference
As an alternative to the Laplace approximation (see previous section), variational inference (Bei et al. 2017 Titsias & Lawrence 2010 can be employed to approximate the posterior p(f'ly') x p(y'l "ly''l')p(("" when the likelihood p(y' is non-Gaussian The idea is to approximate the posterior with the variational distribution (('') ~ N (f'l Eg)) and to learn the parameters Hg: Eg by minimizing the Kullback-Leibler divergence between variational distribution and posterior. In order to ensure positive definiteness, the entries of the covariance matrix Eg are not learned directly, but instead the entries of its Cholesky factor Lq where it holds that Eg L Since exact minimization of the KL divergence is intractable, the evidence lower bound ELBO = log p(y') (T((((")) y')) S log p(y') is maximized in its stead.
The predictive equations in terms of the variational parameters are given by
F m K.K (My -m). cov (f) K,, K (EqK-IT I)K, K., and the ELBO can be rewritten in terms of numerically tractable, one-dimensional integrals
(B.9a)
(B.9b)
ELBO = E, -LXG)///())))-KKX)))))))])))))))
(B.10)
/ (og(uly/1r")) (UApo////))aaauurKK-KL-k))/)))))
Σ / -
Equations B.9 are obtained analogously to [B.7 for derviation of (B.10). see Titsias & Lawrence 2010
The parameters Mg. Eg of the variational distribution and the GP hyperparameters are determined jointly by maximizing the ELBO, which we do by employing gradient descent. Same as for the Laplace approximation, the likelihood for the transformed data p(y' is required when calculating the ELBO and can be obtained via B.Sh.
In our experiments, the variational approach tended overfit the data more than the Laplace approximation. It is possible that different parameterization of the variational covariance matrix, different optimization
20

---

Published in Transactions on Machine Learning Research (11/2022)
Algorithm 3 Constraining the GP Special Case of Constant Task Interdependencies
Input: data mean ma(-) 1; data kernel ka(:,") +); task mean μe; task covariance matrix Σ;;
constraints S): (transformod) data y'y'; points of prediction X.
Output: constrained predictive distribbution f. f ,',,, X.
Note: During hyperparameter optimization X, = {} and hence f" =0
Step 1: Use F, S to calculate constrained μ4, Σ according to (8b)
Step 2:
Construct parameters m, K of the (single (ask) joint prior distribution according to B.1 -omit noise term I
Step 3: Use p4, My, m, K to construct constrained (multi task) Pío Eno according 4
Step 4: Remove entries in plate Lior correppodiing to incomplete measurements as detailed in Section B.2 if Hyperparameter optimization then
Step 5: Calculate the log marginal likelihood according to B.7c
Step 6: Perform optimization step
else if Prediction then
Step 5: Calculate the predictive distribution f[ y' , X. according to B.7a
end if
scheme would manage to yield better results. In the paper, we went with the Laplace approximation over variational inference; for a brief comparison of the two approaches at the example of the harmonic oscillator, see Appendix A.4
B.5 Kernel and Mean
Throughout the paper we use radial basis function (RBF) kernel (also: squared exponential kernel) as data kernel,
KRBF ARBF(x,X') (x,x') oeexp (- ||x -
x'133
(B.11)
2/2
where of is scale factor and l is the length scale. We use the position independent index kernel provided by gpytorch (Gardner et al. 2018),
kg = BBT T + diag(v).
(B.12)
where B is a low-rank matrix and v is a non-negative vector; we chose the rank of B to be equal to the number of tasks the GP in question. The parameters oy, I, B and v are to be learned during the training process. For more examples of possible kernels see c.g. Rasmussen & Williams (2006): MacKay 1998) The
Gram matrix then constructed via the Kronecker product
Kpp X') KRBF (X, X') 8 kt.
(B.13)
We chose constant mean functions for all outputs of the multitask GP. All the models have been implemented in python with the library gpytorch (Gardner et al. 2018).
B.6 Special Case of Constant Constraints
In Section 3.2 in the main paper, detailed the method for incorporating the sum constraint into the GP in the general, non-constant case. Subsequently, in Section 3.2.1 we pointed out the possibility of implementing the sum constraint in a more efficient way for the case, where all of the constraints are constant and where the kernel of the GP factorizes into data and task kernel, as [4] and (B.13). The main ideas are discussed in the main paper, here we summarize the modified procodure in Algorithm A proof that the factorization holds also for the constrained GP is given in the next section.
B.6.1 Proof of factorization
In this section we provide formal proof that the claims made in Section 3.2.1 hold, i.c. that directly constraining the mean and task covariance matrix and subsequently performing the Kronecker product with the data mean and covariance matrix does indeed lead to the constrained GP from Section 3.2
21

---

Published in Transactions on Machine Learning Research (11/2022)
To start, let summarize the objects involved:
Eoo K @ Σt
(C.1)
Prot = m @ 14t
F Foo IN.. @
Soot 1 N X S
Here, K and Σ, denote data and task covariance matrix, whereas m and μr denote data and task mean, respectively. F and define the constraint single point. IN... denotes identity matrix and 1 vector of only ones of dimension Noa. The quantities with the tot subscript give the quantities that correspond to the general approach from Section 3.2.
Due to the requirement of constant constraint and inter-task dependencies, we also need to pick a constant data mean m, with entries a = const. We introduce the new quantity S' S which is used when constraining the task mean and covariance matrix. With equation we find the following:
Diot (Fror F Flot) 1
"1 FF not EJ
(C.2)
(IINun F)(K E, (INir ()T) (IN F) & E.)
(K @ (FE,FT) (K) @ FE!)
(K-1 @ (FE,F") - (K @ FF.)
K @ (FE, F FF,
FE, T
IN & D
A N & IN,
D F Ftot
(C.3)
fot
IN 8 IN: (INiot (INu F)
INSt & IN, IN & D
I N & (IN, D=F)
N 8 A
plo Prot + D Stot
(C.4)
(INemt @ A pe) (IN.. & (1Nre 8S)
m & Ap4 + N & D
m (Ape + DTS
m & μ4
E A Diot AT
(C.5)
(INoa & A)(K E.)(INo A)T
-K (A
K & E
Hence we have shown that pho and Eier of the constrained GP factorize into Kronecker products between the data mean and covariance matrix and the constrained task mean and covariance matrix, respectively.
22

---

Published in Transactions on Machine Learning Research (11/2022)
B.7 Credible Intervals
While standard deviation and variance for the backtransformed outputs f cannot be recovered via a simple backtransformation of the corresponding quantities of f', due the potentially nonlinear and piecewise backtransformation, it is possible to recover credible intervals in this way: for the transformed outputs f", we generate the upper and lower bounds of the 2o credible interval; subsequently those bounds can be backtransformed in the same way we do for the mean of the GP. That means that the posterior of the constrained f can be a bit skewed, i.e. the mean may not lie exactly in the middle between upper and lower credible interval. When auxiliary outputs are involved in the backtransformation, their respective means should be used also when recovering the credible intervals, for the results be consistent with the posterior means.
B.8 Training Procedure
The models have been trained using the Adam optimizer provided by gpytorch. For each experiment, the corresponding learning rate (Ir), number of iterations (iter) and (il applicable) scheduler settings are given in Table 7 The scheduler multiplies the learning rate with s-factor after s-steps iterations. The two different scheduler parameters given for the double pendulum correspond to the constrainod and the unconstrained GP, respectively.
During the training of all datasets, we checked for errors in the Cholesky decomposition, which can happen when matrix becomes singular due to numerical errors; when that happened, hyperparameter optimization was restarted with of the constrained GP proved to be less stable than for the other datasets. To counteract the issue, we tested for two further failure modes of the GP. First, we checked the learned lengthscale of the GP; if it was unreasonably small (smaller than 0.1), the training was repeated. Very small lengthscales typically correspond to the case where the GP learns an almost constant function with spikes towards all of the training points. Sccondly, we confirmed that gradient descent had actually converged during training: to this end, we took the loss values over the last 40 iterations and checked, whether the standard deviation smaller than 0.1. If either of the two checks failed, the training was repeated with newly initialized hyperparameters.
random initialization. For the non-square nonlinearity (logsin) experiment, training
iter
s-steps
s-factor
HO (Sec.
0.1
200
100
0.5
Triangle
Sec.
4.2
0.1
2000
800
0.2
DP (Sec.
4.3
0.1
2000
800/500
0.2/0.5
Free fall
A.1
0.1
200
100
0.5
Damped
HO
A.2
0.1
200
100
0.5
Non-square
A.3
0.1
200
100
0.5
Table 7: Training paramcters for the experiments.
B.9 Computing Power Available for the Experiments
All experiments have been conducted on a system with NVIDIA GTX 1060, 6GB GPU, Intel Core i7 7700-K @ 4.2GHz CPU and 16GB RAM. The creation average values for one set of parameters as displayed in typically took between 15 minutes and three hours.
Tables 4
C Details on Simulated Datasets
In this section we provide information on how the data used in the different simulation experiments was generated.
23

---

Published in Transactions on Machine Learning Research (11/2022)
C.1 Harmonic Oscillator
The data for the harmonic oscillator toy problem generated from
(t) 2g sin((w)))
(C.1a)
(t) 2gw/ cos (wot).
(C.1b)
The energy is given by
E 2 =z(t)²+ + m =tt)²=
(C.2)
= 2 0
2 sin" (wntt) + m
2 z cos (wot)
k _2-
2 2
(C.3)
We have chosen E = 0,8 m = 1 kg, wo 1s-1 and it holds that kt = mwj and = /2E/k
Training data has been generated by evaluating the function on the equally spaced grid E linspace [s]. Subsequently, random noise -N(0, was added the data and output components were omitted random with probability fai the values for on and fa are given in Table I in the main text. Test data has been generated on the grid E Vinspace(-0.1110,100 [s].
C.2 Damped Harmonic Oscillator
The data for the damped harmonic oscillator was generated from
(t) = (t) sin(wit,
(C.4a)
(t) = =((()) cos(wt) 2t( 2m
sin(wt),
(C.4b)
where =b{(t) = exp
-bt 2m
and w wug b
)² 2m
The energy is given by
Ett
k 2t)
m 2 u(t)²
(C.5)
2
which is now time dependent and no longer yields a constant expression.
We have chosen E J, m 1 kg, wa 1s b 0.1 0gs and it holds that k = mwg and 20 V2E/k.
Training data has been generated by evaluating the function on the equally spaced grid f E 1mspace(0,100,00 [s]. Subsequently, random noise ~ "0 was added to the data and output components were omitted at random with probability fa; the values for on and fa are given in Table 4 Test data has been generated on the grid e linspace(-O.1,10.11
[s].
C.3
Free Fall
The data for the free fall was generated from
(t) = Not ,
(C.6a)
v(l) = oo gl.
(C.6b)
The energy is given by
E mgz
m utt) = 171
2 Un.
(C.7)
We have chosen E 200JJ 1k and it holds that Un /2E/m, and the gravitational acceleration earth is m Training data has been generated by evaluating the function on the equally spaced grid t E linspacc(0.6,20) [s]. Subsequently, random noisc ~ (0) was added to the data and oulput components were omitied at random with probability fai the values for о, and fa are given in Table To ensure good visibility and learnability, we scaled the data y with a factor a = 20: + y/a. Both in Figure and in Table 3 the results are given terms of the rescaled data (and noise values); the results terms of the original scale can be obtained by multiplying with Test data has been generated on the grid E linspace(-0.1.. 1,6,100)) [s).
24

---

Published in Transactions on Machine Learning Research (11/2022)
C.4 Son-square Nonlinearity
The data for the experiment with non-square nonlinearities was generated from
f1(2) 2e 1( + + 0.2,
(C.8)
x²
fz(2)
(C.9)
Training data has been generated on the equally spaced grid linspace(-1...2220.. Subsequently, random noise ~ (0, of) was added to the data and output components were omitted at random with probability fa; the values for and fa are given in Table 5 Test data has been generated on the grid & E linspace(-1.2,2,1000
C.5 Triangle in the Plane
In terms of the parameter the trajectory that we used for the triangle in the plane in Section 4.2 is given
by
4 8
8.4
Zo 4 4
(C.10a)
6
Z, Zo + d(a))
(C.10b)
Z = ( + d(()).
(C.10c)
where each column of the matrix contains the coordinates of one corner point of the triangle, and where d(a) = 1/2 c2s(2a) and R(a) a rolation matrix. Subsequently, random noise ~ ~N(0, was added to Z; the values for on are given in Table 2 the main text. We then added an auxiliary point of known position (4, 4) to each datapoint, which will be important for the backtransformation:
Z =
21z
23z
(C.11)
Z1y
22y
23g
4
Following the approach from Salzmann & Urtasun (2010a). we constructed the matrix Q = ZZ and used the upper triangular elements of Q as transformed outputs for the constrained GP:
y' [Q11, Q12, Q13. Q14. Q22, Q23, Q21, Q33, Q34 Qua).
(C.12)
Then the matrix F and the corresponding vector S encoding the length constraints for all the edges of the triangle become
2
0
0
1
0
0
0
2
0
0
0
0
F
(C.13)
0
0
0
0
1
2
0
1
0
0
0
0
0
0
0
0
0
0
0
S =
L2
L
L2
L3
T
(C.14)
where Lyy denote the distances between the points and j. The the last row of F corresponds to the constraint on the distance between the auxiliary point and the origin the coordinate system. Note, that Q14, Q24 and Q34 could in principle be learned separately from the remaining transformed outputs, as they do not enter into any of the constraints and the corresponding columns in (C.13) are zero. Furthermore, Qua could be omitted from the learning process entirely, as the value is known.
After training the constrained GP, the predicted values f' are rearranged into the (symmetric) matrix Q, analogously to (C.12). Then the matrix Z is recovered via singular value decomposition (SVD) of Q This decomposition is not unique and the auxiliary point comes into play: we compare the learned with the known position and determine the angle between them, which enables us to rotate the learned coordinates to their true positions.
Training data was generated on the grid a E [0, 5), consisting of 20 uniformly spacod points. Subsoquently, random noise e (O( of) was added the data; the values for on are given in the main text. Test data generated over the same range (0, 5), although this time with the grid divided into 100 points.
25

---

Published in Transactions on Machine Learning Research (11/2022)
400
500 Hz
Ekin
Epot
Figure 8: Comparison of the energy of the double pendulum for different frame rates of the camera. It is apparent that frame rate of 400 Hz incompatible with the principle of energy conservation; while the energy is decreasing in the long term due to friction, it should never increase. No choice of mass ratio mb/mg was able to resolve this issue. On the other hand, a frame rate of 500 Hz together with the mass ratio ma/ma- = 6.5 is compatible with energy conservation, within the bounds of error.
D Details on the Double Pendulum Dataset
D.1 Parameters
In Section 4.3 we demonstrated the applicability of our approach to the 'Double Pendulum Chaotic' dataset. A description of the dataset can be found in Asseman et al. (2018); to prevent confusion, we should mention that the blue and the green marker in our paper correspond to the green and the blue marker in Asseman et al. (2018), respectively (i.e. the colors have been exchanged). The lengths of the two pendula are given as To 91 mm and Ig = 70 mm, where the subfix refers to the pendulum with blue marker and g to the one with green marker. However, in order calculate the energy (up to a constant factor), knowledge of the masses, at least of the ratio ma/m, y is required. From information given by the authors of the paper and the manufacturer of the double pendulum, together with some experimentation of our own we estimated this ratio as ma/n m x 6.5. Note that in our description of the double pendulum that it consists of two point masses, which is only approximately true.
we made the assumption
Another quantity of interest the frame rate of the camera that was used to create the dataset; enters into the model when calculating the velocitios of the masses. In their paper (Asseman et al.] 2018), the authors state frame rate of 400 Hz. However, our experiments with the dataset and keeping the energy constraint in mind strongly indicate frame rate of 500 Hz; for 400 Hz there are segments of the motion where the total energy E clearly increases which violates the principle of energy conservation (see Figure 8
The 'Double Pendulum Chaotic' dataset was published under the "Community Data License Agreement
Sharing - Version 1.0".
D.2 Implementation Details
The 'Double Pendulum Chaotic' dataset provides data in the form of annotated positions of the masses attached to the ends of the two pendula (together with the position of the top of the apparature holding the pendulum which does not change and which we therefore omitted). We now have the positions as points on an equally spaced grid; in terms of the camera frame rate r the spacing between two adjacent points is given by 1/r. To obtain the velocities we numerically take the gradient of the positions on the grid and we receive the data which we use for our GP, with outputs
f [Zbar Zby: Egas Zgys Ubzi Ubys Ugas Ugy]" T
(D.1)
To obtain positions and velocities with comparable absolute values, which enhances the performance of the GP and which makes the quantities easier to compare in plots, we scaled positions by a factor of 20 and velocities by a factor of v10; the time t was scaled by factor of 5.
As outlined in Scction 4.3 we obtain training data, be used during hyperparameter optimization, and test. data, to evaluate the quality of predictions, by picking random interval of 200 datapoints from the
26

---

Published in Transactions on Machine Learning Research (11/2022)
second half of the trajectories provided by the dataset; out of those use 15 points training data and the rest test data. Note that the value E received by evaluating 11] and averaging over the training data will in general be less accurate estimate than the value of the energy Eg received when averaging over all datapoints in the interval, since the average is taken over fewer points in the former case. Hence, when determining the accuracy of the constraiint fulfillment, the results in Section 4.3 have been compared Eg-
For the double pendulum, we receive the transformed outputs
v = |2by Zgy. Unay by gz gy
(D.2)
with corresponding
[mag. mgg, 2 2 2
1ib
17th
mg
mg
(D.3)
2
The auxiliary outputs are
faux
[Thes Rgas Uba May Ugz: Ugul
(D.4)
note, that the outputs The and Ege are not actually auxiliary outputs, but since they are not involved in the constraint (11) (i.c. the corresponding entry in F would be zero), they can be learned separately from the constrained putputs, logether with the auxiliary outputs. Same as for the harmonic oscillator, we created virtual measurements for th bby Ugz 22 at zero crossings of the auxiliary outputs ubzi Ubys Uga Ugy-
E
Comparison to Jidling et al. (2017)
In this section we will investigate the parallels between the method of Jidling et al. (2017) and our own. While they only consider homogeneous constraints in the paper, in the context of constant linear sum constraints it is simple extend the method to affine constraints, as we will see below.
In the approach of Jidling al. (2017), vectors spanning the nullspace of the constraint F are used to construct the task covariance matrix. Given sum constraint F(f) Ff = S, and vectors h, spanning the mullspace (ie. Fh;- = 0), can define the matrix
G G h2... hn..
(E.1)
where n denotes the dimension of the nullspace. A suitable task covariance matrix can then be constructed via kg GGT, and in order accommodate the non-zero right hand side of the sum constraint, the task mean m, is chosen such that Fm, = S. Then we obtain the multivariate Gaussian N (m,, kt), samples of
which obey the constraint F. Note that k, is the projector the mullspace of the constraint; in Matthews et al. (017)) the relationship between constrained multivariate Gaussian distributions and the nullspace of the corresponding linear operator discussed. Subsequently, the full mean and covariance matrix of the GP can be constructed according to (4.
To make this more conerete, we consider again the example of the harmonic oscillator from Section 3.1.2 Here, the constraint given by F [k/2,m/3 O, 0] and E. Then the corresponding matrix G can be constructed (the choice of null vectors not unique)
0
0
Vm
R2
0
0
G
=
Vm
+
(E.2)
1
0
0
1
The task mean and task covariance matrix become
m
ma
0
0
m
k, GGT
mik
0
m;
(E.3)
1
0
0
1.
27

---

Published in Transactions on Machine Learning Research (11/2022)
Now if we approach the problem from the other side, as it turns out, starting with the identity matrix as task covariance matrix and then conditioning it according to (8) leads to the same k, obtained in [E.3 Here we note one advantage of our approach: it is straightforward to include additional correlations into the task covariance matrix, e.g. correlations between constrained outputs and those not involved in the constraint.
For example, when introducing an additional correlation between tasks one and three we obtain (after setting m 1)
0.5
0
0.5
-0.5
0.25
0
0
0
0
-0.5
0.5
-0.25
0
kg
- k4
(F.4)
0.5
0
1
0
0.25
0.25
0.875
0
0
0
1
0
0
1
where kg and ke are the unconstrained and the constrained task covariance matrix, respectively. So the constraint alters correlations between outputs involved the constraint and other outputs. In the approach of Jidling et al. (2017), correlations between the nullspace dimensions could be added by introducing a non-diagonal matrix between GT and G in E.3 However, that method would allow to introduce arbitrary correlations between the tasks demonstrated in (E.4). Throughout the work, we have used the structure given in B.12 as the starting point for our task covariance matrices.
28


Reviews:
[{"review_id": "oaURFPpfwU", "reviewer": "Reviewer_68RS", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strength\n\n* This paper addresses an important problem of Gaussian process modelling, especially in the multi-output setting. As demonstrated by examples from physics, vision, and control domains, real-world modelling tasks often require imposing equality constraints over multiple outputs of a function. Efficient inference schemes must be developed in order to use GPs as priors for such functions. Therefore, work that addresses this problem will likely be of interest to a wide community. \n\n* The idea is simple (reparameterizing nonlinear constraints as linear ones and fit gps to the transformed function values). It works well together with a clever trick that solves non-monotonic nonlinearity. \n\n* The method was tested using a variety of experiments in real-world scenarios. \n\n## Weaknesses\n\n* I found several claims in the paper clueless:\n   * Page 5, \"In case of a position dependent constraint, it is important to note that the values of the functions c(x) and a(x) must be known at all N_{tot} points.\" I don't think this is necessary. Correct me if I am wrong, but I believe it is possible to only condition the GP on the constraints over some subset of points. \n   * Section 3.2.1, \"It now suffices to enforce the constraints F on the task mean and covariance matrix\". I don't see how the constraints on the original GP translate to constraints on the task mean and covariance matrix. A proof should be given here if it exists. This is a major concern I have on the correctness of the proposed method.\n\n* The complexity of the proposed method seems high (impractical). How does it apply to large datasets? Could it be combined with sparse variational inference approaches?\n\n* There is artifacts in the uncertainty due to GP values might fall outside of the output domain of the variable (e.g., v^2 must be non-negative). This seems a major drawback of the proposed method. ", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This work propose a method to impose linear and nonlinear constraints on multiple outputs of a Gaussian process. It is based on the framework of multitask GPs, which decomposes the kernel into a kronecker product of task kernel and data kernel. Nonlinear constraints are transformed into linear ones through reparameterization. In order to address non-monotonic nonlinearities, an unconstrained GP is fit to the data to determine the region of the outputs. The paper also proposes to regularize the constrained GP with results from the unconstrained GP. Results on several toy and real-world datasets show that incorporating constraints often improves the prediction accuracy and reduce the credible intervals. ", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "I would like to see my comment on the two claims addressed. And I also believe the artifact problem is a big concern and should either be solved or by providing evidence that it does not affect practical utility of the approach. ", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "eV7mD9iWai", "reviewer": "Reviewer_JvJ3", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- The method is principled and appears to have a sensible derivation. While perhaps obvious, it is nice to see it written down all in one place so that people can refer to the details.\n\nWeaknesses\n- The various scenarios under which GP posteriors can be inferred from observations of linear transformations of GPs is well-known, and follows directly from the basic properties of finite-dimensional Gaussian vectors. The ``multi-task\" element does not really add much to the difficulty of the problem. The Laplace approximation is also well-known, as are its limitations. I don't think any reader will find the results here surprising or particularly enlightening, as the result could be conveyed with the summary I wrote above.\n\nAll-in-all, I am leaning very slightly towards acceptance, because while the results might not be surprising, it is nice to see all the details collected together in one place. ", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors consider the problem of inference using vector-valued Gaussian processes when certain equality constraints of the GP must be satisfied. These are, in order of difficulty, linear constraints (where the linear combination of a GP must be equal to some value), monotonic nonlinear constraints (where the linear combination of a monotonic nonlinear function of a GP is equal to some value), and arbitrary nonlinear constraints. Linear constraints are easily handled by the fact that evaluations of linear operators of Gaussian random variables are jointly Gaussian. Monotonic nonlinearities are handled via a Laplace approximation, and the underlying GP can be found by inverting the monotonic transformation. Arbitrary nonlinearities are handled in a similar way, but a pseudo-inverse is sought rather than an inverse. Three tasks are considered: a synthetic harmonic oscillator, a synthetic pose estimation problem, and a real double pendulum.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "As the authors mention in the Appendix,\n``There is no guarantee that Newton’s method will determine the correct maximum ˆf in case of multimodal\ndistributions, or that the resulting Gaussian distribution will constitute a good approximation of the true\nposterior. For these reasons, it has to be decided on a case by case basis whether the Laplace approximation\nshould be employed or not\"\n\nPerhaps I missed it, but I did not see anywhere in the examples a discussion on whether the Laplace approximation should be employed or not. Can you provide one? What are the advantages of Laplace approximation over other methods? Is it possible to try other methods for at least one of the problems you tried?\n\n", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "oL4dxFvHYT", "reviewer": "Reviewer_aAJV", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n\n- Incorporating prior knowledge in probabilistic modeling is relevant and of interest to the TMLR community. \n- The paper is easy-to-follow and well-written.  The paper uses a harmonic oscillator example throughout the text, which provides a nice demonstration of the ideas.\n- The proposed constrained multi-task GP outperforms alternatives in terms of constraint violations, while achieving comparable or better RMSE\n- Some limitations are discussed (discontinuities, potentially for overfitting due to virtual observations)\n\nWeaknesses\n\n- Depending on the transformation, the transformed outputs may not be well modeled with a standard stationary GP kernel. An RBF kernel is used in all experiments (for constrained and unconstrained GPs), but it seems unlikely that this is the ideal choice for both. To what extent is the choice of kernel impactful on this evaluation?\n- The Laplace approximation requires that the likelihood of the transformed observations be known. How common is this expected to be? Arbitrary non-linear transformations of a Gaussian distribution would not follow well known probability distributions. Is the proposed method applicable in such cases?\n - The quality of the Laplace approximation (i.e. the ability to approximate the likelihood with a Gaussian) will depend on the transformation. The Laplace approximation is used in all examples in this paper. The authors note that the approximation should be evaluated on a case by case basis. Is the proposed technique amenable to alternative approximations if the Laplace approximation is poor? \n- For arbitrary non-linear constraints, it seems that appropriate transformations may not be obvious. Are there cases where such transformation is not possible or consist of a very large pieces?", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This work tackles the problem of leveraging prior knowledge regarding sum constraints in multi-task gaussian process (GP) regression. These problems arise in physics domains, where often multiple possibly correlated outputs are known to obey sum constraints. This paper demonstrates how sum constraints that are non-linear functions of the outputs can be transformed into linear sum constraints by transforming the outputs. The paper then provides a procedure for fitting multi-task GPs under sum constraints and using them for prediction. The method is empirically evaluated on two synthetic and one real data example, where the proposed approach achieves lower constraint violation on out of sample prediction tasks than multi-task GPs that do not account for constraints and comparable or better RMSE.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "Important\n- In Sec 3.2.1, how are the constraints $F$ and $S$ changed when they are only applied to the inter-task covariance matrix directly?\n- It is unclear at several points in the paper if the setting considered is where all outputs are observed for each input. In section 2.1, a Kronecker product $\\otimes$ is used (as in 3.2.1), but it is only a Kronecker product if all outputs are observed for each input. In the harmonic oscillator example in figure 1, not all outputs are observed for each input. Hence, the multi-task kernel should be a Hadamard product, not a Kronecker product. Please unify and clarify this in the paper.\n\nStrengthen:\n\n- Is the auxiliary variable always simply $f$? Are there examples where other auxiliary variables should be used? \n- How are the inputs corresponding to virtual observations obtained? Are they found by finding where the posterior mean is zero?\n- In figure 1, it would be nice to show the unconstrained GP after applying the transformation to its outputs (e.g. by Monte-Carlo sampling in cases where analytic representation is not tractable).\n- What are the plus or minus terms in the tables? Are these numbers two standard errors of the mean? Please state this.\n- Sec 4.3, what is the training data (how many points and from which trajectories) for the double pendulum problem?\n- A metric that evaluates the quality of the uncertainty estimation such as log-likelihood would be nice to add.\n\nMinor Grammar:\n- There are a couple instances where hyphens are used instead of em dashes (e.g. Sec 2.1).\n- Sec 2.1: “and observations” -> and observation vector\n- Sec 4: missing space after period\n- Sec B2: Shouldn’t $\\boldsymbol f’$ be used instead of $\\boldsymbol f$ throughout this section? \n- Using $N_\\text{tot}$ for the number of training points and the number of training + test points is confusing\n- Figure 3: the right subplot is confusing. Is the left inset the unconstrained GP and the right inset the constrained GP? The caption is a little confusing, so a clearer description or legend would be helpful", "broader_impact_concerns": "n/a", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Almost all questions/concerns raised during the discussion phase were addressed and clarified in the latest version. The reviewers noted that only a restricted class of non-linear constraints has been considered, the occasional invalid outputs produced by the approximation and the ideas seem to be straightforward. However, based on the paper's correctness and relevance to the TMLR community, an acceptance decision is recommended.
Decision: Accept (poster)