Merge pull request #151 from GeoStat-Framework/pseudo_inv

Use pseudo inverse to prevent singular matrices

.travis.yml

@@ -1,3 +1,6 @@
+# do not build pull request from base repo twice (only build PRs of forks)
+if: type != pull_request OR fork = true
+
 language: python
 python: 3.8
 

pykrige/core.py

@@ -25,10 +25,15 @@
 import numpy as np
 from scipy.spatial.distance import pdist, squareform, cdist
 from scipy.optimize import least_squares
+import scipy.linalg as spl
+
 
 eps = 1.0e-10  # Cutoff for comparison to zero
 
 
+P_INV = {"pinv": spl.pinv, "pinv2": spl.pinv2, "pinvh": spl.pinvh}
+
+
 def great_circle_distance(lon1, lat1, lon2, lat2):
     """Calculate the great circle distance between one or multiple pairs of
     points given in spherical coordinates. Spherical coordinates are expected
@@ -674,7 +679,13 @@ def _calculate_variogram_model(
 
 
 def _krige(
-    X, y, coords, variogram_function, variogram_model_parameters, coordinates_type
+    X,
+    y,
+    coords,
+    variogram_function,
+    variogram_model_parameters,
+    coordinates_type,
+    pseudo_inv=False,
 ):
     """Sets up and solves the ordinary kriging system for the given
     coordinate pair. This function is only used for the statistics calculations.
@@ -694,6 +705,11 @@ def _krige(
     coordinates_type: str
         type of coordinates in X array, can be 'euclidean' for standard
         rectangular coordinates or 'geographic' if the coordinates are lat/lon
+    pseudo_inv : :class:`bool`, optional
+        Whether the kriging system is solved with the pseudo inverted
+        kriging matrix. If `True`, this leads to more numerical stability
+        and redundant points are averaged. But it can take more time.
+        Default: False
 
     Returns
     -------
@@ -755,15 +771,23 @@ def _krige(
     b[n, 0] = 1.0
 
     # solve
-    res = np.linalg.solve(a, b)
+    if pseudo_inv:
+        res = np.linalg.lstsq(a, b, rcond=None)[0]
+    else:
+        res = np.linalg.solve(a, b)
     zinterp = np.sum(res[:n, 0] * y)
     sigmasq = np.sum(res[:, 0] * -b[:, 0])
 
     return zinterp, sigmasq
 
 
 def _find_statistics(
-    X, y, variogram_function, variogram_model_parameters, coordinates_type
+    X,
+    y,
+    variogram_function,
+    variogram_model_parameters,
+    coordinates_type,
+    pseudo_inv=False,
 ):
     """Calculates variogram fit statistics.
     Returns the delta, sigma, and epsilon values for the variogram fit.
@@ -782,6 +806,11 @@ def _find_statistics(
     coordinates_type: str
         type of coordinates in X array, can be 'euclidean' for standard
         rectangular coordinates or 'geographic' if the coordinates are lat/lon
+    pseudo_inv : :class:`bool`, optional
+        Whether the kriging system is solved with the pseudo inverted
+        kriging matrix. If `True`, this leads to more numerical stability
+        and redundant points are averaged. But it can take more time.
+        Default: False
 
     Returns
     -------
@@ -810,6 +839,7 @@ def _find_statistics(
                 variogram_function,
                 variogram_model_parameters,
                 coordinates_type,
+                pseudo_inv,
             )
 
             # if the estimation error is zero, it's probably because

pykrige/lib/cok.pyx

@@ -35,7 +35,21 @@ cpdef _c_exec_loop(double [:, ::1] a_all,
 
     cdef double [::1] variogram_model_parameters = np.asarray(pars['variogram_model_parameters'])
 
-    cdef double [::1,:] a_inv = scipy.linalg.inv(a_all)
+    cdef double [::1,:] a_inv = np.asfortranarray(np.empty_like(a_all))
+
+
+    if pars['pseudo_inv']:
+        if pars['pseudo_inv_type'] == "pinv":
+            a_inv = np.asfortranarray(scipy.linalg.pinv(a_all))
+        elif pars['pseudo_inv_type'] == "pinv2":
+            a_inv = np.asfortranarray(scipy.linalg.pinv2(a_all))
+        elif pars['pseudo_inv_type'] == "pinvh":
+            a_inv = np.asfortranarray(scipy.linalg.pinvh(a_all))
+        else:
+            raise ValueError('Unknown pseudo inverse method selected.')
+    else:
+        a_inv = np.asfortranarray(scipy.linalg.inv(a_all))
+
 
     for i in range(npt):   # same thing as range(npt) if mask is not defined, otherwise take the non masked elements
         if mask[i]:

pykrige/ok.py

@@ -15,7 +15,7 @@
 ----------
 .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in
     Hydrogeology, (Cambridge University Press, 1997) 272 p.
-.. [2] N. Cressie, Statistics for spatial data, 
+.. [2] N. Cressie, Statistics for spatial data,
    (Wiley Series in Probability and Statistics, 1993) 137 p.
 
 Copyright (c) 2015-2020, PyKrige Developers
@@ -31,6 +31,7 @@
     _initialize_variogram_model,
     _make_variogram_parameter_list,
     _find_statistics,
+    P_INV,
 )
 import warnings
 
@@ -143,19 +144,32 @@ class OrdinaryKriging:
         coordinates, both given in degree. Longitudes are expected in
         [0, 360] and latitudes in [-90, 90]. Default is 'euclidean'.
     exact_values : bool, optional
-        If True, interpolation provides input values at input locations. 
-        If False, interpolation accounts for variance/nugget within input 
-        values at input locations and does not behave as an 
+        If True, interpolation provides input values at input locations.
+        If False, interpolation accounts for variance/nugget within input
+        values at input locations and does not behave as an
         exact-interpolator [2]. Note that this only has an effect if
         there is variance/nugget present within the input data since it is
         interpreted as measurement error. If the nugget is zero, the kriged
         field will behave as an exact interpolator.
+    pseudo_inv : :class:`bool`, optional
+        Whether the kriging system is solved with the pseudo inverted
+        kriging matrix. If `True`, this leads to more numerical stability
+        and redundant points are averaged. But it can take more time.
+        Default: False
+    pseudo_inv_type : :class:`str`, optional
+        Here you can select the algorithm to compute the pseudo-inverse matrix:
+
+            * `"pinv"`: use `pinv` from `scipy` which uses `lstsq`
+            * `"pinv2"`: use `pinv2` from `scipy` which uses `SVD`
+            * `"pinvh"`: use `pinvh` from `scipy` which uses eigen-values
+
+        Default: `"pinv"`
 
     References
     ----------
     .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in
        Hydrogeology, (Cambridge University Press, 1997) 272 p.
-    .. [2] N. Cressie, Statistics for spatial data, 
+    .. [2] N. Cressie, Statistics for spatial data,
        (Wiley Series in Probability and Statistics, 1993) 137 p.
     """
 
@@ -186,7 +200,14 @@ def __init__(
         enable_statistics=False,
         coordinates_type="euclidean",
         exact_values=True,
+        pseudo_inv=False,
+        pseudo_inv_type="pinv",
     ):
+        # config the pseudo inverse
+        self.pseudo_inv = bool(pseudo_inv)
+        self.pseudo_inv_type = str(pseudo_inv_type)
+        if self.pseudo_inv_type not in P_INV:
+            raise ValueError("pseudo inv type not valid: " + str(pseudo_inv_type))
 
         # set up variogram model and parameters...
         self.variogram_model = variogram_model
@@ -335,6 +356,7 @@ def __init__(
                 self.variogram_function,
                 self.variogram_model_parameters,
                 self.coordinates_type,
+                self.pseudo_inv,
             )
             self.Q1 = core.calcQ1(self.epsilon)
             self.Q2 = core.calcQ2(self.epsilon)
@@ -506,6 +528,7 @@ def update_variogram_model(
             self.variogram_function,
             self.variogram_model_parameters,
             self.coordinates_type,
+            self.pseudo_inv,
         )
         self.Q1 = core.calcQ1(self.epsilon)
         self.Q2 = core.calcQ2(self.epsilon)
@@ -619,7 +642,11 @@ def _exec_vector(self, a, bd, mask):
         zero_index = None
         zero_value = False
 
-        a_inv = scipy.linalg.inv(a)
+        # use the desired method to invert the kriging matrix
+        if self.pseudo_inv:
+            a_inv = P_INV[self.pseudo_inv_type](a)
+        else:
+            a_inv = scipy.linalg.inv(a)
 
         if np.any(np.absolute(bd) <= self.eps):
             zero_value = True
@@ -650,7 +677,11 @@ def _exec_loop(self, a, bd_all, mask):
         zvalues = np.zeros(npt)
         sigmasq = np.zeros(npt)
 
-        a_inv = scipy.linalg.inv(a)
+        # use the desired method to invert the kriging matrix
+        if self.pseudo_inv:
+            a_inv = P_INV[self.pseudo_inv_type](a)
+        else:
+            a_inv = scipy.linalg.inv(a)
 
         for j in np.nonzero(~mask)[
             0
@@ -876,6 +907,8 @@ def execute(
                     "variogram_model_parameters",
                     "variogram_function",
                     "exact_values",
+                    "pseudo_inv",
+                    "pseudo_inv_type",
                 ]
             }
 

pykrige/ok3d.py

@@ -14,7 +14,7 @@
 ----------
 .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in
     Hydrogeology, (Cambridge University Press, 1997) 272 p.
-.. [2] N. Cressie, Statistics for spatial data, 
+.. [2] N. Cressie, Statistics for spatial data,
 (Wiley Series in Probability and Statistics, 1993) 137 p.
 
 Copyright (c) 2015-2020, PyKrige Developers
@@ -30,6 +30,7 @@
     _initialize_variogram_model,
     _make_variogram_parameter_list,
     _find_statistics,
+    P_INV,
 )
 import warnings
 
@@ -154,19 +155,32 @@ class OrdinaryKriging3D:
     enable_plotting : bool, optional
         Enables plotting to display variogram. Default is False (off).
     exact_values : bool, optional
-        If True, interpolation provides input values at input locations. 
-        If False, interpolation accounts for variance/nugget within input 
-        values at input locations and does not behave as an 
+        If True, interpolation provides input values at input locations.
+        If False, interpolation accounts for variance/nugget within input
+        values at input locations and does not behave as an
         exact-interpolator [2]. Note that this only has an effect if
         there is variance/nugget present within the input data since it is
         interpreted as measurement error. If the nugget is zero, the kriged
         field will behave as an exact interpolator.
+    pseudo_inv : :class:`bool`, optional
+        Whether the kriging system is solved with the pseudo inverted
+        kriging matrix. If `True`, this leads to more numerical stability
+        and redundant points are averaged. But it can take more time.
+        Default: False
+    pseudo_inv_type : :class:`str`, optional
+        Here you can select the algorithm to compute the pseudo-inverse matrix:
+
+            * `"pinv"`: use `pinv` from `scipy` which uses `lstsq`
+            * `"pinv2"`: use `pinv2` from `scipy` which uses `SVD`
+            * `"pinvh"`: use `pinvh` from `scipy` which uses eigen-values
+
+        Default: `"pinv"`
 
     References
     ----------
     .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in
        Hydrogeology, (Cambridge University Press, 1997) 272 p.
-    .. [2] N. Cressie, Statistics for spatial data, 
+    .. [2] N. Cressie, Statistics for spatial data,
        (Wiley Series in Probability and Statistics, 1993) 137 p.
     """
 
@@ -199,7 +213,14 @@ def __init__(
         verbose=False,
         enable_plotting=False,
         exact_values=True,
+        pseudo_inv=False,
+        pseudo_inv_type="pinv",
     ):
+        # config the pseudo inverse
+        self.pseudo_inv = bool(pseudo_inv)
+        self.pseudo_inv_type = str(pseudo_inv_type)
+        if self.pseudo_inv_type not in P_INV:
+            raise ValueError("pseudo inv type not valid: " + str(pseudo_inv_type))
 
         # set up variogram model and parameters...
         self.variogram_model = variogram_model
@@ -332,6 +353,7 @@ def __init__(
             self.variogram_function,
             self.variogram_model_parameters,
             "euclidean",
+            self.pseudo_inv,
         )
         self.Q1 = core.calcQ1(self.epsilon)
         self.Q2 = core.calcQ2(self.epsilon)
@@ -514,6 +536,7 @@ def update_variogram_model(
             self.variogram_function,
             self.variogram_model_parameters,
             "euclidean",
+            self.pseudo_inv,
         )
         self.Q1 = core.calcQ1(self.epsilon)
         self.Q2 = core.calcQ2(self.epsilon)
@@ -604,7 +627,11 @@ def _exec_vector(self, a, bd, mask):
         zero_index = None
         zero_value = False
 
-        a_inv = scipy.linalg.inv(a)
+        # use the desired method to invert the kriging matrix
+        if self.pseudo_inv:
+            a_inv = P_INV[self.pseudo_inv_type](a)
+        else:
+            a_inv = scipy.linalg.inv(a)
 
         if np.any(np.absolute(bd) <= self.eps):
             zero_value = True
@@ -635,7 +662,11 @@ def _exec_loop(self, a, bd_all, mask):
         kvalues = np.zeros(npt)
         sigmasq = np.zeros(npt)
 
-        a_inv = scipy.linalg.inv(a)
+        # use the desired method to invert the kriging matrix
+        if self.pseudo_inv:
+            a_inv = P_INV[self.pseudo_inv_type](a)
+        else:
+            a_inv = scipy.linalg.inv(a)
 
         for j in np.nonzero(~mask)[
             0