# Covariance kernels for random processes with multidimensional output do not output symmetric matrices

**URL:** <https://openturns.discourse.group/t/covariance-kernels-for-random-processes-with-multidimensional-output-do-not-output-symmetric-matrices/55>\
**Category:** Development\
**Tags:** kriging, distribution, bug\
**Created:** [November 6, 2020, 4:15pm UTC](https://openturns.discourse.group/t/covariance-kernels-for-random-processes-with-multidimensional-output-do-not-output-symmetric-matrices/55 "2020-11-06T16:15:58Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![josephmure](https://yyz2.discourse-cdn.com/free1/user_avatar/openturns.discourse.group/josephmure/32/114_2.png) [@josephmure](https://openturns.discourse.group/u/josephmure)\
**Post date:** [November 6, 2020, 4:15pm UTC](https://openturns.discourse.group/t/covariance-kernels-for-random-processes-with-multidimensional-output-do-not-output-symmetric-matrices/55/1 "2020-11-06T16:15:58Z")

</div>

I think there is a mistake in the way covariance kernels for random processes with multidimensional output are handled by the library.

In the following code block, I do not think `operator()` should return a `CovarianceMatrix` because that matrix is not necessarily symmetric.

> <https://github.com/openturns/openturns/blob/6de2a00118458512c4cbbf179648073c92a44d6d/lib/src/Base/Stat/openturns/TensorizedCovarianceModel.hxx#L60>

To prove this point, consider a stationary Gaussian process Y with 1D input and 1D output. Let k:\mathbb{R} \rightarrow \mathbb{R} be its (stationary) covariance kernel.

Now let h be a real number and let Z be the Gaussian process defined by Z(t) = Y(t+h).

The process (Y,Z)^\intercal is Gaussian, stationary, with 1D input and 2D output. Let K be its covariance kernel. For any \tau \in \mathbb{R},

K(\tau) = \begin{pmatrix} Cov(Y(0), Y(\tau)) & Cov(Y(0), Z(\tau)) \\ Cov(Z(0), Y(\tau)) & Cov(Z(0), Z(\tau)) \end{pmatrix} ;

K(\tau) = \begin{pmatrix} Cov(Y(0), Y(\tau)) & Cov(Y(0), Y(\tau+h)) \\ Cov(Y(h), Y(\tau)) & Cov(Y(h), Y(\tau+h)) \end{pmatrix} ;

K(\tau) = \begin{pmatrix} k(\tau) & k(\tau + h) \\ k(\tau-h) & k(\tau) \end{pmatrix} .

I will try to propose a fix shortly.

---

<div class="post-metadata">

**Author:** ![josephmure](https://yyz2.discourse-cdn.com/free1/user_avatar/openturns.discourse.group/josephmure/32/114_2.png) [@josephmure](https://openturns.discourse.group/u/josephmure)\
**Post date:** [November 23, 2020, 7:47pm UTC](https://openturns.discourse.group/t/covariance-kernels-for-random-processes-with-multidimensional-output-do-not-output-symmetric-matrices/55/2 "2020-11-23T19:47:17Z")

</div>

Let \boldsymbol{Y} = (Y\_1,...,Y\_d)^T be a stationary Gaussian Process with one-dimensional input and d-dimensional output. Let K be its (stationary) covariance kernel. For any \tau \in \mathbb{R}, the (i,j)-th element of K(\tau) is equal to Cov(Y\_i(0), Y\_j(\tau)). And the (j,i)-th element of K(-\tau) is

Cov(Y\_j(0), Y\_i(-\tau)) = Cov(Y\_i(-\tau), Y\_j(0)) = Cov(Y\_i(0), Y\_j(\tau)):

it is equal to the (i,j)-th element of K(\tau).  
Therefore the matrix K(-\tau) is the transpose of the matrix K(\tau).

In the `UserDefinedStationaryCovarianceModel` class, we only need the user to specify the matrices K(t\_k) for nonnegative timestamps t\_k, but we need to transpose the matrix K(t\_k) when the `operator()` is called on -t\_k.