A Mutually-Dependent Hadamard Kernel for Modelling Latent Variable Couplings
Sami Remes Markus Heinonen Samuel Kaski


The revival of the Baldwin Effect
    Mauro Santos
    December 30, 2023
=================================



 
We introduce a novel kernel that models input-dependent couplings across multiple latent processes.
The pairwise joint kernel measures covariance along inputs and across different latent signals in a mutually-dependent fashion. 
A latent correlation Gaussian process (LCGP) model combines these non-stationary latent components into multiple outputs by an input-dependent mixing matrix. Probit classification and support for multiple observation sets are derived by Variational Bayesian inference. Results on several datasets indicate that the LCGP model can recover the correlations between latent signals while simultaneously achieving state-of-the-art performance. We highlight the latent covariances with an EEG classification dataset where latent brain processes and their couplings simultaneously emerge from the model.
 


Gaussian process, non-stationary kernel, cross-covariance, latent variable modelling





§ INTRODUCTION



Gaussian processes (GP) are Bayesian non-parametric models that explicitly characterize the uncertainty in the learned model by describing distributions over functions <cit.>. These models assume a prior over functions, and subsequently the function posterior given the data can be derived. The prior covariance plays the key roles of both regularising the model by determining its smoothness properties, and characterising how the underlying function varies in the input space.


Recently, there has been interest in deriving non-stationary covariance kernels, where the general signal variances or the intrinsic kernel parameters – such as the lengthscales in the squared exponential or Matérn kernels – are input-dependent <cit.>. For instance, in geostatistical applications, a non-stationary kernel can both model a difference in the covariance along or across geological formations <cit.>. Input-dependent, heteroscedastic noise models have also been studied in single-task <cit.> and in multi-task settings <cit.>. 


In multi-task learning Gaussian processes are utilized by modeling the output covariances between possibly several latent functions <cit.>. In latent function models[Coined as linear models of coregionalisation (LCM) in geostatistics literature <cit.>.] the outputs are linear combinations of multiple underlying latent functions <cit.>. In Gaussian Process Regression Networks (GPRN) the mixing coefficients of multiple independent latent signals are input-dependent Gaussian processes as well, leading to a general multi-task framework that adaptively combines latent signals into outputs along the input space <cit.>.


The main contribution of this paper is to introduce a mutually-dependent Hadamard product kernel that combines a covariance structure between the latent signals that depends on the inputs, with an input kernel that depends on the latent signal indices. The signal and input kernels are interdependent, conditional on each other. This is in contrast to earlier Kronecker-based joint kernels where inputs and latent signals would be assumed independent. The kernel generalizes Wishart processes <cit.> into cross-covariances for input-dependent correlation structure, and a non-stationary Gaussian kernel <cit.> for measuring input-space correlations at specific latent signals. We deploy this kernel to extend the GPRN framework by a non-stationary cross-covariance function for the latent signals.

Furthermore, the proposed latent correlation Gaussian process (LCGP) incorporates multiple latent signals that are linearly combined into multiple outputs in an input-dependent fashion. The latent signals have a structured Wishart-Gibbs model that leads to non-stationary signal variances. We account for both regression, and Probit-based classification. Finally, the model is extended for multiple observation sets, where each observation is modeled by a separate latent model with shared latent correlations. In such a model, the latent correlations effectively regularize the latent models of each observation. Variational Bayesian inference with whitened gradients is derived for a scalable implementation. 

We highlight the model with several datasets where interesting latent signal covariance models emerge, while retaining or improving the state-of-the-art regression and classification performance. Multi-observation classification is demonstrated on EEG data from a large set of scalp measurements from several subjects, where the model is able to learn the covariance model between the underlying brain processes. In simulation studies, we show that our model is capable of accurately learning the latent variable correlations.




§ LATENT CORRELATION GAUSSIAN PROCESS


We consider M-dimensional observations (x) ∈^M over N data points (x_1, …, x_N). We denote vectors with boldface symbols, matrices with capital symbols and block matrices with boldface capital symbols. In this section we first construct the multi-output regression model for (x), and then develop a novel kernel for latent variables in such a model as our main contribution. Section <ref> further extends the framework into a classification setup.



 §.§ Multi-output regression


Following <cit.>, we model the M-dimensional outputs (x) ∈^M as an input-dependent mixture of Q latent signals (̆x) ∈^Q via a mixing matrix B(x) ∈^M × Q,

    (x)    = (x) +  = B(x)((̆x) + ) + ,

where = (x) is zero-mean M-dimensional Gaussian observation noise and = (x) is zero-mean Q-dimensional latent noise

    ∼(, ω_f^-1I⃗),   ω_f ∼( α_f, β_f), 
       ∼(, ω_u^-1I⃗),    ω_u ∼(α_u, β_u)   .


We model both the latent variables $̆ as well as the elements of the mixing matrixBas Gaussian processes. A GP priorϕ(x) ∼(μ(x),k(x,x'))defines a distribution over functionsϕ(x)with expectation𝔼[ϕ(x)] = μ(x)and covariance of the values between pointsxandx'is[ϕ(x),ϕ(x')] = k(x,x'). A set of function valuesϕ⃗= (ϕ(x_1), …, ϕ(x_N))^Tfollows a Gaussianϕ⃗∼(μ, K)withK_ij = k(x_i,x_j)andμ_i = μ(x_i).

The mixing matrixB(x)is anM ×Qmatrix of independent Gaussian processes over outputsmand latent signalsq,

    B_mq(x) ∼(0, k_b(x,x')).

The kernelk_b(x,x')between two input pointsxandx'determines how mixing of latent signals into outputs evolves along the input space. For instance, with temporal data the mixing matrix allows time-dependent linear combinations of the outputs. The full model is depicted in Figure <ref>. Next, we proceed to derive a kernel for the latent variables$̆.



 §.§ Wishart-Gibbs Hadamard Product Kernel




The latent signals u_p(x) are functions of the signal index p and input x. We propose to encode the latent signals (̆x) as mutually dependent Gaussian processes over pairs of inputs (x, x') and signals (p, q),

    u_p(x)    ∼(0, A_xx'(p,q) K_pq(x,x') ),

such that the joint covariance [ u_p(x), u_q(x')] = A_xx'(p,q) K_pq(x,x') is a product of signal and input similarities. Both similarities depend on each other to produce a non-stationary joint covariance. 


The pairwise, mutually dependent Hadamard kernel k(x,x',p,q) = A_xx'(p,q) K_pq(x,x') encodes a rich similarity between input x of latent signal p and x' of latent signal q as the product of the two conditional kernels. The kernel A_xx'(p,q) encodes signal similarity between inputs x and x', while the kernel K_pq(x,x') denotes input similarity at latent signals p and q. Since the two kernels depend on each other, a simple model such as Kronecker kernel product <cit.> is not suitable. Both kernels can be interpreted as cross-covariances. The Gibbs kernel K_pq restricts the flexibility of the Wishart kernel A_xx' (See Figure <ref>).

For instance, in EEG data the kernels could signify correlations A_tt'(p,q) between latent brain processes p and q at two time points t and t', while K_pq(t,t') is a smooth temporal kernel that connects events that occur at similar time points. In geospatial applications, the correlations A_'(p,q) can encode similarity between two latent ore functions p and q at two locations ,' ∈^2, for instance between cadmium and zinc concentrations <cit.>. The location kernel K_pq(, ') could encode a smooth spatial proximity. A conventional Kronecker kernel k(,',p,q) = A(p,q) K(,') would assume – in contrast – that (i) the same spatial proximity K(,') applies to all ore functions p, and (ii) two ore concentrations would correlate similarly independent of the location . 

We start forming the joint kernel by considering a non-stationary Gaussian kernel for the inputs x  <cit.>,

    K_pq(x,x') = √(2ℓ_pℓ_q/ℓ_p^2+ℓ_q^2)exp(-(x-x')^2/ℓ_p^2+ℓ_q^2),

which encodes specific lengthscales ł_1, …, ł_Q for each latent signal. 
The kernel within a single latent signal K_pp reduces into a standard Gaussian kernel, while the cross-covariance similarity K_pq measures similarity of two inputs with different associated lengthscales.

We base our construction of the mutually dependent covariance structure A_xx'(p,q) on Wishart processes. A Generalized Wishart Process (GWP) prior on a covariance matrix, that depends on a single variable x, is <cit.>

    A(x)    = ∑_r=1^ν L _r(x) _r(x)^T L^T ∼(V,ν,K_z),

where V = LL^T and all z_pr(x) ∼(0, K_z(x,x') )
are independent Gaussian processes for p = 1, …, Q and r = 1, …, ν. 
The kernel K_z determines the change of A(x) in the input space. From this formulation we define our joint kernel, such that we preserve the GWP marginal for A(x) by extending the GWP into cross-covariances of two variables, as

    A_xx'(p,q) = _p(x)^T _q(x'),

where we have for each element of _p(x) ∈^ν a GP prior (See Figure <ref>). With this choice the prior expectation of the covariance is the identity matrix.

The resulting covariance of u_p(x) is then a product of covariance between inputs at signals p and q, and a covariance between signals at inputs x and x'. 
This covariance can be seen marginally from two perspectives, 

    _̆p    ∼_N(, A(p,p) ∘ K_pp) 
    (̆x)    ∼_Q(, A(x) ∘ P),

where _̆p ∈^N is a single latent signal that follows a Normal distribution weighted by variances A(p,p), and (̆x) ∈^Q contains all Q latent signals at input x and follows a Normal distribution with generalized Wishart process prior, scaled by the matrix P_pq = √(2ℓ_pℓ_q/ℓ_p^2+ℓ_q^2). The element-wise, or Hadamard, product of x and y is denoted by x ∘ y.

The joint covariance over the concatenated column vector of all latent signals ∈̆ℝ^QN is a block matrix

    (,̆)̆   = ( Z_i Z_j^T ∘ K_ij)_i,j=1^N + _u = ^T ∘_Q + _u
 
where Z_i is a Q ×ν matrix, and = (Z_1, …, Z_N)^T. 

The Q × Q kernel K_ij = (K_pq(x_i,x_j))_p,q=1^Q gives the signal similarities at inputs x_i and x_j, the block matrix _Q = (K_ij)_i,j=1^N collects them into (N × N) blocks of kernel values, and finally the noise matrix is _u = ω_u^-1I⃗_QN, introducing the latent noise directly into the covariance of the $̆'s. The resulting joint input-output covariance^T ∘_Qconsists ofN ×Nblock matrices of sizeQ ×Q. See Figure <ref> for a visualisation. 
The kernel matrix_Qis positive semi-definite (PSD) as an outer product, and the Gaussian kernel is PSD as well <cit.>. The Hadamard product_Q ∘_Qretains this property. We refer to this kernel as the Wishart-Gibbs cross-covariance (WGCC).

The proposed latent correlation Gaussian process (LCGP) model is a flexible Bayesian regression model that simultaneously infers the latent signals and their mixing to match the output processes, while learning the underlying correlation structure of the latent space using the WGCC kernel. The latent correlations are parameterised by two terms that characterise the input and signal similarities with Gaussian and Wishart functions, respectively. A key feature of the model is the ability adaptively couple and decouple latent processes along the input space.








§ CLASSIFICATION WITH MULTIPLE OBSERVATIONS


We further suppose that we haveSobservations or samples^(s)(x)associated with a class label, or response,r^(s), and assume that all these observations share their latent space. We then learn separate latent functions^̆(s)for each sample, while keeping the mixing modelB(x), latent correlations(x)and_Q, and the noise precisionsω_fandω_ushared. The noiseless sample is then reconstructed as^(s)(x) = B(x) ^̆(s)(x),which results in the same likelihood as in eq. (<ref>).

We build a classifier in the latent signal space as
a Probit classification model over all latent signals^T ^̆(s)with Gaussian-Gamma priors, where∈^NQis a concatenated column vector of linear weights_p ∈^Nfor theQlatent signals. This allows us to reduce the data dimensionality for classification, asMcan potentially be very large. The classifier is then 

    r^(s)|, ^̆(s)   ∼(Φ(^T ^̆(s) + b)), 
    _p|λ_w    ∼(0, λ_w^-1),  λ_w ∼( α_w, β_w) 
    
        b|λ_b    ∼(0, λ_b^-1),  λ_b ∼(α_b, β_b),

where we index the observations withs, andandbare the classifier weights and bias, respectively. The Gaussian CDF is denoted byΦ(·). We additionally assume, for notational clarity, that all data are observed at the same input pointsx_1, …, x_N.  

Essentially, our model now has two likelihoods for the two types of data, one defined for the output data in eq. (<ref>) and one for the class labels related to the outputs in eq. (<ref>).









§ INFERENCE


 §.§ Variational Bayes

For inference in our Bayesian model we adopt the Variational Bayesian (VB) approach <cit.>, which is
based on maximising a lower bound on the log marginal likelihood of the data with respect to a distributionq(Θ), whereΘrepresents all model parameters. The lower bound is of an easier form than the true posterior distributionp(Θ|), where= (Y^(s),r^(s))_s=1^SandY^(s) ∈^M ×N. The lower bound is obtained
by Jensen's inequality

    log p()    = log∫ q(Θ) p(, Θ)/q(Θ) dΘ≥∫ q(Θ) logp(, Θ)/q(Θ) dΘ≡Ł(q)   .

Typically, a factorised approximationq(Θ) = ∏_i q(θ_i)is used, whereθ_iare some disjoint subsets of the variablesΘ. It can be
shown that the optimal solution that maximizesŁ(q)is

    q(θ_i) ∝exp(⟨log p(Θ, ) ⟩_θ_-i),

in which the expectation is taken with respect to all variables exceptθ_i.
The VB algorithm consists of iterating through updating each factorq(θ_i).



 §.§ VB for the LCGP classification model


We employ the following factorization 

    q(Θ)    = ∏_s q(^̆(s)) q(h^(s)) ∏_m q(_m) q(ω_f) q(,b) q(λ_w) q(λ_b) q(),

whereq(_m)factorizes the mixing matrixrow-wise.
Most factors have standard distributions, the update formulas are shown in Table <ref>. The VB inference procedure is summarised in Algorithm <ref>. LCGP can be run with or without the classification part of the model; without classification the parameters involved are ignored (see Figure <ref>).



Auxiliary variableshare introduced to make the variational inference tractable for Probit classification <cit.>,

    h |, ∼̆(^T+̆ b, 1).

Class labels depend on the sign ofh, i.e.r = +1ifh > 0. Integrating outhrecovers the Probit likelihood 

    p(r|,)̆ = ( r | Φ(^T+̆b)).

The posteriorq(h)is a truncated Gaussian <cit.>, which has analytical formulas for first and second moments.

Finally,q()is updated by optimising the lower boundŁ()with respect to. We optimiseŁ()using L-BFGS in whitened domain employing a change of variables= L⃗^-1with the Cholesky of the kernel_z = L⃗ L⃗^Tto make the optimization more efficient <cit.>, see the Supplementary for details.

Predictions to new inputsx^*can be made by applying standard GP formulas to obtain(x^*),(̆x^*)and(x^*)based on the optimized variational posteriorq(Θ). For new observation with unkown class labelr^(s^*), we can apply the update forq(^̆(s^*))without classification related terms.




§ RELATED WORK

In semiparametric latent factor models (SLFM) the signal(x) = B (̆x)overMoutputs is a linear combination ofQindependent latent Gaussian process signals(̆x)with a fixed mixing matrixB ∈^M ×Q, with appropriate hyperparameter learning <cit.>. A Gaussian process regression network (GPRN) <cit.> extends this model by considering a mixing matrixB(x)where each elementB_pi(x)is an independent Gaussian process alongx.

In geostatistics vector-valued regression with Gaussian processes is called cokriging<cit.>. In linear coregionalization models (LCM) latent Gaussian processes are mixed from latent signalsu_p(x)andu_q(x')that are independent. In contrast to SLFM, each signalu_p(x)is an additional mixture ofR_Qsignals with separate shared covariancesK_q(x,x'). In the intrinsic coregionalization model (ICM) only a single (Q=1) latent mixture with a single shared kernel is used, while in SLFM there are multiple latent singleton (R_Q=1) signals. In spatially varying LCMs (SVLCM) the mixing matrices are input-dependent, similar to GPRNs <cit.>. <cit.> used non-orthogonal latent signalsu_p(x)andu_q(x')with fixed covariances.

Multi-task Gaussian processes employ structured covariances that combine a task covariance with an input covariance. Simple Kronecker products between the covariances assume that task and input covariances are independent functions  <cit.>. This is computationally efficient <cit.>, but it does not take into account interactions between the tasks and inputs.

In Generalised Wishart Processes <cit.> an input-dependent covariance matrixΣ(x) = ∑_n=1^ν_n(x) _n(x)^Tis a sum ofνouter products. 
The random variablesz_ni(x) ∼(0, K(x,x'))are all independent Gaussian processes. Copula processes also describe dependencies of random variables by Gaussian processes <cit.>. In Bayesian nonparametric covariance regression, covariances of multiple predictors share a common dictionary of Gaussian processes <cit.>.

Finally, Gaussian process dynamical or state-space systems are a general class of discrete-time state-space models that combine the latent state into time-dependent outputs as Markov processes <cit.>. In Gaussian process factor analysis the outputs are described as factors that have GP priors, however not modeling the factor dependencies <cit.>.







§ EXPERIMENTS


In the first experiment we show that our model[Our Matlab implementation can be found on <https://github.com/sremes/wishart-gibbs-kernel>] can recover the true latent correlations in a simple simulated-data case, and compare our method with GPRNs, which is a state-of-the-art multi-output Gaussian process regression model <cit.>. We employ the mean-field variational inference implementation of GPRN by <cit.>. Second, we apply our method to the Jura geospatial dataset to elucidate latent ore concentration process couplings. Finally, we demonstrate our full modelling framework on an EEG single-trial classification task, outperforming state-of-the-art regularised LDA in classification and additionally recovering an interesting latent representation that we further evaluate in a simulation study.
Results from the experiments are summarised in Table <ref>.





 §.§ Simulated Data Experiments


  §.§.§ Wishart-Gibbs kernel in multi-output GP


We simulated a dataset that contains a clear switch in the coupling of three outputs in the middle of an interval[-1, 1]. A traditional Kronecker multi-output kernel of form= A ⊗K, withA = ∑_k _k_k^TandKa Gaussian kernel, cannot model this, but our proposed Wishart-Gibbs kernel can adapt to this switch point. The data and posterior fit with both kernels are shown in Figure <ref>. Our kernel obtains an MSE of 0.44, and Kronecker kernel 0.60, with a baseline of 1.92 (predicting zero).



  §.§.§ Recovering latent covariance with LCGP

We use simple toy data to show that we are able to recover known latent correlation structure.
We generated data with a varying number of latent componentsQ = 2, …, 5and amount of samplesS = 1, …, 20. The mixing matrix was binary such that one output maps to one latent variable. For simplicity, we only consider LCGP without the classifier.



To assess the accuracy, we measured the correlation between the elements of the true covariance matrix to the one estimated. With GPRN we computed the empirical covarianceΣ̂ = ∑_s ^̆(s)^̆(s)  ^Tof the latent variables. As the order of the recovered latent variables is not identifiable, we computed the correlations over all permutations and report the best. Rotations of the latent space are not accounted for, however. The results in  Figure <ref> show that our model can recover the true underlying latent covariance with high correlations.




 §.§ Jura


The Jura dataset[Data available at <https://sites.google.com/site/goovaertspierre/pierregoovaertswebsite/publications/book>.] 
consists of measurements of cadmium, nickel and zinc concentrations in a region of the Swiss Jura <cit.>.
For training we are given the concentrations measured at 259 locations and for validation the measurements at 100 additional locations.
We set hyperparameters for both our model and GPRN asℓ_u = 0.5andℓ_b = 1, 
and for our model the parameter for the latent correlation lengthscale toℓ_z = 1.
We learned the models withQ=2latent variables, which resulted in the best model performance. We report both the mean squared and absolute errors for the predicted concentrations 
in Table <ref>. Our model performs at the same level as the state-of-the-art competitor GPRN, with slightly better performance in absolute errors.

Figure <ref> shows the inferred model. The latent variables are 2D spatial surfaces on which the measurement points are indicated as black points. The two latent variables learn different geological processes that have an interesting two-pronged correlation pattern that indicates two kinds of negative correlations (the scatter plot). By explicitly modelling the latent covariance, we are able to see the regions of the input space that contribute to this pattern; the latent covariances indicate the combined covarianceC_pq(,') = A_'(p,q) K_pq(, '). The diagonal plots of Figure <ref>a show the variances of the two latent signals, while the off-diagonal covariance plot indicates the two-pronged negative correlation model between the geological processes. Finally, the mixing matrices of the two latent components reconstruct the three ore observation surfaces. 









 §.§ EEG


Our main motivation for developing the present model was in modelling EEG data. 
We demonstrate LCGP on data from a P300 study <cit.>, where
the subjects were shown either a target or non-target stimulus, specifically a green or a red LED 
flashing, respectively. The classification task is to classify the stimulus based on the brain measurements.
Additionally, we evaluate our modelling approach in a simulation study.



  §.§.§ Classification Results


We evaluate the classification performance using a Monte Carlo cross-validation scheme where in each fold we randomly sample training and test sets ofS=1000trials from the full dataset consisting in total of 7351 trials from 16 subjects. A single trial is the continuous voltage measurement ofM=19channels in an EEG cap for 800ms withN=89after filtering and downsampling the time series <cit.>. We report the average area under the ROC curve (AUC) statistic over 100 folds, and compare our method to the state-of-the-art regularised LDA method implemented in the BBCI toolbox <cit.>. Results in Table <ref> show that our method performs better than RLDA (p<0.05).

An example visualization of the model from one of the cross-validation folds is depicted in Figure <ref> for the three first latent signals. Panel (a) indicates the shared variances and covariances of the latent signals along time. The first and third latent signals have a monotonically increasing covariance coupling, while the first and second latent signals have a periodicity in the covariance. The average latent variables of the target and non-target trials are shown in panel (b). The third latent variable captures a strong dynamic between time points[0.3, 0.5], which coincides with the expected P300 activity approximately 300 ms after the stimulus representation. The first two variables show peaks also at approximately 300 ms. 
In general the positive trials have a remarkably different latent representations than the negative trials. Panel (c) shows the classifier weights_pfor the three latent signals with the average classification plotted. Finally, a subset of the EEG channels are shown in panel (d), highlighting the differences in the channel dynamics. In panel (e) the components are found to be discriminative also when plotted on the scalp map.



  §.§.§ Finding the Latent Correlations


In addition to the classification results, we evaluated our modelling approach in a simulation study to test whether we can find the latent correlations correctly using data that resembles the real EEG as closely as possible, but where we know the ground truth. To this end, we simulated datasets from the fitted models, and learn a new model on the simulated data.
We repeated this for varying number of latent variables (Q = 2,3,4) with 10 simulations done for each value ofQ. For each simulation we computed the empiricalp-value with the hypothesis that the accuracy of our model is greater than using randomly simulated covariances from the model. Thep-values from the simulations were combined using the Fisher's method <cit.>, with results reported in Table <ref>. We use again the correlation-based score for assessing the accuracy as with the toy data case.



§ DISCUSSION


The LCGP is a flexible framework for multi-task learning. We demonstrate in two experiments that our model can robustly learn the latent variable correlations. The model also achieves state-of-the-art performance in both regression and classification. The added modeling of the correlations of the underlying latent processes both improves model interpretability, and regularises the model especially with multiple observations. The novel Wishart-Gibbs cross-covariance kernel encodes mutually-dependent covariances between latent signals and inputs in a parameterised way without being too flexible.

In place of the non-stationary Gaussian kernel other non-stationary kernels are possible. <cit.> propose a class of non-stationary convolution kernels containing, for instance, a non-stationary Matérn kernel. For future work coupling the spectral kernels <cit.> with Wishart correlations is another highly interesting avenue for a general family of dependent, structured kernels. The mutually-dependent Hadamard kernel would also be interesting to study in context of structured multi-task learning to model dependent input-output relations.

 This work was supported by the Finnish Funding Agency for Innovation (project Re:Know) and Academy of Finland (COIN CoE, and grants 299915, 294238 and 292334). The research leading to this results has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement no 611570. We acknowledge the computational resources provided by the Aalto Science-IT project.

§ OPTIMISATION OF KERNEL PARAMETERS


The factorq()corresponding to the GP parameters of the proposed Hadamard kernel is updated by finding point estimates by maximizing the variational lower bound. The relevant part of the bound is given by

    Ł()    =  ∑_s ⟨log p(^̆(s)|) ⟩ + log p()
        = -(Slog|| + ∑_s ⟨^̆(s)T^-1^̆(s)⟩ +  ( ^T_z^-1 ))

where= ^T ∘_Q + , and its gradient by

    Ł_ij = ( [ ^-1(∑_s ⟨^̆(s)^̆(s)T⟩)^-1 - S^-1] _ij) -  [_z^-1]_ij

where_z = K_z ⊗I_Qis a block matrix of full size(QN ×QN), and

    _ij = (^T ∘_Q + )_ij   = (1⃗_ij^T + 1⃗_ij^T) ∘_Q.

The cost function in the whitened domain can be evaluated asŁ(L⃗ )and gradient as 

    Ł = L = L⃗^T Ł .

We can similarly optimize noise precisionω_uand lengthscalesℓ_u.