From 1fd9e42dc8f88c4432c8d223c4a274865c3d3da0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BAlio=20Hoffimann?= Date: Tue, 7 Jul 2026 10:07:05 -0300 Subject: [PATCH 1/3] Use SVD factorization in the presence of missing data --- src/krig.jl | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/src/krig.jl b/src/krig.jl index 142e516..e5873e2 100644 --- a/src/krig.jl +++ b/src/krig.jl @@ -67,7 +67,7 @@ function fit(model::KrigingModel, data) nfun, miss = setlhs!(model, LHS, data) # factorize LHS - FHS = lhsfactorize(model, LHS) + FHS = lhsfactorize(model, LHS, miss) # record Kriging state state = KrigingState(data, LHS, RHS, FHS, nfun, miss) @@ -139,24 +139,26 @@ function setlhs!(model::KrigingModel, LHS, data) end # choose appropriate factorization of LHS -lhsfactorize(model::GeoStatsModel, LHS) = _lhsfactorize(model.fun, LHS) - -# enforce Bunch-Kaufman factorization for symmetric functions -_lhsfactorize(::Variogram, LHS) = bunchkaufman!(Symmetric(LHS), check=false) -_lhsfactorize(::Covariance, LHS) = bunchkaufman!(Symmetric(LHS), check=false) +lhsfactorize(model::GeoStatsModel, LHS, miss) = _lhsfactorize(model.fun, LHS, miss) # enforce SVD factorization for rank-deficient matrices # use QRIteration to avoid LAPACK bug: https://github.com/Reference-LAPACK/lapack/issues/672 -_lhsfactorize(::Transiogram, LHS) = svd!(LHS, alg=QRIteration()) +_lhsfactorize(::Transiogram, LHS, miss) = svd!(LHS, alg=QRIteration()) # choose appropriate factorization for other functions (e.g., CompositeFunction) -function _lhsfactorize(fun::GeoStatsFunction, LHS) - if issymmetric(fun) - # enforce Bunch-Kaufman factorization - bunchkaufman!(Symmetric(LHS), check=false) +function _lhsfactorize(fun::GeoStatsFunction, LHS, miss) + if isempty(miss) + if issymmetric(fun) + # enforce Bunch-Kaufman factorization + bunchkaufman!(Symmetric(LHS), check=false) + else + # fallback to LU factorization + lu!(LHS, check=false) + end else - # fallback to LU factorization - lu!(LHS, check=false) + # fallback to SVD factorization for rank-deficient matrices + # use QRIteration to avoid LAPACK bug: https://github.com/Reference-LAPACK/lapack/issues/672 + svd!(LHS, alg=QRIteration()) end end From 6eab8ac1f4c5223c19a1efc13b4dd04ec8289c1d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BAlio=20Hoffimann?= Date: Tue, 7 Jul 2026 10:16:30 -0300 Subject: [PATCH 2/3] Add missing deps to benchmarks --- benchmark/Project.toml | 1 + benchmark/benchmarks.jl | 1 + 2 files changed, 2 insertions(+) diff --git a/benchmark/Project.toml b/benchmark/Project.toml index 62ff9a1..13608a4 100644 --- a/benchmark/Project.toml +++ b/benchmark/Project.toml @@ -1,4 +1,5 @@ [deps] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +GeoStatsFunctions = "6771c435-bc22-4842-b0c3-41852a255103" GeoTables = "e502b557-6362-48c1-8219-d30d308dcdb0" Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa" diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index 11925dd..e92daf5 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -1,5 +1,6 @@ using BenchmarkTools using GeoStatsModels +using GeoStatsFunctions using GeoTables using Meshes From 7d8d33ff86c1f38542e4ee76071903ad13e4236d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BAlio=20Hoffimann?= Date: Tue, 7 Jul 2026 11:33:25 -0300 Subject: [PATCH 3/3] More refactoring --- src/krig.jl | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/krig.jl b/src/krig.jl index e5873e2..b207ace 100644 --- a/src/krig.jl +++ b/src/krig.jl @@ -12,16 +12,15 @@ abstract type KrigingModel <: GeoStatsModel end Base.range(model::KrigingModel) = range(model.fun) """ - KrigingState(data, LHS, RHS, FHS, nfun, miss) + KrigingState(data, FHS, RHS, nfun, miss) A Kriging state stores information needed to perform estimation at any given geometry. """ -mutable struct KrigingState{D<:AbstractGeoTable,L,R,F} +mutable struct KrigingState{D<:AbstractGeoTable,F,R} data::D - LHS::L - RHS::R FHS::F + RHS::R nfun::Int miss::Vector{Int} end @@ -67,10 +66,10 @@ function fit(model::KrigingModel, data) nfun, miss = setlhs!(model, LHS, data) # factorize LHS - FHS = lhsfactorize(model, LHS, miss) + FHS = lhsfactorize!(model, LHS, miss) # record Kriging state - state = KrigingState(data, LHS, RHS, FHS, nfun, miss) + state = KrigingState(data, FHS, RHS, nfun, miss) FittedKriging(model, state) end @@ -139,14 +138,14 @@ function setlhs!(model::KrigingModel, LHS, data) end # choose appropriate factorization of LHS -lhsfactorize(model::GeoStatsModel, LHS, miss) = _lhsfactorize(model.fun, LHS, miss) +lhsfactorize!(model::GeoStatsModel, LHS, miss) = _lhsfactorize!(model.fun, LHS, miss) # enforce SVD factorization for rank-deficient matrices # use QRIteration to avoid LAPACK bug: https://github.com/Reference-LAPACK/lapack/issues/672 -_lhsfactorize(::Transiogram, LHS, miss) = svd!(LHS, alg=QRIteration()) +_lhsfactorize!(::Transiogram, LHS, miss) = svd!(LHS, alg=QRIteration()) # choose appropriate factorization for other functions (e.g., CompositeFunction) -function _lhsfactorize(fun::GeoStatsFunction, LHS, miss) +function _lhsfactorize!(fun::GeoStatsFunction, LHS, miss) if isempty(miss) if issymmetric(fun) # enforce Bunch-Kaufman factorization