Use pseudo inverse to prevent singular matrices

[MuellerSeb]

This is adressing #87, #150 and others.
The pseudo-inverse eventually solves the system by the least square error solution of the kriging equation. This helps to solve over-/under-determined systems.

If there are redundant points, their values will be averaged in the resulting krige field.

This feature can be pluged in via a switch in the init routines of Ordinary and Universal kriging in 2D/3D.

[MuellerSeb]

And another note:
The example uses a linear model, that is not a valid one in 2D.

[rth]

Are we confident that the matrix we are trying to invert is always Hermetian/symmetric?

Also at least in _exec_vector it looks like we could solve the linear system directly with scipy.linalg.solve which should be faster than inverting the matrix, or am I missing something?

[MuellerSeb]

Yepp, the symmetry of the kriging matrix is a property we can use. pinvh is a lot faster than pinv for big matrices.

We could also solve it, but we should use scipy.linalg.lstsq, to guarantee a unique solution (should coincide with the pseudo-inverse solution)

I could check, if that solves the problem above.

[MuellerSeb]

Problem is:
linalg.solve and linalg.lstsq don't support masked arrays. So this would need more refactoring.

[MuellerSeb]

For some strange reason, one test is failing with pinvh but succeeds with pinv... so there are some numerical errors occurring. I would stick to pinv then.

@rth since the inversion in _exec_vector and _exec_loop is only done once, I think it would be ok to not use linalg.solve.

I would add a switch to use the pseudo-inverse, so the current behavior stays the same and user can enable the pseudo-inverse, if they need it.

[MuellerSeb]

I will implement a switch inv in the execute method to select the used method.
Either:

• inv="inv" (default) use the inverse as before
• inv="pinv" use pinv from scipy which uses lstsq to compute the pseudo-inverse
• inv="pinv2" use pinv2 from scipy which uses singular value decomposition to compute the pseudo-inverse
• inv="pinvh" use pinvh from scipy which uses the eigenvalues to compute the pseudo-inverse

Then everything will stay as it is and the users can select other routines for the inversion.

[MuellerSeb]

I now decided to add two new keywords to the Kriging class __init__ routines, following GeoStat-Framework/GSTools#97:

    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:`int`, optional
        Here you can select the algorithm to compute the pseudo-inverse matrix:

            * `1`: use `pinv` from `scipy` which uses `lstsq`
            * `2`: use `pinv2` from `scipy` which uses `SVD`
            * `3`: use `pinvh` from `scipy` which uses eigen-values

So the current behavior is conserved and the usage of the pseudo-inverse can be opted in.

@rth : Are you OK with that?

[MuellerSeb]

@rth ping. :-)

[MuellerSeb]

I think, I'll just merge this by next week!

[rth]

Thanks @MuellerSeb , sorry for late response. Generally LGTM but I haven't done a very thorough review.

One comment is that it might be better to use string options for pseudo_inv_type, as UniversalKrigging(pseudo_inv_type='pinvh') is more explicit/readable than UniversalKrigging(pseudo_inv_type=3).