Multivariate Data Functions

This page describes AMA Spline Library functions which compute spline approximations and interpolants of independent variable data ${\bf X}\in{\rm R}^n$ and dependent variable data ${\bf Z}\in{\rm R}^1$.

The data is given as $({\bf X}_\ell,z_\ell,w_\ell,\epsilon_\ell)$ for $\ell=1,\ldots,N$ and the spline approximation is given as

\[ s({\bf X}) = \sum_{j_n=1}^{m_n}\cdots\sum_{j_1=1}^{m_1}\alpha_{j_1,\ldots,j_n}B_{d_1,j_1}(x_1|{\bf\Lambda}_1)\cdots B_{d_n,j_n}(x_n|{\bf\Lambda_n}) \]

where $n$ is the number of independent variables for the spline and ${\bf X}\in{\rm R}^n$ is a vector of independent variables ${\bf X} = (x_1,x_2,\ldots,x_n)$, where the $x_k$'s, for $1\le k\le n$; are the spline's independent variables. The $\alpha_{j_1,\ldots,j_n}$'s are the coefficients of the tensor product B-splines

\[ \Phi_{j_1,\ldots,j_n}({\bf X}) = B_{d_1,j_1}(x_1|{\bf\Lambda}_1)\cdots B_{d_n,j_n}(x_n|{\bf\Lambda}_n) \]

where the $B_{d_k,j_k}(x_k|{\bf\Lambda}_k)$ 's are the $m_k$ univariate B-splines of degree $d_k$ defined by the knot vector

\[ {\bf\Lambda}_k = ( \lambda_{k,1}, \cdots, \lambda_{k,d_k+1}, \lambda_{k,d_k+2}, \cdots,\lambda_{k,m_k}, \lambda_{k,m_k+1}, \cdots, \lambda_{k,m_k+d_k+1})^T. \]

The number of spline coefficients is $M=\prod_{k=1}^nm_k$ and the number of univariate B-spline's $m_k$ and the knots $\lambda_{k,i}$, for $i=1,\ldots,m_k+d_k+1$ and $1\le k\le n$, depend on the method employed to approximate the data. The five methods available in the AMA Spline Library are Least Squares Approximation, Approximation, Interpolation, Monotonic Approximation and Monotonic Interpolation. If Interpolation or Monotonic Interpolation is being performed, then the independent variable data ${\bf X}_\ell\in{\rm R}^n$ and the dependent variable data $z_\ell$, for $\ell=1,\ldots,N$, must be provided. If a Least Squares Approximation is being performed, then the weights $w_\ell$, for $\ell=1,\ldots,N$, and the knot vectors ${\bf\Lambda}_k$, for $1\le k\le n$, must be provided along with the independent and dependent variable data. And, if Approximation or Monotonic Approximation is being performed, then the approximation tolerances $\epsilon_\ell$, for $\ell=1,\ldots,N$, must be provided in addition to the independent and dependent variable data.

The AMA Spline Library provides functions for computing approximations and interpolants of multivariate data where the independent variable data distribution is either random; lies upon an underlying rectilinear grid with some data missing, the so called gridded data with holes distribution; or defines a rectilinear grid.

The independent and dependent variable data for randomly distributed data and gridded data with holes is given as $({\bf X}_\ell,z_\ell)$, for $\ell=1,\ldots,N$. For gridded data with holes the independent variable data is a subset of the rectilinear grid $ (x_{1,1},x_{1,2},\cdots,x_{1,N^g_1-1},x_{1,N^g_1})\times\cdots\times(x_{n,1},x_{n,2},\cdots,x_{n,N^g_n-1},x_{n,N^g_n}), $ that is, $N < \prod_{k=1}^nN^g_k$. If the data is randomly distributed, then the independent variable data distribution has no discernable pattern. In either case, a Least Squares Approximation requires the weights $w_\ell$, for $\ell=1,\ldots,N$; while Approximation and Monotonic Approximation require the approximation tolerances $\epsilon_\ell$, for $\ell=1,\ldots,N$.

For gridded data the independent variable data is given as $x_{k,\ell_k}$, for $\ell_k=1,\ldots,N^g_k$ and $1\le k\le n$, and it defines the rectilinear grid $ (x_{1,1},x_{1,2},\cdots,x_{1,N^g_1-1},x_{1,N^g_1})\times\cdots\times(x_{n,1},x_{n,2},\cdots,x_{n,N^g_n-1},x_{n,N^g_n}). $ The dependent variable data lies on the grid and is given as $z_{\ell_1,\cdots,\ell_n}$ for $\ell_1=1,\ldots,N^g_1,\cdots,\ell_n=1,\ldots,N^g_n$; that is, there are $N=\prod_{k=1}^nN^g_k$ dependent variable values. Also, a Least Squares Approximation requires the weights $w_{\ell_1,\cdots,\ell_n}$ while Approximation and Monotonic Approximation require the approximation tolerances $\epsilon_{\ell_1,\cdots,\ell_n}$, for $\ell_1=1,\ldots,N^g_1,\cdots,\ell_n=1,\ldots,N^g_n$.

Examples of these three data distributions which are based on the river current function presented in Zermelo's Navigation are shown in Figure Data Distributions. The figure consists of three panels in which the data is represented as red circles and the yellow curves represent the river's banks. An example of randomly distributed data is illustrated in the left hand panel, the center panel illustrates gridded data with holes and the right hand panel illustrates gridded data. For the gridded data and gridded data with holes the underlying independent variable grid consists of 31 $\times$ 31 uniformly distributed points defined over the region $[0.0,2\pi]\times[-2.0,2.0]$. For the gridded data all 961 data points are included while for the gridded data with holes the onshore data is omitted, leaving 465 data points within the river's banks. The randomly distributed data consists of 961 points within the region $[0.0,2\pi]\times[-2.0,2.0]$ where the onshore points are omitted, leaving a total of 481 data points within the river's banks.

Data Distributions
Random Data Gridded Data with Holes Gridded Data

The AMA Spline Library functions for computing spline approximations and interpolants of multivariate data are separated into two categories. Those which compute approximations and interpolants of randomly distributed data or gridded data with holes and those which compute approximations and interpolants of gridded data. These functions are presented in the following two sections.

Multivariate Random Data Functions

This section presents functions which compute spline approximations and interpolants of independent variable data ${\bf X}\in{\rm R}^n$ and dependent variable data ${\bf Z}\in{\rm R}^1$ where the independent variable data is randomly distributed or has a gridded data with holes distribution. The five methods available in the AMA Spline Library for approximating multivariate data which is either randomly distributed or has a gridded data with holes distribution are Least Squares Approximation, Approximation, Interpolation, Monotonic Approximation and Monotonic Interpolation. These five approximation and interpolation methods are provided by the following AMA Spline Library functions:

Because the knot vectors are specified for the least squares function AMA_MltvLstsqr() it can be used on large and extremely randomly distributed data sets. However, the approximation and interpolation functions define knot vectors based on the underlying grid upon which the data resides. An example of this situation is illustrated in Figure AMA_MltvApprox() and AMA_MltvInterp() Grid. Hence, for large randomly distributed data sets, the approximation and interpolation functions define an extremely large number of tensor product spline coefficients and they require large amounts of memory and unacceptably long computation times to compute their approximations and interpolants. Therefore, the approximation and interpolation functions are better suited for gridded data with holes distributions and small random data sets.

Although these functions are provided to compute spline approximations and interpolants of randomly distributed data or gridded data with holes distributions, they can also be used on gridded data sets. However, they do not utilize the efficient computational algorithm employed by the Multivariate Gridded Data Functions; but, they do include second order cross partials in the penalty term. Examples of the approximations produced by the Multivariate Random Data Functions are given in the following sections:

The examples are based on a specific instance of Zermelo's Problem. Ernst Zermelo was a German mathematician who first presented the problem that now bears his name. Zermelo navigation has been used to describe the motion of many things including aircraft, ships, birds and robots. Byson and Ho describe the problem as follows.

A ship must travel through a region of strong currents. The magnitude and direction of the currents are known functions of position $u(x_1,x_2)$ and $v(x_1,x_2)$ where $(x_1,x_2)$ are rectangular coordinates and $(u(x_1,x_2),v(x_1,x_2))$ are the velocity components of the current in the $x_1$ and $x_2$ directions, respectively. The magnitude of the ship's velocity relative to the water is a constant $V$. The problem is to steer the ship in such a way as to minimize the time necessary to go from a point A to a point B. The equations of motion are

\begin{eqnarray*} {dx_1\over dt} &= V\cos\theta(t) + u(x_1,x_2) \cr & \cr {dx_2\over dt} &= V\sin\theta(t) + v(x_1,x_2) \cr \end{eqnarray*}

where $\theta(t)$ is the heading angle of the ship's axis relative to the (fixed) coordinate axes and $(x_1,x_2)$ represents the position of the ship. The Zermelo's Problem considered herein involves computing the minimum time path for a ship to cross a river. The river is sinusoidal and its center line is given by

\[ x_2(x_1) = \sin(x_1) \]

and its current flows tangent to the curve

\[ x_2(x_1) = x^c_2(x_1) + w \]

where $w$ is the north-south distance from the center line. For a given maximum distance $\bar w$ the north and south shorelines are defined as

\begin{eqnarray*} x^n_2(x_1) = x^c_2(x_1) + \bar w, \cr x^s_2(x_1) = x^c_2(x_1) - \bar w. \cr \end{eqnarray*}

If the magnitude of the current velocity is given as $R(w) = \bar Rd(w)$ where $\bar R$ is constant and

\[ d(w) = {\rm exp}{\left[-\left(w\over \bar w\right)^2\right]} \]

is the current distribution across the river, then the current velocity components are

\begin{eqnarray*} u(x_1,x_2) = {\bar R \over \sqrt{1+\cos^2(x_1)}} {\rm exp}{\left[-\left(x_2 - \sin(x_1)\over \bar w\right)^2\right]} \cr v(x_1,x_2) = {\bar R \cos(x_1) \over \sqrt{1+\cos^2(x_1)}} {\rm exp}{\left[-\left(x_2 - \sin(x_1)\over \bar w\right)^2\right]}. \cr \end{eqnarray*}

In what follows the AMA Spline Library's Multivariate Random Data Functions are used to compute cubic spline approximations and interpolants of data extracted from the $u(x_1,x_2)$ component of the river current function. Approximations of both $u(x_1,x_2)$ and $u(x_1,x_2)+E(x_1,x_2)$ are considered, where $E(x_1,x_2)=U(x_1,x_2)u(x_1,x_2)$ and $U(x_1,x_2)$ is an uniformly distributed random number in the interval $[-0.01,0.01]$; that is, a maximum 1% randomly distributed absolute error is added to $u(x_1,x_2)$. The sum of squares error

\[ \sum_{\ell=1}^{N}(E(x_{1,\ell},x_{2,\ell}))^2 \]

and the maximum error

\[ \max_\ell\vert E(x_{1,\ell},x_{2,\ell})\vert \]

associated with a maximum 1% randomly distributed absolute error for the three data distributions are given in Table River Current Errors. They are provided for comparison with the maximum and sum of squares errors associated with the approximations produced by the multivariate functions discussed below.

River Current Errors
Data Distribution Random Data Gridded Data with Holes Gridded Data
Maximum Error $2.585\times 10^{-2}$ $2.716\times 10^{-2}$ $2.631\times 10^{-2}$
Sum of Squares Error $6.351\times 10^{-2}$ $5.699\times 10^{-2}$ $5.688\times 10^{-2}$

The program exampleAMA_Mltv.c was used to compute the spline approximations and interpolants illustrated in this section and to generate the VisIt files used to produce the figures.

AMA_MltvLstsqr() - Least Squares Approximation of Multivariate Data

The least squares approximation function AMA_MltvLstsqr() requires the independent variable data ${\bf X}_\ell=(x_{1,\ell},\cdots,x_{n,\ell})$, the dependent variable data $z_\ell$ and the weights $w_\ell$, $\ell = 1,\ldots,N$, along with either the knot vector ${\bf\Lambda}_k$ or ${\bf\Lambda}^o_k$, for $1\le k\le n$; see AMA_MltvLstsqr(). For this example $n=2$ and the independent variable data ${\bf X}_\ell=(x_{1,\ell},x_{2,\ell})$ is the randomly distributed data illustrated in the left hand panel of Figure Data Distributions, the dependent variable data is $z_{\ell}=u(x_{1,\ell},x_{2,\ell})$, for $\ell=1,\ldots,481$, and the weights are $w_\ell=1.0$, for $\ell=1,\ldots, 481$. The knot vectors ${\bf\Lambda}^o_1$ and ${\bf\Lambda}^o_2$ consist of 31 knots uniformly distributed over the intervals $[0.0,2\pi]$ and $[-2.0,2.0]$, respectively. These knots are chosen because they represent the knots defined by AMA_MltvApprox() and AMA_MltvInterp() to approximate the river current data which lies on the gridded data with holes distribution illustrated in the center panel of Figure Data Distributions. Hence, they help to facilitate comparing the approximations produced by the various Multivariate Random Data Functions.

In Figure Tensor Product Grid and Data Distributions the tensor product spline grid defined by the 31 uniform knots is illustrated in conjunction with the random and gridded data with holes distributions. In both distributions the regions below and above the river's banks contain no data and the tensor product spline coefficients associated with tensor product B-splines whose domains lie entirely within those regions are not unique. However, the algorithm employed by AMA_MltvLstsqr() defines those coefficients by keeping them fixed and equal to an initial value of zero. A nonzero initial value for the spline coefficients can be specified with AMA_OptionsSetCoefficients(). Within the river's banks both distributions contain areas where the tensor product spline coefficients are not uniquely defined by the data. For the gridded data with holes distribution areas along the river's banks and at the left and right hand boundaries of the data region may have non-unique coefficients. For the randomly distributed data areas of non-unique coefficients are scattered throughout the river. In particular, consider the highlighted elliptical region along the left hand boundary of the randomly distributed data shown in the left hand panel of Figure Tensor Product Grid and Data Distributions. At the bottom of the ellipse and near the lower river boundary there lies two data points which nearly align with two horizontal grid lines. And, at the top of the ellipse there are another two data points which nearly align with two horizontal grid lies. However, in the middle of the ellipse there are three horizontal grid lines which do not have any data nearby. Therefore, the tensor product spline coefficients associated with the tensor product B-splines which lie within the elliptical region are not uniquely defined by the data.

The function AMA_MltvLstsqr() has two options which are designed to alleviate the problems associated with non-unique tensor product spline coefficients. The first option is the minimum norm optimization which is performed subsequent to the least squares approximation and does not degrade the sum of squares error. The second option is the penalty term which balances between approximating the data and providing a smooth approximation. The next two examples explore these two options. The first example illustrates the effect of the minimum norm optimization on an approximation of the randomly distributed river current data $u(x_1,x_2)$. The second example demonstrates the influence of the penalty term on an approximation of the randomly distributed river current data $u(x_1,x_2)+E(x_1,x_2)$.

Tensor Product Grid and Data Distributions
Random Data Gridded Data with Holes

Least squares approximations of the randomly distributed river current data $u(x_1,x_2)$ are illustrated in Figure Least Squares Approximation of Random Data. The figure consists of three panels in which the data is represented as red spheres and the approximation is represented as a gray surface. The left hand panel illustrates a least squares approximation with the minimum norm optimization disabled. The approximation exhibits two prominent characteristics caused by non-unique spline coefficients. The first characteristic is the flat region towards the front and left of the figure which corresponds to the data less region in the lower left hand corner of Figure Tensor Product Grid and Data Distributions. Since there is no data in this region the least squares approximation remains fixed and equal to its default initial value of zero. Since the region lies outside the river's banks, where there is no current, this characteristic is reasonable. The second characteristic lies within the highlighted elliptical region shown in the left hand panel of Figure Tensor Product Grid and Data Distributions which corresponds to the spike in the approximation to the front and left in the figure. This spike is a result of the two data points which lie at the bottom of the ellipse, where the approximation rises to fit the data, and the lack of data in the center of the ellipse, where the approximation remains equal to an initial value of zero. This characteristic is not acceptable since the approximation suggests the current is zero within the river banks. Hence, any solution to the river crossing problem based on this approximation is probably unrealistic.

Least Squares Approximation of Random Data
Minimum Norm Optimization Disabled Minimum Norm Optimization Enabled Bounded Minimum Norm Approximation

A least squares approximation of the river current data based on a minimum norm optimization produces a more reasonable approximation. Such an approximation is illustrated in the center panel of Figure Least Squares Approximation of Random Data. Along with the river current data and its approximation the figure contains the surface $u(x_1,x_2)=0.0$, shown as a transparent blue green plane. Examination of the approximation in the elliptical region highlighted in Figure Tensor Product Grid and Data Distributions shows it fits the data at the top and bottom of the ellipse as well as the approximation without the minimum norm optimization; and, within the center of the ellipse it assumes behavior consistent with the data. In fact, within the entire data region the approximation behaves in a manner consistent with the data. However, outside the river's banks the approximation becomes negative as the trend established by the data within the river's banks is maintained.

In addition to the minimum norm optimization option, AMA_MltvLstsqr() has the ability to impose lower and upper bounds on an approximation. Since the river current data is everywhere positive it seems appropriate to impose a lower bound of zero. An approximation subject to such a lower bound is illustrated in the right hand panel of Figure Least Squares Approximation of Random Data along with the plane $u(x_1,x_2)=0.0$. This approximation maintains the desirable behavior of the unbounded approximation within the data region and is everywhere positive. Since positivity is a fundamental characteristic of the river current data, a lower bound of zero is imposed on all of the approximations presented throughout the remainder of this section.

The sum of squares error

\[ \sum_{\ell=1}^{481}(s(x_{1,\ell},x_{2,\ell}) - u(x_{1,\ell},x_{2,\ell}))^2 \]

and the maximum error

\[ \max_\ell\vert s(x_{1,\ell},x_{2,\ell}) - u(x_{1,\ell},x_{2,\ell})\vert \]

for the three approximations is given in Table Least Squares Approximation Errors given below.

Least Squares Approximation Errors
Minimum Norm Disabled Enabled Enabled and Bounded
Maximum Error $1.014\times 10^{-4}$ $9.189\times 10^{-5}$ $7.456\times 10^{-5}$
Sum of Squares Error $1.714\times 10^{-7}$ $1.233\times 10^{-7}$ $1.152\times 10^{-7}$

These errors indicate that, for all practical purposes, 31 uniformly distributed knots are sufficient to "interpolate" the randomly distributed river current data $u(x_1,x_2)$ and they show the minimum norm optimization does not degrade the sum of squares error. They also show imposing the lower bound of zero on the approximation does not degrade the errors. Because of their corresponding errors and the quality of the two minimum norm approximations, it's likely they will produce realistic solutions to the river crossing problem. It should be noted that the river current data $u(x_1,x_2)$ lends itself nicely to the minimum norm optimization criteria and this example is not necessarily indicative of minimum norm optimization results in general.

Besides performing the minimum norm optimization AMA_MltvLstsqr() provides an option to include a penalty term on a least squares approximation. While the minimum norm optimization is performed subsequent to the least squares approximation and does not degrade the sum of squares error, the penalty term balances between approximating the data and producing a smoother surface. Hence, the penalty term provides a means to mitigate the effect of error within the data on the approximation. Three examples of approximating the river current data $u(x_1,x_2)+E(x_1,x_2)$ utilizing the penalty term option are presented in Figure Penalty Term Approximation of Random Data with Error. All the approximations are subject to a lower bound of zero. Since the penalty term and minimum norm optimization are mutually exclusive, the minimum norm optimization is disabled. The left hand panel illustrates an approximation with the penalty term weight $\theta = 0.0$, the approximation in the center panel corresponds to $\theta = 0.01$ and the approximation in the right hand panel is based on $\theta = 0.1$.

Penalty Term Approximation of Random Data with Error
Penalty Term Weight $\theta = 0.0 $ Penalty Term Weight $\theta = 0.01$ Penalty Term Weight $\theta = 0.1 $

The approximation for $\theta = 0.0$ exhibits the same two prominent characteristics as the approximation shown in the left hand panel of Figure Least Squares Approximation of Random Data. As before, they are caused by nonunique spline coefficients associated with the distribution of the data with respect to the knots and by initializing the spline coefficients to zero. However, the approximation of $u(x_1,x_2) + E(x_1,x_2)$ shown in the left hand panel of Figure Penalty Term Approximation of Random Data with Error also has multiple spikes along the data boundary which are not evident in the approximation of $u(x_1,x_2)$. That behavior is a result of the nonunique spline coefficients accentuating the error within the data. Although the maximum and sum of squares errors given in Table Penalty Term Approximation Errors indicate the approximation is acceptable(the maximum error is less then the maximum error of $2.492\times 10^{-2}$ within the data), the magnitude of the spikes along the data boundary and within the river banks deem it unusable for the purpose of solving the river crossing problem.

The center and right hand panels of Figure Penalty Term Approximation of Random Data with Error illustrate that a nonzero penalty term weight mitigates the influence of the error within the data on the approximation; and, they show the approximation becomes smoother with increasing penalty term weight. Also, from Table Penalty Term Approximation Errors, the maximum and sum of squares errors for $\theta=0.01$ are smaller then the corresponding errors for gridded data with holes shown in Table River Current Errors; while, the maximum error for the $\theta=0.1$ approximation is larger. Continuing to increase the penalty term weight produces smoother approximations with a corresponding degradation of the maximum and sum of squares errors. In fact, a value of $\theta = 1.0$ produces a planar approximation of the data with maximum and sum of squares errors of $1.520$ and $2.206\times 10^2$, respectively.

Penalty Term Approximation Errors
Theta $0.00$ $0.01$ $0.10$
Maximum Error $2.129\times 10^{-2}$ $2.486\times 10^{-2}$ $2.836\times 10^{-2}$
Sum of Squares Error $3.934\times 10^{-3}$ $1.347\times 10^{-2}$ $2.475\times 10^{-2}$

The errors shown in Table Least Squares Approximation Errors indicate that the 31 uniform knots provided to AMA_MltvLstsqr() are sufficient to "interpolate" the river current data $u(x_1,x_2)$. But, this is a consequence of the nature of the data and in order to insure satisfaction of the approximation and interpolation constraints imposed by AMA_MltvApprox() and AMA_MltvInterp() those functions must, in general, define knot vectors based on the underlying grid upon which the data lies. With respect to its underlying grid, the randomly distributed river current data can be considered a gridded data with holes distribution with multiple holes scattered throughout the grid. This situation is illustrated in Figure Underlying Grid for Random Data which shows the randomly distributed river current data in conjunction with its underlying grid. This underlying grid is also the tensor product grid defined by AMA_MltvApprox() and AMA_MltvInterp() when they are employed to compute approximations of the randomly distributed river current data.

Underlying Grid for Random Data
AMA_MltvApprox() and AMA_MltvInterp() Grid

The tensor product grid illustrated in Figure Underlying Grid for Random Data defines in access of 230,000 tensor product spline coefficients and the amount of memory required by AMA_MltvApprox() and AMA_MltvInterp() to compute their approximations exceeds the memory limits of most computing systems. Therefore, it is not possible to employ these functions on the randomly distributed river current data. The situation worsens with AMA_MltvMonoApprox() and AMA_MltvMonoInterp() as these routines include additional knots internal to the underlying data grid in order to insure satisfaction of the local monotonicity constraints. Because the approximation and interpolation functions define an excessive number of knots, it is recommended these functions are used only on gridded data with holes distributions or small randomly distributed data sets. Hence, the examples presented in the next four sections are based on the gridded data with holes distribution illustrated in the center panel of Figure Data Distributions.

AMA_MltvApprox() - Approximation of Multivariate Data

The approximation function AMA_MltvApprox() requires the independent variable data ${\bf X}_\ell=(x_{1,\ell},\cdots,x_{n,\ell_n})$, the dependent variable data $z_\ell$ and the approximation tolerances $\epsilon_\ell$, $\ell = 1,\ldots,N$. For this example $n=2$ and the independent variable data ${\bf X}_\ell=(x_{1,\ell},x_{2,\ell})$ is the gridded data with holes distribution illustrated in the center panel of Figure Data Distributions; the dependent variable data is the river current data $u(x_{1,\ell},x_{2,\ell})$, for $\ell=1,\ldots,465$; or the data $u(x_{1,\ell},x_{2,\ell})+E(x_{1,\ell},x_{2,\ell})$ where the approximation tolerances are either $\epsilon_\ell=0.0$ or $\epsilon_\ell=E(x_{1,\ell},x_{2,\ell})$, for $\ell=1,\ldots,465$. In general, a statistical analysis of the process which produced the data would provide nonzero approximation tolerances, but $E(x_1,x_2)$ serves to illustrate the influence of the approximation tolerances on the approximation produced by AMA_MltvApprox().

Based on the data AMA_MltvApprox() defines the knot vectors ${\bf\Lambda}_0$ and ${\bf\Lambda}_1$ which lie on the independent variable data's underlying rectilinear grid. They are illustrated in the right hand panel of Figure Tensor Product Grid and Data Distributions. These knot vectors define sufficient spline coefficients to insure satisfaction of the approximation constraints $z_\ell - \epsilon_\ell \le s(x_{1,\ell},x_{2,\ell}) \le z_\ell + \epsilon_\ell$, for $\ell=1,\ldots, 465$.

Three approximations of the river current data whose independent variable data defines a gridded data with holes distribution are illustrated in Figure Approximation. The left hand and center panels illustrate approximations of $u(x_1,x_2)$ and $u(x_1,x_2)+E(x_1,x_2)$, respectively; with the approximation tolerances $\epsilon_\ell=0.0$ for all $\ell=1,\ldots,465$. The right hand panel shows an approximation of $u(x_1,x_2)+E(x_1,x_2)$ with the approximation tolerances $\epsilon_\ell=E(x_{1,\ell},x_{2,\ell})$, for $\ell=1,\ldots,465$.

Approximation
Approximation of $u(x_1,x_2)$ with $\epsilon=0.0$ Approximation of $u(x_1,x_2)+E(x_1,x_2)$ with $\epsilon=0.0$ Approximation of $u(x_1,x_2)+E(x_1,x_2)$ with $\epsilon=E(x_1,x_2)$

The sum of squares and maximum errors for the three approximations are given in Table Approximation Errors. They show the spline approximations of the river current data $u(x_1,x_2)$ and the data with a random error $u(x_1,x_2)+E(x_1,x_2)$ for $\epsilon=0.0$ are, as expected, spline interpolants of the data. In the latter case the error within the data causes ripples within the surface which are not evident in the former case, where the absence of error within the data results in a smoother surface.

Approximation Errors
Data Approximated $u(x_1,x_2)$ with $\epsilon=0.0$ $u(x_1,x_2)+E(x_1,x_2)$ with $\epsilon=0.0$ $u(x_1,x_2)+E(x_1,x_2)$ with $\epsilon=E(x_1,x_2)$
Maximum Error $4.356\times 10^{-13}$ $6.741\times 10^{-13}$ $2.692\times 10^{-2}$
Sum of Squares Error $4.146\times 10^{-25}$ $1.033\times 10^{-24}$ $4.570\times 10^{-2}$

However, as the right hand panel of Figure Approximation illustrates, including the approximation tolerances $\epsilon_\ell=E(x_{1,\ell},x_{2,\ell})$ also produces a smooth approximation. And, from Table Approximation Errors, the maximum and sum of squares errors for the approximation are smaller then the corresponding errors for gridded data with holes shown in Table River Current Errors. This is a direct consequence of the approximation constraints imposed by AMA_MltvApprox() which insure the absolute errors are within the approximation tolerances.

AMA_MltvInterp() - Interpolation of Multivariate Data

The interpolation function AMA_MltvInterp() requires the independent variable data ${\bf X}_\ell=(x_{1,\ell},\cdots,x_{n,\ell_n})$ and the dependent variable data $z_\ell$, for $\ell = 1,\ldots,N$. For this example $n=2$ and the independent variable data ${\bf X}_\ell=(x_{1,\ell},x_{2,\ell})$ is the gridded data with holes distribution illustrated in the center panel of Figure Data Distributions and the dependent variable data is either the river current data $u(x_{1,\ell},x_{2,\ell})$ or the data $u(x_{1,\ell},x_{2,\ell})+E(x_{1,\ell},x_{2,\ell})$, for $\ell=1,\ldots,465$.

Based on the data AMA_MltvInterp() defines the knot vectors ${\bf\Lambda}_0$ and ${\bf\Lambda}_1$ which lie on the independent variable data's underlying rectilinear grid. They are illustrated in the right hand panel of Figure Tensor Product Grid and Data Distributions. These knot vectors define sufficient spline coefficients to satisfy the interpolation constraints $s(x_{1,\ell},x_{2,\ell}) = z_\ell$, for $\ell=1,\ldots, 465$, and the boundary conditions ${\partial s^{2}(0.0,x_2)\over \partial x^2_1} = {\partial s^{2}(2\pi,x_2)\over \partial x^2_1} = 0.0$ and ${\partial s^{2}(x_1,-2.0)\over \partial x^2_2} = {\partial s^{2}(x_1,2.0)\over \partial x^2_2} = 0.0$, see AMA_MltvInterp() for a description of the boundary conditions it imposes.

Two interpolants of the river current data whose independent variable data defines a gridded data with holes distribution are illustrated in Figure Interpolation. The interpolant for $u(x_1,x_2)$ is shown in the left hand panel and the interpolant for $u(x_1,x_2)+E(x_1,x_2)$ is shown in the right hand panel.

Interpolation
Interpolation of $u(x_1,x_2)$ Interpolation of $u(x_1,x_2)+E(x_1,x_2)$

The sum of squares and maximum errors for the two interpolants given Table Interpolation Errors verify the function AMA_MltvInterp() does produce spline interpolants of the data. As with the approximations produced by AMA_MltvApprox() when the approximation tolerances equal zero, the error within the data causes ripples within the interpolant of $u(x_1,x_2)+E(x_1,x_2)$ which are not evident in the interpolant of $u(x_1,x_2)$, where the absence of error within the data results in a smoother surface.

Interpolation Errors
Data Approximated $u(x_1,x_2)$ with $\epsilon=0.0$ $u(x_1,x_2)+E(x_1,x_2)$
Maximum Error $5.004\times 10^{-13}$ $6.543\times 10^{-13}$
Sum of Squares Error $3.043\times 10^{-25}$ $9.824\times 10^{-25}$

With the exception of the natural boundary conditions imposed by AMA_MltvInterp(), its algorithm is equivalent to the algorithm employed by AMA_MltvApprox() when the approximation tolerances are equal to zero. Therefore, AMA_MltvApprox() can be used with the approximation tolerances equal to zero when interpolation without natural boundary conditions is desired.

AMA_MltvMonoApprox() - Monotonic Approximation of Multivariate Data

The monotonic approximation function AMA_MltvMonoApprox() requires the independent variable data ${\bf X}_\ell=(x_{1,\ell},\cdots,x_{n,\ell_n})$, the dependent variable data $z_\ell$ and the approximation tolerances $\epsilon_\ell$, $\ell = 1,\ldots,N$. For this example $n=2$ and the independent variable data ${\bf X}_\ell=(x_{1,\ell},x_{2,\ell})$ is the gridded data with holes distribution illustrated in the center panel of Figure Data Distributions; the dependent variable data is the river current data $u(x_{1,\ell},x_{2,\ell})$, for $\ell=1,\ldots,465$; or the data $u(x_{1,\ell},x_{2,\ell})+E(x_{1,\ell},x_{2,\ell})$ where the approximation tolerances are either $\epsilon_\ell=0.0$ or $\epsilon_\ell=E(x_{1,\ell},x_{2,\ell})$, for $\ell=1,\ldots,465$.

Based on the data AMA_MltvMonoApprox() defines knot vectors ${\bf\Lambda}_0$ and ${\bf\Lambda}_1$ which provide sufficient spline coefficients to satisfy the approximation constraints $z_\ell - \epsilon_\ell \le s(x_{1,\ell},x_{2,\ell}) \le z_\ell + \epsilon_\ell$, for $\ell=1,\ldots, 465$, and the local monotonicity constraints, see AMA_MltvMonoApprox() for a definition of the local monotonicity constraints it imposes.

Three monotonic approximations of the river current data whose independent variable data defines a gridded data with holes distribution are illustrated in Figure Monotonic Approximation. The left hand and center panels illustrate monotonic approximations of $u(x_1,x_2)$ and $u(x_1,x_2)+E(x_1,x_2)$, respectively; with the approximation tolerances $\epsilon_\ell=0.0$ for all $\ell=1,\ldots,465$. The right hand panel shows a monotonic approximation of $u(x_1,x_2)+E(x_1,x_2)$ with the approximation tolerances $\epsilon_\ell=E(x_{1,\ell},x_{2,\ell})$, for $\ell=1,\ldots,465$.

Monotonic Approximation
Monotonic Approximation of $u(x_1,x_2)$ with $\epsilon=0.0$ Monotonic Approximation of $u(x_1,x_2)+E(x_1,x_2)$ with $\epsilon=0.0$ Monotonic Approximation of $u(x_1,x_2)+E(x_1,x_2)$ with $\epsilon=E(x_1,x_2)$

The sum of squares and maximum errors for the three approximations are given in Table Monotonic Approximation Errors. As was the case with AMA_MltvApprox() the spline approximation of the river current data $u(x_1,x_2)$ and the approximation of the data with a random error $u(x_1,x_2)+E(x_1,x_2)$, for $\epsilon=0.0$, are actually spline interpolants of the data. The effect of the monotonicity constraints is most noticeable at the local maximum highlighted with an ellipse in the left hand panel. The ellipse contains two data points which, for the sake of discussion, are denoted $u(x_{1,\ell},x_{2,\ell})$ and $u(x_{1,\ell+1},x_{2,\ell+1})$. Note that $x_{1,\ell} < x_{1,\ell+1}$ and $x_{2,\ell} = x_{2,\ell+1}$. The river current data at these two points is such that $u(x_{1,\ell},x_{2,\ell}) = u(x_{1,\ell+1},x_{2,\ell+1})$. In this case a zero partial derivative constraint ${\partial s(x_1,x_2)\over\partial x_1}=0.0$ is imposed at $x_{2,\ell}$ over the interval $x_1\in[x_{1,\ell},x_{1,\ell+1}]$. Hence, instead of obtaining a local maximum which exceeds the data, the zero partial derivative constraint forces the approximation to maintain a local maximum consistent with the data - causing the approximation to flatten at the local maximum. Additionally, the knot vectors defined by AMA_MltvMonoApprox(), which insure sufficient degrees of freedom to interpolate the data and to satisfy the monotonicity constraints, cause the flat region to extend beyond the interval.

Monotonic Approximation Errors
Data Approximated $u(x_1,x_2)$ with $\epsilon=0.0$ $u(x_1,x_2)+E(x_1,x_2)$ with $\epsilon=0.0$ $u(x_1,x_2)+E(x_1,x_2)$ with $\epsilon=E(x_1,x_2)$
Maximum Error $8.881\times 10^{-16}$ $8.881\times 10^{-16}$ $2.716\times 10^{-2}$
Sum of Squares Error $3.993\times 10^{-29}$ $4.179\times 10^{-29}$ $3.901\times 10^{-2}$

For the approximation of $u(x_1,x_2)+E(x_1,x_2)$ given in the center panel the situation differs as, due to error, the data within the ellipse is $u(x_{1,\ell},x_{2,\ell})+E(x_{1,\ell},x_{2,\ell})$ and $u(x_{1,\ell+1},x_{2,\ell+1})+E(x_{1,\ell+1},x_{2,\ell+1})$, which satisfies $u(x_{1,\ell},x_{2,\ell})+E(x_{1,\ell},x_{2,\ell})< u(x_{1,\ell+1},x_{2,\ell+1})+E(x_{1,\ell+1},x_{2,\ell+1})$. Now, a positive partial derivative constraint ${\partial s(x_1,x_2)\over\partial x_1}> 0.0$ is imposed at $x_{2,\ell}$ over the interval $x_1\in[x_{1,\ell},x_{1,\ell+1}]$. As before, the positive partial derivative constraint forces the approximation to maintain a local maximum consistent with the data, causing the approximation to obtain a local maximum at $(x_{1,\ell+1},x_{2,\ell+1})$. Also note that in spite of the monotonicity constraints the random error within the data causes ripples in the approximation of $u(x_1,x_2)+E(x_1,x_2)$, which are not evident in the approximation of $u(x_1,x_2)$.

The approximation of $u(x_1,x_2)+E(x_1,x_2)$ using the approximation tolerances $\epsilon_\ell = E(x_{1,\ell},x_{2,\ell})$ is similar to the approximation of $u(x_1,x_2)$. In this case the data at $(x_{1,\ell},x_{2,\ell})$ and $(x_{1,\ell+1},x_{2,\ell+1})$ is considered to lie within the intervals $[u(x_{1,\ell},x_{2,\ell})-E(x_{1,\ell},x_{2,\ell}),u(x_{1,\ell},x_{2,\ell})+E(x_{1,\ell},x_{2,\ell})]$ and $[u(x_{1,\ell+1},x_{2,\ell+1})-E(x_{1,\ell+1},x_{2,\ell+1}),u(x_{1,\ell+1},x_{2,\ell+1})+E(x_{1,\ell+1},x_{2,\ell+1})]$, respectively. Since these intervals intersect a zero partial derivative constraint ${\partial s(x_1,x_2)\over\partial x_1}=0.0$ is imposed at $x_{2,\ell}$ over the interval $x_1\in[x_{1,\ell},x_{1,\ell+1}]$. This constraint causes the approximation to maintain a local maximum which lies within the aforementioned data intervals and, because of the knots defined by AMA_MltvMonoApprox(), the local maximum extends beyond the interval $x_1\in[x_{1,\ell},x_{1,\ell+1}]$.

AMA_MltvMonoInterp() - Monotonic Interpolation of Multivariate Data

The monotonic interpolation function AMA_MltvMonoInterp() requires the independent variable data ${\bf X}_\ell=(x_{1,\ell},\cdots,x_{n,\ell_n})$ and the dependent variable data $z_\ell$, for $\ell = 1,\ldots,N$. For this example $n=2$ and the independent variable data ${\bf X}_\ell=(x_{1,\ell},x_{2,\ell})$ is the gridded data with holes distribution illustrated in the center panel of Figure Data Distributions and the dependent variable data is either the river current data $u(x_{1,\ell},x_{2,\ell})$ or the data $u(x_{1,\ell},x_{2,\ell})+E(x_{1,\ell},x_{2,\ell})$, for $\ell=1,\ldots,465$.

Based on the data AMA_MltvMonoInterp() defines knot vectors ${\bf\Lambda}_0$ and ${\bf\Lambda}_1$ which provide sufficient spline coefficients to satisfy the interpolation constraints $s(x_{1,\ell},x_{2,\ell}) = z_\ell$, for $\ell=1,\ldots, 465$, and the local monotonicity constraints, see AMA_MltvMonoInterp() for a definition of the local monotonicity constraints it imposes.

Two monotonic interpolants of the river current data whose independent variable data defines a gridded data with holes distribution are illustrated in Figure Monotonic Interpolation. The monotonic interpolant for $u(x_1,x_2)$ is shown in the left hand panel and the monotonic interpolant for $u(x_1,x_2)+E(x_1,x_2)$ is shown in the right hand panel.

Monotonic Interpolation
Monotonic Interpolation of $u(x_1,x_2)$ Monotonic Interpolation of $u(x_1,x_2)+E(x_1,x_2)$

Since AMA_MltvMonoInterp() does not impose boundary conditions on its approximation, the approximation it produces is identical to the approximation produced by AMA_MltvMonoApprox() when its approximation tolerances are set equal to zero. Hence, AMA_MltvMonoInterp() is redundant but, nonetheless, is provided for convenience.

Multivariate Gridded Data Functions

This section presents functions which compute spline approximations and interpolants of independent variable data ${\bf X}\in{\rm R}^n$ and dependent variable data ${\bf Z}\in{\rm R}^1$ where the independent variable data lies on a grid.

For gridded data the independent variable data is given as $x_{k,\ell}$, for $\ell=1,\ldots,N^g_k$ and $1\le k\le n$, and it defines the rectilinear grid $ (x_{1,1},x_{1,2},\cdots,x_{1,N^g_1-1},x_{1,N^g_1})\times\cdots\times(x_{n,1},x_{n,2},\cdots,x_{n,N^g_n-1},x_{n,N^g_n}). $ The dependent variable data lies on the grid and is given as $z_{\ell_1,\cdots,\ell_n}$ for $\ell_1=1,\ldots,N^g_1,\cdots,\ell_n=1,\ldots,N^g_n$; that is, there are $N=\prod_{k=1}^nN^g_k$ dependent variable values. A Least Squares Approximation requires the weights $w_{\ell_1,\cdots,\ell_n}$ while Approximation and Monotonic Approximation require the approximation tolerances $\epsilon_{\ell_1,\cdots,\ell_n}$, for $\ell_1=1,\ldots,N^g_1,\cdots,\ell_n=1,\ldots,N^g_n$.

These functions employ an efficient algorithm for computing spline approximations and interpolants of multivariate gridded data. This algorithm consists of multiple invocations of the function's corresponding Univariate Data Functions to compute a multivariate spline which satisfies the approximation and interpolation constraints. Because these functions employ the Univariate Data Functions, they do not minimize the cross partial terms of the penalty term utilized by the Multivariate Random Data Functions. Therefore, with the exception of least squares approximation under certain circumstances, the Multivariate Random Data Functions produce different spline approximations of gridded data sets than their gridded data counterparts.

The five methods available in the AMA Spline Library for approximating multivariate gridded data are Approximation, Interpolation, Least Squares Approximation, Monotonic Approximation and Monotonic Interpolation. These methods are provided in the following AMA Spline Library functions:

Examples of the approximations produced by these functions are given in the following sections:

The examples are based on a specific instance of Zermelo's Problem. Ernst Zermelo was a German mathematician who first presented the problem that now bears his name. "Zermelo navigation" has been used to describe the motion of many things including aircraft, ships, birds and robots. See Zermelo's Navigation for a description of the problem provided by Byson and Ho.

In what follows the AMA Spline Library's Multivariate Gridded Data Functions are used to compute cubic spline approximations and interpolants of data extracted from the $u(x_1,x_2)$ component of the river current function described in Zermelo's Navigation. The data is defined over the rectilinear grid illustrated in the right hand panel of Figure Data Distributions. As in Multivariate Random Data Functions, approximations of both $u(x_1,x_2)$ and $u(x_1,x_2)+E(x_1,x_2)$ are considered, where $E(x_1,x_2)=U(x_1,x_2)u(x_1,x_2)$ and $U(x_1,x_2)$ is an uniformly distributed random number in the interval $[-0.01,0.01]$; that is, a maximum 1% randomly distributed absolute error is added to $u(x_1,x_2)$.

The program exampleAMA_MltvGrd.c was used to compute the spline approximations and interpolants illustrated in this section and to generate the VisIt files used to produce the figures.

AMA_MltvGrdLstsqr() - Least Squares Approximation of Multivariate Gridded Data

The least squares approximation function AMA_MltvGrdLstsqr() requires the independent variable data $x_{k,\ell_k}$, for $\ell_k=1,\ldots,N^g_k$ and $1\le k\le n$; the dependent variable data $z_{\ell_1,\cdots,\ell_n}$ and the weights $w_{\ell_1,\cdots,\ell_n}$, for $\ell_1=1,\ldots,N^g_1,\cdots,\ell_n=1,\ldots,N^g_n$; along with either the knot vector ${\bf\Lambda}_k$ or ${\bf\Lambda}^o_k$, for $1\le k\le n$; see AMA_MltvGrdLstsqr(). For this example $n=2$ and the independent variable data is the gridded data distribution illustrated in the right hand panel of Figure Data Distributions, the dependent variable data is $z_{\ell_1,\ell_2} = u(x_{1,\ell_1},x_{2,\ell_2})$ and the weights are $w_{\ell_1,\ell_2} = 1$, for $\ell_1=1,\ldots,31$ and $\ell_2=1,\ldots,31$.

When approximating gridded data with tensor product splines it is always possible to define knot vectors which provide sufficient spline coefficients to interpolate the data and to define an unique approximation. In the case of cubic spline approximation the knot vectors ${\bf\Lambda}^o_1 = (x_{1,1},x_{1,3},\cdots,x_{1,N^g_1-2},x_{1,N^g_1})$ and ${\bf\Lambda}^o_2 = (x_{2,1},x_{2,3},\cdots,x_{2,N^g_2-2},x_{2,N^g_2})$ are such knots. Note, the knot vectors equal the independent variable data with the second and next to last independent variable data points removed.

Figure Least Squares Approximation of Gridded Data illustrates least squares spline approximations of gridded river current data $u(x_1,x_2)$ based on these knots. In the left hand and center panels the data is represented as red spheres and the spline approximation as a continuous gray surface. The least squares approximation produced by AMA_MltvGrdLstsqr() is shown in the left hand panel, the least squares approximation produced by AMA_MltvLstsqr() is shown in the center panel, and a comparison of the two approximations is shown in the right hand panel where the blue surface represents the approximation computed by AMA_MltvLstsqr().

Least Squares Approximation of Gridded Data
Approximation Using Gridded Data Function Approximation Using Random Data Function Comparison of Approximations

From the right hand panel of Figure Least Squares Approximation of Gridded Data it can be seen that the two approximations are similar. It can be shown that the algorithms employed by AMA_MltvGrdLstsqr() and AMA_MltvLstsqr() are mathematically equivalent and, in the absence of bounds, they produce equivalent approximations. The difference illustrated in the right hand panel of the figure is attributed to the imposition of a lower bound. Therefore, for small data sets, either function can be used. However, for large gridded data sets the memory requirements of AMA_MltvLstsqr() probably exceeds the memory available on most computing systems. Since the function AMA_MltvGrdLstsqr() employs an algorithm which invokes the univariate function AMA_UnvLstsqr() multiple times, its memory requirements are significantly less and it can compute least squares spline approximations of extremely large multivariate gridded data sets.

The above example states, in the absence of bounds, that when the knots vectors ${\bf\Lambda}_k$ or ${\bf\Lambda}^o_k$, for $1\le k\le n$, are such that a least squares spline approximation of gridded data has an unique solution; then the functions AMA_MltvLstsqr() and AMA_MltvGrdLstsqr() produce equivalent approximations. However, this is only true if the penalty term option is not enabled. Consider the least squares approximations of the river current data $u(x_1,x_2)+E(x_1,x_2)$ based on the aforementioned knots with penalty term weight $\theta = 0.1$ which are illustrated in Figure Penalty Term Approximation of Gridded Data with Error. As before, the least squares approximation produced by AMA_MltvGrdLstsqr() is shown in the left hand panel, the least squares approximation produced by AMA_MltvLstsqr() is shown in the center panel, and a comparison of the two approximations is shown in the right hand panel where the blue surface represents the approximation computed by AMA_MltvLstsqr().

Penalty Term Approximation of Gridded Data with Error
Approximation Using Gridded Data Function Approximation Using Random Data Function Comparison of Approximations

The difference between the two approximations is due to AMA_MltvApprox() including the cross partial term

\[ F^{(2)}(s(x_1,x_2)) = 2.0\int_R\left(\partial^2 s({\bf X})\over\partial x_1 \partial x_2\right)^2dR \]

in its penalty term; which is omitted by AMA_MltvGrdLstsqr(), since it employs AMA_UnvLstsqr() to compute its approximation. The inclusion of the cross partial term does produce a smoother surface at the expense of longer computation times. For this example the differences are almost neglectable; but, in general, they can be more substantial.

Another situation where the functions AMA_MltvGrdLstsqr() and AMA_MltvLstsqr() produce different least squares spline approximations is when the knot vectors are such that an unique approximation does not exist and the minimum norm optimization is enabled. For instance, consider the knot vectors ${\bf\Lambda}^o_1 = (x_{1,1},x_{1,2},\cdots,x_{1,N^g_1-1},x_{1,N^g_1})$ and ${\bf\Lambda}^o_2 = (x_{2,1},x_{2,2},\cdots,x_{2,N^g_2-1},x_{2,N^g_2})$ which are equivalent to the independent variable data. It should also be noted that they are equivalent to the knots defined by AMA_MltvGrdApprox() and AMA_MltvGrdInterp() to approximate the gridded river current data with a cubic spline.

Based on the knot vectors defined above a cubic spline has $(N^g_1+2)(N^g_2+2)>N^g_1N^g_2$ spline coefficients. These knots provide sufficient spline coefficients to interpolate the data but, since the number of spline coefficients exceeds the number of data points, the least squares approximation of the gridded river current data is not unique. Therefore, the minimum norm optimization option must be enabled in order to insure an unique approximation. Figure Minimum Norm Approximation of Gridded Data illustrates least squares spline approximations of the gridded river current data $u(x_1,x_2)$ based on these knots with the minimum norm optimization enabled. As before, the least squares approximation produced by AMA_MltvGrdLstsqr() is shown in the left hand panel, the least squares approximation produced by AMA_MltvLstsqr() is shown in the center panel, and a comparison of the two approximations is shown in the right hand panel where the blue surface represents the approximation computed by AMA_MltvLstsqr().

Minimum Norm Approximation of Gridded Data
Approximation Using Gridded Data Function Approximation Using Random Data Function Comparison of Approximations

From the right hand panel of Figure Minimum Norm Approximation of Gridded Data it can be seen that the two approximations are similar in the interior of the data region where there is sufficient data to insure uniqueness of the spline coefficients. However, they differ in the boundary regions where extra degrees of freedom are caused by the inclusion of the second and next to last independent variable data points into the knot vectors. It is in these regions where the minimum norm optimization has its greatest influence. Since AMA_MltvLstsqr() includes the spline's cross partial terms in the minimum norm optimization's penalty term, it produces a different minimum norm approximation than AMA_MltvGrdLstsqr() within the boundary regions.

AMA_MltvGrdApprox() - Approximation of Multivariate Gridded Data

The approximation function AMA_MltvGrdApprox() requires the independent variable data $x_{k,\ell_k}$, for $\ell_k=1,\ldots,N^g_k$ and $1\le k\le n$; the dependent variable data $z_{\ell_1,\cdots,\ell_n}$ and the approximation tolerances $\epsilon_{\ell_1,\cdots,\ell_n}$, for $\ell_1=1,\ldots,N^g_1,\cdots,\ell_n=1,\ldots,N^g_n$. For this example $n=2$ and the independent variable data is the gridded data distribution illustrated in the right hand panel of Figure Data Distributions, the dependent variable data is the river current data $u(x_{1,\ell_1},x_{2,\ell_2})$ or the data $u(x_{1,\ell_1},x_{2,\ell_2})+E(x_{1,\ell_1},x_{2,\ell_2})$ where the approximation tolerances are either $\epsilon_{\ell_1,\ell_2}=0.0$ or $\epsilon_{\ell_1,\ell_2}=E(x_{1,\ell},x_{2,\ell})$, for $\ell=1,\ldots,31$ and $\ell_2=1,\ldots,31$. In general, a statistical analysis of the process which produced the data would provide nonzero approximation tolerances, but $E(x_1,x_2)$ serves to illustrate the influence of the approximation tolerances on the approximation produced by AMA_MltvGrdApprox().

Based on the data AMA_MltvGrdApprox() defines the knot vectors ${\bf\Lambda}_0$ and ${\bf\Lambda}_1$ which lie on the independent variable data's underlying rectilinear grid. They are illustrated in the right hand panel of Figure Tensor Product Grid and Data Distributions. These knot vectors define sufficient spline coefficients to insure satisfaction of the approximation constraints $z_{\ell_1,\ell_2} - \epsilon_{\ell_1,\ell_2} \le s(x_{1,\ell},x_{2,\ell}) \le z_{\ell_1,\ell_2} + \epsilon_{\ell_1,\ell_2}$, for $\ell_1=1,\ldots,31$ and $\ell_2=1,\ldots,31$.

Three approximations of the river current data whose independent variable data defines a gridded data distribution are illustrated in Figure Multivariate Gridded Approximation. The left hand and center panels illustrate approximations of $u(x_1,x_2)$ and $u(x_1,x_2)+E(x_1,x_2)$, respectively; with the approximation tolerances $\epsilon_{\ell_1,\ell_2}=0.0$ for all $\ell_1=1,\ldots,31$ and $\ell_2=1,\ldots,31$. The right hand panel shows an approximation of $u(x_1,x_2)+E(x_1,x_2)$ with the approximation tolerances $\epsilon_{\ell_1,\ell_2}=E(x_{1,\ell},x_{2,\ell})$, for $\ell_1=1,\ldots,31$ and $\ell_2=1,\ldots,31$. As before in each panel the data is given as red spheres and the approximation produced by AMA_MltvGrdApprox() is shown as a gray surface. Also, for comparison, the AMA_MltvApprox() approximation is shown as a blue surface.

Multivariate Gridded Approximation
Approximation of $u(x_1,x_2)$ with $\epsilon=0.0$ Approximation of $u(x_1,x_2)+E(x_1,x_2)$ with $\epsilon=0.0$ Approximation of $u(x_1,x_2)+E(x_1,x_2)$ with $\epsilon=E(x_1,x_2)$

The sum of squares and maximum errors for the three approximations are given in Table Multivariate Gridded Approximation Errors. They show the spline approximations of the river current data $u(x_1,x_2)$ and the data with a random error $u(x_1,x_2)+E(x_1,x_2)$ with $\epsilon=0.0$ are, as was the case with AMA_MltvApprox(), spline interpolants of the data. Also, comparing the errors for the approximation of $u(x_1,x_2)+E(x_1,x_2)$ with $\epsilon=E(x_1,x_2)$ with those provided for gridded data in Table River Current Errors indicates the approximation satisfies its approximation constraints.

Multivariate Gridded Approximation Errors
Data Approximated $u(x_1,x_2)$ with $\epsilon=0.0$ $u(x_1,x_2)+E(x_1,x_2)$ with $\epsilon=0.0$ $u(x_1,x_2)+E(x_1,x_2)$ with $\epsilon=E(x_1,x_2)$
Maximum Error $1.332\times 10^{-15}$ $1.776\times 10^{-15}$ $2.481\times 10^{-2}$
Sum of Squares Error $6.889\times 10^{-29}$ $1.065\times 10^{-29}$ $4.119\times 10^{-2}$

The multivariate gridded data approximation functions employ their univariate counterparts to solve a sequence of approximation problems along the data's underlying grid to produce a multivariate approximation. Since AMA_MltvGrdApprox() uses the univariate function AMA_UnvApprox(), it is not possible for it to include the second order cross partial terms in the penalty term. The effect of this is clearly seen in the left hand and center panels of Figure Multivariate Gridded Approximation. When the approximation tolerances are nonzero; as shown in the right hand panel of the figure, then the approximation constraints are imposed only while computing the univariate approximations in the first independent variable. This is another reason an AMA_MltvGrdApprox() approximation differs from an AMA_MltvApprox() approximation. It also causes the approximation to depend on the independent variable order.

AMA_MltvGrdInterp() - Interpolation of Multivariate Gridded Data

The interpolation function AMA_MltvGrdInterp() requires the independent variable data $x_{k,\ell_k}$, for $\ell_k=1,\ldots,N^g_k$ and $1\le k\le n$, and the dependent variable data $z_{\ell_1,\cdots,\ell_n}$, for $\ell_1=1,\ldots,N^g_1,\cdots,\ell_n=1,\ldots,N^g_n$. For this example $n=2$ and the independent variable data is the gridded data distribution illustrated in the right hand panel of Figure Data Distributions, the dependent variable data is based on the river current data $u(x_{1,\ell_1},x_{2,\ell_2})$, for $\ell=1,\ldots,31$ and $\ell_2=1,\ldots,31$.

Based on the data AMA_MltvGrdInterp() defines the knot vectors ${\bf\Lambda}_0$ and ${\bf\Lambda}_1$ which lie on the independent variable data's underlying rectilinear grid. They are illustrated in the right hand panel of Figure Tensor Product Grid and Data Distributions. These knot vectors define sufficient spline coefficients to satisfy the interpolation constraints $s(x_{1,\ell_1},x_{2,\ell_2}) = z_{\ell_1,\ell_2}$, for $1\le\ell_1\le 31$ and $1\le\ell_2\le 31$, and the boundary conditions ${\partial s^{2}(0.0,x_2)\over \partial x^2_1} = {\partial s^{2}(2\pi,x_2)\over \partial x^2_1} = 0.0$ and ${\partial s^{2}(x_1,-2.0)\over \partial x^2_2} = {\partial s^{2}(x_1,2.0)\over \partial x^2_2} = 0.0$, see AMA_MltvGrdInterp() for a description of the boundary conditions it imposes.

Two interpolants of the river current data whose independent variable data defines a gridded data distribution are illustrated in Figure Multivariate Gridded Interpolation. The left hand panel illustrates an interpolant of $u(x_1,x_2)$ and the right hand panel shows an interpolant of $u(x_1,x_2)+E(x_1,x_2)$. In both panels the interpolant produced by AMA_MltvGrdInterp() is shown as a gray surface and, for comparison, the AMA_MltvInterp() interpolant is shown as a blue surface.

Multivariate Gridded Interpolation
Interpolation of $u(x_1,x_2)$ Interpolation of $u(x_1,x_2)+E(x_1,x_2)$

Figure Multivariate Gridded Interpolation indicates the two interpolants produced by AMA_MltvGrdInterp() and AMA_MltvInterp() are equivalent. In fact, it can be shown that the algorithms employed by AMA_MltvGrdInterp() and AMA_MltvInterp() are mathematically equivalent and, when boundary conditions are imposed, the interpolant is unique and the two functions produce identical interpolants. Therefore, for small data sets, either function can be used. However, for large gridded data sets the memory requirements of AMA_MltvInterp() probably exceeds the memory available on most computing systems. Since the function AMA_MltvGrdInterp() employs an algorithm which invokes the univariate function AMA_UnvInterp() multiple times, its memory requirements are significantly less and it can compute least squares spline interpolants of extremely large multivariate gridded data sets.

With the exception of the boundary conditions the interpolation algorithm is equivalent to the approximation algorithm when the approximation tolerances are equal to zero. Therefore, if interpolation without boundary conditions is desired, then AMA_MltvGrdApprox() can be used with the approximation tolerances set equal to zero.

AMA_MltvGrdMonoApprox() - Monotonic Approximation of Multivariate Gridded Data

The monotonic approximation function AMA_MltvGrdMonoApprox() requires the independent variable data $x_{k,\ell_k}$, for $\ell_k=1,\ldots,N^g_k$ and $1\le k\le n$; the dependent variable data $z_{\ell_1,\cdots,\ell_n}$ and the approximation tolerances $\epsilon_{\ell_1,\cdots,\ell_n}$, for $\ell_1=1,\ldots,N^g_1,\cdots,\ell_n=1,\ldots,N^g_n$. For this example $n=2$ and the independent variable data is the gridded data distribution illustrated in the right hand panel of Figure Data Distributions, the dependent variable data is based on the river current data $u(x_{1,\ell_1},x_{2,\ell_2})$ and the approximation tolerances are either $\epsilon_{\ell_1,\ell_2}=0.0$ or $\epsilon_{\ell_1,\ell_2}=E(x_{1,\ell},x_{2,\ell})$, for $\ell=1,\ldots,31$ and $\ell_2=1,\ldots,31$.

Based on the data AMA_MltvGrdMonoApprox() defines knot vectors ${\bf\Lambda}_0$ and ${\bf\Lambda}_1$ which provide sufficient spline coefficients to satisfy the approximation constraints $z_{\ell_1,\ell_2} - \epsilon_{\ell_1,\ell_2} \le s(x_{1,\ell},x_{2,\ell}) \le z_{\ell_1,\ell_2} + \epsilon_{\ell_1,\ell_2}$, for $\ell_1=1,\ldots,31$ and $\ell_2=1,\ldots,31$, and the local monotonicity constraints; see AMA_MltvGrdMonoApprox() for a definition of the local monotonicity constraints it imposes.

Three monotonic approximations of the river current data whose independent variable data defines a gridded data distribution are illustrated in Figure Multivariate Gridded Monotonic Approximation. The left hand and center panels illustrate monotonic approximations of $u(x_1,x_2)$ and $u(x_1,x_2)+E(x_1,x_2)$, respectively; with the approximation tolerances $\epsilon_{\ell_1,\ell_2}=0.0$ for all $\ell_1=1,\ldots,31$ and $\ell_2=1,\ldots,31$. The right hand panel shows a monotonic approximation of $u(x_1,x_2)+E(x_1,x_2)$ with the approximation tolerances $\epsilon_{\ell_1,\ell_2}=E(x_{1,\ell},x_{2,\ell})$, for $\ell_1=1,\ldots,31$ and $\ell_2=1,\ldots,31$. In both panels the approximation produced by AMA_MltvGrdMonoApprox() is shown as a gray surface and, for comparison, the AMA_MltvMonoApprox() approximation is shown as a blue surface.

Multivariate Gridded Monotonic Approximation
Monotonic Approximation of $u(x_1,x_2)$ with $\epsilon=0.0$ Monotonic Approximation of $u(x_1,x_2)+E(x_1,x_2)$ with $\epsilon=0.0$ Monotonic Approximation of $u(x_1,x_2)+E(x_1,x_2)$ with $\epsilon=E(x_1,x_2)$

The sum of squares and maximum errors for the three approximations are given in Table Multivariate Gridded Monotonic Approximation Errors. They show the monotonic spline approximations of the river current data $u(x_1,x_2)$ and the data with a random error $u(x_1,x_2)+E(x_1,x_2)$ with $\epsilon=0.0$ are, as was the case with AMA_MltvMonoApprox(), spline interpolants of the data. Also, comparing the errors for the monotonic approximation of $u(x_1,x_2)+E(x_1,x_2)$ with $\epsilon=E(x_1,x_2)$ with those provided for gridded data in Table River Current Errors indicates the monotonic approximation satisfies its approximation constraints.

Multivariate Gridded Monotonic Approximation Errors
Data Approximated $u(x_1,x_2)$ with $\epsilon=0.0$ $u(x_1,x_2)+E(x_1,x_2)$ with $\epsilon=0.0$ $u(x_1,x_2)+E(x_1,x_2)$ with $\epsilon=E(x_1,x_2)$
Maximum Error $1.332\times 10^{-13}$ $1.425\times 10^{-13}$ $2.481\times 10^{-2}$
Sum of Squares Error $2.130\times 10^{-25}$ $1.420\times 10^{-24}$ $4.385\times 10^{-2}$

Because AMA_MltvGrdMonoApprox() produces an interpolant by invoking AMA_UnvMonoApprox() multiple times to compute univariate monotonic approximations along the underlying data grid, the monotonicity constraints are imposed separately along the data grid; instead of collectively, as in they are by AMA_MltvMonoApprox(). This, along with the absence of the cross partials from the penalty term, causes the differences between the approximations illustrated in the left hand and center panels of the figure. When the approximation tolerances are nonzero; as shown in the right hand panel of the figure, then the approximation constraints are imposed only while computing the univariate approximations in the first independent variable. This is another reason an AMA_MltvGrdMonoApprox() approximation differs from an AMA_MltvMonoApprox() approximation. It also causes the approximation to depend on the independent variable order.

AMA_MltvGrdMonoInterp() - Monotonic Interpolation of Multivariate Gridded Data

The monotonic interpolation function AMA_MltvGrdMonoInterp() requires the independent variable data $x_{k,\ell_k}$, for $\ell_k=1,\ldots,N^g_k$ and $1\le k\le n$, and the dependent variable data $z_{\ell_1,\cdots,\ell_n}$, for $\ell_1=1,\ldots,N^g_1,\cdots,\ell_n=1,\ldots,N^g_n$. For this example $n=2$ and the independent variable data is the gridded data distribution illustrated in the right hand panel of Figure Data Distributions and the dependent variable data is based on the river current data $u(x_{1,\ell_1},x_{2,\ell_2})$, for $\ell=1,\ldots,31$ and $\ell_2=1,\ldots,31$.

Based on the data AMA_MltvGrdMonoInterp() defines knot vectors ${\bf\Lambda}_0$ and ${\bf\Lambda}_1$ which provide sufficient spline coefficients to satisfy the interpolation constraints $s(x_{1,\ell_1},x_{2,\ell_2}) = z_{\ell_1,\ell_2}$, for $\ell_1=1,\ldots, 31$ and $\ell_2=1,\ldots, 31$, and the local monotonicity constraints, see AMA_MltvGrdMonoInterp() for a definition of the local monotonicity constraints it imposes.

Two monotonic interpolants of the river current data whose independent variable data defines a gridded data distribution are illustrated in Figure Multivariate Gridded Monotonic Interpolation. The monotonic interpolant for $u(x_1,x_2)$ is illustrated in the left hand panel and the monotonic interpolant for $u(x_1,x_2)+E(x_1,x_2)$ is shown in the right hand panel. In both panels the interpolant produced by AMA_MltvGrdMonoInterp() is shown as a gray surface and, for comparison, the AMA_MltvMonoInterp() interpolant is shown as a blue surface.

Multivariate Gridded Monotonic Interpolation
Monotonic Interpolation of $u(x_1,x_2)$ Monotonic Interpolation of $u(x_1,x_2)+E(x_1,x_2)$

It was shown in Section AMA_MltvGrdInterp() - Interpolation of Multivariate Gridded Data that when boundary conditions are imposed the interpolants produced by AMA_MltvInterp() and AMA_MltvGrdInterp() are equivalent. However, since AMA_MltvGrdMonoInterp() produces an interpolant by invoking AMA_UnvMonoInterp() multiple times to compute univariate monotonic interpolants along the underlying data grid, the monotonicity constraints are imposed separately along the data grid; instead of collectively, as in they are by AMA_MltvMonoInterp(). This, along with the absence of the cross partials from the penalty term, causes the differences between the interpolants illustrated above.

Since AMA_MltvGrdMonoInterp() does not impose boundary conditions on its approximation, the approximation it produces is identical to the approximation produced by AMA_MltvGrdMonoApprox() when its approximation tolerances are set equal to zero. Hence, AMA_MltvGrdMonoInterp() is redundant but, nonetheless, is provided for convenience.