Didier Rullière
email: didier.rulliere@emse.fr

Marc Grossouvre
email: marcgrossouvre@urbs.fr

A Joint Kriging Model with Application to Constrained Classification

Keywords: Multi-output Kriging, Cokriging, Constrained classification, Spatial Prediction, multi-task Gaussian Process regression

Interpolating or predicting data is of utmost importance in machine learning, and Gaussian Process Regression is one of the numerous techniques that is often used in practice. This paper considers the case of multi-input and multi-output data. It proposes a simple Joint Kriging model where common combination weights are applied to all output variables at the same time. This dramatically reduces the number of hyperparameters to be optimized, while keeping nice interpolating properties. An original constraint on predicted values is also introduced, useful for considering external information or adverse scenarios. Finally, it is shown that applied to membership degrees, the model is especially helpful for fuzzy classification problems. In particular, the model allows for prescribed average percentages of each class in predictions. Numerical illustrations are provided for both simulated and real data, and show the importance of the constraint on predicted values. The method also competes with state-of-the-art techniques on an open real world data set.

Introduction

Interpolating data is widely used in many fields of computer experiments, it is especially useful to predict values of one or several variables of interest, in the context of time-consuming or costly computer experiments. One considers here a Kriging interpolation problem on several output variables, with specific constraints on predicted values, and with possible applications to classification. Let us detail the need to deal with such problem.

Kriging on several outputs. Kriging or Gaussian Process Regression is a method of interpolation, especially suited when there are only few observations that have to be interpolated. It is widely used in many fields of Machine Learning, originally for geostatistical studies and spatial interpolation, but also for computer experiments in many domains (finance, industry, environment, etc.). The most basic Kriging theory aims a predicting a single real-valued quantity of interest, the output (for instance gold concentration in the ground), depending on some explanatory variables that are referred to as input values or locations (for instance latitude, longitude and depth). From a statistical point of view, the Kriging method is based on a best linear unbiased combination of observed outputs, with an assumption that observations are random variables whose correlation depends on locations. From a Gaussian random field point of view, in a Gaussian setting, the interpolation is the mean of a conditional Gaussian random field, with confidence bands derived from the variance of the conditional random field (see e.g. [START_REF] Williams | Gaussian processes for machine learning[END_REF], for an in-depth review).

The method has several advantages. First, it is interpretable: the prediction is a weighted average of observations, with quite logical behavior of weights. Second, the method fully interpolates the data, that is, predicts exactly an observed output if one uses the same input values. And third, it not only gives a prediction, but also confidence intervals for this prediction.

Among limitations of the method, and proposed extensions in the literature, one can cite the difficulty to handle a large number of observations (see e.g. [START_REF] Cressie | Fixed Rank Kriging for Very Large Spatial Data Sets[END_REF][START_REF] Banerjee | Efficient gaussian process regression for large datasets[END_REF]Rullière et al., 2018, and references therein), the difficulty to specify the covariance model and to estimate its hyperparameters (see e.g. [START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyper-parameters of gaussian processes with model misspecification[END_REF], the difficulty to treat multi-valued output (see [START_REF] Furrer | Aggregation-cokriging for highly multivariate spatial data[END_REF]).

In the present paper, we mainly consider this latter multi-valued output problem, which is clearly of practical interest. Co-Kriging techniques are built to treat several outputs, but there is usually one main output and others are used to improve the prediction of the main considered output. Furthermore, all cross-covariances between outputs at different locations have to be modelled, which creates O(p 2 ) covariance models, where p is the number of outputs (see e.g. [START_REF] Goovaerts | Ordinary cokriging revisited[END_REF][START_REF] Ver Hoef | Multivariable spatial prediction[END_REF][START_REF] Furrer | Aggregation-cokriging for highly multivariate spatial data[END_REF].

While highly parametrized models are useful in many situations, the prediction quality rely on the proper specification of the model and on the estimation of its parameters. A fine model with wrong parameters can some-times be less efficient than a simpler model with more control on few parameters (see [START_REF] Rasmussen | Occam's razor[END_REF]. In this paper, one considers a model where Kriging is applied to multi-valued outputs in R p , but with a specific simplification leading to a single covariance function to tune instead of O(p 2 ). This way, the prediction is feasible with a limited number of hyperparameters, interpretable and easy to explain.

One originality of the paper is that this kind of multi-valued interpolation is also applied to membership degrees in a classification setting. This allows to benefit from Kriging advantages (interpretation/explainability, interpolation, confidence intervals). Moreover, the proposed simplification of the model is especially useful since it keeps the property of membership degrees summing to one. At last, another novelty is to consider a specific constraint on predicted values. As detailed below, it will allow for proportion constraints in a classification setting.

Constraints on predicted values. There is a well known joke on actuaries: How much is two plus two? An actuary will ask "What do you want it to equal?". At first glance, it seems dishonest to require constraints on predicted values, especially if these constraints are very precise. But such constraints can be useful in some situations, when having external information, for adverse modelling or for homogenization needs, as illustrated below. The constraint we consider in this paper focuses on the average of predicted values.

It can be very helpful to prescribe a specific value for the average of predicted values. Let us instantiate some examples: due to an industrial accident, one wishes to measure the pollution in the soil for different chemical products. Measures are done at some spatial places, but the number of measures is limited. One would like to infer the quantity of all chemical product everywhere in the soil. Knowing stockpiles of products before the accident, the total quantity of lost chemicals may be known for every chemical products. While Gaussian Process Regression is especially suited to predict one product dosage in the soil, it has difficulty to handle jointly a lot of products, as it needs to model many cross-covariances. Furthermore it cannot handle at all constraints like prescribing the average of predicted values. Another example is the case where one needs to build a prediction under an adverse scenario: even if the total quantity of lost chemical is unknown, it can be useful to get an idea of the distribution of pollutants in an adverse case of massive loss. Finally, in other investigations, there may be external knowledge to consider: for example a regional study might want to be in line with some given national statistics, if there is no reason that the regional statistics differs in average. One can observe data due to an exceptional situation (e.g. COVID), and one may want to use it knowing that the situation has returned to normal. Or one might want predictions over different years or over different regions to coincide at least on average: for instance one may want that some disease incidence prediction does not differ, on average, over different medical centers. Another situation is the following one: imagine that one knows, under arbitrage free setting, that some predicted stock returns must be zero on average, or imagine that the regulator wants to force a prediction under specific shock scenarios. Fairness constraints can also be introduced to limit unfairness in algorithmic decision making (see e.g. [START_REF] Zafar | Fairness constraints: A flexible approach for fair classification[END_REF].

In a nutshell, prescribing the average value of predictions is useful in many contexts, be it external information (known quantity of chemical, national statistic...), adverse modelling (regulation, simulation under specific scenarios...), or need to homogenize results (over different regions, observed years, fairness constraints...).

Constrained Kriging and classification. Once Kriging is adapted to multi-output prediction, applying jointly Kriging on membership degrees has several advantages for classification.

One advantage is that the interpolation property can be preserved, which is not necessarily the case for other classification or clustering techniques like kNN: even very near from one observation, if the chosen number of neighbors is greater than one, kNN can propose another class than the observed one. When applying Kriging, we will see that some constraints on predicted values, or alternatively the use of a nugget effect, may break this interpolation property. But Kriging still offers the choice between prescribing an interpolation constraint, or prescribing more specific constraint on predicted values. Another advantage of using a Joint Kriging model on membership degrees is to get an estimation of the uncertainty relying on the prediction. For instance, at a specific location one may predict 10 percent of chance that the class is one, but one can also give a confidence interval for this quantity.

Applied to classification, a constraint on average predicted membership degrees is also useful: imagine that one interpolates membership degrees of different classes, depending on some explanatory variables. One may want to predict the membership degrees for many possible values of explanatory variables. But due to some external knowledge, to adverse scenarios modelling or to homogenization needs, one may want to prescribe the proportion in each class over all predicted values. Predicted class (or more precisely predicted membership degree) may differ depending on explanatory variables, but one may want for sure that in average, the proportion in the different considered classes is given. Some works can be found in the literature about clustering or classification under constraints: a survey on constrained classification can be found in [START_REF] Gordon | A survey of constrained classification[END_REF], see references therein. Some researches are treating size constraints for clustering [START_REF] Bradley | Constrained k-means clustering[END_REF]; [START_REF] Höppner | Clustering with size constraints[END_REF]; [START_REF] Ganganath | Data clustering with cluster size constraints using a modified k-means algorithm[END_REF], some are treating the problem of fuzzy clustering with weights (membership degrees), as in the present paper, see for example [START_REF] Benatti | Theoretical analysis of classic and capacity constrained fuzzy clustering[END_REF]. Fairness constraints are also considered in [START_REF] Zafar | Fairness constraints: A flexible approach for fair classification[END_REF].

Regarding Kriging and classification, some works on classification using Gaussian settings can be found in a dedicated Chapter 3 in the book [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF]. In particular, for binary classification, membership probabilities can be approximated by a sigmoid-like transformations of some latent Gaussian Process. The approach can be generalized to multiclass problems, and bayesian inference can be conducted using analytic approximations of integrals, or solutions based on Monte Carlo sampling (see [START_REF] Williams | Bayesian classification with gaussian processes[END_REF]Rasmussen et al., 2006, and references therein). Other recent approaches involving Multi-task Gaussian process, using several latent Gaussian processes and bayesian inference with approximations or sampling, can be found in [START_REF] Dahl | Grouped gaussian processes for solar power prediction[END_REF]; [START_REF] Panos | Large scale multi-label learning using gaussian processes[END_REF].

Among works closer to what is proposed in the present paper, Indicator Kriging (IK) aims at determining the cdf of an underlying random field at an unknown location, as a weighted average of indicators. It uses as well linear combinations of transformed observations, but relies on a direct link between indicators and underlying random field, using thresholds. Hence, it seems not directly suited to classification for non ordinal data, without any hierarchy between classes. It also requires the observation of the latent process that generates the indicators (see e.g. [START_REF] Journel | Nonparametric estimation of spatial distributions[END_REF][START_REF] Meer | Classification of remotely-sensed imagery using an indicator kriging approach: application to the problem of calcite-dolomite mineral mapping[END_REF][START_REF] Goovaerts | Auto-ik: A 2d indicator kriging program for the automated non-parametric modeling of local uncertainty in earth sciences[END_REF]Chiang et al., 2013, and references therein). Extensions like indicator co-Kriging require a large number of cross-covariances (see [START_REF] Agarwal | Copula-based multiple indicator kriging for non-gaussian random fields[END_REF]. In the present paper, the proposed method can be applied to non ordinal data, and does not require a specific model and thresholds between indicators of membership and an underlying real random field; furthermore in a simplified setting, the whole method can also rely on a single covariance function.

Proposal. To the best of our knowledge, the use of Kriging on several outputs with application to classification under constraints on predicted values has not been developed yet.

By itself, as previously detailed, Joint Kriging on several outputs has advantages for interpretability, for its capacity to interpolate data, for the uncertainty measurement associated to each prediction, for the limited number of hyperparameters and the simplicity of their estimation, while allowing specific model characteristics as periodicity for example. The present approach yields directly closed form formulas, without the need of conditional density approximations or sampling. With constraints on predicted values, Kriging is useful for using external information, for adverse modelling, for the need to homogenize results or to introduce fairness constraints. Such original constraints are not typically addressed by Kriging or Gaussian Process regression, but we show that with suitable adaptations, Kriging is an excellent choice for this task.

Applied to classification, the very same advantages hold: it allows for capacity constraints, as well as model uncertainty, interpretability, simple covariance modelling, and sometimes interpolation. Contrarily to Bayesian inference methods relying on integral approximations or Monte-Carlo sampling, this provides directly closed-form formulas for the predictor, using only basic linear algebra. It does not require creating a model where classes are derived from underlying Gaussian Processes and thresholds, so that it is also suited for non ordinal classes. It seems to us that using multi-outputs Kriging on classification offers many modelling perspectives, as well as practical results and performance: we will see that the proposed model compete with best available methods on an open data set.

The paper is structured as follows. In Section 2, one defines a simplified Kriging model that is suited for multi-valued outputs. The model is detailed in different cases: without specific constraint, similarly to Simple Kriging; with weights summing to 1, similarly to Ordinary Kriging; with constraints on weights summing to 1 and on average predicted values. In each case we derive optimal weights together with prediction mean and variance. An extension using an affine prediction is also developed. In Section 3, the proposed interpolation technique is applied to membership degrees, and it is shown that it preserves useful basic properties for the prediction. In Section 4, numerical applications of the proposed interpolation technique are given. One considers in particular a minimal application on a toy example, an illustration on a multi-valued time series on a real data set, and a more detailed real-world application on a classification problem. A conclusion closes the paper.

Appendix and supplementary material.

For more readability, all proofs are gathered in a Section A provided in Appendix. A list of notations is given in a dedicated section B. All illustrations in the paper are generated with notebooks that are available as supplementary material at https:// gitlab.emse.fr/marc.grossouvre/jointkrigingsupplementary/, in modifiable and executable format .Rmd and in already executed directly readable .html format (see [START_REF] Grossouvre | Supplementary material to: A Joint Kriging Model with Application to Constrained Classification[END_REF]. Hence the results are fully reproducible, and all specifications for drawing figures are easy to retrieve.

Joint Kriging model

Let us consider a multi-valued random field Y(x) := (Y 1 (x), . . . , Y p (x)) ⊤ ∈ R p , x ∈ X where X is a metric set of input points, typically X = R d .

For the sake of clarity and using analogy with geostatistics, we will refer to x as locations, but of course X may contain any explanatory variable. The components Y 1 (.), . . . , Y p (.) will be referred to as the p considered targets.

Components of the random field Y(x) can be dependent. Furthermore Y or its components are not necessarily Gaussian. However one assumes that first and second order moments exists. One considers here that Y(x) ∈ R p and X = R d , but other metric spaces would be possible, as soon as expectation and covariances between Y(x) and Y(x ′ ) can be derived.

Given n observations of Y(x 1 ), . . . , Y(x n ), we aim at predicting the values of the random field at some unobserved locations x ⋆ 1 , . . . , x ⋆ q , i.e. we aim at giving a predictor of Y(x ⋆ 1 ), . . . , Y(x ⋆ q ). At an unobserved location x ⋆ , we define the Joint Kriging predictor as a predictor M(x) = (M 1 (x), . . . , M p (x)) ⊤ depending linearly on observations, where real coefficients apply jointly to all components of the observations:

M(x ⋆ ) := n i=1 α i (x ⋆ )Y(x i ) , (1)
where each α i (x ⋆ ) ∈ R. These weights α(x ⋆ ) := (α 1 (x ⋆ ), . . . , α n (x ⋆ )) ⊤ are optimized in order to minimize some error that we will detail later on, under various possible constraints. Now, defining the p × n matrix

Y := [Y(x 1 ), . . . , Y(x p )],
Equation (1) also writes in a compact way

M(x ⋆ ) = Yα(x ⋆ ) . (2)
The main assumption behind this Joint Kriging predictor is that the weights are impacting all components the same way: the first component M 1 (x ⋆ ) is a linear combination of the observed first components, namely Y 1 (x 1 ), . . . , Y 1 (x n ), the second component M 2 (x ⋆ ) is the same linear combination of the observed second components Y 2 (x 1 ), . . . , Y 2 (x n ), etc. In other words, the weights affect jointly (or simultaneously) all the components of observed Y(x i ), i = 1, . . . , n, hence the chosen name of Joint Kriging model. This assumption can be quite logical, and easier to explain, to a decision maker: all targets are interpolated given observed targets, and the interpolation weights do not depend on the considered target. We will see in Section 3 that this key assumption is very useful, especially for classification under constraints. It would be technically possible to release this assumption, e.g. replacing weights α i (x) by some p × p matrix for i = 1, . . . , n: one would get closer to some general co-Kriging model with O(p 2 ) covariance models, but it is not the purpose of the present paper. Let us define the prediction error associated to a vector of weights α(x ⋆ ), at a prediction location x ⋆ . This loss is defined as the scalar value

∆(x ⋆ ) := E ∥M(x ⋆ ) -Y(x ⋆ )∥ 2 W . (3)
where ∥v∥ 2 W := v ⊤ Wv is a squared norm with W a given symmetrical positive-definite matrix of real weights. For instance, if one changes the unit of the first target, say multiply it by 100, then it sounds logical that the resulting norm is unchanged. Thus, some weights matrix may seem reasonable: an inverse covariance matrix as in the Mahalanobis distance, or a diagonal matrix of inverse variances, etc. For simplicity, the reader may imagine that all p targets are already scaled and that W is the p × p identity matrix.

The main difficulty is to derive the optimal weights α(x ⋆ ) under the various constraints one would like to consider. At all prediction locations x ⋆ 1 , . . . , x ⋆ q one thus aims at determining the optimal weights, gathered in a n × q matrix A := α(x ⋆ 1 ), . . . , α(x ⋆ q ) . This is performed in the three following subsections, under different constraints.

Optimal weights without constraints

In this subsection, one defines optimal weights minimizing the prediction error, without supplementary constraints.

The following Proposition expresses the weights such that M(x ⋆ ) is a best linear unbiased predictor (BLUP) of Y(x ⋆ ), in the sense of minimizing the loss (3). The result looks exactly the same as in the simple Kriging model, but the components in the symmetric positive semidefinite matrix K and in the vector h(x ⋆ ) here aggregate the values of all p observed targets. One retrieves usual Simple Kriging equations in the case where p = 1 and W is the identity matrix.

Proposition 1 (Simple Joint Kriging weights). The optimal weights α(x ⋆ ) minizing the loss of Equation (3) are given by the n × 1 vector:

α(x ⋆ ) = K -1 h(x ⋆ ) , (4)
or equivalently, using a matrix expression,

A = K -1 H , ( 5
)
where the n × n matrix K := E Y ⊤ WY is assumed to be invertible, the n × 1 vector h(x ⋆ ) := E Y ⊤ WY(x ⋆ ) , and the n×q matrix H := h(x ⋆ 1 ), . . . , h(x ⋆ q ) . If furthermore for any target j = 1, . . . , p, for any location x ∈ X , E [Y j (x)] = 0, then M(x ⋆ ) is unbiased.

Proof. The proof is postponed to Appendix. Note that the matrix K is necessarily a covariance matrix, since it is symmetric positive semidefinite. In particular, this appears clearly when E [Y j (x)] = 0 for any target j = 1, . . . , p and any location x ∈ X .

Optimal weights summing to one

In this section, one considers an additional constraint. This constraint raises naturally when Y i (x) are not centered, and leads to weights summing to one, as in Ordinary Kriging, see [START_REF] Cressie | Spatial prediction and ordinary kriging[END_REF], namely for all x ⋆ ,

α ⊤ (x ⋆ )1 n = 1 . ( 6
)
Where

1 n is a n × 1 vector of ones. It is clear that if for any x ∈ X , E [Y(x)] = µ then constraint (6) will imply that E [M(x ⋆ )] = µ, where µ = (µ 1 , . . . , µ p ) ⊤ , so that E [M(x ⋆ )] = E [Y(x ⋆ )].
The reverse implication is also true if all µ i are distinct from zero. Hence constraint ( 6) is a very natural constraint. It does not imply however that M(x ⋆ ) is a convex combination of all Y(x i ), because some weights can still be negative.

Under this constraint of weights summing to one, the following proposition gives the optimal weights. One retrieves exactly the same formulae as in ordinary Kriging, but the involved elements in matrices K and H are different: they are computed taking into account all p targets over all observations.

Proposition 2 (Ordinary Joint Kriging weights). Under the constraint of Equation (6), the optimal weights α(x ⋆ ) minimizing the loss of Equation (2) are given by the n × 1 vector α(x ⋆ ), together with the scalar λ(x ⋆ ):

α(x ⋆ ) = K -1 (h(x ⋆ ) + λ(x ⋆ )1 n ) λ(x ⋆ ) = 1 δ 1 -1 n ⊤ K -1 h(x ⋆ ) (7)
Equivalently, using matrix expressions, one gets

A = K -1 H + 1 n λ ⊤ λ ⊤ = 1 δ 1 q ⊤ -1 n ⊤ K -1 H (8)
where

K := E Y ⊤ WY , h(x ⋆ ) := E Y ⊤ WY(x ⋆ )
, and with scalar δ := 1 n ⊤ K -1 1 n . For matrix expressions, λ := (λ(x ⋆ 1 ), . . . , λ(x ⋆ q )) ⊤ , and the n × q covariance matrix H := h(x ⋆ 1 ), . . . , h(x ⋆ q ) . If furthermore, for all targets i = 1, . . . , p, for all locations

x ∈ X , E [Y i (x)] = µ i , then M(x ⋆ ) is unbiased.
Proof. The proof is postponed to Appendix.

An originality of the presentation of the result is that matrices can be expressed indifferently with compact expressions (using K and h(x ⋆ )) or with more classical covariances (using K and h(x ⋆ )), as stated in the next remark.

Remark 1 (Covariance matrices). Let us define the true unknown values of Y at all prediction points by

Y ⋆ := Y(x ⋆ 1 ), . . . , Y(x ⋆ q )
. Assume E [Y(x)] = µ for all x ∈ X . Furthermore, assume that either weights sum to one that is to say α(x ⋆ ) ⊤ 1 n = 1, or µ = 0 p .

Then the matrices K, H and the vector h(x ⋆ ) can be replaced by

     K = E Y ⊤ WY -E Y ⊤ W E [Y] H = E Y ⊤ WY ⋆ -E Y ⊤ W E [Y ⋆ ] h(x ⋆ ) = E Y ⊤ WY(x ⋆ ) -E Y ⊤ W E [Y(x ⋆ )] (9)
everywhere in Proposition 2, without changing the optimal weights α(x ⋆ ).

The next remark shows that the model can be built from a single covariance function k(x, x ′ ) between an implicitly weighted sum of all targets.

Remark 2 (Using correlation functions). Assume that there exists a positive definite matrix W, such that the covariances

k(x, x ′ ) := E Y(x) ⊤ WY(x ′ ) -E Y(x) ⊤ W E Y(x ′ )
depend only on some distance between x and x ′ .

Then the components of the covariances matrices K and H can be derived from the covariance function k : X × X → R by setting:

K ij = k(x i , x j ) H ik = k(x i , x ⋆ k ) (10)
In practice, with one hyperparameter for the variance σ 2 > 0, and d positive hyperparameters θ 1 , . . . , θ d , usually referred to as characteristic length-scales (see Rasmussen et al., 2006, bottom of p.14), one can set for example

k(x, x ′ ) = σ 2 r 0 ∥x -x ′ ∥ θ , (11)
where r 0 a unit correlation function and ∥x -

x ′ ∥ 2 θ = d i=1 x i -x ′ i θ i 2
corresponds to a rescaled Euclidean norm. Notice that this expression does not depend on W any more, so that when using the above assumption we do not have to estimate W.

In Remark 2, k(x, x) = σ 2 for all x ∈ X . It would be possible to consider a non homogeneous variance, this is not treated here for the sake of simplicity.

As a consequence of previous Remark 2, for a given matrix W (e.g. once the components variance is set to one, say), a noticeable advantage of the Joint Kriging method is the limited number of hyperparameters to optimize. Despite the multivariate output of the p × 1 response vector Y(x), x ∈ X , there is only a few hyperparameters to choose for defining the covariances: for instance σ 2 , θ, and the covariance family. This is quite different from co-Kriging techniques where all cross-covariances between components Y i (x) and Y j (x ′ ) should be defined, i, j ∈ {1, . . . , p}, x, x ′ ∈ X , which ends up in an order of O(p 2 ) covariance functions, and many associated hyperparameters.

Optimal weights with constraint on predictions

The constraint we consider here is more original than the previous one: we would like that, given observations Y, the average of predicted values has some prescribed value. Formally

E [M(X ⋆ ) | Y] = m . ( 12
)
where X ⋆ is a random variable taking values in prediction locations {x ⋆ 1 , . . . , x ⋆ q }, and where m is a given p × 1 vector. This constraint is quite original at it relies on predicted values for a given set of observations. The idea is to force the optimal weights to take into account this a posteriori constraint.

Notice the importance of the conditioning by Y, otherwise if all Y i (x) are centered, say, then the constraint would not be possible to satisfy in general, since all M(x ⋆ ) would be centered. We will see that this kind of constraint is particularly useful for classification, when one wishes to force classes proportions whatever the observed values.

Gathering all predictors in a p × q matrix M = M(x ⋆ 1 ), . . . , M(x ⋆ q ) , one can write M = YA. By conditioning on the value of X ⋆ , denoting

π x ⋆ = P [X ⋆ = x ⋆ ], x ⋆ ∈ {x ⋆ 1 , . . . , x ⋆ q }, Equation(12) simply writes YAπ = m . ( 13
)
where π is the q × 1 column vector π = (π x ⋆ 1 , . . . , π x ⋆ q ), and where the n × q matrix A gather all vectors of weights: A := α(x ⋆ 1 ), . . . , α(x ⋆ q ) . One specificity is that the resulting weights in matrix A will have to be solved all at once, for all q prediction locations. This is quite different from usual Kriging settings where prediction locations can be treated separately, if desired.

In the following proposition, we give the optimal weights matrix A when both constraints are considered at the same time: the constraint ( 12) on predicted values, and the constraint (6) on weights summing to one.

Proposition 3 (Joint Kriging weights under predicted values constraint).

The Joint Kriging weights minizing the loss of Equation (3) under the constraint of weights summing to one of Equation (6), and prescribed average predicted values of Equation (12) write:

A = K -1 H + 1 n λ ⊤ + Y ⊤ λ ′ π ⊤ (14)
with Lagrange multipliers, provided that

1 δ uu ⊤ -YK -1 Y ⊤ is invertible, λ ′ = γ -1 1 δ uu ⊤ -YK -1 Y ⊤ -1 YK -1 Hπ + 1 δ u 1 -1 n ⊤ K -1 Hπ -m λ = δ -1 1 q -H ⊤ K -1 1 n -πλ ′ ⊤ u where u := YK -1 1 n , γ := π ⊤ π ∈ R and δ := 1 n ⊤ K -1 1 n ∈ R, and
H := h(x ⋆ 1 ), ..., h(x ⋆ q ) . Proof. The proof is postponed to Appendix.
Note that it is also possible to compute a model with the only constraint YAπ = m, without requiring weights summing to one. In view of further classification applications, we do not develop it here and keep both constraints.

Again, an originality is that the previous result can be expressed using compact expressions for K and H or more classical covariances, as stated in the following remark. As a result, under suitable assumptions, all these covariance matrices can be filled using a covariance function k(., .), as detailed in Remark 2.

Remark 3 (Covariance matrices with two constraints). Under the assumptions of Remark 1 and using the same notations, the matrices K and H can be replaced by K and H everywhere in Proposition 3, without changing the optimal weights A.

Proof. The proof is postponed to Appendix.

Optimal weights with affine extension

Up to this point, one has only considered linear predictors, where a predictor is a linear combination of observed responses Y(x 1 ), . . . , Y(x n ), under various constraints. In this subsection, one considers the case where the prediction involves one supplementary term.

Remind the requirement on predicted values: they should be m in average. The hidden idea behind the knowledge of vector m is that there is some external information giving an hint on the prediction. This information may come for instance from some known overall statistics on the territory, some expert knowledge or from an expectancy estimator. Let us denote by Z the p × 1 random vector containing this external source of information.

With this in mind, it is quite natural to define an affine prediction:

M + (x ⋆ ) := α 0 (x ⋆ )Z + n i=1 α i (x ⋆ )Y(x i ) , (15)
Given Z = m, a constant term is included in the sum, hence the name affine prediction.

The sum of weights constraint on the new vector α + = (α 0 (x ⋆ ), . . . , α n (x ⋆ )) can be written

1 n+1 ⊤ α + (x ⋆ ) = 1 . ( 16
)
This way, if the p components of m and Y(x i ), i = 1, . . . , n are percentages summing to one, then the p components of the predictor M(x ⋆ ) will also sum to one.

For the second constraint on average predicted values, previously detailed in Equation ( 12), there is an implicit conditioning by Z = m. this constraint may write, with X ⋆ a r.v. defined on {x ⋆ 1 , . . . , x ⋆ q }:

E M + (X ⋆ ) Z = m, Y = m . (17)
Finally, provided covariances between Z and Y(x) are given, x ∈ X , then the setting is absolutely the same as in previous Propositions 1, 2 and 3, excepts one observation Z = m is added in the vectors of observations Y(x 1 ), . . . , Y(x n ). The covariance matrices are also updated. This is detailed in the following Proposition.

Proposition 4 (Affine version of predictors). Assume that the following covariance vectors are given

     P ⊤ := E Z ⊤ WY -E Z ⊤ W E [Y] Q ⊤ := E Z ⊤ WY ⋆ -E Z ⊤ W E [Y ⋆ ] σ 2 Z := E Z ⊤ WZ -E Z ⊤ W E [Z] (18)
Then affine predictors corresponding to the simple unconstrained case, to the ordinary case with one constraint, to the case with two constraints can be obtained by replacing Y, K, H by

Y + = m Y , K + = σ 2 Z P ⊤ P K , H + = Q ⊤ H , (19)
in Propositions 1, 2 and 3 respectively.

Proof. The proof is straightforward, hence not appearing in the Appendix.

Notice that the previous Proposition 4 can be easily extended to several sources of information Z 1 , Z 2 , . . .. For the sake of simplicity, this is not developed here.

Recall that, from Remark 2, that matrices K, H can be derived from simple correlation functions. Now it remains to derive one expression for P and Q.

Remark 4 (Extra covariances for affine prediction). Let P = (P 1 , . . . , P n ), Q = (Q 1 , . . . , Q q ) and Z = (Z 1 , . . . , Z p ). To fix ideas, in the case where W is the identity matrix, then for i = 1, . . . , n, j = 1, . . . , q

P i = p k=1 Cov [Z k , Y k (x i )] Q j = p k=1 Cov Z k , Y k (x ⋆ j )
This is easy to adapt when W is not the identity matrix, by considering the components of W 1/2 Z and W 1/2 Y(x), x ∈ X . Assuming that the above covariance does not depend on x i nor on x ⋆ j (i.e. the general source of information gives hints on the whole process, not on a particular location), then one can propose

P i = ρσσ Z , i = 1, . . . n Q j = ρσσ Z , j = 1, . . . q (20)
This corresponds for example to a model where

Y k (x) σ = ρ Z k σ Z + G k (x), k = 1, . . . , p, where all G k (x) are independent from all Z k .
Other assumptions can be chosen, leading to different vectors P and Q.

In the previous remark 4, the parameter ρ ∈ [-1, 1] measures how redundant is the information provided by Z, and can even be set to 0 if one considers that the external information source is completely independent from observations. The parameter σ Z measures how certain is the external information: when σ Z is high, the added information cannot be trusted and one retrieves the linear predictor, when σ Z is low, the added information is trustable, so that far from observed locations, M(x) gets nearer to m. In practice one can set 0 < σ Z ≪ σ to see the maximal difference with the linear predictor. One can even optimize this parameter σ Z to smoothly switch from a linear to an affine model.

Joint Kriging Mean and Variance

In this subsection we derive the mean predictor and the prediction error, assuming the optimal weights have been calculated with chosen constraints, as detailed in previous subsections.

Consider M(x ⋆ ) and α(x ⋆ ) a Joint Kriging predictor and the associated weights with or without constraints. In the following, we call Joint Kriging mean the value of the the predictor M(x ⋆ ) and Joint Kriging variance the value of the quadratic error ∆(x ⋆ ). Let us recall that:

M(x ⋆ ) := Yα(x ⋆ ) , ∆(x ⋆ ) := E ∥M(x ⋆ ) -Y(x ⋆ )∥ 2 W ,
where Y = [Y(x 1 ), . . . , Y(x p )] is the p × n matrix of observations. If p = 1 and if W is the identity matrix, Joint Kriging mean and Joint Kriging variance are exactly the Kriging mean and the Kriging variance usually known in Kriging.

The following Proposition gives a closed formula to compute the Joint Kriging variance.

Proposition 5 (Joint Kriging variance with arbitrary weights). Let α(x ⋆ ) be any vector of weights, possibly satisfying supplementary constraints.

The associated Joint Kriging variance writes:

∆(x ⋆ ) = α(x ⋆ ) ⊤ Kα(x ⋆ ) -2α(x ⋆ ) ⊤ h(x ⋆ ) + v(x ⋆ ) , (21)
or using a matrix expression, denoting ∆ := (∆(x ⋆ 1 ), . . . , ∆(x ⋆ 1 )) ⊤ , we get

∆ = diag A ⊤ KA -2 diag A ⊤ H + diag[K ⋆ ] ,
where the n×n matrix

K := E Y ⊤ WY , the q×q matrix K ⋆ := E Y ⋆⊤ WY ⋆ , the n×1 vector h(x ⋆ ) := E Y ⊤ WY(x ⋆ ) , and the scalar v(x ⋆ ) := E Y(x ⋆ ) ⊤ WY(x ⋆ )
are assumed to be known. diag [.] is the vector whose entries are the diagonal of the considered matrix.

Proof. The proof is postponed to Appendix.

Note that the above Proposition 5 can be directly adapted to the affine case of Proposition 4, by replacing Y, K, H by Y + , K + , H + , v being unchanged: one can interpret the predictor to be a linear predictor with one more observation, with correct covariances.

As previously stated in Remarks 1 and 3, and using the same notations, one can replace K, H, h with K, H, h, provided that the following new quantities are defined:

   K ⋆ = E Y ⋆⊤ WY ⋆ -E Y ⋆⊤ W E [Y ⋆ ] v(x ⋆ ) = E Y(x ⋆ ) ⊤ WY(x ⋆ ) -E Y(x ⋆ ) ⊤ W E [Y(x ⋆ )] .
The result is stated in Remark 5 below. Hence, in practice all these covariances can be filled using a given covariance function k(x, x ′ ), under suitable assumptions, as detailed in Remark 2.

Remark 5 (Covariance matrices in Joint Kriging mean and variance). Under the assumptions of Remark 1 and using the same notations, the matrices K, H, K ⋆ , the vector h(x ⋆ ) and the scalar v(x ⋆ ) can be replaced by K, H, K, h(x ⋆ ) and v(x ⋆ ) everywhere in Proposition 5, without changing the Joint Kriging mean and variance.

Proof. The proof is postponed to Appendix. Now, remark that Proposition 5 gives only an overall error

∆(x ⋆ ) := E ∥M(x ⋆ ) -Y(x ⋆ )∥ 2 W
which is a weighted sum of errors over all components of M(x ⋆ ).

This is a strength of the method, since the quantity to optimize is realvalued, which allows using standard covariance functions as detailed in Remark 2. This is also an important limitation, because in practice, one surely needs prediction errors for each component of M(x):

δ i (x ⋆ ) := E ∥M i (x ⋆ ) -Y i (x ⋆ )∥ 2 , i = 1, . . . , p .
The following Proposition 6 shows that one can get this error δ i (x ⋆ ) for each component i = 1, . . . , p. It relies on a supplementary assumption on the matrix W, but this assumption is only useful for determining the confidence bands for each component of the predictor M(x), not for computing M(x) itself.

Proposition 6 (Variance sharing). Assume that transformed observations Y(x) := W 1/2 Y(x) are such that components of Y are uncorrelated and bear the same share of the covariance function k, that is to say:

Cov Y i (x), Y j (x ′ ) = 1 p k(x, x ′ )1 {i=j} , i, j ∈ {1, . . . , p}, x, x ′ ∈ X , (22)
Assume also that E [Y(x)] = µ for all x ∈ X . Furthermore, assume that either the weights sum to one, or µ = 0 p . Then the local errors write

δ i (x ⋆ ) = σ 2 i σ 2 ∆(x ⋆ ) , i = 1, . . . , p . ( 23
)
where

σ 2 i := Var [Y i (x)
] is the variance of the component Y i (x), assumed to be constant over x.

Proof. The proof is postponed to Appendix.

The result of Proposition 6 is quite logical, it simply states that for a well chosen matrix W, the error δ i (x ⋆ ) is proportional to the unit global error σ -2 ∆(x ⋆ ): one simply has to apply the variance σ 2 i of the component instead of the variance σ 2 of the aggregated weighted components.

Application to constrained classification

In this section, we now apply multi-output prediction to membership degrees for fuzzy classification. We show that the Joint Kriging predictor, together with constraints on weights and predicted values, is especially suited to this task since it remains very simple, and constraints are very logical in a classification setting.

Consider a classification problem with p possible labels. Labels are depending on some explanatory variables x ∈ X , so that one may observe labels ℓ(x 1 ), . . . , ℓ(x n ) taking values in {1, . . . , p}.

At an unobserved location x ⋆ , for a predictor L(x ⋆ ) ∈ {1, . . . , p} of ℓ(x ⋆ ), the reader may convince himself that, for given probabilities p j , j ∈ {1, . . . , p}, constraints on predicted classification such as

P [L(X ⋆ ) = j | observed labels] = p j , (24)
are not so easy to handle, even if X ⋆ is a uniformly distributed random variable over prediction points x ⋆ 1 , . . . , x ⋆ q . This is because such constraints usually act in a non-linear way on the predictor L(X ⋆ ), and the predictor L(.) itself can be some complicated function of observed labels. Existing predictors, such as Indicator Kriging, may be unable to deal with such constraints. Furthermore, they may not be appropriate in cases where the considered labels do not correspond to ordinal classes.

The originality here is to use Joint Kriging to propose a fuzzy classification under both constraints: membership degrees summing to one, and prescribed average of predicted degrees. Indeed, in a practical context, it is natural to require that over all predictions, predicted values are distributed like the observed ones, for instance. Or the sampling bias of observations is known so that the expected label percentages are known. Or sometimes an external source of information gives the expected label percentages. It can be the case for a regional study, knowing some statistics at a national level. At last, it can be used for modelling adverse scenarios.

Technically, in a classification problem, each label ℓ ∈ {1, . . . , p} can be converted into a p × 1 vector of indicator functions, namely

Y := 1 {j=ℓ} j=1,...,p .
This transformation is well known in the machine learning community as label binarization (see also dummy variables or one-hot encoding), and is implemented in many languages. Obviously using this simple representation,

1 p ⊤ Y = 1.
In practice, it is common to observe true label values depending on some explanatory variables x ∈ X . But it may also happen that one observes uncertain labels: multiple and distinct observed labels for a same x ∈ X , uncertainty relying on the value of x, etc. To handle this problem, one generalizes slightly the previous label binarization: one assumes here that observations consist in a distribution of possible labels, so that one observes n vectors Y(x 1 ), . . . , Y(x n ), such that the components of each vector are summing to one:

1 p ⊤ Y(x i ) = 1, i = 1, . . . , n.
In other words, the p components of Y(x i ) represent the membership degrees of the p possible classes, at an observed location x i , i = 1, . . . , n. Using previous notations, recall that Y = [Y(x 1 ), . . . , Y(x n )], so that observed membership degrees satisfy

1 p ⊤ Y = 1 n ⊤ . ( 25
)
Now, using a Joint Kriging model, one can infer the membership degree of an unobserved location x ⋆ , using the predictor of Equation (1):

M(x ⋆ ) := n i=1 α i (x ⋆ )Y(x i ) (26)
The next remark details the impact of both constraints, weights summing to one and prescribed average predicted values, in the particular setting of membership degrees that are summing to one. Despite being very straightforward, it aims in particular at recalling the previous results in the classification context.

Remark 6 (Constraints impact). Consider the two previous constraints on weights and predicted values, namely constraint of Equation (6) and Equation (12) . Consider also the membership degree assumption given in Equation (25), 1 p ⊤ Y = 1 q ⊤ . Then the Joint Kriging model imply that

• Predicted membership degrees are summing to one:

1 p ⊤ M(x ⋆ ) = 1 ,
for any prediction point x ⋆ ∈ X . In particular

1 p ⊤ M = 1 q ⊤ .
• Average membership degree over prediction points can be chosen:

E [M(X ⋆ ) | Y] = m ,
where m is a prescribed average of predicted membership degrees of each class, with 1 p ⊤ m = 1, and X ⋆ a random variable over all prediction points.

Proof. The proof is postponed to Appendix.

As noticed before, although predicted membership degrees are summing to one, no guarantee of positivity is forced for the predicted values. In practice, as it is the case with numerous machine learning methods, a posttreatment of results may be required in specific cases involving negative membership degrees (see e.g. the use of softmax function with neural networks).

Recall that, computing the Joint Kriging predictor requires the knowledge of cross moments of Y(x) and Y(x ′ ), namely E Y(x) ⊤ WY(x ′ ) for all observed and predicted locations, i.e. for x, x ′ ∈ {x 1 , . . . , x n } ∪ {x ⋆ 1 , . . . , x ⋆ q }, where W is a given definite-positive matrix. In practice, the Remark 2 can be used to build all covariance matrices from a single covariance function: it is enough to choose a covariance function k(x, x ′ ) with few associated hyperparameters to get all needed covariance matrices.

Other techniques could be used for deriving covariances. Ordinal labels may derive from, say, an underlying 1D Gaussian Process U (x), x ∈ X , by setting Y(x) = 1 {g(U (x))=j} j=1,...,p , where g : R → {1, . . . , p} is a given function. The derivation of all covariances from the ones of U (.) is feasible, as with Indicator Kriging, but out of the scope of the present paper. And finally, it also ends up in the choice of a covariance function k(x, x ′ ).

Numerical illustrations

In this section one considers different numerical illustrations, for both prediction and classification. The first illustration focuses on the impact of constraints with one output, the second one on the behavior of the predictor with multiple outputs. The third illustration gives an application to classification and a benchmark with numerous competitors. All the illustrations are created in R markdown notebooks, one per subsection, available as a supplementary material at https://gitlab.emse.fr/marc.grossouvre/ jointkrigingsupplementary/ [START_REF] Grossouvre | Supplementary material to: A Joint Kriging Model with Application to Constrained Classification[END_REF]. Notebooks are given in both an executable format and an executed html format. The presented figures are directly created from the notebooks and results are fully reproducible.

A simplified toy example

One considers here the very simple case where there is a single target: the output Y(x) is belonging to R p , with p = 1. The interest for testing the Joint Kriging with one single target is to discuss the impact of the constraint on predicted values, and the impact of the affine prediction.

For p = 1, Simple Joint Kriging and Ordinary Joint Kriging are identical to common Simple Kriging and Ordinary Kriging, but the constraint on predicted values leads to a new original predictor. We keep here the vector bold font for vectors Y(x) ∈ R p and m ∈ R p , even though p = 1, in order to keep the very same notations as in the rest of the paper.

Let us consider that the process Y(x) aims at approximating an hidden function, with say a = 1 and b = 4, f (x) := a + sin(x/b) .

Observed locations x 1 , . . . , x n are randomly chosen with a uniform distribution over the interval [-10, 5], and q prediction locations x ⋆ 1 , . . . , x ⋆ q are chosen regularly spaced over the interval [-3, 10]. Both intervals are purposely shifted so that some prediction points are far from observations, and vice-versa.

Observed responses in R p , with p = 1, are Y(x i ) = f (x i ), i = 1, . . . , n with n = 10. Prediction is made over a set of q = 100 points. One defines X ⋆ a discrete uniform random variable over all prediction points.

The purpose here is not to interpolate as precisely as possible the hidden function f given few observations, but only to illustrate the differences between various possible interpolators, and the impact of requiring a prescribed average values for predicted values.

The prescribed value for m ∈ R p , with p = 1, is the scalar m = 1.5. The covariances between Y(.) are modelled as prescribed in Remark 2, from a single covariance function, using a squared exponential kernel. One could also pick a kernel that reflects f periodicity. However, the purpose is not to do the best possible prediction but, rather, to understand the impact of various constraints.

Cov Y(x), Y(x ′ ) = k(x, x ′ ) = σ 2 exp - (x -x ′ ) 2 2θ 2 .
We set σ 2 = 0.6 mainly for the visibility of the confidence band in presented figures, and θ = 1.2.

In the Figure 1, one exclusively considers the constraint of sum of weights, which is assumed to be one:

1 n ⊤ α = 1.
The predictor M(x) appears in red thick line, together with confidence intervals built from the variance ∆(x).

The Panel 1a presents the result of ordinary Kriging exposed in Proposition 2. As is well known, when the location x is large (and far from observed locations), the ordinary Kriging mean tends to return to the estimated mean of the observations. The average value of the Kriging mean E [M(X ⋆ ) | Y] ≃ 1.12 is quite different from the value m = 1.5 (horizontal dashed line), which is natural as this constraint has not been taken into account yet.

The Panel 1b uses the Proposition 4 to add a supplementary affine term to the linear combination, while preserving the sum of weights being one. The affine term is derived from a random variable Z, and we choose σ Z = σ/10, so that this external information is assumed trustable (small variance). Given Z = m, the consequence is that, far from observed locations, the prediction tends to put all weight on this external source of information, so that the prediction gets closer to m, as one can see at the extreme right of this Figure 1b. This also makes the average E [M(X ⋆ ) | Y] ≃ 1.37 closer to m = 1.5, but the values of those two quantities remain distinct. Another consequence of the affine term is the reduction of the confidence band width, as a new source of information has been added.

In Figure 2, one considers both the constraint of sum of weights, which is assumed to be one: 1 n ⊤ α = 1, together with the prescribed average of predicted values E [M(X ⋆ ) | Y] = m. The predictor M(x) appears in thick blue line, together with confidence intervals built from the variance ∆(x).

The Panel 2a presents the result of ordinary Kriging exposed in Proposition 3. The average value of the Kriging mean E [M(X ⋆ ) | Y] = 1.5 is exactly the prescribed one m = 1.5 (horizontal dashed line), which is natural as this constraint has been taken into account during the joint optimization of all α(x ⋆ j ), j = 1, . . . , q. However, the predictor is no more interpolating. This is logical: if q = 1, one has one only prediction point

x ⋆ 1 , the constraint E [M(X ⋆ ) | Y] = m becomes M(x ⋆ 1 ) = m, which is distinct to an observation Y(x n ), even if x ⋆
1 gets closer to x n . Another example: if on the one hand observation points and prediction points are the same, if on the other hand m is not the average value of observations, then at least one prediction must be different from the associated observation to satisfy the constraint. The Panel 2b uses Proposition 4 to add a supplementary affine term to the previous linear predictor of Panel 2a, while preserving the sum of weights being equal to one. The affine term is derived from a random variable Z, and we choose as previously σ Z = σ/10, so that this external information is assumed trustable. As above, the average of predicted values is exactly the prescribed one, by construction. Again, given Z = m, the consequence is that, far from observed locations, the prediction tends to put all weights on this external source of information, so that the prediction gets closer to m, as one can see at the extreme right of this Figure 2b. Another consequence of the affine term is the reduction of the confidence band width, as a new source of information has been added. With the prescribed average of predicted value, the predictor is not interpolating, but adding the affine term helps the prediction to get closer to observations. In this simple toy example, one can check numerically that each prediction satisfies the constraints that it should. One can also clearly visualize the impact of the specific constraint on average predicted values, and the behavior of the predictor when adding an affine term.

The illustrations that have been presented in this subsection are available in the supplementary material notebook Application1D.

A multi-output time series example

In the previous example, we have illustrated the impact of constraints on the prediction of a one-dimensional output. Hence, the joint aspect of the estimation was not discussed. In the present example, one considers a multioutput data, so as to illustrate the specificity of the single hyperparameter estimation with multiple outputs. We choose one dimensional inputs in R to facilitate the interpolation representation, but considering more general inputs in R d , d > 1, would be easy. It would only change the number of hyperparameters to estimate, d instead of 1.

Imagine the following situation: a city wants to infer the history of some pollutants concentration at a particular crossroad based on a small series of measurements. This simple problem requires a model that takes time as input and multiple concentrations as output. Obviously, the end purpose would be to have a model with space and time as input but this is out of this illustration's framework.

Using the data air quality (see [START_REF] Vito | Air Quality. UCI Machine Learning Repository[END_REF], one tries to infer the concentration of several pollutants, from only few values. Studied pollutants in the data where chosen arbitrarily: CO, C6H6, NOx and NO2. The time range of learning data has been selected so that visually there is not too much missing data in the period (sensor stuck to an inferior bound or missing), it corresponds to hourly measurements from 23/04/2004 18.00.00 to 28/04/2004 17.00.00. Missing values are tagged with -200 values in this data, they have been all filtered before the study, as if they were not informative at all. The challenge is to predict all hourly measurements in the selected period from only n = 10 values.

The purpose here is not to give specific conclusions about the measured pollution, but only to illustrate the capacity of the Joint Kriging model to handle complex multi-valued data, with very few hyperparameters to optimize. The idea is to create a joint model that would be as simple as possible. Many refinements of the model could be suggested, but this is not the purpose of this example.

Let us model the covariances between components of Y(.) using Remark 2. The proposed method does not require the definition of each crosscorrelation between a pollutant concentration at one location and a different pollutant concentration at a different location. It just takes one covariance function k(x, x ′ ) between an implicitly weighted sum of all targets. We use the multiplication of two covariance kernels (hence it is positive semidefinite): a periodic kernel with period of one day, and a kernel of the Matérn 3/2 family (see Chapter 4 and Equation (4.31) in [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF].

k(x, x ′ ) = σ 2 exp -sin 2 (π|x -x ′ |) 1 + |x -x ′ | θ exp - |x -x ′ | θ .
(27) The parameterization has been simplified, e.g. factors √ 3 in Matérn covariance expressions are not used here: they have the same effect as a rescaling of the characteristic length-scale θ. Notice that despite the p dimensional output where p = 4 is the number of studied pollutants, the kernel k(x, x ′ ) in Equation ( 27) depends only on two hyperparameters θ and σ 2 . Since σ 2 impacts the uncertainty measurement but not the prediction itself, it is set to σ 2 = 1.

Let us consider first one single constraint: the sum of weights should be one. It corresponds to the Joint Ordinary Kriging predictor.

Figure 3 shows the optimization of the single length-scale hyperparameter θ. As this study does not aim at comparing the prediction accuracy with other methods, we did not use a separate test sample, but only a validation sample, keeping in mind that it may lead to an overfit. The validation data used for this single hyperparameter estimation is set to all hourly measurements in the selected period. For the hyperparameter optimization, a specific error has been chosen, where one optimizes the worst standardized mean absolute error over all p = 4 series: the errors have been standardized in order to make them unitless and scale invariant. The best estimation is θ ≃ 1.4, it is kept for all other illustrations of the subsection.

The optimization here depends quite heavily on the chosen observation locations, so that in practice, an averaged error on several training and validation data would probably be more stable: in many situations on real data, the error function is monotonic, either increasing and leading to extremely small optimized hyperparameter θ (the prediction then tends to return quickly to an average value), or either decreasing, leading to a very large value of θ (the prediction then tends to smooth data a lot). Classical co-Kriging strategies that define a large number of cross-covariance hyperparameters would probably worsen the situation, highlighting the utility of a small number of hyperparameters.

Figure 4 presents the simultaneous predictions of the four pollutant concentrations with Joint Kriging, the only constraint being that weights sum to 1. Confidence band associated with a given pollutant is represented proportional to the standard deviation of this pollutant's concentration as detailed in Proposition 6. Pollutant concentrations have very different orders of magnitude but when applying Proposition 6, the obtained confidence bands looks quite comparable between series, as desired.

With very few hyperparameters and with a rough covariance model, the result has a lot of room for improvement. Nevertheless, despite the single model hyperparameter θ, and considering the limited number of observations n = 10, the predictions of the p = 4 concentrations seem quite reasonable. By construction each prediction is a combination of observed values of the considered pollutant, with weights summing to one. No other Considering this single hyperparameter model, with a basic covariance model, the results also seem reasonable when using two constraints. In the left panels, the average of predicted values is exactly set to 30% more than the observed average of pollutant concentration, which is a lot. However the visual differences between true and predicted sequences looks surprisingly moderate, even in this adverse scenario: despite satisfying all constraints, the model still offers a good fit with observations. The goal of this numerical experiment is to demonstrate the Joint Kriging model's ability to handle complex multi-valued data. It also illustrates the advantage of having a limited number of hyperparameters. One sees here that with a quite simple model, in a difficult problem (predicting four quite erratic time series from 10 observations), the model performs reasonably well. Furthermore, it allows for introducing some constraints, like setting an adverse scenario of 30% increase of the pollutant concentration.

The illustrations that have been presented in this subsection are available in the supplementary material notebook ApplicationAirQuality.

A constrained classification example

We present in this subsection the specific case of multi-dimensional outputs derived from a classification problem. As presented in section 3, Joint Kriging can be implemented for fuzzy classification. Different modalities of a classification variable are regarded as multiple output variables with values in [0, 1].

Imagine the case of an event that may occur at a given location of a territory, with a measurable intensity. We are interested in classifying the intensity of this event, if it occurs, into multiple classes, depending on some thresholds. In the following, this event is an earthquake and its intensity is its Richter magnitude.

The Quake data set given in [START_REF] Simonoff | Smoothing Methods in Statistics[END_REF], visualized in Figure 6, describes 2 178 earthquakes with their latitude, longitude, focal depth and magnitude. A given location x has coordinates latitude, longitude and focal depth. For a single observation at location x, the target

Y(x) = (Y 1 (x), Y 2 (x)) ⊤ is equal to (1, 0) ⊤ if
an earthquake is occurring here with a magnitude above the data set average magnitude, or (0, 1) ⊤ otherwise. If a location x is observed repeatedly, the membership degrees at x are averaged out over observations. It makes sense to impose that membership degrees are summing to one, so that 1 p ⊤ Y(x) = 1. Extensions with more thresholds would be easy to conduct, see Figure 10, we keep here p = 2 for comparison to existing benchmarks. The binarized data is available at www.openml.org/search?type=data&id=772, on openML website (see [START_REF] Bischl | OpenML benchmarking suites[END_REF]. Figure 6: Earthquakes observations. An earthquake is a point with coordinates latitude, longitude and focal depth (given by the color). Triangles represent earthquakes which magnitude is above average. Circles represent earthquakes which magnitude is below average.

The purpose here is to compare the performance of Joint Kriging with a set of 69 other models' performances. The study available at www. openml.org/search?type=task&id=4516 compares models, called flows in openML, performing 10 times a 10-fold cross-validation and computing the predictive accuracy as performance indicator (see tab Analysis, measure predictive_accuracy, the user must click on "fetch next 100 runs" several times in order to see all benchmarked models).

Remember from Remark 2 that although we are constructing a bivariate model, we need a single covariance kernel. The latter should be periodic with respect to the latitude and longitude, not with respect to focal depth. A simple way to define an admissible kernel is to multiply 3 kernels associated with the 3 dimensions (see [START_REF] Williams | Gaussian processes for machine learning[END_REF]:

k(x, x ′ ) = σ 2 exp -2 sin 2 ((x 1 -x ′ 1 )/2) θ 2 1 -2 sin 2 ((x 2 -x ′ 2 )/2) θ 2 2 -2 (x 3 -x ′ 3 ) 2 θ 2 3
The hyperparameters estimation has been treated separately on other train/test splits to avoid overfitting the data. The resulting values for θ are 2.3 for latitude, 0.9 for longitude and 196.8 for focal depth.

In order to visualize the algorithm's behaviour, we predict on a grid of latitude, longitude and focal depth values. In addition to imposing the sum of membership degrees to be 1, we set the output mean expectation to be the same as in the data set. Predictions on a grid of latitude, longitude and focal depth are presented on Figure 7. One can observe that maps representing membership degrees (first two rows) can be deduced from each other by y = 1 -x. The third row shows a segmentation of the plane into areas where membership degree for "P: magnitude is greater than average" is greater than 0.5, and areas where the converse is true. This segmentation depends on the focal depth: small focal depth on the top row (21km, first quartile) and greater one on the bottom row (68km, third quartile). For instance, looking at the bottom left corner of the map, which is around the Fiji archipelago, one can predict that earthquakes with small focal depth are more likely to be of large magnitude than deep earthquakes. However the converse is true in the south Atlantic area (bottom center part of the map). Moreover, the predictor achieves a circular coherence along longitude due to the periodicity of covariance. Periodicity along latitude is more difficult to observe because is covers only 180 • . Performances are evaluated using the Predictive Accuracy: it is the percentage of instances that are classified correctly. It is measured on binarized predicted membership degrees, on a 10 times 10-fold cross-validation, as in the OpenML benchmark, in order to get comparable results. Prior to that, the hyperparameters optimization has been treated separately on other train/test splits in order not to overfit the data.

Figure 8 presents from top to bottom: two results found in openML i.e. the best recorded model which is kernel logistic regression with RBFK kernel and Random Forest for reference, below are presented results of Joint Kriging models i.e. simple model without constraint, with weights summing to 1, with constraint on the prediction and weights summing to 1, affine with weights summing to 1 and affine with constrained output.

For the 10 runs, each diagram shows the Predictive Accuracy of each run (colored points), the minimum, first quartile, median, 3rd quartile and maximum, as well as the mean value materialized by a cross. Although the runs' performances stay in the range of those observed for Random Forest and Kernel Logistic Regression, the average values obtained with Joint Kriging are greater: the average is 0.556 ± 0.018 for the best model in OpenML benchmark, and 0.5661 ± 0.0038 for the best Joint Kriging model. The latter was even slightly greater, 0.5669, during hyperparameter optimisation, due to a slight overfit that has been reduced when using different train/test splits. Benchmark being based on this average value, it means that Joint Kriging has a better performance than the 69 models tested in the OpenML benchmark.

One can expect the multiplication of constraints to have an adverse effect on performance as a constrained optimization has less degree of freedom than an unconstrained one. On the other hand, injecting useful information through constraints may improve the performance. Figure 8 shows that overall, the performance is improved, especially when adding the constraint on the output. The whisker plots show that the Joint Kriging performance is less dispersed than that of presented OpenML competitors. This is partly due to the fact that lengthscales have been optimized separately, so that the dispersion of their estimators is not taken into account: the variability rely on different train/test splitting. This prior optimisation slightly reduces the average measured performance of Joint Kriging models, due to different train/test split, but it also reduces its dispersion. From top to bottom: membership degree of "P: magnitude is above average", membership degree of "N: magnitude is below average", binarized prediction (1 if membership degree of P is greater than 0.5). Left: focal depth of 21km. Right: focal depth of 68km.

In Figure 9, one uses the affine version of Joint Kriging with two constraints, weights summing to one and prescribed average prediction. In left panels, an adverse scenario forces the average predicted membership degrees of the first class (large magnitude events) to be equal to 65%. In the right panels, this percentage is set to the observed percentage of large magnitude events, 55%. This illustrates the usefulness of the constraint for adverse modeling.

In order to compare the results with existing benchmarks, we studied above the p = 2 binary classification problem. But the method can handle more classes as well. As an example, in Figure 10, we give a prediction for p = 4 classes. Observations have been converted into four classes, using three Richter magnitude thresholds 5.85, 5.95, 6.15. Specific thresholds have been chosen for this illustration, in order to get enough observations in each classes (at least 17% observations), but a seismology study might focus on other thresholds. Once again, the predictor achieves a circular coherence along longitude, and one can observe complex patterns that would be difficult to catch with classification trees. The presented classification task was constructed from indicators deriving from an underlying real value, the Richter magnitude, and from thresholds, thus creating ordinal classes. But the prediction can be derived as well for observations of non ordinal class labels, without any underlying process or thresholds. From top to bottom: membership degrees of "1: magnitude is smaller than 5.85", "2: magnitude is between 5.85 and 5.95", "3: magnitude is between 5.95 and 6.15", "4: magnitude is greater than 6.15" and class of greatest membership degree in the 5th row colored by increasing magnitude from dark to light blue. From left to right: focal depth of 21km, focal depth of 68km.

In this classification example, we used a direct implementation of the model with a single covariance family. Nevertheless the average performance is the best one among the whole OpenML benchmark, fine tuning of the model would surely lead to performance improvements. The illustration above aims at demonstrating that with very basic assumptions, the method is competitive with an open benchmark that has numerous competitors, as shown in Figure 8. It also aims at showing that it can model adverse scenarios, as in Figure 9, or multiple classes, as in Figure 10.

The illustrations that have been presented in this subsection are available in the supplementary material ApplicationClassificationQuake.

Conclusion

A Joint Kriging model on multiple outputs has been presented, where at each prediction location, the same weights apply to all outputs. This simplification allows for an easy covariance modelling, with very few hyperparameters even though the number of outputs p is large. Still, the model benefits from Kriging advantages: interpretability, ability to interpolate data, uncertainty measurement associated to each prediction, specific covariance modelling. As any simplification, the model can surely be improved and may have some limitations compared to heavily parameterized models: e.g. coKriging with many cross-covariance functions might be more flexible for dealing with time series with different regularities, or models with parameterized distortions of locations might be convenient for dealing with non-stationarities. However, the limited number of hyperparameters and the simplicity of their estimation is an asset of the model, while allowing specific model characteristics as periodicity. Furthermore the model is not limited to Gaussian Processes as it only relies on the existence of moments of order one and two.

An original constraint on predicted values was also introduced. It appears to be useful for using external information, for adverse modelling, for homogenizing results or for considering fairness constraints. To handle this constraint, all weights of predicted points need to be computed at the same time, unlike usual Kriging techniques. But the resulting predictor itself is quite simple to derive since it is given by a closed formula. Some extensions using an affine term were also proposed, allowing to account for an external information, and providing more control on the behavior of the predictor far from observations.

Ultimately, an application to classification was developed. Applying a multi-outputs Kriging model on classification is feasible through the prediction of membership degrees. Even without constraints, it is in itself interesting: it allows for interpretability, modelling uncertainty estimation and interpolating data. Using Joint Kriging with the proposed constraints eas-ily ensures that membership degrees sum to one, and allows for prescribed percentages of each predicted class. The simplified covariance model eases a lot the hyperparameters estimation. At the same time, with Joint Kriging, classification tasks benefit from the diversity of covariance kernels including periodicity. The resulting classification performs especially well in the investigated practical case: in the quake numerical example, the model competes with the best provided approaches on an open data set with numerous competitors.

Multiple extensions to the model can be imagined. For instance, the model with constrained predicted values does not guarantee continuous interpolation so that further work may fix this problem. A specific estimation procedure of the underlying joint covariance structure could also be of interest. Moreover, once applied to classification, membership degrees summing to one do not imply the combinations to be convex: some weights can still be negative or greater than 1. Ensuring the combinations to be convex could also be an improvement.

A Proofs

Proof of Proposition 1 Simple Joint Kriging weights. The proof is very similar to the geo-statistical proof of Simple Kriging model. It does not rely on any Gaussian assumption, but just on existing moments of order two. Recall that W is a symmetrical positive definite matrix, so that W = W ⊤ . Let us calculate the gradient of ∆(x ⋆ ) with respect to α(x ⋆ ):

∇ α(x ⋆ ) E ∥M(x ⋆ ) -Y(x ⋆ )∥ 2 W = ∇ α(x ⋆ ) E (M(x ⋆ ) -Y(x ⋆ )) ⊤ W (M(x ⋆ ) -Y(x ⋆ )) = ∇ α(x ⋆ ) E (Yα(x ⋆ ) -Y(x ⋆ )) ⊤ W (Yα(x ⋆ ) -Y(x ⋆ )) = ∇ α(x ⋆ ) E α(x ⋆ ) ⊤ Y ⊤ WYα(x ⋆ ) -2α(x ⋆ ) ⊤ Y ⊤ WY(x ⋆ ) + Y(x ⋆ ) ⊤ WY(x ⋆ ) = 2 E Y ⊤ WY α(x ⋆ ) -2 E Y ⊤ WY(x ⋆ ) . Thus, ∇ α(x ⋆ ) ∆(x * ) = 2Kα(x ⋆ ) -2h(x ⋆ ) . (28)
Where 4) when the gradient is zero. The matrix expression A = K -1 H of Equation ( 5) is obtained by binding column vectors of Equation ( 4), for all prediction locations. Remark that, under assumption that E [Y(x)] = 0 p for all x ∈ X , it is clear that E [M(x ⋆ )] = E [Y(x ⋆ )] = 0, so that the predictor is unbiased.

K := E Y ⊤ WY is a n × n matrix and h(x ⋆ ) := E Y ⊤ WY(x ⋆ ) is a n × 1 vector, thus leading to a n × 1 gradient. Hence α(x ⋆ ) = K -1 h(x ⋆ ) in Equation (
In that case, the (i, j) component of the matrix

Y ⊤ Y is (E Y ⊤ Y ) ij = n k=1 E [Y i (x k )Y j (x k )] = n k=1 Cov [Y i (x k ), Y j (x k )] .
Hence Y ⊤ Y is a symmetric positive semi-definite matrix. The same holds for K:

writing K = W 1/2 Y ⊤ W 1/2 Y , it is clear that for any vector v, v ⊤ Kv = v ⊤ v ≥ 0,
where the vector v := W 1/2 Yv. Thus K is a symmetric semi-definite positive matrix, i.e. a covariance matrix.

Proof of Proposition 2 Ordinary Joint Kriging weights. Under the constraint (6), and using a Lagrange multiplier λ ∈ R, the loss to minimize is ∆ 1 (x ⋆ ) := ∆(x ⋆ ) -2λ(x ⋆ ) α(x ⋆ ) ⊤ 1 n -1

Using Equation ( 28), the gradient of ∆ 1 (x ⋆ ) with respect to α(x ⋆ ) is

∇ α(x ⋆ ) ∆ 1 (x ⋆ ) = 2 E Y ⊤ WY α(x ⋆ ) -2 E Y ⊤ WY(x ⋆ ) -2λ(x ⋆ )1 n (29)
Setting this ∇ α(x ⋆ ) ∆ 1 (x ⋆ ) to be zero for all of its p components, we get

E Y ⊤ WY α(x ⋆ ) = E Y ⊤ WY(x ⋆ ) + λ(x ⋆ )1 n .
and finally, Kα(x ⋆ ) = h(x ⋆ ) + λ(x ⋆ )1 n . Once α(x ⋆ ) is written as a function of λ(x ⋆ ), one easily gets the value of λ(x ⋆ ) by setting 1 n ⊤ α(x ⋆ ) = 1. Hence the result. Matrix expressions are obtained by binding column vectors for all x ⋆ in {x ⋆ 1 , . . . , x ⋆ q }

Proof of Remark 1 Covariance matrices. Recall that K = E Y ⊤ WY and h(x ⋆ ) = E Y ⊤ WY(x ⋆ ) . Under the chosen mean assumption, both E [Y(x ⋆ )] = µ and E [Y] = µ1 n ⊤ . Thus, under the given constraint α(x ⋆ ) ⊤ 1 n = 1, or when µ = 0 p ,

E Y ⊤ W E [Y] α(x ⋆ ) = E Y ⊤ W E [Y(x ⋆ )] = 1 n µ ⊤ Wµ .
Hence the gradient of ∆(x ⋆ ) in Equation ( 28) also writes ∇ α(x ⋆ ) ∆(x * ) = 2Kα(x ⋆ ) -2h(x ⋆ ) = 2 Kα(x ⋆ ) -2 h(x ⋆ ) .

As a consequence, the gradient of ∆ 1 (x ⋆ ) in Equation ( 29) is unchanged when replacing both (K, h) by ( K, h). Thus, one can freely replace both (K, h) by ( K, h) in the rest of the proof of Proposition 2, without changing the result.

Proof of Proposition 3 Joint Kriging weights under predicted values constraint. Under both constraints, the quantity to minimize is

∆ 2 (x ⋆ ) := ∆(x ⋆ ) -2λ(x ⋆ ) α(x ⋆ ) ⊤ 1 n -1 -2λ ′ ⊤ (YAπ -m) ,
where λ ′ is a p × 1 vector of Lagrange multipliers. The gradient of the last term, with respect to α(x ⋆ ) is

∇ α(x ⋆ ) 2λ ′ ⊤ (E [M(X ⋆ ) | Y] -m) = ∇ α(x ⋆ ) 2λ ′ ⊤ (P [X ⋆ = x ⋆ ] E [M(x ⋆ ) | Y] + P [X ⋆ ̸ = x ⋆ ] E [M(X ⋆ ) | X ⋆ ̸ = x ⋆ , Y] -m) = ∇ α(x ⋆ ) 2λ ′ ⊤ (P [X ⋆ = x ⋆ ] E [M(x ⋆ ) | Y] -m) + 0 = ∇ α(x ⋆ ) 2λ ′ ⊤ (P [X ⋆ = x ⋆ ] E [Yα(x ⋆ ) -m) | Y] = ∇ α(x ⋆ ) 2λ ′ ⊤ (π x ⋆ Yα(x ⋆ ) -m) = 2π x ⋆ Y ⊤ λ ′
Hence, provided the matrix factor is invertible,

λ ′ = γ -1 1 δ uu ⊤ -YK -1 Y ⊤ -1 YK -1 Hπ + 1 δ u 1 -1 n ⊤ K -1 Hπ -m
Once λ ′ computed, one gets for λ

1 q ⊤ = 1 n ⊤ K -1 H + δλ ⊤ + 1 n ⊤ K -1 Y ⊤ λ ′ π ⊤ δλ ⊤ = -1 n ⊤ K -1 H -1 n ⊤ K -1 Y ⊤ λ ′ π ⊤ + 1 q ⊤
And finally, using u = YK -1 1 n ,

λ = δ -1 1 q -H ⊤ K -1 1 n -πλ ′ ⊤ u
Proof of Remark 3 Covariance matrices with two constraints. The proof is similar to the one of Remark 1 and uses the fact that, under chosen assumptions and for any prediction point x ⋆ ,

Kα(x ⋆ ) -h(x ⋆ ) = Kα(x ⋆ ) -h(x ⋆ ) .
Hence the gradient of ∆ 2 (x ⋆ ) in Equation ( 30) is unchanged when replacing K and h(x ⋆ ) by K and h(x ⋆ ), and all further expressions follows the same way in the proof of Proposition 3.

Proof of Proposition 5 Joint Kriging variance with arbitrary weights. The first equation is a simple vector rewriting of Equation (1). For the prediction error, one simply write, whatever the weights α(x ⋆ ),

∆(x ⋆ ) = E ∥M(x ⋆ ) -Y(x ⋆ )∥ 2 W = E (M(x ⋆ ) -Y(x ⋆ )) ⊤ W (M(x ⋆ ) -Y(x ⋆ )) = E (Yα(x ⋆ ) -Y(x ⋆ )) ⊤ W (Yα(x ⋆ ) -Y(x ⋆ )) = E α(x ⋆ ) ⊤ Y ⊤ WYα(x ⋆ ) -2α(x ⋆ ) ⊤ Y ⊤ WY(x ⋆ ) + Y(x ⋆ ) ⊤ WY(x ⋆ ) .
Hence the result.

Proof of Remark 5 Covariance matrices in Joint Kriging mean and variance.

The case where µ = 0 p is straightforward, as in that case K = K, h(x ⋆ ) = h(x ⋆ ) and v(x ⋆ ) = v(x ⋆ ), whatever the weights α(x ⋆ ). It remains the case where weights are summing to one. As in previous remarks, under chosen assumptions one gets Kα(x ⋆ ) -h(x ⋆ ) = Kα(x ⋆ ) -h(x ⋆ ) , and moreover one can show that

-α(x ⋆ ) ⊤ h(x ⋆ ) + v(x ⋆ ) = -α(x ⋆ ) ⊤ h(x ⋆ ) + v(x ⋆ ) .
Hence the result.

Proof of Proposition 6 Variance sharing. The difficulty here is to derive the cross-covariance k ij (x, x ′ ) = Cov [Y i (x), Y j (x ′ )] from the expression of k(x, x ′ ) that is detailed in Remark 2

k(x, x ′ ) := E Y(x) ⊤ WY(x ′ ) -E Y(x) ⊤ W E Y(x ′ )
Denoting Y(x) := W 1/2 Y(x), x ∈ X , this scalar covariance writes

k(x, x ′ ) = E Y(x) ⊤ Y(x ′ ) -E Y(x) ⊤ E Y(x ′ ) (31)
One would like to compute the p × p cross-covariance matrix between Y(x) and Y(x ′ ), using Y(x) = W -1/2 Y(x), x ∈ X :

K Y (x, x ′ ) := E Y(x)Y(x ′ ) ⊤ -E [Y(x)] E Y(x ′ ) ⊤ = W -1/2 E Y(x) Y(x ′ ) ⊤ -E Y(x) E Y(x ′ ) ⊤ W -1/2 ⊤ = W -1/2 K Y (x, x ′ )W -1/2 ⊤ (32)
where one defines

K Y (x, x ′ ) := E Y(x) Y(x ′ ) ⊤ -E Y(x) E Y(x ′ ) ⊤ .

Now assume that:

Cov Y i (x), Y j (x ′ ) = 0 whenever i ̸ = j, x, x ′ ∈ X .

This implies that W 1/2 is proportional to a whitening transformation, so that all components of Y 1 (x), . . . , Y p (x) are uncorrelated. Assume furthermore that:

Cov Y 1 (x), Y 1 (x ′ ) = . . . = Cov Y p (x), Y p (x ′ ) , x, x ′ ∈ X .

Then one easily sees from Equation (31) that the scalar k(x, x ′ ) satisfies

k(x, x ′ ) = p i=1
Cov Y i (x), Y i (x ′ ) = p Cov Y j (x), Y j (x ′ ) , j = 1, . . . , p

Hence under these assumptions, denoting I p the p × p identity matrix,

K Y (x, x ′ ) = 1 p k(x, x ′ )I p .
As a consequence, from Equation ( 32),

K Y (x, x ′ ) := E Y(x)Y(x ′ ) ⊤ -E [Y(x)] E Y(x ′ ) ⊤ = 1 p k(x, x ′ )W -1 (33) Cov Y i (x), Y j (x ′ ) = 1 p k(x, x ′ )(W -1 ) ij (34)
Now from this, one can derive the local cross errors

δ ij (x, x ′ ) := E (M i (x) -Y i (x)) M j (x ′ ) -Y j (x ′ )
Let us denote by Y i. the ith row vector of the matrix Y. We get

δ ij (x, x ′ ) = E (Y i. α(x) -Y i (x)) Y j. α(x ′ ) -Y j (x ′ ) = α(x) ⊤ E Y i. ⊤ Y j. α(x ′ ) -α(x) ⊤ E Y i. ⊤ Y j (x ′ ) -E Y i (x) ⊤ Y j. α(x ′ ) + E Y i (x)Y j (x ′ )
Now assume E [Y(x)] = µ for all x ∈ X . Then from Equation ( 34),

E Y i (x) ⊤ Y j (x ′ ) -µ i µ j = 1 p k(x, x ′ )(W -1 ) ij
which implies, using the matrix K defined in Equation ( 10):

E Y i. ⊤ Y j. -µ i µ j 1 n 1 n ⊤ = 1 p (W -1 ) ij K and E Y i. ⊤ Y j (x ′ ) -µ i µ j 1 n = 1 p (W -1 ) ij h(x ′ ) .
Furthermore, assume that either weights sum to one or µ = 0 p , then terms in µ i µ j vanish and one gets:

δ ij (x, x ′ ) = 1 p (W -1 ) ij α(x) ⊤ Kα(x ′ ) -α(x) ⊤ h(x ′ ) -h ⊤ (x)α(x ′ ) + k(x, x ′ ) .
In particular from Proposition 5, using Remark 2 and Remark 5,

δ i (x ⋆ ) = 1 p (W -1 ) ii ∆(x ⋆ ) . (35)
From Equation (33), when k(x, x) = σ 2 for all x, one can write

K Y (x, x) = 1 p σ 2 W -1 .
Using σ 2 i := Var [Y i (x)], assumed to be constant over x,

1 p (W -1 ) ii = (K Y (x, x)) ii σ 2 = σ 2 i σ 2 .
Hence from Equation ( 35),

δ i (x ⋆ ) = σ 2 i σ 2 ∆(x ⋆ ) .
Proof of Remark 6 Constraints impact. The result is a very straightforward rewriting and interpretation of constraints ( 6) and ( 12). From YAπ = m one derives 1 p ⊤ m = 1 p ⊤ YAπ = 1 n ⊤ Aπ = 1 q ⊤ π = 1, hence the natural constraint on prescribed average membership degrees in m, that must sum to one.

Figure 1 :

 1 Figure1: Prediction with one constraint: weights summing to one. Left: linear predictor, and right: affine predictor. In both cases, the average of predicted value is distinct from the prescribed value m = 1.5 (horizontal dashed line). Observations are black dots, the thin dotted blue line is the underlying function. In the right panel, one applies the assumption in Remark 4 with ρ = 0 and σ Z = σ/10.

Figure 2 :

 2 Figure2: Prediction with two constraints: weights summing to one, and average of predicted values set to m = 1.5 (horizontal dashed line). The average of predictions is equal to this value m = 1.5 in both cases. Observations are black dots, the thin dotted blue line is the underlying function. In the right panel, one applies the assumption in Remark 4 with ρ = 0 and σ Z = σ/10.

Figure 3 :

 3 Figure 3: Optimization of the single correlation hyperparameter θ for the four selected pollutants, data extracted from Air quality data set.

Figure 4 :

 4 Figure 4: Joint Kriging interpolation: using the joint ordinary model with weights summing to one, with very few data points (black dots) and a single optimized length-scale hyperparameter obtained in Figure 3. Upper left: CO, upper right: C6H6, lower left: NOx, lower right: NO2. Predictions are in thick solid lines, true values are in thin black solid lines.

FocalFigure 7 :

 7 Figure7: Joint Kriging with 2 constraints, earthquakes magnitude prediction into 2 classes. From top to bottom: membership degree of "P: magnitude is above average", membership degree of "N: magnitude is below average", binarized prediction (1 if membership degree of P is greater than 0.5). Left: focal depth of 21km. Right: focal depth of 68km.

FocalFigure 9 :Figure 10 :

 910 Figure9: Adverse scenario: predicted membership degrees of earthquakes magnitude using Joint Kriging with two constraints. Top panels: adverse scenario, first class output average constrained to be 65%. Bottom panels: regular scenario, output average constrained to 55.5%. Left: focal depth of 21km. Right: focal depth of 68km.

  Distribution of performances for 10 runs of two OpenML models which are kernel logistic regression (best model in benchmark) and Random Forest, and the 5 different types of Joint Kriging model . The whisker plots give minimum, first quartile, median, third quartile and maximum. The dark red cross indicates the average predictive accuracy, the higher the better. Average is 0.566 1 ± 0.003 8 for the best Joint Kriging model, 0.556 ± 0.018 for the best model in OpenML benchmark.

	Best model in benchmark (openML)	seed
		Random Forest	1235
	(from benchmark)	1236
		Simple	1237
	Model	With weights summing to 1	1238 1239 1240
	With constrained
		output	1241
		Affine with	1242
	weights summing to 1	1243
		Affine with	1244
	constrained output
		0.50 0.51 0.52 0.53 0.54 0.55 0.56 0.57 0.58 0.59
		Predictive accuracy
	Figure 8:

constraint is added on Figure 4, so that the average of predictions does not correspond at all to a specific prescribed value.

Figure 5 presents the simultaneous predictions of the four pollutant concentrations with Joint Kriging on which both constraints on the weights and on the predicted values are imposed using the affine model of Remark 4. Left panels show an adverse scenarios where the average of predictions is set to 130% of the true average, right panels are normal scenarios where the average of predictions is set to 100% of the true average. Using this setting, the interpolation property is lost, as seen in the previous example of Section 4.1, but the n = 10 observations still have a large influence, and the global shape of the prediction is preserved. By construction, the average of predicted values (thick solid line) is exactly the prescribed one (horizontal thick dashed dark green line).

Hence using the gradient of ∆(x ⋆ ) in Equation ( 28), one gets

Setting the gradient to be equal to a n × 1 vector of zeros, we get for all prediction locations

As optimal weights are gathered in the n × q matrix A := α(x ⋆ 1 ), ..., α(x ⋆ q ) , if one defines the n × q matrix H := h(x ⋆ 1 ), ..., h(x ⋆ q ) , then the previous system can be written, by binding columns for all x ⋆ ∈ {x ⋆ 1 , . . . , x ⋆ q }:

with q × 1 Lagrange multiplier λ, and p × 1 Lagrange multiplier λ ′ . If K is invertible, then the first equation writes

Injecting this value of A into the first constraint 1 n ⊤ A = 1 q ⊤ , denoting γ := π ⊤ π ∈ R and δ := 1 n ⊤ K -1 1 n ∈ R one gets:

Now injecting the value of A into the second constraint YAπ = m, and using the last equation, denoting the p × 1 vector u := YK -1 1 n , one gets

and finally, the vector λ ′ must satisfies

B Notations

locations X set of locations (inputs/design points). n, q number of observed locations, of prediction locations.

x any location. x 1 , . . . , x n are all observed locations.

x ⋆ any prediction location. x ⋆ 1 , . . . , x ⋆ q are all prediction locations. X ⋆ a random variable over prediction locations. π = (π x ⋆ 1 , . . . , π x ⋆ q ) the q × 1 distribution of X ⋆ over prediction locations. γ = π ⊤ π an intermediate real value used in calculations. targets p number of targets (i.e. number of outputs). Y(x) the p × 1 vector of targets at location x.

1 ), . . . , Y(x ⋆ q ) all p × q unknown targets at prediction locations.

. . , M(x ⋆ q ) the p × q matrix of all predictions. α(x ⋆ ) the n × 1 linear weights for the prediction in x ⋆ . A = (α(x ⋆ 1 ), . . . , α(x ⋆ q ) the n × q matrix of weights for all predictions. m a given constant p × 1 vector of prescribed mean predicted values. ∆(x ⋆ ), ∆ 1 (x ⋆ ), ∆ 2 (x ⋆ ) losses to be minimized for finding M(x ⋆ ). λ a q × 1 vector of Lagrange multipliers (relative to sum of weights) λ ′ a p × 1 vector of Lagrange multipliers (relative to predicted values) u = YK -1 1 n an intermediate p × 1 vector in calculations. Z an additional p × 1 factor for affine predictions. covariances W a given symmetric positive definite matrix for computing norms.

), H other covariances using centred expressions. δ = 1 n ⊤ K -1 1 n an intermediate real value in calculations. P additional n × 1 covariance vector between Z and Y(x i ) Q additional q × 1 covariance vector between Z and Y(x ⋆ j ) miscellaneous v a generic vector for defining norm or checking psd characteristic.

1 n , 1 p , 1 q a vector of ones of size n, p, q respectively. 0 n , 0 p , 0 q a vector of zeros of size n, p, q respectively.