BUG: HSGP set_boundary causes unhandled divide by zero and inf eigenvalues when feature range is zero

Is there an existing issue for this?

• I have searched the existing issues

What happened?

In pymc.gp.hsgp_approx, set_boundary(X, c) calculates the domain radius $S$ via:

S = (pt.max(X, axis=0) - pt.min(X, axis=0)) / 2.0
L = (c * S).eval()

If an input matrix X contains any feature dimension with zero variance (e.g. constant covariate column or single-point coordinate slice where pt.max(X) == pt.min(X)), S evaluates to 0.0, setting L[d] = 0.0.

Downstream in calc_eigenvalues(L, m):

return np.square((np.pi * S_arr) / (2 * L))

Division by zero occurs (RuntimeWarning: divide by zero encountered in divide), populating eigvals with infs and corrupting basis matrix phi without raising an explicit ValueError.

Reproduction Code

import numpy as np
import pymc as pm
from pymc.gp.hsgp_approx import set_boundary, calc_eigenvalues

# Input matrix where 1st feature dimension has zero spatial range
X = np.array([[1.0, 5.0], [1.0, 10.0], [1.0, 15.0]])
c = 1.5
m = [5, 5]

L = set_boundary(X, c)            # L[0] evaluates to 0.0
eigvals = calc_eigenvalues(L, m)   # Raises RuntimeWarning and populates infs

Expected Behavior

set_boundary / calc_eigenvalues should validate that L > 0 for all feature dimensions, or raise a descriptive ValueError(f"Dimension {d} in X has zero spatial range. Boundary length L must be positive.") instead of silently propagating inf and NaN values into the Laplacian eigenvector basis.

Record

• I agree to follow this project's Code of Conduct
• I want to work on this issue