48  Data-driven functions

As early as Chapter 3 of this book, we noted that a function can be described as a table where each row stands for one set of input values together with a corresponding output value. We did not, however, make much use of the table-as-function concept. Instead, we used data tables to motivate the choice of parameters, as in linear combinations of the basic modeling functions. We called this function fitting: constructing a function that stays near the data values. (We will say more about function fitting in Block 5 where we introduce new tools for treating functions as geometrical objects.)

This chapter introduces yet another important method for constructing functions that match with data. What’s different here is that each data point will be a mandate; the function is required to go through each and every data point exactly.

48.1 Generating smooth motion

As a motivating example, consider the programming of robotic arms of the sort displayed in Figure 48.1:

Figure 48.1: Robotic arms on the factory floor.

Since this isn’t a robots course, we will simplify. The arm has a resting position. When a car frame comes into place, the arm moves so that its welding electrodes are at a specific, known place in space near the car body. Then it moves in sequence to other places where a weld is required, perhaps passing through waypoints to avoid obstacles.

The problem of converting the discrete list of weld and waypoints into a continuous signal for the actuator is an instance of a mathematical process called interpolation. In real robot arms, there are multiple joints that need to be controlled simultaneously. For our illustration, we will use a simple setup (see Figure 48.2) where the robot hand rolls along a set of rails in the y-direction and another x-rail running crosswise to the y direction.

Figure 48.2: A CNC (Computer Numerical Control) carving machine. The cutting head (which plays the role of our “hand” in the narrative) is at the bottom of the rig mounted vertically. It can ove left-to-right along the bar labeled “XCARVE”. That bar, in turn, can move to-and-fro along the two parallel panels.

The task for our example robot will be to visit the points shown in Figure 48.3 in order, taking 15 seconds to traverse the whole path.

Figure 48.3: The waypoints on a path the robot hand is supposed to follow. All the action is taking place in roughly 1x1 meter area.

Figure 48.3 shows a continuous path in $(x,y)$ coordinates together with discrete labels indicating when each waypoint is to be reached. Note that the path is not a function $y(x)$. Mathematical functions are required to be single valued, meaning that for each value of the input (in the function domain) there can be only one, unique output value. The path in Figure 48.3 often involves two or more different $y$ values for a single $x$ value. There is even a small domain of $x$ near $x=900$ where the path at any given $x$ crosses six different $y$-values.

Even so, functions can be a useful way of describing the $(x,y)$ path. The key is the plural: functions. For the path in Figure 48.3 we need two quantities varying with time in a coordinated way. One approach, familiar to navigators, is to specify direction of movement and velocity at each instant of time. Perhaps not as familiar, but more fundamental mathematically, is to specify $x$ as a function of time and, separately, $y$ as a function of time. Using this formalism, the trajectory of the robot arm will consist of two functions, $x(t)$ and $y(t)$. To build those functions, we start with the waypoints stored in the data frame Robot_stations.

t	x	y
1	496	191
2	1037	138
3	1251	191
... 16 rows in total ...
15	928	432
16	737	240

The data frame `Robot_stations”

The $x(t)$ and $y(t)$ functions in this table aren’t complete enough to operate the robot. We need to provide the $x,y$-location data in the form of two continuous functions of $t$ so that the robot, at any time $t$, can look up where it is supposed to be, what its velocity should be, and how that velocity should be changing in time (acceleration).

One strategy is to construct the functions as piecewise linear functions of $t$, like this:

Figure 48.4: Two functions $x(t)$ and $y(t)$ which describe the path shown in Figure 48.3.

It can be difficult at first glance to see the relationship between the $x(t)$ and $y(t)$ functions and the path shown in Figure 48.3. As an exercise, look specifically at the segment $9\le t\le 10$. In Figure 48.4, this is the segment connecting points 9 and 10. In the path view, you can see that on this segment $x$ changes a lot while $y$ changes only a little. Correspondingly, in the function view, ${\partial }_{t}x(t)$ is large in magnitude compared to ${\partial }_{t}y(t)$.

Each functions in Figure 48.4 is an interpolating function. You’re entitled to think of the $x(t)$ function as connecting with lines the sequence of $x$ versus $t$ coordinates from the table and similarly for $y(t)$. Each of the two functions is continuous. But, based on your work in Blocks 1 through 3, you have a richer set of concepts for interpreting those two functions.

For instance, let’s look at ${\partial }_{t}y(t)$. Since $y$ is a position along the cross rail, ${\partial }_{t}y(t)$ is the velocity in that direction. Figure 48.5 shows the velocity versus time for both the $x$ and $y$ components of the movement.

Figure 48.5: Velocity versus time time along the path defined by $x(t)$ and $y(t)$ as shown in Figure Figure 48.4.

The speed of the robot arm maxes out at about 600 mm-per-second. You can get a sense for this by moving your finger two feet in 1 second: a normal human speed of movement.

Since the original $x(t)$ and $y(t)$ functions are piecewise linear, it makes sense that the derivatives with respect to time are piecewise constant. But the robot hand is a physical thing; it has to have a velocity at every instant in time. It cannot instantaneously have an undefined velocity.

Think about what it is that causes the change from one velocity step to another. There is a motor that is spinning and changing its rate of spin, perhaps using a pulley and a belt to move the robot hand to the right position at any instant of time. Changing the velocity requires a force to create an acceleration. We can differentiate the velocity to see what the acceleration must be to create the simple piecewise linear function shown in Figure 48.4.

Figure 48.6: Acceleration in the x and the y direction for the piecewise linear path. There is a huge change in velocity at the juncture between linear segments.

Mathematically, the second derivatives ${\partial }_{tt}x(t)$ and ${\partial }_{tt}y(t)$ do not exist, because ${\partial }_{t}x(t)$ and ${\partial }_{t}y(t)$ are discontinuous. There is no physical amount of force that will change the velocity in an instant.

As an accommodation to the physical existence of the robot hand, we’ve softened the transition between consecutive velocity segments to allow it to take 0.2 seconds, ramping up from zero force 0.1 second before the hand reaches the station, to maximum force at the station, then back down to zero 0.1 second after the hand reaches the station. Consequently the actual motion is smoother and the maximum acceleration is about half that of gravity. Figure 48.7 shows the resulting trajectory which can be likened to that of a baseball player rounding a base.

Figure 48.7: A smoothed x-trajectory near station 2. The position of the station is marked with a dot.

A consequence of smoothing the trajectory is that the robot hand comes near, but does not touch the station. It misses by about 2 mm. For many human tasks that might be good enough, but for precision manufacturing a miss by 2 mm is about 1000 times more than allowed.

If you like working with practical problems, you might find a simple solution to the problem. For instance, we could have aimed the robot hand 2 mm further to the right than the actual station. In falling short by 2mm, the hand would miss the new target but cross right over the originally intended station.

Solutions like this are sometimes called ad hoc, meaning that they are so specifically tailored to one situation that they do not generalize well to slightly different problems. The next section introduces a superior approach that is much more general.

48.2 Piecewise but smooth

The approach we will take to smoothly connect the points on the path is based on ideas of derivatives and on the construction of low-order polynomials. In Block 2, we emphasized low-order polynomials up to the square term, and we will pick that up again here for demonstration purposes. For this example, we will construct only the $y(t)$ function. Constructing $x(t)$ would be done using the same procedure.

Our task is to find a function $y(t)$ to interpolate discrete points such as those shown in Figure 48.8. The discrete points are called knots¹ in the language of interpolating functions. Each knot is a coordinate pair $(t_i,y(t_i))$ shown as an orange dot in Figure 48.8.

The piecewise linear interpolating function is easily constructed and is shown as a dotted curve. As we saw in the previous section, such a function has a discontinuous first derivative. We would like something smoother, with a continuous first derivative. A curve such as the one we seek is shown as the multi-colored function.

Figure 48.8: Two interpolating functions of the four discrete points (orange). One is piecewise linear (dotted curve), the other is piecewise quadratic (multi-color curve).

The framework we will adopt for the smooth interpolating function is piecewise quadratic segments between adjacent knots. There are four knots, requiring three segments. We will call the segment $p_1(y)$ connecting the first knot to the second, with $p_2(y)$ connecting the second to the third knot and $p_3(y)$ connecting the third to the fourth knot. Each of those segments will be a second-order polynomial. To keep things organized, we will use coefficient names systematically:

$$p_1(t)\equiv a_1+b_1[t-t_1]+c_1{[t-t_1]}^2$$ $$p_2(t)\equiv a_2+b_2[t-t_2]+c_2{[t-t_2]}^2$$ $$p_3(t)\equiv a_3+b_3[t-t_3]+c_3{[t-t_3]}^2$$ The four knots are $$\left[\begin{array}{c}(t_1,x_1)\\ (t_2,y_2)\\ (t_3,y_3)\\ (t_4,y_4)\end{array}\right]$$ which you can think of as two columns of a data frame, one with the $t$-coordinates of the knots and the other with the $y$-coordinates. For the knots in Figure 48.8 the data table is

t	y
1	0.0
2	2.0
3	0.5
4	1.7

Constructing the interpolating function is a matter of making good choices for $a_1,$ $a_2,$ $a_3,$ $b_1,$ $b_2,$ $b_3,$ $c_1,$ $c_2,$ and $c_3$.

We require these things of each of the interpolating polynomials:

1. It passes exactly through the two knots marking the segment’s endpoints. That is $p_1(t_1)=y_1$ and $p_1(t_2)=y_2=p_2(t_2)$ and $p_2(t_3)=y_3=p_3(t_3)$ and, finally, $p_3(t_4)=y_4$. Note that at the interior knots where two polynomials join, the left-hand polynomial and the right-hand polynomial should exactly match the function value and each other.
2. The derivative (with respect to $t$) should match where the segments join. That is, ${\partial }_t{p}_1(t_1)={\partial }_t{p}_2(t_2)$ and ${\partial }_t{p}_2(t_3)={\partial }_t{p}_3(t_3)$. Thus, the function we want to build will be $C^1$, that is, have a continuous first derivative.

How to accomplish (1) and (2)?

Notice first that because we wrote each of the polynomials in the style of Taylor polynomials, we can read the values of $a_1$, $a_2$, and $a_3$ directly from the data table:

$$p_1(t_1)=a_1=y_1$$ $$p_2(t_2)=a_2=y_2$$ $$p_3(t_3)=a_3=y_3$$

We can find other coefficients from the requirement that the right side of each segment pass through the knot on that side. This gives:

$$p_1(t_2)=y_2=a_1+b_1[t_2-t_1]+c_1{[t_2-t_1]}^2$$ $$p_2(t_c)=y_3=a_2+b_2[t_3-t_2]+c_2{[t_3-t_2]}^2$$ $$p_3(t_c)=y_4=a_3+b_3[t_4-t_3]+c_3{[t_4-t_3]}^2$$ (Notice that $t_2-t_1$ and the like are simply numbers that can be computed from the known knot points.)

Another two conditions are that the derivatives of the polynomials from either side of each interior knot point must match at the knot point. Finding the derivatives of the segments is a simple exercise:

$${\partial }_t{p}_1(t)=b_1+2c_1[t-t_1]$$ $${\partial }_t{p}_1(t)=b_2+2c_2[t-t_2]$$ $${\partial }_t{p}_1(t)=b_3+2c_3[t-t_3]$$

Matching these derivatives at the $t_2$ and $t_3$ knot points—the interior knots where two segments come together—gives two more equations: $${\partial }_t{p}_1(t_2)=b_1+2c_1[t_2-t_1]=b_2={\partial }_t{p}_2(t_2)$$ $${\partial }_t{p}_2(t_3)=b_2+2c_2[t_3-t_2]=b_3={\partial }_t{p}_3(t_3)$$ All together, we have five equations in six unknowns: $b_1,b_2,b_3$ and $c_1,c_2,c_3$.

Plugging in the specific values $t_1$ through $t_4$, and $x_1$ through $x_4$ from the data table translates the equations for the polynomial values and derivatives gives this system of equations:

$$b_1+c_1=x_2-x_1=2$$ $$b_2+c_2=x_3-x_2=-1.5$$ $$b_3+c_3=x_4-x_3=1.2$$ $$b_1+2c_1-b_2=0$$ $$b_2+2c_2-b_3=0$$

This is not the place to go into the details of solving the five equations to find the six unknowns. (Block 5 introduces the mathematics of such things, which turns out to the same math used to find model parameters to “fit” data.) But there are some simple things to say about the task.

First, you may recall being told in high-school mathematics that to find six unknowns you need six equations. We have only five equations to work with. But it is far from true that there is no solution for six unknowns with five equations. There are in fact an infinite number of solutions. (Again, Block 5 will show the mathematics behind this statement.) Essentially, all we need to do is make up a sixth equation to identify a particular one of the infinite number of solutions. It is nice if this made-up equation reflects something interpretable about the curve.

We will choose to have the sixth equation specify what the derivative of the interpolating function should be at the far right end of the graph. That right-most derivative value will be $${\partial }_t{p}_3(t_4)=b_3+2c_3[t_4-t_3].$$ We can set this value to anything we like. For instance, in Figure 48.8 the right-most derivative is set to zero; you can see this from the curve being flat at the right-most knot point.

## Warning: All aesthetics have length 1, but the data has 101 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
##   a single row.

## Warning: All aesthetics have length 1, but the data has 101 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
##   a single row.

## Warning: All aesthetics have length 1, but the data has 101 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
##   a single row.

## Warning: All aesthetics have length 1, but the data has 101 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
##   a single row.

Figure 48.9: Four different $C^1$ piecewise quadratic functions that interpolate the knot points. The functions have different values of the derivative at the right end of the domain.

Keeping in mind the piecewise nature of the interpolating polynomial, it may seem surprising that changing the slope at $t_4$ leads to a change in value of the function almost everywhere. Yet the stiffness of the parabolic segments means that conditions in one segment have an impact on adjacent segments. In turn, the segments adjacent to these also change, a change that percolates down to every segment in turn.

Try it! 48.1

Try it! 48.1 Splining through the robot stations

Quadratic spline functions can be created with the R/mosaic qspliner() function. The second argument is a data frame giving the knot locations. The first argument is a tilde expression specifying the variables to use from the data frame.

1

2

xfun <- qspliner(x ~ t, data = Robot_stations)

yfun <- qspliner(y ~ t, data = Robot_stations)

Having constructed the functions, we can plot them in the ordinary way, here adding orange dots at the $x$-values of the knots.

Active R chunk 48.1

1

2

3

4

5

## Plot the x time series and the knots

slice_plot(xfun(t) ~ t, domain(t = 1:16), npts=500) |>

  gf_point(x ~ t, data = Robot_stations, 

           color="orange", alpha=0.5, size=4) |>

  gf_labs(y="x(t)")

The graph produced by Active R chunk 48.1 quadratic spline interpolation of the $x$-coordinates of the robot-path knots. The data hardly speak for themselves, since the interpolationg function tends to alternate between concave up and concave down in adjacent segments between the knots.

Putting together the $x(t)$ and $y(t)$ interpolating functions, each of which has that extremum between knot points, leads to an absurdly complicated path, as seen in Figure 48.10.

Figure 48.10: Connecting the robot-path knots with a piecewise quadratic polynomial, constructed to be $C^1$. The path is pretty perhaps, but absurd.

Quadratic splines are rarely used in practice. (A cubic spline provides helpful flexibility. See Section 48.3.) In Figure 48.9 you can see one of the reasons: the quadratic form is so stiff that the interpolating function tends to shift from concave up to concave down (or vice versa) at each knot point. This results in the interpolating function tending to have a local minimum or maximum between adjacent knots, even if the data themselves to not indicate such a structure.

48.3 C² smooth functions

In the previous section, we arranged the functions $x(t)$ and $y(t)$ composed from the piecewise quadratic segments to be $C^1$ smooth. (Recall that $C^1$ smooth means that the derivatives ${\partial }_{t}x(t)$ and ${\partial }_{t}y(t)$ are continuous.) We established this continuity by make sure that each segment has a value of the derivative at its end-point know that matches the derivative of the adjacent segment.

To arrange $C^2$ continuity requires that the segment include a new parameter. Most commonly, this is done by moving from quadratic segments to cubic segments. This can be done by an approach similar to that of the previous section but somewhat more elaborate. Such a $C^2$ interpolating function is called a cubic spline. Cubic splines are very commonly encountered in applications requiring interpolation.

With the ability to match piecewise cubic polynomials to a set of knots, we can easily construct the smooth path to connect the knots in Figure 48.3. Figure 48.11 shows a $C^2$ path connecting the knots. The path is constructed by plotting simultaneously the output of two functions, $x(t)$ and $y(t)$, with the input $t$ on the domain $1\le t\le 16$.

Cubic splines have smoother joints

Cubic spline functions can be created with the R/mosaic spliner() function. The second argument is a data frame giving the knot locations. The first argument is a tilde expression specifying the variables to use from the data frame.

xfun <- spliner(x ~ t, data = Robot_stations)
yfun <- spliner(y ~ t, data = Robot_stations)

Figure 48.11: Connecting the robot-path knots with a piecewise cubic polynomial, constructed to be $C^2$. This is a much smoother path than produced by interpolation with quadratic polynomials.

Figure 48.12: A cubic spline interpolation of the $x$-coordinates of the robot-path knots. The cubic spline respects the monotonicity of consecutive knot points.

Application area 48.1 —Continuous second derivatives keep motion smooth.

Application area 48.1 Avoiding infinite force

We saw that using a line-segment interpolation produces discontinuity in the derivative of the function. Mathematically, discontinuity in the velocity can be thought of as an infinite acceleration, requiring an infinite force. In the physical world, accelerations must be finite. Even if a force is large, there is often slack in connections between components and the components are not perfectly rigid.

Figure 48.13 shows motion of a robotic dog. At the start, the motors in the robot are being asked to make a straight-line transition between waypoints. The result is vibration and a tremor-like movement.

Then the motors are given a smoothly interpolated signal. (Much of the video shows the programming involved. You can skip directly to time 13:02 to see the motion with smoothly interpolated waypoints.) The smooth interpolation produces a gentle and vibration-free movement.

Figure 48.13: Link to entire video by James Bruton.

48.4 Bézier splines

The sort of interpolating functions described in the previous two sections were designed to be smooth at the $C^1$ level (quadratic spline) or the $C^2$ level (cubic spline). Such smoothness makes sense for, say, robotic motion where we want at all times to keep the force on each robot joint small.

Not all path-related design problems require such smoothness. Indeed, in some settings, non-smoothness is called for. For instance, Figure 48.14 shows the outline of a familiar shape

Figure 48.14: The outline of a letter in a computer font is often specified by a series of knot points (red dots). The path passes smoothly through some of the knot points, but has a discontinuous derivative at others.

For both quadratic and the more commonly used cubic splines, matching the derivatives on either side of a knot is essential to constructing the function. For Bézier splines, each segment is mathematically independent from every other segment. It is up to the human designer of the curve to determine whether the curves derivative should be continuous or discontinuous at the knot point between two segments.

The shape of a Bézier spline segment is established by four independent points. The first and last points determine the endpoints of the segment. Each endpoint is associated with a control point that sets the angle and “speed” with which the path leaves or enters the endpoint. You can interact with the graph in Figure 48.15 to develop an intuition.

Figure 48.15: A single Bézier segment is defined by two endpoints and two control points. Drag the control points to see how the shape of the curve is defined by them.

The curve for a given segment is gratifyingly smooth. The real power of Bézier splines stems from how segments can be connected in various ways. Figure 48.16 shows two Bézier segments that have been initialized to have a smooth junction at endpoints 4 and 5. The smoothness is set by the corresponding control points (marked 3 and 6). So long as those four points (3, 4, 5, 6) are colinear and in order, the junction will be smooth. You can alter control points 2 and 7 in any way you like; the junction will remain smooth.

Figure 48.16: Two Bézier segments can be arranged in to create a smooth or non-smooth junction between them.

Consider the path followed by a Bézier curve as it leaves one of the endpoints. Figure 48.16 has been initialized so that the tangent to the curve is horizontal at the right endpoint and almost vertical at the left endpoint. The further the control point is from the endpoint, the longer the Bézier curve will stay close to the tangent line. Another way to think of this is that the position of a control point has little impact on the shape of the curve near the opposite endpoint. You can observe this on the canvas by, say, moving control point 2 and observing the relatively little change near endpoint 4.

```

Figure 48.17: A Bézier curve leaves each endpoint in a direction that is tangent to the line drawn between the endpoint and its control point.

Algebraically, each Bézier segement is a pair of cubic functions, $x(t)$ for the x-coordinate and $y(t)$ for the y-coordinate. The input $t$ varies between 0 and 1 for each segment. The coordinate pair $(x(0),y(0))$ is one endpoint of the curve, while $(x(1),y(1))$. Each intermediate value of $t$ corresponds to a point on the interior of the curve.

The $x(t)$ and $y(t)$ functions have the same form, the difference between the functions being only the values of the end values ($x_1$ and $x_4$ for the $x(t)$ function, and similarly $y_1$ and $y_4$ for the $y(t)$ function), as well as the control-point values ($x_2$ and $x_3$ for one function, $y_2$ and $y_3$ for the other.)

$$x(t)=(1-t{)}^3{x}_1+3(1-t{)}^{2}t{x}_2+3(1-t)t^2{x}_3+t^3{x}_4$$ and $$y(t)=(1-t{)}^3{y}_1+3(1-t{)}^{2}t{y}_2+3(1-t)t^2{y}_3+t^3{y}_4$$

Calculus history—Splines in wood and metal

Before the advent of digital design and manufacturing, smooth curves were described by clay or wooden models hand-crafted by skilled workers. Material was removed to conform to the models by machine tools directed by cams running over the models, by hand sanding and polishing, as shown in this video of propeller manufacture during World War II.

Video 48.1: Machining in the pre-digital world: airplane propellers

Spline functions and digital actuators have largely replaced such analog models.

---

1. Called such possibly because the curves are tied together at each of the knots.