Научная статья на тему 'Универсальные алгоритмы для обобщенных дискретных матричных уравнений Беллмана с симметрической матрицей теплица'

Универсальные алгоритмы для обобщенных дискретных матричных уравнений Беллмана с симметрической матрицей теплица Текст научной статьи по специальности «Математика»

CC BY
98
45
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
UNIVERSAL ALGORITHMS / MATRIX BELLMAN EQUATIONS / TOEPLITZ MATRIX / SEMIRINGS / IDEMPOTENT SEMIRINGS

Аннотация научной статьи по математике, автор научной работы — Сергеев С. Н.

This paper presents two universal algorithms for generalized discrete matrix Bellman equations with symmetric Toeplitz matrix. The algorithms are semiring extensions of two well-known methods solving Toeplitz systems in the ordinary linear algebra.

i Надоели баннеры? Вы всегда можете отключить рекламу.
iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

UNIVERSAL ALGORITHMS FOR GENERALIZED DISCRETE MATRIX BELLMAN EQUATIONS WITH SYMMETRIC TOEPLITZ MATRIX

This paper presents two universal algorithms for generalized discrete matrix Bellman equations with symmetric Toeplitz matrix. The algorithms are semiring extensions of two well-known methods solving Toeplitz systems in the ordinary linear algebra.

Текст научной работы на тему «Универсальные алгоритмы для обобщенных дискретных матричных уравнений Беллмана с симметрической матрицей теплица»

MSC 68W30, 15A80, 15B05, 16Y60

Universal algorithms for generalized discrete matrix

Bellman equations with symmetric Toeplitz matrix

l

© S. N. Sergeev

Ecole Polvtechnique, Paris, France

This paper presents two universal algorithms for generalized discrete matrix Bellman equations with symmetric Toeplitz matrix. The algorithms are semiring extensions of two well-known methods solving Toeplitz systems in the ordinary linear algebra

Keywords: universal algorithms, matrix Bellman equations, Toeplitz matrix, semirings, idempotent semirings

1. Introduction

As observed by B, A, Carré and E, C, Backhouse [2], [4], the Gaussian elimination without pivoting can be viewed as a prototype for some algorithms on graphs, M, Gondran [8] and G, Rote [21] made this observation precise by proving that the Gausssian elimination, under certain conditions, can be applied to the linear systems of equations over semirings. The notion of universal algorithm over semiring was introduced by G, L, Litvinov, V, P, Maslov and E, V, Maslova in [10], [15], These papers are to be considered in the framework of publications [9], [11], [12], [17], [18], [20] of the Russian idempotent school, and more generally in the framework of idempotent and tropical mathematics, see [1], [5], [14], [16] and references therein. Essentially, an algorithm is called universal if it does not depend on the computer representation of data and on a specific realization of algebraic operations involved in the algorithm [15], Linear algebraic universal algorithms include generalized bordering method, LU- and LDM-decompositions for solving matrix equations. These methods are basically due to B, A, Carré, see also [15],

It was observed in [13], [15] that universal algorithms can be implemented by means of objective-oriented programming supported by C++, MATLAB, Scilab, Maple and other computer systems and languages. Such universal programs can be instrumental in many areas including the problems of linear algebra, optimization theory, and interval analysis over positive semirings, see [13], [17], [18], [19],

1This work is partially supported by the RFBR-CRNF grant 11-01-93106

1751

This paper presents new universal algorithms based on the methods of Durbin and Levinson, see [7], Sect, 4,7, These algorithms solve systems of linear equations with symmetric Toeplitz matrices. Our universal algorithms have the same complexity O(n2) as their prototypes which beats the complexity O(n3) of the LDM-decomposition method. All algorithms are described as MATLAB-programs (following the style of [7]), meaning that they can be actually implemented.

The author is grateful to G, L, Litvinov and A, N, Sobolevskii for drawing his attention to this problem and for valuable discussions.

§ 2. Semirings and universal algorithms

A set S equipped with addition © and 'multiplication © is a .semiring (with unity) if the following axioms hold:

1) (S, ©) is a commutative semigroup with neutral element 0;

2) (S, ©) is a semigroup with neutral element 1 = 0;

3) a © (b © c) = (a © b) © (a © c), (a © b) © c = (a © c) © (b © c) for all a,b,c G S (distributivitv);

4) 0 © a = a © 0 = 0 for all a G S,

©

The semiring S is called idempotent if a © a = a for any a G S. In this case © induces the canonical partial order relation

a -< b ^ a © b = b. (1)

The semiring S is called complete (cf. [6]), if any subset {xa} C S is summable and the infinite distributivitv

c © (© xa) = ©(c © Xa),

a a

(© Xa) © c = 0(Xa © c),

aa

holds for all c G S and {xa} C S, This property is natural in idempotent semirings and also in the theory of partially ordered spaces (cf, G, Birkhoff [3]) with partial

a

Consider the closure operation

<X>

a* = © a%. (2)

i=0

1752

In the complete semirings it is defined for all elements. The property

a* = 1 © aa* = 1 © a*a,

(3)

reveals that the closure operation is a natural extension of (1 — a)-1.

We give some examples of semirings living on the set of reals R totally ordered by <: the semiring R+ with customary operations © = +, © = ■ and neutral elements

0 = ^d 1 = 1; the semiring Rmax = R U { — ro} with operations © = max, © = + and neutral elem ents 0 = — ro, 1 = 0; the semir ing R max = Rmax U {ro}, which is a completion of Rmax with the element ro satisfying a © ro = ro for all a a ©ro = ro© a = ro for a = 0 and 0 ©ro = ro© 0 = 0; the semiring Rmax,min = R U {ro} U { — ro} with © = max © = min, 0 = —ro, and 1 = ro.

Consider operation (2) for the examples above. In R+ the closure a* equals (1 — a)-1 if a < 1 and is undefined otherwise; in Rmax it equaIs 1 if a < 1 and is undefined otherwise; in Rmax we have a* = 1 for a < ^d a* = ro for a > 1; in Rmax,min we have a* = 1 for all a. Note that Rmax and Rmax,min are a-eomplete, so the closure is defined for any element of these semirings,

©©

linear algebra. Denote by Matmn(S) the set of all m x n matrices over the semiring

S. By In we denote the n x n unity matrix, that is, matrix with 1 on the diagonal and 0 off the diagonal. As usual, we have AIn = InA = A and A0 = In for any A G Matnn(S). The set Matnn(S) of all n x n matrices is a semiring. Its unity is In and its zero is 0n, the square matrix with all entries equal to 0. If S (and hence Matnn(S)) is complete, the closure A* is defined for any matrix A and it satisfies (3), Note that if S is partially ordered, then Matnn(S) is ordered elementwise: A ^ B iff Aij ^ Bj for all i = 1,... ,m and j = 1,..., n. If S is idempotent and canonically ordered (1), then the elementwise order of Matnn(S) also satisfies (1),

The closure operation of matrices is important for the (discrete stationary) matrix Bellman equations

a-complete idempotent semirings the matrix A*B is the least solution of equation (4) with respect to the canonical order (1),

Since A* is a generalization of (I — A)-1, the known universal algorithms for A* are generalizations of the methods for matrix inverses, and the known algorithms for Bellman equations are generalizations of the methods for AX = B, Further we consider the generalized bordering method (see also [2,4,8,15,21]),

Let A be a square matrix. Closures of its main submatrices Ak can be found inductively. The base of induction is A1, the closure of the first diagonal entry. Generally we represent Afc+1 as

assuming that we have found the closure of Ak. In this represent at ion, and hk are

X = AX © B.

(4)

If the closure of A exists and (3) holds, then X = A*B is a solution of (4), In

1753

columns with k entries and ak+1 is a scalar. We also represent Ak+1 as

W Mfc+1

Uk vk

Using (3) we obtain that

uk+1 = (hfc Akgk © ak+1)* , vk Ak gk uk+1, .

wTk = uk+1hTAk, ( j

Uk = Ak gk uk+1hT Ak © Ak.

Consider the bordering method for finding the solution x = A*b to (4), where X = x and B = b are column vectors. Firstly, we have x(1) = A^, Let x(k) be the (k — 1)

z = x(k) © Ak gk xk+1.

We have to compute Akgk. In the next section we show that this can be done

A

We also note that the bordering method described by (5) and (6) is valid more generally over Conway semirings, see [6] for definition.

§ 3. Universal algorithms for Toeplitz linear systems

Formally, a matrix A G Matnn(S) is called Toeplitz if there exist scalars ri;

i = —n +1,..., 0,..., n — 1, such th at Aij- = rj-i for a 11 ¿and j. Informally, Toeplitz matrices are such that their entries are constant along any line parallel to the main diagonal (and along the main diagonal itself). For example,

is Toeplitz, Such matrices are not necessarily symmetric. However, they are always persymmetric, that is, symmetric with respect to the inverse diagonal. This property

Using (5) we obtain that

Xk+1 = Uk+1(hTx(k) © bk+1),

(6)

/ ro r1 r2 rs\

r-1 ro r1 r 2

r-2 r-1 ro r1

\r_3 r-2 r-1 ro/

1754

is algebraically expressed as A = EnATEn, where En = [en,..., ei]. By e* we denote the column whose ith entry is 1 and other entries are 0. The property E^ = (where is th e nxn identity matrix) implies that the product of two persymmetric matrices is persymmetric. Hence any degree of a persymmetric matrix is persymmetric, and so is the closure of a persymmetric matrix. Thus, if A is persymmetric, then

E„A* = (A*)T En. (7)

Further we deal only with symmetric Toeplitz matrices. Consider the equation y = Tny © r(n), where r(n) = (r1,..., rn)T, and Tn is defined by the scalars r0, r1,..., rn-1 so that Tj = r|j-i| for a 11 i and j, This is a generalization of the Yule-Walker problem [7]. Assume that we have obtained the least solution y(;) to the svstem y = T;y © r(;) for some k such th at 1 < k < n — 1, where T; is the ma in k x k submatrix of Tn, We write Tk+1 as

Tk

E;r(;)\ k+1 = U(;)T Efc ro )

We also write y(k+1^d r(k+1) as

„(fc)

y(k+1)=i ¿j ■ r(k+1)=

Using (6), (7) and the identity T;r(;) = y(;), we obtain that

a; = (ro © r(;)Ty(k))* (r(k)TE;y(;) © r;+1), z = E;y(;)a; © y(;).

Denote 3; = r0 © r(;)Ty(;), The following argument shows that 3; can be found recursively if (3;_ 1)-1 exists.

3; = ro © [r(;-1)T r;] p;-1y(;-1)a;-1 © y(;-1)^

V a;-1 J

k-1

= ro © r(;-1)Ty(;-1) © (r(;-1)TE;-1y(;-1) © r;)«;- = ^

= £;_1 © (^^-1)-1 © «2-1.

Existence of (^k-1 )-1 is not universal, and this will make us write two versions of our

algorithm, the first one involving (8), and the second one not involving it. We will

write these two versions in one program and mark the expressions which refer only

to the first version or to the second one by the MATLAB-style comments %1 and

%2, respectively. Collecting the expressions for 3;, a; and z we obtain the following

y(k)

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

3; = ro © r(;)Ty(;), %2

3; = 3;-1 © (^fc-1)-1 © a2-1, %1

a; = (/3;)* © (r(;)TE;y(;) © r;+1), (9)

y(;+1) = (E;y(;)a; © y(;r

V a;

1755

Recursive expression (9) is a generalized version of the Durbin method for the Yule-Walker problem [7]. Using this expression we obtain the following algorithm.

Algorithm 1. The Yule-Walker problem for the Bellman equations with symmetric Toeplitz 'matrix.

function y = durbin(ro, r) n = size(r) + 1

y(1) = ro © r(1)

3 = ro %1

a = ro © r(1) for k = 1 : n — 1

3 = ro © r(1 : k) © y(1 : k) %2

3 = 3 © (3*)-1 © a2 %1

a = 3* © (r(k : —1:1) © y(1 : k) © r(k + 1)) z(1 : k) = y(1 : k) © y(k : —1:1) © a

y(1 : k) = z(1 : k)

y(k + 1) = a end

Now we consider the problem of finding x(n) = T*6(n) where Tn is as above and

6(n) = (61,..., bn) is arbitrary. We also introduce the column vectors y(;) which solve the Yule-Walker problem: y(;) = T;r(;), The main idea is to find the expression for x(;+1) = Tfc*+16(;+1) involving x(;) and y(;). We write i(;+1^d b(;+1) as

Making use of the persymmetry of Tfc* and of the identities T£bk = x(k) and Tfc*rk =

y( k), we specialize expressions (6) and obtain that

^ = (ro © r(k)Ty(k))* © (r(k)TEfcx(fc) © bfc+i), v = fc © x(k).

The coefficient r0 © r( k)Ty( k) = Ak can be expressed again as Ak = Ak-1 © (^^-1)-1 ©

(ak-1)2, if the closure (Ak-1)* is invertible. Using this we obtain the following recursive expression:

Ak = ro © r(k)Ty(k), %2

Ak = Ac-1 © (Afc-1)-1 © «2-1, %1

^ k = ft © (r(k)TEfcx( k)© bk+1), (10)

x( k+1) = (Eky(k)^ k © x(k)^

V ^k /

Expressions (9) and (10) yield the following generalized version of the Levinson algorithm for solving linear symmetric Toeplitz systems [7]:

Algorithm 2. Bellman system with symmetric Toeplitz matrix

1756

function y = levinson(ro, r, 6) n = size(b)

y(1) = ro © r(1); x(1) = ro © 6(1);

3 = ro %1

a = ro © r(1) for k = 1 : n — 1

3 = ro © r(1 : k) © y(1 : k) %2

3 = 3 © (3*)-1 © a2 %1

^ = 3* © (r(k : —1:1) © x(1 : k) © 6(k + 1))

v(1 : k) = x(1 : k) © y(k : —1:1) © ^

x(1 : k) = v(1 : k)

x(k + 1) = ^

if k < n — 1

a = 3* © (r(k : —1:1) © y(1 : k) © r(k + 1)) z(1 : k) = y(1 : k) © y(k : —1:1) © a

y(1 : k) = z(1 : k)

y(k + 1) = a

end

end

The computational complexity of all methods described in this section is O(n2),

References

1, F, L, Baeeelli, G, Cohen, G, J, Olsder, and J, P. Quadrat, Synchronization and Linearity: an Algebra for Discrete Event Systems, Wiley, 1992,

2, E, C, Backhouse and B, A, Carré, Regular algebra applied to path-finding

problems, J, Inst, Math, and Appl,, 1975, vol. 15, 161-186,

3, G, Birkhoff, Lattice theory. Amer, Math, Soe,, 1967,

4, B.A, Carré, An algebra for network routing problems, J, Inst, Math, and Appl,,

1971, vol. 7.

5, E. A. Cuninghame-Green, Minimax Algebra, Lecture Notes in Economics and Mathematical Systems, vol. 166, Springer, Berlin, 1979,

6, J, Golan, Semirings and their applications, Kluwer, 2000,

7, G, H, Golub and C, van Loan, Matrix Computations, John Hopkins Univ. Press, Baltimore and London, 1989,

8, M, Gondran, Path algebra and algorithms, In: B, Eov (ed), Combinatorial programming: methods and applications, Eeidel, Dordrecht, 1975, 137-148,

9, V, N, Kolokoltsov and V, P. Maslov, Idempotent analysis and its applications, Kluwer Academic Pub,, 1997,

1757

10. G, L, Litvinov and V, P. Maslov, The correspondence principle for idempotent calculus and some computer applications, In: J, Gunawardena (ed), Idempotenev, Cambridge Univ. Press, Cambridge, 1998, 420-443,

11. G, L, Litvinov, V, P, Maslov, and G, B, Shpiz, Idempotent functional analysis: An algebraic approach, Mat, Zametki, 2001, vol. 69, No, 5, 758-797,

12. G, L, Litvinov and G, B, Shpiz, Dequantization of mathematics and group representations, Vestnik Tambov Univ. Ser: Est, tekhn, nauki, 2005, vol. 10, issue

4, 390-411.

13. G. L. Litvinov and V. P. Maslov. The correspondence principle for idempotent calculus and some computer applications. In: J. Gunawardena (ed), Idempotenev, Cambridge Univ. Press, 1994, 420-443. E-print arXiv: math/0101021.

14. G. L. Litvinov and V. P. Maslov (Eds). Idempotent mathematics and mathematical physics. Contemporary Mathematics, AMS, 2005, vol. 377.

15. G. L. Litvinov and E. V. Maslova. Universal numerical algorithms and their software implementation. Programming and Computer Software, 2000, vol. 26, No.

5, 275-280. E-print arXiv: math.SC/0102114.

16. G. L. Litvinov and S. N. Sergeev (Eds). Tropical and Idempotent Mathenatics, Contemporary Mathematics, AMS, 2009, vol. 495.

17. G. L. Litvinov and A. N. Sobolevskii, Exact interval solutions of the discrete Bellman equation and polynomial complexity in interval idempotent linear algebra, Dokladv Mathematics, 2000, vol. 62, No. 2,199-201. E-print arXiv math.LA/010141.

18. G. L. Litvinov and A. N. Sobolevskii. Idempotent interval analysis and optimization problems, Reliable Computing, 2001, vol. 7, No.5, 353-377.

19. P. Loreti and M. Pedicini. An object oriented approach to idempotent analysis: Integral equations as optimal control problems. In: G. L. Litvinov and V. P. Maslov (Eds), Idempotent Mathematics and Mathematical Physics, Contemporary Mathematics, AMS, 2005, vol. 377, 187-208.

20. V. P. Maslov. Operator methods, M,: Mir, 1987.

21. G, Rote, A systolic array algorithm for the algebraic path problem, Computing, 1985, vol. 34, 191-219,

1758

i Надоели баннеры? Вы всегда можете отключить рекламу.