Refactor factorization in the presence of missing data

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"

benchmark/benchmarks.jl

@@ -1,5 +1,6 @@
 using BenchmarkTools
 using GeoStatsModels
+using GeoStatsFunctions
 using GeoTables
 using Meshes
 

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)
+  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,24 +138,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