More Efficient Matrix Inversion for Spline

[cmdty]

Linear system (18) of the maximum smoothness spline document involves a block matrix whose sub-matrices are sparse. This structure potentially lends itself to a much more efficient (in terms of both memory and CPU usage) solution strategy than the naive approach in the current implementation. A more efficient approach would be something like:

1. Formulate the inversion of the block matrix in (18) as a 2 x 2 block matrix with sub-matrices a function of the sub-matrices and inverses of these sub-matrices of the original matrix. The fact that the bottom right sub-matrix is a zero matrix should hopefully simplify things.
2. Use the block-diagonal sparseness of the sub-matrices to implement efficient inversions, multiplications, subtractions and additions to evaluate the matrix inverse formulated in the previous step.

[cmdty]

Looks like this isn’t as straight forward as previously thought due to the matrix H being singular, due to the columns of zeros it contains.