gpmp.num module

The gpmp.num module is the numerical backend interface used by GPmp. It exposes the array operations needed by GPmp models, kernels, diagnostics, and samplers while keeping the backend choice in one place.

Backend-independent GPmp code usually imports this module as gpmp.num as gnp. Mean functions, covariance functions, and selection criteria should use gnp operations rather than direct backend-library calls unless the code is intentionally backend-specific.

This page documents the portable subset of gpmp.num used by GPmp code. The module may expose additional backend compatibility names, but code meant to run with both backends should rely on the functions described here.

Configuration and import order

The gpmp.config module owns the process-wide configuration used by GPmp. It records the requested backend, the floating-point dtype policy, the global random seed, the logger, and small internal caches. The gpmp.num module reads that configuration when it is loaded, imports the selected backend implementation, resolves the backend dtype object, and re-exports the numerical functions documented below.

This separation matters because not all settings have the same lifetime. Backend and dtype determine which functions and array classes are installed in gpmp.num; they must therefore be fixed before gpmp.num is loaded. The seed and log level remain ordinary runtime state and can be changed after import.

Backend and dtype are fixed when gpmp.num is imported. In ordinary scripts, choose the backend with environment variables before importing gpmp:

import os

os.environ["GPMP_BACKEND"] = "numpy"  # or "torch"
import gpmp.num as gnp

From a shell, the equivalent is:

GPMP_BACKEND=numpy python script.py

If GPMP_BACKEND is not set, GPmp uses the torch backend when PyTorch is available and otherwise uses NumPy. The dtype is fixed to float64.

The seed is runtime state. Use set_seed() after import to set the random generators used by GPmp:

import gpmp.num as gnp

gnp.set_seed(1234)

After import, the current configuration can be inspected with gpmp.config.get_config(). Runtime settings such as the seed and log level can still be changed. Backend and dtype changes require a new Python process.

Backend array contract

GPmp functions operate on arrays from the active numerical backend:

• With GPMP_BACKEND=numpy, backend arrays are NumPy ndarray objects.

• With GPMP_BACKEND=torch, backend arrays are PyTorch Tensor objects.

Inputs passed to GPmp should be arrays from the active backend or objects convertible by asarray(). Outputs use the active backend array type unless a function explicitly documents conversion to NumPy. Use to_np() at plotting, reporting, or external-library boundaries when a NumPy array is required.

GPmp uses float64 for floating-point computations. The resolved backend dtype is available through get_dtype(), and GPmp constructors default to it where applicable.

Shape conventions

The constructor convention is intentionally explicit.

• Deterministic constructors take the shape as one object: gnp.zeros((2, 3)), gnp.ones((n, d)), gnp.empty((p, p)).

• Random constructors take dimensions as positional arguments: gnp.rand(2, 3) and gnp.randn(n, d).

• One-dimensional shapes should still be written explicitly when using deterministic constructors: gnp.zeros((n,)).

This convention avoids exposing backend-specific constructor signatures such as separated dimensions for deterministic constructors.

Constructors and conversion

gpmp.num.asarray(x, dtype=None)

Convert x to an array of the active backend.

Floating-point inputs are converted to GPmp’s default floating dtype unless dtype is provided. Integer and boolean inputs keep an integer or boolean dtype when possible. Existing backend arrays are returned without a copy when no dtype conversion is needed.

gpmp.num.array(x, dtype=None)

Create an array of the active backend from x.

This is the constructor form of asarray(). Use it when a new array is expected. Use asarray() when accepting user input that may already be a backend array.

gpmp.num.empty(shape, dtype=None)

Return an uninitialized array with the given shape.

shape is an integer or a tuple of integers. Floating arrays use get_dtype() unless dtype is provided.

gpmp.num.zeros(shape, dtype=None)

Return an array of zeros with the given shape.

shape must be passed as one object, for example zeros((2, 3)).

gpmp.num.ones(shape, dtype=None)

Return an array of ones with the given shape.

shape must be passed as one object, for example ones((n, 1)).

gpmp.num.full(shape, fill_value, dtype=None)

Return an array with the given shape and constant fill_value.

gpmp.num.eye(n, m=None, k=0, dtype=None)

Return a two-dimensional identity-like array.

n is the number of rows. If m is None, the result is square. k selects the diagonal offset.

gpmp.num.triu(x, k=0)

Return the upper-triangular part of a two-dimensional backend array.

k selects the diagonal offset. Entries below that diagonal are set to zero.

gpmp.num.arange(*args, **kwargs)

Return evenly spaced values over an integer-like interval.

This follows the usual start, stop, step convention. Use it for index arrays and simple grids.

gpmp.num.linspace(start, stop, num=50, endpoint=True, retstep=False, dtype=None, axis=0)

Return num evenly spaced values between start and stop.

The result uses the active backend. Floating values use GPmp’s floating dtype unless dtype is provided.

gpmp.num.logspace(start, stop, num=50, endpoint=True, base=10.0, dtype=None, axis=0)

Return values spaced evenly on a logarithmic scale.

gpmp.num.meshgrid(*arrays, **kwargs)

Return coordinate arrays from one-dimensional coordinate vectors.

The return type follows the active backend. Use this for plotting grids and tensor-product designs.

gpmp.num.zeros_like(x)

Return an array with the same shape and backend as x.

gpmp.num.ones_like(x)

Return an array of ones with the same shape and backend as x.

gpmp.num.empty_like(x)

Return an uninitialized array with the same shape and backend as x.

gpmp.num.full_like(x, fill_value)

Return an array with the same shape and backend as x and constant fill_value.

gpmp.num.to_np(x)

Convert x to a NumPy array when needed.

With the NumPy backend this returns x. With the torch backend, tensors are detached, moved to CPU, and converted to NumPy arrays. Use this function at plotting or external-library boundaries, not inside backend-independent numerical code.

gpmp.num.to_scalar(x)

Convert a scalar backend object to a Python scalar.

Use this for logging, printing, and calls to external APIs that require a Python scalar. Do not use it inside differentiable computations.

gpmp.num.get_dtype()

Return the resolved floating-point dtype of the active backend.

The value is numpy.float64 with the NumPy backend and torch.float64 with the torch backend.

gpmp.num.asdouble(x)

Convert x to the backend double-precision floating dtype.

gpmp.num.asint(x)

Convert x to a backend integer array.

gpmp.num.isarray(x)

Return whether x is an array object of the active backend.

gpmp.num.copy(x)

Return a backend array copy of x.

gpmp.num.array_equal(x, y)

Return whether two arrays have the same shape and entries.

Elementwise mathematical functions

Elementwise functions apply component by component and return an array from the active backend. Inputs may be backend arrays or objects accepted by asarray(). Standard backend broadcasting rules apply.

gpmp.num.exp(x)

Return the elementwise exponential of x.

gpmp.num.log(x)

Return the elementwise natural logarithm of x.

gpmp.num.log10(x)

Return the elementwise base-10 logarithm of x.

gpmp.num.log1p(x)

Return the elementwise value of log(1 + x).

gpmp.num.gammaln(x)

Return the elementwise logarithm of the absolute gamma function.

This is used by Matérn covariance implementations and returns an array from the active backend.

gpmp.num.compute_gammaln(up_to_p)

Return a cached table of gammaln(k) values for integer k.

The result is a one-dimensional backend array. The cache is stored in gpmp.config and grows when a larger table is requested.

gpmp.num.sqrt(x)

Return the elementwise square root of x.

gpmp.num.abs(x)

Return the elementwise absolute value of x.

gpmp.num.sin(x)

Return the elementwise sine of x.

gpmp.num.cos(x)

Return the elementwise cosine of x.

gpmp.num.tan(x)

Return the elementwise tangent of x.

gpmp.num.tanh(x)

Return the elementwise hyperbolic tangent of x.

gpmp.num.floor(x)

Return the elementwise floor of x.

gpmp.num.ceil(x)

Return the elementwise ceiling of x.

gpmp.num.minimum(x1, x2)

Return the elementwise minimum of x1 and x2.

gpmp.num.maximum(x1, x2)

Return the elementwise maximum of x1 and x2.

gpmp.num.clip(x, min=None, max=None, out=None)

Clip values below min or above max.

gpmp.num.where(condition, x, y)

Select entries from x or y according to condition.

gpmp.num.isnan(x)

Return a boolean array indicating NaN entries.

gpmp.num.isinf(x)

Return a boolean array indicating infinite entries.

gpmp.num.isfinite(x)

Return a boolean array indicating finite entries.

gpmp.num.isclose(x, y, **kwargs)

Return elementwise closeness tests between x and y.

gpmp.num.allclose(x, y, **kwargs)

Return whether all entries of x and y are close within tolerances.

gpmp.num.logical_not(x)

Return the elementwise logical negation of x.

gpmp.num.logical_and(x, y)

Return the elementwise logical conjunction of x and y.

gpmp.num.logical_or(x, y)

Return the elementwise logical disjunction of x and y.

gpmp.num.nan_to_num(x, copy=True, nan=0.0, posinf=None, neginf=None)

Replace NaN and infinite values by finite values.

gpmp.num.inftobigf(x, bigf=None)

Replace infinite values by a large finite value.

This helper is mainly used in optimization and diagnostic code that must pass finite values to external routines.

Reductions and statistics

Reduction functions use NumPy-style argument names where possible. In particular, use axis rather than torch’s dim argument in backend-independent code.

gpmp.num.sum(x, axis=None, **kwargs)

Sum entries of x along axis.

If axis is None, all entries are reduced. Extra keyword arguments are forwarded to the active backend where supported.

gpmp.num.mean(x, axis=None, **kwargs)

Compute the arithmetic mean along axis.

gpmp.num.prod(x, axis=None, **kwargs)

Multiply entries of x along axis.

gpmp.num.cumsum(x, axis=None, **kwargs)

Return cumulative sums along axis.

gpmp.num.all(x, axis=None, **kwargs)

Return whether all entries along axis are true.

gpmp.num.any(x, axis=None, **kwargs)

Return whether any entry along axis is true.

gpmp.num.var(x, axis=None, ddof=1, keepdims=False)

Compute the variance along axis.

By default, GPmp computes the unbiased empirical variance with ddof=1, so the divisor is N - 1. More generally, ddof is the delta degrees of freedom and the divisor is N - ddof. Pass ddof=0 for the population convention. Use keepdims=True to keep reduced axes with length one.

gpmp.num.std(x, axis=None, ddof=1, keepdims=False)

Compute the standard deviation along axis.

The ddof and keepdims arguments have the same meaning as in var().

gpmp.num.cov(m, y=None, rowvar=True, bias=False, ddof=None, fweights=None, aweights=None, dtype=None)

Estimate a covariance matrix.

The API follows NumPy’s cov convention. If rowvar=True, each row is a variable and columns are observations. If rowvar=False, columns are variables. Frequency weights and analytic weights are backend-dependent and should be checked before use with the torch backend.

gpmp.num.percentile(x, q, axis=None, method='linear', keepdims=False, **kwargs)

Return percentiles of x.

q is expressed in percent. The portable interpolation method is "linear".

gpmp.num.norm(x, axis=None, ord=2)

Return a vector or matrix norm.

Use axis for vector norms along an axis. Matrix norm support depends on the active backend.

gpmp.num.min(x, axis, keepdims=False)

Return minimum values along axis.

Pass axis explicitly in backend-independent code.

gpmp.num.max(x, axis, keepdims=False)

Return maximum values along axis.

Pass axis explicitly in backend-independent code.

gpmp.num.argmin(x, axis=None)

Return indices of minimum values.

gpmp.num.argmax(x, axis=None)

Return indices of maximum values.

Array composition

These functions preserve the active backend and follow NumPy-style naming. They are useful when writing backend-independent covariance and mean functions.

gpmp.num.reshape(x, shape)

Return x with a new shape.

gpmp.num.expand_dims(x, axis)

Insert a length-one axis at position axis.

gpmp.num.concatenate(arrays, axis=0)

Join arrays along an existing axis.

gpmp.num.hstack(arrays)

Stack arrays horizontally.

gpmp.num.vstack(arrays)

Stack arrays vertically.

gpmp.num.stack(arrays, axis=0)

Join arrays along a new axis.

gpmp.num.tile(x, reps)

Construct an array by repeating x according to reps.

gpmp.num.split(x, indices_or_sections, axis=0)

Split an array along axis.

gpmp.num.transpose(x, dim0, dim1)

Swap two axes of x.

gpmp.num.diag(x, k=0)

Extract a diagonal from a two-dimensional array, or construct a diagonal array from a one-dimensional input.

Linear algebra

Linear-algebra functions expect finite arrays from the active backend. They may raise backend linear-algebra exceptions when a matrix is singular, not positive definite, or has incompatible shape. Higher-level optimization wrappers may catch some of these failures and convert them to infinite criterion values.

gpmp.num.safe_inf()

Return a backend-compatible positive-infinity sentinel.

This is mainly for GPmp internals that need a scalar failure value compatible with the active backend.

gpmp.num.safe_neginf()

Return a backend-compatible negative-infinity sentinel.

This is mainly for GPmp internals that need a scalar failure value compatible with the active backend.

gpmp.num.solve(A, B, **kwargs)

Solve the linear system A X = B.

A has shape (n, n). B has shape (n,) or (n, k). The result has the same backend as A.

gpmp.num.solve_triangular(A, B, trans=0, lower=False, unit_diagonal=False, **kwargs)

Solve a triangular linear system.

A is triangular with shape (n, n). B has shape (n,) or (n, k). Set lower=True when A is lower triangular. trans may request the non-transposed, transposed, or conjugate-transposed system.

gpmp.num.inv(A)

Return the inverse of a square matrix.

A has shape (n, n). Singular matrices raise a backend linear-algebra exception.

gpmp.num.cholesky(A)

Return a Cholesky factor of a positive-definite matrix.

A has shape (n, n). Non-positive-definite inputs raise a backend linear-algebra exception.

gpmp.num.cho_factor(A, lower=False, **kwargs)

Compute a Cholesky factorization of a positive-definite matrix.

Returns (C, lower). The factor orientation follows lower.

gpmp.num.cho_solve(c_and_lower, b, **kwargs)

Solve a positive-definite system from a Cholesky factorization returned by cho_factor().

gpmp.num.cholesky_solve(A, b)

Solve A x = b by Cholesky factorization.

Returns (x, L), where L is the Cholesky factor used internally.

gpmp.num.cholesky_inv(A)

Return the inverse of a positive-definite matrix using a Cholesky-based computation where available.

gpmp.num.logdet(A)

Return the logarithm of the determinant of a square matrix.

This is intended for positive-definite covariance matrices. Non-positive or singular determinants raise an exception or return backend-specific non-finite values.

gpmp.num.svd(A, full_matrices=True, hermitian=True)

Compute a singular value decomposition.

Return values follow the active backend. Use this for diagnostics rather than core GP solves when a Cholesky or triangular solve is available.

gpmp.num.qr(A, mode='reduced')

Compute a QR factorization of a two-dimensional array.

gpmp.num.cond(A)

Return a condition-number estimate for a matrix.

Distance functions

The distance functions implement GPmp’s anisotropic lengthscale convention. loginvrho stores -log(rho_j) for each input coordinate, so exp(loginvrho) is the inverse lengthscale vector.

gpmp.num.scaled_distance(loginvrho, x, y)

Compute the full matrix of scaled Euclidean distances.

x has shape (n, d), y has shape (m, d), and loginvrho has shape (d,) or is broadcastable to the coordinate dimension. The result has shape (n, m).

gpmp.num.scaled_squared_distance(loginvrho, x, y, zero_diagonal=False)

Compute the full matrix of squared scaled Euclidean distances.

If y is None, the function uses y = x. Set zero_diagonal=True to force a zero diagonal for square self-distance matrices.

gpmp.num.scaled_distance_elementwise(loginvrho, x, y)

Compute row-by-row scaled distances.

x and y must have the same shape (n, d) unless y is None or x is y. The result has shape (n,).

gpmp.num.scaled_squared_distance_elementwise(loginvrho, x, y)

Compute row-by-row squared scaled distances.

Matérn correlation primitive

gpmp.num.matern_correlation_from_squared_distance(q, nu, *, nodes=160, t_max=20.0, large_nu_method='auto')

Evaluate the normalized Matérn correlation from a nonnegative squared distance q and a positive regularity nu.

q and nu follow the active backend’s broadcasting rules. The result has their broadcast shape and uses the active backend. At q == 0, the value is exactly one.

large_nu_method accepts "auto", "quadrature", and "uniform_asymptotic" with both backends. The uniform expansion requires scalar nu >= 70. Under "auto", both backends use this expansion in that range. Below it, eligible CPU float64 Torch evaluations use either a recurrence or quadrature according to the calibrated dispatch profile, while NumPy uses SciPy. nodes and t_max control the base Torch quadrature rule and have no effect with NumPy. See Numerical evaluation of the Matérn covariance for the dispatch policy, scaling convention, and Torch derivatives.

Random numbers

Random functions use GPmp’s process-wide random generators. Call set_seed() to reset them. With the torch backend, set_seed() also resets the NumPy-compatible generator returned by random_generator(), which is used by SciPy/NumPy utilities such as design generation.

Random constructors accept positional dimensions or a shape tuple:

x = gnp.rand(10, 2)
z = gnp.randn(10)

gpmp.num.set_seed(seed)

Reset GPmp’s random generators to seed.

gpmp.num.random_generator()

Return the NumPy-compatible random generator controlled by set_seed().

This is intended for NumPy/SciPy utilities that need an explicit generator.

gpmp.num.rand(*shape)

Draw independent uniform random numbers on [0, 1).

Dimensions can be passed as positional integers or as one tuple. rand(2, 3) and rand((2, 3)) both return an array with shape (2, 3). rand() returns a scalar backend object.

gpmp.num.randn(*shape)

Draw independent standard normal random numbers.

Dimensions follow the same convention as rand().

gpmp.num.choice(a, size=None, replace=True, p=None)

Draw random samples from a one-dimensional population a.

p gives optional sampling probabilities. The return type follows the active backend.

gpmp.num.permutation(x)

Return a random permutation.

If x is an integer, the result is a permutation of 0, ..., x - 1. If x is an array, the result is a permutation along the first axis.

class gpmp.num.normal

Univariate normal distribution helper with backend-returning methods.

cdf(x, loc=0.0, scale=1.0)

Return the normal cumulative distribution function at x.

logcdf(x, loc=0.0, scale=1.0)

Return the logarithm of the normal cumulative distribution function.

pdf(x, loc=0.0, scale=1.0)

Return the normal probability density function at x.

logpdf(x, loc=0.0, scale=1.0)

Return the logarithm of the normal probability density function.

class gpmp.num.multivariate_normal

Multivariate normal distribution helper.

rvs(mean=0.0, cov=1.0, n=1)

Draw n samples. cov is either a scalar variance or a square covariance matrix.

logpdf(x, mean=0.0, cov=1.0)

Return log-density values. This method is available with the NumPy backend.

cdf(x, mean=0.0, cov=1.0)

Return cumulative probabilities. This method is available with the NumPy backend.

Differentiation

Differentiation support depends on the active backend.

• With the torch backend, gradients use automatic differentiation. The function being differentiated must return a scalar tensor connected to its input.

• With the NumPy backend, gradients use finite differences and are intended for low-dimensional scalar objectives.

gpmp.num.grad(f)

Return a function that evaluates the gradient of a scalar function f.

The returned function accepts one array argument and returns an array with the same shape. Under the NumPy backend this uses finite differences. Under the torch backend this uses automatic differentiation.

gpmp.num.value_and_grad(f, x)

Return (value, gradient) for the scalar function f at x.

value is scalar-like and gradient has the same shape as x. With the torch backend, non-finite scalar values return a zero gradient.

class gpmp.num.DifferentiableSelectionCriterion(criterion, x, z)

Wrap a scalar selection criterion criterion(param, x, z) for use by SciPy optimizers and GPmp parameter-selection routines.

evaluate(param)

Evaluate the criterion at param.

evaluate_no_grad(param)

Evaluate the criterion without preparing a gradient computation.

evaluate_pre_grad(param)

Evaluate the criterion and store the state needed by gradient().

gradient(param)

Return the gradient from the last evaluate_pre_grad() call.

With the torch backend, calling this method before evaluate_pre_grad() raises an error. With the NumPy backend, selection gradients are usually supplied by finite-difference wrappers.

class gpmp.num.BatchDifferentiableSelectionCriterion(criterion, loader, reduction='mean', batches_per_eval=0)

Wrap a scalar selection criterion evaluated over batches of observations.

loader yields (x_batch, z_batch) pairs. reduction is "mean" or "sum". If batches_per_eval is zero, each evaluation uses the full loader. If it is positive, each evaluation uses that many batches and cycles through the loader.

evaluate(param)

Evaluate the batch-reduced criterion.

evaluate_no_grad(param)

Evaluate the batch-reduced criterion without preparing a gradient.

evaluate_pre_grad(param)

Evaluate the criterion and prepare the gradient when supported by the backend.

gradient(param)

Return the gradient prepared by evaluate_pre_grad() when available.