A Parallel Modified Lagrangian Method for an Optimal Control Problem of a Linear Distributed Stochastic System

Thomas Pohl

Johann Wolfgang Goethe - Universität Frankfurt am Main, Fachbereich Mathematik, Postfach 111932, D-60054 Frankfurt, Germany pohl@math.uni-frankfurt.de

Wilfried Grecksch

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

Holger Blaar

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

Abstract

We consider an optimal control problem for a Hilbert space valued linear stochastic evolution equation with additive noise and a quadratic goal function in integral form. We solve an infinite dimensional linear stochastic regulator problem by applying an augmented Lagrangian algorithm. Because of high amounts of computation time required by numerical solution of such problems we investigate a parallelization of this algorithm by distributing the computations of space stages over several processor nodes of a parallel computer. We obtain an efficient parallel algorithm with low communication cost by using a ring topology.

Key Words: Stochastic Optimization, Control, Parallel Algorithms

AMS Classification: 65Y05, 93E20, 93E23

1 Introduction

A high amount of computation time is often required to solve or simulate problems by mathematical models. Parallel computers are well suited to perform such expensive computations in less time by distributing necessary computations to several processor nodes, i. e. several computations or several algorithm parts can be carried out at the same time. Communication and synchronization costs of the parallel algorithm usually 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. We consider an optimal control problem for a Hilbert space valued linear stochastic evolution equation with additive noise and a quadratic goal function in integral form. Linear parabolic Ito equations with respect to white noise in time and space are examples for such evolution equations. Infinite dimensional stochastic linear regulator problems are studied by many authors. For example, stochastic maximum principles are proved in [BGR92], [G87]. It is known that the calculation of adjoint variables of the maximum principle is difficult. The dynamic programming is applied for example in [KD92], [CP78]. An optimal feedback control is determined by the solution of infinite dimensional Riccati systems. A third possibility soluting stochastic control problems consists in the interpretation as an infinite dimensional optimization problem. Then it is possible with the help of the general optimization method to construct ε-optimal methods. For example, the authors have determined ε-optimal controls by applying parallel variants of the stochastic quasigradient method [PGB99]. Here we solve an infinite dimensional linear stochastic regulator problem by applying an augmented Lagrangian algorithm. The advantage is that we can devide the problem in minimizing the controls and the states respectively. Consequently two sequences of optimization problems are constructed. 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 a ring topology can be mapped easily to current parallel computer architectures. In the second chapter we define an optimal control problem for an abstract linear parabolic Ito equation. A sequential modified Lagrangian algorithm is formulated in the third chapter. The convergence proof is given in chapter 4. In chapter 5 we propose an efficient parallel algorithm to solve the problem numerically. Finally, chapter 6 contains some numerical results and a discussion of available speedup and efficiency.

2 An Optimal Control Problem for a Linear Parabolic Ito equation

Let (V, H, V^{∗}) be a triple of rigged Hilbert spaces. A denotes a linear operator from V in V^{∗} so that

< Av, v > +α∥v∥_{V}^{2} ≤ 0  (1)

for all v ∈ V, where < Av, v > is the value of the linear continuous functional Av ∈ V^{∗} over V. Then we can say that A is on

D(A) = {v ∈ V : Av ∈ H}  (2)

an infinitesimal operator of a semigroup of contractions in H [PZ92]. D(A) is with respect to the norm ∥v∥^{2} = ∥v∥_{H}^{2} + ∥Av∥_{H}^{2} a separable Hilbert space. Next we consider on a complete probability space (Ω, F, P) a Wiener process (Z(t), F_{t})_{t∈[0,T]} with values in a separable Hilbert space K and a kernel operator Q as covariance operator. We assume that f(t), (t ∈ [0, T]) is a function with values in the space of linear operators from Q^{1/2}K in H so that f(t)Q^{1/2} is a Hilbert Schmidt operator from K in H and ∫_{0}^{T} ∥f(t)Q^{1/2}∥_{HS}^{2}dt < ∞, and the Ito integral ∫_{0}^{t} f(s)dZ(s) defines a martingale with continuous paths in D(A) with respect to F_{t} for t ∈ [0, T]. Let (U(t))_{t∈[0,T]} be an adapted D(A)-valued process with E ∫_{0}^{T} ∥U(t)∥_{D(A)}^{2}dt ≤ c_{1} where c_{1} is a given positive constant. The set of all these processes is called the set U_{ad} of admissible controls. Finally, X_{0} denotes a F_{0}-measurable D(A)-valued variable with E∥X_{0}∥_{D(A)}^{2} < ∞. Under these assumptions the controlled problem

dX(t) = (AX(t) + U(t))dt + f(t)dZ(t), X(0) = X_{0}  (3)

has a unique (with probability 1) strong solution, that is:

X(t) ∈ D(A), P − a.s., for all t ∈ [0, T],  (4)
∫_{0}^{T} ∥X(t)∥_{D(A)}^{2}dt < ∞,  (5)
X(t) = X_{0} + ∫_{0}^{t} AX(s)ds + ∫_{0}^{t} U(s)ds + ∫_{0}^{t} f(s)dZ(s),  (6)

where the equation (6) holds for all t ∈ [0, T], P -a.s. The paths of X(t) are continuous in H. This fact follows from Theorem 2.3 on page 76 in [GT95]. The strong solution of (3) can be approximated by solutions of problems with finite dimensional Wiener processes. It follows from the definition of a Hilbert space valued Wiener process that

Z(t) = ∑_{j=1}^{∞} \sqrt{λ_{j}}Z_{j}(t)e_{j}  (7)

where λ_{j} and e_{j} are the eigenvalues and eigenvectors of Q and Z_{j}(t) = \frac{1}{\sqrt{λ_{j}}} (Z(t), e_{j})_{K}. The convergence in (7) is in the sense of mean square (Proposition 4.1 on page 87 in [PZ92]). Let X^{n} be the strong solution of

dX^{n}(t) = (AX^{n}(t) + U(t))dt + f(t)dZ^{n}(t), X(0) = X_{0}  (8)

where Z^{n}(t) = ∑_{j=1}^{n} Z_{j}(t)e_{j}.

Applying the Ito formula for the norm square ∥·∥_{H}^{2} [KR79] and Gronwall lemma we can easily prove:

Lemma 2.1 It holds

lim_{n→∞}( sup_{s∈[0,t]} E∥X^{n}(s) − X(s)∥_{H}^{2} + E ∫_{0}^{t} ∥X^{n}(s) − X(s)∥_{V}^{2} ds) = 0

for all t ∈ [0, T].

The problem of optimal control of (3) consists in the minimization of

J(U, X) = E\{\frac{1}{2} ∫_{0}^{T} ∥X(t) − X_{d}∥_{H}^{2}dt + \frac{ν}{2} ∫_{0}^{T} ∥U(t)∥_{H}^{2}dt\}  (9)

with respect to U ∈ U_{ad} where X_{d} : [0, T] → H with ∫_{0}^{T} ∥X_{d}(t)∥_{H}^{2}dt < ∞ and ν > 0 are given. X denotes the solution of (3). The optimal control problem has a unique solution U ∈ U_{ad}. This fact follows from the general optimization theory (see for example, Theorem 1.3 in [LP66]).

3 A Sequential Modified Lagrangian Algorithm

The idea of the algorithm is that the restriction “state equation” is introduced by a “penalty function” into the goal function. We choose ε = 1/r (r = 1, 2, ...) and consider the problem

min_{U,X} E \{\frac{1}{2} ∫_{0}^{T} ∥X(t) − X_{d}(t)∥_{H}^{2} dt + \frac{ν}{2} ∫_{0}^{T} ∥U(t)∥_{H}^{2} dt

+ \frac{1}{2ε}E ∫_{0}^{T} ∥X(t) − X_{0} − ∫_{0}^{t} AX(s) ds − ∫_{0}^{t} f(s) dZ(s) − ∫_{0}^{t} U(s) ds∥_{H}^{2} dt\}
= min_{U,X} J^{ε}(X, U),  (10)

where the minimum is defined with respect to the admissible controls U and the adapted D(A)- valued processes (X(t))_{t∈[0,T]} with E ∫_{0}^{T} ∥X(t)∥_{D(A)}^{2} < ∞ (denotation: X ∈ H^{2}). Without restriction of the generality we assume X_{0} = 0. In the following we omit the time paramter for stochastic processes. We define the process

q^{ε} = \frac{1}{ε} (X − ∫_{0}^{t} AX ds − ∫_{0}^{t} f dZ − ∫_{0}^{t} U ds)

for admissible controls U and X ∈ H^{2}. The optimal solution of (10) is denoted by U^{ε}, X^{ε}. We now introduce the following sequential algorithm: Sequential algorithm

S1 Set initial values q_{0} ∈ H and u_{−1} ∈ U_{ad},

S2 Calculation of X_{k} ∈ H^{2}

min_{X∈H^{2}} \{\frac{1}{2}E ∫_{0}^{T} ∥X − X_{d}∥_{H}^{2} dt
+ E ∫_{0}^{T} (q_{k} , X − ∫_{0}^{t} AX ds − ∫_{0}^{t} f dZ)_{H} dt
+ \frac{1}{2ε}E ∫_{0}^{T} ∥X − ∫_{0}^{t} AX ds − ∫_{0}^{t} f dZ − ∫_{0}^{t} u_{k−1} ds∥_{H}^{2} dt\}

S3 Calculation of u_{k} ∈ U_{ad}

min_{U∈U_{ad}} \{\frac{ν}{2} E ∫_{0}^{T} ∥u∥_{H}^{2} dt + E ∫_{0}^{T} (q_{k}, u)_{H} dt
+ \frac{1}{2ε}E ∫_{0}^{T} ∥X_{k} − ∫_{0}^{t} AX_{k} ds − ∫_{0}^{t} f dZ − ∫_{0}^{t} U ds∥_{H}^{2} dt\}

S4 IF X_{k} − ∫_{0}^{t} AX_{k} ds − ∫_{0}^{t} f dZ − ∫_{0}^{t} u_{k} ds = 0 THEN STOP ELSE Calculation of q_{k+1}

q_{k+1} := q_{k} + p̃(X_{k} − ∫_{0}^{t} AX_{k} ds − ∫_{0}^{t} f dZ − ∫_{0}^{t} u_{k} ds)

GOTO S2

with 0 ≤ p_{0} ≤ p̃ ≤ r (r … penalty factor).

Theorem 3.1 Let (U^{∗}, X^{∗}) be the optimal solution of the original optimal control problem (3), (9). Then the following convergence properties hold:

lim_{k→∞} E ∫_{0}^{T} ∥U_{k} − U^{∗}∥_{H}^{2}dt = 0,

lim_{k→∞} E ∫_{0}^{T} ∥X_{k} − X^{∗}∥_{H}^{2}dt = 0

and

lim_{k→∞} E ∫_{0}^{T} ∥q_{k+1} − q_{k}∥_{H}^{2}dt = 0.

The convergence proof is given in the next section. The algorithm will be implemented by the following example: Example We choose V = Ŵ_{2}^{1}[0, 1] = {v ∈ W_{2}^{1}[0, 1] : v(0) = v(1) = 0}, H = L^{2}[0, 1], V^{∗} = W_{2}^{−1}[0, 1] and A = Δ = \frac{d^{2}}{dx^{2}} and consider

dy(x, t) = (Δy(x, t) + u(x, t))dt + f(x, t)dW (x, t)  (11)

with

y(x, 0) = y_{0}(x), y(0, t) = y(1, t) = 0,
J(u, y) = E \{\frac{1}{2} ∫_{0}^{1} ∫_{0}^{1}(y(x, t) − y_{d}(x, t))^{2}dxdt + ν ∫_{0}^{1} ∫_{0}^{1} u^{2}(x, t)dxdt\} .  (12)

Here (W (x, t))_{x,t∈[0,1]} is a Wiener field, that is:

• W : Ω × [0, 1] × [0, 1] → R^{1}, W (0, x) = W (t, 0) = W (0, 0) = 0.

• For all t_{0}, …, t_{m}, x_{0}, …, x_{n} with t_{0} < … < t_{m}, x_{0} < … < x_{n} the random variables R_{i,j} = W (x_{i+1}, t_{j+1}) − W (x_{i}, t_{j+1}) − W (x_{i+1}, t_{j}) + W (x_{i}, t_{j}) (i = 0, …, n − 1, j = 0, …, m − 1) are independent Gaussian distributed variables with expectation zero and variance (x_{i+1} − x_{i})(t_{j+1} − t_{j}).

The admissible controls are real F_{x,y}-measurable random fields u(x, r) with values in a compact set, where F_{x,y} = σ {W (s, τ) : s ≤ x, τ ≤ t}, and further, y_{d} ∈ L^{2}([0, 1] × [0, 1]). We can interprete the Wiener field (W (x, t))_{x,t∈[0,1]} as Hilbert space valued Wiener process. The eigenfunctions e_{i}(x) = \sqrt{2} sin(πix) of Δ define an orthonormal basis in L^{2}[0, 1] and the functions \overline{e}_{i}(x) = πie_{i}(x) (i = 1, 2, …) an orthormal basis in W_{2}^{−1}[0, 1]. The processes

(∫_{0}^{t} ∫_{0}^{1} e_{i}(x)dW (s, x))_{t∈[0,1]} (i = 1, 2, …)

define independent real Wiener processes (w_{i}(t))_{t∈[0,T]}. The stochastic differential in (11) is defined by a two parameter Ito integral (see [P71]). Then

Z(·, t) = ∑_{n=1}^{∞} \frac{1}{πn}w_{n}(t)\overline{e}_{n}(·)

defines a K = W_{2}^{−1}-valued Wiener process since ∑_{n=1}^{∞} \frac{1}{π^{2}n^{2}} < ∞. This fact follows again from Proposition 4.1 on page 87 in [PZ92]. Δ obviously fulfills the conditions for A. Consequently we can interprete (11) as an example for (3) if f(x, t) is such a function so that ∫_{0}^{t} f(·, t)dZ(·, t) has paths in D(Δ). Another interpretation of (11) is given by

y(x, t) = ∫_{0}^{1} y_{0}(x)G(x − z, t)dx
+ ∫_{0}^{t} ∫_{0}^{1} u(z, s)G(x − z, t − s)dzds
+ ∫_{0}^{t} ∫_{0}^{1} f(t, s)G(x − z, t − s)dW (z, s)  (13)

where G(x, t) is the Green function of the corresponding deterministic problem. It follows from [J91] that X(·, t) = y(·, t) (P -a.s., for all t ∈ [0, 1]), where X is the strong solution of (11). We approximate the solution of (11) by a time and space discretization with t_{j+1} − t_{j} = h for all j by the solutions of the difference equations

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

with

y_{−1,m} = y_{n−1,m} = 0, y_{n,m} = y_{0,m} = 0, y_{i,0} = y_{0}(\frac{i}{n}).

It follows from Theorem 1.2 in [DG95].

Lemma 3.2 Let assume 4n^{2}h ≤ 1. Then there is a positive constant C > 0 so that

E|y_{ij} − y(\frac{i}{n}, jh)|^{2} ≤ \frac{C}{n}

for i = 1, …, n − 1 and j = 1, …, m − 1.

The above algorithm will be implemented by time and space discretization of (11) with (14). The integrals in (10) are approximated by Riemann sums. Then we obtain for the steps S2 and S3:

S2 Calculation of y^{k} :

min_{y_{i,j}} EJ_{1}^{ε}(y_{i,j}, u_{i,j})

= min_{y_{i,j}} \frac{1}{nm} \{\frac{1}{2}E ∑_{i=1}^{n} ∑_{j=1}^{m}(y_{i,j} − y_{d}(x_{i}, t_{j}))^{2}
+ E ∑_{i=1}^{n} ∑_{j=1}^{m} q_{i,j}^{k} (y_{i,j+1} − y_{i,j} − n^{2}h(y_{i+1,j} + y_{i−1,j} − 2y_{i,j}) − nf_{i,j}R_{i,j})
+ \frac{1}{2ε}E ∑_{i=1}^{n} ∑_{j=1}^{m} (y_{i,j+1} − y_{i,j} − n^{2}h(y_{i+1,j} + y_{i−1,j} − 2y_{i,j}) − hu_{i,j}^{k} − nf_{i,j}R_{i,j})^{2}\}

S3 Calculation of F_{x_{i},t_{j}}-measurable U-valued random variables u_{i,j} from L^{2}(Ω):

min_{u_{i,j}} J_{2}^{ε}(y_{i,j}, u_{i,j})

= min_{u_{i,j}} \frac{1}{nm} \{\frac{ν}{2} E ∑_{i=1}^{n} ∑_{j=1}^{m} u_{i,j}^{2} + E ∑_{i=1}^{n} ∑_{j=1}^{m} q_{i,j}^{k}u_{i,j}
+ \frac{1}{2ε}E ∑_{i=1}^{n} ∑_{j=1}^{m} (y_{i,j+1}^{k} − y_{i,j}^{k} − n^{2}h(y_{i+1,j}^{k} + y_{i−1,j}^{k} − 2y_{i,j}^{k}) − hu_{i,j} − nf_{i,j}R_{i,j})^{2}\} .

The optimization problems in S2 ans S3 can be solved, for example, by a stochastic quasigradient method [E76].

4 Convergence Proof

At first we will show that the sequences (U^{ε}), (X^{ε}), (q^{ε}) are bounded in the following sense:

Lemma 4.1 There are k_{0} > 0 and C > 0 so that for a given ε_{0} we have

sup_{0<ε≤ε_{0}} E ∥ U^{ε} ∥_{H}^{2}≤ k_{0},
sup_{0<ε≤ε_{0}} E ∥ X^{ε} ∥_{H}^{2}≤ k_{0},
sup_{0<ε≤ε_{0}} E \frac{1}{ε} ∥ X^{ε} − ∫_{0}^{t} AX^{ε} ds − ∫_{0}^{t} f dZ + ∫_{0}^{t} U^{ε} ds ∥_{H}^{2}≤ k_{0}
sup_{0<ε≤ε_{0}} E ∥ q_{ε} ∥_{H}^{2}≤ C.

Proof It holds

0 ≤ J(X^{ε}, U^{ε}) ≤ J^{ε}(X^{ε}, U^{ε}) ≤ J^{ε}(X^{∗}, U^{∗}) = J(X^{∗}, U^{∗}) < ∞,  (15)

where the last equation is valid since we have q^{ε} = 0 for X = X^{∗}. Consequently, we get

E ∫_{0}^{1} ∥ X^{ε} − X_{d} ∥_{H}^{2} dt + νE ∫_{0}^{1} ∥ U^{ε} ∥_{H}^{2} dt
+\frac{1}{ε}E ∫_{0}^{1} ∫_{0}^{1} ∥ X^{ε} − ∫_{0}^{t} X^{ε} ds − ∫_{0}^{t} f dZ + ∫_{0}^{t} U^{ε} ds ∥_{H}^{2} dt ≤ J^{∗}.  (16)

Obviously the first three a-priori estimates follow from (16). To prove the fourth a-priori estimate we choose χ ∈ L_{H}^{2}[0, 1] : ∫_{0}^{1} ∥ χ ∥_{H}^{2}≤ 1, p̃ > 0, X_{χ} = T (u − p̃χ) ∈ L_{H}^{2}([0, 1] × Ω), where T : L_{H}^{2}([0, 1] × Ω) → L_{H}^{2}([0, 1] × Ω) denotes the operator which is defined by the right hand side of the stochastic differential equation (3). It follows from the necessary optimality criterion

J^{ε′}(X^{ε}, U^{ε})(X − X^{ε}, u − U^{ε}) ≥ 0,
∀X ∈ H^{2}, ∀u ∈ U_{ad},

where

J^{ε′}(X^{ε}, U^{ε})(X − X^{ε}, u − U^{ε})
= E ∫_{0}^{1} (X^{ε} − X_{d} , X − X^{ε})_{H} dt + νE ∫_{0}^{1} (U^{ε} , u − U^{ε})_{H} dt
− E ∫_{0}^{1} (q_{ε} , ∫_{0}^{t}(u − U^{ε})ds)_{H} dt + E ∫_{0}^{1} (q_{ε}, X − X^{ε} − ∫_{0}^{t} A(X − X^{ε})ds)_{H} dt

and it holds for the control u_{χ} = u − p̃χ

E ∫_{0}^{1} (X^{ε} − X_{d} , X_{χ} − X^{ε})_{H} dt + νE ∫_{0}^{1} (U^{ε} , u_{χ} − U^{ε})_{H} dxdt
− E ∫_{0}^{1} (q^{ε}, ∫_{0}^{t}(u_{χ} − U^{ε})ds)_{H} dt + E ∫_{0}^{1} (q^{ε} , X_{χ} − X^{ε} − ∫_{0}^{t} A(X_{χ} − X^{ε})ds)_{H} dt
= E ∫_{0}^{1} (X^{ε} − X_{d} , X_{χ} − X^{ε})_{H} dt + νE ∫_{0}^{1} (U^{ε} , u_{χ} − U^{ε})_{H} dt
− E ∫_{0}^{1} (q^{ε} , ∫_{0}^{t}(u_{χ} − U^{ε})ds)_{H} dt + E ∫_{0}^{1} (q^{ε}, ∫_{0}^{t}(u_{χ} − p̃χ) ds − ∫_{0}^{t} U^{ε} ds − εq^{ε})_{H} dt
= E ∫_{0}^{1} (X^{ε} − X_{d} , X_{χ} − X^{ε})_{H} dt + νE ∫_{0}^{1} (u^{ε} , u_{χ} − U^{ε})_{H} dt
+ E ∫_{0}^{1} (q^{ε}, ∫_{0}^{t}(−p̃χ) ds − εq^{ε})_{H} dt ≥ 0.

Consequently

E ∫_{0}^{1} (q^{ε} , p̃∫_{0}^{t} χ ds)_{H} dt ≤ νE ∥ U^{ε} ∥_{L^{2}}∥ u_{χ} − U^{ε} ∥_{L^{2}}
+ E ∥ X^{ε} − X_{d} ∥_{L^{2}}∥ X_{χ} − X^{ε} ∥_{L^{2}} −εE ∫_{0}^{1} ∥q^{ε}∥_{H}^{2} ≤ C,

where C is a constant independently from ε. If we consider the supremum over χ then the definition of the norm of a linear continuous functional over a Hilbert space gives the assertion. • Proof of Theorem 3.1 It follows from the necessary optimality criterion for the problem (9)

J′(X^{∗}, U^{∗})(X − X^{∗}, U − U^{∗}) ≥ 0,
∀X ∈ H^{2}, ∀U ∈ U_{ad},

where

J′(X^{∗}, U^{∗})(X − X^{∗}, U − U^{∗})
= E ∫_{0}^{1}(X^{∗} − X_{d} , X − X^{∗})_{H}dt + Eν ∫_{0}^{1}(U^{∗} , U − U^{∗})_{H}dt.

If we note

X^{∗} − ∫_{0}^{t} AX^{∗}ds − ∫_{0}^{t} fdZ − ∫_{0}^{t} U^{∗}ds = 0,

then the necessary optimality criterion gives the following inequalities, where the process q^{∗} is the process q for U = U^{∗}, X = X^{∗}:

E ∫_{0}^{1}(X^{∗} − X_{d} , X − X^{∗})_{H}dt + E ∫_{0}^{1} (q^{∗} , X − X^{∗} − ∫_{0}^{t} A(X − X^{∗})ds)_{H} dt
+ \frac{1}{ε}E ∫_{0}^{1} (X^{∗} − ∫_{0}^{t} AX^{∗}ds − ∫_{0}^{t} fdZ − ∫_{0}^{t} U^{∗}ds , X − X^{∗} − ∫_{0}^{t} A(X − X^{∗})ds)_{H} dt ≥ 0 (17)

and

νE ∫_{0}^{1}(U^{∗} , U − U^{∗})_{H}dt − E ∫_{0}^{1} (q^{∗} , ∫_{0}^{t}(U − U^{∗})ds)_{H} dt
− \frac{1}{ε}E ∫_{0}^{1} (X^{∗} − ∫_{0}^{t} AX^{∗}ds − ∫_{0}^{t} fdZ − ∫_{0}^{t} U^{∗}ds , ∫_{0}^{t}(U − U^{∗})ds)_{H} dt ≥ 0. (18)

Calculating the gradients of the functionals of steps S2 and S3 of the continuous sequential algorithm in the points X − X_{n} and U − u_{n}, then

E ∫_{0}^{1}(X_{n} − X_{d} , X − X_{n})_{H}dt + E ∫_{0}^{1} (q_{n} , X − X_{n} − ∫_{0}^{t} A(X − X_{n})ds)_{H} dt

+ \frac{1}{ε}E ∫_{0}^{1} (X_{n} − ∫_{0}^{t} AX_{n}ds − ∫_{0}^{t} fdZ − ∫_{0}^{t} u_{n−1}ds, X − X_{n} − ∫_{0}^{t} A(X − X_{n})ds)_{H} dxdt ≥ 0

(19)

and

νE ∫_{0}^{1}(u_{n} , U − u_{n})_{H}dt − E ∫_{0}^{1} (q_{n} , ∫_{0}^{t}(U − u_{n})ds)_{H} dt
− \frac{1}{ε}E ∫_{0}^{1} (X_{n} − ∫_{0}^{t} AX_{n}ds − ∫_{0}^{t} fdZ − ∫_{0}^{t} u_{n}ds , ∫_{0}^{t}(U − u_{n})ds)_{H} dt ≥ 0. (20)

Next we substitute X = X_{n} in (17), X = X^{∗} in (19), U = u_{n} in (18), and U = U^{∗} in (20). Then we mutiply these inequalities by −1 and add these inequalities, and we get

E ∫_{0}^{1}(−X^{∗} + X_{d} , X_{n} − X^{∗})_{H}dt − E ∫_{0}^{1} (q^{∗} , X_{n} − X^{∗} − ∫_{0}^{t} A(X_{n} − X^{∗})ds)_{H} dt

− \frac{1}{ε}E ∫_{0}^{1} (X^{∗} − ∫_{0}^{t} AX^{∗}ds − ∫_{0}^{t} fdZ − ∫_{0}^{t} U^{∗}ds , X_{n} − X^{∗} − ∫_{0}^{t} A(X_{n} − X^{∗})ds)_{H} dt
+ E ∫_{0}^{1}(X_{n} − X_{d} , −X^{∗} + X_{n})_{H}dt + E ∫_{0}^{1} (q_{n}, −X^{∗} + X_{n} − ∫_{0}^{t} A(−X^{∗} + X_{n})ds)_{H} dt
+ \frac{1}{ε}E ∫_{0}^{1} (X_{n} − ∫_{0}^{t} AX_{n}ds − ∫_{0}^{t} fdZ − ∫_{0}^{t} u_{n−1}ds , −X^{∗} + X_{n} − ∫_{0}^{t} A(−X^{∗} + X_{n})ds)_{H} dt
− Eν ∫_{0}^{1}(U^{∗}, u_{n} − U^{∗})_{H}dt + E ∫_{0}^{1} (q^{∗} , ∫_{0}^{t}(u_{n} − U^{∗})ds)_{H} dt
− \frac{1}{ε}E ∫_{0}^{1} (X^{∗} − ∫_{0}^{t} AX^{∗}ds − ∫_{0}^{t} fdZ − ∫_{0}^{t} U^{∗}ds , ∫_{0}^{t}(−u_{n} + U^{∗})ds)_{H} dt
+ Eν ∫_{0}^{1}(u_{n} , −U^{∗} + u_{n})_{H}dt − E ∫_{0}^{1} (q_{n} , ∫_{0}^{t}(−U^{∗} + u_{n})ds)_{H} dt
+ \frac{1}{ε}E ∫_{0}^{1} (X_{n} − ∫_{0}^{t} AX_{n}ds − ∫_{0}^{t} fdZ − ∫_{0}^{t} u_{n}ds , ∫_{0}^{t}(u^{∗} − u_{n})ds)_{H} dt ≤ 0.

Elementary calculations of the last inequality give

E ∫_{0}^{1} ∥ X_{n} − X^{∗} ∥_{H}^{2} +νE ∫_{0}^{1} ∥ u_{n} − U^{∗} ∥_{H}^{2}
+ E ∫_{0}^{1} (q_{n} − q^{∗} , X_{n} − X^{∗} − ∫_{0}^{t} A(X_{n} − X^{∗})ds − ∫_{0}^{t}(u_{n} − U^{∗})ds)_{H} dt
+ \frac{1}{ε}E ∫_{0}^{1} ∥ X_{n} − X^{∗} − ∫_{0}^{t} A(X_{n} − X^{∗})ds ∥_{H}^{2} dt + \frac{1}{ε}E ∫_{0}^{1} ∥ ∫_{0}^{t}(U^{∗} − u_{n})ds ∥_{H}^{2} dt
+ \frac{1}{ε}E ∫_{0}^{1} (X_{n} − X^{∗} − ∫_{0}^{t} A(X_{n} − X^{∗})ds , ∫_{0}^{t}(2U^{∗} − u_{n−1} − u_{n})ds)_{H} dt ≤ 0. (21)

Obviously

E ∫_{0}^{1} (X_{n} − X^{∗} − ∫_{0}^{t} A(X_{n} − X^{∗})ds , ∫_{0}^{t}(2U^{∗} − u_{n−1} − u_{n})ds)_{H} dxdt
= E ∫_{0}^{1} ∥ X_{n} − X^{∗} − ∫_{0}^{t} A(X_{n} − X^{∗})ds − ∫_{0}^{t}(u_{n} − U^{∗})ds ∥_{H}^{2} dt
− E ∫_{0}^{1} ∥ X_{n} − X^{∗} − ∫_{0}^{t} A(X_{n} − X^{∗})ds ∥_{H}^{2} dt − E ∫_{0}^{1} ∥ ∫_{0}^{t}(U^{∗} − u_{n})ds ∥_{H}^{2} dt
+ E ∫_{0}^{1} (X_{n} − X^{∗} − ∫_{0}^{t} A(X_{n} − X^{∗})ds , ∫_{0}^{t}(u_{n} − u_{n−1})ds)_{H} dt.

If we substitute the last term in the last summand of the left hand side of (21) then

E ∫_{0}^{1} ∥ X_{n} − X^{∗} ∥_{H}^{2} dt + νE ∫_{0}^{1} ∥ u_{n} − U^{∗} ∥_{H}^{2} dt
+ E ∫_{0}^{1} (q_{n} − q^{∗} , X_{n} − X^{∗}) − ∫_{0}^{t} A(X_{n} − X^{∗})ds − ∫_{0}^{t}(u_{n} − U^{∗})ds)_{H} dt
+ \frac{1}{ε}E ∫_{0}^{1} ∥ X_{n} − X^{∗} − ∫_{0}^{t} A(X_{n} − X^{∗})ds − ∫_{0}^{t}(u_{n} − u^{∗})ds ∥_{H}^{2} dt
+ \frac{1}{ε}E ∫_{0}^{1} (X_{n} − X^{∗} − ∫_{0}^{t} A(X_{n} − X^{∗})ds , ∫_{0}^{t}(u_{n} − u_{n−1}) ds)_{H} dt ≤ 0.  (22)

One obtains with the substitutions X̃_{j} := X_{j} − X^{∗} and ũ_{j} := u_{j} − U^{∗} (j = 1, 2, …, n)

E ∫_{0}^{1} (X_{n} − X^{∗} − ∫_{0}^{t} A(X_{n} − X^{∗})ds , ∫_{0}^{t}(u_{n} − u_{n−1})ds)_{H} dt

= E ∫_{0}^{1} (X̃_{n} − X̃_{n−1} − ∫_{0}^{t} A(X̃_{n} − X̃_{n−1})ds , ∫_{0}^{t}(ũ_{n} − ũ_{n−1})ds)_{H} dt
+ E ∫_{0}^{1} (X̃_{n−1} − ∫_{0}^{t} AX̃_{n−1}ds − ∫_{0}^{t} ũ_{n−1}ds , ∫_{0}^{t}(ũ_{n} − ũ_{n−1})ds)_{H} dt
+ \frac{1}{2}E ∫_{0}^{1} [∥ ∫_{0}^{t} ũ_{n}ds ∥_{H}^{2} − ∥ ∫_{0}^{t} ũ_{n−1} ds ∥_{H}^{2} − ∥ ∫_{0}^{t}(ũ_{n} − ũ_{n−1}) ds ∥_{H}^{2}] dt. (23)

First we substitute in (20) n by n − 1 and u by u_{n} and second we substitute in (20) u by u_{n−1} and then we add these inequalities. We obtain with easy transformations

νE ∫_{0}^{1} ∥ u_{n} − u_{n−1} ∥_{H}^{2} dt − E ∫_{0}^{1} (q_{n} − q_{n−1} , ∫_{0}^{t}(u_{n} − u_{n−1}) ds)_{H} dt
+ \frac{1}{ε}E ∫_{0}^{1} ∥ ∫_{0}^{t}(u_{n} − u_{n−1}) ds ∥_{H}^{2} dt
− \frac{1}{ε}E ∫_{0}^{1} (X_{n} − X_{n−1} − ∫_{0}^{t} A(X_{n} − X_{n−1})ds , ∫_{0}^{t}(u_{n} − u_{n−1}) ds)_{H} dt ≤ 0.

Consequently

\frac{1}{ε}E ∫_{0}^{1} (X̃_{n} − X̃_{n−1} − ∫_{0}^{t} A(X̃_{n} − X̃_{n−1})ds , ∫_{0}^{t}(ũ_{n} − ũ_{n−1})ds)_{H} dt

≥ νE ∫_{0}^{1} ∥ ũ_{n} − ũ_{n−1} ∥_{H}^{2} dt + \frac{1}{ε}E ∫_{0}^{1} ∥ ∫_{0}^{t}(ũ_{n} − ũ_{n−1}) ds ∥_{H}^{2} dt
− p̃E ∫_{0}^{1} (X_{n−1} − ∫_{0}^{t} Ay_{n−1}ds − ∫_{0}^{t} fdZ − ∫_{0}^{t} u_{n−1} ds , ∫_{0}^{t}(ũ_{n} − ũ_{n−1}) ds)_{H} dt (24)

where the difference q_{n} − q_{n−1} was represented as in step S4. Applying in (23) the inequality (24) and

X^{∗} − ∫_{0}^{t} AX^{∗}ds − ∫_{0}^{t} fdZ − ∫_{0}^{t} U^{∗}ds = 0

and the definitions of X̃_{n−1} and ũ_{n−1}, then

\frac{1}{ε}E ∫_{0}^{1} (X_{n} − X^{∗} − ∫_{0}^{t} A(X_{n} − X^{∗})ds , ∫_{0}^{t}(u_{n} − u_{n−1}) ds)_{H} dt

≥ νE ∫_{0}^{1} ∥ u_{n} − u_{n−1} ∥_{H}^{2} dt + \frac{1}{2ε}E ∫_{0}^{1} ∥ ∫_{0}^{t}(u_{n} − u_{n−1}) ds ∥_{H}^{2} dt
+ \frac{1}{2ε}E [∫_{0}^{1} ∥ ∫_{0}^{t}(u_{n} − u^{∗}) ds ∥_{H}^{2} dt − ∫_{0}^{1} ∥ ∫_{0}^{t}(u_{n−1} − u^{∗}) ds ∥_{H}^{2} dt]
+ (\frac{1}{ε} − p̃)E ∫_{0}^{1} (X_{n−1} − ∫_{0}^{t} AX_{n−1}ds − ∫_{0}^{t} fdZ − ∫_{0}^{t} u_{n−1} ds , ∫_{0}^{t}(ũ_{n} − ũ_{n−1}) ds)_{H} dt
≥ νE ∫_{0}^{1} ∥ u_{n} − u_{n−1} ∥_{H}^{2} dt + \frac{p̃}{2}E ∫_{0}^{1} ∥ ∫_{0}^{t}(u_{n} − u_{n−1}) ds ∥_{H}^{2} dt
+ \frac{1}{2ε}E ∫_{0}^{1} [∥ ∫_{0}^{t}(u_{n} − u^{∗}) ds ∥_{H}^{2} − ∥ ∫_{0}^{t}(u_{n−1} − u^{∗}) ds ∥_{H}^{2}] dt
+ \frac{p̃ − \frac{1}{ε}}{2} E ∫_{0}^{1} ∥ X_{n−1} − ∫_{0}^{t} AX_{n−1}ds − ∫_{0}^{t} fdZ − ∫_{0}^{t} u_{n−1} ds ∥_{H}^{2} .

Consequently one obtains from (22)

E ∫_{0}^{1} ∥ X_{n} − X^{∗} ∥_{H}^{2} dt + ν ∫_{0}^{1} E ∥ u_{n} − u^{∗} ∥_{H}^{2} dt

+ E ∫_{0}^{1} (q_{n} − q^{∗} , X_{n} − ∫_{0}^{t} AX_{n} ds − ∫_{0}^{t} u_{n} ds − ∫_{0}^{t} f dZ)_{H} dt
+ \frac{1}{ε}E ∫_{0}^{1} ∥X_{n} − ∫_{0}^{t} AX_{n} ds − ∫_{0}^{t} u_{n} ds − ∫_{0}^{t} f dZ∥_{H}^{2} dt
+ νE ∫_{0}^{1} ∥ u_{n} − u_{n−1} ∥_{H}^{2} dt + \frac{p̃}{2} E ∫_{0}^{1} ∥ ∫_{0}^{t}(u_{n} − u_{n−1}) ds ∥_{H}^{2} dt
+ \frac{1}{2ε} E ∫_{0}^{1} ∥ ∫_{0}^{t}(u_{n} − U^{∗}) ds ∥_{H}^{2} dt − \frac{1}{2ε} E ∥ ∫_{0}^{t}(u_{n−1} − U^{∗}) ds ∥_{H}^{2} dt
+ \frac{p̃ − \frac{1}{ε}}{2} E ∫_{0}^{1} ∥ X_{n−1} − ∫_{0}^{t} AX_{n−1}ds − ∫_{0}^{t} fdZ − ∫_{0}^{t} X_{n−1} ds ∥_{H}^{2} dt ≤ 0.

Let’s introduce

γ_{n} = \frac{1}{2ε} ∫_{0}^{1} ∥ ∫_{0}^{t}(u_{n−1}−U^{∗}) ds ∥_{H}^{2} dt+ \frac{\frac{1}{ε} − p̃}{2} ∫_{0}^{1} ∥ X_{n−1}−∫_{0}^{t} AX_{n−1}ds−∫_{0}^{t} fdZ−∫_{0}^{t} u_{n−1} ds ∥_{H}^{2},

then

E ∫_{0}^{1} ∥ X_{n} − X^{∗} ∥_{H}^{2} dt + νE ∫_{0}^{1} ∥ u_{n} − U^{∗} ∥_{H}^{2} dt
+ \frac{\frac{1}{ε} + p̃}{2} E ∫_{0}^{1} ∥X_{n} − ∫_{0}^{t} AX_{n} ds − ∫_{0}^{t} u_{n} ds − ∫_{0}^{t} f dZ∥_{H}^{2} dt
+ νE ∫_{0}^{1} ∥ u_{n} − u_{n−1} ∥_{H}^{2} dt + \frac{p̃}{2} E ∫_{0}^{1} ∥ ∫_{0}^{t}(u_{n} − u_{n−1}) ds ∥_{H}^{2} dt
≤ E γ_{n} − E γ_{n+1} − E ∫_{0}^{1} (q_{n} − q^{∗} , X_{n} − ∫_{0}^{t} AX_{n} ds − ∫_{0}^{t} u_{n} ds − ∫_{0}^{t} f dZ)_{H} dt. (25)

The definition of q in step S4 shows

E ∫_{0}^{1} (q_{n} − q^{∗} , X_{n} − ∫_{0}^{t} AX_{n} ds − ∫_{0}^{t} u_{n} ds − ∫_{0}^{t} f dZ)_{H} dt

= \frac{1}{p̃} E ∫_{0}^{1}(q_{n} − q^{∗} , q_{n+1} − q_{n})_{H}dt,
q_{n+1} = \frac{1}{p̃} E [−\frac{1}{2} ∫_{0}^{1} ∥q_{n}∥_{H}^{2}dt + \frac{1}{2} ∫_{0}^{1} ∥q_{n+1}∥_{H}^{2}dt − \frac{1}{2} ∫_{0}^{1} ∥q_{n} − q_{n+1}∥_{H}^{2}dt
− ∫_{0}^{1}(q^{∗} , q_{n+1})_{H}dt + ∫_{0}^{1}(q^{∗} , q_{n})_{H}dt] .

Substituting

μ_{n} = γ_{n} + \frac{1}{2p̃} ∫_{0}^{1} ∥ q_{n}∥_{H}^{2}dt − ∫_{0}^{1}(q^{∗}q_{n})_{H}dt

then it follows from (25)

E ∫_{0}^{1} ∥ X_{n} − X^{∗} ∥_{H}^{2} dt + νE ∫_{0}^{1} ∥ u_{n} − u^{∗} ∥_{H}^{2} dt

+ \frac{\frac{1}{ε} + p̃}{2} E ∫_{0}^{1} ∥X_{n} − ∫_{0}^{t} AX_{n} ds − ∫_{0}^{t} u_{n} ds − ∫_{0}^{t} f dZ∥_{H}^{2} dt
≤ E μ_{n} − E μ_{n+1} + \frac{p̃}{2}E ∫_{0}^{1} ∥X_{n} − ∫_{0}^{t} AX_{n} ds − ∫_{0}^{t} u_{n} ds − ∫_{0}^{t} f dZ∥_{H}^{2} dt.  (26)

Finally one obtains

0 ≤ E ∫_{0}^{1} ∥ X_{n} − X^{∗} ∥_{H}^{2} dt + νE ∫_{0}^{1} ∥ u_{n} − U^{∗} ∥_{H}^{2} dt

+ \frac{1}{2ε} E ∫_{0}^{1} ∥X_{n} − ∫_{0}^{t} AX_{n} ds − ∫_{0}^{t} u_{n} ds − ∫_{0}^{t} f dZ∥_{H}^{2} dt  (27)
≤ E μ_{n} − E μ_{n+1}.

Consequently, the sequence (E μ_{n}) is monotone non-increasing. Because of γ_{n} ≥ 0 and

E μ_{n} ≥ \frac{1}{2p̃}E [∫_{0}^{1} ∥ q_{n}∥_{H}^{2}dt − 2 ∫_{0}^{1}(q^{∗} , q_{n})_{H}dt] ≥ −\frac{1}{2p̃} E ∫_{0}^{1} ∥q^{∗}∥_{H}^{2}dt

it follows E (μ_{n}−μ_{n+1}) → 0 for n → ∞. The inequality (27) shows that the theorem is proved. •

Remark 4.1 The convergence Theorem 3.1 is also true if we consider the equation (11) in the interpretation (13). Then we substitute in the proof above AX, AX_{u}, AX^{∗} by ∫_{0}^{1} y_{0}(x)G(x − t, t)dx, ∫_{0}^{t} u(s)ds by ∫_{0}^{t} ∫_{0}^{1} u(z, s)G(x−z, t−s)dzds, ∫_{0}^{t} u_{n}(s)ds by ∫_{0}^{t} ∫_{0}^{1} u(z, s)G(x−z, t−s)dzds, ∫_{0}^{t} u^{∗}(s)ds by ∫_{0}^{t} ∫_{0}^{1} u^{∗}(z, s)G(x − z, t − s)dzds and ∫_{0}^{t} f(s)dW (s) by ∫_{0}^{t} ∫_{0}^{1} u(z, s)G(x − z, t − s)dW (dz, ds). Then we especially see that the convergence properties also hold if we consider the case f(x, t) ≡ 1. Obviously in this case there does not exist a strong solution in the sense of (3) but a solution in the sense of (13).

5 A Parallel Algorithm

The algorithm above is numerically expensive. Simulating Wiener fields, solving stochastic differential equations, and computing integrals require a high amount of computing time. Thus, a parallelization could be of advantage. Our goal is to distribute computations on several processors in such a way that the communication overhead will be small, i. e. after long computing parts only small amounts of data communication should occur. This is possible by distributing the computations of n space stages to several processors. Each of p processors computes a new state in space if p = n, or, if p < n, each processor has to compute several new states in space. Without loss of generality we can assume that n = pl, l ∈ N holds. Our parallel algorithm is designed for parallel computers with distributed memory (so called message passing systems). Here the minimization of communication cost of an algorithm is especially important. The nature of our problem shows that a parallelization over the space requires only data communication between neighbors owing to space discretization. Each processor P_{i}, i = 1, ..., p, has to change data only with two other processors, that compute the space states i − 1 and i + 1 (in the case of p = n, accordingly if n = pl, l ∈ N). Only in the case of y processor P_{i} has to send and to receive values of space states i − 2 and i + 2 in addition. This is easy to carry out by a second communication between neighbors after the first one, i. e. after the values of states i − 1 and i + 1 has reached their neighbors. Hence, a ring (see Figure 1) is a very suitable topology here, and embedding a ring topology into real networks of parallel computers is quite easy. A parallel algorithm, the same for each processor, can be formulated as follows.

Parallel Algorithm on processor P_{i}

S1 Initialize control and state variables

For j:=1 to m do

q[i][j] := q̃; u[i][j] := ũ; y[i][j] := ỹ iteration := true

While (iteration = true)

S2 Send (right, q_{i}); Recv(left, q_{i−1})

Send (left, q_{i}); Recv(right, q_{i+1}) Send (right, u_{i}); Recv(left, u_{i−1}) Send (left, u_{i}); Recv(right, u_{i+1}) For k:=1 to N do

For j:=1 to m do

R[i][j] := N(0, \frac{h}{n}) Send (right, R_{i}); Recv (left, R_{i−1}) Send (left, R_{i}); Recv (right, R_{i+1}) Send (right, y_{i}); Recv (left, y_{i−1}) Send (left, y_{i}); Recv (right, y_{i+1}) Send (right, y_{i−1}); Recv (left, y_{i−2}) Send (left, y_{i+1}); Recv (right, y_{i+2}) For j:=0 to m do

y[i][j] := y[i][j] − ρ_{k}grad_{y[i][j]}J^{ε}

S3 Send (right, y_{i}); Recv(left, y_{i−1})

Send (left, y_{i}); Recv(right, y_{i+1}) For k:=1 to N do

For j:=1 to m do

R[i][j] := N(0, \frac{h}{n}) Send (right, u_{i}); Recv(left, u_{i−1}) Send (left, u_{i}); Recv(right, u_{i+1}) For j:=1 to m do

u[i][j] := u[i][j] − ρ_{k}grad_{u[i][j]}J^{ε}

Send (right, u_{i}); Recv(left, u_{i−1}) Send (left, u_{i}); Recv(right, u_{i+1})

S4 If (y[i][j] − (y[i + 1][j] + y[i − 1][j] − 2y[i][j]) − hũ[i][j] − nf[i][j]R[i][j]) < eps)

iteration := false else

q[i][j] := q[i][j]+p̃(y[i][j +1]−y[i][j]−n^{2}h(y[i+1][j]+y[i−1][j]−2y[i][j])− hũ[i][j] − nf[i][j]R[i][j])

with N = iteration number, n = number of space stages, and m = number of time steps. Send and Recv declare data exchanges between neighboring processors within the ring. R_{i} denotes the vector R[i][j], j = 1, 2, ..., m, y_{i}, u_{i}, and q_{i} correspondingly.

processor P_{1}

processor P_{2}

…

processor P_{i}

processor P_{i+1}

processor P_{p}

processor P_{p−1}

…

processor P_{i+3} processor P_{i+2}

1.

1.

2.

2.

Figure 1: Ring topology

Before looking at numerical experiments 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} = nm(Q_{init} + Y_{init} + U_{init}) + (M − 1) [Nnm(2R + Y_{new} + U_{new}) + nmA + nmQ_{new}]
+ Nnm(2R + Y_{new} + U_{new}) + nmA,

with

Q_{init}, U_{init}, Y_{init}: cost of initializing the variables in step S1,

R: cost of generation of a Gaussian random variable,

U_{new}, Y_{new}: cost of computing new values of u and y by one iteration with a stochastic gradient method,

Q_{new}: cost of computing a new value of q

A: computation cost of termination condition

M: number of running the while-loop (S2, S3, S4), until termination condition is true

For the parallel algorithm we obtain, in consideration of data transport, a parallel time of T_{p} = \frac{n}{p}m(Q_{init} + Y_{init} + U_{init}) + (M − 1) [N · \frac{n}{p}m(2R + Y_{new} + U_{new}) + \frac{n}{p}mA + \frac{n}{p}mQ_{new}]

+ N · \frac{n}{p}m(2R + Y_{new} + U_{new}) + \frac{n}{p}mA + 8MD_{t_{m}} + 8MND_{t_{m}} with

p : number of processors with n = pl, l ∈ N,

D_{t_{m}}: cost required for exchanging a block of m real numbers between direct connected processors (D_{t_{m}} = t_{s} + mt_{w}, t_{s} = startup time, t_{w} = per-word transfer time),

Therefore we have S_{p} = \frac{(Q_{init}+Y_{init}+U_{init})+MN(2R+Y_{new}+U_{new})+MA+(M−1)Q_{new}}{\frac{1}{p} (Q_{init}+Y_{init}+U_{init})+ \frac{MN}{p} (2R+Y_{new}+U_{new})+ \frac{M}{p} A+ \frac{M−1}{p} Q_{new}+ \frac{8M}{n} ( \frac{ts}{m} +t_{w})(1+N)},

E_{p} = _{\frac{(Q_{init}+Y_{init}+U_{init})+MN(2R+Y_{new}+U_{new})+MA+(M−1)Q_{new}}{(Q_{init}+Y_{init}+U_{init})+MN(2R+Y_{new}+U_{new})+MA+(M−1)Q_{new}+ \frac{8Mp}{n} ( \frac{ts}{m} +t_{w})(1+N)}} .

Now we want to have a look at the behavior of speedup and efficiency of our parallel algorithm if we increase the iteration number N and the discretization granularity n. To increase n may be necessary if higher computing accuracy of the solution is required. By increasing n we obtain from S_{p} and E_{p} above

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

We see that speedup and efficiency increase by increasing the number of discretization intervals and approach to the ideal amounts p and 1 respectively. But this result is only interesting if high computing accuracy is required (T_{1} and T_{p} increase absolutely, too). If higher accuracy of the stochastic gradient method is required we consider

lim_{N→∞} S_{p} = \frac{M(2R + Y_{new} + U_{new})}{\frac{M}{p} (2R + Y_{new} + U_{new}) + \frac{8M}{n} (\frac{t_{s}}{m} + t_{w})}, n, m fixed,
lim_{N→∞} E_{p} = \frac{M(2R + Y_{new} + U_{new})}{M(2R + Y_{new} + U_{new}) + \frac{8Mp}{n}(\frac{t_{s}}{m} + t_{w})}, n, m fixed.

We can notice that speedup and efficiency increase slightly by increasing N, but the communication cost cannot be reduced relatively in such a way as by increasing n. (In steps S2 and S3 in the parallel algorithm some data exchanges are required in every loop run independent of N.) By using p = n processors we obtain increasing speedup but decreasing efficiency with increasing M or N (it follows from the formulas to determine S_{p} and E_{p}). Increasing M has only a very small influence on speedup and efficiency, because all algorithm steps except initialization step S1 are carried out M times.

6 Numerical Results

We discuss some numerical experiments for several examples of our optimal control problem (9) with

y_{d}(x, t) = sin(πx)e^{t}

f(x, t) ≡ 1
y_{0}(x) = sin(πx).

The algorithm was testet by using a two-dimensional mesh of transputers (T805). The ring topology, used by the parallel algorithm as communication topology, can easily be embedded into a real mesh topology. Some results obtained with p = 1 and p = 8 processors are shown in Table 1. We have obtained acceptable amounts of speedup and efficiency. We notice that the influence of the iteration number N is very small as discussed in the previous chapter. That means already for a small iteration number an acceptable effect of our parallelization was achieved.

N Time (sec) Speedup Efficiency
T_{1}  T_{p}  S_{p}  E_{p}
1  1,85 0,40  4,63  0,58
5  8,86 1,70  5,21  0,65
10 17,61 3,32  5,30  0,66
20 35,12 6,57  5,35  0,67
50 87,65 16,31  5,37  0,67
100 175,21 32,55  5,38  0,67

Table 1: Experimental results (p = 8)

A very important fact to get an efficient implementation of a parallel algorithm is to minimize the number of necessary data exchanges, i. e. the number of sending routine calls. This can be achieved by sending rather few long messages than many short messages. The reason is the high startup time compared with the lower word-transfer time in message passing systems. Thus in our parallel algorithm complete vectors are sent. Some tests showed that sending of separate values (instead of complete vectors) decreases efficiency down to 0, 2. 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 those with shared memory.

References

[BGR92] A. Bensoussan, R. Glowinski, A. Rascanu: Approximation of some stochastic differential equations by the splitting up method, Appl. Math. Opt. 25, No. 1, 81-106, 1992

[CP78] R. F. Curtain, A. J. Pritchard: Infinite dimensional linear systems theory, Lect. Notes in Contr. Inf. Sci. 8, 1978 [DG95] A. M. Davie, J. G. Gaines: Convergence of Implicit Schemes for Numerical Solution

of Parabolic Stochastic Partial Differential Equations, Preprint of the University of Edinburgh, Edinburgh, 1995.

[E76] M. Ermolev: Methods of stochastic programming (in Russ.), Glawnaja Redakzija Fisiko-Matematičskoi Literatiri, Moskva, 1976.

[G87] W. Grecksch: Stochastische Evolutionsgleichungen und deren Steuerung, Teubner-Texte zur Mathematik 98, Teubner, Leipzig, 1987

[GT95] W. Grecksch, C. Tudor: Stochastic Evolution Equations: A Hilbert Space Approach, Akademie Verlag, Berlin, 1995.

[J82] G. Jetschke: Different approaches to stochastic parabolic differential equations, Proceedings of the 10th Winter School, Supplemento ai Rendiconti del Circolo Matematico di Palermo, Serie II, No. 2, 161-169, 1982.

[J91] G. Jetschke: Lattice approximation of a nonlinear stochastic partial differential equation with white noise, International Series of Numerical Mathematics, Vol. 102, 107-126, Birkhäuser Verlag, 1991.

[KD92] H. J. Kushner, P. G. Dupuis: Numerical Methods for Stochastic Control Problems in Continuous Time, Springer, New York, 1992.

[KR79] N. V. Krylov, B. L. Rozovskij: Ob evoljucionnych stochastičeskich uravnenijach V, R. V. Gamkrelidze (Glavn. Red.): Itogi nauki i techniki, serija sovremennych problemy matematiki, tom 14, VINITI, Moskva, 71-146, 1979.

[KGGK94] V. Kumar, A. Grama, A. Gupta, G. Karypis: Introduction to Parallel Computing. Design and Analysis of Algorithms, Benjamin/Cummings Publishing Comp., Redwood City, 1994

[LP66] E. S. Levtin, B. T. Poljak: Metody Minimizacija c (in Russ.), Journal Vyčisli. Matematiki i Matematičeskoj Fiziki, 787-823, 1966.

[P71] L. L. Ponomarento: Stochastic integrals with multiparameter Brownian motion and corresponding stochastic equations, Teor. ver. i mat. stat., No. 2, 100-109, 1971.

[PZ92] G. Da Prato, J. Zabczyk: Stochastic Equations in Infinite Dimensions, Cambridge University Press, 1992.

[PGB99] T. Pohl, W. Grecksch, H. Blaar: A Parallel Application of a Quasigradient Method in Stochastic Control Theory, to appear in Optimization, 1999.

[Q94] M. J. Quinn: Parallel Computing. Theory and Practice, 2nd ed., McGraw-Hill, New York, 1994