Closed MuellerSeb closed 3 years ago
And another note: The example uses a linear model, that is not a valid one in 2D.
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?
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.
Problem is:
linalg.solve
and linalg.lstsq
don't support masked arrays. So this would need more refactoring.
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.
I will implement a switch inv
in the execute
method to select the used method.
Either:
inv="inv"
(default) use the inverse as beforeinv="pinv"
use pinv
from scipy which uses lstsq
to compute the pseudo-inverseinv="pinv2"
use pinv2
from scipy which uses singular value decomposition to compute the pseudo-inverseinv="pinvh"
use pinvh
from scipy which uses the eigenvalues to compute the pseudo-inverseThen everything will stay as it is and the users can select other routines for the inversion.
I now decided to add two new keywords to the Kriging class __init__
routines, following https://github.com/GeoStat-Framework/GSTools/pull/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?
@rth ping. :-)
I think, I'll just merge this by next week!
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.