Performance optimization with LU decomposition

[rth]

Making this @bsmurphy comment a separate issue,

BTW, @rth (and anyone else who thinks about matrices, I don't think about them that much admittedly), I've been wondering if we could use LU decomposition on the kriging matrix to speed up the solution. Any thoughts? I haven't thought about this that much, but since the LHS matrix is the same (just the RHS vector that's changing), we might be able to leverage this for the looping backend...

[rth]

@bsmurphy Thanks for this idea!

The current situation (at least with ordinary kriging) is that to solve the linear kriging system we use,

1. matrix inversion in _exec_loop (loop backend), _exec_vector (vectorized backend) and _c_exec_loop (C loop backend) .

• linalg.solve in _exec_loop_moving_window (loop + moving window backend) and _c_exec_loop_moving_window (C loop backend + mooving window)

The linalg.solve uses a LU decomposition internally , so I guess the questions are,

1. would it be faster to get rid of matrix inversion, and use LU decomposition instead?
2. in the cases when linalg.solve is used repeatedly inside a loop could we precompute the LU decomposition outside of the loop and solve a simpler system in the loop?

Both sound interesting and would be worth exploring (strarting from the pure Python implementations that are easier to change). In any case we need benchmarks to see how it would impact performance and memory use (cf. PR #36) ..

[MuellerSeb]

We should use the properties of the kriging matrix to optimize that.
The kriging equation can always be reformulated to use the covariance-matrix instead of the semi-variogram-matrix (as done now). Then, this part of the matrix is a symmetric and positive definite matrix (if we use a valid covariance model), which could be tackled by the Cholesky decomposition. The full inversion could then be computed by block matrix-inversion as stated here: #52 (comment)