1: $t (:)$ – double array: The dimension of the array t must be at least $\max (1, n)$

$t (i + 1)$ must contain the value $τ_{i}$ , for $i = 0, 1, \dots, n - 1$ .

Constraint: $t (1) > 0.0$ . Note that if this is not true, then the Toeplitz matrix cannot be positive definite.
2: $b (:)$ – double array: The dimension of the array b must be at least $\max (1, n)$

The right-hand side vector $b_{n}$ .
3: $x (:)$ – double array: The dimension of the array x must be at least $\max (1, n)$

With $n > 1$ the ( $n - 1$ ) elements of the solution vector $x_{n - 1}$ as returned by a previous call to nag_linsys_real_toeplitz_update (f04mf). The element $x (n)$ need not be specified.
4: $work (:)$ – double array: The dimension of the array work must be at least $\max (1, 2 \times n - 1)$

With $n > 2$ the elements of work should be as returned from a previous call to nag_linsys_real_toeplitz_update (f04mf) with ( $n - 1$ ) as the argument n.

Optional Input Parameters

1: $n$ – int64int32nag_int scalar: Default: the dimension of the arrays t, b, x.
The order of the Toeplitz matrix $T$ .

Constraint: $n \geq 0$ . When $n = 0$ , then an immediate return is effected.

Output Parameters

1: $x (:)$ – double array

The dimension of the array x will be

\max (1, n)

The solution vector

x_{n}

2: $p$ – double scalar

The reflection coefficient

p_{n - 1}

. (See Further Comments.)

3: $work (:)$ – double array

The dimension of the array work will be

\max (1, 2 \times n - 1)

The first (

n - 1

) elements of work contain the solution to the Yule–Walker equations

T_{n - 1} y_{n - 1} = - t_{n - 1},

where

t_{n - 1} = {(τ_{1} τ_{2} \dots τ_{n - 1})}^{t}

4: $ifail$ – int64int32nag_int scalar

ifail = 0

unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

$ifail = - 1$

On entry,	$n < 0$ ,
or	$t (0) \leq 0.0$ .

W $ifail = 1$: The Toeplitz matrix $T_{n}$ is not positive definite to working accuracy. If, on exit, p is close to unity, then $T_{n}$ was probably close to being singular.

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.

$ifail = - 999$: Dynamic memory allocation failed.

Accuracy

The computed solution of the equations certainly satisfies

r = T_{n} x_{n} - b_{n},

where

{‖r‖}_{1}

is approximately bounded by

{‖r‖}_{1} \leq c ε C (T_{n}),

c

being a modest function of

n

ε

being the machine precision and

C (T)

being the condition number of

T

with respect to inversion. This bound is almost certainly pessimistic, but it seems unlikely that the method of Levinson is backward stable, so caution should be exercised when

T_{n}

is ill-conditioned. The following bound on

T_{n}^{- 1}

holds:

\max (\frac{1}{\prod_{i = 1}^{n - 1} (1 - p_{i}^{2})}, \frac{1}{\prod_{i = 1}^{n - 1} (1 - p_{i})}) \leq {‖T_{n}^{- 1}‖}_{1} \leq \prod_{i = 1}^{n - 1} (\frac{1 + |p_{i}|}{1 - |p_{i}|}) .

(See Golub and Van Loan (1996).) The norm of

T_{n}^{- 1}

may also be estimated using function nag_linsys_real_gen_norm_rcomm (f04yd). For further information on stability issues see Bunch (1985), Bunch (1987), Cybenko (1980) and Golub and Van Loan (1996).

Further Comments

The number of floating-point operations used by this function is approximately

8 n

y_{i}

is the solution of the equations

T_{i} y_{i} = - {(τ_{1} τ_{2} \dots τ_{i})}^{T},

then the reflection coefficient

p_{i}

is defined as the

i

th element of

y_{i}

Example

This example finds the solution of the equations

T_{k} x_{k} = b_{k}

k = 1, 2, 3, 4

, where

T_{4} = (\begin{array}{l} 4 & 3 & 2 & 1 \\ 3 & 4 & 3 & 2 \\ 2 & 3 & 4 & 3 \\ 1 & 2 & 3 & 4 \end{array}) and b_{4} = (\begin{array}{l} 1 \\ 1 \\ 1 \\ 1 \end{array}) .

Open in the MATLAB editor: f04mf_example

function f04mf_example


fprintf('f04mf example results\n\n');

t = [4; 3; 2; 1];
b = [1; 1; 1; 1];
x = [0];
work = zeros(9,1);
fprintf('  order  refl. coeff   solution\n');
for k=1:4
  [x, p, work, ifail] = f04mf( ...
                               t(1:k), b(1:k), x, work);
  if k > 1
    fprintf('%6d%11.4f    ', k, p);
  else
    fprintf('%6d%11.4s    ', k, ' ');
  end
  fprintf('%8.4f',transpose(x));
  fprintf('\n');
  if k < 4
    x = [x; 0]; % Extend x by one element
  end
end

f04mf example results

  order  refl. coeff   solution
     1                 0.2500
     2    -0.7500      0.1429  0.1429
     3     0.1429      0.1667  0.0000  0.1667
     4     0.1667      0.2000  0.0000 -0.0000  0.2000

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox