MuellerSeb commented 4 years ago

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:
$R=(\alpha_1,\ldots,\alpha_M)$
where M is the number of 2D sub-spaces (from main axis):
$M=\binom{N}{2}=\frac{N\cdot(N-1)}{2}$
each angle describes the rotation in the x_i-x_j sub-space where the indices are given in the following order:
$\begin{align*} I_{1} & =(1,2)\\ I_{2} & =(1,3)\\ \vdots\\ I_{N-1} & =(1,N)\\ I_{N} & =(2,3)\\ \vdots\\ I_{M} & =(N-1,N) \end{align*}$
the rotation-matrix for a given angle in the i-j plane is given by:
$\left(\text{Rot}_{\alpha}^{(i,j)}\right)_{mn}=\begin{cases} \cos\alpha & m=n=i\\ \cos\alpha & m=n=j\\ -\sin\alpha & m=i,n=j\\ \sin\alpha & m=j,n=i\\ 1 & m=n\\ 0 & \text{else} \end{cases}$
consequently, the rotation matrix to bring unrotated coordinates into the desired rotation is given by:
$\begin{align*} \text{Rot}_{R} & =\text{Rot}_{\alpha_{M}}^{I_{M}}\cdot\ldots\cdot\text{Rot}_{\alpha_{1}}^{I_{1}}\\ \, & =\overset{\curvearrowleft}{\prod_{i=1}^{M}}\text{Rot}_{\alpha_{i}}^{I_{i}} \end{align*}$
since we are interested in derotating given points, we calculate the derotation matrix with:
$\begin{align*} \text{DeRot}_{R} & =\left(\text{Rot}_{R}\right)^{-1}\\ \, & =\text{Rot}_{-\alpha_{1}}^{I_{1}}\cdot\ldots\cdot\text{Rot}_{-\alpha_{M}}^{I_{M}}\\ \, & =\overset{\curvearrowright}{\prod_{i=1}^{M}}\text{Rot}_{-\alpha_{i}}^{I_{i}} \end{align*}$

Anisotropy

Anisotropy is given by N-1 anisotropy factors

$A=(e_{x_2},\ldots,e_{x_N})$

We can also formulate a transformation matrix by:

$\begin{align*} \text{AnIso}_{A} & =\text{diag}(1,e_{x_{2}},\ldots,e_{x_{N}})\\ \text{Iso}_{A} & =\text{AnIso}^{-1}\\ & =\text{diag}(1,e_{x_{2}}^{-1},\ldots,e_{x_{N}}^{-1})\\ e_{x_{i}} & =\frac{\ell_{x_{i}}}{\ell}\\ \ell & =\ell_{x_{1}}\text{ main len-scale} \end{align*}$

normalized lag

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

$\tilde{\underline{h}}=\text{Iso}_{A}\cdot\text{DeRot}_{R}\cdot\underline{h}$

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 commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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: https://github.com/GeoStat-Framework/GSTools/issues/77

LSchueler commented 4 years ago

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

MuellerSeb commented 4 years ago

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 commented 3 years ago

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

MuellerSeb commented 3 years ago

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

GeoStat-Framework / PyKrige

[Proposal] Rotation and Anisotropy in N-D #138

Rotation

Anisotropy

normalized lag

Conclusion