#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 . Must satisfy $2\le$ nind $\le$ AMA_MXNIND.
ng	[in] Array of size nind containing the number of points where ng[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 where degree[k] . Must satisfy $1\le$ degree[k] $\le 5$ .
mlamda	[in] Array of size nind containing the number of knots where mlamda[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 .
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

Functions

Function Documentation