Proposal: unify kriging routines

[MuellerSeb]

Since all kriging routines could be reformulated in order to use the covariance matrix instead of the variogram matrix, we could unify the kriging base class to serve as a swiss knife for kriging.

All other kriging routines are then just a shortcut with a limited input parameter list.

[MuellerSeb]

It would also be good to use the pseudo-inverse, if the kriging matrix is singular.
We could check this as suggested here:
https://stackoverflow.com/a/9155489/6696397

When using the Covariance matrix, we could also use SVD to invert the matrix as suggested here:
https://www.quora.com/How-do-I-get-the-inverse-of-a-matrix-using-SVD-in-Python/answer/Samuel-S-Watson-1

The fully inverted kriging matrix could then be calculated by using the inversion rules for block matrices:
https://en.wikipedia.org/wiki/Block_matrix#Block_matrix_inversion

[LSchueler]

Great idea! And nice research, the solution from the stackoverflow link is pretty funny.

What about performance? - Do you have any idea at what expense the flexibility will come?

[MuellerSeb]

The try-except block doesn't include any overhead, when the current implementation succeeds. The pinv takes about 5 times more time to compute the inverse (if we have numpy 1.17, which brings the hermitian keyword):

import numpy as np

a = np.random.random((1000, 1000))
a += a.T

import timeit

%timeit np.linalg.inv(a)
94.8 ms ± 5.18 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit np.linalg.pinv(a)
1.17 s ± 76.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit np.linalg.pinv(a, hermitian=True)
452 ms ± 25.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

[LSchueler]

Actually, I mean the fully flexible approach, not just the pseudo inverse. That's a thing we should definitely do.

[MuellerSeb]

Sorry for getting you wrong. The flexible approach is not including any overhead. At the time of coding I was just not aware, that kriging can always be formulated using the covariance function instead of the variogram function (in case of finite sill and stationary models, which we always assume). So we get rid of the different approaches for creating the kriging matrix, since every method only extends the matrix:

simple -> ordinary -> (ext)drift

And I propose to simply use this hierarchy and extend the matrix if needed.

[MuellerSeb]

As seen in #89, sometimes a LinAlgError isn't even detected, so we should stick to the pinvfrom the beginning.
@LSchueler, do you think it would be ok, to increase the minimal required versions:

• numpy to v1.17.0 (from v1.14.5)
• Cython to v0.29.11 (from v0.28.3)
• scipy to v1.2.2 (from v1.1.0)
  They all support py>=3.5.

[MuellerSeb]

I just noticed, that scipy 1.1.0 already provides a pseudo inverse algorithm for symmetric matrices:
https://docs.scipy.org/doc/scipy-1.1.0/reference/generated/scipy.linalg.pinvh.html#scipy.linalg.pinvh

So we could use this one, without increasing any required versions.

[MuellerSeb]

This was resolved with #97