diff --git a/README.md b/README.md
index e552e36a..26ea9c57 100644
--- a/README.md
+++ b/README.md
@@ -125,7 +125,7 @@ result, and prediction interfaces.
 
 | Formula surface    | Supported terms                                                                                      |
 | ------------------ | ---------------------------------------------------------------------------------------------------- |
-| Univariate smooths | `s(..., bs='bs')`, `cr`, `cs`, `cc`, `cp`, `ds`, `ps`, `tp`, `ts`                                    |
+| Univariate smooths | `s(..., bs='bs')`, `cr`, `cs`, `cc`, `cp`, `ds`, `gp`, `ps`, `tp`, `ts`                              |
 | Structured smooths | random effects `re`, factor smooths `fs`, sum-to-zero factor smooths `sz`                            |
 | Tensor products    | `te(...)` and `ti(...)` over supported numeric marginals                                             |
 | Parametric terms   | numeric and factor terms, supported interactions, intercept policies, and formula offsets            |
diff --git a/docs/generate_notebooks.py b/docs/generate_notebooks.py
index 9b160c09..ea428518 100644
--- a/docs/generate_notebooks.py
+++ b/docs/generate_notebooks.py
@@ -1514,6 +1514,7 @@ def gam_notebook() -> dict:
 | `cr`, `cs` | cubic regression spline; `cs` adds null-space shrinkage |
 | `cc` | cyclic cubic spline for periodic covariates |
 | `ps` | P-spline with difference penalties |
+| `gp` | low-rank Gaussian-process smooth with spherical, power-exponential, or Matérn covariance |
 | `tp`, `ts` | thin-plate regression spline; `ts` adds shrinkage |
 | `te(...)` | scale-invariant tensor product including main-effect directions |
 | `ti(...)` | tensor interaction with marginal main-effect directions removed |
@@ -1535,6 +1536,9 @@ def gam_notebook() -> dict:
     "shrinkage_cubic": GAM(formula="demand ~ s(temperature, bs='cs', k=10)"),
     "cyclic": GAM(formula="demand ~ s(hour, bs='cc', k=8)"),
     "p_spline": GAM(formula="demand ~ s(temperature, bs='ps', k=10)"),
+    "gaussian_process": GAM(
+        formula="demand ~ s(temperature, humidity, bs='gp', k=20, m=[3, 1.0])"
+    ),
     "thin_plate": GAM(formula="demand ~ s(temperature, bs='tp', k=10)"),
     "shrinkage_thin_plate": GAM(formula="demand ~ s(temperature, bs='ts', k=10)"),
     "tensor_surface": GAM(
diff --git a/docs/notebooks/01_gam.ipynb b/docs/notebooks/01_gam.ipynb
index 1442d1b3..8c93cec3 100644
--- a/docs/notebooks/01_gam.ipynb
+++ b/docs/notebooks/01_gam.ipynb
@@ -629,6 +629,7 @@
     "| `cr`, `cs` | cubic regression spline; `cs` adds null-space shrinkage |\n",
     "| `cc` | cyclic cubic spline for periodic covariates |\n",
     "| `ps` | P-spline with difference penalties |\n",
+    "| `gp` | low-rank Gaussian-process smooth with spherical, power-exponential, or Mat\u00e9rn covariance |\n",
     "| `tp`, `ts` | thin-plate regression spline; `ts` adds shrinkage |\n",
     "| `te(...)` | scale-invariant tensor product including main-effect directions |\n",
     "| `ti(...)` | tensor interaction with marginal main-effect directions removed |\n",
@@ -652,6 +653,7 @@
        "['cubic',\n",
        " 'cyclic',\n",
        " 'factor_smooth',\n",
+       " 'gaussian_process',\n",
        " 'p_spline',\n",
        " 'random_effect',\n",
        " 'shrinkage_cubic',\n",
@@ -675,6 +677,9 @@
     "    \"shrinkage_cubic\": GAM(formula=\"demand ~ s(temperature, bs='cs', k=10)\"),\n",
     "    \"cyclic\": GAM(formula=\"demand ~ s(hour, bs='cc', k=8)\"),\n",
     "    \"p_spline\": GAM(formula=\"demand ~ s(temperature, bs='ps', k=10)\"),\n",
+    "    \"gaussian_process\": GAM(\n",
+    "        formula=\"demand ~ s(temperature, humidity, bs='gp', k=20, m=[3, 1.0])\"\n",
+    "    ),\n",
     "    \"thin_plate\": GAM(formula=\"demand ~ s(temperature, bs='tp', k=10)\"),\n",
     "    \"shrinkage_thin_plate\": GAM(formula=\"demand ~ s(temperature, bs='ts', k=10)\"),\n",
     "    \"tensor_surface\": GAM(\n",
diff --git a/nampy/gam/compiler/factory.py b/nampy/gam/compiler/factory.py
index c1739629..8cfa2bfe 100644
--- a/nampy/gam/compiler/factory.py
+++ b/nampy/gam/compiler/factory.py
@@ -18,6 +18,7 @@
 from ..smooths.univariate.bs import DerivativeBSplineTerm1D
 from ..smooths.univariate.cr import CubicSplineTerm
 from ..smooths.univariate.ds import DuchonSplineTerm
+from ..smooths.univariate.gp import GaussianProcessTerm
 from ..smooths.univariate.ps import PSplineTerm1D
 from ..specs import LinearPredictorSpec, PenaltyGroupSpec, TermSpec
 from ..specs.smooth import (
@@ -27,6 +28,7 @@
     DerivativeBSplineSmoothSpec,
     DuchonSplineSmoothSpec,
     FactorSmoothInteractionSpec,
+    GaussianProcessSmoothSpec,
     PSplineSmoothSpec,
     RandomEffectSmoothSpec,
     ShapeConstrainedSmoothSpec,
@@ -207,6 +209,25 @@ def instantiate_term(term_like: TermSpec | Any):
             metadata=metadata,
         )
 
+    if isinstance(smooth_spec, GaussianProcessSmoothSpec):
+        return GaussianProcessTerm(
+            feature=features,
+            k=smooth_spec.k,
+            m=smooth_spec.m,
+            label=label,
+            term_id=term_like.term_id,
+            smoothing_id=smoothing_id,
+            by=by,
+            sp=smooth_spec.sp,
+            select=smooth_spec.select,
+            fixed=smooth_spec.fx,
+            constraint_mode=smooth_spec.constraint_mode,
+            pc=smooth_spec.pc,
+            knots=smooth_spec.knots,
+            xt=smooth_spec.xt,
+            metadata=metadata,
+        )
+
     if isinstance(smooth_spec, ShapeConstrainedSmoothSpec):
         if len(features) == 2:
             return BivariateShapePSplineTerm(
diff --git a/nampy/gam/compiler/linked_basis.py b/nampy/gam/compiler/linked_basis.py
index d8c930e6..ea9aff0e 100644
--- a/nampy/gam/compiler/linked_basis.py
+++ b/nampy/gam/compiler/linked_basis.py
@@ -92,6 +92,17 @@ def attach_shared_basis_metadata(predictor_specs, X, feature_names):
                 "different feature sets; mgcv 1.9-4 fails while constructing "
                 "the shared tensor basis."
             )
+        has_scalar_pc = any(
+            term.smooth_spec is not None
+            and str(term.smooth_spec.special).lower() == "s"
+            and getattr(term.smooth_spec, "pc", None) is not None
+            for term in group_terms
+        )
+        if has_scalar_pc and len(feature_tuples) > 1:
+            raise NotImplementedError(
+                "pc= is not supported across id-linked s() terms with different "
+                "feature sets; mgcv 1.9-4 fails while constructing the shared basis."
+            )
         for term in group_terms[1:]:
             _clone_linked_smooth_spec(base_term, term)
 
diff --git a/nampy/gam/smooths/__init__.py b/nampy/gam/smooths/__init__.py
index 44c5b6e4..07ab6dae 100644
--- a/nampy/gam/smooths/__init__.py
+++ b/nampy/gam/smooths/__init__.py
@@ -22,6 +22,7 @@
 from .univariate.bs import DerivativeBSplineTerm1D
 from .univariate.cr import CubicSplineTerm
 from .univariate.ds import DuchonSplineTerm
+from .univariate.gp import GaussianProcessTerm
 from .univariate.ps import PSplineTerm1D
 from .univariate.tp import ThinPlateSplineTerm
 
@@ -33,6 +34,7 @@
 bs = DerivativeBSplineTerm1D
 cr = cs = cc = CubicSplineTerm
 ds = DuchonSplineTerm
+gp = GaussianProcessTerm
 cp = ps = PSplineTerm1D
 tp = ts = ThinPlateSplineTerm
 fs = FSmoothInteractionTerm
@@ -60,6 +62,7 @@
     "build_penalty_definition",
     "CubicSplineTerm",
     "DuchonSplineTerm",
+    "GaussianProcessTerm",
     "DerivativeBSplineTerm1D",
     "PSplineTerm1D",
     "ThinPlateSplineTerm",
@@ -76,6 +79,7 @@
     "cs",
     "cc",
     "ds",
+    "gp",
     "cp",
     "ps",
     "tp",
diff --git a/nampy/gam/smooths/categorical/fs.py b/nampy/gam/smooths/categorical/fs.py
index f674866d..484dcb63 100644
--- a/nampy/gam/smooths/categorical/fs.py
+++ b/nampy/gam/smooths/categorical/fs.py
@@ -18,6 +18,7 @@
 from ..univariate.bs import DerivativeBSplineTerm1D
 from ..univariate.cr import CubicSplineTerm
 from ..univariate.ds import DuchonSplineTerm
+from ..univariate.gp import GaussianProcessTerm
 from ..univariate.ps import PSplineTerm1D
 from .categorical_utils import (
     as_object_1d,
@@ -109,7 +110,7 @@ def _build_base_smooth_term(
     Build the per-level base smooth used inside fs/sz.
 
     Supported base smooth classes in the current codebase:
-    bs, cr, cs, cc, cp, ds, ps, tp, ts
+    bs, cr, cs, cc, cp, ds, gp, ps, tp, ts
     """
     base_bs = str(base_bs).lower()
     metric_features = list(metric_features)
@@ -120,22 +121,23 @@ def _build_base_smooth_term(
     if mode == "fs" and base_bs in {"cs", "ts"}:
         raise NotImplementedError(_fs_full_rank_base_error(base_bs))
 
-    if len(metric_features) > 1 and base_bs not in {"ds", "tp", "ts"}:
+    if len(metric_features) > 1 and base_bs not in {"ds", "gp", "tp", "ts"}:
         raise NotImplementedError(
             f"Current {mode} implementation supports multivariate base smooths only "
-            f"for bs in {{'ds','tp','ts'}}, got base bs={base_bs!r}."
+            f"for bs in {{'ds','gp','tp','ts'}}, got base bs={base_bs!r}."
         )
 
     if xt_rest is not None and base_bs not in {
         "bs",
         "cp",
         "ds",
+        "gp",
         "ps",
         "tp",
         "ts",
     }:
         raise NotImplementedError(
-            "Extra xt options are currently only supported for bs/cp/ds/ps/tp/ts "
+            "Extra xt options are currently only supported for bs/cp/ds/gp/ps/tp/ts "
             "base smooths, "
             f"got xt={xt_rest!r} with base bs={base_bs!r}."
         )
@@ -220,6 +222,24 @@ def _build_base_smooth_term(
             metadata=metadata,
         )
 
+    if base_bs == "gp":
+        return GaussianProcessTerm(
+            feature=metric_features,
+            k=k,
+            m=outer_m,
+            label=label,
+            smoothing_id=None,
+            by=by,
+            sp=None,
+            select=bool(select),
+            fixed=bool(fixed),
+            constraint_mode=str(constraint_mode),
+            pc=None,
+            knots=knots,
+            xt=xt_rest,
+            metadata=metadata,
+        )
+
     if base_bs in {"tp", "ts"}:
         return make_smooth_term(
             base_bs,
@@ -242,7 +262,7 @@ def _build_base_smooth_term(
 
     raise NotImplementedError(
         f"Current {mode} implementation supports base bs in "
-        f"{{'bs','cr','cs','cc','cp','ds','ps','tp','ts'}}, got {base_bs!r}."
+        f"{{'bs','cr','cs','cc','cp','ds','gp','ps','tp','ts'}}, got {base_bs!r}."
     )
 
 
@@ -251,6 +271,8 @@ def _penalty_rank_from_base_term(base_term, basis_matrix, penalty_matrix) -> int
         return int(base_term._setup.ranks[0])
     if isinstance(base_term, DuchonSplineTerm):
         return int(base_term._setup.rank)
+    if isinstance(base_term, GaussianProcessTerm):
+        return int(base_term._setup.rank)
     if isinstance(base_term, PSplineTerm1D) and len(base_term.penalties) > 0:
         if str(base_term.basis_name).lower() == "cp":
             return int(base_term._setup.rank)
diff --git a/nampy/gam/smooths/tensor/marginals.py b/nampy/gam/smooths/tensor/marginals.py
index a25906f5..afdf833d 100644
--- a/nampy/gam/smooths/tensor/marginals.py
+++ b/nampy/gam/smooths/tensor/marginals.py
@@ -10,11 +10,12 @@
 from ..univariate.bs import DerivativeBSplineTerm1D
 from ..univariate.cr import CubicSplineTerm
 from ..univariate.ds import DuchonSplineTerm
+from ..univariate.gp import GaussianProcessTerm
 from ..univariate.ps import PSplineTerm1D
 from ..univariate.tp import ThinPlateSplineTerm
 
 TENSOR_MARGINAL_BASES = frozenset(
-    {"bs", "cr", "cs", "cc", "cp", "ds", "ps", "tp", "ts"}
+    {"bs", "cr", "cs", "cc", "cp", "ds", "gp", "ps", "tp", "ts"}
 )
 
 
@@ -130,6 +131,22 @@ def make_tensor_marginal_term(
             metadata=metadata,
         )
 
+    if basis == "gp":
+        return GaussianProcessTerm(
+            feature=marginal_features,
+            k=k,
+            m=m,
+            xt=xt,
+            label=str(feature),
+            smoothing_id=None,
+            by=None,
+            select=False,
+            fixed=False,
+            constraint_mode=constraint_mode,
+            knots=knots,
+            metadata=metadata,
+        )
+
     if basis in {"tp", "ts"}:
         return ThinPlateSplineTerm(
             feature=marginal_features,
diff --git a/nampy/gam/smooths/tensor/te.py b/nampy/gam/smooths/tensor/te.py
index 3562ee0d..146dbc56 100644
--- a/nampy/gam/smooths/tensor/te.py
+++ b/nampy/gam/smooths/tensor/te.py
@@ -114,6 +114,12 @@ def __init__(
 
         self._by_state = None
 
+    @property
+    def expected_linked_penalty_count(self):
+        if self.select:
+            return None
+        return int(np.sum(~np.asarray(self.fixed_flags, dtype=bool)))
+
     def fit(self, X, feature_names):
         marginal_shared_setups = self._linked_id_marginal_setups(self.feature)
         marginals, _, _ = build_tensor_marginal_terms(
diff --git a/nampy/gam/smooths/tensor/ti.py b/nampy/gam/smooths/tensor/ti.py
index aca048d9..51877be7 100644
--- a/nampy/gam/smooths/tensor/ti.py
+++ b/nampy/gam/smooths/tensor/ti.py
@@ -117,6 +117,12 @@ def __init__(
 
         self._by_state = None
 
+    @property
+    def expected_linked_penalty_count(self):
+        if self.select:
+            return None
+        return int(np.sum(~np.asarray(self.fixed_flags, dtype=bool)))
+
     def fit(self, X, feature_names):
         self._set_by_state(X, feature_names)
 
diff --git a/nampy/gam/smooths/univariate/__init__.py b/nampy/gam/smooths/univariate/__init__.py
index 87b0c1c7..290e0a69 100644
--- a/nampy/gam/smooths/univariate/__init__.py
+++ b/nampy/gam/smooths/univariate/__init__.py
@@ -1,12 +1,14 @@
 from .bs import DerivativeBSplineTerm1D
 from .cr import CubicSplineTerm
 from .ds import DuchonSplineTerm
+from .gp import GaussianProcessTerm
 from .ps import PSplineTerm1D
 from .tp import ThinPlateSplineTerm
 
 bs = DerivativeBSplineTerm1D
 cr = cs = cc = CubicSplineTerm
 ds = DuchonSplineTerm
+gp = GaussianProcessTerm
 cp = ps = PSplineTerm1D
 tp = ts = ThinPlateSplineTerm
 
@@ -14,7 +16,8 @@
     "DerivativeBSplineTerm1D",
     "CubicSplineTerm",
     "DuchonSplineTerm",
+    "GaussianProcessTerm",
     "PSplineTerm1D",
     "ThinPlateSplineTerm",
 ]
-__all__ += ["bs", "cr", "cs", "cc", "cp", "ds", "ps", "tp", "ts"]
+__all__ += ["bs", "cr", "cs", "cc", "cp", "ds", "gp", "ps", "tp", "ts"]
diff --git a/nampy/gam/smooths/univariate/gp.py b/nampy/gam/smooths/univariate/gp.py
new file mode 100644
index 00000000..458ca809
--- /dev/null
+++ b/nampy/gam/smooths/univariate/gp.py
@@ -0,0 +1,243 @@
+"""Gaussian-process smooth term (``bs='gp'``)."""
+
+from __future__ import annotations
+
+import numpy as np
+
+from ...constraints.absorption import (
+    fit_single_penalty_with_constraint_policy,
+    fit_single_penalty_with_setup_basis,
+)
+from ...penalties.algebra import penalty_rescale_factor, scale_penalty
+from ...splines.univariate.gp import (
+    build_gaussian_process_setup,
+    predict_gaussian_process,
+)
+from ..registry import register_smooth
+from ..smooth_base import BaseSmoothTerm, _resolve_feature, columns_as_float_matrix
+
+
+@register_smooth("gp")
+class GaussianProcessTerm(BaseSmoothTerm):
+    term_type = "smooth"
+    basis_name = "gp"
+    supports_tensor_marginal = False
+
+    def __init__(
+        self,
+        feature,
+        k=-1,
+        m=None,
+        label=None,
+        term_id=None,
+        smoothing_id=None,
+        by=None,
+        sp=None,
+        select=False,
+        fixed=False,
+        constraint_mode="auto",
+        pc=None,
+        knots=None,
+        xt=None,
+        null_penalty_tol=1e-10,
+        metadata=None,
+    ):
+        features = list(feature) if not isinstance(feature, (str, int)) else [feature]
+        super().__init__(
+            feature=features,
+            label=label or f"s({', '.join(map(str, features))})",
+            term_id=term_id,
+            smoothing_id=smoothing_id,
+            by=by,
+            sp=sp,
+            metadata=metadata,
+        )
+        self.k = int(k)
+        self.m = m
+        self.select = bool(select)
+        self.fixed = bool(fixed)
+        self.constraint_mode = str(constraint_mode).lower()
+        self.pc = pc
+        self.knots = knots
+        self.xt = xt
+        self.null_penalty_tol = float(null_penalty_tol)
+
+        if self.constraint_mode not in {"auto", "factor_by", "always", "never"}:
+            raise ValueError(
+                "constraint_mode must be one of "
+                "{'auto', 'factor_by', 'always', 'never'}."
+            )
+
+        self._feature_indices = None
+        self._feature_names = None
+        self._by_state = None
+        self._basis_train = None
+        self._penalties = None
+        self._setup = None
+
+    @property
+    def expected_linked_penalty_count(self):
+        return None if self.select else 1
+
+    def fit(self, X, feature_names):
+        feature_indices = []
+        resolved_names = []
+        for feature in self.feature:
+            index, name = _resolve_feature(feature, feature_names)
+            feature_indices.append(index)
+            resolved_names.append(name)
+
+        values = columns_as_float_matrix(X, feature_indices)
+        self._set_by_state(X, feature_names)
+        self._feature_indices = feature_indices
+        self._feature_names = resolved_names
+        self._set_resolved_features(resolved_names)
+
+        shared_X = self._linked_id_setup_matrix(feature_names)
+        if shared_X is not None:
+            setup_values = columns_as_float_matrix(shared_X, feature_indices)
+        else:
+            setup_values = values
+        self._setup = build_gaussian_process_setup(
+            setup_values,
+            k=self.k,
+            m=self.m,
+            knots=self.knots,
+            xt=self.xt,
+        )
+        setup_base = np.asarray(self._setup.basis_train, dtype=np.float64)
+        base = (
+            setup_base
+            if shared_X is None
+            else np.asarray(
+                predict_gaussian_process(values, self._setup), dtype=np.float64
+            )
+        )
+        raw_penalty = np.asarray(self._setup.penalty, dtype=np.float64)
+        penalty = scale_penalty(setup_base, raw_penalty)
+        self._set_penalty_rescale_factors(
+            [penalty_rescale_factor(setup_base, raw_penalty)]
+        )
+
+        if self.pc is not None:
+            constrained_basis, constrained_penalties, transform, _ = (
+                self._apply_point_constraint(
+                    base,
+                    [penalty],
+                    self.pc,
+                    feature_names=self._feature_names,
+                    point_basis_fn=lambda points: predict_gaussian_process(
+                        points, self._setup
+                    )[0],
+                    fixed=self.fixed,
+                )
+            )
+            self._basis_train = np.asarray(constrained_basis, dtype=np.float64)
+            self._penalties = constrained_penalties
+            self._record_constraint_result("pc", transform, absorbed_by="runtime")
+            return self
+
+        auto_constrain = bool(self._by_state.is_constant)
+        if shared_X is None:
+            result = fit_single_penalty_with_constraint_policy(
+                base,
+                penalty,
+                self._by_state,
+                constraint_mode=self.constraint_mode,
+                fixed=self.fixed,
+                auto_constrain_when=auto_constrain,
+            )
+        else:
+            result = fit_single_penalty_with_setup_basis(
+                base,
+                setup_base,
+                penalty,
+                self._by_state,
+                constraint_mode=self.constraint_mode,
+                fixed=self.fixed,
+                auto_constrain_when=auto_constrain,
+            )
+        self._basis_train = result.basis_train
+        self._penalties = result.penalties
+        self._record_constraint_result(
+            result.constraint_kind,
+            result.constraint_transform,
+            absorbed_by=(
+                "runtime" if result.constraint_transform is not None else None
+            ),
+        )
+        return self
+
+    def get_penalty_definitions(self):
+        self._require_fitted()
+        if not self.penalties:
+            return []
+        metadata = {
+            "term_type": self.term_type,
+            "basis_name": self.basis_name,
+            "feature": list(self.feature),
+            "label": self.label,
+            "by": self.by,
+            "by_name": self._by_state.feature_name,
+            "by_is_constant": bool(self._by_state.is_constant),
+            "constraint_mode": self.constraint_mode,
+            "constraint_kind": self.constraint_kind,
+            "pc": self.pc,
+            "knots": self.knots,
+            "xt": self.xt,
+            "m": self.m,
+            "gp_definition": np.asarray(self._setup.definition).tolist(),
+            "original_null_space_dim": self._setup.null_space_dim,
+            "fixed": bool(self.fixed),
+        }
+        selection_metadata = {**metadata, "is_selection_penalty": True}
+        metadata = self._penalty_metadata_with_scale(metadata, penalty_index=0)
+        return self._build_penalty_block(
+            self.penalties[0],
+            rank=int(self._setup.rank),
+            smooth_metadata=metadata,
+            selection_metadata=selection_metadata,
+        )
+
+    def transform_new(self, X_new):
+        self._require_fitted()
+        values = columns_as_float_matrix(X_new, self._feature_indices)
+        basis = predict_gaussian_process(values, self._setup)
+        return self._apply_constraint_transform_and_by(basis, X_new)
+
+    def tensor_marginal_fit_matrices(
+        self, *, centered=False, apply_np=False, x_train=None
+    ):
+        del apply_np, x_train
+        self._require_fitted()
+        setup_base = np.asarray(self._setup.basis_train, dtype=np.float64)
+        setup_penalty = np.asarray(self._setup.penalty, dtype=np.float64)
+        if centered:
+            if (
+                self._linked_id_setup() is not None
+                and self.constraint_transform is not None
+            ):
+                transform = np.asarray(self.constraint_transform, dtype=np.float64)
+                scaled = scale_penalty(setup_base, setup_penalty)
+                return (
+                    np.asarray(setup_base @ transform, dtype=np.float64),
+                    np.asarray(transform.T @ scaled @ transform, dtype=np.float64),
+                    None,
+                )
+            return super().tensor_marginal_fit_matrices(centered=True)
+        return setup_base, setup_penalty, None
+
+    def tensor_marginal_predict_matrix(
+        self, X_new, *, centered=False, np_transform=None
+    ):
+        if centered:
+            basis = np.asarray(self.transform_new(X_new), dtype=np.float64)
+        else:
+            values = columns_as_float_matrix(X_new, self._feature_indices)
+            basis = predict_gaussian_process(values, self._setup)
+        if np_transform is not None:
+            basis = basis @ np.asarray(np_transform, dtype=np.float64)
+        return np.asarray(basis, dtype=np.float64)
+
+
+__all__ = ["GaussianProcessTerm"]
diff --git a/nampy/gam/specs/__init__.py b/nampy/gam/specs/__init__.py
index 38c68879..b8d9d702 100644
--- a/nampy/gam/specs/__init__.py
+++ b/nampy/gam/specs/__init__.py
@@ -11,6 +11,7 @@
     DerivativeBSplineSmoothSpec,
     DuchonSplineSmoothSpec,
     FactorSmoothInteractionSpec,
+    GaussianProcessSmoothSpec,
     PSplineSmoothSpec,
     RandomEffectSmoothSpec,
     ShapeConstrainedSmoothSpec,
@@ -34,6 +35,7 @@
     "DerivativeBSplineSmoothSpec",
     "DuchonSplineSmoothSpec",
     "FactorSmoothInteractionSpec",
+    "GaussianProcessSmoothSpec",
     "PSplineSmoothSpec",
     "RandomEffectSmoothSpec",
     "ShapeConstrainedSmoothSpec",
diff --git a/nampy/gam/specs/build.py b/nampy/gam/specs/build.py
index c73b1958..a93add5c 100644
--- a/nampy/gam/specs/build.py
+++ b/nampy/gam/specs/build.py
@@ -477,8 +477,8 @@ def _fs_by_without_factor_feature_base_spec(smooth_spec: FactorSmoothInteraction
 
     if base_bs == "ps":
         kwargs["m"] = None if xt_rest is None else xt_rest.get("m", None)
-    elif base_bs in {"ds", "tp", "ts"}:
-        if base_bs == "ds":
+    elif base_bs in {"ds", "gp", "tp", "ts"}:
+        if base_bs in {"ds", "gp"}:
             kwargs["m"] = smooth_spec.m
         kwargs["xt"] = xt_rest
     elif xt_rest:
diff --git a/nampy/gam/specs/modeling.py b/nampy/gam/specs/modeling.py
index b881fd27..371ded20 100644
--- a/nampy/gam/specs/modeling.py
+++ b/nampy/gam/specs/modeling.py
@@ -42,7 +42,7 @@ def make_predictor_specs(model, feature_names, *, knots=None):
                     metadata={},
                 )
             )
-        elif basis in {"bs", "ds", "ps", "cp"}:
+        elif basis in {"bs", "ds", "gp", "ps", "cp"}:
             main_terms.append(
                 TermSpec(
                     kind="smooth",
@@ -105,7 +105,7 @@ def make_predictor_specs(model, feature_names, *, knots=None):
         else:
             raise NotImplementedError(
                 "Automatic main-effect construction currently supports "
-                "{'bs','cr','cs','cc','cp','ds','ps','tp','ts','re'}, "
+                "{'bs','cr','cs','cc','cp','ds','gp','ps','tp','ts','re'}, "
                 f"got {model.basis!r}."
             )
 
diff --git a/nampy/gam/specs/smooth.py b/nampy/gam/specs/smooth.py
index 6590e60d..a8164142 100644
--- a/nampy/gam/specs/smooth.py
+++ b/nampy/gam/specs/smooth.py
@@ -67,6 +67,15 @@ class DuchonSplineSmoothSpec(BaseSmoothSpec):
     pc: Any = None
 
 
+@dataclass(frozen=True)
+class GaussianProcessSmoothSpec(BaseSmoothSpec):
+    bs: str = "gp"
+    m: Any = None
+    xt: Any = None
+    constraint_mode: str = "auto"
+    pc: Any = None
+
+
 @dataclass(frozen=True)
 class ShapeConstrainedSmoothSpec(BaseSmoothSpec):
     """SCAM SCOP-spline specification for a named shape basis code."""
@@ -142,6 +151,7 @@ class TensorInteractionSmoothSpec(BaseSmoothSpec):
     CyclicCubicRegressionSmoothSpec,
     DerivativeBSplineSmoothSpec,
     DuchonSplineSmoothSpec,
+    GaussianProcessSmoothSpec,
     CubicShrinkageSmoothSpec,
     PSplineSmoothSpec,
     ShapeConstrainedSmoothSpec,
diff --git a/nampy/gam/specs/smooth_build.py b/nampy/gam/specs/smooth_build.py
index ef1da39e..40420e35 100644
--- a/nampy/gam/specs/smooth_build.py
+++ b/nampy/gam/specs/smooth_build.py
@@ -15,6 +15,7 @@
     DerivativeBSplineSmoothSpec,
     DuchonSplineSmoothSpec,
     FactorSmoothInteractionSpec,
+    GaussianProcessSmoothSpec,
     PSplineSmoothSpec,
     RandomEffectSmoothSpec,
     ShapeConstrainedSmoothSpec,
@@ -129,6 +130,21 @@ def _build_s_ds(opts) -> DuchonSplineSmoothSpec:
     )
 
 
+def _build_s_gp(opts) -> GaussianProcessSmoothSpec:
+    return GaussianProcessSmoothSpec(
+        special="s",
+        k=opts["k"],
+        fx=opts["fx"],
+        select=opts["select"],
+        sp=opts["sp"],
+        knots=opts["knots"],
+        m=opts["m"],
+        xt=opts["xt"],
+        constraint_mode=opts["constraint_mode"],
+        pc=opts["pc"],
+    )
+
+
 def _build_s_shape(opts) -> ShapeConstrainedSmoothSpec:
     return ShapeConstrainedSmoothSpec(
         special="s",
@@ -219,6 +235,7 @@ def _build_s_sz(opts) -> SumToZeroFactorSmoothSpec:
     "cp": _build_s_ps,
     "bs": _build_s_bs,
     "ds": _build_s_ds,
+    "gp": _build_s_gp,
     "tp": _build_s_tp,
     "ts": _build_s_ts,
     "re": _build_s_re,
@@ -318,6 +335,7 @@ def _is_vector_fx(fx) -> bool:
     "cr",
     "cs",
     "ds",
+    "gp",
     "ps",
     "tp",
     "ts",
@@ -337,7 +355,7 @@ def _dispatch_smooth_spec_from_options(opts) -> SmoothSpec:
             raise NotImplementedError(
                 f"pc= is not supported for s(..., bs={merged['bs']!r}); "
                 "point constraints are only supported for bs in "
-                "{'bs', 'cc', 'cp', 'cr', 'cs', 'ds', 'ps', 'tp', 'ts'}."
+                "{'bs', 'cc', 'cp', 'cr', 'cs', 'ds', 'gp', 'ps', 'tp', 'ts'}."
             )
         return builder(merged)
     if has_pc and special_key not in {"te", "ti"}:
@@ -535,10 +553,10 @@ def _default_k_for_smooth(kind, basis, features, default_k):
         # mgcv::te()/ti() default k to 5^d per marginal. The current
         # Python tensor surface supports one feature per marginal, so d = 1.
         return [5] * len(features)
-    if str(basis).lower() in {"bs", "ds", "tp", "ts"}:
+    if str(basis).lower() in {"bs", "ds", "gp", "tp", "ts"}:
         # mgcv/R/smooth.r::s() leaves k = -1. The basis constructor then
         # resolves its dimension-dependent default (TP/TS use M + 8/27/100;
-        # DS uses M + 10/30/100). A flat default here would be wrong in more
+        # DS uses M + 10/30/100; GP uses d + 1 + 10/30/100). A flat default here would be wrong in more
         # than one dimension.
         return -1
     return _default_k_for_basis(basis, default_k)
diff --git a/nampy/gam/splines/univariate/__init__.py b/nampy/gam/splines/univariate/__init__.py
index 1ab6b528..650f5474 100644
--- a/nampy/gam/splines/univariate/__init__.py
+++ b/nampy/gam/splines/univariate/__init__.py
@@ -26,6 +26,15 @@
     normalize_duchon_orders,
     predict_duchon_spline,
 )
+from .gp import (
+    GaussianProcessSetup,
+    build_gaussian_process_setup,
+    default_gp_k,
+    gp_kernel,
+    gp_polynomial_basis,
+    normalize_gp_definition,
+    predict_gaussian_process,
+)
 from .ps import (
     PSplineBasisSetup,
     bspline_design_matrix,
@@ -58,6 +67,13 @@
     "duchon_polynomial_basis",
     "normalize_duchon_orders",
     "predict_duchon_spline",
+    "GaussianProcessSetup",
+    "build_gaussian_process_setup",
+    "default_gp_k",
+    "gp_kernel",
+    "gp_polynomial_basis",
+    "normalize_gp_definition",
+    "predict_gaussian_process",
     "add_full_rank_shrinkage",
     "bspline_design_matrix",
     "cyclic_cubic_bd",
diff --git a/nampy/gam/splines/univariate/gp.py b/nampy/gam/splines/univariate/gp.py
new file mode 100644
index 00000000..428e3bd0
--- /dev/null
+++ b/nampy/gam/splines/univariate/gp.py
@@ -0,0 +1,278 @@
+"""Low-rank Gaussian-process smooth construction matching mgcv ``bs='gp'``."""
+
+from __future__ import annotations
+
+import warnings
+from dataclasses import dataclass
+
+import numpy as np
+
+from .ds import (
+    _duchon_setup_locations,
+    _duchon_unique_rows,
+    _normalize_duchon_knots,
+)
+from .tp import _top_eigensystem
+
+
+def normalize_gp_definition(m) -> tuple[np.ndarray, bool]:
+    """Normalize mgcv's ``m=(signed type, range, power)`` definition."""
+    if m is None:
+        values = []
+    elif np.isscalar(m):
+        values = [m]
+    else:
+        values = list(np.asarray(m, dtype=object).ravel())
+
+    missing = not values
+    if len(values) == 1:
+        try:
+            missing = bool(np.isnan(float(values[0])))
+        except (TypeError, ValueError):
+            missing = False
+
+    if missing:
+        signed_type = 3.0
+        stationary = False
+    else:
+        try:
+            first = float(values[0])
+            if not np.isfinite(first):
+                raise ValueError
+            gp_type = abs(int(np.rint(first)))
+        except (TypeError, ValueError, OverflowError) as exc:
+            raise ValueError("incorrect arguments to GP smoother") from exc
+        signed_type = float(np.sign(first) * gp_type)
+        stationary = bool(first < 0.0)
+
+    try:
+        rho = float(values[1]) if len(values) > 1 else -1.0
+        power = float(values[2]) if len(values) > 2 else 1.0
+    except (TypeError, ValueError) as exc:
+        raise ValueError("incorrect arguments to GP smoother") from exc
+    return np.asarray([signed_type, rho, power], dtype=np.float64), stationary
+
+
+def gp_polynomial_basis(x, definition):
+    """Port ``gpT``: constant tail for stationary GP, linear tail otherwise."""
+    values = np.asarray(x, dtype=np.float64)
+    if values.ndim == 1:
+        values = values.reshape(-1, 1)
+    definition = np.asarray(definition, dtype=np.float64).ravel()
+    constant = np.ones(values.shape[0], dtype=np.float64)
+    if definition[0] < 0.0:
+        return constant.reshape(-1, 1)
+    return np.asarray(np.column_stack([constant, values]), dtype=np.float64)
+
+
+def gp_kernel(x, knots, definition=None):
+    """Port ``gpE`` and return both the covariance and resolved definition."""
+    values = np.asarray(x, dtype=np.float64)
+    setup_knots = np.asarray(knots, dtype=np.float64)
+    if values.ndim == 1:
+        values = values.reshape(-1, 1)
+    if setup_knots.ndim == 1:
+        setup_knots = setup_knots.reshape(-1, 1)
+    if values.shape[1] != setup_knots.shape[1]:
+        raise ValueError("GP data and knots must have the same dimension.")
+
+    differences = values[:, None, :] - setup_knots[None, :, :]
+    distances = np.sqrt(np.sum(differences * differences, axis=2))
+    normalized, _ = normalize_gp_definition(definition)
+    signed_type, rho, power = map(float, normalized)
+    gp_type = abs(int(np.rint(signed_type)))
+    if np.isnan(rho) or np.isnan(power):
+        raise ValueError("missing value where TRUE/FALSE needed")
+    if rho <= 0.0:
+        rho = float(np.max(distances))
+    with np.errstate(divide="ignore", invalid="ignore", over="ignore"):
+        scaled = distances / rho
+
+    if gp_type not in {1, 2, 3, 4, 5} or power > 2.0 or power <= 0.0:
+        raise ValueError("incorrect arguments to GP smoother")
+
+    if gp_type == 1:
+        covariance = (1.0 - 1.5 * scaled + 0.5 * scaled**3) * (scaled <= 1.0)
+    elif gp_type == 2:
+        covariance = np.exp(-(scaled**power))
+    else:
+        exponential = np.exp(-scaled)
+        if gp_type == 3:
+            covariance = (1.0 + scaled) * exponential
+        elif gp_type == 4:
+            covariance = exponential + (scaled * exponential) * (
+                1.0 + scaled / 3.0
+            )
+        else:
+            covariance = exponential + (scaled * exponential) * (
+                1.0 + 0.4 * scaled + scaled**2 / 15.0
+            )
+
+    resolved = np.asarray([signed_type, rho, power], dtype=np.float64)
+    return np.asarray(covariance, dtype=np.float64), resolved
+
+
+def default_gp_k(dimension: int) -> int:
+    """Return the literal mgcv default ``d + 1 + c(10,30,100)[d]``."""
+    dimension = int(dimension)
+    if dimension < 1:
+        raise ValueError("Gaussian-process smooths require at least one covariate.")
+    if dimension > 3:
+        raise ValueError(
+            "An omitted k for bs='gp' is undefined upstream above three dimensions; "
+            "supply k explicitly."
+        )
+    return dimension + 1 + (10, 30, 100)[dimension - 1]
+
+
+def _parse_gp_xt(xt):
+    max_knots = 2000
+    seed = 1
+    if xt is None:
+        return max_knots, seed
+    if not isinstance(xt, dict):
+        raise NotImplementedError(
+            "For bs='gp', xt must be None or a dict with optional keys "
+            "{'max.knots', 'seed'}."
+        )
+    if xt.get("max.knots") is not None:
+        max_knots = int(xt["max.knots"])
+    if xt.get("seed") is not None:
+        seed = int(xt["seed"])
+    if max_knots < 1:
+        raise ValueError("For bs='gp', xt['max.knots'] must be positive.")
+    return max_knots, seed
+
+
+@dataclass
+class GaussianProcessSetup:
+    shift: np.ndarray
+    knots: np.ndarray
+    UZ: np.ndarray
+    definition: np.ndarray
+    null_space_dim: int
+    rank: int
+    bs_dim: int
+    basis_train: np.ndarray
+    penalty: np.ndarray
+    used_supplied_knots: bool
+    used_subsampling: bool
+
+
+def build_gaussian_process_setup(X, *, k=-1, m=None, knots=None, xt=None):
+    """Port ``smooth.construct.gp.smooth.spec`` and retain prediction state."""
+    values = np.asarray(X, dtype=np.float64)
+    if values.ndim == 1:
+        values = values.reshape(-1, 1)
+    if values.ndim != 2 or values.shape[1] < 1:
+        raise ValueError("GP smooth data must be a non-empty numeric matrix.")
+
+    n_obs, dimension = values.shape
+    raw_definition, stationary = normalize_gp_definition(m)
+    requested_k = int(k)
+
+    unique = _duchon_unique_rows(values)
+    if requested_k >= 0 and unique.shape[0] < requested_k:
+        raise ValueError(
+            "A term has fewer unique covariate combinations than specified "
+            "maximum degrees of freedom"
+        )
+
+    supplied = _normalize_duchon_knots(knots, dimension)
+    if supplied is not None and supplied.shape[0] == 0:
+        supplied = None
+    if supplied is not None and supplied.shape[0] > n_obs:
+        warnings.warn(
+            "more knots than data in an ms term: knots ignored.",
+            stacklevel=2,
+        )
+        supplied = None
+
+    shift = np.mean(values, axis=0)
+    if supplied is not None:
+        supplied = supplied - shift[None, :]
+    max_knots, seed = _parse_gp_xt(xt)
+    setup_knots, used_subsampling = _duchon_setup_locations(
+        values,
+        shift,
+        supplied,
+        max_knots=max_knots,
+        seed=seed,
+    )
+
+    covariance, definition = gp_kernel(setup_knots, setup_knots, raw_definition)
+    bs_dim = requested_k if requested_k >= 0 else default_gp_k(dimension)
+    if bs_dim < dimension + 2:
+        bs_dim = dimension + 2
+        warnings.warn("basis dimension reset to minimum possible", stacklevel=2)
+
+    null_space_dim = 1 if stationary else dimension + 1
+    rank = int(bs_dim - null_space_dim)
+    n_knots = int(setup_knots.shape[0])
+    if n_knots < rank:
+        raise ValueError(
+            "Gaussian-process smooth requires at least as many knot locations as "
+            "penalized basis coefficients."
+        )
+
+    penalty = np.zeros((bs_dim, bs_dim), dtype=np.float64)
+    if rank < n_knots:
+        eigenvalues, eigenvectors = _top_eigensystem(
+            covariance,
+            rank,
+            tolerance_exponent=0.5,
+        )
+        penalty[:rank, :rank] = np.diag(eigenvalues)
+    else:
+        eigenvectors = np.eye(rank, dtype=np.float64)
+        penalty[:rank, :rank] = covariance
+
+    setup = GaussianProcessSetup(
+        shift=np.asarray(shift, dtype=np.float64),
+        knots=np.asarray(setup_knots, dtype=np.float64),
+        UZ=np.asarray(eigenvectors, dtype=np.float64),
+        definition=np.asarray(definition, dtype=np.float64),
+        null_space_dim=int(null_space_dim),
+        rank=int(rank),
+        bs_dim=int(bs_dim),
+        basis_train=np.zeros((n_obs, bs_dim), dtype=np.float64),
+        penalty=np.asarray(penalty, dtype=np.float64),
+        used_supplied_knots=bool(supplied is not None),
+        used_subsampling=bool(used_subsampling),
+    )
+    setup.basis_train = predict_gaussian_process(values, setup)
+    return setup
+
+
+def predict_gaussian_process(X_new, setup: GaussianProcessSetup):
+    """Port ``Predict.matrix.gp.smooth`` including its knot-sized chunks."""
+    values = np.asarray(X_new, dtype=np.float64)
+    if values.ndim == 1:
+        values = values.reshape(-1, 1)
+    shifted = values - setup.shift[None, :]
+    n_obs = int(shifted.shape[0])
+    n_knots = int(setup.knots.shape[0])
+
+    def _block(block):
+        covariance, _ = gp_kernel(block, setup.knots, setup.definition)
+        tail = gp_polynomial_basis(block, setup.definition)
+        return np.asarray(np.column_stack([covariance @ setup.UZ, tail]))
+
+    if n_obs <= n_knots:
+        return np.asarray(_block(shifted), dtype=np.float64)
+    out = np.empty((n_obs, setup.bs_dim), dtype=np.float64)
+    for start in range(0, n_obs, n_knots):
+        stop = min(start + n_knots, n_obs)
+        out[start:stop, :] = _block(shifted[start:stop, :])
+    return out
+
+
+__all__ = [
+    "GaussianProcessSetup",
+    "build_gaussian_process_setup",
+    "default_gp_k",
+    "gp_kernel",
+    "gp_polynomial_basis",
+    "normalize_gp_definition",
+    "predict_gaussian_process",
+]
diff --git a/pyproject.toml b/pyproject.toml
index de0443e1..67e09538 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -208,6 +208,7 @@ markers = [
     "smooth_bs: tests covering integrated-derivative B-spline smooths",
     "smooth_cp: tests covering cyclic P-spline smooths",
     "smooth_ds: tests covering Duchon regression spline smooths",
+    "smooth_gp: tests covering Gaussian-process smooths",
     "smooth_ps: tests covering P-spline smooths",
     "smooth_tp: tests covering thin plate smooths",
     "smooth_ts: tests covering shrinkage thin plate smooths",
diff --git a/tests/SUBSYSTEM_COVERAGE.md b/tests/SUBSYSTEM_COVERAGE.md
index d5f3c782..6232deb2 100644
--- a/tests/SUBSYSTEM_COVERAGE.md
+++ b/tests/SUBSYSTEM_COVERAGE.md
@@ -15,7 +15,7 @@ local development notes rather than duplicated here.
 | Subsystem | Primary owner(s) | Primary tests | Notes |
 | --- | --- | --- | --- |
 | Formula/spec parsing | `nampy/gam/formula/`, `nampy/gam/specs/` | `tests/parity/test_mgcv_formula_parse_parity.py` | Direct formula parity vs `mgcv`. |
-| Smooth constructors / raw basis owners | `nampy/gam/smooths/`, `nampy/gam/splines/` | `tests/smooths/test_mgcv_raw_constructor_parity.py`, `tests/smooths/test_mgcv_smoothcon_parity.py`, `tests/parity/test_gam_spec_build_owner_contracts.py`, `tests/parity/test_mgcv_cp_combinations_parity.py`, `tests/parity/test_mgcv_bs_combinations_parity.py`, `tests/parity/test_mgcv_ds_combinations_parity.py` | Basis, penalty, constructor, and smoothCon surfaces, including cyclic P-splines (`cp`), integrated-derivative B-splines (`bs`), and multivariate Duchon splines (`ds`) across prediction, selection, linked bases, tensor/factor-smooth combinations, and the upstream tensor-`m` wrong-length warning and zero fallback. |
+| Smooth constructors / raw basis owners | `nampy/gam/smooths/`, `nampy/gam/splines/` | `tests/smooths/test_mgcv_raw_constructor_parity.py`, `tests/smooths/test_mgcv_smoothcon_parity.py`, `tests/parity/test_gam_spec_build_owner_contracts.py`, `tests/parity/test_mgcv_cp_combinations_parity.py`, `tests/parity/test_mgcv_bs_combinations_parity.py`, `tests/parity/test_mgcv_ds_combinations_parity.py`, `tests/parity/test_mgcv_gp_combinations_parity.py` | Basis, penalty, constructor, and smoothCon surfaces, including cyclic P-splines (`cp`), integrated-derivative B-splines (`bs`), multivariate Duchon splines (`ds`), and five-family stationary/nonstationary Gaussian-process smooths (`gp`) across prediction, selection, linked bases, tensor/factor-smooth combinations, and the upstream tensor-`m` wrong-length warning and zero fallback. |
 | `pc=` / linked `id=` routing | smooth metadata + linked basis owners | `tests/smooths/test_mgcv_pc_id_parity.py` | Localizes shared-smoothing and point-constraint issues, including `te`/`ti` point constraints. |
 | Design / pre-fit assembly | `nampy/gam/compiler/`, `nampy/gam/fit/penalized_system.py` | `tests/optimization/test_mgcv_gam_setup_assembly_parity.py`, `tests/optimization/test_mgcv_preoptimization_blocks_parity.py`, `tests/optimization/test_mgcv_preoptimization_reparam_parity.py` | Setup, blocks, and reparameterization parity, including one global shared-component block with overlapping linear-predictor indices. |
 | Term wrapping / by-variable / offset routing | predictor wrapping + compiled term owners | `tests/optimization/test_gam_term_wrapping_owner_contracts.py` | Localizes wrapped predictor blocks, offset routing, and general-family block ownership before broader prediction parity. |
@@ -47,7 +47,7 @@ closes the former seven-stage under-tested backlog.
 | Pipeline stage | Combination coverage now owned |
 | --- | --- |
 | 1. Formula parsing and canonical specs | Formula lists/shared components, transformed covariates, supported numeric/factor interactions, factor-by smooths, formula and fit offsets, intercept policies, weights, knots, `min_sp`, `drop_intercept`, and tensor-`m` warning/fallback behavior; the combined interaction recipe is rebuilt on newdata and compared with a committed `mgcv` reference fixture. |
-| 2. Runtime terms and low-level bases | Train/newdata pairing for every supported basis (`bs/cr/cs/cc/cp/ds/ps/tp/ts/re/fs/sz/te/ti`), including univariate, multivariate Duchon, and structured boundary response/SE parity. Row-permutation contracts cover linked and identified-FS terms, PS/TP/TS bases, `te`/`ti`, SZ, and factor-by smooths. Non-unique constructor representations use column-space projectors and penalized response operators instead of arbitrary coefficient orientation. |
+| 2. Runtime terms and low-level bases | Train/newdata pairing for every supported basis (`bs/cr/cs/cc/cp/ds/gp/ps/tp/ts/re/fs/sz/te/ti`), including univariate, multivariate Duchon/GP, and structured boundary response/SE parity. Row-permutation contracts cover linked and identified-FS terms, PS/TP/TS bases, `te`/`ti`, SZ, and factor-by smooths. Non-unique constructor representations use column-space projectors and penalized response operators instead of arbitrary coefficient orientation. |
 | 3. Constructed/wrapped terms | Numeric-by plus linked `id=`, factor-by plus linked `id=`, tensor-by, and mixed fixed/free/select penalty ownership, numerical coefficient-map composition, and fitted response/SE parity. |
 | 4. Predictor/model compilation | Multi-predictor layouts with unequal intercept and offset policies, overlapping shared-component coefficient indices, and three-term linked, fixed/free, select, and rank-deficient assembly; supported `gaulss` and `gammals` layouts are also fitted and compared with `mgcv`. |
 | 5. Side conditions and identifiability | Repeated, nested, reverse-formula-order, tensor/main-effect, three-way, both identified near-rank regimes, zero-width, no-intercept, ordered factor-by, linked, general-family, two-predictor, SZ, and exempt random/factor smooth cases. Deletion rank and behavior replace raw pivot-column identity where QR/eigen choices are non-unique. |
diff --git a/tests/TAXONOMY.md b/tests/TAXONOMY.md
index d5a6a46c..c86ea87d 100644
--- a/tests/TAXONOMY.md
+++ b/tests/TAXONOMY.md
@@ -12,7 +12,7 @@ The GAM test suite is intentionally overlapping. The goal is fast subset runs an
 - `tests/`: shared helpers, marker inference, taxonomy registry, static reference fixtures, and parity-generation R scripts
 
 ## Taxonomy Axes
-- `smooth_<name>`: `bs`, `cr`, `cs`, `cc`, `cp`, `ds`, `ps`, `tp`, `ts`, `te`, `ti`, `fs`, `sz`, `re`
+- `smooth_<name>`: `bs`, `cr`, `cs`, `cc`, `cp`, `ds`, `gp`, `ps`, `tp`, `ts`, `te`, `ti`, `fs`, `sz`, `re`
 - `family_<name>`: `gaussian`, `binomial`, `poisson`, `gamma`, `negbin`, `gaulss`, `gammals`, `general`
 - `method_<name>`: `fixed`, `reml`, `ml`, `laml`, `gcv`, `ubre`
 - `link_<name>`: `identity`, `log`, `inverse`, `logit`, `probit`, `cloglog`, `cauchit`, `sqrt`
diff --git a/tests/_taxonomy_registry.py b/tests/_taxonomy_registry.py
index ebe8d7b8..ae41327c 100644
--- a/tests/_taxonomy_registry.py
+++ b/tests/_taxonomy_registry.py
@@ -20,6 +20,7 @@ def _leaf(leaf_id: str, *nodeid_parts: str) -> LeafCoverageExpectation:
     "cc": "smooth_cc",
     "cp": "smooth_cp",
     "ds": "smooth_ds",
+    "gp": "smooth_gp",
     "ps": "smooth_ps",
     "tp": "smooth_tp",
     "ts": "smooth_ts",
@@ -143,6 +144,11 @@ def _leaf(leaf_id: str, *nodeid_parts: str) -> LeafCoverageExpectation:
         "tests/smooths/test_mgcv_raw_constructor_parity.py",
         "tests/smooths/test_mgcv_smoothcon_parity.py",
     ),
+    "smooth_gp": (
+        "tests/parity/test_mgcv_gp_combinations_parity.py",
+        "tests/smooths/test_mgcv_raw_constructor_parity.py",
+        "tests/smooths/test_mgcv_smoothcon_parity.py",
+    ),
     "smooth_cs": (
         "tests/smooths/test_mgcv_pc_id_parity.py",
         "tests/smooths/test_mgcv_raw_constructor_parity.py",
diff --git a/tests/mgcv_invariant_policy.py b/tests/mgcv_invariant_policy.py
index 581a3c89..097555bd 100644
--- a/tests/mgcv_invariant_policy.py
+++ b/tests/mgcv_invariant_policy.py
@@ -309,6 +309,17 @@ def _canonicalize_duchon_raw_state(state):
     return state
 
 
+def _canonicalize_gp_raw_state(state):
+    extra = state["extra"]
+    extra.pop("used_supplied_knots", False)
+    extra.pop("used_subsampling", False)
+    extra.pop("pure_knot", False)
+    state["S"] = [penalty_spectrum(S) for S in state["S"]]
+    state["X"] = matrix_self_gram(state["X"])
+    extra["UZ"] = stable_column_space_projector(extra["UZ"])
+    return state
+
+
 def _canonicalize_cs_raw_state(state):
     state["S"] = [penalty_spectrum(S) for S in state["S"]]
     return state
@@ -365,6 +376,8 @@ def canonicalize_raw_representation_state(state: dict[str, Any]) -> dict[str, An
         return _canonicalize_tprs_raw_state(state)
     if class_name == "duchon.spline":
         return _canonicalize_duchon_raw_state(state)
+    if class_name == "gp.smooth":
+        return _canonicalize_gp_raw_state(state)
     if class_name == "fs.interaction":
         return _canonicalize_fs_raw_state(state)
     if class_name == "sz.interaction":
diff --git a/tests/mgcv_parity_utils.py b/tests/mgcv_parity_utils.py
index d2e79681..85a4c65f 100644
--- a/tests/mgcv_parity_utils.py
+++ b/tests/mgcv_parity_utils.py
@@ -940,10 +940,11 @@ def _normalize_snapshot_payload(payload):
     if cached is not None:
         return _normalize_snapshot_payload(cached)
 
+    formula_for_r = _normalize_python_formula_text(formula)
     if optimizer is not None:
         result = _run_mgcv_snapshot_single(
             data,
-            formula,
+            formula_for_r,
             family_token,
             method,
             select=select,
@@ -954,7 +955,7 @@ def _normalize_snapshot_payload(payload):
         try:
             result = _run_mgcv_snapshot_batched(
                 data,
-                formula,
+                formula_for_r,
                 family_token,
                 method,
                 select=select,
@@ -963,7 +964,7 @@ def _normalize_snapshot_payload(payload):
         except Exception:
             result = _run_mgcv_snapshot_single(
                 data,
-                formula,
+                formula_for_r,
                 family_token,
                 method,
                 select=select,
@@ -1631,6 +1632,12 @@ def _run_mgcv_raw_constructor(
       shift = pack_vector(sm$shift, "numeric"),
       p_order = pack_vector(sm$p.order, "numeric")
     ),
+    "gp.smooth" = list(
+      knt = pack_matrix(sm$knt),
+      UZ = pack_matrix(sm$UZ),
+      shift = pack_vector(sm$shift, "numeric"),
+      gp_defn = pack_vector(sm$gp.defn, "numeric")
+    ),
     "random.effect" = list(
       C = pack_constraint(sm$C),
       random = isTRUE(sm$random),
diff --git a/tests/parity/test_mgcv_gp_combinations_parity.py b/tests/parity/test_mgcv_gp_combinations_parity.py
new file mode 100644
index 00000000..98d7806c
--- /dev/null
+++ b/tests/parity/test_mgcv_gp_combinations_parity.py
@@ -0,0 +1,293 @@
+"""Integrated parity coverage for Gaussian-process smooths (``bs='gp'``)."""
+
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from nampy.gam import GAM
+from nampy.gam.splines.univariate.gp import build_gaussian_process_setup
+from tests.mgcv_parity_utils import (
+    _assert_basic_mgcv_parity,
+    _fit_nampy_model,
+    _fit_nampy_snapshot,
+    _run_mgcv_snapshot,
+)
+
+
+def _gp_data(seed=301, n=180):
+    rng = np.random.default_rng(seed)
+    x0 = rng.uniform(-2.0, 2.0, size=n)
+    x1 = rng.uniform(-1.5, 1.5, size=n)
+    x2 = rng.uniform(-1.8, 1.8, size=n)
+    x3 = rng.uniform(-1.2, 1.7, size=n)
+    z = 0.8 + rng.uniform(-0.4, 0.7, size=n)
+    f = np.asarray(["a", "b", "c"], dtype=object)[np.arange(n) % 3]
+    f1 = np.asarray(["u", "v"], dtype=object)[np.arange(n) % 2]
+    y = (
+        0.2
+        + z * np.sin(1.2 * x0)
+        + 0.3 * x1**2
+        - 0.25 * np.cos(x2)
+        + 0.15 * x3
+        + 0.2 * (f == "b")
+        - 0.15 * (f1 == "v")
+        + rng.normal(scale=0.12, size=n)
+    )
+    return pd.DataFrame(
+        {
+            "y": y,
+            "x0": x0,
+            "x1": x1,
+            "x2": x2,
+            "x3": x3,
+            "z": z,
+            "f": f,
+            "f1": f1,
+        }
+    )
+
+
+def _assert_snapshot_fit(actual, expected, *, atol=3e-7):
+    for key in ("response", "link"):
+        np.testing.assert_allclose(
+            actual["predictions"][key],
+            expected["predictions"][key],
+            atol=atol,
+            rtol=atol,
+        )
+    np.testing.assert_allclose(
+        actual["fit"]["edf_total"],
+        expected["fit"]["edf_total"],
+        atol=atol,
+        rtol=atol,
+    )
+
+
+def test_gp_numeric_by_select_true_matches_mgcv():
+    data = _gp_data(seed=302)
+    formula = 'y ~ s(x0, x1, by=z, bs="gp", k=10, m=[4,.8])'
+    actual = _fit_nampy_snapshot(data, formula, "gaussian", "REML", select=True)
+    expected = _run_mgcv_snapshot(data, formula, "gaussian", "REML", select=True)
+    assert len(actual["fit"]["smoothing_params"]) == 2
+    _assert_basic_mgcv_parity(
+        actual,
+        expected,
+        pred_atol=5e-7,
+        pred_rtol=5e-7,
+        sp_log_atol=8e-6,
+        criterion_atol=3e-7,
+    )
+
+
+def test_gp_stationary_select_has_no_extra_null_penalty_and_matches_mgcv():
+    data = _gp_data(seed=303)
+    formula = 'y ~ s(x0, x1, bs="gp", k=10, m=[-5,.9])'
+    model = _fit_nampy_model(data, formula, "gaussian", "REML", select=True)
+    assert len(model.smoothing_params) == 1
+    _assert_basic_mgcv_parity(
+        model.parity_snapshot(X=data, include_covariances=True),
+        _run_mgcv_snapshot(data, formula, "gaussian", "REML", select=True),
+        pred_atol=5e-7,
+        pred_rtol=5e-7,
+        sp_log_atol=8e-6,
+        criterion_atol=3e-7,
+    )
+
+
+def test_gp_factor_by_fixed_sp_matches_mgcv():
+    data = _gp_data(seed=304)
+    formula = 'y ~ s(x0, x1, by=f, bs="gp", k=10, m=[2,.7,1.4], sp=.7)'
+    actual = _fit_nampy_snapshot(data, formula, "gaussian", "fixed")
+    expected = _run_mgcv_snapshot(data, formula, "gaussian", "REML")
+    _assert_snapshot_fit(actual, expected)
+
+
+def test_gp_linked_terms_clone_first_definition_but_retain_each_xt():
+    data = _gp_data(seed=305, n=150)
+    formula = (
+        'y ~ s(x0, x1, bs="gp", k=10, m=[4,.8], '
+        'xt={"max.knots":14,"seed":7}, id="shared_gp", sp=[.7,.9])'
+        ' + s(x2, x3, bs="gp", k=12, m=[1,.5], '
+        'xt={"max.knots":16,"seed":9}, id="shared_gp")'
+    )
+    model = _fit_nampy_model(data, formula, "gaussian", "fixed", select=True)
+    assert len(model.smoothing_params) == 2
+    runtimes = [
+        term.predict_fn.__self__
+        for term in model.gam_result_.compiled_model.compiled_terms
+        if term.term_type == "smooth"
+    ]
+    assert [runtime._setup.knots.shape[0] for runtime in runtimes] == [14, 16]
+    for runtime in runtimes:
+        np.testing.assert_allclose(runtime._setup.definition, [4.0, 0.8, 1.0])
+    _assert_snapshot_fit(
+        model.parity_snapshot(X=data, include_covariances=True),
+        _run_mgcv_snapshot(data, formula, "gaussian", "REML", select=True),
+        atol=8e-7,
+    )
+
+
+def test_gp_point_constraint_and_fixed_basis_match_mgcv():
+    data = _gp_data(seed=306)
+    formula = 'y ~ s(x0, x1, bs="gp", k=10, m=[3,.8], pc=[.2,-.3], sp=.8)'
+    actual = _fit_nampy_snapshot(data, formula, "gaussian", "fixed")
+    expected = _run_mgcv_snapshot(data, formula, "gaussian", "REML")
+    _assert_snapshot_fit(actual, expected)
+
+    fixed_formula = 'y ~ s(x2, x3, bs="gp", k=10, m=[-3,.8], fx=True)'
+    model = _fit_nampy_model(data, fixed_formula, "gaussian", "fixed", select=True)
+    assert model.gam_result_.compiled_model.compiled_penalties == ()
+    _assert_snapshot_fit(
+        model.parity_snapshot(X=data, include_covariances=True),
+        _run_mgcv_snapshot(
+            data, fixed_formula, "gaussian", "REML", select=True
+        ),
+    )
+
+
+@pytest.mark.parametrize(
+    "formula",
+    [
+        'y ~ te(x0, x2, bs=["gp","cr"], k=[8,5], '
+        'm=[[2,.7,1.5],None], sp=[.6,.8])',
+        'y ~ ti(x0, x1, x2, d=[2,1], bs=["gp","cr"], k=[10,5], '
+        'm=[[-5,.9],None], sp=[.6,.8])',
+    ],
+    ids=["te_univariate", "ti_multivariate_stationary"],
+)
+def test_gp_tensor_margins_fixed_sp_match_mgcv(formula):
+    data = _gp_data(seed=307, n=150)
+    actual = _fit_nampy_snapshot(data, formula, "gaussian", "fixed")
+    expected = _run_mgcv_snapshot(data, formula, "gaussian", "REML")
+    _assert_snapshot_fit(actual, expected, atol=1e-6)
+
+
+@pytest.mark.parametrize(
+    "formula",
+    [
+        'y ~ s(f, x0, x1, bs="fs", xt="gp", k=10, m=[-1,.6], '
+        'sp=[.7,.9])',
+        'y ~ s(f, f1, x0, x1, bs="sz", xt="gp", k=10, '
+        'm=[1,.6], id="shared_gp_sz", sp=.7)',
+    ],
+    ids=["fs", "sz"],
+)
+def test_gp_multivariate_factor_smooth_base_fixed_sp_matches_mgcv(formula):
+    data = _gp_data(seed=308, n=150)
+    actual = _fit_nampy_snapshot(data, formula, "gaussian", "fixed")
+    expected = _run_mgcv_snapshot(data, formula, "gaussian", "REML")
+    _assert_snapshot_fit(actual, expected, atol=1e-6)
+
+
+def test_gp_persistence_preserves_chunked_extrapolation(tmp_path):
+    data = _gp_data(seed=309, n=130)
+    formula = (
+        'y ~ s(x0, x1, bs="gp", k=10, m=[-1,.6], '
+        'xt={"max.knots":14,"seed":7}, sp=.7)'
+    )
+    model = _fit_nampy_model(data, formula, "gaussian", "fixed")
+    newdata = pd.DataFrame(
+        {
+            "x0": np.linspace(-4.0, 4.0, 31),
+            "x1": np.linspace(-3.0, 3.5, 31),
+        }
+    )
+    expected = model.predict(newdata, type="link", block_size=1)
+    path = tmp_path / "gp.pkl"
+    model.save_model(path)
+    restored = GAM.load_model(path)
+    np.testing.assert_allclose(
+        restored.predict(newdata, type="link", block_size=1), expected
+    )
+    _assert_snapshot_fit(
+        model.parity_snapshot(X=data, include_covariances=True),
+        _run_mgcv_snapshot(data, formula, "gaussian", "REML"),
+    )
+
+
+def test_gp_array_api_with_explicit_k_builds_one_term_per_feature():
+    data = _gp_data(seed=312, n=100)
+    features = data[["x0", "x1"]]
+    model = GAM(
+        family="gaussian",
+        basis="gp",
+        k=10,
+        optimize_smoothing=False,
+        smoothing_params=[0.5, 0.8],
+    ).fit(X=features, y=data["y"].to_numpy(dtype=np.float64))
+    assert len(model.smoothing_params) == 2
+    assert model.predict(features.iloc[:7]).shape == (7,)
+
+
+def test_gp_definition_warnings_and_validation():
+    x = np.linspace(-1.0, 1.0, 30)
+    X = np.column_stack([x, np.sin(x)])
+    rounded = build_gaussian_process_setup(X, k=8, m=[1.6, 0.7, 1.5])
+    np.testing.assert_allclose(rounded.definition, [2.0, 0.7, 1.5])
+
+    with pytest.warns(UserWarning, match="basis dimension reset"):
+        assert build_gaussian_process_setup(X, k=2, m=[-3, 0.7]).bs_dim == 4
+    with pytest.warns(UserWarning, match="more knots than data in an ms term"):
+        ignored = build_gaussian_process_setup(
+            X, k=8, m=[3, 0.7], knots=[np.arange(31), np.arange(31)]
+        )
+    assert not ignored.used_supplied_knots
+
+    with pytest.raises(ValueError, match="incorrect arguments"):
+        build_gaussian_process_setup(X, k=8, m=[0, 0.7])
+    with pytest.raises(ValueError, match="incorrect arguments"):
+        build_gaussian_process_setup(X, k=8, m=[1, 0.7, 2.1])
+    with pytest.raises(ValueError, match="same length"):
+        build_gaussian_process_setup(
+            X,
+            k=8,
+            m=[3, 0.7],
+            knots=[np.linspace(-1, 1, 10), np.linspace(-1, 1, 9)],
+        )
+    with pytest.raises(ValueError, match="fewer unique covariate combinations"):
+        build_gaussian_process_setup(np.repeat(X[:4], 5, axis=0), k=8)
+    with pytest.raises(ValueError, match="at least as many knot locations"):
+        build_gaussian_process_setup(
+            X,
+            k=10,
+            m=[3, 0.7],
+            knots=[np.arange(6), np.arange(6)],
+        )
+    with pytest.raises(ValueError, match="supply k explicitly"):
+        build_gaussian_process_setup(np.column_stack([X, X]), k=-1)
+
+
+def test_gp_linked_pc_different_features_and_derivative_are_explicitly_unsupported():
+    data = _gp_data(seed=310, n=80)
+    with pytest.raises(NotImplementedError, match="different feature sets"):
+        GAM(
+            formula=(
+                'y ~ s(x0, bs="gp", k=10, pc=[0], id="gp_pc")'
+                ' + s(x1, bs="gp", k=10, pc=[0], id="gp_pc")'
+            ),
+            optimize_smoothing=False,
+        ).fit(data=data)
+
+    model = GAM(
+        formula='y ~ s(x0, bs="gp", k=10, sp=.7)',
+        optimize_smoothing=False,
+    ).fit(data=data)
+    with pytest.raises(NotImplementedError, match="derivative provider"):
+        model.derivative(data, smooth_number=1)
+
+
+def test_gp_tensor_stationarity_requires_nested_m():
+    data = _gp_data(seed=311, n=90)
+    with pytest.raises(ValueError, match="incorrect arguments to GP smoother"):
+        GAM(
+            formula='y ~ te(x0, x1, bs=["gp","gp"], k=[8,8], m=[-3,-3])',
+            optimize_smoothing=False,
+        ).fit(data=data)
+
+    nested = GAM(
+        formula='y ~ te(x0, x1, bs=["gp","gp"], k=[8,8], m=[[-3],[-3]])',
+        optimize_smoothing=False,
+    ).fit(data=data)
+    assert nested.predict(data).shape == (len(data),)
diff --git a/tests/reference_fixtures/mgcv/03afb257b606dc2aabe12601f681cceaf7f8fa09e91550cfb53fb98f3c4c408c.json.gz b/tests/reference_fixtures/mgcv/03afb257b606dc2aabe12601f681cceaf7f8fa09e91550cfb53fb98f3c4c408c.json.gz
new file mode 100644
index 00000000..3e7a07ff
Binary files /dev/null and b/tests/reference_fixtures/mgcv/03afb257b606dc2aabe12601f681cceaf7f8fa09e91550cfb53fb98f3c4c408c.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/0da72acd5c05004ca793a5e0d9b653afd1180deb4b49281bf189b3ff384b70d1.json.gz b/tests/reference_fixtures/mgcv/0da72acd5c05004ca793a5e0d9b653afd1180deb4b49281bf189b3ff384b70d1.json.gz
new file mode 100644
index 00000000..fa17d3a5
Binary files /dev/null and b/tests/reference_fixtures/mgcv/0da72acd5c05004ca793a5e0d9b653afd1180deb4b49281bf189b3ff384b70d1.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/22b839f0dfd05d2f290ca76a6632a6e037575b4bbb89216a9201fc04d5159aa1.json.gz b/tests/reference_fixtures/mgcv/22b839f0dfd05d2f290ca76a6632a6e037575b4bbb89216a9201fc04d5159aa1.json.gz
new file mode 100644
index 00000000..105e3b67
Binary files /dev/null and b/tests/reference_fixtures/mgcv/22b839f0dfd05d2f290ca76a6632a6e037575b4bbb89216a9201fc04d5159aa1.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/3f00065b333c6ad95f3ae3cc48d5ef9acde68612eec29933ffb5f6fde03ef0f7.json.gz b/tests/reference_fixtures/mgcv/3f00065b333c6ad95f3ae3cc48d5ef9acde68612eec29933ffb5f6fde03ef0f7.json.gz
new file mode 100644
index 00000000..4aca27cc
Binary files /dev/null and b/tests/reference_fixtures/mgcv/3f00065b333c6ad95f3ae3cc48d5ef9acde68612eec29933ffb5f6fde03ef0f7.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/40c55686c1fe6ca835ff7a764c6aa5c3c1abadde3f3e6eb8e643dc57bf27e07c.json.gz b/tests/reference_fixtures/mgcv/40c55686c1fe6ca835ff7a764c6aa5c3c1abadde3f3e6eb8e643dc57bf27e07c.json.gz
new file mode 100644
index 00000000..951e5feb
Binary files /dev/null and b/tests/reference_fixtures/mgcv/40c55686c1fe6ca835ff7a764c6aa5c3c1abadde3f3e6eb8e643dc57bf27e07c.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/5b146985163b8faf499298269d13581ccae2e4f7688778fcc74b581985c15309.json.gz b/tests/reference_fixtures/mgcv/5b146985163b8faf499298269d13581ccae2e4f7688778fcc74b581985c15309.json.gz
new file mode 100644
index 00000000..f370e500
Binary files /dev/null and b/tests/reference_fixtures/mgcv/5b146985163b8faf499298269d13581ccae2e4f7688778fcc74b581985c15309.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/64760b0e53015494fe31e9277d071d2a9bbd32c2eecf6d235277f37157faab08.json.gz b/tests/reference_fixtures/mgcv/64760b0e53015494fe31e9277d071d2a9bbd32c2eecf6d235277f37157faab08.json.gz
new file mode 100644
index 00000000..1a8b8cb8
Binary files /dev/null and b/tests/reference_fixtures/mgcv/64760b0e53015494fe31e9277d071d2a9bbd32c2eecf6d235277f37157faab08.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/6cca1ccfcd8e3bf978422430903f21631f46b34f21e1cef8bc83ec692e2917a9.json.gz b/tests/reference_fixtures/mgcv/6cca1ccfcd8e3bf978422430903f21631f46b34f21e1cef8bc83ec692e2917a9.json.gz
new file mode 100644
index 00000000..2c6f9766
Binary files /dev/null and b/tests/reference_fixtures/mgcv/6cca1ccfcd8e3bf978422430903f21631f46b34f21e1cef8bc83ec692e2917a9.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/7251f8748a8d91949dbf0e4f89ff2c51ef1667aed45fa93415ce2a267e2aa09f.json.gz b/tests/reference_fixtures/mgcv/7251f8748a8d91949dbf0e4f89ff2c51ef1667aed45fa93415ce2a267e2aa09f.json.gz
new file mode 100644
index 00000000..a61f8b7a
Binary files /dev/null and b/tests/reference_fixtures/mgcv/7251f8748a8d91949dbf0e4f89ff2c51ef1667aed45fa93415ce2a267e2aa09f.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/7e22726ace7ed90c3c518ddb7a1ced68ec59ee50c490ad621cc2aa3946e953c5.json.gz b/tests/reference_fixtures/mgcv/7e22726ace7ed90c3c518ddb7a1ced68ec59ee50c490ad621cc2aa3946e953c5.json.gz
new file mode 100644
index 00000000..4b74e7dd
Binary files /dev/null and b/tests/reference_fixtures/mgcv/7e22726ace7ed90c3c518ddb7a1ced68ec59ee50c490ad621cc2aa3946e953c5.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/84e69fa5486b8ae184e01507742c7c67c0ad342c1ab63ea9fb0cd988ddb230b9.json.gz b/tests/reference_fixtures/mgcv/84e69fa5486b8ae184e01507742c7c67c0ad342c1ab63ea9fb0cd988ddb230b9.json.gz
new file mode 100644
index 00000000..a9b1e8a4
Binary files /dev/null and b/tests/reference_fixtures/mgcv/84e69fa5486b8ae184e01507742c7c67c0ad342c1ab63ea9fb0cd988ddb230b9.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/8dc5f1aa43317d040b3c4b888d715d653ae2a8dad2369a67deb2a9fbaad39feb.json.gz b/tests/reference_fixtures/mgcv/8dc5f1aa43317d040b3c4b888d715d653ae2a8dad2369a67deb2a9fbaad39feb.json.gz
new file mode 100644
index 00000000..4e68a63b
Binary files /dev/null and b/tests/reference_fixtures/mgcv/8dc5f1aa43317d040b3c4b888d715d653ae2a8dad2369a67deb2a9fbaad39feb.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/9109ade814a36a953c0c3393b50ace0309af3e847fafc31a6682eda8fdb24de5.json.gz b/tests/reference_fixtures/mgcv/9109ade814a36a953c0c3393b50ace0309af3e847fafc31a6682eda8fdb24de5.json.gz
new file mode 100644
index 00000000..ab649441
Binary files /dev/null and b/tests/reference_fixtures/mgcv/9109ade814a36a953c0c3393b50ace0309af3e847fafc31a6682eda8fdb24de5.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/9a3275d48edcf5bdaec3733fd0c0d8514a0f6a5a93dd534495c2572180b56fd7.json.gz b/tests/reference_fixtures/mgcv/9a3275d48edcf5bdaec3733fd0c0d8514a0f6a5a93dd534495c2572180b56fd7.json.gz
new file mode 100644
index 00000000..38ff8896
Binary files /dev/null and b/tests/reference_fixtures/mgcv/9a3275d48edcf5bdaec3733fd0c0d8514a0f6a5a93dd534495c2572180b56fd7.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/a979cc59c412790592955bcddd64fc385e9ba98e564dbb5660f83d62d005be9c.json.gz b/tests/reference_fixtures/mgcv/a979cc59c412790592955bcddd64fc385e9ba98e564dbb5660f83d62d005be9c.json.gz
new file mode 100644
index 00000000..7ff2628c
Binary files /dev/null and b/tests/reference_fixtures/mgcv/a979cc59c412790592955bcddd64fc385e9ba98e564dbb5660f83d62d005be9c.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/be6bcf4b6042c290c1220b518469149f51b0972393ea81e345d1f515c1739837.json.gz b/tests/reference_fixtures/mgcv/be6bcf4b6042c290c1220b518469149f51b0972393ea81e345d1f515c1739837.json.gz
new file mode 100644
index 00000000..3ce8c533
Binary files /dev/null and b/tests/reference_fixtures/mgcv/be6bcf4b6042c290c1220b518469149f51b0972393ea81e345d1f515c1739837.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/c2e3e1399ee9dd12141beb621b9b65ab530cffdee89a4ad762a8dddc7ab7d83e.json.gz b/tests/reference_fixtures/mgcv/c2e3e1399ee9dd12141beb621b9b65ab530cffdee89a4ad762a8dddc7ab7d83e.json.gz
new file mode 100644
index 00000000..ba55d478
Binary files /dev/null and b/tests/reference_fixtures/mgcv/c2e3e1399ee9dd12141beb621b9b65ab530cffdee89a4ad762a8dddc7ab7d83e.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/c2e91b7dca1828a2deb9fb677895c561a01d28865f894a37da620477085d2ace.json.gz b/tests/reference_fixtures/mgcv/c2e91b7dca1828a2deb9fb677895c561a01d28865f894a37da620477085d2ace.json.gz
new file mode 100644
index 00000000..c0dbbcc0
Binary files /dev/null and b/tests/reference_fixtures/mgcv/c2e91b7dca1828a2deb9fb677895c561a01d28865f894a37da620477085d2ace.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/c48025d8992bba6d93601d5dc4b31d6bd94104aa51fd38a63408dc8ca7a91eb8.json.gz b/tests/reference_fixtures/mgcv/c48025d8992bba6d93601d5dc4b31d6bd94104aa51fd38a63408dc8ca7a91eb8.json.gz
new file mode 100644
index 00000000..5feeec47
Binary files /dev/null and b/tests/reference_fixtures/mgcv/c48025d8992bba6d93601d5dc4b31d6bd94104aa51fd38a63408dc8ca7a91eb8.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/d11e343500c1464cc7fdc18571c348bcfd8bb8317e7d9e6f201bab11dd71729e.json.gz b/tests/reference_fixtures/mgcv/d11e343500c1464cc7fdc18571c348bcfd8bb8317e7d9e6f201bab11dd71729e.json.gz
new file mode 100644
index 00000000..40cc0090
Binary files /dev/null and b/tests/reference_fixtures/mgcv/d11e343500c1464cc7fdc18571c348bcfd8bb8317e7d9e6f201bab11dd71729e.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/da77452e3c5c4dcc0fd6cffba9777fa8916fbf91a59e349bb25fcb8083a1a07a.json.gz b/tests/reference_fixtures/mgcv/da77452e3c5c4dcc0fd6cffba9777fa8916fbf91a59e349bb25fcb8083a1a07a.json.gz
new file mode 100644
index 00000000..da28bb7e
Binary files /dev/null and b/tests/reference_fixtures/mgcv/da77452e3c5c4dcc0fd6cffba9777fa8916fbf91a59e349bb25fcb8083a1a07a.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/dc2eaedbdd661f04f69e780147fabb72617ffa801e819f230837d5b5416ba2bf.json.gz b/tests/reference_fixtures/mgcv/dc2eaedbdd661f04f69e780147fabb72617ffa801e819f230837d5b5416ba2bf.json.gz
new file mode 100644
index 00000000..fba60403
Binary files /dev/null and b/tests/reference_fixtures/mgcv/dc2eaedbdd661f04f69e780147fabb72617ffa801e819f230837d5b5416ba2bf.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/e4527392ccbc1d1fda5b3c5768e5892b28d138c5c2c2734474bd811d363ab0bd.json.gz b/tests/reference_fixtures/mgcv/e4527392ccbc1d1fda5b3c5768e5892b28d138c5c2c2734474bd811d363ab0bd.json.gz
new file mode 100644
index 00000000..d16267dc
Binary files /dev/null and b/tests/reference_fixtures/mgcv/e4527392ccbc1d1fda5b3c5768e5892b28d138c5c2c2734474bd811d363ab0bd.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/e8d1b70c1495f8264328110cc06975e8ad9679894cb1e01a85855bfc5e5d593a.json.gz b/tests/reference_fixtures/mgcv/e8d1b70c1495f8264328110cc06975e8ad9679894cb1e01a85855bfc5e5d593a.json.gz
new file mode 100644
index 00000000..2af943e8
Binary files /dev/null and b/tests/reference_fixtures/mgcv/e8d1b70c1495f8264328110cc06975e8ad9679894cb1e01a85855bfc5e5d593a.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/fbb60a2328883171e218e4fc674beb6adb5a71ea557f5e9e71db63dbd39ea724.json.gz b/tests/reference_fixtures/mgcv/fbb60a2328883171e218e4fc674beb6adb5a71ea557f5e9e71db63dbd39ea724.json.gz
new file mode 100644
index 00000000..d6a594ff
Binary files /dev/null and b/tests/reference_fixtures/mgcv/fbb60a2328883171e218e4fc674beb6adb5a71ea557f5e9e71db63dbd39ea724.json.gz differ
diff --git a/tests/smooths/test_mgcv_raw_constructor_parity.py b/tests/smooths/test_mgcv_raw_constructor_parity.py
index 72b3af25..f16ced60 100644
--- a/tests/smooths/test_mgcv_raw_constructor_parity.py
+++ b/tests/smooths/test_mgcv_raw_constructor_parity.py
@@ -29,6 +29,7 @@
 from nampy.gam.smooths.univariate.bs import DerivativeBSplineTerm1D
 from nampy.gam.smooths.univariate.cr import CubicSplineTerm
 from nampy.gam.smooths.univariate.ds import DuchonSplineTerm
+from nampy.gam.smooths.univariate.gp import GaussianProcessTerm
 from nampy.gam.smooths.univariate.ps import PSplineTerm1D
 from nampy.gam.smooths.univariate.tp import ThinPlateSplineTerm
 from nampy.gam.specs.build import build_formula_model
@@ -639,6 +640,86 @@ def _build_duchon_case_matrix():
     ]
 
 
+def _build_gp_case_matrix():
+    return [
+        _case(
+            "gp_1d_default_k",
+            _factory(_make_univariate_data, seed=160, n=90),
+            'y ~ s(x, bs="gp")',
+            atol=2e-8,
+        ),
+        _case(
+            "gp_2d_default_k",
+            _factory(_make_gaussian_data, seed=161, n=100),
+            'y ~ s(x0, x1, bs="gp")',
+            atol=2e-8,
+        ),
+        _case(
+            "gp_2d_spherical",
+            _factory(_make_gaussian_data, seed=162, n=90),
+            'y ~ s(x0, x1, bs="gp", k=10, m=[1, .7])',
+            atol=2e-8,
+        ),
+        _case(
+            "gp_2d_stationary_spherical",
+            _factory(_make_gaussian_data, seed=163, n=90),
+            'y ~ s(x0, x1, bs="gp", k=10, m=[-1, .7])',
+            atol=2e-8,
+        ),
+        _case(
+            "gp_2d_power_exponential",
+            _factory(_make_gaussian_data, seed=164, n=90),
+            'y ~ s(x0, x1, bs="gp", k=10, m=[2, .8, 1.5])',
+            atol=2e-8,
+        ),
+        _case(
+            "gp_2d_matern_35",
+            _factory(_make_gaussian_data, seed=165, n=90),
+            'y ~ s(x0, x1, bs="gp", k=10, m=[5, .9])',
+            atol=2e-8,
+        ),
+        _case(
+            "gp_2d_supplied_truncated",
+            _factory(_make_gaussian_data, seed=166, n=90),
+            'y ~ s(x0, x1, bs="gp", k=10, m=[4, .8])',
+            atol=2e-8,
+            knots_factory=_observed_row_knots(["x0", "x1"], 14),
+        ),
+        _case(
+            "gp_2d_supplied_pure_knot",
+            _factory(_make_gaussian_data, seed=167, n=90),
+            'y ~ s(x0, x1, bs="gp", k=10, m=[3, .8])',
+            atol=2e-8,
+            knots_factory=_observed_row_knots(["x0", "x1"], 7),
+        ),
+        _case(
+            "gp_2d_stationary_pure_knot",
+            _factory(_make_gaussian_data, seed=168, n=90),
+            'y ~ s(x0, x1, bs="gp", k=10, m=[-3, .8])',
+            atol=2e-8,
+            knots_factory=_observed_row_knots(["x0", "x1"], 9),
+        ),
+        _case(
+            "gp_2d_max_knots_xt",
+            _factory(_make_gaussian_data, seed=169, n=70),
+            'y ~ s(x0, x1, bs="gp", k=10, xt={"max.knots": 14, "seed": 7})',
+            atol=2e-8,
+        ),
+        _case(
+            "gp_3d_basic",
+            _factory(_make_gaussian_data_3col, seed=170, n=100),
+            'y ~ s(x0, x1, x2, bs="gp", k=15)',
+            atol=5e-8,
+        ),
+        _case(
+            "gp_3d_default_k",
+            _factory(_make_gaussian_data_3col, seed=173, n=120),
+            'y ~ s(x0, x1, x2, bs="gp")',
+            atol=2e-6,
+        ),
+    ]
+
+
 def _build_re_case_matrix():
     penalty_multi = {
         "S": [
@@ -731,6 +812,18 @@ def _build_factor_smooth_case_matrix():
             'y ~ s(f1, f2, x, bs="sz", k=6, id="shared")',
             atol=2e-8,
         ),
+        _case(
+            "fs_base_gp_multivariate",
+            _factory(_make_factorized_gaussian_data, seed=171, n=96),
+            'y ~ s(f, x0, x1, bs="fs", xt="gp", k=10, m=[1, .6])',
+            atol=2e-8,
+        ),
+        _case(
+            "sz_base_gp_multivariate",
+            _factory(_make_sz_metric2d_data, seed=172, n=90),
+            'y ~ s(f1, f2, x0, x1, bs="sz", xt="gp", k=10, m=[1, .6])',
+            atol=2e-8,
+        ),
     ]
     fs_xt_cases = [
         ("cr", "cr"),
@@ -930,6 +1023,7 @@ def _build_tensor_case_matrix():
     *_build_bs_case_matrix(),
     *_build_tprs_case_matrix(),
     *_build_duchon_case_matrix(),
+    *_build_gp_case_matrix(),
     *_build_re_case_matrix(),
     *_build_factor_smooth_case_matrix(),
     *_build_tensor_case_matrix(),
@@ -1147,6 +1241,27 @@ def _serialize_duchon_raw(term):
     )
 
 
+def _serialize_gp_raw(term):
+    setup = term._setup
+    basis = np.asarray(setup.basis_train, dtype=np.float64)
+    return _common_raw_state(
+        "gp.smooth",
+        basis,
+        [np.asarray(setup.penalty, dtype=np.float64)],
+        rank=int(setup.rank),
+        null_space_dim=int(setup.null_space_dim),
+        extra={
+            "knt": np.asarray(setup.knots, dtype=np.float64),
+            "UZ": np.asarray(setup.UZ, dtype=np.float64),
+            "shift": np.asarray(setup.shift, dtype=np.float64),
+            "gp_defn": np.asarray(setup.definition, dtype=np.float64),
+            "used_supplied_knots": bool(setup.used_supplied_knots),
+            "used_subsampling": bool(setup.used_subsampling),
+            "pure_knot": bool(setup.knots.shape[0] == setup.rank),
+        },
+    )
+
+
 def _serialize_re_raw(term):
     B = np.asarray(term._basis_train, dtype=np.float64)
     q = int(B.shape[1])
@@ -1374,6 +1489,8 @@ def _serialize_term_raw(term, X):
         return _serialize_tprs_raw(term)
     if isinstance(term, DuchonSplineTerm):
         return _serialize_duchon_raw(term)
+    if isinstance(term, GaussianProcessTerm):
+        return _serialize_gp_raw(term)
     if isinstance(term, RandomEffectTerm):
         return _serialize_re_raw(term)
     if isinstance(term, FSmoothInteractionTerm):