Spline Library in DRIP


SplineLibrary provides the functionality for building, calibrating, and evaluating different kinds of splines for use in latent state representation. It implements the functionality behind spline design, spline constructions, customization, calibration, and evaluation of a wide variety of spline types and basis functions.


SplineLibrary achieves its design goals by implementing its functionality over several packages the perform the following:

·        Univariate Function Package: The univariate function package implements the individual univariate functions, their convolutions, and reflections.

·        Univariate Calculus Package: The univariate calculus package implements univariate difference based arbitrary order derivative, implements differential control settings, implements several integrand routines, and multivariate Wengert Jacobian.

·        Spline Parameters Package: The spline parameters package implements the segment and stretch level construction, design, penalty, and shape control parameters.

·        Splinbe Basis Function Set Package: The spline basis function set package implements the basis set, parameters for the different basis functions, parameters for basis set construction, and parameters for B Spline sequence construction.

·        Spline Segment Package: The spline segment package implements the segment’s inelastic state, the segment basis evaluator, the segment flexure penalizer, computes the segment monotonicity behavior, and implements the segment’s complete constitutive state.

·        Spline Stretch Package: The spline stretch package provides single segment and multi segment interfaces, builders, and implementations, along with custom boundary settings.

·        Spline Grid/Span Package: The spline grid/span package provides the multi-stretch spanning functionality. It specifies the span interface, and provides implementations of the overlapping and the non-overlapping span instances. It also implements the transition splines with custom transition zones.

·        Spline PCHIP Package: The spline PCHIP package implements most variants of the local piece-wise cubic Hermite interpolating polynomial smoothing functionality. It provides a number of tweaks for smoothing customization, as well as providing enhanced implementations of Akima, Preuss, and Hagan-West smoothing interpolators.

·        Spline B Spline Package: The spline B Spline package implements the raw and the processed basis B Spline hat functions. It provides the standard implementations for the monic and the multic B Spline Segments. It also exports functionality to generate higher order B Spline Sequences.

·        Tension Spline Package: The tension spline package implements closed form family of cubic tension splines laid out in the basic framework outlined in Koch and Lyche (1989), Koch and Lyche (1993), and Kvasov (2000).


Download Source and the binary from the Downloads area.






Framework Symbology and Terminology


1.      Predictor Ordinates: The segment independent/input values.

2.      Response Values: The segment dependent/output values.

3.      C0, C1, and C2 Continuity: C0 refers to base function continuity. C1 refers to the continuity in the first derivative, and C2 refers to continuity in the second.

4.      Local Piece-wise Parameterized Splines: Here the space formulation is in the local variate space that spans 0 to 1 within the given segment – this is also referred to as piece-wise parameterization.

5.      Bias: This is the first term in the Spline Objective Function – essentially measures the exactness of fit.

6.      Variance: This refers to the second and the subsequent terms in the Spline Objective Function – essentially measures the curvature/roughness.





1.      Definition: “Spline is a sufficiently smooth polynomial function that is piecewise-defined, and possesses a high degree of smoothness at the places where the polynomial pieces connect (which are known as knots). [Spline (Wiki), Judd (1998), Chen (2009)]

2.      Drivers:

a.       Lower degree, gets rid of oscillation associated with the higher degrees [Runge’s phenomenon (Wiki)]

b.      Easy, accurate higher degree smoothness specification

3.      Basic Spline: Covered in [Spline (Wiki), Bartels, Beatty, and Barsky (1987), Judd (1998), De Boor (2001), Fan and Yao (2005), Chen (2009), Katz (2011)].

4.      History: Schoenberg (1946), Ferguson (1964), Epperson (1998).



Document Layout


1.      Overview

2.      Calibration Framework SKU

3.      Spline Builder SKU

4.      B Splines

5.      Polynomial Spline Stretches

6.      Local Spline Stretches

7.      Spline Calibration

8.      Spline Jacobian

9.      Shape Preserving Splines

10.  Koch-Lyche-Kvasov Tension Splines

11.  Multi-dimensional Splines

12.  References

13.  Figures

14.  Spline Library Software Components

Calibration Framework






1.      Definition: Calibration is the process of inferring the latent state elastic properties from the specified inputs.

2.      Classes of static fields: Elastic and inelastic

3.      Calibrator Creation: On creation, objects acquire specific values for the constitutive inelastic fields. Volatile Latent State elastic fields may as yet be undefined.

4.      Calibration is Inference: Since calibrated parameters are used for eventual prediction, calibration is essentially inference. Bayesian classification (an alternate, generalized calibration exercise) is inference too.

5.      Calibration and entity-variate focus:

6.      Latent State Construction off of hard/soft signals: Hard Signals are typically the truthness signals. Typically reduce to one calibration parameter per hard observation, and they include the following:

·        Actual observations => Weight independent true truthness signals

·        Weights => Potentially indicative of the truthness hard signal strength

Soft signals are essentially signals extracted from inference schemes. Again, typically reduce to one calibration parameter per soft inference unit, and they include the following:

·        Smoothness signals => Continuity, first, second, and higher-order derivatives match – one parameter per match.

·        Bayesian update metrics => Inferred using Bayesian methodologies such as maximum likelihood estimates, variance minimization, and error minimization techniques.

7.      Directionality Bias: Directionality “bias” is inherent in calibration (e.g., left to right, ordered sequence set, etc:) – this simplifies the solution space significantly, as it avoids simultaneous non-linearity. Therefore, the same directional bias also exists in the calibration nodal sequence.

8.      Truthness/Smoothness vs. Information Propagation: In segment-by-segment calibration challenges associated with inferring a composite Latent State, if the inference is based purely off of the truthness measurements, the information directionality/propagation/flow is irrelevant. Notions such as  are important primarily owing to the smoothness axioms.

9.      Head Node Calibration: Calibration of the head node is typically inherently different from the other nodes, because the inputs needed/used by it could be different. The other nodes use continuity/smoothness parameters, which the head node does not.

10.  Parameter Space Explosion: Generally not a problem as long as it is segment-localized (in linear systems parlance, as long the Latent State matrix is tri-diagonal, or close to it), i.e., local information discovery does not affect far away nodes/segments.

·        Also maybe able to use optimization techniques to trim them.

11.  Live Calibrated Parameter Updating: Use automatic differentiation to:

·        Estimate parametric Jacobians (or sub-coefficient micro-Jacobians) to the observed product measures.

·        Re-adjust the shifts using the hard-signal strength.

·        Update the parameters from the calculated shifts.

·        Re-construct the curve periodically (for increments, as opposed to a full re-build).

12.  Spline Segment Calibrator: For a given segment, its calibration depends only on the segment local value set – the other inputs come from the prior segments (except in the case of “left-most” segment, whose full set of inputs will have to be extraneously specified).

13.  Block Diagonalization vs. Segment/Segment Calibration: Since segment-by-segment span calibration occurs through  transmission, it will effectively be block diagonal if it is a linear system, and hence computationally efficient (this also allows for explicitly setting segment level controls). In other words, although the local matrix is dense, the span-level matrix is still sparse,

14.  Calibration Perspective of Supervised vs. Unsupervised:

·        Supervised – Alternate View => Since supervised learning depends upon a training step, post supervised activity may be viewed as a post-inference systems, where the inference/calibration has already occurred during the training step (and parsimonization into the parameter step). So, post-training, the supervised machine is only used for prediction.

·        Unsupervised – Alternate View => Here, inference and prediction happen simultaneously.

·        Hybrid Supervised/Unsupervised => Clearly no TRUE unsupervised are possible (as there are prior views on linearity, Markov nature, error process, basis spline representation choice, etc;). So unsupervised systems are typically mostly hybrid systems.



Latent State Specification


1.      Latent State Specification: The latent-state here refers to the state whose dependent response values we wish to calibrate/infer as a function of the predictor ordinate. For e.g., in fixed income finance, the Prevailing Interest Rate, the survival, the forward  rate, FX spot/forward would each constitute a potentially separate latent state that needs to be inferred.

2.      Latent State Quantification Metric: A given latent state may be described by one/more alternate, mutually exclusive quantification metric. Again, the discounting latent state may be quantified using a discount factor, zero rate etc:

3.      Manifest Measure: The Latent State may be inferred using a variety of external experimental measurements each of which produces a convolved signal of the latent state. Each such signal is referred to as a manifest metric. Again, the discounting latent state may be inferred from the cash instrument manifest metric, the swap rate manifest metric etc:

Spline Builder Setup



Design Objectives behind Interpolating Splines


1.      Symbols and Definition:

2.      Monotonicity:  increases with increase in  (and vice versa).

3.      Convexity:  should also be convex wherever  is convex (and vice versa).

4.      Smoothness: Smoothness (also called shape-preserving) corresponds to the least curvature. Even C0 can be “smooth”, and so is Ck.

5.      Locality: Locality means that the dependence of  is primarily only on  and . This is advantageous to schemes that locally modify/insert the points.

6.      Approximation Order: Approximation Order indicates the smallest polynomial degree by which  departs from  as the density of  increases. More formally, it is the m in , where .

7.      Other Desired Criteria:

8.      Assessment of Monotonicity and Convexity: An individual segment can be assessed to be monotone/convex etc:, but from the data PoV, peaks, valleys, and inflection occur only at the knots. These can be assessed only at the span level.



Base Formulation


1.      Base Mathematical formulation:

·        , therefore .

·        From known nodes  and , we can draw the 2 linear equations for :



·        From known nodal derivatives , where , we can draw the following r linear equations for :

o        where

2.      Linear of Segment Coefficients to the Response Values (): In all the spline formulations, the Jacobian  is constant (i.e., independent of the response values or their nodal derivative inputs).

3.      Span Boundary Specification:

·        “Natural” Spline – Energy minimization problem – Second Derivative is Zero at either of the extreme nodes.

·        “Financial” Spline – Second Derivative at the left extreme is zero, but first derivative at the right extreme is zero.

·        Clamped Boundary Conditions: , and .

·        Not-A-Knot Boundary Conditions:  and .

4.      Discrete Segment Mesh vs. Inserted Knots: Inserting knot point is similar to discretizing the segment into multiple grids, with one key difference:

·        Discretization uses the same single spline across all the grid units of the segment.

·        Inserted knots introduce additional splines – on between each knot pair.

5.      Segment Inelastics: These are effectively the same as shape controller, i.e., the following are the shape controlling inelastic parameter set:

·        Tension

·        Number of basis

·        Continuity

·        Optimizing derivative set order






1.      Motivation: As postulated by De Boor et. al. (see De Boor (2001)), B Splines have a geometric interpolant context – thereby with the correspondingly strong CADG/curve/surface construction focus. Smoothening occurs as a natural part of this.

2.      kth Order B-Spline Interpolant: Higher order B-Splines are defined by the recurrence  where  and  if , and  otherwise.

3.      Recursive Interpolant Scheme: B Spline formulation is recursively interpolant, i.e., the order k spline is interpolant over the order  splines on nodes  and  - this formulation automatically ensures  nodal continuity.

4.      B-Spline Order Relationships: Assuming no coincident knots, the following statements are all EQUIVALENT/TRUE:

5.      Expository Formulation:

6.      Spline Coefficient Partition of Unity: Using the earlier formulation , it is easy to show that . This simply follows from the recursive nodal interpolation property.

7.      Smoothness Multiplicity Order Linker: # smoothness conditions at knot + the multiplicity at the knot = B-Spline Order.

8.      Starting Node de-biasing: Left node is always weighted by  in the interpolation scheme, but the left node asymmetry is maintained because the denominator in , i.e.,  increases in length.

9.      Other Single B-Spline Properties:

10.  Formulation off of Starting Node and Starting Order: Given the starting node  and the starting order , the contribution to the node  (i.e., m nodes after the start) and the order  (i.e., n nodes after start) can be “series”ed as

11.  Cardinal B-Spline Knot Sequence: Knot sequence  => Uniformly spaced knots, simplifying the interpolant/recursive analysis significantly - .



·                    Cardinal B-Spline Order 3:



12.              Non-coinciding B Spline Segment Relations:



13.              Bernstein B-Spline Knot Sequence: Knot sequence  -  occurs  times, and  occurs  times.

·                     for .

·                     has  derivatives at , and  derivatives at  - this is also referred to as  smoothness conditions at , and  smoothness conditions at .

14.              B-Spline vs. Spline: B-Spline is just a single basis polynomial that is valid across a set of knots. “Spline” is a linear combination of such basis B Splines – i.e., the set of all the ’s.

15.              Spline Definition:  where . ’s are the coefficients – or nodal points  - that can be interpolated.



B Spline Derivatives


1.      B-Spline Derivative Formulation:


2.      B Spline Order 3 Nodal Slopes: The slopes match across the left and the right segment, as shown below, thereby making   continuous.


Left Slope

Right Slope





3.      B Spline Continuity Condition: From the B Spline derivative formulation it is clear that if both  and  are  continuous, then  will be  continuous. Given that B Spline order 3 is  continuous, by induction,  is  continuous.

Polynomial Spline Basis Function



Plain Polynomial Spline Basis Function


1.      Most direct polynomial spline fit is the Lagrange polynomial that passes through the sequence of given points.

2.      Young (1971) was one of first to apply shape-preserving polynomial using diminished Lagrange Polynomials (Lagrange Polynomials (Wiki)), showing that co-monotone interpolant with an upper bound on the polynomial degree exists (Raymon (1981)).

3.      Knot Insertion and Control Techniques: Careful knot insertion can produce:



Bernstein Polynomial Basis Function


1.      Bernstein Polynomial of degree n, and term :  where .

2.      Derivative of the Bernstein Polynomial: .

3.      Bernstein Recurrence: .

4.      Reduction of B-Splines to Bernstein’s Polynomial: From the recurrence relation, it is clear that this is exactly the same recurrence as that for B-splines, except that it happens over repeating knots at  and .

Local Spline Stretches



Local Interpolating/Smoothing Spline Stretches


1.                  Hermite Cubic Splines: The “local information” here takes the form of user specified left/right slopes + calibration points.

a.                   2 User Specified local slopes + 2 points => 4 sets of equations. Solve for the coefficients.

b.                  C1 continuity is maintained, and C2 continuity is not.

c.                   Segment control is completely local. Both the head and non-head calibration are identical/analogous for this reason.

2.                  C1 Hermite Formal Definition: For C1, the Hermite polynomial of degree  is given as , where  and  are expressed in terms of the jth Lagrange coefficient of degree ,  as



3.                  Catmull-Rom Cubic Splines (Catmull and Rom (1974)): Instead of explicitly specifying the left/right segment slopes, they are inferred from the “averages” of the prior and the subsequent points, i.e., , and . Here refers to the slope vector, and  to the point vector.

a.                   Again, C1 continuity is maintained, and C2 continuity is not.

b.                  Segment control is not completely local, but still local enough – it only depends on the neighborhood of 3 points.

4.                  Cardinal Cubic Splines: This is a generalization of the Catmull-Rom spline with a tightener coefficient , i.e., , and .  corresponds to tightening, and  corresponds to loosening.

a.                   Again, C1 continuity is maintained, and C2 continuity is not.

b.                  Segment control is “local” in the Catmull-Rom sense - it only depends on the neighborhood of 3 points.

5.                  Catmull-Rom/Cardinal Splines as Interpolation Splines: As interpolating splines, both Catmull-Rom and Cardinal are primarily useful in heuristic knot-insertions – Catmull-Rom as linear in the gaps, and Cardinal as tense linear gap knots.

·        The local knot point insertion may be generalized as follows: The targeted knot insertions follow the formulation paradigm , where  is the set of the neighborhood points. Similar formulation (with potentially different function forms, of course) may be used for each of the  derivatives. Catmull-Rom and cardinal use 1D, strictly neighboring adjacencies, as well as tense linear averaging.

6.                  Hermite-Bessel Splines: These splines use 4 basis functions per segment, therefore they are cubic polynomial, but . The first are set at each node  as the first derivative of the quadratic that passes through , , and  (the edges are handled slightly differently, as shown below). Specifically:



·         for

7.                  Hyman’s Monotone Preserving Cubic Spline:

·        Hyman (1983) applies the stringent conditions to preserve monotonicity by applying the de Boor-Schwarz criterion.

·        Define . If locally monotone (i.e., ), then set . If non-monotone (i.e., ), then set .

·        Put another way (Iwashita (2013)): For cubic polynomial splines, the first derivative should be in the range  where .

·        Adjustment for Spurious Extrema => To ensure that no spurious extrema is introduced in the interpolant,  if , and  if .

8.                  Hyman89 Extension to Hyman83: Doherty, Edelman, and Hyman (1989) relax the Hyman83 stringency posited for monotonicity preservation. Define the following:





·        If , , , , and  all have the same sign, then .

·        If , , , , and  all have the same sign, then .

·        Finally, set  if , and  otherwise.

9.                  Hyman’s Monotone Preserving Quintic Spline:

·        For quintic polynomial splines, the first derivative should be in the range  where .

·        The constraint on the second derivative is: .

·        Monotonicity Preserving Quintic Spline => Enhancement of the criterion established by de Boor and Schwartz (1977) (Hyman (1983), Doherty, Edelman, and Hyman (1989)).

o       Set  if , and  otherwise. Then

o       If , then , AND:

o       If , then .

·        Second Derivative Tests for Monotonicity Preserving Quintic Spline =>

o       Define the following constants:

§         .

§         .

§          if , and otherwise.

§          if , and otherwise.

o       Define Ranges  and  as:

§         , AND

§         .

If  and  overlap, then  should lie in their common range. If they do not overlap, i.e., , reset  as . Setting this  ensures that  and  overlap, so the other tests aspects may now continue.

10.              Harmonic Splines: Introduced by Fritsch and Butland (1984) as:

·         if ,  if . Boundary Conditions are:



·        Harmonic Spline Monotonicity Filter =>

o        if  AND  AND

o        if

o        if

o        if  AND  AND

o        if

o        if

·        Continuous Limiters => For harmonic splines, as the predictor ordinate widths become identical , setting , we get . This is the Van Leer limit (Van Leer (1974)). Huynh (1993) reviews several such limiters.

·        Shortcoming of these Limiters => Since they rely on min/max/modulus functions, by definition they are not smooth close to transition edge. This is rectified by Le Floc’h (2013), who defines a new limiter based on rational functions:  for , and  otherwise. This produces a stable  interpolator.

11.  Akima Cubic Interpolator (Akima (1970)):

12.  Kruger’s Constrained Cubic Interpolant (Kruger (2002)):

·         if .

·         if .

·         and  at the end points.

13.  Shape-Preserving Knot-based  Cubic: Ideas are taken from the awesome paper by Pruess (1993). The basic idea is to take the interval , and partition it into 3 parts by locating the two knots at , and . Evolve the criteria for the selection of and  (and, of course, their corresponding responses) such that the local spline has shape-preserving feature, and avoids being global (i.e., preserves locality).

a.       Using the above notations, the basic equations are:

·         for .

·         for .

·         for .

b.      The corresponding  maintenance solution then becomes (in terms of , , , and , whose specification will then complete the inference):

·        .

·        .

·        .

·        .

·        .

·        .

·        .

·        .

·        .

c.       Choice of  and :  and  may be generated using typical generation schemes (e.g., using the Fritsch and Butland (1984) algorithm).

d.      The Preuss Inequalities: It is specified as follows. Set

·        ; ;

·        If  and , you are done, since the chosen  and  also preserve convexity – you can go and set the second derivatives, and set  and .

e.       Mismatch in the Preuss Inequalities: If the Preuss inequalities are not met,  and  need to be modified such that  where  are obtained using the double sweep algorithm below.

f.        Preuss (1993) Double Sweep: First find  and  from the following regimes:

·        If , ; .

·        If , ; .

·        If , ; .

·        Finally the  and  initializations are set from:

·        If , ; .

·        If , ; .

·        Preuss (1993) Backward Sweep for : First set , then set  using the following:

o       If , then , and .

o       If , then , and .

o       If , then , and .

·        Preuss (1993) Setting 2nd Derivatives: .

·        Preuss (1993) Final Step – Setting  and : Setting , verify the following inequalities:



o       If these inequalities are satisfied, then you have your , , , and . Otherwise, reduce  and  till they are satisfied.



Space Curves and Loops


1.                  Space Curve Reproduction: Here is one way to construct loops that are not possible using the ordered variates, i.e., .

·                    If the ordering  is switched out in favor of the DAG , where the DAG vertices correspond to the loop trace, normal splines may be used to represent space curve loops.

2.                  Second Degree Parameterization: However, on using a second-degree parameterization of x and y such as  and , it may be possible to enforce the order . The corresponding control points are  and .

a.                   Side effect of this – is that you need to work on two pairs of splines – one each for  and .

b.                  This can offer additional customization and freedom in the design of the surface, at the expense of computing additional splines.

3.                  Closed Loops: Further, if the start/end points coincide, this corresponds to a closed loop that satisfies the C2 continuity criterion.

a.                   This also implies that no extra head/tail C1 slope specifications are required.

Spline Segment Calibration





1.                  Spline Segment Calibrator: For a given segment, its calibration depends only on the segment local value set – the other inputs come from the prior segments (except in the case of “left-most” segment, whose full set of inputs will have to be extraneously specified).

2.                  Bayesian Techniques in Spline Calibration: Frequentist and Bayesian techniques such as MLE and MAP regression ought to be possible in the calibration of the spline segment/span coefficients.

3.                  Main Calibration Inputs Modes: Here we consider the following segment calibration input modes: a) Smoothing Best fit Splines, and b) Segment Best Fit Response Inputs with Constraints.



Smoothing Best Fit Splines


1.      Definition: Here the treatment is limited to within a segment. In this, the segment coefficients are calibrated to the following inputs:

2.      Nomenclature:

3.      Spline Set Setup:

4.      Best Fit Penalizer Setup:

5.      Curvature Penalizer Setup:

6.      Second Derivative:

a.       !

7.      Joint Linearized Minimizer Setup:

8.      Number of Equations/Unknowns Review:

a.       For intermediate segments, the following equations determine the unknowns:

                                                   i.      Number of Continuity Constraints =>

                                                 ii.      Number of Left/Right Node Values => 2

                                                iii.      Number of Roughness Penalizer Constraint => at least 1

                                               iv.      Thus, minimal number of degrees of freedom on a per-segment basis: . This will be the number of “free” parameters we will use for to extract for each segment.

b.      For left most segments, the following equations determine the unknowns:

                                                   i.      Number of Left/Right Node Values => 2

                                                 ii.      Number of Roughness Penalizer Constraint => at least 1

                                                iii.      Thus, for the set of  parameters, the number of undetermined parameters:

c.       For the span as a whole, the number of degrees of freedom/undetermined parameters is . You may determine:

                                                   i.      The right-most second derivative, AND

                                                 ii.      Possibly, the left-most second derivative

                                                iii.      For , this will complete the set of undetermined coefficients.



Segment Best Fit Response with Constraint Matching


1.      Purpose: Here we assume that a linear transform exists the hidden state quantification metric and the measurement manifest metric.

2.      Caveat with the Segment-wise Representation: Optimizing on certain constraints (such as multi-segment constraints) now ends up producing a highly non-sparse, dense matrix. This is simply a reflection on the multi-segment spanning nature of the constraint and the eventual optimization.

3.      Constraint and Least-Square Spec:  where  and . Note that when the hidden state quantification metric is identically the same as the measurement manifest metric,  and . This corresponds to computing the least-square minimization over the observations.

4.      Constraint Formulation Development: , or , where . The parallel between this and the original least-squares formulation can now be extended in a straightforward manner.

5.      Weighted Constrained LSM:

6.      Optimization of :

Spline Jacobian





1.                  Chain Rule vs. Matrix Operations of Linear Basis Function Combination: When it comes to extracting variate Jacobian of coefficients from boundary inputs, these are absolutely equivalent – in fact, the coefficient Matrix is in reality a Jacobian itself.

·        Matrix entry as a Jacobian => Every entry of the Matrix  where  is actually a Jacobian entry, i.e., .

2.                  Self-Jacobian: Given an ordered pair  that needs to be interpolated/splined across, the self-Jacobian is defined as the vector . More generally, the self-Jacobian may be defined as  where  is an input.

o                   Self-Jacobian tells you the story of sensitivity/perturbability of the interpolant (y) on non-local points through , since the  transmission occurs through . Within a single segment, quadratic and greater splines cause fairly non-banded, dispersed Jacobians, indicating that the impact is non-local; linear splines produce simple banded/tri-diagonal Jacobians; and tension splines produce a combination of the two depending on the tension parameter (and therefore dense within the segment).

o                   Obviously Jacobian of any function  is going to be dependent on the self-Jacobian  because of the chain rule.



Optimizing Spline Basis Function Jacobian


1.      Coefficient- Value Micro-Jacobian:  where

·        A is the matrix of the basis coefficients

·        Y is the matrix (column valued) of the values (RHS). In particular, it is the boundary segment calibration nodal values in the following order:

·        F is the matrix of the coefficients of the basis function values and their derivatives. It is the following 2D Matrix:




o        where .

2.      Coefficient-Value Micro-Jacobian: Given , the coefficient-value micro-Jack is .



Spline Input Quote Sensitivity Jacobian


1.      Segment Quote Jacobian: Formulation for quote Jacobian is different than those for the coefficient edge value Jacobian, since the former automatically figures in the design matrix in the sensitivity matrix extractor pseudo-calibration stage. Thus, the quote sensitivities are effectively external sensitivity constraints transmitted via the design matrix quote sensitivities.

2.      Quote Jacobian Matrix: The Quote Sensitivity coefficients are calibrated identically to that of the base coefficient sensitivities. This simply follows from the linearity of the quote sensitivity formulation. The only area where there is non-linearity is in the product term , and that appears only at the constraint equations. Others are identically the same.

3.      Latent State Quote Sensitivity: Spline Formulation of the Latent State automatically implies that the quote sensitivity of the latent state is restricted by the above, and is therefore also a spline in itself. This further implies that the boundary formulation is subject to the similar edge conditions as before.

4.      Terminology and Nomenclature:

    1.  => Input “Calibration Quotes”
    2.  => Number of explicit Input Node Value Matches
    3.  => Number of explicit Input Constraint Value Matches
    4.  => Number of explicit Input Derivative Value Matches
    5.  => Number of explicit Input Penalty Value Matches
  1. Explicit Input-to-Response Match: This emanates from the  node match continuity criterion:  for .
  2. Explicit Input-to-Constraint Match:  for . This may be re-written as .
  3. Explicit Input-to-Derivative Match: This emanates from the  continuity criterion:  for .
  4. Explicit Input-to-Penalizer Match: This emanates from the criterion for the curvature and length penalties:  for .

Shape Preserving Spline



Shape Preserving Tension Spline


1.      Integrated vs. Partitioned Shape Controller: Integrated Shape Controllers apply shape control on a basis function-by-basis function basis (certain basis functions such as flat/linear polynomial functions need no shape controller applied on them). Partitioned shape controllers apply shape control on a segment-by-segment basis.

2.      Shape Controller Parameter Types:

3.      Shape Control as part of Basis Function formulation:

4.      Drawbacks of Shape Control as part of Basis Function formulation:

5.      Shape Control using over-determined Basis Function Set:

6.      Drawbacks of Shape Control using Over-determined Basis Function formulation:

7.      Potentially Best of Both – Partitioned Basis and Shape Control:

8.      Drawbacks of Using Partitioned Basis Functions:

9.      Partitioned vs. Integrated Tension Splines: Partitioned splines are designed such that the interpolant functional and the shape control functional are separated by formulation (e.g., rational splines). Integrated tension splines are formulated such that the shape preservation is an inherent consequence of the formulation, and there is no separation between the interpolant and the shape control functionality.

10.  Explicit Shape Preservation Control in Partitioned Splines: , where  is the interpolant, and  is the shape controller. Typically  is determined (among other things) by the continuity criterion , and  contains an explicit design parameter for shape control (for e.g.,  in the case of rational splines).

11.  Shape Control Design: Asymptotically, depending on the shape design parameter ,  should switch between linear and polynomial (i.e., typically cubic – Qu and Sarfraz (1997)). Further, design  such that , so that  and .

12.  Rational Cubic Spline Formulation:

13.  Rational Cubic Spline Coefficients:

14.  Rational Cubic Spline Derivatives:

15.  Designing  for the Segment Inflection/Extrema Control:

16.  Co-convex choice for : A similar analysis can be done to make the spline co-convex, but the corresponding formulation requires a non-linear solution for .

17.  Generalized Shape Controlling Interpolator: Given a pair of points , a  spline , and a  spline , we define a shape controlling interpolator spline  by , with the constraints .

18.  Generically Partitioned Spline Derivative:

19.  Partitioned Interpolating Spline Coefficient: Given ,

20.  Interpolating Polynomial Splines of Degree n: Given ,

                                                   i.      Mathematical simplicity

                                                 ii.      Completeness.

21.  “Derivative Completeness” Nature of Polynomial Basis Function: One big advantage for polynomial basis functions is that they are “derivative complete” in the local as well as global sense, i.e., the  basis polynomials are sufficient to uniquely determine the  continuity constraint. This is not true of non-polynomial basis functions (exponential basis functions, for e.g., need an infinite number of derivatives for completely derivative coefficient determination), therefore their shape needs global determination.

22.  Polynomial Interpolating Spline Coefficient micro-Jack:

·        ,,

·        ,,

·        ,,

23.  Curvature Design in Integrated Tension Splines: Cubic spline is interpolant on  across the nodes, and linear spline is interpolant on y. Thus,  (the tension spline interpolant) offers the tightness vs. curvature smoothness trade-off.

·        Tightness vs. Smoothness Generalization =>  is linear in x, given k is even. Of course, for  this describes a tension spline (hyperbolic or exponential). Schweikert (1966) used to improve the shape preservation characteristics.

24.  Basis Function Interpolant:

·         that is linear in x is satisfiable only by hyperbolic and exponential basis splines.

·         that is linear in x is satisfiable by hyperbolic, exponential, or sinusoidal basis splines.

·        More generally,  that is linear in x, and where  and  is satisfied only by hyperbolic and exponential splines.

·         that is linear in x, and where  and  is satisfied only by hyperbolic, exponential, or sinusoidal splines.

25.  Integrated Tension Spline Types: Sets containing both exponential and hyperbolic basis splines and a linear spline satisfy .

26.  Exponential Basis Functions:

·        Base Segment Formulation =>



·        Global <-> Local =>





·        Local <-> Global =>





·        Co-efficient Calibration =>





·        Coefficient to Input Sensitivity Grid =>





·        Local Derivatives =>





·        Global <-> Local Derivatives =>




27.  Hyperbolic Basis Functions:

·        Base Segment Formulation =>



·        Global <-> Local =>







·        Coefficient to Input Sensitivity Grid =>





·        Local Derivatives =>




o                    if r is even.

o                    if r is odd.

·        Global <-> Local Derivatives =>




28.  Alternate specifications of the segment interpolation (Trojand (2011)). Renka (1987) provides techniques for setting σ under several circumstances:

o                   Finding σ when f is bound.

§                     To get the minimum tension factor required we need to find the zeros of f’ (Renka (1987)).

o                   Finding σ when f’ is bound.

§                     To get the minimum tension factor required we need to find the zeros of f’ (Renka (1987)).

o                   Finding σ from the bound values of convexity/concavity (Renka (1987)).

29.  Problems with Hyperbolic/Tension Splines:

·        Hyperbolic and exponential functions are time consuming to compute (Preuss (1976)), Lynch (1982)).

·        They are somewhat unstable to wide parameter ranges (Spath (1969), Sapidis, Kaklis, and Loukakis (1988)).

·        In certain cases, reasonable alternatives have been provided by  splines (Nielson (1974)) and rational splines.



Shape Preserving  Splines


1.      Generic  Spline Formulation: Approach here is somewhat similar to Foley (1988), although different language/symbology.

·        p-set Basis Splines per each Segment.

·        n Data Points

·        Penalty of degree m

·         Continuity Criterion

·        Data Point Set:

·        Spline Objective Function:

2.      Number of Unknowns Analysis: In the above, , and .

·        Number of equations from the end points per segment => 2.

·        Number of equations from the coefficients determined by the  Continuity Criterion: .

·        Number of equations from the Shape Optimization Formulation: .

·        Total number of equations: .

·        Number of coefficients per segment => .

3.      Node matching constraints: Given that we are examining shape preserving splines, on applying the node match criterion  to  formulated earlier, we get  where  is the node matched Spline Objective Function.

4.      Generic Curvature Optimization Formulation: Using the above, the curvature optimization for spline basis function inside a local segment i corresponds to .

5.      Generic Curvature Optimization Minimizer: Given the basis function set , .

6.      Generic Coefficient Constrained Optimization Setup:  where .

7.      Polynomial Formulation for : For the set of polynomial basis functions, we set  on a segment-by-segment basis.

·        We also seek to optimize  on a per-segment basis by re-casting  to .

8.      : .

9.      : .

10.  Minimization of : , .

·        Here .

·        .

·        Since  and ,  corresponds to the minimum of .

·        Thus, if  and ,  becomes .

11.  Polynomial  Splines – Number of unknowns:

·        Number of coefficients (unknown) =>

·        Number of Nodal Start/End Values (known) => 2

·        Number of Calibrated coefficients from the  criterion (known):

·        Net number of unknowns: .

12.  Ordered Unknown Coefficient Set in Polynomial  Splines: Given that ,  through , as well as , are known.

·         where  are the unknown coefficients.

·        For e.g., for  cubic polynomial spline, the number of unknowns are

13.  Maximum number of equations available from Optimizing  Splines: Number of equations available from the optimization is .

·        Determinacy criterion => Thus if , or , there are no solutions!

·        Alternatively, for completeness, derive m from k as  for completeness.

·        Finally, if , optimizing run is needed.

14.  Advantage of Basis Curve Optimizing Formulation: This formulation can readily/easily incorporate linearized constraints in an automatic manner – as long as the explicit constraints are re-cast to be specified with the current segment.



Alternate Tension Spline Formulations


1.      Kaklis-Pandelis Tension Spline: As described in Kaklis and Pandelis (1990), here , where , and  is the Kaklis-Pandelis shape-controlling tension polynomial exponent.

·         corresponds to the cubic spline interpolant on .

·         corresponds to linear interpolant on .

2.      Manni’s Tension Spline: The methodology is explained in detail in Manni (1996a), Manni and Sablionniere (1997), and Manni and Sampoli (1998). Here,  on  where and are cubic polynomials. Further, is strictly increasing in , so that  is well defined (Manni (1996b)).

·        The boundary conditions are: ; further, we impose that , , , and (see Manni (2001)). The claim is that if , , thus becomes cubic. Also if ,  reduces to linear.

3.      KLK Splines: Next section is completely devoted to this.

Koch-Lyche-Kvasov Tension Splines





1.      Exponential B-Spline Specification: Expounded in detail in Cline (1974), Koch and Lyche (1989), Koch and Lyche (1993), and Kvasov (2000). First extend the knot set with 6 new points , ¸¸¸, and  such that  and , but arbitrary otherwise.

2.      Exponential Hat Functions:

3.      Properties of the  Splines:  as defined above is the basis on top of which all the higher order splines are built.  is non-zero only for , where .

4.      Layout of Base Monic Setup: With reference to figure 9, the monic basis function  may be estimated from the corresponding primitive hat functions  and  (referred to as A and B respectively in Figure 9) as:

5.      Monic B-Spline Normalizer:

6.      Monic B-Spline Cumulative Normalized Integrand:

7.      Monic B-Spline Scaled Integrand:

8.      Monic B-Spline Scaled Integrand: .

9.      Quadratic and Cubic Exponential Tension Splines: Higher order splines are recursively defined from  where:

Here . Further,  and  correspond to quadratic and cubic tension splines, respectively.

10.  Similarities between Exponential Tension B-Splines and Polynomial B-Splines: Notice the similarities, the iterative higher-order definitions, and the partition of unity as well.

11.  Cubic Exponential Tension B-Spline: This corresponds to the  case, i.e., , with validity in the interval .

12.  Piecewise Cubic Interpolant Expansion: Remember that, no matter what the basis tension functions are, for piecewise tension  continuity, they are expected to satisfy  for , i.e., this entity varies linearly across the segment.

13.  Tension Spline Curvature Penalizing Norm: The pure curvature penalizer may now be altered to become a curvature + length penalizer. Thus, . Notice that both  and  (i.e., separated squares) are included individually in the set up.

14.  Constrained Optimizer Estimate for : If the RMS best-fit error is to be limited to  where  is an extraneously specified closeness of fit metric, the constraint may be expressed as .

15.  Drawbacks of the above method: This involves yet another non-linear root extractor. Other non-linear root extraction parts in curve building are:

Thus, the stability of the precision norm technique outlined above is riddled with challenges.

16.  Parallel with Hagan-West Forward Interpolator: In Hagan-West (Hagan and West (2006)) minimalist quadratic interpolator, the segment length is incorporated in a slightly different way – as a minimizer of the  jumps at the knots (i.e., if , minimizing  at the knots).

17.  The other Tension Splines: They all have the property that the tension parameter moves smoothly from cubic (low tension) to linear (high tension), and have different forms for  and . These forms may make them computationally less expensive too.

o       Non-uniform Rational Cubic Tension Spline with linear Denominator =>



o       Non-uniform Rational Cubic Tension Spline with Quadratic Denominator =>



o       Non-uniform Exponential Rational Spline =>



18.  Tension Implied by the Basis Function Set: Given the tension  interpolant relation  inside a given segment , we can infer  from .

19.  Caveat for using KLK-type Splines for Local State Shape Proxying: Often (and this may be true for other B Splines as a whole too), the B Spline basis choice may produce segment node edge values and their corresponding derivatives of zero. In this case, you may have a singular calibration matrix that does not calibrate. In particular, this is the case for iterated B Splines that are constructed to vanish and fade rapidly at the edges. This poses for problems for segment-local splines (that may span between 0 and 1 within a given segment).

Smoothing Splines



Penalty Minimization Risk Function


1.      Penalty Minimizer Estimator Metric: Choice of the “normalized curvature area” shown in figures 5) and 6) are two possible penalty estimator choices. Obviously, closer the area is to zero, the better the penalizing match is.

2.      Dimensionless Penalizing Fit Metric: Choosing the representation in 5), and recognizing that the segment is set in the flat base , we can derive the representation in 7).

3.      Dimensionless Penalty Estimator (DPE): Using Figure 7), we now define DPE as

4.      Pros/Cons of the above Representation of DPE: If the basis functions have near-delta functional forms (Figure 8), DPE will still remain , and the metric is not very meaningful in that case. Fortunately, such delta-type basis functions are rare.

5.      Aggregate DPE Measure: Need a consolidated DPE metric that spans across all the segments in a span, i.e., the span DPE.



Smoothing Splines Setup


1.                  Process Control using Weights: Dimensionless units (such as Reynolds’ number) can effectively account for the ratio of competing natural forces. Similar use can be done for process control to be able to guide/control between 2 or more competing objectives. For example in the instance of the smoothing spline:

·        First Objective => Closeness of match using the most faithful reproducer, or curve fit.

·        Second Objective => Smoothest curve through the given points, without necessarily fitting them – of course, “smoothest” possible “curve” is a straight line.

2.                  Penalizing Smootheners: Penalizing smootheners are the consequence of Bayes estimation applied on the Quadratic Penalties with Gaussian Priors (also referred to with maxim “The Penalty is the Prior”).

3.                  Smoothing Spline Formulation: Given , and the function  that fits the points  from  (see Hastie and Tibshirani (1990) and Smoothing Spline (Wiki)). The smoothing spline estimate  is the minimizer

·         is the Spline Objective Function.

·         is needed to the left term to make it finite as , otherwise  will also have to be infinite.

·        The derivative “k” corresponds to what makes  linear. Thus, for cubic splines, k = 2.

4.                  Bias Curvature/Variance Fit Trade-off: Smaller the , the more you will fit for bias (low curvature penalty). Bigger the , more you fit for curvature/roughness penalty.

5.                  Curvature Penalty Minimizer Spline: It can be theoretically shown that the curvature penalty minimizer spline is a cubic spline. Here is how.

·        First, notice that any spline of degree >= 0 can reproduce the knot inputs.

·        By default, curvature corresponds to k = 2. Thus,  varies linearly inside a segment, thus this becomes the least possible curvature.

·        Higher order splines will have a non-linear curvature.

·        Lesser order (spline order less than 3) will violate the C2 continuity constraint.

6.                  Bias Curvature/Variance Fit Trade-off: Smaller the , the more you will fit for bias (low curvature penalty). Bigger the , more you fit for curvature/roughness penalty.

7.                  Smoothing Output Criterion:

·        Speed of Fitting

·        Speed of Optimization

·        Boundary Effects

·        Sparse, Computationally Efficient Designs

·        Semi-Parametric Models

·        Non-normal Data

·        Ease of Implementation

·        Parametrically determinable Limits

·        Specialized Limits

·        Variance Alteration/inflation

·        Adaptive Flexibility Possible

·        Adaptive Flexibility Available

·        Compactness of Results

·        Conservation of data distribution moments

·        Easy Standard Errors

8.                  Smoothing vs. Over-fitting: Since  is a control parameter, it can always be attained by a parametric specification. To estimate optimal value of  against over-fitting, use one of the following other additional criteria to penalize the extra parameters used in the fit, such as the following. Each one of them comes with its own advantages/disadvantages.

·        Cross-validation

·        Global Cross-Validation

·        Akaike Information Criterion

·        Bayesian Information Criterion

·        Deviance

·        Kullback-Leibler Divergence metric

9.                  Segment Stiffness Control:  may also be customized to behave as a segment stiffener or a penalty/stiffness controller, thus providing extra knobs for the optimization control.

10.              Relation of Lagrangian to Smoothing Spline:

·        Lagrangian objective function is used to optimize a multi-variate function  to incorporate the constraint  as . Here  is the Lagrange multiplier.

·        Optimized formulation of the smoothing spline is given by minimizing the spline objective function (a form of optimization) . Here  is the spline objective function shadow price of curvature penalty, i.e., , where c is the constraint constant defined analogous to the constraint constant in the Lagrangian objective function: .



Least Squares Best Fit + Curvature + Segment Length Penalty Formulation


1.      Nomenclature:

2.      The Formulation:  where .

3.      Segment-Level Decomposition: Segment-level decomposition ensures optimal segment coefficient formulation to within the boundaries of a segment (from the least squares fit point of view) – however, not necessarily global optimum. Further, these optimal constraints provide an extra degree of freedom at the segment level, and not necessarily at the stretch/span level.

4.      Segment-Level Decomposition Formulation:  where , , and .

5.      Least Squares Minimization Review: From earlier,

6.      Curvature Penalty Minimization: Again, from earlier,

7.      Segment Length Penalty Minimization: Similar to the curvature penalty, we get,

8.      Combining it all: .



Alternate Smootheners


1.      Compendium of Smoothing Methods:

·        Kernel Smoothing with or without binning.

·        Local Regression with or without binning.

·        Smoothing Splines with or without band solvers.

·        Regression splines with fixed/adaptive knots.

·        Penalizing B Splines.

·        Density Smoothing.

2.      Kernel Bandwidth Selector: Kernel bandwidth selection is analogous to the optimal knot point selection employed in the regression spline schemes.

·        Remember that the kernel methods essentially use the periodic functions as their basis functions.

3.      Regression Splines: Here the data is simply fit to a (hugely) reduced set of basis spline functions, typically using least squares, without any smoothness penalty.

·        Penalized regression splines are pretty much the same as regression splines in that they do use a reduced set of basis splines. However, they do impose a roughness penalty. Penalized regression splines are also referred to as smoothing splines.

·        Polynomial Regression Splines do curve fitting/regression analysis using selective insertion/removal of knots. Knots are added according to the Rao criterion, and removed according to the Wald criterion.

1.      Log Splines are a customization of the polynomial regression splines targeted for density estimation. The log of the density is modeled as a cubic spline.

·        Tension Regression Splines => In addition to the curvature penalty and the least squares fit penalty, tension regression splines also penalize the segment length.

4.      Base Density Smoothing Formulation: Log-likelihood density smoothing is analogous to maximizing the multinomial likelihood histogram , where  is the empirical observation count, and  is the probability of finding an observation in the cell .

Multi-dimensional Splines



1.      Non-symmetrical multi-dimensional Variates: Again, considering 2D as an example, it makes sense to use the basis splines separately across both , as in .

2.      Bi-polynomial 2D Spline: For the 2D Segment Range  and , working in the local variate space  and , we transform the spline basis on to the local representation basis as .

3.      Bi-linear 2D Spline: This produces a  surface. Here , therefore the first derivatives (and on) are discontinuous on the grid boundaries. From the observation set , we get from the boundary match the following values for :





4.      Bi-Cubic Interpolation: This produces a  surface. Here , therefore the first, second, and the first cross derivatives (i.e., , , , , , and ) are continuous across the grid boundaries. From the observation set , their first derivatives, and their cross derivatives, we get from the boundary match the following values for  as before. The common way is to cast these as a sequence of 2D relations by unraveling the continuity constraints, and thereby linearizing the formulation.

5.      Symmetrical Multi-dimensional variates: The trivial univariate ordering  needs revising in the context of certain multivariates, e.g., symmetrical multivariates (Smith, Price, and Lowser (1974), Graham (1983), and Lee (1989)).

·        A general “distance from focal node”  makes to more sense to set in the ascending order. Thus , where  are the multivariate nodes corresponding to the focal node.

·        Alternatively, the distance from the prior node parametrizer  may also work.

·        Use Cartesian/radial/axial basis functions to formulate the segments in terms of the surface vector coefficients in “symmetrical variate” situations.

6.      Surface Energy Minimization: Surface energy minimization using the “sigma” tension parameter – formulate equation.

·        Thin plate splines are an effective way to achieve surface energy minimization, i.e., for a 2D surface, the smoothing spline surface may be created by the minimization of the following least squares surface spline objective function . Again, apparently this is more appropriate if  are symmetrical.

7.      Elastic Maps Method for Manifold Learning: This method combines the least squared penalty for the approximation error with the bending/torsional-stretching penalty for the proxy manifold. It then uses a coarse discretization to extract the solution for the optimization problem.




·        Akima, H. (1970): A New Method of Interpolation and Smooth Curve Fitting based on Local Procedures, J. Association for Computing Machinery 17 (4): 589-602.

·        Andersen, L. (2005): Discount Curve Construction with Tension Splines, SSRN eLibrary.

·        Bartels, Beatty, and Barsky (1987): An introduction to Splines for use in Computer Graphics and Geometric Modeling.

·        Butland, J. (1980): A method of interpolating reasonable-shaped curve through any data, Proc. Computer Graphics 80: 409-422.

·        Catmull, E., and R. Rom (1974): A Class of local Interpolating Spline, in Computer Aided Geometric Design, R. E. Barnhill and Reisenfled (Eds.), Academic Press, New York.

·        Chen, W. (2009): Feedback, Nonlinear, and Distributed Circuits. CRC Press.

·        Cline, A. K. (1974): Scalar and Planar Valued Curve Fitting using Splines under Tension Communications of the ACM 17: 218-223.

·        Costantini, P., and R. Morandi (1984): Monotone and convex cubic spline interpolation, Calcolo 21: 281-294.

·        Craven, P., and G. Wahba (1979): Smoothing Noisy Data with Spline Function: Estimating the correct Degree of Smoothness by the Method of Generalized Cross Validation, Numerische Mathematik. 31: 377-403.

·        Delbourgo, R. (1989): Shape preserving interpolation to convex data to rational functions with quadratic numberator and linear denominator, IMA J. Numer. Anal. 9: 123-136.

·        Delbourgo, R., and J. Gregory (1983):  rational spline quadratic interpolation to monotonic data, IMA J. Numer. Anal. 3: 141-152.

·        Delbourgo, R., and J. Gregory (1985a): The determination of the derivative parameters for a monotonic rational quadratic interpolant, IMA J. Numer. Anal. 5: 397-406.

·        Delbourgo, R., and J. Gregory (1985b): Shape preserving piecewise rational interpolation, SIAM J. Sci. Stat. Comput. 6: 967-976.

·        De Boor, C. (2001): Practical Guide to Splines (Revised Edition), Springer.

·        De Boor, C. and B. Schwartz (1977): Piecewise Monotone Interpolation, Journal of Approximation Theory 21: 411-416.

·        Dougherty, R., A. Edelman, and J. Hyman (1989): Non-negativity, Monotonicity, and Convexity Preserving Cubic and Quintic Hermite Interpolation, Mathematics of Computation 52 (186): 471-494.

·        Duchon, J. (1976): Interpolation des Fonctions de deux Variables suivant le Principe de a Minces, RAIRO Analyse Numerique 10: 5-12.

·        Duchon, J. (1977): Splines minimizing Rotation-Invariant Semi-norms in Sobolev Spaces, Lecture Notes in Mathematics (ZAMP) 57: 85-100.

·        Eilers, P. H. C., and B. D. Marx (1996): Flexible Smoothing with B-Splines and Penalties, Statistical Science 11 (2): 89-121.

·        Epperson (1998): History of Splines, NA Digest, 98 (26).

·        Fan, K., and Q. Yao (2005): Non-linear time series: parametric and non-parametric methods. Springer.

·        Ferguson, J. (1964): Multi-variable curve interpolation, J ACM 11 (2): 221-228.

·        Foley, T. (1988): A Shape preserving Interpolant with Tension Controls, Computer Aided Geometric Design 5.

·        Fritsch, F. N., and Butland, J. (1984): A method for constructing local monotone piecewise cubic interpolants. SIAM J. Sci. Stat. Comput. 5: 300-304.

·        Fritsch, F. N., and Carlson, R. E. (1980): Monotone piecewise cubic interpolation. SIAM J. Numer. Anal. 17: 238-246.

·        Goodman, T. (2002): Shape preserving interpolation by curves.

·        Graham, N. Y. (1983): Smoothing with Periodic Cubic Splines.

·        Gregory, J. (1984): Shape preserving rational spline interpolation, in Rational Approximation and Interpolation, Graves-Morris, Saff, and Varga (Eds.), Springer-Verlag, 431-441.

·        Gregory, J. (1986): Shape preserving spline interpolation, Computer Aided Design 18: 53-58.

·        Hagan, P., and G. West (2006): Interpolation Methods for Curve Construction, Applied Mathematical Finance 13 (2): 89-129.

·        Hastie, T. J., and R. J. Tibshirani (1990): Generalized Additive Models, Chapman and Hall.

·        He, X., and P. Ng (2006): COBS – Qualitatively Constrained Smoothing via Linear Programming Computational Statistics 14 (3): 315-337.

·        Huynh, H. (1993): Accurate Monotone Cubic Interpolation, SIAM Journal on Numerical Analytisis 30 (1): 57-100.

·        Hyman, J. M. (1983): Accurate Monotonicity Preserving Cubic interpolation, SIAM Journal of Scientific and Statistical Computing 4 (4): 645-654.

·        Iwashita, Y. (2013): Piecewise Polynomial Interpolations, Open Gamma Technical Report.

·        Judd, K. (1998): Numerical Methods in Economics. MIT Press.

·        Kaklis, P. D., and D. G. Pandelis (1990): Convexity preserving polynomial splines of non-uniform degree, IMA J. Numer. Anal. 10: 223-234.

·        Katz, M. (2011): Multi-variable Analysis: A Practical Guide for Clinicians and Public Health Researchers. Cambridge University Press.

·        Koch, P. E., and T. Lyche (1989): Exponential B-Splines in Tension, in Approximation Theory VI: Proceedings of the Sixth International Symposium on Approximation Theory, vol. II, C. K. Chui et. al. (eds.), Academic Press, Boston 361-364.

·        Koch, P. E., and T. Lyche (1993): Interpolation with Exponential Splines in Tension, in Geometric Modeling, Computing/Supplementum 8, G. Farin et. al. (eds.), Springer Verlag, Wien 173-190.

·        Kruger, C. J. C. (2002): Constrained Cubic Spline Interpolation for Chemical Engineering Applications.

·        Kvasov, B. (2000): Methods of Shape-Preserving Spline Approximation World Scientific Publishing Singapore.

·        Lagrange Polynomial (Wiki): Wikipedia Entry for Lagrange Polynomial.

·        Lamberti, P., and C. Manni (2001): Shape preserving  functional interpolation via parametric cubics, IMA J. Numer. Anal. 28: 229-254.

·        Lee, E. T. Y. (1989): Choosing Nodes in Parametric Curve Interpolation.

·        Le Floc’h, F. (2013): Stable Interpolation for the Yield Curve, Calypso Technology Working Paper Series.

·        Lynch, R. (1982): A Method for Choosing a Tension Factor for a Spline under Tension Interpolation. M. Sc, University of Texas, Austin.

·        Manni, C. (1996a):  comonotone Hermite interpolation via parametric cubics, J. Comp. App. Math. 69: 143-157.

·        Manni, C. (1996b): Parametric shape preserving Hermite interpolation by piecewise quadratics, in Advanced Topics in Multi-variate Approximation, Fontanella, Jetter, and Laurent (Eds.), World-Scientific, 211-228.

·        Manni, C. (2001): On shape preserving  Hermite interpolation, BIT. 14: 127-148.

·        Manni, C., and P. Sablonniere (1997): Monotone interpolation of order 3 by  cubic splines, IMA J. Numer. Anal. 17: 305-320.

·        Manni, C., and M. L. Sampoli (1998): Comonotone parametric Hermite interpolation, in Mathematical Methods for Curves and Surfaces II, Daehlen, Lyche, and Schumaker (Eds.), Vanderbilt University Press, 343-350.

·        McAllister, D., E Passow, and J Roulier (1977): Algorithms for computing shape preserving spline interpolation to data. Math. Comp. 31: 717-725.

·        McAllister, D., and J Roulier (1981a): An Algorithm for computing shape-preserving osculating quadratic splines. ACM Trans. Math. Software. 7: 331-347.

·        McAllister, D., and J Roulier (1981b): Shape-preserving osculating Splines. ACM Trans. Math. Software. 7: 384-386.

·        Nielson, G. (1974): Some Piecewise Polynomial Alternatives to Splines under Tension, in Computer Aided Geometric Design, R. Barnhill, R. Reisenfeld (Eds.), Academic Press, 209-235.

·        Passow, E., and J Roulier (1977): Monotone and convex interpolation. Society for Industrial and Applied Mathematics, J. Numer. Anal. 14: 904-909.

·        Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery (1992, 2007), Numerical Recipes in C: the Art of Scientific Computing 2nd Edition, Cambridge University Press.

·        Pruess, S. (1976): Properties of splines in tension, J. Approx. Theory. 17: 86-96.

·        Pruess, S. (1993): Shape Preserving  Cubic Spline Interpolation, IMA Journal of Numerical Analysis 13 (4): 493-507.

·        Qu, R., and M. Sarfraz (1997): Efficient method for curve interpolation with monotonicity preservation and shape control, Neural, Parallel, and Scientific Computations 5: 275-288.

·        Raymon, L. (1981): Piecewise Monotone Interpolation in Polynomial Type, SIAM J. Math. Anal. 12: 110-114.

·        Reinsch, C. H. (1967). Smoothing by Spline Functions.

·        Renka, R. (1987): Interpolator tension splines with automatic selection of tension factors. Society for Industrial and Applied Mathematics, J. ScL. Stat. Comput. 8 (3): 393-415.

·        Rentrop, P. (1980): An Algorithm for the Computation of Exponential Splines Numerische Matematik 35 81-93.

·        Runge’s phenomenon (Wiki): Wikipedia Entry for Runge’s phenomenon.

·        Ruppert, D., M. P. Wand, and R. J. Carroll (2003): Semiparametric Regression, Cambridge University Press.

·        Sapidis, N., P Kaklis, and T Loukakis (1988): A Method for computing the Tension Parameters in Convexity preserserving Spline-in-Tension Interpolation, Numer. Math 54: 179-192.

·        Schoenberg, I. (1946): Contributions to the problem of approximation of equi-distant data by analytic functions, Quart. Appl. Math. 4: 45-99, and 112-141.

·        Schaback, R (1973): Spezielle rationale Splinefunktionen. J. Approx. Theory. 7: 281-292.

·        Schaback, R. (1988): Adaptive Rational Splines, NAM-Bericht Nr. 60, Universitat Gottingen.

·        Schweikert, D. G. (1966): An Interpolation Curve using a Spline in Tension Journal of Mathematics and Physics 45: 312-313.

·        Schumaker, L. L. (1983): On shape-preserving quadratic spline interpolation, SIAM J.. Numer. Anal. 20: 854-864.

·        Smith Jr., R. E., J. M. Price, and L. M. Howser (1974): A Smoothing Algorithm Using Cubic Spline Functions.

·        Smoothing Spline (Wiki): Wikipedia Entry for Smoothing Spline.

·        Spath, H. (1969): Exponential Spline Interpolation, Computing 4: 225-233.

·        Spath, H. (1974): Spline Algorithms for Curves and Surfaces. Utilitas Mathematica Pub. Inc. Winnipeg.

·        Spline (Wiki): Wikipedia Entry for Spline.

·        Tanggaard, C. (1997): Non-parametric Smoothing of Yield Curves Review of Quantitative Finance and Accounting 9: 251-267.

·        Trojand, D. (2011): Splines Under Tension.

·        Utreras, F. I., and V. Celis (1983): Piecewise Cubic Monotone Interpolation: A Variational Approach. Departamento de Matematicas, Universidad de Chile, Tech. Report MA-83-B-281.

·        Van Leer, B. (1974): Towards the Ultimate Conservative