5. Latent networks

Everything so far was one GP behind one variable. The package’s actual model is a network of latent nodes, a directed graph the data flows through, and the reason is not architectural ambition. The network is where geological hypotheses go. “The grade rides on a regional trend plus a local structure”, “these seven elements share three underlying processes”, “this orebody is folded, so distances should be measured in an unfolded space”: each of these is a network shape, trained end to end by the same ELBO.

5.1 The nodes

Sources. BasicInput turns coordinates into the network’s first latent variables and holds the inducing points. It is deterministic by design, because coordinates are facts, and its transform is the door where anisotropy, projections or a fault step in. That is chapter 2’s ellipsoid, made composable. GaussianInput is the same node for inputs that are not facts: a mean and a variance per coordinate, which is what a GaussianData container holds (chapter 10). It is built for a high-dimensional input with missing entries, each given the mean and variance it could have, and it serves an uncertain location the same way. The variance rides through the transform, and a BasicGP above takes the moments of the mixture over it: the expected kernel (5.1, depth) for the mean, and for the variance the spread of the posterior’s mean over the input as well, in closed form for the Gaussian kernel and over 64 points of the input for the others. So such a row is used for what it says rather than dropped or imputed. Before 0.9.0 BasicGP read an uncertain input through an inflated kernel, measured to understate the spread several times over, and UncertainInputGP, now deprecated, integrated over it by quadrature. For a location error a noise term already absorbs, telling the model nothing was measured the better choice. The variance of the mean reads the kernel at 64 points of every uncertain input against every inducing point, so it costs time and memory in proportion; dividing the inducing points among experts does not shorten it, but training them one at a time (train_by_expert) holds a fraction of it in memory.

Workers. BasicGP(parent, size=k) is the GP node of chapter 3. The size argument gives it \(k\) latent columns that share one kernel and one set of inducing points, which is cheap: one covariance structure serves all of them.

It is worth being precise about what those \(k\) columns are, because the shorthand “one node, several variables” invites the wrong reading. The columns are independent fields. They share a kernel, so they have the same range, the same shape at the origin and the same anisotropy, but knowing one of them tells you nothing about the others. Correlation is not something a BasicGP produces, and it has to be introduced deliberately, in one of two places:

In the network, by passing the node’s output through a Linear node, whose trainable matrix mixes the independent columns into correlated outputs. This is the package’s reading of the linear model of coregionalization, with the loadings learned rather than declared.
On the data side, through a mixing warping. PCA, RobustPCA and Rotation decorrelate the measured variables before they reach the field, so independent latent columns come back as correlated variables after the back-transform. Chapter 16 takes this route.

Combiners, each a hypothesis:

Add is superposition. A long-range node plus a short-range node is the trainable version of a nested variogram, with the split between structures estimated rather than declared.
Multiply is interaction, one field modulating another’s amplitude.
Linear maps its parent’s columns onto however many outputs are wanted, through a trainable matrix. Close to the root it rotates coordinates, in the middle it is a bottleneck, and at the end it is where correlation between outputs comes from.
LinearCombination mixes several parents of the same size with positive weights that sum to one, keeping the output’s variance under control.
Concatenate stacks parents side by side into a joint input for a deeper node.
Bias is a trainable constant, for when the mean is not the data’s.

Depth. A BasicGP whose parent is another BasicGP receives a distribution rather than a point. The parent’s uncertainty rides along, which is what makes the stack a deep GP instead of two models glued together. The outer node averages its kernel over that uncertainty, the expected kernel, and what decides the average is how uncertain the inner layer is about the difference between two locations: two places whose warped positions move together stay correlated, two whose positions are uncertain independently lose correlation. So each node hands the next its covariance between locations along with its variances. The practical reading is simpler: the inner layers warp space. A stationary kernel in the warped space is non-stationary in the real one, so folded veins and curved orebodies stop being kernel problems and become network problems, and the network is trainable.

How that uncertainty is carried is not a detail. Before 0.9.0 each location’s variance travelled alone and widened the outer node’s range by half of it, which barely registered: on a folded synthetic section a two-layer network’s predictive variance came out at a few thousandths of what its own posterior implied, training bought a tiny noise variance with that confidence, and new drillholes scored 25 nats each where the expected kernel scores 9 (GPOptions(propagation=...), whose "marginal" keeps the old rule for models saved under it). The outer node of a deep network takes the Gaussian, exponential, Matérn or rational quadratic kernel: the expected kernel needs a kernel that is a mixture of Gaussians.

Depth comes with one habit that is close to mandatory. Concatenate the inner node with the original coordinates before feeding the next layer. An inner GP is free to map two distant regions onto the same place, and if the outer layer sees only the inner node’s output it has no way to tell them apart, so the space collapses and points that are far away become artificially correlated. Keeping the coordinates in the joint input costs a couple of columns and removes the failure mode, and it is the recommended construction of every deep network in this manual.

5.2 Choosing a shape

Start with BasicInput → BasicGP and earn every addition:

variance left at short range that the kernel cannot take: Add a short-range node;
several variables that should share structure: one BasicGP(size=k) followed by a Linear, or a decorrelating warping on the data side;
mapped geometry the ellipsoid cannot express: depth, an inner GP or an SDE node warping space, which is the subject of the 2025 paper;
more data than one inducing set carries: experts (chapter 3).

The network mirrors belief, and it also documents it. to_dot() renders the graph, and the model’s repr prints every node with its parameters, worth reading before training the way one reads a variogram model before kriging with it.

5.3 A deep model, drawn

Two layers on Walker Lake: a two-column inner GP warping space, and a one-column outer GP reading the warped coordinates alongside the real ones. The inducing points are chapter 3’s grid of experts, and the warping is the positive chain of chapter 4, so the predictions cannot come back negative.

import geoml
import numpy as np

geoml.set_seed(1234)
walker, walker_grid = geoml.datasets.walker()

experts = geoml.data.inducing.grid_experts(walker_grid, 10.0, block=8)

root = geoml.latent.BasicInput(
    experts,
    transform=geoml.transform.Isotropic(50))

# the inner layer moves the ground: two smooth fields that the outer layer
# will read as if they were coordinates
inner = geoml.latent.BasicGP(
    root,
    size=2,
    kernel=geoml.kernels.Gaussian())

# the real coordinates travel alongside the warped ones, so that two
# distant places cannot end up at the same address
deep_input = geoml.latent.Concatenate(root, inner)

# a GP reading an uncertain input takes a kernel that is a mixture of
# Gaussians: the Matern32 for the spherical's roughness at the origin
outer = geoml.latent.BasicGP(
    deep_input,
    size=1,
    kernel=geoml.kernels.Matern32())

warping = geoml.warping.ChainedWarping(
    geoml.warping.BoxCox(1, shift=1.0),
    geoml.warping.ZScore(1))

model = geoml.models.VGPNetwork(
    walker, "V",
    geoml.likelihood.Gaussian(warping),
    outer,
    options=geoml.models.GPOptions(verbose=False))

model.train_full(max_iter=100)

print([node.name for node in [outer] + outer.get_unique_parents()])
print(model.to_dot().splitlines()[0])

The node list is the graph the ELBO trains, and the DOT text renders with any Graphviz viewer. Prediction confirms what the warping promised:

subset = geoml.data.PointData.from_array(
    np.asarray(walker_grid.coordinates)[::20])
model.predict(subset, n_sim=20)

print("lowest prediction:",
      round(float(np.min(subset.values("V/prediction"))), 1))

The value is positive, and it is positive by construction rather than by luck. Chapter 2’s model could return a negative grade wherever it was unsure. This one cannot, at any location and in any realization, because the chain’s inverse power returns a zero grade for anything the latent field puts below where a zero lands.

Depth costs iterations. A two-layer model settles more slowly than a flat one, and it repays the wait only where the geometry is genuinely curved, which Walker Lake’s is only mildly. The folded vein of chapter 17 is the honest showcase.

In the code. Every node class lives in geoml.latent, and viz/graphviz.py renders any model or network through to_dot(). The expert propagation of chapter 3 keys off this same graph, so what to_dot draws is exactly what runs.

References

Damianou, A., & Lawrence, N. D. (2013). Deep Gaussian processes. Proceedings of the 16th International Conference on Artificial Intelligence and Statistics (AISTATS), 207–215.

Lu, C.-K., & Shafto, P. (2021). Conditional deep Gaussian processes: empirical Bayes hyperdata learning. Entropy, 23(11), 1545. https://doi.org/10.3390/e23111545

Gonçalves, Í. G. et al. (2025). Uncertainty propagation in deep Gaussian process networks. Mathematical Geosciences. https://doi.org/10.1007/s11004-025-10187-4

Gonçalves, Í. G. et al. (2026). Scalable variational Gaussian process framework for implicit geological modelling and compositional grade interpolation. Artificial Intelligence in Geosciences. https://doi.org/10.1016/j.aiig.2026.100218

Paciorek, C. J., & Schervish, M. J. (2003). Nonstationary covariance functions for Gaussian process regression. Advances in Neural Information Processing Systems 16.

5. Latent networks

5.1 The nodes

5.2 Choosing a shape

5.3 A deep model, drawn

Further reading

References