[Proposal] Rotation and Anisotropy in N-D

[MuellerSeb]

In order to provide N-D kriging, we should think about, what rotation and anisotropy means in higher dimensions.

My proposal is:

Rotation

• in N-D we can describe rotation with M angles:
  
  where M is the number of 2D sub-spaces (from main axis):
  

• each angle describes the rotation in the x_i-x_j sub-space where the indices are given in the following order:
  

• the rotation-matrix for a given angle in the i-j plane is given by:
  

• consequently, the rotation matrix to bring unrotated coordinates into the desired rotation is given by:
  

• since we are interested in derotating given points, we calculate the derotation matrix with:
  

Anisotropy

Anisotropy is given by N-1 anisotropy factors

We can also formulate a transformation matrix by:

normalized lag

The normalized lag (Isotropified and derotated) is now given by:

Conclusion

The proposed formulation coincides with the current descriptions of rotation and anisotropy in 1D (none), 2D (rotation in x-y plane, 1 aniso-factor) and 3D (Tait–Bryan angles, 2 aniso-factors) and generalizes them to N-D.

What do you think? @rth, @bsmurphy, @LSchueler

[rth]

It is indeed a valid point that I haven't considered before. I'm not competent on this subject, so I can't really give you useful input on this proposal.

In order to provide N-D kriging, we should think about, what rotation and anisotropy means in higher dimensions.

As a side note, we don't necessarily have to resolve this for N-D krigging refactoring. Having one interface for 1D, 2D, 3D would already be great and address most common use cases. In a first approach one could say that N > 3 are not supported, and then in a second phase add support for it if there is interest.

[MuellerSeb]

With the above formulation, it should be quite easy to implement. So we can get rid of different classes for different dimensions.

PS: I updated the issue to fix some wrong calculations.

[bsmurphy]

This looks good to me, and shouldn't be too hard to implement for arbitrary dimensions. At one point, @rth had suggested refactoring the rotation/scaling into a general 'anisotropy transformation' that could be implemented as part of a pipeline, similar to how things are done in scikit-learn. Maybe worth separating the anisotropy stuff out from the whole kriging code (currently the anisotropy adjustments are pretty well baked into the main PyKrige code), would at the very least make code maintenance somewhat easier.

[MuellerSeb]

I am planing to include the transformation matrices in the covariance model in GSTools, so we can simply use it in PyKrige later on. See: GeoStat-Framework/GSTools#77

[LSchueler]

I really like the idea. I've started a GSTools PR to get things started.

[MuellerSeb]

In order to be in line with the Tait-Bryan angles, the signs of the sin functions should alter with the rotation planes to preserve orientation. (meaning the rotation in xz plane is clockwise in 3D).

[MuellerSeb]

Some links:
https://en.wikipedia.org/wiki/Givens_rotation
https://en.wikipedia.org/wiki/Davenport_chained_rotations

[MuellerSeb]

This is now solved in GSTools and will be used in PyKrige in the future:

• CovModel: Update and Refactoring GSTools#109
• [Field] refactoring and update GSTools#112