IEEE TRANSACTIONS ON IMAGE PROCESSING,VOL.10,NO.4,APRIL 2001 493

Fast Image Recovery Using Dynamic Load

Balancing in Parallel Architectures,by Means of

Incomplete Projections

Francisco J.González-Castaño,Ubaldo M.García-Palomares,José L.Alba-Castro,Associate Member,IEEE,

and José M.Pousada-Carballo

Abstract This paper formulates an incomplete projection algo-

rithmthat is applied to the image recovery problem.The algorithm

allows an easy implementation of dynamic load balancing for par-

allel architectures.Furthermore,the local computation - commu-

nication load ratio can be adjusted,since each processor performs

a finite number of iterations of any projection-type technique,and

this number can be provided as a parameter of the algorithm.Nu-

merical results compare favorably with those obtained by the ex-

trapolated method of parallel subgradient projections.

Index Terms Load balancing,parallel algorithms,recovery,

restoration.

I.I

NTRODUCTION

T

HE image recovery problemis the estimation of an image

that has been degraded (blurred) by a process we have

partial information about [1].In convex set-theoretic image re-

covery [2] the solution space is a Hilbert space,but since we are

interested in strong convergence we restrict our solution space

to a finite dimensional space.The image is described solely by

a family of convex constraints

arising

from observed data.If we define the convex sets

then the recovery problem reduces to the convex inequality

problem (CIP) of finding

in the closed convex set

.We refer the reader to the bibliography

survey in [1] for related and recent work.Other interesting

references are [3],[4],[9],[10],and [17].

Sequential and parallel projection techniques are robust

methods for solving the CIP,and literature on their properties

and convergence abounds.A recent and comprehensive survey

is [5].Given an estimate

,the key idea is to simultaneously

compute

,the projections of

on closed convex sets

and generate the next estimate

as

(1)

Manuscript received February 22,1999;revised January 12,2001.The as-

sociate editor coordinating the review of this manuscript and approving it for

publication was Prof.Timothy J.Schulz.

F.J.González-Castaño,J.L.Alba-Castro,and J.M.Pousada-Carballo are

with the Departamento de Tecnologías de las Comunicaciones,E.T.S.I.Teleco-

municación,Campus Universitario,36200 Vigo,Spain.

U.M.García-Palomares is with Departamento de Procesos y Sistemas,Uni-

versidad Simón Bolívar,Caracas 1080-A,Venezuela.

Publisher Item Identifier S 1057-7149(01)02729-4.

The notation in (1) will be detailed later.In the sequel,we

commit a minor abuse in notation and denote

as the

sequence generated by an infinite loop of a programming lan-

guage.

It is well known that,under suitable conditions,the sequence

,generated by the basic algorithmgivenbelow,converges

to some

[5].

Basic algorithm

Given an initial estimate

,

REPEAT Choose the set

,the relaxation factor

and the

weights

,and let

UNTIL Convergence

END.

Convergence results for the previous algorithm exist for

and

[5].Recent research shows that

convergence is preserved,and in fact improved,for relaxation

parameters

out of the conventional (0,2) interval [1],[9],

[12],[16].Within the basic algorithms framework,Combettes

suggested the extrapolated method of parallel subgradient

projections (EMOPSP),which outperforms previous strategies

like Successive Projection [6],SIRT [7] and Simultaneous

Projection [8].These strategies require exact projections on

,whereas EMOPSP carries out (easy) projections on affine

hyperplanes.This feature makes the performance of EMOPSP

look superior.

EMOPSP can be viewed as an acceleration of the scheme

presented by Bauschke and Borwein [5].It is worth mentioning

that

(2)

and a sufficient convergence condition is

(3)

Suitable conditions to ensure (3) are well known [1],[5].

Since our implementation slightly differs from Combettes,

we describe below one iteration of EMOPSP,for the special

case of equal

weights.Further details and generalization to

a Hilbert space can be found in [1].Given

and

,

10577149/01$10.00 © 2001 IEEE

494 IEEE TRANSACTIONS ON IMAGE PROCESSING,VOL.10,NO.4,APRIL 2001

EMOPSP ITERATION (Generate

)

1.Determine an index set

of violated constraints,i.e.,

.

2.Compute subgradients

,and define

.

3.Compute

and

define

.

4.Define

.

END OF ITERATION.

In a multiprocessor environment,Combettes suggests making

,the number of processors.Note that steps 2 and

3 can be carried out in parallel,whereas the other steps need

communication among processors.Unless the subgradient eval-

uation

involves a large computational load,communi-

cation costs could surpass the computational costs,which ren-

ders a parallel method inefficient.

Themainobjectiveof this paper is toproposeamethodthat im-

proves the local computation versus communication load ratio,

and hopefully improve the speedup.The next section formulates

the incomplete projection algorithm (IPA),which replaces the

exact projections

in (1) by points that satisfy

the Fejér condition (FC) stated below.Section III explains how

to apply IPA in order to 1) produce a complete load balance

among processors and 2) increase the processors computational

load.Numerical results inSection IVshowthat IPAis better than

EMOPSP as a fast initial approximation to the solution.

II.I

NCOMPLETE

P

ROJECTION

Inequality (2),usually coined as the Fejér property,holds for

all projection-type methods,and is endowed with nice conver-

gence properties [5].Given

we say that

in (1) sat-

isfies the Fejér condition (FC) if

We next prove a lemma that will enable us to develop the

incomplete projection algorithm (IPA).Its main feature is that

points satisfying (FC) substitute exact projections.This is the

key to implement algorithms with a better computation-com-

munication cost ratio.

Lemma 1:If

satisfies (FC),then

Proof:If there exists

such that

we obtain

a contradiction

Theorem 1:Let

be a nonfeasible point (

),

let

be points in

satisfying (FC),let

,and let

.

Then

Proof:By the previous lemma we have

,we can generate,either sequentially or in parallel,points

by incomplete projection,which is namely,

a finite number of projection steps (or any other technique

that ensures that (FC) is fulfilled).Then we form a linear

nonnegative combination of the directions

to

generate the next estimate

,which also satisfies the Fejér

condition (FC).

Note that

(4)

If we want to obtain optimal

values,we may solve either

the following quadratic programming problem:

Maximize

(5)

or the following fractional programming problem:

Maximize

(6)

which can be solved by a sequence of quadratic programming

problems.Problem (6) was also suggested by Kiwiel [12].We

do not proceed with this matter any further,because an effi-

cient implementation of a quadratic programming algorithm is

required.Here we report results with

.

Let us mention that if exact projections are used,i.e.,

,very close results are obtained [16].In particular,

GONZÁLEZ-CASTAÑO et al.:FAST IMAGE RECOVERY USING DYNAMIC LOAD BALANCING 495

(7)

III.D

YNAMIC

L

OAD

B

ALANCING IN

P

ARALLEL

A

RCHITECTURES BY

M

EANS OF

I

NCOMPLETE

P

ROJECTION

Parallel architectures and,in particular,clusters of work-

stations [13],[14] have been widely used in recent years.In

distributed memory parallel architectures,the scalability of a

method is especially sensitive to relative communication load.

We have stated that EMOPSP presents a large communication

load as compared to local computation load in a parallel

environment.This is mainly because each processor solely

completes one easy projection at step 3 of the EMOPSP

iteration.We intend to replace that step by a finite number

of iterations of a projection-type method,which,as we have

pointed out,generates a

that satisfies (FC).This keeps the

desirable effect of having the processors carry out homoge-

neous tasks.In order to do this effectively we need to provide a

dynamic load balance to the processors.If this is possible,the

effective execution time per task will be the same.

In general,the load balancing problemis a very difficult one.

We need to redistribute the data at the beginning of a parallel

iteration.Moreover,fast dynamic load balancing requires the

data to be easily available to all processors.In practice,in a

distributed memory architecture,this supposes that all the data

shouldbeknownandreadytobeusedbyall processors.This may

not be possible inlarge problems duetomemoryrestrictions,and

a static distribution must be imposed [15].If limited memory is

not an issue,we can use the easy dynamic load balance given in

step 1 of the dynamic IPAalgorithmbelow( DIPA).We believe

that the joint effect of the dynamic load balance and efficient

local projection methods,like those proposed in [1],[11],[18],

makes our algorithm amenable for solving the general convex

feasibility problem.In short,we do the following.

We impose some control over the local computation,

thus improving the computation versus communication

ratio.We generate

in a finite number of

iterations of a projection method.

We dedicate some communication load to the search of

a good balance load among processors,providing the

dynamic load balancing we are looking for.

We have implemented the incomplete projection algorithm

with dynamic load distribution (DIPA) given below,with ad-

ditional peculiarities that will be explained in the next section.

Remarks that follow the description of the algorithm attempt

to further clarify its implementation.We do not claim that this

is the best algorithm,because it is certainly problem-depen-

dent. We are assuming that we have

processors and

con-

straints.

Incomplete Projection Algorithm with Dynamic Load

Distribution (DIPA).

Given an estimate

of the original image and a tolerance

,

REPEAT

1) Dynamic load:

(a) Determine an index set of violated constraints:

.

(b) Distribute the set

evenly among processors:Define

such that

,and

.

2) EMOPSP iteration:

(a)EachprocessorperformsoneEMOPSPiterationwith

.

(b) Let

.

3) Combination:

(a) Define

.

UNTIL

Remark:Optimal

values were not considered.

Remark 1a:To ensure convergence (3) must hold,that is

Remark 1b:

is desirable.All proces-

sors would theoretically employ the same execution time if

they compute alike projections at step 2.Since we do not im-

pose

[1],[5],[12],and the usual conver-

gence test is

.In some cases,it is

advisable to state a less stringent condition,namely,

Max

,where

is the

direction generated by processor

at the

of the squared distance of

to the set

,namely,

,is known,the algorithm can

stop as soon as [12]

496 IEEE TRANSACTIONS ON IMAGE PROCESSING,VOL.10,NO.4,APRIL 2001

Fig.1.Original image.

This stopping condition is indeed a warning of infeasibility and

arises from (7) with

at all iterations.

IV.N

UMERICAL

R

ESULTS

This section applies IPA with dynamic load balancing to

a standard digital image restoration problem.We performed

numerical tests and compared them with EMOPSP,as imple-

mented in [1].All images are 8-bit gray-scale,and have

pixels,with

.Problemlayout is similar to [1].

Original image

(Fig.1) is degraded by convo-

lutional blur with an uniform

kernel.Edge effects

are taken into account.

Uniform white noise

,and

,for

,where

is the

th component of the degraded image,and

is the

row of the matrix

.We recall that

,

where

.Note that

,

but the projection of a point on

is a difficult task.On the

contrary,projections on

are straightforward [1].We

highlight

,the projection on

,because we will especially

consider it in the algorithm

(8)

where

.We used

will take significantly longer processing time

than

.Hence,we decided to avoid the computation

of

at each iteration and implemented the dynamic load bal-

ancing among

.The proposed algo-

rithm is as follows.

GONZÁLEZ-CASTAÑO et al.:FAST IMAGE RECOVERY USING DYNAMIC LOAD BALANCING 497

Adapted DIPA

Given the iteration counter

of the

original image,

REPEAT Set the iteration counter

2)

3) Generate

:

(a) Set

(b) Parallel stage.Dynamic load:

Determine the index set of violated constraints:

4) Combination:

(a) Define

.

(b) Let

.

UNTIL

or execution

time limit is reached.

END of Adapted DIPA

In order to verify the effectiveness of the dynamic load dis-

tribution,we also implemented a static load.At the onset,we

evenly split the sets

among all pro-

cessors,that is,step 3b is omitted.We must mention that in all

cases both the static and the dynamic algorithms reached the ex-

ecution time limit.

To compare results,we defined the comparison functions

and

(9)

,

and

498 IEEE TRANSACTIONS ON IMAGE PROCESSING,VOL.10,NO.4,APRIL 2001

Fig.5.Dynamic IPA after 4 min,eight processors.

V.C

ONCLUSIONS

We have formulated the (Dynamic) Incomplete Projection

Algorithm,which theoretically does not require exact projec-

tions on convex supersets.It allows a dynamic load distribu-

tion in parallel computation and an efficient local computation

to communication ratio.DIPA can exploit the same EMOPSP

strategies and it is extremely flexible as well.Numerical results

showthat DIPAcan obtain a fast approximation to a discretized

image.Furthermore,it admits roomfor improvements by means

of strategies like set selection,weight optimization,or appro-

priate choice of projection-type methods.DIPA can be indeed

considered as an enhancement of EMOPSP,or any other projec-

tion type technique,in parallel environments.

A

CKNOWLEDGMENT

The authors thank Prof.R.R.Rubio for helping them with

the programming interface to the ORIGIN 2000.They also ac-

knowledge the careful reading by anonymous referees,who sug-

gested meaningful improvements to this paper.

R

EFERENCES

[1] P.L.Combettes,Convex set theoretic image recovery by extrapolated

iterations of parallel subgradient projections, IEEE Trans.Image Pro-

cessing,vol.6,pp.493506,Apr.1997.

[2] P.L.Combettes,The foundations of set theoretic estimation, Proc.

IEEE,vol.81,pp.182208,Feb.1993.

[3] Y.Censor,M.D.Altschuler,and W.D.Powlis,On the use of Cim-

minos simultaneous projections method for computing a solution of the

inverse problemin radiation therapy treatment planning, Inv.Prob.,vol.

4,pp.607623,1988.

[4] G.T.Herman and L.B.Meyer,Algebraic reconstruction techniques

can be made computationally efficient, IEEE Trans.Med.Imag.,vol.

12,pp.600609,1993.

[5] H.H.Bauschke and J.M.Borwein,On projection algorithms for

solving convex feasibility problems, SIAM Rev.,vol.38,no.3,pp.

367426,1996.

[6] L.G.Gurin,B.T.Polyak,and E.V.Raik,The method of projections for

finding the common point of convex sets, USSR Comput.Math.Math.

Phys.,vol.7,no.6,pp.124,1967.

[7] P.Gilbert,Iterative methods for the three-dimensional reconstruction of

an object fromprojections, J.Theor.Biol.,vol.36,no.1,pp.105117,

1972.

[8] P.L.Combettes,Inconsistent signal feasibility problems:Least-squares

solutions in a product space, IEEE Trans.Signal Processing,vol.42,

pp.29552966,Nov.1994.

[9] H.Ozaktas,M.C.Pinar,M.Akgul,T.Kurc,and C.Aykanat,An algo-

rithm with long steps for the simultaneous block projections approach

for the linear feasibility problem, J.Discrete Math.Sci.Cryptogr.,to

be published.

[10]

,Restoration of space-variant global blurs caused by severe

camera movements and coordinate distortions, J.Optics,vol.29,pp.

303310,1998.

[11] U.M.García-Palomares,Parallel projected aggregation methods for

solving the convex feasibility problem, SIAM J.Optim.,vol.3,pp.

882900,1993.

[12] K.C.Kiwiel,Block-iterative surrogate projection methods for convex

feasibility problems, Linear Algebra Applicat.,vol.215,pp.225259,

1995.

[13] T.M.Warschko,J.M.Blum,and W.F.Tichy,The parastation project:

Using workstations as building blocks for parallel computing, in Proc.

Int.Conf.Parallel Distributed Processing,Techniques,Applications,

1996.

[14] D.Ridge,D.Becker,and P.Merkey,Beowulf:Harnessing the power

of parallelism in a pile-of-PCs, Proc.IEEE Aerospace,1997.

[15] F.J.González-Castaño,Contribución al encaminamiento óptimo en

redes de circuitos virtuales con supervivencia,mediante arquitecturas de

memoria distribuida, Ph.D.dissertation,Universidad de Vigo,Spain,

1998.

[16] U.M.García-Palomares,Relaxation in projection methods, in

Encyclopedia of Optimization,P.M.Pardalos and C.A.Floudas,

Eds.Norwell,MA:Kluwer.

[17] Y.Censor,T.Elfving,and G.T.Herman,Linear algebra in image recon-

struction fromprojections, Linear Algebra Applicat.,vol.130,1990.

[18] U.M.García-Palomares,A superlinearly convergent projection algo-

rithm for solving the convex inequality problem, Oper.Res.Lett.,vol.

22,pp.97103,1998.

[19] M.I.Sezan and H.Stark,Image restoration by the method of convex

projections:Part 2Applications and numerical results, IEEE Trans.

Med.Imag.,vol.MI-1,pp.95101,1982.

Francisco J.González-Castaño received the M.Sc.

and Ph.D.degrees (cum laude) in telecommunica-

tions engineering from the Universidad de Santiago,

Spain,in 1990 and the Universidad de Vigo,Spain,

in 1998,respectively.

He has written ten papers in international refereed

journals in diverse fields of computer science

and telecommunications and has participated in

several R&D projects.He is currently a Profesor

Titular with the Departamento de Tecnologías de

las Comunicaciones,Universidad de Vigo,and has

been a Visiting Assistant Professor with the Computer Science Department,

University of Wisconsin,Madison.He has one patent.

Ubaldo M.García-Palomares was born in Caracas,

Venezuela.He received the B.S.degree in electrical

engineering in 1966 from Universidad Central de

Venezuela,Caracas,and the Ph.D.degree from

University of Wisconsin,Madison,in 1973.

He retired as Full Professor from the Universidad

Simón Bolívar,Venezuala,in 1993,but is still

engaged in research and teaching at that university

and at Universidad de Vigo,Spain.His main

interests include the theory and applications of

convex feasibility and optimization problems.He has

written over 40 reviewed papers on several nonlinear programming subjects.

GONZÁLEZ-CASTAÑO et al.:FAST IMAGE RECOVERY USING DYNAMIC LOAD BALANCING 499

José L.Alba-Castro (A96) received the M.Sc.and

Ph.D.degrees (cumlaude) intelecommunications en-

gineering from the Universidad de Santiago,Spain,

in 1990 and the Universidad de Vigo,Spain,in 1997,

respectively.

His research interests include neural networks

for classification applications,image segmentation,

statistical pattern recognition,automatic speech,

and speaker recognition and biometrics.He is a

Professor of discrete signal processing and image

processing at the Universidad de Vigo.

José M.Pousada-Carballo received the M.Sc.and

Ph.D.degrees (cum laude) in telecommunications

engineering from the Universidad Politécnica de

Madrid,Spain,in 1991 and the Universidad de Vigo,

Spain,in 1999,respectively.

He has written eight papers in international ref-

ereed journals in diverse fields of computer science

and telecommunications,and has participated in sev-

eral R&D projects.He is currently a Profesor Titular

with the Departamento de Tecnologías de las Comu-

nicaciones,Universidad de Vigo.He has one patent.

## Σχόλια 0

Συνδεθείτε για να κοινοποιήσετε σχόλιο