An error-free Levison algorithm to solve integer Toeplitz system

An Error-Free Levison Algorithm to Solve Integer Toeplitz System Miroslav MorhiE Institute of Physics Slovak Academy of Sciences 842 28 Bratislava, Czechoslovakia

Transmitted by John Casti

ABSTRACT This paper nonsymmetrical

presents

an error-free

when using classic algorithms, of original

form

Toeplitz system. It removes

Levison

deteriorate

algorithms.

of Levison

algorithm

to solve integer

roundoff errors. Their accumulation

can,

or destroy the solution. It retains the speed

Its modular

form

can

be easily implemented

for

parallel computer.

1.

INTRODUCTION

The standard method to solve the system of N linear equations with N unknowns, AZ = q, consists in finding of inverse matrix A-’ and its multiplication with vector g or in using a method like Gauss’ elimination or any other from direct or iterative algorithms. The methods of these types have computation complexity proportional to N 3. The internal special structure can be utilized in some cases to construct faster algorithm. One such a case is Toeplitz system with the matrix specified numbers ak, k = -N + 1, . . . , - IO, 1, . . . , N - 1 as

r

a0,a_l,a_2,...,a-N+2,a-N+l

a,, a,, a_l,.

APPLIED

MATHEMATICS

. . , a-N+3’

:!

by 2N -

1

Xl a-N+2

AND COMPUTATION

0 Elsevier Science Inc., 1994 655 Avenue of the Americas, New York, NY 10010

x2 = xN-1 XN

61:135-149

(1994)

135 0096-3003 /94 /$7.01)

M. MORHie

136 Computation

task here is to find vector

E from vector

g and matrix

practice, N may be large (e.g., 1,000 or more), and therefore ment of the fast algorithms is becoming more urgent.

A. In

the develop-

In past years, different kinds of fast algorithms appeared. The well-known Levison algorithm [I] can be employed upon the condition that the Toeplitz matrix is symmetrical. The Berlekamp-Massey algorithm [2] is applicable for a Toeplitz

system of the form AS = [-UN,

-UN+1 )...)

-Q+J.

(2)

We remind you also of the methods based on Euclid’s algorithm [l], and the methods applicable for upper or lower triangular Toeplitz systems, for Toeplitz-plus-Hankel systems [3], for Trench algorithm [l, The problem to solve Toeplitz system originates in including spectral estimation, linear prediction, error codes, and deconvolution. Mainly, the problems of deconvolution ill conditioned. The rounding-off errors, when ence very negatively the solution of the system necessity to look for an algorithm to eliminate Toeplitz systems. The aim of the paper is to derive an effective Toeplitz systems. The algorithm is designed for matrix A and vector Ij consists of integers), not a restrictive factor in any way. 2.

MATHEMATICAL

BACKGROUND

41, etc. many applications signal processing, are almost always

solving large systems, influ(1). This fact emphasizes the rounding-off errors in solving error-free algorithm to solve integer Toeplitz systems (the

but for many applications

OF THE

this is

ALGORITHM

The Levison’s method described in [l, Sj is well known. Its generalization to nonsymmetrical Toeplitz system with nonzero principal diagonal minors is presented in [6]. We shall make use of it to derive the error-free algorithm to solve integer nonsymmetrical Toeplitz system. Equation (1) can be written as

ai_jxj

f

= yi,

i = 1,2 )...,

iv.

(3)

j=l

Levison derived the algorithm dure solving the system

E

aj_jvjM)

j=l

for fast solution of (I) using recursive

= yi,

i = 1,2,...,

M

proce-

(4

Levinson Algorithm for Toeplitz Systems

for M = 1,2,.

. . , N . The vector

M = N, ijcN) equals j;. The regular integer Toeplitz

137

GCM) denotes

the result in Mth step. For

system (1) can be written in the form

(5) where D, is the determinant of the matrix A. The elements of the vector X and the determinant D,, are integers. Our task consists in finding the vector Z and the determinant D,. Then dividing th ese values one easily obtains the solution of (1). Analogously,

one can write integer form of the relation (4)

fai_jg =yi,

i = 1,2 ,..a> M.

(6)

j=l

Recursing

from step M to step M + 1 gives

where i = 1,2,.

. . , A4 + 1. By elimination

of yi in relations (6) and (7), one

has

and then

I=

'i-(M+l),

i = 1,2 ,...,

M.

(9)

M. MORHii:

138 Lettingi+M+l-iandj-+M+l-jyields

i = 1,2,...,

M

(10)

where

D M+l’M+l-jCM)

D,vLM,:?~

-

From (ll), v(M+r)

M+ l-j

j = 1,2,.

E b!M’ I ’

(M+l)

'M-t

1

.*>M.

(11)

it follows that = 1

D[M

CM) DM+l _ ‘M+l(M+l)b(M) j

‘M+l-j

j

17

=

1,2 ,...,

M.

(12)

To calculate components of the vector UCM“) in the M + 1st iteration step we need to know the quantities D, ,vb">] _j, bjM, from the Mth iteration step plus D,+ 1 and vL”+:r) from M + 1st iteration step. Substituting

i =

M +

1 into (7) gives f

uM+

l-j

(13)

YM+l

j=l

and hence

D M+l

(M+O=YM+lvM+l

_

_

a0

Letting

aM+l-jvj

a0

(M+l)

(14)

.

j=l

j --) M + 1 - j in (12) we get

v!M+l) .I

= _

1

D[M

After substituting

#)D~+~ I

j = 1,2 ,*.*>

_ vLM+;l)b(M) a4+1-j I )

(15) into (la),

the sought quantity

v&y:r’

M.

(15)

can be expressed

M YM+IDM

-

C"M+,-jvjM) j=l

(M+l)=

%f+

1

M

aoDM

-

C ‘M+,-jbLF-j j=1

D M+l

(16)

139

L.evinson Algorithm for Toeplitz Systems

It remains to express recursion to compute the vector 6 and determinant

D. Let us return to the beginning. Analogously to (6) the recursive procedure to solve Toeplitz system can be written

E

aj_,

g

=

i = 1,2 ,...,

yi,

M.

(17)

j=l The matrix of the set of linear equations (17) is the transposed matrix of system (6) and therefore the determinants D, of both systems in the Mth iteration step are equal. Again it holds that the vector GcN) in Nth iteration step represents the solution of (5). Then 1)(N) = fjj(W

=

(18)

x. >

i.e., the solution can be calculated in two ways. Employing the same procedure as for (6) one obtains the relation analogous to (10)

i = 1,2 ,...,

M,

(1%

where

&M) = I

CM) D M+lWM+l-j

-

D,w$~~$

j = 1,2 >**.> M.

(M+l)

(29)

wM+l

Comparing (19) to (6) one sees that the same resulting relations hold for the values c! M, as for the values v.(M). The only exception is that yi is replaced by ii. Then according to (l&) obviously, it holds that

uM+~DM -

EaM+l_jcj”)

j=l

cF++ll) =

M %DM

-

C aM+,-jbL"+'l-j j=l

D M+l

(21)

M. MORHLE

140 Making use of the same analogy in opposite

direction

yields

M a_,_,D, b(M+l) hi+1

-

CUj-M_&j(M) j=l

= ?

all&l-

j=l

(22) *

aj-M-lc(MM!l-j

D M+l To express cjMfl)

(see (15)) and bj”+l),

results in final relations .(M+l)

J

we apply the same procedure.

This

for recursions = -

1

DM [

CjM)DM+l

_ C(MMt+ll)bW)

$M)D~+~

_ b$M+;l)CW)

M+l-j

I

(23)

I *

(24)

and b!M+U = 1

I

DM [

M+l-j

The last problem that must be solved is the recursion to calculate determiLet us assume the value of the determinant in the step M is nant D,,,,. known. Let us divide the matrix ACM+I) to submatrices ..a

A

alTa2p.-.TaM

a0

ACM+])

=

a-1

a,,

a-2

a-1,

'-M

a l,*.*,“M-l

. . , aM-2

a,,.

.

.

.

.

+aT =

i

a0> -a,

ACM'

1 ’

a_M+1,a-M+2,...?a0

+a=[al,az,...,aM]’ -ii = In accordance with (lo),

[a_1,a_2,..-,a_M]T-

using matrix notation it holds

(24a)

Levinson Algorithm for Toeplitz Systems

141

Let us add the linear combination of columns of the matrix ACM+I)

+a

1

--

g(M)

D, [ ACM)1 to its first column. Then the matrix ACM+‘) has the form

A(M+l)

a0

=

1

--

+aT .&CM)

DM

Then the determinant of the matrix ACM+‘) is

--

D

1

+zT .5(M)

.D

M

D,

(25)

or

D Mtl

= DM*ao

-

f

j=l

ajbiM) = D,*a,

-

t

aM+,_jbLM>l_j.

(26)

j=l

Analogously to this derivation we can, starting from (19), come to the relation M

D M+l

=

DMUo

-

c

aj-M_lt$(!?-j.

(27)

j=l

Denominators of the expressions (23) and (24) equal 1, so we have

CM+11 ‘M+ 1

_ -

aM+lDM -

:

aM+l-jcj”)

(28)

j=l

and

b(M+U Mfl

=

a_M_lDM

-

j~laj-M-lb~M’.

(2%

M. MORHie

142

Relations (26), (161, (281, and (29) to determine A4 + 1, relations_(l5), (23), and (24) to determine the vectors G, Z, b, and initial conditions

form the basis of the Toeplitz system. 3.

(1) =

q(l)=yl,

Dl =a,,

ILLUSTRATIVE

effective

with the index components of

by) = a_,,

a,,

Cl

error-free

quantities remaining

algorithm

(30)

to solve the

integer

EXAMPLE

To make some facts, following from the above given relations, introduce the solution of integer Toephtz system for N = 3

clear, we

(31)

We start from initial values according

to (30) Cl(1) =

“Il) = y1,

Dl = ao,

al

b',')

,

=

a_

1.

We set M = 1, and using (291, (281, and (16), we calculate

b(,2)= a_,D,

-

~‘22) = a,D, v$‘) =

yz D,

-

a_,bl') a1c1 (I)

=

=

a-2ao

-

a2ao - a:,

a,v(,‘) = yzao - a, yl.

at,,

(32) (33) (34)

Using (26) or (27) we obtain

D, = D,a, - a,b(,‘) = at - ala_l

(35)

or

D, = D,a,

- a_,~(,‘)

= af, -

ala_,.

(36)

Levinson Algorithm for Toeplitz Systems

143

The resulting values are the same, no matter which of two relations we use. Now we continue

in calculation

by’ = +

of b$

[b’,“&

according

to (24)

- bb2’cr”] = a,~_,

- CZ__~U~,

(37)

1

cy) according

to (23)

cp’ =

+D,

- cf'b',"] =

aOal

u~u_~,

-

(38)

1

and u$

according

to (15)

VP =

&‘D, - ~$2’b~1’] = a,,~,

- y2u_,.

1

The

components

of vectors

previous assumptions. In the second iteration

5,

i?, and 5 are integers

that correspond

to

step for M = 2, again using (291, (281, and (161,

one obtains bk3) = a_,D,

- u_,bi2)

= u;u_3 ~6

-

u,a_,a_,

= a3 D, - u2ct2) 2 UOU,

=

- u_,bi2)

- u_lulu3

-

2u,u_,u_,

+ a,a”_,

+ a?,,

(40)

(2)

UlC2 -

2u,u,u,

+ u_iu;

+ a;,

(41)

where u3 = u_~ = 0, and (3) = 03

=

y3

D, - u2vi2) - u,ui2)

y&4

-

a0a2)

-

Y2(%%

- 82@-1) + Y3(4

- %a-,).

Comparing (40) to (40, one can see that the indices of both relations opposite signs. We continue in calculation of the determinant using (26) or (27)

(42) have

D, = a, D, - u2bi2) - u,bjz) =u

3 0 -

2u,u,u_,

- u2u0u_2

+ u,u2,

+ u$z_2

(43)

M. MORHk

144 or

D, = a,D,

- a_,ci2)

Now employing

- a2uou_2

2alaou_l

=a i -

u_,c(,~)

-

+ u,a2_,

+ ufu-2.

(44)

(241, (23), and OS), we calculate

b(3) = 1

@2'D,

_

b$+.$2' I

= af)a-l b(23) = ;

-

-

QlaOa-2

[ bi2)D,

-

a,&,

+ u2u-1a-2,

(45)

b’3’c;2’]

2 =+_,

c(3) = 1

=

L

-

u,u~_,

42)D,

-

+

ulu_la_z

u,aY,,

(46)

1

c$3Q,$2)

D[2

uou;- a-

-

,aoa2

(47)

u_1aT + a_2ulu2,

-

= aiu2 - aouf + a_lalu2

- a-24,

(48)

and

(3) 01

=

;

[

vi2)D3 - t~$~‘bi~‘]

2

=

(3) = 02

- a1a_I)

Y&4

- y2(aoa_l

- ala-z>

+

y3(u?,

-

w-2>7

(49)

[ oi2)D3 - t$)b(12)]

$ 2

=

-Yh%

-

u~u_~)

+

y2(ai

-

u,u-,)

-

y3(aoa-,

-a1u-2).

(50)

Levinson

The

Algorithm

resulting

for Toeplitz Systems

vector

145

x = l/D,[uy),

up), z$)lT

gives

the

solution

Toeplitz system (31). Its correctness is obvious. The presented accordance with assumptions given in Section 2 shows that l l l

4.

of the

example

in

all quantities are integers, determinants calculated according to (26) and (27) are equal, and indices of quantities “a” in corresponding elements have opposite signs.

ALGORITHM OF TOEPLITZ

OF ERROR-FREE SYSTEM

SOLUTION

From the above given example, one can see that all quantities are integers. After some iteration steps, however, their magnitudes increase rapidly, which complicates the realization of the calculation. To overcome this problem, we can proceed l

to perform

in several ways: the calculation

with these long integers

what means to watch

overflows; the calculation becomes more time consuming and inconvenient; . to perform the calculation in floating point form, which means to give up the accuracy of the algorithm; for ill-conditioned systems, this can give rise to the well-known problems with stability of the algorithm; and l

to perform the calculation in modulo arithmetic; one does not need watch overflows, and also the roundoff errors are removed.

When using modulo arithmetic, we carry out the calculation prime modulo classes m,, i = 1,2, . . . , r, so that it holds [7] as

to

in different

(51)

where

P must fulfill the condition

P > 2.max{NN/‘-M(

A)N,

N(N

1)‘N-1)‘2.M(

-

A)Np”M(Y)},

(52)

where

M(A)

= maxla,),

M(Y)

=

mmlfJjl,

iE(--N+l,N-l),

j

E

(1,

N).

(53)

(54)

The calculation in the given modulo class is independent of other modular classes, and thus this model of implementation is well suited for parallel computers with SIMD structure.

M. MORHx%

146 Now we introduce the complete basic algorithm Let us start with variables initialized as follows up =[a,,a,,a, an

= [U-l,

,..., u-2,

for one modulo class mk.

~,_~,qr,

u-3 ,...,

a_.+,,O]T,

3 = [Yl> YZ> Y37*-*r Yr.L constant A = a,, determinant

D =a,, and modulus

S=mk. We shall need working vectors

6 = [b,,b, ,.../ b‘?J, c=

[c,,c,

5 =

[u1,u2 ,..., UJ,

CJJT,

)...,

and auxiliary variables r, q, DP, and DI. Further let b, = u_~, c1 = a,, DP = D, and the iteration step M = 1. The procedure is then as = Yip

01

follows: (a) b[M

+ 1] =

D.un[M +

l](mod

S)

c[ M + 11 = D.up[ M + l](mod S) V[ M + 11 = D. y[M + l](mod S> DP = DP.A(mod S). (b) For j = 1,2,. . . , M calculate b[ M + l] = b[ M + I]- un[ M + 1 -j].b[j](mod S> c[ M + 11 = c[ M + 11- up[M + 1 -j].c[j](mod S) u[M + 11 = u[ M + I]- up[M + I - j].v[jl(mod S) DP = DP - up[M + 1 - j].b[ M + 1 - j](mod S). (c) Calculate inverse value of D (mod S> and assign DZ = D-‘(mod S),

d=

DP.

Levinson Algorithm for Toeplitz Systems

147

cd> For j = 1,2, . . . , M calculate v[j] u[j]

= u[j].D - u[M + l].b[M + 1 -jKmod S) = v[j].DZ(mod S). (e) If M = N - 1 finish the calculation, if not, continue

(f) For j = 1,2,. ?- = b[j] q=c[M+l-j]

in the step J

. . , M calculate

b[j]

= r. D - b[ M + ll.q(mod S) c[M + 1 -jl = q.D - c[M + ll.r(mod

= b[j]. DZ(mod S) c[M + 1 -jl = c[M + 1 -jl.DZ(mod ($ Increment iteration step

S>

b[j]

S).

M=M+l. (h) Repeat The

resulting

determinant

from step (a) on. vector

in the

given

modulo

class

is stored

in 6 and the

in D. Let us denote them Z;, and D,,

k = 1,2 ,...,

r.

By inverse conversion from residual representation, employing the Chinesse theorem [l, 81, the resulting vector X and determinant D is obtained according to (5). Division of vector jT: by D gives the solution of system (1).

EXAMPLE. Find exact solution of the Toephtz system

[~2~~q[r;]=

[ -iI.

We initialize the vectors and the variables F

= [3,2,01T,

an

= [-1,2,0]7

Ij = [-1,3,11r, 6 = [ -1,

-,

E = [3, -, 6 = [-l,-)-IT,

-I’, -IT,

148

M. MORHk

A = 1 and D = 1. With respect

to (51)

(52), (53)

and (54)

we choose the

moduli sr = 5, sp = 7, and sa = 11. We apply the above presented algorithm for given residual classes. Carrying out the first iteration step yields

z(mod5)

= [3, I, -IT;

c(mod 5) = [0,3,

-IT;

G(mod5)

-IT;

D(mod5)

= [2,1, = 4;

g(mod 7) = [O, 1, -IT;

&(mod 11) = 14, I, -I*

E(mod 7) = [5,0, Z(mod 7) = [2,6,

?(mod 11) = [5,4, G(mod 11) = [2,6, D(mod 11) = 4.

@mod

-IT; -IT;

7) = 4;

-IT -]r

After the second iteration step (M = N - 1 = 2) the calculation is finished. The resulting values of vectors 6, C, U and determinant D are then as follows:

Rmod 5) = [3,1, OIT; c(mod 5) = [O, 3,4]r;

Rmod 7) = [O, 1, llT; c(mod 7) = [5,0, 41T;

Z(mod 5) = [l, 3, 21T; D(mod5) = 3;

G(mod 7) = [2,3, D(mod 7) = 2;

By the conversion

of the vector

31T;

Z(mod 11) = [4,1,4]’ Xmod 11) = [5,4,O]r Zj(mod 11) = [5,3, 41T D(mod 11) = 1.

U and the determinant

D from residue class

code, one obtains U = [l6,3,367]r,

D = 23.

Let us suppose that positive numbers are represented in the range (0, (I’ 1)/2) and negative numbers in the range ((P - 1)/2 + 1, P - 1). Then the negative values are, for the values being greater than (P - 1)/2, calculated by subtracting P. In our case, the third element of resulting vector 367 - 385 = - 18 (367 > 5.7.11/2).Then the final solution is

if = ;[16,3,

u(3) =

-181T.

In the whole algorithm, there does not occur any operation of division so that the rounding-off errors are excluded of the solution. The only operation worth noting is the finding of inverse number D-‘(mod s). As the modulus is supposed to be prime, we can apply the Euclid’s algorithm described, e.g., in [l]. The vector of inverse numbers can be calculated beforehand in form of table. Then, in the run-time, we need just to assign for a value D E (1, s - 1) the appropriate D-‘(mod s).

149

Levinson Algorithm for Toeplitz Systems 5.

COMPLEXITY

OF THE

ALGORITHM

AND CONCLUSION

The complexity of the algorithm can be expressed in the number of necessary multiplications and additions. For the above presented algorithm, the number of needed modulo multiplications is

N, = T.N.(N

and number

of needed

-

I) + 4.c~

-

I)

module additions is

N, = ~.N.(N

When using a classic computer,

-

I).

the entire calculation

must be repeated

for

all residual classes, which increases the calculation time. The individual processes are independent of each other, and thus they can be carried out in parallel. Using parallel computer makes it possible, simultaneously with precision,

to achieve

also a high-speed

calculation

of the integer

Toeplitz

system. REFERENCES R. E. Blahut,

Fast Algorithms for Digital Signal Processing,

IBM

Corporation,

Owego, NY, 1985. E. R. Berkelamp,

Algebraic Coding Theory, McGraw-Hill,

G. Heinig, P. Janowski, Hankel matrices, W. F. Trench, 12:512-522 E.

Numer.

Math. 52:665-682

An algorithm for the inversion of finite Toeplitz matrices,

1. SIAM

and

S. Treibel,

Geophysical

Signal Analysis,

Prentice-Hall,

Cliffs, NJ, 1980.

H. W. Press et al., Numerical Recipes, Cambridge York/New

of Toeplitz-plus-

(1988).

(1964).

A. Robinson

Englewood

New York, 1968.

and K. Rost, Fast inversion algorithms

Rochelle/Melbourne/Sydney,

M. Newman,

Solving equations

Univ. Press, Cambridge/New

1986.

exactly,

Nut.

Bureau

Standards

71B:171-179

(1967). H.

J.

Nussbaumer,

Springer-Verlag,

Fast

Berlin,

Fourier

1981.

Transform

and

Convolution

Algorithms,