pylops.signalprocessing.InterpCubicSpline

class pylops.signalprocessing.InterpCubicSpline(dims, iava, bc_type='natural', axis=-1, dtype='float64', name='S')[source]

Cubic spline interpolation operator.

Interpolate a regularly sampled input vector along axis at the fractional positions iava using a cubic spline, i.e., a \(C^{2}\)-continuous piecewise third order polynomial.

Currently, only cubic splines with natural boundary conditions are supported, i.e., the second derivatives at the first and last sampling point are both zero.

Note

The vector iava should contain unique values only. If the same fractional index is present multiple times, an error will be raised. Elements that exceed the last index dims[axis] - 1 are clipped to the closest float value right below dims[axis] - 1 to avoid extrapolation.

Parameters:
dimslist or int

Number of samples for each dimension. A cubic spline requires dims[axis] > 2.

iavalist or numpy.ndarray

Floating indices of the locations to which the spline should interpolate.

bc_typestr, optional

The type of boundary condition. Currently, only "natural" is supported.

axisint, optional

Axis along which interpolation is applied. By default, the interpolation is carried out over the last axis.

dtypenumpy.dtype-like, optional

The data type of the input and output arrays. For complex input, both the real and the imaginary parts are interpolated separately. Only double precision versions of numpy.inexact are supported, i.e., either "float64" (default) or "complex128". Multiplication of the operator with data with less precise data types will result in a type promotion.

namestr, optional

Name of operator (to be used by pylops.utils.describe.describe).

Attributes:
dimstuple

Shape of the array after the adjoint, but before flattening. For example, x_reshaped = (Op.H * y.ravel()).reshape(Op.dims).

dimsdtuple

Shape of the array after the forward, but before flattening. For example, y_reshaped = (Op * x.ravel()).reshape(Op.dimsd).

shapetuple

Operator shape.

Raises:
ValueError

If dims[axis] < 2.

ValueError

If the iava contains duplicate values.

NotImplementError

If bc_type != "natural".

TypeError

If dtype is neither numpy.float64 nor numpy.complex128.

See also

pylops.signalprocessing.Interp

General interpolation operator

scipy.interpolate.CubicSpline

An equivalent implementation of the forward operator \(\mathbf{S}\mathbf{x}\)

Notes

Cubic spline interpolation of an \(\left(N\times 1\right)\)-vector \(\mathbf{x}\) at the \(L\) fractional positions iava can be represented as \(\mathbf{y}=\mathbf{S}\mathbf{x}\) where \(\mathbf{S}\) is the cubic spline operator.

The \(\left(N\times 1\right)\)-grid-point vector at which \(\mathbf{x}\) was equidistantly sampled is denoted by the “knot” vector \(\mathbf{k}\), i.e., \(\mathbf{x}_{j} = f\left(\mathbf{k}_{j}\right)\). Furthermore, the vector iava is denoted by \(t\) in the following mathematical expressions.

When iava[i] (\(t_{i}\)) is mapped to the \(j\)-th knot-to-knot-interval \(\mathbf{k}_{j}\le t_{i}\lt\mathbf{k}_{j + 1}\), the corresponding polynomial \(s_{j}\left(t\right)\) can be expressed as

\[s_{j}\left(t\right) = p_{j,0}\left(t\right)\cdot\mathbf{x}_{j} + p_{j,1}\left(t\right)\cdot\mathbf{x}_{j + 1} + p_{j,2}\left(t\right)\cdot\mathbf{m}_{j} + p_{j,3}\left(t\right)\cdot\mathbf{m}_{j + 1}\]

where the \(\left(N\times 1\right)\)-vector \(\mathbf{m}\) holds, the second order derivatives of the cubic spline at its knots and the base polynomials are given by

  • \(p_{j,0}\left(t\right) = \mathbf{k}_{j + 1} - t\)

  • \(p_{j,1}\left(t\right) = t - \mathbf{k}_{j}\)

  • \(p_{j,2}\left(t\right) = \frac{1}{6}\cdot\left(\left(\mathbf{k}_{j + 1} - t\right)^{3} - \left(\mathbf{k}_{j + 1} - t\right)\right)\)

  • \(p_{j,3}\left(t\right) = \frac{1}{6}\cdot\left(\left(t - \mathbf{k}_{j}\right)^{3} - \left(t - \mathbf{k}_{j}\right)\right)\)

Applying this 4-term linear combination to all points in iava leads to the matrix representation

\[\begin{split}\mathbf{S}\mathbf{x} = \mathbf{P}\begin{bmatrix} \mathbf{x} \\ \mathbf{m} \end{bmatrix}\end{split}\]

where each row of the very sparse \(\left(L\times 2N\right)\)-base-polynomial-matrix \(\mathbf{P}\) has only 4 non-zero entries:

\[\mathbf{P}_{i} = \begin{bmatrix} \dots & 0 & \dots & p_{j,0}\left(t_{i}\right) & p_{j,1}\left(t_{i}\right) & \dots & 0 & \dots & p_{j,2}\left(t_{i}\right) & p_{j,3}\left(t_{i}\right) & \dots & 0 & \dots \end{bmatrix}\]

at the indices \(j\), \(j + 1\), \(j + N\), and \(j + N + 1\), respectively.

\(\mathbf{P}\) forms the first linear operator of \(\mathbf{S}\) that carries out the linear combination and requires the vectors \(\mathbf{x}\) and \(\mathbf{m}\) as combination weights. Those need to be obtained from the second linear operator \(\mathbf{F}\)

\[\begin{split}\mathbf{F}\mathbf{x} = \begin{bmatrix} \mathbf{x} \\ \mathbf{m} \end{bmatrix}\end{split}\]

For cubic splines, \(\mathbf{F}\) is given by

\[\begin{split}\mathbf{F} = \begin{bmatrix} \mathbf{I}_{N} \\ \mathbf{B}^{-1}\mathbf{D}_{2} \end{bmatrix}\end{split}\]

\(\mathbf{I}_{N}\) is the \(\left(N\times N\right)\)-identity matrix. The \(\left(N\times N\right)\)-tridiagonal or -near-tridiagonal matrix \(\mathbf{B}\) and the \(\left(N\times N\right)\)-second-order-finite- difference matrix \(\mathbf{D}_{2}\) originate from the linear system

\[\mathbf{B}\mathbf{m}=\mathbf{D}_{2}\mathbf{x}\]

that needs to be solved for \(\mathbf{m}\), the second order derivatives of the cubic spline at its knots, by

\[\mathbf{m}=\mathbf{B}^{-1}\mathbf{D}_{2}\mathbf{x}\]

Assuming \(\mathbf{x}\) was sampled at equidistant knots, \(\mathbf{B}\) simplifies to

\[\begin{split}\mathbf{B} = \frac{1}{6}\cdot\begin{bmatrix} \mu_{0} & \lambda_{0} & 0 & 0 & 0 & \dots & 0 & 0 & 0 & 0 & \theta_{0} \\ 1 & 4 & 1 & 0 & 0 & \dots & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 4 & 1 & 0 & \dots & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 4 & 1 & \dots & 0 & 0 & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & 0 & \dots & 1 & 4 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \dots & 0 & 1 & 4 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & \dots & 0 & 0 & 1 & 4 & 1 \\ \theta_{N-1} & 0 & 0 & 0 & 0 & \dots & 0 & 0 & 0 & \lambda_{N-1} & \mu_{N-1} \end{bmatrix}\end{split}\]

The special values \(\mu_{i}\), \(\lambda_{i}\), \(\theta_{i}\) for the first (\(i = 0\)) and last row (\(i = N - 1\)) are determined by the boundary conditions, i.e., the behaviour prescribed at the first and last knot. See below for different boundary conditions and their corresponding special values.

Similarly, the second order finite difference matrix \(\mathbf{D}_{2}\) can be reduced to

\[\begin{split}\mathbf{D}_{2} = \begin{bmatrix} \ & \ & \ & \ & \ & \mathbf{d}_{0} & \ & \ & \ & \ \\ 1 & -2 & 1 & 0 & 0 & \dots & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 & \dots & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & -2 & 1 & \dots & 0 & 0 & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & 0 & \dots & 1 & -2 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \dots & 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & \dots & 0 & 0 & 1 & -2 & 1 \\ \ & \ & \ & \ & \ & \mathbf{d}_{N-1} & \ & \ & \ & \ \end{bmatrix}\end{split}\]

Again, the special rows \(\mathbf{d}_{i}\) depend on the boundary conditions.

For the “natural” boundary condition for which the second derivatives of the spline are exactly zero at the boundaries, the special values and rows are

  • \(\mu_{0} = \mu_{N - 1} = 6\)

  • \(\lambda_{0} = \lambda_{N - 1} = 0\)

  • \(\theta_{0} = \theta_{N - 1} = 0\)

  • \(\mathbf{d}_{0} = \mathbf{d}_{N - 1}\) are all zeroes

and \(\mathbf{B}\) is thus truly tridiagonal.

In summary, the cubic spline interpolation can be decomposed as

\[\mathbf{y}=\mathbf{S}\mathbf{x}=\mathbf{P}\mathbf{F}\mathbf{x}\]

The adjoint operator \(\mathbf{S}^{H}\) can be derived by rearranging the involved operators. All of them are purely real and consequently, a transpose is sufficient. This yields

\[\mathbf{S}^{H}=\mathbf{S}^{T} = \mathbf{F}^{T}\mathbf{P}^{T} = \begin{bmatrix} \mathbf{I}_{N} & {\mathbf{D}_{2}}^{T} \mathbf{B}^{-T} \end{bmatrix} \mathbf{P}^{T}\]

where \(\mathbf{B}^{-T} = \left(\mathbf{B}^{-1}\right)^{T} = \left(\mathbf{B}^{T}\right)^{-1}\).

Computationally,

  • \(\mathbf{D}_{2}\) and \({\mathbf{D}_{2}}^{T}\) can be computed efficiently via finite differences with appropriate boundary handling

  • \(\mathbf{B}\) and \(\mathbf{B}^{T}\) are individually factorized via partially pivoted tridiagonal \(LU\)-decompositions that are kept in memory

  • \(\mathbf{B}^{-1}\mathbf{z}\) and \(\mathbf{B}^{-T}\mathbf{z}\) then rely on highly optimized backward- and forward-substitutions with those factorizations to solve the linear systems

  • \(\mathbf{P}\) is represented as a scipy.sparse.csr_matrix which involves a tiny bit of overhead during operator initialization, but allows for a simple transpose.

References

Methods

__init__(dims, iava[, bc_type, axis, dtype, ...])

adjoint()

apply_columns(cols)

Apply subset of columns of operator

cond([uselobpcg])

Condition number of linear operator.

conj()

Complex conjugate operator

div(y[, niter, densesolver])

Solve the linear problem \(\mathbf{y}=\mathbf{A}\mathbf{x}\).

dot(x)

Matrix-matrix or matrix-vector multiplication.

eigs([neigs, symmetric, niter, uselobpcg])

Most significant eigenvalues of linear operator.

matmat(X)

Matrix-matrix multiplication.

matvec(x)

Matrix-vector multiplication.

reset_count()

Reset counters

rmatmat(X)

Matrix-matrix multiplication.

rmatvec(x)

Adjoint matrix-vector multiplication.

todense([backend])

Return dense matrix.

toimag([forw, adj])

Imag operator

toreal([forw, adj])

Real operator

tosparse()

Return sparse matrix.

trace([neval, method, backend])

Trace of linear operator.

transpose()

Examples using pylops.signalprocessing.InterpCubicSpline

Restriction and Interpolation

Restriction and Interpolation