A Parallel Version of a Quasigradient Method in Stochastic Control Theory

Thomas Pohl, Wilfried Grecksch

Martin-Luther-Universität Halle-Wittenberg, FB Mathematik und Informatik, Institut für Optimierung und Stochastik, Theodor-Lieser-Str. 5, D-06099 Halle, Germany pohl@mathematik.uni-halle.de, grecksch@mathematik.uni-halle.de

Holger Blaar

Martin-Luther-Universität Halle-Wittenberg, FB Mathematik und Informatik, Institut für Informatik, Kurt-Mothes-Str. 1, D-06099 Halle, Germany blaar@informatik.uni-halle.de

Abstract: We solve an optimal control problem for controlled parabolic Ito equations by a stochastic quasigradient method. Because of high amounts of computation time required by numerical solution of such problems we investigate the parallelization of the algorithm. We distribute the computations of space stages over several processor nodes of a parallel computer. We obtain an efficient algorithm with low communication cost by using a ring topology.

Key Words: Stochastic Optimization, Control, Parallel Algorithms

AMS Classification: 65Y05, 93E20, 93E23

1 Introduction

Simulation by mathematical models often requires high amount of computation time. Parallel computers are well suited to perform such computations in less time by distributing necessary computations to several processor nodes. So several computations or several algorithm parts can be carried out at the same time. Communication and synchronization costs of the parallel algorithm can produce overhead and can reduce the advantage of parallel computation. In general, it is not a trivial problem to transform a usual sequential algorithm into an efficient parallel algorithm or to create a new, powerful parallel algorithm. In addition, the efficiency of a parallel implementation depends on the type of parallel computer used. One of the first steps of parallel algorithm design should be a theoretical performance evaluation by determining the complexity of both the fastest sequential and a parallel algorithm as well as solving the same problem. Especially optimal deterministic and stochastic control problems are expensive tasks of this nature. The evolution of optimal conditions of dynamic programing or maximum principles leads to the problem of solving equations of Riccati type, if the problems considered are linear regulators or are close to the linear regulator case [FR75]. Consequently, there are papers which construct parallel algorithms to solve Ricatti systems (see for example the deterministic case and the stochastic case [ZS93]). Here we choose another point of view. A stochastic control problem can be interpreted by an optimization problem in a Hilbert space. In the finite dimensional case, optimization problems are solved by random search methods, for example in [G79], and by quasigradient methods [E76]. Methods of stochastic approximation are often examples of stochastic quasigradient methods. Stochastic approximation problems are discussed in the infinite dimensional case in [Y92]. Here we consider an optimal control problem for a controlled parabolic Ito equation with respect to white noise in time and space and determine ε-optimal control by time and space discretization of the state equation and the objective functional. These problems are optimization problems of high dimensions and the optimization problems are solved by a stochastic quasigradient method. The parallelization idea is to distribute the computation of space stages to several processors of a message passing system, that is a parallel computer without global memory. The communication between processors is solved by using a ring topology. This simple topology gives a sufficient efficient parallel solution, and ring topology algorithms can be easily mapped to current parallel computer architectures.

2 A Stochastic Quasigradient Method for an Optimization Problem on Hilbert Space

We define a convex functional F on a closed and convex set G of a separable Hilbert space H with the scalar product ⟨·, ·⟩. Let us consider the minimization problem

min{F (u) : u ∈ G ⊂ H}.  (1)

We assume, that a solution of (1) exists. We don’t know the gradient of F, but we have the possibility to get observations. All random variables are defined on a complete probability space (Ω, F, P). E(·) denotes the mathematical expectation with respect to P. The stochastic quasigradient method defines a sequence of H-valued random variables X^{s} by

X^{s+1} := Pr_{G} (X^{s} − ρ_{s}γ_{s}Y^{s}) , s = 0, 1, 2, ….  (2)

Y^{s} describes a random variable with values in H satisfying

E {Y^{s}|X^{0}, …, X^{s}} = a_{s}F̂(X^{s}) + b_{s},  (3)

where b_{s} is a H-valued random variable and a_{s} is a positive real-valued random variable. The random variable b_{s} is measurable with respect to σ-algebra B_{s} := σ(X^{0}, …, X^{s}) and X^{0} is an element of G. Let us write F̂(X^{s}) for the subgradient of F in X^{s}, that is: F̂(X^{s}) is a B_{s}-measurable H-valued variable with

F̂(X^{s})(X − X^{s}) ≤ F (X) − F (X^{s})

for all X ∈ H. If the gradient exists, then F̂(X^{s}) is the gradient (see [BP86] page 107). We denote the step size by ρ_{s} , γ_{s} is a normalizing factor and P r_{G} is the projection operator from H onto the set G. Here it is shown that in the infinite dimensional case, X^{s} defines a minimizing sequence with probability 1. We can transform the proof of the finite dimensional case to the above infinite dimensional case. The concept of the Féjer sequence is useful for the proof [E76].

Definition 2.1 ([E76]) A sequence of H-valued random variables (Z^{s}) ⊂ H with ||Z^{0}||^{2} < ∞ called generalized Féjer-Sequence with respect to G, if we obtain for any y ∈ G:

E {||y − Z^{s+1}||^{2} | Z^{0}, …, Z^{s}} ≤ ||y − Z^{s}|| + W_{s}, s = 0, 1, …,

where W_{s} are real non-negative random variables, measurable with respect to the σ-algebra B_{s} := σ(Z^{0}, …, Z^{s}) and ∑_{s=0}^{∞} EW_{s} < ∞.

Theorem 2.1 If (Z^{s}) is a generalized Féjer-Sequence in G, then we obtain:

1. lim_{s→∞} ||y − Z^{s}||^{2} P-a.s., ∀ y ∈ G.

2. The set of limit points of Z^{s}(ω) is not empty for almost all ω ∈ Ω.

3. If Z′(ω), Z′′(ω) are two limit points of (Z^{s}(ω)), which are not in G, then the set G is in a hyperplane, whose points have the same distance from Z′ and Z′′.

Proof If we want to prove the second statement, it is sufficient to prove the first statement. If

Z̃_{s} = ||y − Z^{s}||^{2} + ∑_{k=s}^{∞} W_{k}

then Z̃_{s}

E(Z̃_{s+1}|Z̃_{0}, …, Z̃_{s}) ≤ Z̃_{s}, s = 0, 1, ….

Therefore {Z̃_{s}} is a supermartingale. The first statement follows from a known convergence theorem on supermartingales (see for example [E76], page 20). It follows from 1. that all limit points of (Z^{s}(ω)) lie on a sphere with this radius. If Z′(ω), Z′′(ω) are two limit points then we obtain for any fixed y ∈ G

||y − Z′(ω)||^{2} = ||y − Z′′(ω)||^{2}.

Thus

2(y, Z′′(ω) − Z′(ω)) + ||Z′(ω)||^{2} − ||Z′′(ω)||^{2} = 0. •

Proposition 2.2 If Z(ω) is a limit point of Z^{s}(ω) in G for any ω, then Z(ω) is limit of Z^{s}(ω).

Proof Let Z′(ω) be an limit point in G for any ω. We assume that there is a second limit point Z′′(ω) with Z′′(ω) ≠ Z′(ω). We obtain from (1) of Theorem 2.1

||y − Z′(ω)||^{2} = ||y − Z′′(ω)||^{2}.

If we choose y = Z′(ω), then

||Z′(ω) − Z′′(ω)||^{2} = 0,

in contradiction of our assumption. Thus Z′(ω) is the limit of the sequence. •

Theorem 2.3 Let η_{s} ≥ δ > 0 be random variables which are measurable with respect to B_{s}. For all L > 0 we have an C_{L} > 0 such that

E {||Y^{s}||^{2} |X^{0}, …, X^{s}} ≤ η_{s}^{2} ≤ C_{L} for ||X^{k}|| ≤ L, k = 0, 1, …, s.  (4)

Further we have numbers \underline{γ}, \overline{γ} > 0 such that for every s the equation

0 < \underline{γ} ≤ γ_{s} (η_{s} + τ_{s}||X^{s}||) ≤ \overline{γ} < ∞  (5)

is true where

τ_{s} = \{ \begin{matrix}1 : & ||b_{s}|| > 0 \\ 0 : & ||b_{s}|| = 0. \end{matrix}

Let ρ_{s} ≥ 0 and a_{s} ≥ 0 and

∑_{s=0}^{∞} E (ρ_{s}||b_{s}|| + ρ_{s}^{2}) < ∞.  (6)

Then we have:

1. The sequence (X^{s}), defined by (2), is a generalized Féjer-Sequence with respect to G^{∗}, where G^{∗} denotes the set of solutions of the optimization problem.

2. If ∑_{s=0}^{∞} ρ_{s}a_{s} = ∞ (f.s.), then (X^{s}(ω)) converges to (X^{∗}(ω)) ∈ G^{∗} for almost every ω ∈ Ω.

The theorem in the case of finite dimension is proven in [E76]. We can use the main idea of this proof for the case of infinite dimension. Proof of Theorem 2.3 From (5) we obtain

γ̃_{1} ≤ γ_{s} ≤ γ̃_{2}.

Let us assume X^{∗} be from G^{∗}. The projection property gives

||X^{∗} − X^{s+1}||^{2} ≤ ||X^{∗} − X^{s} + ρ_{s}γ_{s}Y^{s}||^{2}
= (X^{∗} − X^{s}, X^{∗} − X^{s} + ρ_{s}γ_{s}Y^{s}) + (ρ_{s}γ_{s}Y^{s}, X^{∗} − X^{s} + ρ_{s}γ_{s}Y^{s})
= ||X^{∗} − X^{s}||^{2} + ρ_{s}^{2}γ_{s}^{2}||Y^{s}||^{2} + 2ρ_{s}γ_{s} (Y^{s}, X^{∗} − X^{s}) .

Further we obtain from (3)

E {||X^{∗} − X^{s+1}||^{2}|B_{s}}
≤ ||X^{∗} − X^{s}||^{2} + 2ρ_{s}γ_{s}a_{s} (F̂(X^{s}), X^{∗} − X^{s})
+2ρ_{s}γ_{s} (b_{s}, X^{∗} − X^{s}) + ρ_{s}^{2}γ_{s}^{2}E (||Y^{s}||^{2}|B_{s}) .  (7)

Using the definition of subgradients and the optimality of X^{∗} we have (F̂(X^{∗}), X^{∗} − X^{s}) ≤ F (X^{∗}) − F (X^{s}) ≤ 0.

The Schwarz inequality and (5) give

2ρ_{s}γ_{s} (b_{s}, X^{∗} − X^{s}) ≤ 2ρ_{s}γ_{s}||b_{s}|| ||X^{∗} − X^{s}|| ≤ 2ρ_{s}||b_{s}|| (γ̃_{2}||X^{∗}|| + \overline{γ})

and consequently,

E {||X^{∗} − X^{s+1}||^{2} |B_{s}}
≤ ||X^{∗} − X^{s}||^{2} + 2ρ_{s}||b_{s}|| (γ̃_{2}||X^{∗}|| + \overline{γ}) + \overline{γ}^{2}ρ_{s}^{2}.  (8)

It follows then from the assumptions

∑_{s=0}^{∞} E (ρ_{s}||b_{s}|| (γ̃_{2}||X^{∗}|| + \overline{γ}) + \overline{γ}^{2}ρ_{s}^{2}) < ∞.

Therefore (X^{s}) is a generalized Féjer-Sequence with respect to G^{∗}. To prove the second part of the theorem, we show that one of the limit points is in G^{∗} for almost all ω. Then Proposition 2.2 gives the second statement of the theorem. Using the estimate of the first part of the proof recursively we then have

||X^{∗} − X^{s+1}||^{2} ≤ ||X^{∗} − X^{0}||^{2} + 2 ∑_{k=0}^{s} ρ_{k}γ_{k}(Y^{k}, X^{∗} − X^{k})
+ ∑_{k=0}^{s} ρ_{k}^{2}γ_{k}^{2}||Y^{k}||^{2}.

By the properties of the conditional expectation and (7), (8) we obtain the following inequality:

0 ≤ E||X^{∗} − X^{s+1}||^{2} ≤ E||X^{∗} − X^{0}|| + 2E ∑_{k=0}^{s} ρ_{k}a_{k}γ_{k}(F̂(X^{k}), X^{∗} − X^{k})
+ 2(γ̃_{2}||X^{∗}|| + \overline{γ}) ∑_{k=0}^{s} Eρ_{k}||b_{k}|| + \overline{γ}^{2} ∑_{k=0}^{s} Eρ_{k}^{2}.

Therefore we have

E||X^{∗} − X^{0}|| + 2E ∑_{k=0}^{∞} ρ_{k}a_{k}γ_{k}(F̂(X^{k}), X^{∗} − X^{k})

+2(γ̃_{2}||X^{∗}|| + \overline{γ}) ∑_{k=0}^{∞} Eρ_{k}||b_{k}|| + \overline{γ}^{2} ∑_{k=0}^{∞} Eρ_{k}^{2} ≥ 0

and by (6) we obtain

E ∑_{k=0}^{∞} ρ_{k}a_{k}γ_{k}(F̂(X^{k}), X^{∗} − X^{k}) > −∞.

From ∑_{k=0}^{∞} ρ_{k}a_{k} = ∞ a. s. we know that a subsequence (s_{k}) = (s_{k}(ω)) exists for almost all ω with the property

γ_{s_{k}}(F̂(X^{s_{k}}), X^{∗} − X^{s_{k}}) ⟶^{k→∞} 0. P − a.s.

It follows from the assumptions (4) and (5) that γ_{s_{k}} is bounded. Hence we have

(F̂(X^{s_{k}}), X^{∗} − X^{s_{k}}) ⟶^{k→∞} 0.

By

F (X^{∗}) − F (X^{s_{k}}) ≥ (F̂(X^{s_{k}}), X^{∗} − X^{s_{k}}) ⟶^{k→∞} 0

we obtain

0 ≤ F (X^{s_{k}}) − F (X^{∗}) ≤ (F̂(X^{s_{k}}), X^{s_{k}} − X^{∗}) ⟶^{k→∞} 0.

Thus (X^{s_{k}}) is a minimizing sequence. Therefore each limit point X′ of the sequence (X^{s_{k}}) satifies the relation F (X′) = F (X^{∗}). From Proposition 2.2 we know that the whole sequence converges with probability 1. • Remark: We obtain similar results in the case of the optimization problem (1) without restrictions, that is: G = H.

3 An Optimal Stochastic Control Problem for a Parabolic Ito Equation

In this section we consider a stochastic control problem with a stochastic parabolic differential equation as state equation. Let W (x, t) (x ∈ [0, 1], t ∈ [0, t]) be a two parameter Wiener field, that is:

• P {W (0, t) = 0} = P {W (x, 0) = 0} = 1 for all t ∈ [0, T], x ∈ [0, 1].

• The random variables R_{ij} = W (x_{i+1}+h, t_{j+1}+k)−W (x_{i}+h, t_{j++1}+k)−W (x_{i+1}+h, t_{j} +k)+W (x_{i}+h, t_{j} +k) are Gaussian distributed with expectation 0 and variance (x_{i+1} − x_{i})(t_{j+1} − t_{j}) for all h, k > 0 with x_{i+1} + h ≤ 1, t_{i+1} + k ≤ T, i = 0, …, n − 1, j = 0, … m − 1, 0 = x_{0} < x_{1} < … < x_{n} = 1, 0 = t_{0} < t_{1} < … < t_{m} = T.

Let F_{xt} be the σ-algebra wich is generated by W (u, s), u ∈ [0, x], s ∈ [0, t]. We denote a two parameter process y(x, t) adapted if y(x, t) is F_{xt} measurable for all x ∈ [0, 1], t ∈ [0, T]. We introduce now a optimal control problem:

min_{u} E \{\frac{1}{2} ∫_{Q}(y(x, t) − z_{d}(x, t))^{2} dxdt + \frac{M}{2} ∫_{Q} u(x, t)^{2} dxdt\} = min_{u} J(y, u)  (9)

where

dy(x, t) = \frac{∂^{2}y(x, t)}{∂x^{2}}dt + u(x, t)dt + σ(x, t)dW_{1}
+ σ_{2}(x, t)u(x, t)dW_{2} + σ_{3}(x, t)y(x, t)dW_{3} in Q =]0, 1[×]0, T],
y(x, 0) = y_{0}(x), x ∈]0, 1[, P-a.s.,
y(0, t) = y(1, t) = 0, P-a.s..  (10)

W_{1}, W_{2}, W_{3} denote independent two parameter Wiener fields. The admissible controls are adapted processes, where the values u(ω, x, t) are in a closed bounded convex set U. We denote U_{ad} as set of admissible controls. We suppose M > 0 and z_{d}, σ, σ_{2}, σ_{3} :]0, 1[×]0, T] → R^{1} are bounded functions. There are a lot of solution definitions for stochastich partial differential equations (see for example, [GRTU95]). Here we define the solution of (10) by an integral equation

y(x, t) = ∫_{0}^{1} y_{0}(z)p(x − z, t) dz + ∫_{0}^{t} ∫_{0}^{1} u(z, s)p(x − z, t − s) dsdz
+ ∫_{0}^{t} ∫_{0}^{1} σ(z, s)p(x − z, t − s) dW_{1}(s, z)
+ ∫_{0}^{t} ∫_{0}^{1} σ_{2}(z, s)p(x − z, t − s)u(z, s) dW_{2}(s, z)
+ ∫_{0}^{t} ∫_{0}^{1} σ_{3}(z, s)y(z, s)p(x − z, t − s) dW_{3}(s, z)

where p(x, t) is the Green function

p(x, t) = ∑_{r=0}^{∞} e^{−π^{2}r^{2}t} sin(πrx).

The case u(z, s) = 0, σ(z, s) = 0 is discussed in [DG95]. The problem (10) has a unique adapted solution y with sup{e^{−μt}E|y(x, t)|^{2} : x ∈ [0, 1], t ∈ [0, T]} < ∞. This fact can be proved by standard arguments as in [DG95], since the right hand side of (9) defines for fixed u a Lipschitz continuous function Ψ with respect to y in the Banach space of all adapted two parameter processes y with the norm

∥y∥ = sup{e^{−μt}(E|y(x, t)|^{2})^{1/2} : x ∈ [0, 1], t ∈ [0, t]}

for every conatant μ > 0. One can choose μ so that Ψ is a contraction in this space.

4 Calculation of an ε-optimal Control by a Parallel Algorithm

First we consider a general stochastic optimization problem

F (u) = Ef(u, ω) = min_{u∈G} .  (11)

We assume that the square residual term of Taylor formula of F (·) is bounded for all admissible u. We introduce for (11):

u^{k+1} := P r_{G} (u^{k} − ρ_{k} ∑_{i=1}^{N} \frac{f(u^{k} + α_{k}e_{i}, ω^{k}) − f(u^{k}, ω^{k})}{α_{k}}e_{i}) ,  (12)

where (e_{i}) is an orthogonal system in H, N a natural number, α_{k} a sufficiently small positive number and ω^{k} an independent realization of ω. The random functions

Y_{s} = ∑_{i=1}^{N} \frac{f(u^{k} + α_{k}e_{i}, ω^{k}) − f(u^{k}, ω^{k})}{α_{k}}e_{i}

fulfil (3), since Taylor formula gives

E (ρ_{k} ∑_{i=1}^{N} \frac{f(u^{k} + α_{k}e_{i}, ω^{k}) − f(u^{k}, ω^{k})}{α_{k}}e_{i}|u^{k}) = ρ_{k} ∑_{i=1}^{N} \frac{F (u^{k} + α_{k}e_{i}) − F (u^{k})}{α_{k}}e_{i}
= ρ_{k}∇F (u^{k}) + v^{k}α_{k},

where v^{k} denotes the square residual term of Taylor formula for u ∈ U_{ad}. If we have

∑_{k=0}^{∞} ρ_{k} = ∞, ∑_{k=0}^{∞}(ρ_{k})^{2} < ∞, ∑_{k=0}^{∞} ρ_{k}α_{k} < ∞,

then the algorithm is convergent in the sense of Theorem 2.3. Under these conditions we can formulate a sequential algorithm for the stochastic control problem. We approximate (9), (10) with 0 = t_{0} < t_{1} < … < t_{m} = T, 0 = x_{0} < x_{1} < … < x_{n} = 1 by the following discrete system:

min_{u} E \{\frac{1}{2} ∑_{i=1}^{n} ∑_{j=1}^{m}(y_{ij} − z_{d_{ij}})^{2} + \frac{M}{2} ∑_{i=1}^{n} ∑_{j=1}^{m} u_{ij}^{2}\} = min_{u} J_{N}(y_{··}, u)  (13)

where

y_{i,j+1} = y_{i,j} + n^{2}h (y_{i+1,j} + y_{i−1,j} − 2y_{i,j}) + hu_{ij} + nσ(x_{i}, t_{j})R_{i,j}^{1}
+ nσ_{2}(x_{i}, t_{j})u_{ij}R_{i,j}^{2} + nσ_{3}(x_{i}, t_{j})y_{ij}R_{i,j}^{3} i = 0, 1, … n − 1,  (14)

(n = 2, 3, ⋯) with y_{−1,m} = y_{n−1,m}, y_{n,m} = y_{0,m} and the initial condition y_{i,0} = y_{0}(\frac{j}{n}) and let us define R_{i,j}^{k} := W_{k}(x_{i+1}, t_{j+1}) − W_{k}(x_{i}, t_{j+1}) − W_{k}(x_{i+1}, t_{j}) + W_{k}(x_{i}, t_{j}) (k = 1, 2, 3). Let h be the time-step size. The solutions of (14) approximate the solution of (10) in the following sense:

Theorem 4.1 (see [G95]) We assume that for y_{0}(x) there exists the first derivative with respect to x. Then there is a constant k > 0 independent of L, σ(0), u_{0} and T such that for 4n^{2}h ≤ 1 it holds

E|y_{j,m} − y(j/n, mh)|^{2} ≤ kn^{−1}, 0 ≤ j ≤ n − 1, 0 ≤ mh ≤ T.

Definition 4.1 Let ε be any positive number. An adapted control u_{ε} is called ε-optimal of (9), (10) if J (y(u_{ε}), u_{ε}) − J (y(u^{∗}), u^{∗}) < ε, where y(u_{ε}) is solution of (10) for the control u = u_{ε}, y(u^{∗}) is solution of (10) for u = u^{∗}, and u^{∗} is an optimal control of (9), (10).

Theorem 4.2 Let u^{∗} be an optimal control for the original problem (9), (10) and u_{nm}^{∗} = (v_{ij}^{∗})_{i=0,...,n−1; j=0,...,m−1} an optimal control for the discrete problem (13), (14). Then ũ_{nm}^{∗} with ũ_{nm}^{∗}(x, t) = v_{ij}^{∗} for (x, t) ∈ [x_{i}, x_{i+1}[×[t_{j}, t_{j+1}[ is ε-optimal with respect to the original problem (9), (10) for every ε > 0 if n, m are sufficiently large.

Proof We have

|J (y(ũ_{nm}^{∗}), ũ_{nm}^{∗}) − J (y(u^{∗}), u^{∗})| ≤ |J (y(ũ_{nm}^{∗}), ũ_{nm}^{∗}) − J_{N} (y_{··}(u_{nm}^{∗}), u_{nm}^{∗})|
+ |−J (y(u^{∗}), u^{∗}) + J_{N} (y_{··}(u_{nm}^{∗}), u_{nm}^{∗})| .

Because of Theorem 4.1 and the fact that convergence of second moments follows from the convergence of mean square we can choose m and n in such a manner that the first summand of the last inequality is less than \frac{ε}{2}. If

−J_{N}(y_{··}(u_{nm}^{∗}), u_{nm}^{∗}) + J(y(u^{∗}), u^{∗}) ≥ 0,

then we obtain with the same consideration

J(y(u^{∗}), u^{∗}) − J_{N}(y_{··}(u_{nm}^{∗}), u_{nm}^{∗}) ≤ J(y(ũ_{nm}^{∗}), ũ_{nm}^{∗}) − J_{N}(y_{··}(u_{nm}^{∗}), u_{nm}^{∗}) ≤ \frac{ε}{2}.  (15)

If

−J_{N}(y_{··}(u_{nm}^{∗}), u_{nm}^{∗}) + J(y(u^{∗}), u^{∗}) < 0,

then we have for sufficiently large n and m

0 > J(y(u^{∗}), u^{∗}) − J_{N}(y_{··}(u_{nm}^{∗}), u_{nm}^{∗}) > J(y(u^{∗}), u^{∗}) − J_{N}(y_{··}(u^{∗}), u^{∗}) > −\frac{ε}{2},  (16)

where we define y_{ij}(u^{∗}) := y_{ij}(u^{∗}(x_{i}, t_{j})). The inequalities (15), (16) show the theorem. • (14) is a system of linear equations for y_{ij}. The random variables y_{ij} are linear functions from u_{ij}. Consequently, the goal function in (13) fulfils the conditions for (11). Then Theorem 3 can be applied and we can formulate: Sequential Algorithm

S1 Initialize control and state variables

For i:=1 to n do u[i][0] := u

y[i][0] := ỹ[i][0] := y

For k:=1 to N do

S2 Generate Gaussian random numbers For j := 1 to m do

For i := 1 to n do
R_{ij} := N(0, h\frac{1}{n})

S3 Solve the parabolic Ito equation with a finite-difference approximation

For i := 1 to m do
For j := 1 to n do

y_{j,i+1} := y_{j,i} + n^{2}h(y_{j+1,i} + y_{j−1,i} − 2y_{j,i}) + hu_{j,i} + nσ(x_{j}, t_{i})R_{j,i}^{1}

+nσ_{2}(x_{j}, t_{i})u_{ji}R_{j,i}^{2} + nσ_{3}(x_{j}, t_{i})y_{ji}R_{j,i}^{3}.

S4 Compute J_{N}(y, u)

S5 Compute new control variables For ĩ:= 1 to n do

For j̃ := 1 to m do
ũ[ĩ][j̃] := u[ĩ][j̃] + α_{k}
For i := 1 to m do
For j := 1 to n do
ỹ_{j,i+1} := ỹ_{j,i} + n^{2}h(ỹ_{j+1,i} + ỹ_{j−1,i} − 2ỹ_{j,i}) + hũ_{j,i} + nσ(x_{j}, t_{i})R_{j,i}^{1}
+nσ_{2}(x_{j}, t_{i})ũ_{ji}R_{j,i}^{2} + nσ_{3}(x_{j}, t_{i})ỹ_{ji}R_{j,i}^{3}.
Compute J_{N}(ỹ, ũ)
u[ĩ][j̃] := P r_{G} (u[ĩ][j̃] − ρ_{k} \frac{J_{N} (y,u)−J_{N} (ỹ,ũ)}{α_{k}})

S6 Perform steps S2, S3, and S4

Our goal is to design a parallel algorithm with minimized communication cost. A suitable way here is parallel computation on a ring topology (see for instance [KGGK94], [Q94]). We choose an arbitrary number p of processors denoted by P_{1}, ..., P_{p}. The processors are connected into a ring by bidirectional links, so each processor can send and receive data to and from its two direct neighbors at the same time. Processor P_{1} should be the master processor. The master distributes and collects data, and synchronizes the work of the other processors by communication routines. The simple parallelization idea is to distribute the state levels in such a way that each space stage is assigned to one worker processor P_{i} (i = 2, 3, ..., n). This is possible without loss of generality. In practice p will be less than n. Then each worker processor will be responsible to \frac{n}{p} = q, q ∈ N states with a suitable chosen n. In step S5 each processor sends a distinct value to every other processor. Thus, a well suited topology here would be a completely-connected network (each processor has a direct communication link to every other processor). But such a communicaton network is very expensive, and networks in real parallel computers are often meshes, trees, hypercubes, or rings. So we choose a ring because embedding a ring into other, real networks is quite easy. In addition, in our problem the influence of the topology chosen decreases with increasing computation cost in comparison with communication cost. The parallel algorithm can be formulated as follows. The step numbers S1, S2, ..., S6 are correseponding to the sequential algorithm above. Parallel Algorithm Master Processor P_{1}

S1 Initialize control and state variables For i:=1 to n do

u[i][0] := u y[i][0] := ỹ[i][0] := y

For k:=1 to N do

S2 Generate Gaussian random numbers For i := 1 to n do For j := 1 to m do

R_{ij} := N(0, h\frac{1}{n})

S3 Solve the parabolic Ito equation with a finite-difference approximation For i := 1 to m do For j := 1 to n do

y_{j,i+1} := y_{j,i} + n^{2}h(y_{j+1,i} + y_{j−1,i} − 2y_{j,i}) + hu_{ji} + nσ(x_{j}, t_{i})R_{j,i}^{1}

+nσ_{2}(x_{j}, t_{i})u_{ji}R_{j,i}^{2} + nσ_{3}(x_{j}, t_{i})y_{ji}R_{j,i}^{3}.

S4 Compute J_{N}(y, u)

Send J_{N}(y, u) to P_{2}, ..., P_{p} along the bidirectional ring Receive the updated control data from P_{2}, ..., P_{p}

S6 Perform S2, S3, and Compute J_{N}(y, u)

Processor P_{l} l ≠ 1 S1 Initialize control and state variables

For i:=1 to n do u[i][0] := u y[i][0] := ỹ[i][0] := y

For k:=1 to N do

S2 Generate Gaussian random numbers For i := 1 to n do For j := 1 to m do

R_{ij} := N(0, h\frac{1}{n})

S5 Update control variables Receive J_{N}(y, u) For j̃ := 1 to m do

ũ[l][j̃] := u[l][j̃] + α_{k}

For i := 1 to m do
For j := 1 to n do

ỹ_{j,i+1} := ỹ_{j,i} + n^{2}h(ỹ_{j+1,i} + ỹ_{j−1,i} − 2ỹ_{j,i}) + hũ_{ji} + nσ(ỹ_{ji})R_{j,i} +nσ_{2}(x_{j}, t_{i})ũ_{ji}R_{j,i}^{2} + nσ_{3}(x_{j}, t_{i})ỹ_{ji}R_{j,i}^{3}.

Compute J_{N}(ỹ, ũ)

u[l][j̃] := P r_{G} (u[l][j̃] − ρ_{k} \frac{J_{N} (y,u)−J_{N} (ỹ,ũ)}{α_{k}})

Data exchange of u[l][j̃] with the other processors along the ring

In the following we want to evaluate this parallelization. Performance of parallel algorithms can be measured by speedup and efficiency (see for instance [KGGK94], [Q94]). Let T_{p} be the run time of a parallel algorithm with p processors. T_{1} is the sequential run time respectively. The speedup, denoted by S_{p}, is formally defined as the ratio of T_{1} of the best sequential algorithm to T_{p} to solve the same problem, therefore S_{p} = \frac{T_{1}}{T_{p}}. Efficiency, denoted by E_{p}, is defined as the ratio of speedup to the number of processors, that is E_{p} = \frac{S_{p}}{p}. Efficiency measures the fraction of time for which processors are usefully employed. The sequential algorithm above requires computation cost in the amount of

T_{1} = n3X_{init} + N (nmR + nmY + nmJ + nmŨ + n^{2}m^{2}Y + n^{2}m^{2}J + nmU)
+ nmR + nmY + nmJ,

with N : iteration number, n : number of stages in space, m : number of time steps, R : cost of generation of Gaussian random variables, X_{init} : cost of initializing the variables in step S1, Y : cost of solving the finite difference equation, J : computation cost of J_{N}(y, u) for one stage of state and one time step, Ũ : computation cost of ũ[l][j̃] := u[l][j̃] + α_{k}, U : computation cost of u[l][j] := P r_{G} (u[l][j] − ρ_{k} \frac{J_{N} (y,u)−J_{N} (ỹ,ũ)}{α_{k}}). For the parallel algorithm we obtain, in consideration of data transport, a parallel time of

T_{p} = n3X_{init} + N (nmR + nmY + nmJ + \frac{n}{p}mŨ + \frac{n}{p}nm^{2}Y + \frac{n}{p}nm^{2}J + \frac{n}{p}mU)
+ nmR + nmY + nmJ + N \frac{n}{2}D_{t} + NnmD_{t},

with p : number of processors with n = pq, q ∈ N, D_{t} : cost required for exchanging one real number between direct connected processors. Therefore we have

S_{p} = \frac{3X_{init} + Nm (R + (1 + nm)(Y + J) + Ũ + U) + m (R + Y + J)}{3X_{init} + Nm (R + (1 + \frac{n}{p}m)(Y + J) + \frac{1}{p}(Ũ + U)) + m (R + Y + J) + N \frac{1}{2}D_{t} + NmD_{t}} ,

E_{p} = \frac{3X_{init} + Nm (R + (1 + nm)(Y + J) + Ũ + U) + m (R + Y + J)}{3pX_{init} + Nm (R + (p + nm)(Y + J) + Ũ + U) + mp (R + Y + J) + Np\frac{1}{2}D_{t} + NpmD_{t}} .

Speedup and efficiency increase with respect to N. The limits of these functions are

lim_{N→∞} S_{p} = \frac{R + (1 + nm)(Y + J) + Ũ + U}{R + (1 + \frac{n}{p}m)(Y + J) + \frac{1}{p}(Ũ + U) + \frac{1}{2m}D_{t} + D_{t}} , p, n, m fixed,

lim_{N→∞} E_{p} = \frac{R + (1 + nm)(Y + J) + Ũ + U}{pR + (p + nm)(Y + J) + Ũ + U + \frac{1}{2m}D_{t} + pD_{t}} , p, n, m fixed

respectively. We notice, that speedup and efficiency increase slightly by increasing N. Speedup and efficiency are dependent on N and the parallelization will already be efficient for small numbers of iterations. If we choose m = 4n^{2} (see Theorem 4.1) we obtain

lim_{n→∞} S_{p} = p, p fixed,
lim_{n→∞} E_{p} = 1, p fixed,

respectively. Speedup and efficiency increase by increasing the number of discretization intervals and approach to the ideal amounts p and 1 respectively. The reason is that the cost of data transport increases with O(nm) and the cost of calculation increases with O(n^{2}m^{2}). That is, computation time dominates communication time by increasing the discretization granularity. By using p = n processors we obtain increasing speedup and efficiency if we increase the number of discretization steps (it follows from the formulas to determine S_{p} and E_{p} above). This property is advantageous for practical problems which require high discretization numbers. Here an efficiency of the parallel computation close to the optimal value E_{p} = 1 is possible.

5 Numerical Results

We discuss numerical experiments for several examples of (9) and (10). Further, we assume |u(x, t)| ≤ 1. The concrete terms of the state equation are given in Table 1 .

Program Differential equation  z_{d}
Problem 1 dy(x, t) = \frac{∂^{2}y(x,t)}{∂x^{2}} dt + u(x, t) dt + dW (x, t)  1
Problem 2 dy(x, t) = \frac{∂^{2}y(x,t)}{∂x^{2}} dt + u(x, t) dt + dW (x, t)  sin(πx)e^{t}
Problem 3 dy(x, t) = \frac{∂^{2}y(x,t)}{∂x^{2}} dt + u(x, t) dt + dW (x, t) + u(x, t)dW (x, t)  1
Problem 4 dy(x, t) = \frac{∂^{2}y(x,t)}{∂x^{2}} dt + u(x, t) dt + dW_{1}(x, t) + u(x, t)dW_{2}(x, t)  1
Problem 5 dy(x, t) = \frac{∂^{2}y(x,t)}{∂x^{2}} dt + u(x, t) dt + dW (x, t) + u(x, t)dW (x, t) + y(x, t)dW  1
Problem 6 dy(x, t) = \frac{∂^{2}y(x,t)}{∂x^{2}} dt + u(x, t) dt + dW_{1}(x, t) + u(x, t)dW_{2}(x, t) + y(x, t)dW_{3} 1

Table 1: Test examples Remark The behavior of convergence depends on the step size and the initial values of control variables. If the step size is small or the initial values of control variables are ”very close” to the optimal value, then it would be possible that new steps of the quasigradient method will be overlapped from white noise and it would be impossible to achieve convergence. In general the choice of step size depends on the concrete problem. The algorithms were tested by using a two-dimensional mesh of transputers. Some results are shown in Table 2. We notice that speedup and efficiency increase with an increasing number of iterations. This effect is caused by the structure of our parallel algorithm. While the master processor performs step S6, the other processors are waiting. This waiting time will decrease relatively with increasing iteration numbers. Thus, speedup and efficiency will increase. The numerical experiments showed, in correspondence to our theoretical results, that an efficient parallelization of this strong sequential problem is possible. We intend to implement and test the parallel algorithm on more powerful parallel computers, also on such with shared memory.

Program N  Time (in sec)  Speedup Efficiency
p = 1  p = 8  S_{8}  E_{8}
Problem 1 1  157,734  28,382  5,56  0,69
5  766,849 135,554  5,66  0,71
10 1528,182 269,518  5,67  0,71
Problem 2 1 1011,017 171,081  5,91  0,74
Problem 3 1 180,4288 31,3859  5,75  0,72
5 880,2988 150,5291  5,85  0,73
10 1755,1166 299,4577  5,86  0,73
Problem 4 1 197,7717 34,3655  5,75  0,72
5 945,6071 159,4892  5,93  0,74
10 1880,3744 315,4892  5,95  0,74
Problem 5 1 202,6135 34,2177  5,92  0,74
5 991,5815 164,6543  6,02  0,75
10 1977,7990 327,7026  6,04  0,75
Problem 6 1 224,8046 40,3402  5,57  0,70
5 1059,0629 183,3727  5,77  0,72
10 2101,8013 362,1579  5,80  0,73

Table 2: Solution by stochastic gradient method

References

[BP86] Barbu, V. and Precupanu, Th. (1986). Convexity and optimization in Banach spaces. Editura Academiei, University of Romania, Bukarest.

[E76] Ermolew, J. M. (1976). Methods of stochastic programming (Russian). Glawnaja Redakzija Fiziko-Matematitscheskoi Literaturi, Moskau.

[FR75] Fleming, W. H. and Rishel, R. W. (1975). Deterministic and stochastic optimal control. Springer-Verlag, New York.

[G95] Gaines, J. G. (1995). Numerical experiments with S(P)DE’s. London Mathematical Society, Lecture Note Series 216(1995), 55-71.

[G79] Gupal, A. M. (1979). Stochastic methods of solution of non smooth extremal problems (Russian). Naukova Dumka, Kiew.

[GRTU95] Grecksch, W., Tudor, C (1995). Stochastic evolution equations. Akademie Verlag, Berlin.

[KGGK94] Kumar, V., Grama, A., Gupta A. and Karypis, G. (1994). Introduction to parallel computing. Design and Analysis of Algorithms, Benjamin/Cummings Publishing Comp., Redwood City.

[Q94] Quinn, M. J. (1994). Parallel computing. Theory and practice, 2nd ed. McGraw-Hill, New York.

[Y92] Yin, G. (1992). On H-valued stochastic approximation: Finite dimensional projections. Stochastic Analysis and Applications, 10(3), 367-377.

[ZS93] Zoran, G. and Xuemin, S. (1993). Parallel algorithms for optimal control of large scale linear systems, Springer-Verlag, London.