diff --git a/src/ellphi/grad.py b/src/ellphi/grad.py
index 35ffcc9..741e57c 100644
--- a/src/ellphi/grad.py
+++ b/src/ellphi/grad.py
@@ -183,11 +183,35 @@ def coef_from_cov_grad(
     if cov_arr.ndim == 2:
         cov_arr = cov_arr[np.newaxis, :, :]
 
+    if centers.shape[0] != cov_arr.shape[0]:
+        raise ValueError("Mismatch between number of centres and covariance matrices")
+    if cov_arr.shape[-1] != cov_arr.shape[-2]:
+        raise ValueError("Covariance matrices must be square")
+
     n, d = centers.shape
+
+    if cov_arr.shape[-1] != d:
+        raise ValueError("Centre dimensionality and covariance size must agree")
     s2 = scale**2
 
-    # Forward pass
-    inv_cov = np.linalg.inv(cov_arr)  # (n, d, d)
+    # Forward pass – handle singular covariances the same way coef_from_cov
+    # does: return NaN coefficients so callers can detect degeneracy.
+    try:
+        inv_cov = np.linalg.inv(cov_arr)  # (n, d, d)
+    except np.linalg.LinAlgError:
+        m = (d + 1) * (d + 2) // 2
+        coefs = np.full((n, m), np.nan, dtype=float)
+
+        def vjp_nan(
+            grad_coefs: np.ndarray,
+        ) -> tuple[np.ndarray, np.ndarray]:
+            return (
+                np.full_like(centers, np.nan),
+                np.full_like(cov_arr, np.nan),
+            )
+
+        return coefs, vjp_nan
+
     A = inv_cov / s2  # (n, d, d)
     b = -np.einsum("nij,nj->ni", A, centers)  # (n, d)
     c = np.einsum("ni,nij,nj->n", centers, A, centers)  # (n,)
diff --git a/tests/test_coef_from_cov_grad.py b/tests/test_coef_from_cov_grad.py
index ff419c2..56bf8d5 100644
--- a/tests/test_coef_from_cov_grad.py
+++ b/tests/test_coef_from_cov_grad.py
@@ -3,6 +3,7 @@
 from __future__ import annotations
 
 import numpy as np
+import pytest
 
 from ellphi import coef_from_cov_grad
 from ellphi.geometry import coef_from_cov
@@ -185,3 +186,40 @@ def test_grad_shapes(rng):
     grad_X, grad_cov = vjp(g)
     assert grad_X.shape == (n, d)
     assert grad_cov.shape == (n, d, d)
+
+
+# ---------------------------------------------------------------------------
+# Input validation
+# ---------------------------------------------------------------------------
+
+
+def test_batch_size_mismatch():
+    """Mismatched centre/covariance batch sizes raise ValueError."""
+    with pytest.raises(ValueError, match="Mismatch"):
+        coef_from_cov_grad(np.zeros((3, 2)), np.eye(2)[np.newaxis])
+
+
+def test_non_square_cov():
+    """Non-square covariance raises ValueError."""
+    with pytest.raises(ValueError, match="square"):
+        coef_from_cov_grad(np.zeros((2, 2)), np.zeros((2, 3, 2)))
+
+
+def test_dimension_mismatch():
+    """Centre/covariance dimension mismatch raises ValueError."""
+    with pytest.raises(ValueError, match="dimensionality"):
+        coef_from_cov_grad(np.zeros((2, 2)), np.eye(3)[np.newaxis].repeat(2, axis=0))
+
+
+def test_singular_cov_returns_nan():
+    """Singular covariance returns NaN coefficients, matching coef_from_cov."""
+    centers = np.zeros((1, 2))
+    cov_sing = np.zeros((1, 2, 2))
+
+    coefs, vjp = coef_from_cov_grad(centers, cov_sing)
+    ref = coef_from_cov(centers, cov_sing)
+    np.testing.assert_array_equal(np.isnan(coefs), np.isnan(ref))
+
+    grad_X, grad_cov = vjp(np.ones_like(coefs))
+    assert np.all(np.isnan(grad_X))
+    assert np.all(np.isnan(grad_cov))