Background: I'm hoping to find (or write) an algorithm to piecewise linear-interpolate large sets of unevenly sampled functions (10s of thousands of arrays of a thousand or so $x$ and $y$ pairs, where each array $x_i$ are not the same). Iterating over each $y_i(x_i)$ and creating interpolants is expensive (see this stackoverflow question for more info). Because polynomial interpolation can be formulated as a (Vandermond) matrix problem it can be sped up significantly by attacking it as a single 3D matrix.
Question: Is there any way to formulate piecewise interpolation as a similar matrix problem, so that it can be sped up? Alternatively, is there a different set of basis-functions (besides polynomials) which would be guaranteed to be monotonic between each pair of points, which could be linearized?