Abstract

Contents

This page follows the manuscript closely. It is written as mathematical exposition rather than as a summary page.

1 Introduction

Polynomial interpolation is a deceptively simple global problem. Given distinct nodes \[x_0<x_1<\cdots<x_{n-1}\] and sampled values \(y_i=f(x_i)\), one asks for a function agreeing with the data at the nodes and behaving sensibly between them. Classical formulas answer this question in several equivalent ways. Lagrange interpolation writes the interpolant in a nodal basis, Newton interpolation organizes the same polynomial through divided differences, barycentric formulas reorganize Lagrange evaluation for improved numerical behavior, and orthogonal bases replace nodal formulas by coefficient expansions [1, 9, 2, 6]. Mathematically these formulas often represent the same unique polynomial. Computationally they are very different.

That difference becomes sharper in a distributed setting. If a representation is first constructed on one machine and only then broadcast to workers, parallel evaluation can be excellent even though the construction stage remains a centralized bottleneck. Conversely, a method with a slightly less classical representation may be preferable if the representation itself can be assembled from partition-level summaries. MapReduce makes this distinction unavoidable: the relevant question is not only whether evaluations can be mapped independently, but whether the object being evaluated can be constructed without gathering all of the data back to a single driver [3, 10].

The first goal of this paper is therefore diagnostic. We place a broad collection of interpolation methods into a common distributed viewpoint and keep separate four costs that are easily conflated: \[\text{construction},\qquad \text{evaluation},\qquad \text{memory},\qquad \text{communication}.\] Global formulas often expose substantial parallelism in the number of query points while retaining construction costs that grow quadratically or cubically with the number of samples. Piecewise methods avoid this difficulty by replacing one global representation with many small local ones, but the resulting continuation beyond the observed interval is necessarily local.

The second goal is constructive. We replace the global polynomial representation by a finite linear recurrence. For an order \(k\) vector \[a=(a_0,\dots,a_{k-1}),\] we seek a relation \[\label{eq:intro-recurrence} y_{i+k}\approx \sum_{j=0}^{k-1}a_jy_{i+j}.\] When \(k\) is fixed independently of \(n\), the representation has constant dimension. The recurrence can be fitted from local rows of a tall least-squares system, each partition can be reduced to a small triangular summary, and those summaries can be merged through a tree reduction. Once the coefficients are known, the recurrence itself is a fixed-dimensional state evolution.

There is also a structural reason to expect [eq:intro-recurrence] to work. Uniformly sampled exponential polynomials satisfy exact finite recurrences. In particular, this includes uniformly sampled polynomials, exponentials, trigonometric sums, and damped trigonometric sums. For an arbitrary \(C^k\) function the relation need not be exact, but the \(k\)th finite difference gives an explicit order-\(k\) recurrence whose one-step residual is bounded by \(h^k\|f^{(k)}\|_\infty\). Thus the recurrence model is exact on a large algebraically structured class and locally high-order accurate on smooth data.

The remaining issue is stability. A small one-step residual need not remain small after many recurrence steps. Writing the recurrence through its companion matrix separates these two questions cleanly: the residual controls the local defect, while powers of the companion matrix control how the defect is amplified during continuation. This gives a direct criterion for stable long-range propagation and explains why recurrence order alone is not a sufficient model-selection criterion.

1.1 Main results

The paper develops five results around this viewpoint.

First, we make the construction/evaluation distinction explicit for the classical methods used in the benchmark. Their formulas may define the same polynomial, but they incur different driver-side preprocessing and per-query costs.

Second, we prove exact finite-recurrence reconstruction for sampled exponential polynomials. This gives a concrete class on which recurrence continuation is not an approximation of the sampled sequence but an exact representation of it.

Third, for uniformly sampled \(f\in C^k\), we prove an \(O(h^k)\) one-step residual bound and transfer it to the fitted least-squares recurrence. The result gives a deterministic approximation guarantee without assuming that the function itself belongs to an exact recurrence class.

Fourth, we derive a companion-matrix propagation identity. This identity isolates the two quantities governing continuation error: the fitted residuals and the growth of powers of the companion matrix.

Finally, we give a distributed weighted-QR fitting algorithm. For \(p\) workers and fixed recurrence order \(k\), its total work is \(O(nk^2+pk^3)\), each partition is represented by an \(O(k^2)\) triangular summary, and the recurrence continuation admits block parallelism through matrix powers. In the fixed-\(k\) regime, the dependence on the number of samples is linear.

The supplied experimental records are consistent with the motivation for this theory. Piecewise and recurrence-based methods remain feasible on the largest recorded problems, whereas several global constructions deteriorate rapidly. The recurrence method also preserves the strongest extrapolation scores in the supplied benchmark tables. We treat those numerical observations as empirical evidence rather than as a substitute for the stability conditions proved below.

[top]

2 Classical interpolation in a distributed model

We begin with the computational distinction that motivates the recurrence construction.

The distinction is particularly important for formulas with expensive global preprocessing. For example, Newton divided differences permit \(O(n)\) evaluation once the divided-difference table has been formed, but constructing those coefficients is still quadratic in \(n\). Orthogonal-basis interpolation can make evaluation clean while moving the main cost into a dense coefficient solve. Aitken–Neville avoids a separate coefficient table but pays quadratic work at each query. Piecewise linear interpolation, in contrast, avoids a global polynomial representation entirely.

The proposition is not intended as a universal lower bound for every possible implementation of these mathematical methods. It records the computational structure of the implementations being compared. That distinction matters: fast multipoint evaluation and specialized structured algorithms can alter asymptotic costs, but the benchmark question here is how the supplied MapReduce implementations behave when placed under the same execution model.

[top]

3 Finite recurrences

We now turn from global polynomial representations to constant-order relations among neighboring samples.

The principal exact class is the class of exponential polynomials.

This theorem explains why the recurrence model is exact for much more than a single periodic signal. Polynomial trends, sums of oscillatory modes, and exponentially damped oscillations all lie in the same finite-dimensional shift-invariant framework.

3.1 Bidirectional reconstruction

The implementation uses two recurrences rather than forcing every interior value to be generated from one boundary.

[top]

4 Approximation and stability

Exact recurrence structure is not available for an arbitrary smooth function. The useful fact is that a uniform grid still produces an explicit approximate recurrence.

The implemented model does not force these binomial coefficients. It fits the recurrence from the data by weighted least squares.

4.1 Companion-matrix propagation

A local residual bound is only half of a continuation argument. The second half is the stability of repeated recurrence application.

Let \[C(a) = \begin{bmatrix} 0&1&0&\cdots&0\\ 0&0&1&\cdots&0\\ \vdots& & &\ddots&\vdots\\ 0&0&0&\cdots&1\\ a_0&a_1&a_2&\cdots&a_{k-1} \end{bmatrix}\] be the companion matrix, and let \[s_i=(y_i,\dots,y_{i+k-1})^T.\] Then the recurrence defect gives \[s_{i+1}=C(a)s_i+e_k r_i(a),\] where \(e_k=(0,\dots,0,1)^T\).

This is the central stability separation. A recurrence can fit the observed samples extremely well and still be a poor extrapolator if its companion dynamics are unstable. Conversely, a modest local residual can remain controlled when the state evolution is stable. Any practical recurrence-order selection rule should therefore account for both fit quality and the spectrum of the companion matrix.

[top]

5 Distributed recurrence fitting

The recurrence design matrix is tall and narrow. This is exactly the geometry in which local QR reduction is useful.

For each admissible recurrence row define the augmented weighted row \[Z_i=(w_iR_i,\;w_it_i)\in\mathbb{R}^{k+1},\] and let \(Z\) be the matrix obtained by stacking these rows.

Only the short sample histories needed to form recurrence rows that cross partition boundaries must be exchanged before the local reductions. Once a partition has formed its triangular factor, the full local block is no longer needed by later reduction stages.

5.1 Parallel recurrence propagation

Ordinary recurrence continuation is sequential inside one block, but the starting state of a distant block can be obtained directly from a power of the companion matrix.

The construction and continuation now share the same basic feature: every global stage operates on fixed-dimensional objects when \(k\) is fixed.

[top]

6 Reconstruction algorithm

The full method uses one recurrence in each direction. Inside the observed interval the two modeled sequences are blended, while outside the interval the nearest directional recurrence is continued.

The approximation guarantees in [thm:smooth-recurrence,cor:ls-residual] apply directly only to uniform sampling. The supplied benchmark drafts also contain experiments on random and Chebyshev node sets. Those tests are useful as stress tests of the implementation, but they should not be read as being covered by [thm:smooth-recurrence,cor:ls-residual] without an additional nonuniform-grid analysis.

[top]

7 Recorded experiments

The supplied drafts contain benchmark records for classical methods, two piecewise methods, and the finite-recurrence method. We reorganize those records here around three questions: how runtime grows with problem size, how the reported interior score changes with node geometry, and what happens under extrapolation.

The test-function descriptions in the newest source draft use \[f_{\mathrm{easy}}(x)=\sin x,\qquad f_{\mathrm{medium}}(x)=e^x,\qquad f_{\mathrm{hard}}(x)=\frac{1}{1+25x^2}.\] The recorded node families are equispaced, random, and Chebyshev nodes.

7.1 Size-up behavior

The small-\(n\) records already show the separation between quadratic-style global methods and the scalable methods. Selected entries are shown in 1.

The purpose of these measurements is not to infer an asymptotic exponent from five noisy wall-clock values. The useful observation is more basic: Aitken and direct Lagrange deteriorate rapidly over the recorded range, whereas the piecewise and recurrence methods remain essentially flat at these sizes.

The large-scale records continue the comparison for the three methods that remained feasible in the source experiments.

All three methods reach the million-point case in the supplied records. The recurrence timing is comparable to the two piecewise baselines and is the smallest of the three at \(10^6\) points in this particular run. Without the missing hardware metadata, we do not attach significance to the small constant-factor differences; the robust conclusion is the common large-scale feasibility.

7.2 Interior score and node geometry

The recurrence method behaves very differently on the three node families in the supplied score table.

The equispaced results are the most directly aligned with the theory in 9. Random and Chebyshev nodes change the physical spacing represented by one recurrence step, so the degradation there is not surprising: the fitted relation is being asked to model a sequence indexed by sample number even though equal increments in that index no longer correspond to equal increments in \(x\).

The same source table records near-\(100\) scores for the two piecewise baselines across all three node geometries. This is consistent with the role of piecewise methods as local interpolation procedures. It also highlights an important limitation of the present recurrence formulation: its strongest theory and strongest recorded interior behavior are both tied to uniform sampling.

7.3 Recorded extrapolation behavior

The supplied extrapolation table shows the clearest empirical separation. For the easy and medium tests, the recurrence method records a score of \(100\) at every listed shift from \(1\) through \(5\), while the two piecewise methods decrease steadily. The hard test is less favorable for every method, but the recurrence scores remain slightly above the corresponding piecewise values in the supplied record.

These observations should be read together with 13. Recurrence continuation is not automatically stable merely because it avoids a high-degree polynomial. Its stability depends on the fitted companion dynamics. The benchmark suggests that the fitted recurrences in the easy and medium experiments had favorable continuation behavior, while the theorem explains the mechanism that must be checked if that behavior is to be guaranteed on a new data set.

[top]

8 Discussion

The distributed comparison reveals three distinct notions of scalability.

The first is query scalability. Once a representation has been constructed, many interpolation formulas allow independent evaluations at different query points. MapReduce handles this case naturally. This is the easiest form of parallelism and is shared by most methods in the comparison.

The second is construction scalability. Here the classical global methods separate. Pairwise weights, divided-difference tables, and dense basis solves must be formed before the independent evaluations can begin. Parallel query evaluation cannot remove a construction step that remains on the critical path. The recurrence fit differs because every sample contributes only a row of a tall, fixed-width least-squares problem, and those rows can be compressed locally by QR.

The third is continuation scalability. A recurrence appears sequential because each new value depends on the previous state, but the companion-matrix formulation makes block parallelism available. A worker can compute the state at the beginning of its block from a matrix power and then generate that block locally. Thus the recurrence model admits distributed structure both before and after the coefficient vector has been learned.

There is, however, an equally important limitation. The clean approximation result is a uniform-grid theorem. A recurrence indexed only by sample number cannot distinguish between a short physical step and a long one when the nodes are irregular. The random- and Chebyshev-node benchmark scores reflect exactly this sensitivity. Extending the theory to nonuniform grids would require either variable-coefficient recurrences, a reparameterization by arc length or cumulative spacing, or a state model in which the node increment enters explicitly.

A second limitation concerns extrapolation. The recurrence framework gives a mechanism for stable extrapolation, not an unconditional guarantee. The fitted residual and companion spectrum must both be controlled. In practical implementations this suggests regularizing or rejecting recurrence fits with unstable companion roots, even if their in-sample residual is small.

Finally, the benchmark metadata need to be completed before the numerical section is publication-ready. The supplied drafts do not record the full execution environment and do not define the percentage-valued accuracy score used in the tables. These omissions do not affect the mathematical results, but they prevent reproducible interpretation of the wall-clock constants and the percentage scores. The next experimental revision should therefore record hardware, Spark configuration, partition counts, recurrence order, train/evaluation intervals, random seeds, and the exact error-to-score transformation.

[top]

9 Conclusion

The main computational difficulty in large-scale interpolation is not always the evaluation of the interpolant. It is often the construction of the representation that must be evaluated. Classical formulas can distribute their queries while retaining global preprocessing, dense memory, or coefficient-construction costs. Piecewise methods avoid those costs by localizing the approximation. Finite recurrences offer a third route: learn a fixed-dimensional global state rule from distributed summaries.

The recurrence viewpoint is exact for uniformly sampled exponential polynomials and locally high-order accurate for general smooth functions on uniform grids. Its continuation error admits an explicit companion-matrix representation, which separates local fitting quality from dynamical stability. Its weighted least-squares fit can be implemented by local QR factorizations and tree reduction, so for fixed recurrence order the total fitting work is linear in the number of samples and communication is independent of the global data size except through the number of partitions.

The recorded benchmarks support the computational premise. The recurrence method remains feasible at the million-point scale in the supplied experiments and records strong continuation scores relative to the piecewise baselines. At the same time, the node-geometry experiments expose the principal limitation of the present theory: the recurrence is naturally a uniform-grid model. This makes nonuniform recurrence models, spectrum-aware order selection, and fully reproducible distributed benchmarks the most immediate directions for further work.

[top]

References