Padding
How to extend a finite-length vector of data to a virtually infinite-length sequence.
Purpose
We wish to build the real function
There, the sequence \(c\) of coefficients is used to parameterize the function \(f\) and gives us a good many degrees of freedom to shape it to our taste. The adaptability of \(c\) makes it a tool of choice to represent sampled data as the continuously defined function \(f.\) Given a sequence \(y=\left(y[k]\right)_{k\in{\mathbb{Z}}}\) of regularly indexed samples \(y[k]\) for \(k\in{\mathbb{Z}},\) the purpose of the so-called interpolation procedure is to yield a sequence \(c\) such that the interpolation condition \(f(k)=y[k]\) is satisfied.
In case \(\varphi\) is a B-spline, we shall soon see that there exist infinitely many sequences of coefficients that build the same one sequence of samples. Among all these sequences of coefficients , however, there is one that is of particular interest because it can be obtained by a very efficient algorithm to get \(c\) out of \(y.\)
Unfortunately, the theoretical derivations of the algorithm rely on \(y\) being a sequence, which means that infinitely many samples are required. Now, one never has access to a sequence of samples in practice, only to a finite-dimensional vector \({\mathbf{y}}=\left(y_{q}\right)_{q=1}^{K}\in{\mathbb{R}}^{K}\) of samples. To take advantage of the theoretical derivations of the efficient algorithm, it is thus an unavoidable necessity that a method be engineered that converts the vector \({\mathbf{y}}\) into the sequence \(y.\) We call padding the operation that consists in the engineering of the infinite-length subsequences \(\left(y[k]\right)_{k\in{\mathbb{Z}}_{<0}}\) to the left and \(\left(y[k]\right)_{k\in{\mathbb{Z}}_{\geq K}}\) to the right of the provided finite-length sequence (or finite-dimensional vector) \(\left(y[k]\right)_{k=0}^{K-1}=\left(y_{q}\right)_{q=1}^{K}={\mathbf{y}}.\)
The engineering of \({\mathbf{y}}\Rightarrow y\) is application-dependent and every method is fair game. Here are a few cases.
The components of the vector \({\mathbf{y}}\) could represent angular data, in which case one would have to cope with angular wrapping.
Or, \({\mathbf{y}}\) could represent intensity data, in which case one would have to discourage negative intensities in the continuously defined function \(f\) being constructed through the steps \({\mathbf{y}}\Rightarrow y\Rightarrow c\Rightarrow f.\)
Or, one could pretend that all unobserved samples do vanish and take the special value \(0.\)
Or, one could assume that the first observed sample \(y[0]\) has indeed the same value as all unobserved samples that came before, while the last observed sample \(y[K-1]\) has the same value as all unobserved samples that folllow.
Many more padding recipes can be devised. Which one is the best in the context of your application is only a matter of convenience.
Uniqueness
In the sequel, we let \(\varphi\) be a polynomial B-spline of nonnegative integer degree \(n,\) which is a real function \(\beta^{n}:{\mathbb{R}}\rightarrow{\mathbb{R}},x\mapsto\beta^{n}(x).\) For any degree \(n\geq2\) and for \(m\in[1\ldots\left\lfloor n/2\right\rfloor],\) it is known that there exist \(\left\lfloor n/2\right\rfloor\) mutually different, real, negative numbers \(z_{n,m}\) in the open interval \((-1,0)\) that satisfy the relation
These numbers are called poles and are such that \(\sum_{k\in{\mathbb{Z}}}\,\beta^{n}(k)\,z_{n,m}^{-k}=0.\) Because B-splines are even-symmetric, the pole reciprocals \(z_{n,m}^{-1}\in{\mathbb{R}}_{<-1}\) satisfy the same relation.
Suppose now we have identified a sequence \(c\) that verifies the interpolation condition \(\forall k\in{\mathbb{Z}}:y[k]=\sum_{k'\in{\mathbb{Z}}}\,c[k']\,\beta^{n}(k-k').\) Then, it turns out that the sequence \(\left(c'[k']\right)_{k'\in{\mathbb{Z}}}=\left(c[k']+\sum_{m=1}^{\left\lfloor n/2\right\rfloor}\,\left(\lambda_{n,m}^{-}\,z_{n,m}^{-k'}+\lambda_{n,m}^{+}\,z_{n,m}^{k'}\right)\right)_{k'\in{\mathbb{Z}}}\) is also such that \(\forall k\in{\mathbb{Z}}:y[k]=\sum_{k'\in{\mathbb{Z}}}\,c'[k']\,\beta^{n}(k-k'),\) for any choice of the arbitrary constants \(\lambda_{n,m}^{-},\lambda_{n,m}^{+}\in{\mathbb{R}}.\) One concludes that the interpolation condition \(\forall k\in{\mathbb{Z}}:y[k]=f(k)\) alone is not sufficient to make the interpolation task well-defined and that additional requirements are needed.
The splinekit library comes equipped with paddings that guarantee the well-defined interpolation of finite-length data.
Available Paddings
The unobserved samples are, well, unobserved. Consequently, every strategy that assigns specific values to them is valid, but some are less practical than others. The splinekit library deals with paddings of low complexity; in particular, we focus on some for which the overall organization of \(c,\) \(y,\) and \(f\) can be made consistent and solves the uniqueness issue. The seven paddings being considered are
Periodic
Narrow Mirror
Wide Mirror
Anti-Mirror
Nega-Periodic
Nega-Narrow Mirror
Nega-Wide Mirror
Except for the anti-mirror padding, all proposed forms honor some sort of periodicity over \(c,\) albeit the length of a period may differ from the number \(K\) of observed samples. The periodicity of \(c\) then implies that \(y\) and \(f\) are likewise periodic; moreover, the trivial choice \(\lambda_{n,m}^{-}=0\) and \(\lambda_{n,m}^{+}=0\) is the only one that results in the alternative versions \(c'\) of \(c\) being periodic, too, which ensures uniqueness and makes the interpolation problem well-defined. In the case of the anti-mirror padding, a similar reasoning also leads to the uniqueness of the solution.
Side Note The periodicity of the coefficients implies the periodicity of the samples, but the periodicity of the samples does not necessarily imply that the the coefficients are periodic. Indeed, the constants \(\lambda_{n,m}^{-},\lambda_{n,m}^{+}\) can still be chosen arbitrarily to parameterize a family of coefficients that satisfy the interpolation constraint, even for periodic samples. In general, only one member of this family will have the same periodicity as that of the samples.
We give now a piece of code that illustrates visually the effect of the various paddings on random splines of a specified degree.
Jupyter Lab notebook
Periodic Padding
An easy, general-purpose padding approach is to engineer the sequence \(c\) of coefficients to be \(K\)-periodic. This implies that the sequence \(y\) of samples has to be the straighforward \(K\)-periodized version of the vector \({\mathbf{y}}\in{\mathbb{R}}^{K}.\) Ultimately, this also implies that the function \(f\) is itself \(K\)-periodic. In summary, under a periodic padding, the relations being satisfied for any \(k\in{\mathbb{Z}}\) and any \(x\in{\mathbb{R}}\) are
Algorithmic Considerations
In the context of a straightforward periodic padding, there are three major algorithmic approaches to the solution of the interpolation constraint \(y[k]=f(k)\) for \(k\in[0\ldots K-1].\)
The linear-algebra approach first establishes an explicit system of \(K\) linear equations. The \(k\)-th equation of the system would be \(y[k]=\sum_{k'=0}^{K-1}\,\left(\sum_{p\in{\mathbb{Z}}}\,\varphi(k-p\,K-k')\right)\,c[k'].\) Tools of linear algebra would then be deployed to solve the system in terms of the unknown vector \({\mathbf{c}}=\left(c[k]\right)_{k=0}^{K-1}\) that represents one period of the periodic sequence \(c,\) taken such that the first vector component \(c_{1}\) is the sequence element at the origin \(c[0].\) The overall computational cost is \({\mathcal{O}}(K^{3})\) when general solvers are used, and the cost would be \({\mathcal{O}}(K^{2})\) for Toeplitz systems. For periodic paddings, the system is circulant and the overall computational cost reduces to \({\mathcal{O}}(K\,\log K)\) with Fourier-based techniques to solve linear-algebra inversion problems.
The discrete-Fourier approach is best described concisely in matrix notations. Let \({\mathbf{F}}\in{\mathbb{C}}^{K\times K}\) be the discrete-Fourier transform, with the \(\nu\)-th row and \(k\)-th column entry given by \({\mathrm{e}}^{-{\mathrm{j}}\,\left(\nu-1\right)\,\frac{2\,\pi}{K}\,\left(k-1\right)}.\) Let \({\mathbf{\phi}}=(\sum_{p\in{\mathbb{Z}}}\,\varphi(p\,K+k))_{k=0}^{K-1}\) be the data-independent vector that contains the samples (at the integers) of the periodized version of the synthesis function \(\varphi.\) Then, one has that \({\mathbf{c}}={\mathbf{F}}^{-1}\,\left(\left({\mathbf{F}}\,{\mathbf{y}}\right)\oslash\left({\mathbf{F}}\,{\mathbf{\phi}}\right)\right),\) where \(\oslash\) is an element-wise division. In practice, the Fourier transformation and its inverse are implemented via the fast-Fourier algorithm, in which case the overall computational cost is \({\mathcal{O}}(K\,\log K).\)
The recursive-filtering approach is the one followed in the
splinekitlibrary. It requires that \(\varphi\) has a finite support, is even-symmetric, and that the poles of the reciprocal of the \(z\)-transform of its samples at the integers are real numbers. These properties are all satisfied by the polynomial B-splines of nonnegative integer degree \(n\in{\mathbb{N}}+2.\) The overall computational cost is now \({\mathcal{O}}(K\,\left\lfloor n/2\right\rfloor).\)
Recursive Filtering
Start the algorithm by letting \({\mathbf{c}}\leftarrow{\mathbf{y}}.\) Then, iteratively for every one of the poles \(z_{n,m}\in(-1,0)\) indexed by \(m\in[1\ldots\left\lfloor n/2\right\rfloor]\) and associated to the degree \(n,\) apply the in-place recursive updates
In practice, some computations can be spared if the sums that appear in the recursive-update equations are truncated at that index \(k\) where the term \(z_{n,m}^{k}\) becomes negligible.
Experimental Performance
The following experiment establishes some simple statistics over the gain in speed achieved by the recursive approach over the fast-Fourier-based one. For each length \(K\) and for each spline degree, we synthesize \(50\) random vectors \({\mathbf{y}}\in{\mathbb{Z}}^{K}\) and let both the Fourier approach and the recursive approach determine the spline coefficients on the same data. We measure the duration of those computations and report by how many times the recursive approach is faster relatively to the Fourier approach. For instance, the number \(2.0\) would mean that the recursive approach would be \(200\%\) times faster or, equivalently, that it runs twice as fast.
Danger
The link below allows you to inspect the notebook. Unfortunately, running it from the browser is meaningless: the timings of the installation-free version are not representative because the Python kernel is WebAssembly-based and does not run natively.
If you want to test for realistic timings on your own computer, then you will have to first install in full the
splinekitlibrary. Only after that will you be able to launch the notebook either as a regular, full-fledged Jupyter Lab or as a module executed by the native Python kernel.The timings reported in the Results Section correspond to those of the native execution.
Jupyter Lab notebook
Hint
The notebook is available for download in compressed form from
here. Decompression is achieved from the terminal with gunzip padding_speed.ipynb.gz.
Results
In the first table, the number \(K\) of observed data are powers of two, a situation that is very much to the advantage of the discrete-Fourier methods. On a desktop computer of 2021, typical resulting tables are as follows.
Degree 2 3 4 5 6 7 8 9
Period
K = 1 0.5 0.5 0.5 0.5 0.4 0.5 0.5 0.5
K = 2 10.2 8.0 9.9 7.0 11.0 8.8 12.7 7.7
K = 4 8.3 8.2 10.1 7.6 11.6 11.6 12.6 10.9
K = 8 8.9 8.8 10.2 8.1 12.5 10.3 12.7 10.5
K = 16 9.8 9.7 7.9 10.0 12.0 12.2 10.7 11.7
K = 32 11.5 11.4 9.8 12.4 13.0 9.9 11.4 13.7
K = 64 15.2 14.8 12.6 14.7 14.5 14.5 12.4 14.5
K = 128 20.4 20.4 15.6 18.2 16.6 16.5 13.7 15.4
K = 256 27.8 23.8 23.1 22.6 19.4 17.4 17.0 16.7
K = 512 40.3 34.9 29.4 28.9 23.3 21.4 19.3 19.0
K = 1024 52.6 49.6 37.0 35.6 26.9 24.9 21.7 21.1
K = 2048 75.6 69.7 44.8 42.9 31.3 29.2 23.3 23.1
K = 4096 96.5 96.6 54.1 53.1 36.8 33.4 26.5 26.5
K = 8192 101.6 99.5 53.6 53.0 36.1 31.8 25.1 25.0
K = 16384 107.4 105.9 55.7 55.2 38.1 32.4 25.8 25.4
In the second table, we repeat the experiment, with the difference that the \(50\) lengths are now chosen uniformly at random in some range. On a desktop computer of 2021, typical resulting tables are as follows.
Degree 2 3 4 5 6 7 8 9
Period Range
1 ≤ K < 2 0.6 0.5 0.5 0.5 0.5 0.5 0.5 0.5
2 ≤ K < 4 6.1 8.3 9.5 9.1 9.4 11.0 12.5 11.5
4 ≤ K < 8 8.6 8.2 9.8 7.6 11.4 11.3 11.5 9.8
8 ≤ K < 16 9.5 9.1 10.4 8.3 9.2 11.7 12.6 10.2
16 ≤ K < 32 10.7 10.5 11.5 9.4 12.5 12.6 13.5 11.2
32 ≤ K < 64 13.4 13.4 13.7 11.3 13.9 14.3 14.2 14.4
64 ≤ K < 128 19.0 18.3 14.3 17.3 16.2 16.2 13.6 15.5
128 ≤ K < 256 25.0 24.8 18.2 21.1 18.7 18.5 15.3 16.7
256 ≤ K < 512 36.3 35.8 25.1 27.6 22.7 22.3 18.2 19.1
512 ≤ K < 1024 57.3 55.9 36.1 37.9 29.5 29.2 23.1 23.8
1024 ≤ K < 2048 80.2 69.8 47.3 45.9 33.7 31.1 26.0 25.7
2048 ≤ K < 4096 106.5 102.1 62.6 61.3 43.3 38.5 31.3 30.9
4096 ≤ K < 8192 121.8 106.9 59.6 59.3 41.0 35.8 28.6 28.3
8192 ≤ K < 16384 123.4 122.7 60.9 63.5 43.5 37.4 29.4 29.5
16384 ≤ K < 32768 133.3 130.2 68.0 67.0 45.5 39.5 31.1 32.1
Discussion
The conclusion of the experiments is unequivocal: The recursive approach is substantially faster than the fast-Fourier approach, at all lengths (except \(K=1\)), ranges of lengths, and all degrees being investigated.
Narrow-Mirror Padding
Under a narrow-mirror padding, one assumes \(\forall k\in{\mathbb{Z}}\) that the spline coefficients satisfy
Then, it holds for any \(k\in{\mathbb{Z}}\) and any \(x\in{\mathbb{R}}\) that
Recursive Filtering
Start the algorithm by letting \({\mathbf{c}}\leftarrow{\mathbf{y}}.\) Then, iteratively for every one of the poles \(z_{n,m}\in(-1,0)\) indexed by \(m\in[1\ldots\left\lfloor n/2\right\rfloor]\) and associated to the degree \(n,\) apply the in-place recursive updates
Wide-Mirror Padding
Under a wide-mirror padding, one assumes \(\forall k\in{\mathbb{Z}}\) that the spline coefficients satisfy
Then, it holds for any \(k\in{\mathbb{Z}}\) and any \(x\in{\mathbb{R}}\) that
Recursive Filtering
Start the algorithm by letting \({\mathbf{c}}\leftarrow{\mathbf{y}}.\) Then, iteratively for every one of the poles \(z_{n,m}\in(-1,0)\) indexed by \(m\in[1\ldots\left\lfloor n/2\right\rfloor]\) and associated to the degree \(n,\) apply the in-place recursive updates
Anti-Mirror Padding
Under an anti-mirror padding, one assumes \(\forall k\in{\mathbb{Z}}\) that the spline coefficients satisfy
Then, it holds for any \(k\in{\mathbb{Z}}\) and any \(x\in{\mathbb{R}}\) that
Recursive Filtering
Start the algorithm by letting \({\mathbf{c}}\leftarrow{\mathbf{y}}.\) Then, iteratively for every one of the poles \(z_{n,m}\in(-1,0)\) indexed by \(m\in[1\ldots\left\lfloor n/2\right\rfloor]\) and associated to the degree \(n,\) apply the in-place recursive updates
Nega-Periodic Padding
Under a nega-periodic padding, one assumes \(\forall k\in{\mathbb{Z}}\) that the spline coefficients satisfy
Then, it holds for any \(k\in{\mathbb{Z}}\) and any \(x\in{\mathbb{R}}\) that
Recursive Filtering
Start the algorithm by letting \({\mathbf{c}}\leftarrow{\mathbf{y}}.\) Then, iteratively for every one of the poles \(z_{n,m}\in(-1,0)\) indexed by \(m\in[1\ldots\left\lfloor n/2\right\rfloor]\) and associated to the degree \(n,\) apply the in-place recursive updates
Nega-Narrow-Mirror Padding
Under a nega-narrow-mirror padding, one assumes \(\forall k\in{\mathbb{Z}}\) that the spline coefficients satisfy
Then, it holds for any \(k\in{\mathbb{Z}}\) and any \(x\in{\mathbb{R}}\) that
Recursive Filtering
Start the algorithm by letting \({\mathbf{c}}\leftarrow{\mathbf{y}}.\) Then, iteratively for every one of the poles \(z_{n,m}\in(-1,0)\) indexed by \(m\in[1\ldots\left\lfloor n/2\right\rfloor]\) and associated to the degree \(n,\) apply the in-place recursive updates
Nega-Wide-Mirror Padding
Under a nega-wide-mirror padding, one assumes \(\forall k\in{\mathbb{Z}}\) that the spline coefficients satisfy
Then, it holds for any \(k\in{\mathbb{Z}}\) and any \(x\in{\mathbb{R}}\) that
Recursive Filtering
Start the algorithm by letting \({\mathbf{c}}\leftarrow{\mathbf{y}}.\) Then, iteratively for every one of the poles \(z_{n,m}\in(-1,0)\) indexed by \(m\in[1\ldots\left\lfloor n/2\right\rfloor]\) and associated to the degree \(n,\) apply the in-place recursive updates