Natural Cubic Spline Interpolation

[yfnaji]

Natural Cubic Spline Interpolation

This Pull Request implements the cubic spline interpolation as one of the polynomial interpolators mentioned in issue #5.

The implementation was benchmarked against SciPy’s CubicSpline.

Mathematical Background

We want to construct a cubic spline $S(x)$ that interpolates the points

$$\left(x_0, y_0\right), \left(x_1, y_1\right), \ldots,\left(x_n, y_n\right)$$

Since we are working with a natural cubic spline, we define

$$ z_i = \frac{d^2S(x_i)}{dx^2}, $$

the second derivative of the spline at each interpolation point and set

$$ z_0 = z_n = 0 $$

The full derivation of the cubic spline formulation is fairly involved, so some details are omitted here. For reference, the derivation can be found in the lecture notes by T. Gambill (2011) Interpolation/Splines (slides 14-27), which served as the main source for this implementation.

The equation of a natural cubic spline on the interval $x\in\left[x_i, x_{i+1}\right)$ is given by:

$$ S_i(x) = \frac{z_i}{6h_i}\left(x_{i+1} - x\right)^3 +\frac{z_{i+1}}{6h_i}\left(x - x_i\right)^3 +\left(\frac{y_{i+1}}{h_i} - \frac{h_i}{6}z_{i+1}\right)\left(x - x_i\right) +\left(\frac{y_{i}}{h_i} - \frac{h_i}{6}z_{i}\right)\left(x_{i+1} - x\right) $$

(Equation (1))

The next step is to determine suitable values for $z_i$ in order to construct the natural cubic spline.

Differentiating $S_i(x)$ yields

$$ S_i'\left(x\right) = -\frac{z_i}{2h_i}\left(x_{i+1}-x\right)^2 +\frac{z_{i+1}}{2h_i}\left(x-x_i\right)^2 +\frac{y_{i+1}}{h_i} -\frac{h_i}{6}z_{i+1} -\frac{y_i}{h_i} +\frac{h_i}{6}z_i $$

At the interpolation point $x_i$, the spline segments must join smoothly. This requires the gradients at the interpolation point to be equal on both sides. Hence, the following condition is enforced:

Substituting these expressions and rearranging gives the relation

$$ h_{i-1}z_{i-1} + 2\left(h_i + h_{i-1}\right)z_i + h_iz_{i+1} = 6\left(b_i - b_{i-1}\right) $$

where

$$ b_i=\frac{y_{i+1} - y_i}{h_i} $$

This leads to an $(n-1)\times(n-1)$ system of equations:

where

$$ u_i = 2\left(h_i + h_{i+1}\right) $$

Now, we can solve the above equation to obtain the values of $z_i$.

All diagonal and off-diagonal entries of $A$ are positive, implying that $A$ is positive definite. This allows the use of an $LDL^T$ decomposition. A stable inverse of $A$ can be computed using the $LDL^T$ decomposition as follows:

$$ A=LDL^T\Rightarrow A^{-1} = \left(LDL^T\right)^{-1} = (L^{-1})^TD^{-1}L^{-1} $$

The idea behind this approach is that it will be easier the compute the inverses of $L$ and $D$ to then in turn compute the inverse of $A$.

The values of $z_i$ can now be obtained, and these can be substituted into equation (1) to construct the spline functions. Since we are dealing with natural cubic splines, the boundary conditions are $z_0 = z_n = 0$. The $(n-1) \times (n-1)$ system described above provides the values of $z_i$ for $i = 1, 2, \ldots, n-1$.

Notes on implementation

In this implementation, all $0$ entries in matrices have been omitted wherever possible. For example, this applies to a lower triangular matrix (in our case, $L^{-1}$ as mentioned above):

$$ A = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 2 & 3 & 0 & 0 \\ 4 & 5 & 6 & 0 \\ 7 & 8 & 9 & 10 \\ \end{bmatrix} $$

would be represented as an array of arrays

[
    [1],
    [2,3],
    [4,5,6],
    [7,8,9,10]
]

Similarly, $A^T$, an upper triangular matrix, would be represented as

[
    [1,2,4,7],
    [3,5,8],
    [6,9],
    [10]
]

For diagonal matrices, an even simpler approach has been used:

$$ \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 \\ 0 & 0 & 3 & 0 \\ 0 & 0 & 0 & 4 \\ \end{bmatrix} $$

is represented as

[1,2,3,4]

A special note regarding the computation of $L$ in $A = LDL^T$: the diagonal entries of $L$ are always equal to 1, and only the subdiagonal entries can be non-zero. That is, $L$ always has the form

$$ L = \begin{bmatrix} 1 & 0 & 0 & 0 \\ a & 1 & 0 & 0 \\ 0 & b & 1 & 0 \\ 0 & 0 & c & 1 \end{bmatrix} $$

Since this structure is consistent, the diagonal entries can be omitted, and the lower tridiagonal matrix $L$ can be represented simply as

[a, b, c]

These optimizations have been noted in the comments above the functions they have been implemented in.

Amendments to the InterpolationIndex trait

Additional traits have been added to the type Delta:

• Add: Delta + Delta = Delta
• Sub: Delta - Delta = Delta
• IntoValue: Converts a Delta into an InterpolationValue. This is required for cubic splines, as the implementation involves arithmetic operations between Delta and InterpolationValue types. Note: The implementation for numeric types differs from that for date types
• Copy: Required when unpacking Delta values from vectors

Amendments to the InterpolationValue trait

• MulAssign: value_1 *= value_2
• num::FromPrimitive: In order to take an f64 value and convert it to the appropriate InterpolationValue value

Removing Unsigned integers from InterpolationIndex implementation

u8, u16, u32, u64, u128, and usize have been removed from the trait implementation of InterpolationIndex because cubic splines can produce negative Delta values. The trait implementation for i8, i16, i32, i64, i128, and isize should be sufficient for cases where integers are used as indices.

[avhz]

Hi @yfnaji :) Thanks for this one, splines have been needed for a while!

I'll take a closer look tomorrow afternoon, and hopefully merge it

[yfnaji]

Hey @avhz - just saw the Clippy warnings, I am working on them now, but still feel free to review and comment on this PR (:

[avhz]

@yfnaji
I did a bit of refactoring of the Curves today, which required adding Send + Sync to the interpolators. Doesn't look like it created any conflicts, so just a heads up :)

[avhz]

A couple points

[avhz]

I'm a bit unsure about fixing a unit like this, why seconds?

[yfnaji]

At some point we need to apply some arithmetic operations to the Delta type and so a complication arises when apply arithmetic operations between an f64 (the ValueType for time::Date ) and time::Duration e.g. in lower_tri_transpose_inv_times_diag_inv() (line 207):

ValueType::one() / diagonal[j].into_value()

where diagonal[j] is a Delta type.

As for converting time::Duration (to f64), the only possible method that provides this is as_seconds_f64().

We could perhaps create our own conversion?

[avhz]

If I understand correctly, we can't do this because then we are setting the units in seconds which would make no sense on a yield curve, for example.

e.g. for a 1 day delta on a curve, the value would be 86,400.

[avhz]

Also unsure this is needed ?

e.g. i8 as i8 is already possible.

[yfnaji]

Amended this, however, the trait is still required as we have implemented this trait for the time:Date IndexType so it must be implemented for all IndexType's