Functions
AMA_MltvGrdLstsqr.c File Reference
#include <AMA.h>

Functions

long int AMA_MltvGrdLstsqr (AMA_OPTIONS *options, long int nind, long int *ng, double **x, double *z, double *wht, long int *degree, long int *mlamda, double **lamda, double theta, AMA_SPLINE **spline)
 Least Squares Approximation of Multivariate Gridded Data More...
 

Function Documentation

long int AMA_MltvGrdLstsqr ( AMA_OPTIONS options,
long int  nind,
long int *  ng,
double **  x,
double *  z,
double *  wht,
long int *  degree,
long int *  mlamda,
double **  lamda,
double  theta,
AMA_SPLINE **  spline 
)

Least Squares Approximation of Multivariate Gridded Data

This function employs cnspla to compute a least squares spline approximation of independent variable data ${\bf X}\in{\rm R}^n$ and dependent variable data ${\bf Z}\in{\rm R}^1$. 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.

For a given set of data this function employs cnspla to compute the spline

\[ 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}) \]

that minimizes

\[ (1.0-\theta)\sum_{\ell_n=1}^{N^g_n}\cdots\sum_{\ell_1=1}^{N^g_1} w_{\ell_1,\cdots,\ell_n}( s(x_{1,\ell_1},\cdots,x_{n,\ell_n}) - z_{\ell_1,\cdots,\ell_n} )^2 + \sum_{k=1}^n\left[\int_{\lambda_{k,1}}^{\lambda_{k,m_k+d_k+1}}\left(\partial^2 s({\bf X})\over\partial^2 x_k\right)^2dx_k\right] \]

for $0.0\le\theta< 1.0$.

In the above definition of $s({\bf X})$ the $\alpha_{j_1,\ldots,j_n}$, for $j_k=1,\ldots,m_k$ and $1\le k\le n$, 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 $B_{d_k,j_k}(x_k|{\bf\Lambda}_k)$ are the $m_k$ univariate B-splines of degree $1 \le d_k\le 5$ defined by the knot vector ${\bf\Lambda}_k\in{\rm R}^{M_k}$. The knot vector ${\bf\Lambda}_k$ is given as

\[ {\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 \]

where the knots $\lambda_{k,i}$ must satisfy the conditions

  • $\lambda_{k,1} = \ldots = \lambda_{k,d_k+1}$;
  • $\lambda_{k,i} \le \lambda_{k,i+1}$, for $1\le i\le m_k+d_k$;
  • $\lambda_{k,i} < \lambda_{i+d_k+1}$, for $1\le i\le m_k$;
  • $\lambda_{k,m_k+1} = \ldots = \lambda_{k,m_k+d_k+1}$.

If the knot vector ${\bf\Lambda}_k$ is represented by its distinct knot vector

\[ {\rm\bf K}^{{\bf\Lambda}_k} = ( \kappa_{k,1}, \kappa_{k,2}, \cdots, \kappa_{k,m^d_k-1}, \kappa_{k,m^d_k} )^T \]

and its knot multiplicity vector

\[ {\rm\bf H}^{{\bf\Lambda}_k} = ( \eta_{k,1}, \eta_{k,2}, \cdots, \eta_{k,m^d_k-1}, \eta_{k,m^d_k} )^T \]

then the aforementioned conditions are

  • $\kappa_{k,i} < \kappa_{k,i+1}$, for $1\le i\le m^d_k - 1$;
  • $\eta_{k,i} \le d_k + 1$, for $2\le i\le m^d_k-1$;
  • $\eta_{k,1} = \eta_{k,m^d_k} = d_k+1$.

Also the first and last distinct knots must satisfy the conditions $\kappa_{k,1}\le x_{k,1}$ and $\kappa_{k,m^d_k}\ge x_{k,n^g_k}$, respectively. In this case the number of knots is $M_k = m_k + d_k + 1$.

For convenience this function does not require the definition of the $d_k+1$ order knots at the left and right hand endpoints and also accepts the knot vector ${\bf\Lambda}^o_k\in{\rm R}^{M_k}$ given as

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

where the knots $\lambda_{k,i}$ satisfy the conditions

  • $\lambda_{k,i} \le \lambda_{k,i+1}$, for $d_k+2\le i\le m_k-1$;
  • $\lambda_{k,i} < \lambda_{i+d_k+1}$, for $d_k+1\le i\le m_k-d_k$;
  • $\lambda_{k,d_k+1} < \lambda_{k,d_k+2}$ and $\lambda_{k,m_k} < \lambda_{k,m_k+1}$.

If the knot vector ${\bf\Lambda}^o$ is represented by its distinct knot vector

\[ {\rm\bf K}^{{\bf\Lambda}^o} = ( \kappa_{k,1}, \kappa_{k,2}, \cdots, \kappa_{k,m^d_k-1}, \kappa_{k,m^d_k} )^T \]

and its knot multiplicity vector

\[ {\rm\bf H}^{{\bf\Lambda}^o} = ( \eta_{k,1}, \eta_{k,2}, \cdots, \eta_{k,m^d_k-1}, \eta_{k,m^d_k} )^T \]

then the aforementioned conditions are

  • $\kappa_{k,i} < \kappa_{k,i+1}$, for $1\le i\le m^d_k - 1$;
  • $\eta_{k,i} \le d_k + 1$, for $2\le i\le m^d_k-1$;
  • $\eta_{k,1} = \eta_{k,m^d_k} = 1$.

As before, the first and last distinct knots must satisfy the conditions $\kappa_{k,1}\le x_{k,1}$ and $\kappa_{k,m^d_k}\ge x_{k,n^g_k}$, respectively. In this case the number of knots is $M_k = m_k - d_k + 1$.

This function employs an efficient algorithm for computing least squares spline approximations of multivariate gridded data. The algorithm consists of computing the splines

\[ \bar s_k({\bf X}) = \sum_{q_n=1}^{\bar m_n}\cdots\sum_{q_1=1}^{\bar m_1}\bar\alpha^k_{q_1,\ldots,q_n}B_{\bar d_1,q_1}(x_1|{\bf\bar\Lambda}_1)\cdots B_{\bar d_n,q_n}(x_n|{\bf\bar\Lambda_n}) \]

for $1\le k\le n$. The coefficients $\bar\alpha^k_{q_1,\ldots,q_n}$ of $\bar s_k({\bf X})$ are the coefficients of the tensor product B-splines

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

where $B_{\bar d_k,q_k}(x_k|{\bf\bar\Lambda}_k)$ are the $\bar m_k$ univariate B-splines of degree $\bar d_k$ defined by the knot vector ${\bf\bar\Lambda}_k\in{\rm R}^{\bar m_k + \bar d_k + 1}$. The degree $\bar d_{k_o}$ is

\[ \bar d_{k_o} = \cases { d_{k_o}, & for $1\le k_o\le k$, \cr 1 , & for $k < k_o\le n $; \cr } \]

the number of coefficients $\bar m_{k_o}$ is

\[ \bar m_{k_o} = \cases { m_{k_o} , & for $1\le k_o\le k$, \cr n^g_{k_o}, & for $k < k_o\le n $; \cr } \]

and the knot vector ${\bf\bar\Lambda}_{k_o}$ is

\[ {\bf\bar\Lambda}_{k_o} = \cases { {\bf\Lambda}_{k_o} , & for $1\le k_o\le k$, \cr ( x_{k_o,1}, x_{k_o,1}, x_{k_o,2}, \cdots, x_{k_o,n^g_{k_o}-1}, x_{k_o,n^g_{k_o}}, x_{k_o,n^g_{k_o}})^T & for $k < k_o\le n $. \cr } \]

The spline $s({\bf X})$ is based on the spline $\bar s_n({\bf X})$ and is defined by setting its coefficents $\alpha_{j_1,\ldots,j_n}=\bar\alpha^n_{j_1,\ldots,j_n}$ for $j_k=1,\ldots,m_k$ and $1\le k\le n$. The spline $\bar s_n({\bf X})$ is defined through multiple invocations of AMA_UnvLstsqr() in the following manner.

First, the spline $\bar s_1({\bf X})$ is defined by invoking AMA_UnvLstsqr() $\prod_{k=2}^nN^g_{k}$ times to compute the univariate spline

\[ \bar s_{q_2,\cdots,q_n}(x_1) = \sum_{q_1=1}^{m_1}\bar\alpha^1_{q_1,\ldots,q_n}B_{d_1,q_1}(x_1|{\bf\Lambda}_1) \]

that minimizes

\[ (1.0-\theta) \sum_{\ell_1=1}^{N^g_1} w_{\ell_1,q_2,\cdots,q_n} ( \bar s_{q_2,\cdots,q_n}(x_{1,\ell_1}) - z_{\ell_1,q_2,\cdots,q_n} )^2 + \theta \int_{\lambda_{1,1}}^{\lambda_{1,m_1+d_1+1}}\left( \bar s^{(2)}_{q_2,\cdots,q_n}(x_1) \right)^2 dx_1. \]

Computing the univariate splines $\bar s_{q_2,\cdots,q_n}(x_1)$, for $q_k=1,\ldots,N^g_k$ and $2\le k\le n$, defines the coefficients $\bar\alpha^1_{q_1,\ldots,q_n}$ of $\bar s_1({\bf X})$.

After computing $\bar s_1({\bf X})$ the splines $\bar s_k({\bf X})$, for $2\le k\le n$, are defined by invoking AMA_UnvLstsqr() $\prod_{k_o=1}^{k-1}m_{k_o}\prod_{k_o=k+1}^nN^g_{k_o}$ times to compute the univariate spline

\[ \bar s_{q_1,\cdots,q_{k-1},q_{k+1},\cdots,q_n}(x_k) = \sum_{q_k=1}^{m_k}\bar\alpha^k_{q_1,\ldots,q_n}B_{d_k,q_k}(x_k|{\bf\Lambda}_k) \]

that minimizes

\[ (1.0-\theta) \sum_{\ell_k=1}^{N^g_k} ( \bar s_{q_1,\cdots,q_{k-1},q_{k+1},\cdots,q_n}(x_{k,\ell_k}) - \bar\alpha^{k-1}_{q_1,\cdots,q_{k-1},\ell_k,q_{k+1},\cdots,q_n} )^2 + \theta \int_{\lambda_{k,1}}^{\lambda_{k,m_k+d_k+1}}\left( \bar s^{(2)}_{q_1,\cdots,q_{k-1},q_{k+1},\cdots,q_n}(x_k) \right)^2 dx_k. \]

Computing the univariate splines $\bar s_{q_1,\cdots,q_{k-1},q_{k+1},\cdots,q_n}(x_k)$, for $q_{k_o}=1,\ldots,m_{k_o}$ when $1\le k_o\le k-1$ and for $q_{k_o}=1,\ldots,N^g_{k_o}$ when $k+1 \le k_o \le n$, defines the coefficients $\bar\alpha^k_{q_1,\ldots,q_n}$ of $\bar s_k({\bf X})$.

In total this algorithm consists of $\sum_{k=1}^n\left(\prod_{k_o=1}^{k-1}m_{k_o}\prod_{k_o=k+1}^nN^g_{k_o}\right)$ invocations of AMA_UnvLstsqr() to compute a least squares spline approximation.

This function does the following:

  • Checks the input parameters for validity, see AMA_MltvGrdInputCheck().
  • Defines the spline $s({\bf X})$ based on the knot vectors ${\bf\Lambda}_k$, for $1\le k\le n$.
  • Defines the splines $\bar s_k({\bf X})$, for $1\le k\le n$, through multiple invocations of AMA_UnvLstsqr().
  • Stores the least squares approximation $s({\bf X})=\bar s_n({\bf X})$ into spline.
Note
By default the spline coefficients are initialized to zero but a different initial value for the spline coefficients can be set with AMA_OptionsSetCoefficients().
By default the minimum norm optimization is not performed but it can be enabled with AMA_OptionsSetMinimumNorm().
Because this function employs the univariate spline approximation function AMA_UnvLstsqr(), it does not minimizes the cross partial terms of the penalty term utilized by AMA_MltvLstsqr().
If the minimum norm optimization option is enabled with AMA_OptionsSetMinimumNorm(), then it is performed during each invocation of AMA_UnvLstsqr().
The penalty term and minimum norm optimization options are mutually exclusive; that is, the minimum norm optimization is not performed if $\theta > 0.0$.
Since the weights are only used during the computation of $\bar s_1({\bf X})$, the least squares approximation depends on the order of the independent variables.

Parameter Note: In the parameters description given below the limits on $k$ are $1\le k\le n$, k = $k-1$ and N $=\prod_{k=1}^nN^g_k$.

Parameters
options[in] Pointer to AMA_OPTIONS. Must be initialized with AMA_Options() prior to calling AMA_MltvGrdLstsqr().
nind[in] The number of independent variables $n$. Must satisfy $2\le$ nind $\le$ AMA_MXNIND.
ng[in] Array of size nind containing the number of points $N^g_k$ where ng[k] $= N^g_k$. Must satisfy ng[k] $\ge 2$.
x[in] Array of size nind containing arrays of size ng[k] where x[k] contains the independent variable data $x_{k,\ell_k}$, for $\ell_k=1,\ldots,N^g_k$. The values of x[k] must be in ascending order.
z[in] Array of size N containing 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$.
wht[in] Array of size N containing the weights $w_{\ell_1,\cdots,\ell_n}$, for $\ell_1=1,\ldots,N^g_1,\cdots,\ell_n=1,\ldots,N^g_n$. Must satisfy $w_{\ell_1,\cdots,\ell_n}\ge 0.0$ for all $\ell_1=1,\ldots,N^g_1,\cdots,\ell_n=1,\ldots,N^g_n$. If wht = NULL, then this function uses $w_{\ell_1,\cdots,\ell_n}=1.0$.
degree[in] Array of size nind containing the degree $d_k$ where degree[k] $= d_k$. Must satisfy $1\le$ degree[k] $\le 5$.
mlamda[in] Array of size nind containing the number of knots $M_k$ where mlamda[k] $= M_k$. Must satisfy mlamda[k] $\ge 2$.
lamda[in] Array of size nind containing arrays of size mlamda[k] where lamda[k] contains the knot vector ${\bf\Lambda}_k$ or ${\bf\Lambda}^o_k$.
theta[in] The penalty term weight $\theta$. Must satisfy $0.0\le$ theta $< 1.0$.
spline[out] Pointer to AMA_SPLINE pointer containing the least squares spline approximation. Must satisfy spline $\ne$ NULL.
Returns
Success/Error Code.

User Callable Function - Documented 103015