NAG Library Function Document

nag_2d_spline_fit_ts_scat (e02jdc)

Note: this function uses optional arguments to define choices in the problem specification and in the details of the algorithm. If you wish to use default settings for all of the optional arguments, you need only read Sections 1 to 10 of this document. If, however, you wish to reset some or all of the settings please refer to Section 11 for a detailed description of the specification of the optional arguments produced by the function.

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

+− 10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

+− 11 Optional Arguments

11.1 Description of the Optional Arguments

1 Purpose

nag_2d_spline_fit_ts_scat (e02jdc) computes a spline approximation to a set of scattered data using a two-stage approximation method.

The computational complexity of the method grows linearly with the number of data points; hence large datasets are easily accommodated.

2 Specification

#include <nag.h>

#include <nage02.h>

void	nag_2d_spline_fit_ts_scat (Integer n, const double x[], const double y[], const double f[], Integer lsminp, Integer lsmaxp, Integer nxcels, Integer nycels, Integer lcoefs, double coefs[], Integer iopts[], double opts[], NagError *fail)

Before calling nag_2d_spline_fit_ts_scat (e02jdc), nag_fit_opt_set (e02zkc) must be called with optstr set to "Initialize = e02jdc". Settings for optional algorithmic arguments may be specified by calling nag_fit_opt_set (e02zkc) before a call to nag_2d_spline_fit_ts_scat (e02jdc).

3 Description

nag_2d_spline_fit_ts_scat (e02jdc) determines a smooth bivariate spline approximation to a set of data points

(x_{i}, y_{i}, f_{i})

, for

i = 1, 2, \dots, n

. Here, ‘smooth’ means

C^{1}

The approximation domain is the bounding box

[x_{\min}, x_{\max}] \times [y_{\min}, y_{\max}]

, where

x_{\min}

(respectively

y_{\min}

) and

x_{\max}

(respectively

y_{\max}

) denote the lowest and highest data values of the

(x_{i})

(respectively

(y_{i})

The spline is computed by local approximations on a uniform triangulation of the bounding box. These approximations are extended to a smooth spline representation of the surface over the domain. The local approximation scheme is by least squares polynomials (Davydov and Zeilfelder (2004)).

The two-stage approximation method employed by nag_2d_spline_fit_ts_scat (e02jdc) is derived from the TSFIT package of O. Davydov and F. Zeilfelder.

Values of the computed spline can subsequently be computed by calling nag_2d_spline_ts_eval (e02jec) or nag_2d_spline_ts_eval_rect (e02jfc).

4 References

Davydov O and Zeilfelder F (2004) Scattered data fitting by direct extension of local polynomials to bivariate splines Advances in Comp. Math. 21 223–271

5 Arguments

1: n – IntegerInput

On entry:

n

, the number of data values to be fitted.

Constraint:

n > 1

2: x[n] – const doubleInput

3: y[n] – const doubleInput

4: f[n] – const doubleInput

On entry: the

(x_{i}, y_{i}, f_{i})

data values to be fitted.

Constraint:

x [j - 1] \neq x [0]

for some

j = 2, \dots, n

and

y [k - 1] \neq y [0]

for some

k = 2, \dots, n

; i.e., there are at least two distinct

x

and

y

values.

5: lsminp – IntegerInput

6: lsmaxp – IntegerInput

On entry: are control parameters for the local approximations.

Each local approximation is computed on a local domain containing one of the triangles in the discretization of the bounding box. The size of each local domain will be adaptively chosen such that if it contains fewer than lsminp sample points it is expanded, else if it contains greater than lsmaxp sample points a thinning method is applied. lsmaxp mainly controls computational cost (in that working with a thinned set of points is cheaper and may be appropriate if the input data is densely distributed), while lsminp allows handling of different types of scattered data.

Setting

lsmaxp < lsminp

, and therefore forcing either expansion or thinning, may be useful for computing initial coarse approximations. In general smaller values for these arguments reduces cost.

A calibration procedure (experimenting with a small subset of the data to be fitted and validating the results) may be needed to choose the most appropriate values for lsminp and lsmaxp.

Constraints:

$1 \leq lsminp \leq n$ ;
$lsmaxp \geq 1$ .

7: nxcels – IntegerInput

8: nycels – IntegerInput

On entry: nxcels (respectively nycels) is the number of cells in the

x

(respectively

y

) direction that will be used to create the triangulation of the bounding box of the domain of the function to be fitted.

Greater efficiency generally comes when nxcels and nycels are chosen to be of the same order of magnitude and are such that n is

O (nxcels \times nycels)

. Thus for a ‘square’ triangulation — when

nxcels = nycels

— the quantities

\sqrt{n}

and nxcels should be of the same order of magnitude. See also Section 9.

Constraints:

$nxcels \geq 1$ ;
$nycels \geq 1$ .

9: lcoefs – IntegerInput

10: coefs[lcoefs] – doubleOutput

On exit: if

fail . code =

NE_NOERROR on exit, coefs contains the computed spline coefficients.

Constraint:

lcoefs \geq (((nxcels + 2) \times (nycels + 2) + 1) / 2) \times 10 + 1

11: iopts[ $\dim$ ] – IntegerCommunication Array

Note: the dimension,

\dim

, of this array is dictated by the requirements of associated functions that must have been previously called. This array MUST be the same array passed as argument iopts in the previous call to nag_fit_opt_set (e02zkc).

On entry: the contents of iopts MUST NOT be modified in any way either directly or indirectly, by further calls to nag_fit_opt_set (e02zkc), before calling either or both of the evaluation routines nag_2d_spline_ts_eval (e02jec) and nag_2d_spline_ts_eval_rect (e02jfc).

12: opts[ $\dim$ ] – doubleCommunication Array

Note: the dimension,

\dim

, of this array is dictated by the requirements of associated functions that must have been previously called. This array MUST be the same array passed as argument opts in the previous call to nag_fit_opt_set (e02zkc).

On entry: the contents of opts MUST NOT be modified in any way either directly or indirectly, by further calls to nag_fit_opt_set (e02zkc), before calling either or both of the evaluation routines nag_2d_spline_ts_eval (e02jec) and nag_2d_spline_ts_eval_rect (e02jfc).

13: fail – NagError *Input/Output

The NAG error argument (see Section 3.6 in the Essential Introduction).

6 Error Indicators and Warnings

NE_ALL_ELEMENTS_EQUAL: On entry, all elements of x or of y are equal.
NE_ALLOC_FAIL: Dynamic memory allocation failed.
NE_BAD_PARAM: On entry, argument $⟨value⟩$ had an illegal value.
NE_INITIALIZATION: Option arrays are not initialized or are corrupted.
NE_INT: On entry, $lcoefs = ⟨value⟩$ .
Constraint: $lcoefs \geq (((nxcels + 2) \times (nycels + 2) + 1) / 2) \times 10 + 1$ .
On entry, $lsmaxp = ⟨value⟩$ .
Constraint: $lsmaxp \geq 1$ .
On entry, $n = ⟨value⟩$ .
Constraint: $n > 1$ .
On entry, $nxcels = ⟨value⟩$ .
Constraint: $nxcels \geq 1$ .
On entry, $nycels = ⟨value⟩$ .
Constraint: $nycels \geq 1$ .
NE_INT_2: On entry, $lsminp = ⟨value⟩$ and $n = ⟨value⟩$ .
Constraint: $1 \leq lsminp \leq n$ .
NE_INTERNAL_ERROR: An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
An unexpected algorithmic failure was encountered. Please contact NAG.
NE_INVALID_OPTION: The value of optional argument $Polynomial Starting Degree$ was invalid.

7 Accuracy

Technical results on error bounds can be found in Davydov and Zeilfelder (2004).

Local approximation by polynomials of degree

d

for

n

data points has optimal approximation order

n^{- (d + 1) / 2}

The approximation error for

C^{1}

global smoothing is

O (n^{- 2})

Whether maximal accuracy is achieved depends on the distribution of the input data and the choices of the algorithmic parameters. The reference above contains extensive numerical tests and further technical discussions of how best to configure the method.

8 Parallelism and Performance

nag_2d_spline_fit_ts_scat (e02jdc) is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.

nag_2d_spline_fit_ts_scat (e02jdc) makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.

Please consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

n

-linear complexity and memory usage can be attained for sufficiently dense input data if the triangulation parameters nxcels and nycels are chosen as recommended in their descriptions above. For sparse input data on such triangulations, if many expansion steps are required (see lsminp) the complexity may rise to be loglinear.

10 Example

The Franke function

\begin{matrix} f (x, y) & = & 0.75 \exp (- ({(9 x - 2)}^{2} + {(9 y - 2)}^{2}) / 4) + \\ 0.75 \exp (- {(9 x + 1)}^{2} / 49 - (9 y + 1) / 10) + \\ 0.5 \exp (- ({(9 x - 7)}^{2} + {(9 y - 3)}^{2}) / 4) - \\ 0.2 \exp (- {(9 x - 4)}^{2} - {(9 y - 7)}^{2}) \end{matrix}

is widely used for testing surface-fitting methods. The example program randomly generates a number of points on this surface. From these a spline is computed and then evaluated at a vector of points and on a mesh.

11 Optional Arguments

Several optional arguments in nag_2d_spline_fit_ts_scat (e02jdc) control aspects of the algorithm, methodology used, logic or output. Their values are contained in the arrays iopts and opts; these must be initialized before calling nag_2d_spline_fit_ts_scat (e02jdc) by first calling nag_fit_opt_set (e02zkc) with optstr set to "Initialize = e02jdc".

Each optional argument has an associated default value; to set any of them to a non-default value, or to reset any of them to the default value, use nag_fit_opt_set (e02zkc). The current value of an optional argument can be queried using nag_fit_opt_get (e02zlc).

The remainder of this section can be skipped if you wish to use the default values for all optional arguments.

The following is a list of the optional arguments available. A full description of each optional argument is provided in Section 11.1.

11.1 Description of the Optional Arguments

For each option, we give a summary line, a description of the optional argument and details of constraints.

The summary line contains:

the keywords;
a parameter value, where the letters $a$ , $i and r$ denote options that take character, integer and real values respectively;
the default value.

Keywords and character values are case insensitive.

For nag_2d_spline_fit_ts_scat (e02jdc) the maximum length of the parameter cvalue used by nag_fit_opt_get (e02zlc) is

6

Averaged Spline

a

Default

='NO'

When the bounding box is triangulated there are 8 equivalent configurations of the mesh. Setting

Averaged Spline ='YES'

will use the averaged value of the

8

possible local polynomial approximations over each triangle in the mesh. This usually gives better results but at (about 8 times) higher computational cost.

Constraint:

Averaged Spline ='YES'

'NO'

Minimum Singular Value LPA

r

Default

= 1.0

A tolerance measure for accepting or rejecting a local polynomial approximation (LPA) as reliable.

The solution of a local least squares problem solved on each triangle subdomain is accepted as reliable if the minimum singular value

σ

of the matrix (of Bernstein polynomial values) associated with the least squares problem satisfies

Minimum Singular Value LPA \leq σ

In general the approximation power will be reduced as

Minimum Singular Value LPA

is reduced. (A small

σ

indicates that the local data has hidden redundancies which prevent it from carrying enough information for a good approximation to be made.) Setting

Minimum Singular Value LPA

very large may have the detrimental effect that only approximations of low degree are deemed reliable.

Minimum Singular Value LPA

will have no effect if

Polynomial Starting Degree = 0

, and it will have little effect if the input data is ‘smooth’ (e.g., from a known function).

A calibration procedure (experimenting with a small subset of the data to be fitted and validating the results) may be needed to choose the most appropriate value for this parameter.

Constraint:

Minimum Singular Value LPA \geq 0.0

Polynomial Starting Degree

i

Default

= 1

The degree to be used in the initial step of each local polynomial approximation.

At the initial step the method will attempt to fit with local polynomials of degree

Polynomial Starting Degree

. If the approximation is deemed unreliable (according to

Minimum Singular Value LPA

), the degree will be decremented by one and a new local approximation computed, ending with a constant approximation if no other is reliable.

Polynomial Starting Degree

is bounded from above by the maximum possible spline degree,

3

The default value gives a good compromise between efficiency and accuracy. In general the best approximation can be obtained by setting

Polynomial Starting Degree = 3

Constraint:

0 \leq Polynomial Starting Degree \leq 3

nag_2d_spline_fit_ts_scat (e02jdc) (PDF version)

e02 Chapter Contents

e02 Chapter Introduction

NAG Library Manual

NAG Library Function Documentnag_2d_spline_fit_ts_scat (e02jdc)

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

11 Optional Arguments

11.1 Description of the Optional Arguments

NAG Library Function Document

nag_2d_spline_fit_ts_scat (e02jdc)