Zero kriging-matrix diagonal so the nugget is not placed on it#2922
Merged
brendancol merged 2 commits intoJun 4, 2026
Conversation
The semivariogram has gamma(0)=0 by definition; vario_func(0) returns the nugget c0 (the one-sided limit as h->0+), not the value at h=0. _build_kriging_matrix set the variogram block to vario_func(D) where D has a zero diagonal, so a non-zero nugget landed on the matrix diagonal. That forced exact interpolation of noisy data and biased the kriging variance downward. Force the variogram-block diagonal to 0 with np.fill_diagonal. Add two regression tests: the matrix diagonal stays 0 under a non-zero nugget, and a non-zero nugget now smooths predictions and lifts the variance above the nugget floor.
brendancol
commented
Jun 4, 2026
brendancol
left a comment
Contributor
Author
There was a problem hiding this comment.
PR Review: Zero kriging-matrix diagonal so the nugget is not placed on it
Blockers
None.
Suggestions
None.
Nits
None.
What looks good
- The fix matches the standard ordinary-kriging formulation: the semivariogram is gamma(0)=0, and the nugget is the one-sided limit as h -> 0+, so it should not sit on the matrix diagonal. PyKrige does the same.
- The change is minimal and well-commented, touching only
_build_kriging_matrix. _build_kriging_matrixruns in numpy on the host before backend dispatch, so all four backends (numpy / cupy / dask+numpy / dask+cupy) consume the correctedK_inv. The fix is uniform without per-backend edits.- The existing post-failure regularizer (
K[:n,:n] += 1e-10 * np.eye(n)) is a separate Tikhonov term applied only wheninvraises; it is unaffected and still correct. - Both new tests were confirmed to fail on the pre-fix code and pass after, so they are genuine regression guards rather than tautologies.
Notes
- The new tests exercise
_build_kriging_matrix/_kriging_predictdirectly in numpy rather than across all four backends. That is appropriate here: the corrected line is shared host code, and the existingtest_cupy_*/test_dask_cupy_*parity tests already cover dispatch. Adding nugget-specific backend variants would not increase coverage of the fix.
Checklist
- Algorithm matches reference (ordinary kriging, gamma(0)=0)
- All implemented backends consume the corrected matrix
- NaN handling unchanged
- Regression tests fail before / pass after the fix
- No premature materialization or unnecessary copies
- Benchmark not needed (bug fix, no perf-sensitive change)
- README feature matrix unchanged (no new API)
- Docstrings/comments present and accurate
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes #2915
_build_kriging_matrixset the variogram block tovario_func(D), whereDhas a zero diagonal.vario_func(0)returns the nuggetc0, so a non-zero nugget ended up on the matrix diagonal. The semivariogram isgamma(0) = 0by definition (the nugget is the one-sided limit ash -> 0+), and standard ordinary-kriging implementations zero that diagonal.The bug only bit when the fitted nugget was non-zero: it forced exact interpolation of noisy data and shrank the kriging variance. The existing tests use trend-dominated data that fits a near-zero nugget, so they never exercised it.
Changes:
np.fill_diagonalin_build_kriging_matrix.Backends:
_build_kriging_matrixis shared by all four backends (numpy / cupy / dask+numpy / dask+cupy), which buildK_invon the host before dispatch, so the fix applies uniformly.Test plan:
pytest xrspatial/tests/test_interpolation.py -k "Kriging or kriging"(16 passed, includes cupy and dask+cupy parity)