# Multi-output spatial gaussian process Matrix Positive Definite Issue

**URL:** <https://discourse.pymc.io/t/multi-output-spatial-gaussian-process-matrix-positive-definite-issue/17405>\
**Category:** v5\
**Created:** [October 26, 2025, 1:17pm UTC](https://discourse.pymc.io/t/multi-output-spatial-gaussian-process-matrix-positive-definite-issue/17405 "2025-10-26T13:17:36Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![kpjmcg](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/kpjmcg/32/6257_2.png) [@kpjmcg](https://discourse.pymc.io/u/kpjmcg)\
**Post date:** [October 26, 2025, 1:17pm UTC](https://discourse.pymc.io/t/multi-output-spatial-gaussian-process-matrix-positive-definite-issue/17405/1 "2025-10-26T13:17:36Z")

</div>

Hello,

I am trying to model a multi-output spatial gaussian process by adapting [this](https://gitlab.com/essicolo/blog_sm/-/blob/master/2022-03-14_spatial/2022-03-14_pymc3.ipynb) older model. I have 5 parameters I am observing at each location.

The current code looks like:

````py
> with pm.Model() as model:
> # GP
> ## Feature covariance
> length_scale = pm.HalfCauchy(“length_scale”, beta = 10)
> amplitude = pm.HalfCauchy(“amplitude”, beta = 10)
> cov_feature = amplitude \*\* 2 \* pm.gp.cov.Matern52(
> input_dim = X_mod.shape\[1\],
> ls = length_scale,
> active_dims = np.arange(1, X_mod.shape\[1\]) # all except index
> )
>
> ```
> ## Coregion covariance
> W = pm.Normal(
> "W", mu = 0, sigma = 5,
> shape = (n_mod, n_mod)
> )
> kappa = pm.HalfCauchy("kappa", beta = 10, shape = n_mod)
> coreg = pm.gp.cov.Coregion(
> input_dim = X_mod.shape[1],
> active_dims = [id_col], # only index
> kappa = kappa,
> W = W
> )
> 
> ## Combined covariance
> #jitter = 1e-6 # Add this for stability
> jitter = pm.HalfNormal("jitter", sigma=1e-3)
> cov_f = coreg * cov_feature + pm.gp.cov.WhiteNoise(sigma=jitter)
> 
> ## GP noise
> σ = pm.HalfCauchy("σ", beta = 10)
> 
> ## Gaussian process
> gp = pm.gp.Marginal(cov_func = cov_f)
> 
> ## Marginal likelihood
> y_ = gp.marginal_likelihood("y", X = X_mod, y = y_mod, sigma = σ)
> 
> ```

````

I get the sampling error:

```auto
SamplingError: Initial evaluation of model at starting point failed!
Starting values:
{'length_scale_log__': array(2.30258509), 'amplitude_log__': array(2.30258509), 'W': array([[0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0.]]), 'kappa_log__': array([2.30258509, 2.30258509, 2.30258509, 2.30258509, 2.30258509]), 'jitter_log__': array(-6.90775528), 'σ_log__': array(2.30258509)}

Logp initial evaluation results:
{'length_scale': np.float64(-1.14), 'amplitude': np.float64(-1.14), 'W': np.float64(-63.21), 'kappa': np.float64(-5.72), 'jitter': np.float64(-0.73), 'σ': np.float64(-1.14), 'y': np.float64(-inf)}
You can call `model.debug()` for more details.

```

model.debug() indicates the parameters of y evaluate such that they violate positive definite requirements.

My whole data data is structured like:

XY\_mod:

| | index | X | Y | variable | value | variable\_id |
| --- | --- | --- | --- | --- | --- | --- |
| 110 | 110 | -0.581407 | 0.192955 | Sa | -0.624682 | 0 |
| 583 | 16 | 1.217071 | -1.641736 | Smr1 | -0.179331 | 3 |
| 860 | 104 | -0.291302 | 0.563685 | Str | -1.192428 | 4 |
| 322 | 133 | -2.328757 | 2.099471 | Sal | -0.060770 | 1 |
| 328 | 139 | -0.433248 | 0.387877 | Sal | -0.931602 | 1 |
| 546 | 168 | -0.844582 | -0.495918 | Spk | -0.961197 | 2 |
| 310 | 121 | -0.778119 | 0.433171 | Sal | 1.240696 | 1 |
| 347 | 158 | -0.133310 | -0.842685 | Sal | -0.213934 | 1 |

```py
> y_mod = XY_m\[‘value’\].values 
>
> n_mod=5, 
>
> X_mod = XY_m\[\[“variable_id”\] + features\].values
>
> id_col=0

```

Any assistance would be greatly appreciated. I’m sure it’s something silly.

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [October 26, 2025, 5:42pm UTC](https://discourse.pymc.io/t/multi-output-spatial-gaussian-process-matrix-positive-definite-issue/17405/2 "2025-10-26T17:42:46Z")

</div>

I extracted the following runnable example from your post:

```py
import pandas as pd
import pymc as pm
import numpy as np

data = [[110, 110, -0.581407, 0.192955 ,'Sa', -0.624682, 0],
[583, 16, 1.217071, -1.641736 ,'Smr1', -0.179331, 3],
[860, 104, -0.291302, 0.563685 ,'Str', -1.192428, 4],
[322, 133, -2.328757, 2.099471 ,'Sal', -0.060770, 1],
[328, 139, -0.433248, 0.387877 ,'Sal', -0.931602, 1],
[546, 168, -0.844582, -0.49591 ,'Spk', -0.961197, 2],
[310, 121, -0.778119, 0.433171 ,'Sal', 1.240696,],
[347, 158, -0.133310, -0.84268 ,'Sal', -0.213934, 1]]

df = pd.DataFrame(data, columns=['id', 'index', "X", 'Y', 'variable', 'value', 'variable_id']).set_index('id')

y_mod = df['value'].values
X_mod = df[['variable_id', 'X', 'Y']].values
n_mod = 5
id_col =0

with pm.Model() as model:
    
    length_scale = pm.HalfCauchy("length_scale", beta = 10)
    amplitude = pm.HalfCauchy("amplitude", beta = 10)
    
    cov_feature = amplitude ** 2 * pm.gp.cov.Matern52(
        input_dim = X_mod.shape[1],
        ls = length_scale,
        active_dims = np.arange(1, X_mod.shape[1]) # all except index
    )

    ## Coregion covariance
    W = pm.Normal(
        "W", mu = 0, sigma = 5,
        shape = (n_mod, n_mod)
    )
    kappa = pm.HalfCauchy("kappa", beta = 10, shape = n_mod)
    coreg = pm.gp.cov.Coregion(
        input_dim = X_mod.shape[1],
        active_dims = [id_col], # only index
        kappa = kappa,
        W = W
    )

    ## Combined covariance
    #jitter = 1e-6 # Add this for stability
    jitter = pm.HalfNormal("jitter", sigma=1e-3)
    cov_f = coreg * cov_feature + pm.gp.cov.WhiteNoise(sigma=jitter)

    ## GP noise
    σ = pm.HalfCauchy("σ", beta = 10)

    ## Gaussian process
    gp = pm.gp.Marginal(cov_func = cov_f)

    ## Marginal likelihood
    y_hat = gp.marginal_likelihood("y", X = X_mod, y = y_mod, sigma = σ)
    idata = pm.sample()

```

This samples for me without issue on the following versions:

```py
import pymc as pm
import pytensor 

print(f'PyMC Version: {pm. __version__ }')
print(f'Pytensor Version: {pytensor. __version__ }')

PyMC Version: 5.26.0
Pytensor Version: 2.35.1

```