Simulia Isight Optimization Medical

Superiorization: An optimization heuristic for medical physics

6

Instituto de Investigaciones en Matemáticas Aplicadas y en Sistemas,

8

Cd. Universitaria, C.P. 04510, Mexico City, Mexico

13

Purpose: To describe and mathematically validate the superiorization methodol-

14

ogy, which is a recently-developed heuristic approach to optimization, and to discuss

15

its applicability to medical physics problem formulations that specify the desired

16

solution (of physically given or otherwise obtained constraints) by an optimization

17

criterion.

18

Methods: The superiorization methodology is presented as a heuristic solver for

19

a large class of constrained optimization problems. The constraints come from the

20

desire to produce a solution that is constraints-compatible, in the sense of meeting

21

requirements provided by physically or otherwise obtained constraints. The underly-

22

ing idea is that many iterative algorithms for finding such a solution are perturbation

23

resilient in the sense that, even if certain kinds of changes are made at the end of

24

each iterative step, the algorithm still produces a constraints-compatible solution.

25

This property is exploited by using permitted changes to steer the algorithm to a

26

solution that is not only constraints-compatible, but is also desirable according to

27

a specified optimization criterion. The approach is very general, it is applicable to

2

many iterative procedures and optimization criteria used in medical physics.

29

Results: The main practical contribution is a procedure for automatically pro-

30

ducing from any given iterative algorithm its superiorized version, which will supply

31

solutions that are superior according to a given optimization criterion. It is shown

32

that if the original iterative algorithm satisfies certain mathematical conditions, then

33

the output of its superiorized version is guaranteed to be as constraints-compatible as

34

the output of the original algorithm, but it is superior to the latter according to the

35

optimization criterion. This intuitive description is made precise in the paper and

36

the stated claims are rigorously proved. Superiorization is illustrated on simulated

37

computerized tomography data of a head cross-section and, in spite of its general-

38

ity, superiorization is shown to be competitive to an optimization algorithm that is

39

specifically designed to minimize total variation.

40

Conclusions: The range of applicability of superiorization to constrained opti-

41

mization problems is very large.

42

producing a superiorization algorithm from an algorithm aimed at only constraints-

43

compatibility; while non-heuristic (exact) approaches need to be redesigned for a new

44

optimization criterion. Thus superiorization provides a quick route to algorithms for

45

the practical solution of constrained optimization problems.

46

Keywords: superiorization, constrained optimization, heuristic optimization, tomography,

49

Optimization is a tool that is used in many areas of Medical Physics. Prime examples are

50

radiation therapy treatment planning and tomographic reconstruction, but there are others

51

such as image registration. Some well-cited classical publications on the topic are1–12 and

52

some recent articles are13–26.

53

In a typical medical physics application, one uses constrained optimization, where the

54

constraints come from the desire to produce a solution that is constraints-compatible, in

55

the sense of meeting the requirements provided by physically or otherwise obtained con-

56

straints. In radiation therapy treatment planning, the requirements are usually in the form

3

of constraints prescribed by the treatment planner on the doses to be delivered at specific

58

locations in the body. These doses in turn depend on information provided by an imaging

59

instrument, typically a Magnetic Resonance Imaging (MRI) or a Computerized Tomogra-

60

phy (CT) scanner. In tomography, the constraints come from the detector readings of the

61

instrument. In such applications, it is typically the case that a large number of solutions

62

would be considered good enough from the point of view of being constraints-compatible;

63

to a large extent, but not entirely, due to the fact that there is uncertainty as to the exact

64

nature of the constraints (for example, due to noise in the data collection). In such a case,

65

an optimization criterion is introduced that helps us to distinguish the “better” constraints-

66

compatible solutions (for example, this criterion could be the total dose to be delivered to

67

the body, which may vary quite a bit between radiation therapy treatment plans that are

68

compatible with the constraints on the doses delivered to individual locations).

69

The superiorization methodology (see, for example,22,27–32) is a recently-developed heuris-

70

tic approach to optimization. The word heuristic is used here in the sense that the process

71

is not guaranteed to lead to an optimum according to the given criterion; approaches aimed

72

at processes that are guaranteed in that sense are usually referred to as exact. Heuristic

73

approaches have been found useful in practical applications of optimization, mainly because

74

they are often computationally much less expensive than their exact counterparts, but nev-

75

ertheless provide solutions that are appropriate for the application at hand33–35.

76

The underlying idea of the superiorization approach is the following. In many applica-

77

tions there exists a computationally-efficient iterative algorithm that produces a constraints-

78

compatible solution for the given constraints. (An example of this for radiation therapy

79

treatment planning is reported in36, its clinical use is discussed in15.) Furthermore, often

80

the algorithm is perturbation resilient in the sense that, even if certain kinds of changes are

81

made at the end of each iterative step, the algorithm still produces a constraints-compatible

82

solution27–30. This property is exploited in the superiorization approach by using such per-

83

turbations to steer the algorithm to a solution that is not only constraints-compatible, but is

84

also desirable according to a specified optimization criterion. The approach is very general,

85

it is applicable to many iterative procedures and optimization criteria.

86

The current paper presents a major advance in the practice and theory of superiorization.

87

The previous publications22,27–32 used the intuitive idea to present some superiorization

88

algorithms, in this paper the reader will find a totally automatic procedure that turns an

4

iterative algorithm into its superiorized version. This version will produce an output that

90

is as constraints-compatible as the output of the original algorithm, but it is superior to

91

that according to an optimization criterion.

92

true for a very large class of iterative algorithms and for optimization criteria in general,

93

typical restrictions (such as convexity) on the optimization criterion are not essential for

94

the material presented below. In order to make precise and validate this broad claim, we

95

present here a new theoretical framework. The framework of29 is a precursor of what we

96

present here, but it is a restricted one, since it assumes that the constraints can be all

97

satisfied simultaneously, which is often false in medical physics applications. There is no

98

such restriction in the presentation below.

99

The idea of designing algorithms that use interlacing steps of two different kinds (in our

100

case, one kind of steps aim at constraints-compatibility and the other kind of steps aim at

101

improvement of the optimization criterion) is well-established and, in fact, is made use of

102

in many approaches that have been proposed with exact constrained optimization in mind;

103

see, for example, the works of Helou Neto and De Pierro37,38, of Nurminski39, of Combettes

104

and coworkers40,41, of Sidky and Pan and coworkers23,42,43 and of Defrise and coworkers44.

105

However, none of these approaches can do what can be done by the superiorization approach

106

as presented below, namely the automatic production of a heuristic constrained optimization

107

algorithm from an iterative algorithm for constraints-compatibility. For example, in37 it is

108

assumed (just as in the theory presented in our29) that all the constraints can be satisfied

109

simultaneously.

110

A major motivator for the additional theory presented in the current paper is to get rid

111

of this assumption, which is not reasonable when handling real problems of medical physics.

112

Motivated by similar considerations, Helou Neto and De Pierro38 present an alternative

113

approach that does not require this unreasonable assumption. However, in order to solve

114

such a problem, they end up with iterative algorithms of a particular form rather than having

115

the generality of being able to turn any constraints-compatibility seeking algorithm into a

116

superiorized one capable of handling constrained optimization. Also, the assumptions they

117

have to make in order to prove their convergence result (their Theorem 15) indicate that

118

their approach is applicable to a smaller class of constrained optimization problems than

119

the superiorization approach whose applicability seems to be more general. However, for

120

the mathematical purist, we point out that they present an exact constrained optimization

5

algorithm, while superiorization is a heuristic approach. Whether this is relevant to medical

122

physics practice is not clear: exact algorithms are not run forever, but are stopped according

123

to some stopping-rule, the relevant questions in comparing two algorithms are the quality

124

of the actual output and the computation time needed to obtain it.

125

Ultimately, the quality of the outputs should be evaluated by some figures of merit

126

relevant to the medical task at hand.

127

involves superiorization is in30 (Section 4.3), which reports on comparing in CT the efficacy

128

of constrained optimization reconstruction algorithms for the detection of low-contrast brain

129

tumors by using the method of statistical hypothesis testing (which provides a P-value that

130

indicates the significance by which we can reject the null hypothesis that the two algorithms

131

are equally efficacious in favor of the alternative that one is preferable). Such studies bundle

132

together two things: (i) the formulation of the constrained optimization task and (ii) the

133

performance of the algorithm in performing that task. The first of these requires a translation

134

of the medical aim into a mathematical model, it is important that this model should be

135

appropriately chosen.

136

The superiorization approach is not about choosing this model, it kicks in once the model

137

is chosen and aims at producing an output that is “good” according to the mathematical

138

specifications of the constraints and of the optimization criterion. Thus superiorization has

139

been used to compare the effects on the quality of the output in CT when the optimization

140

criterion is specified by total variation (TV) versus by entropy28 or versus by the ℓ1-norm

141

of the Haar transform32. However, the current paper is not about discussing how to trans-

142

late the underlying medical physics task into a constrained optimization problem. For our

143

purposes here, we are assuming that the mathematical model has been worked out and

144

concentrate on the algorithmic approach for solving the resulting constrained optimization

145

problem. We claim that the evaluation of such algorithms should not be based on the

146

medical figures of merit mentioned at the beginning of the previous paragraph, but rather

147

on their performance in solving the mathematical problem. If “good” solutions to the con-

148

strained optimization problem are not medically efficacious, that indicates that something is

149

wrong with the mathematical model and not that something is wrong with the algorithmic

150

approach. For this reason, in this paper we will not carry out a careful investigation of the

151

medical efficacy of any algorithm in the manner that we have done in30 (Section 4.3), but will

152

restrict ourselves to a simple illustration of the performance of the superiorization approach

6

as compared to the previously published algorithm of42 that is aimed at performing exact

154

minimization.

155

Examples of such studies already exist. Superiorization was compared in27 with Algorithm

156

6 of40 and in45 with the algorithm of Goldstein and Osher that they refer to as TwIST46 with

157

split Bregman47 as the substep. In both cases the implementation was done by the proposers

158

of the algorithms. In these reported instances superiorization did well: the constraints-

159

compatibility and the value of the function to be minimized were very similar for the outputs

160

produced by the algorithms being compared, but the superiorization algorithm produced its

161

output four times faster than the alternative. It would be unjustified to draw any general

162

conclusions on the mathematical performance and speed of superiorization based on just a

163

few experiments, but the reported results are encouraging.

164

However, the main reason why we advocate superiorization is different from what is

165

discussed above. The reason why we claim it to be helpful in medical physics research is

166

that it has the potential of saving a lot of time and effort for the researcher. Let us consider

167

a historical example. Likelihood optimization using the iterative process of expectation

168

maximization (EM)48 gained immediate and wide acceptance in the emission tomography

169

community. It was observed that irregular high amplitude patterns occurred in the image

170

with a large number of iterations, but it was not until five years later that this problem

171

was corrected49 by the use of a maximum a posteriority probability (MAP) algorithm with

172

a multivariate Gaussian prior. Had we had at our disposal the superiorization approach,

173

then the introduction of an optimization criterion (Gaussian or other) into the iterative

174

expectation maximization (EM) process would have been a simple matter and we would

175

have saved the time and effort spent on designing a special purpose algorithm for the MAP

176

formulation. A TV -superiorization of the EM algorithm is presented in50.

177

Even though our major claim for superiorization is that it provides a quick route to

178

algorithms for the practical solution of constrained optimization problems, before leaving

179

this introduction let us bring up a question that has to do with the performance of the

180

resulting algorithms: Will superiorization produce superior results to those produced by

181

contemporary MAP methods or is it faster than the better of such methods? At this stage

182

we have not yet developed the mathematical notation to discuss this question in a rigorous

183

manner, we return to it in Subsection II F.

184

In the next section we present in detail the superiorization methodology. In the subse-

7

quent section we provide an illustrative example by reporting on reconstructions produced

186

by algorithms applied to simulated computerized tomography data of a head cross-section.

187

In the final section we discuss our results and present our conclusions.

A.

Problem sets, proximity functions and ε-compatibility

190

Although optimization is often studied in a more general context (such as in Hilbert or

191

Banach spaces), in medical physics we usually deal with a special case, where optimization

192

is performed in a Euclidean space RJ (the space of J-dimensional vectors of real numbers,

193

where J is a positive integer). As often appropriate in practice, we further restrict the

194

domain of optimization to a nonempty subset Ωof RJ (such as the nonnegative orthant RJ

195

that consists of vectors all of whose components are nonnegative).

196

We now turn to formalizing the notion of being compatible with given constraints, a

197

notion that we have used informally in the previous section. In any application, there is a

198

problem set T; each problem T ∈T is essentially a description of the constraints in that

199

particular case. For example, for a tomographic scanner, the problem of reconstruction for

200

a particular patient at a particular time is determined by the measurements taken by the

201

scanner for that patient at that time. The intuitive notion of constraints-compatibility is

202

formalized by the use of a proximity function Pr on T such that, for every T ∈T, PrT

203

maps Ωinto R+, the set of nonnegative real numbers; i.e., PrT : Ω→R+. Intuitively we

204

think of PrT (x) as an indicator of how incompatible x is with the constraints of T. For

205

example, in tomography, PrT (x) should indicate by how much a proposed reconstruction

206

that is described by an x in Ωviolates the constraints of the problem T that are provided

207

by the measurements taken by the scanner. For example, if we use b to denote the vector

208

of estimated line integrals based on the measurements obtained by the scanner and by A

209

the system matrix of the scanner, then a possible choice for the proximity function is the

210

norm-distance ∥b −Ax∥, which we will use as an example in the discussions that follow.

211

An alternative legitimate choice for the proximity function is the Kullback-Leibler distance

212

KL(b, Ax), which is the negative log-likelihood of a statistical model in tomography. The

213

special case PrT (x) = 0 is interpreted by saying that x is perfectly compatible with the

8

constraints; due to the presence of noise in practical applications, it is quite conceivable

215

that there is no x that is perfectly compatible with the constraints, and we accept an x

216

as constraints-compatible as long as the value of PrT (x) is considered to be small enough

217

to justify that decision. Combining these two concepts leads to the notion of a problem

218

structure, which is a pair ⟨T, Pr⟩, where T is a nonempty problem set and Pr is a proximity

219

function on T. For a problem structure ⟨T, Pr⟩, a problem T ∈T, a nonnegative ε and an

220

x ∈Ω, we say that x is ε-compatible with T provided that PrT (x) ≤ε.

221

As an example (whose applicability to tomographic reconstruction is illustrated in Section

222

III), consider the problem structure that arises from the desire to find nonnegative solutions

223

of sequences of blocks of linear equations. Then the appropriate choices are Ω= RJ

224

the problem structure is ⟨S, Res⟩, where the problem set S is

(1)

and the proximity function Res on S is defined, for any problem S

I=1

(bi −⟨ai, x⟩)2.

(2)

Note that each element of this problem set S specifies an ordered sequence of W blocks

229

of linear equations of the form ⟨ai, x⟩= bi where ⟨∗, ∗⟩denotes the inner product in RJ

230

(and thus S is an appropriate representation of the so-called “ordered subsets” approach to

231

tomographic reconstruction51, as well as of other earlier-published block-iterative methods

232

that proposed essentially the same idea52–54). The proximity function Res on S is the residual

233

that we get when a particular x is substituted into all the equations of a particular problem

236

We now define the concept of an algorithm in the general context of problem structures.

237

For technical reasons that will become clear as we proceed with our development, we intro-

238

duce an additional set ∆, such that Ω⊆∆⊆RJ. (Both Ωand ∆are assumed to be known

239

and fixed for any particular problem structure ⟨T, Pr⟩.) An algorithm P for a problem

240

structure ⟨T, Pr⟩assigns to each problem T ∈T an operator PT : ∆→Ω. This definition is

241

used to define iterative processes that, for any initial point x ∈Ω, produce the (potentially)

∞

k=0 (that is, the sequence x, PTx, PT (PTx) , · · · ) of points in

243

Ω. We discuss below how such a potentially infinite process is terminated in practice.

Selecting Ω= Rj

+ and ∆= RJ for the problem structure ⟨S, Res⟩of the previous subsec-

(3)

where S is the problem specified above (2) and, for 1 ≤w ≤W, BSw : ∆→∆is defined by

(4)

where ∥a∥denotes the norm of the vector a in RJ, and Q : ∆→Ωis defined by

248

(Qx)j = max {0, xj} , for 1 ≤j ≤J.

(5)

Note that RS : ∆→Ω. This specific algorithm R is a typical example of the so-called block-

249

iterative methods mentioned above. Except for the presence of Q in (3), which enforces

250

nonnegativity of the components, it is identical to an algorithm used and illustrated in31.

251

With the Q absent from the definition of the algorithm, Ωhas to be the whole of RJ; the

252

practical consequence of the presence versus the absence of Q in the tomographic application

253

is illustrated in Subsection III D. We note also that special cases of the presented algorithm

254

include the classical reconstruction methods ART (if ℓw = 1, for 1 ≤w ≤W) and SIRT (if

255

W = 1); see, for example, Chapters 11 and 12 of55.

256

For a problem structure ⟨T, Pr⟩, a T ∈T, an ε ∈R+ and a sequence R =

257

of points in Ω, we use O (T, ε, R) to denote the x ∈Ωthat has the following properties:

258

PrT(x) ≤ε and there is a nonnegative integer K such that xK = x and, for all nonnegative

Xk

> ε. Clearly, if there is such an x, then it is unique. If there is no

10

such x, then we say that O (T, ε, R) is undefined, otherwise we say that it is defined. The

261

intuition behind this definition is the following: if we think of R as the (infinite) sequence of

262

points that is produced by an algorithm (intended for the problem T) without a termination

263

criterion, then O (T, ε, R) is the output produced by that algorithm when we add to it

264

instructions that make it terminate as soon as it reaches a point that is ε-compatible with

267

The notion of a bounded perturbations resilient algorithm P for a problem structure

268

⟨T, Pr⟩has been defined in a mathematically precise manner29. However, that definition

269

is not satisfactory from the point of view of applications in medical physics (or indeed in

270

any area involving noisy data), because it is useful only for problems T for which there is

271

a perfectly compatible solution (that is, an x such that PrT (x) = 0). We therefore extend

272

here that notion as follows. An algorithm P for a problem structure ⟨T, Pr⟩is said to be

276

have that O (T, ε′, R) is defined for every ε′ > ε and for every sequence R =

(6)

where βkvk are bounded perturbations, meaning that the sequence (βk)∞

280

in RJ is bounded and, for all k ≥0, xk + βkvk ∈∆.

281

In less formal terms, the second of these properties says that for a strongly perturbation

282

resilient algorithm we have that, for every problem and any nonnegative real number ε, if it

283

is the case that for all initial points from Ωthe infinite sequence produced by the algorithm

284

contains an ε-compatible point, then it will also be the case that all perturbed sequences

285

satisfying (6) contain an ε′-compatible point, for any ε′ > ε.

11

Having defined the notion of a strongly perturbation resilient algorithm, we next show

287

that this notion is of relevance to problems in medical physics. We illustrate the use of this

288

in tomography in the next section. We first need to introduce some mathematical concepts.

289

Given an algorithm P for a problem structure ⟨T, Pr⟩and a T ∈T, we say that

290

P is convergent for T if, for every x ∈Ω, there exists a unique y (x) ∈Ωsuch that,

291

limk→∞(PT)k x = y (x), meaning that for every positive real number δ, there exist a non-

(Pt)K X −Y (X)

≤δ, for all nonnegative integers k ≥K.

293

If, in addition, there exists a γ ∈R+ such that PrT (y (x)) ≤γ, for every x ∈Ω, then we

294

say that P is boundedly convergent for T.

295

A function f : Ω→R is uniformly continuous if, for every ε > 0 there exists a δ > 0,

296

such that, for all x, y ∈Ω, |f(x) −f(y)| ≤ε provided that ∥x −y∥≤δ. An example of a

297

uniformly continuous function is ResS of (2), for any S ∈S. This can be proved by observing

298

that the right-hand side of (2) can be rewritten in vector/matrix form as ∥b −Ax∥and

299

then selecting, for any given ε > 0, δ to be ε/ ∥A∥, where ∥A∥denotes the matrix norm of

301

An operator O : ∆→Ω, is nonexpansive if ∥Ox −Oy∥≤∥x −y∥, for all x, y ∈∆.

302

An example of a nonexpansive operator is the RS of (3). The proof of this is also simple.

303

It follows from discussions regarding similar claims in27 that the BSw : RJ →RJ of (4) is a

304

nonexpansive operator, for 1 ≤w ≤W, and that the operator Q of (5) is also nonexpansive.

305

Obviously, a sequential application of nonexpansive operators results in a nonexpansive

306

operator and thus RS is nonexpansive.

307

Now we state an important new result that gives sufficient conditions for strong perturba-

308

tion resilience: If P is an algorithm for a problem structure ⟨T, Pr⟩such that,

309

for all T ∈T, P is boundedly convergent for T, PrT : Ω→R is uniformly

310

continuous and PT : ∆→Ωis nonexpansive, then P is strongly perturbation

311

resilient. The importance of this result lies in the fact that the rather ordinary condition

312

of uniform continuity for the proximity function and the reasonable conditions of bounded

313

convergence and nonexpansiveness of the algorithmic operators guarantee that we end up

314

with a strongly perturbation resilient algorithm. The proof of this new result involves some

315

mathematical technicalities and is therefore presented in the Appendix as Theorem 1.

317

Now suppose, as is indeed the case for the constrained optimization problems discussed

318

in the previous section, that in addition to a problem structure ⟨T, Pr⟩we are also provided

319

with an optimization criterion, which is specified by a function φ : ∆→R, with the

320

convention that a point in ∆for which the value of φ is smaller is considered superior (from

321

the point of view of our application) to a point in ∆for which the value of φ is larger. In the

322

tomography context, any of the functions of x that are listed as a “secondary optimization

323

criterion” (an alternative name is a “regularizer”) in Section 6.4 of55 is an acceptable choice

324

for the optimization criterion φ. These include weighted norms, the negative of Shannon’s

325

entropy and total variation. It is the last of these that we discuss in detail in the illustrative

326

example below. The essential idea of the superiorization methodology presented in this paper

327

is to make use of the perturbations of (6) to transform a strongly perturbation resilient

328

algorithm that seeks a constraints-compatible solution into one whose outputs are equally

329

good from the point of view of constraints-compatibility, but are superior according to the

330

optimization criterion. We do this by producing from the algorithm another one, called its

331

superiorized version, by making sure not only that the βkvk are bounded perturbations, but

Xk

, for all k ≥0.

333

In order to ensure this we introduce a new concept (closely related to the concept of a

334

“descent direction” that is widely used in optimization). Given a function φ : ∆→R and a

335

point x ∈∆, we say that a vector d ∈RJ is nonascending for φ at x if ∥d∥≤1 and

There Is A Δ > 0 Such That For All Λ ∈[0, Δ] ,

(x + λd) ∈∆and φ (x + λd) ≤φ (x) .

(7)

Note that irrespective of the choices of φ and x, there is always at least one nonascending

337

vector d for φ at x, namely the zero-vector, all of whose components are zero. This is a useful

338

fact for proving results concerning the guaranteed behavior of our proposed procedures.

339

However, in order to steer our algorithms toward a point at which the value of φ is small,

340

we need to find a d such that φ (x + λd) < φ (x) rather than just φ (x + λd) ≤φ (x) as in

341

(7). In some earlier papers on superiorization27–31 it was assumed that ∆= RJ and that φ

342

is a convex function. This implied that, for any point x ∈∆, φ had a subgradient g ∈RJ at

343

the point x. It was suggested that if there is such a g with a positive norm, then d should

344

be chosen to be −g/ ∥g∥, otherwise d should be chosen to be the zero vector. However,

13

there are approaches (not involving subgradients) to selecting an appropriate d; an example

346

can be found in32 in which d is found without using subgradients for the case when φ is the

347

ℓ1-norm of the Haar transform. The method we used for selecting a nonascending vector in

348

the experiments reported in this paper is specified at the end of Subsection III A.

350

We now make precise the ingredients needed for transforming an algorithm into its su-

351

periorized version. Let Ωand ∆be the underlying sets for a problem structure ⟨T, Pr⟩

352

(Ω⊆∆⊆RJ, as discussed at the beginning of Subsection II B), P be an algorithm for

353

⟨T, Pr⟩and φ : ∆→R. The following description of the Superiorized Version of Algorithm

355

all k ≥0, (6) is satisfied. We show this to be true, for any algorithm P, after the description

356

of the Superiorized Version of Algorithm P. Furthermore, since the sequence RT is steered

357

by Superiorized Version of Algorithm P toward a reduced value of φ, there is an intuitive

358

expectation that the output of the superiorized version is likely to be superior (from the

359

point of view of the optimization criterion φ) to the output of the original unperturbed

360

algorithm. This last statement is not precise and so it cannot be proved in a mathematical

361

sense for an arbitrary algorithm P; however, that should not stop us from applying the

362

easy procedure given below for automatically producing the Superiorized Version of P and

363

experimentally checking whether it indeed provides us with outputs superior to those of the

364

original algorithm. The well-demonstrated nature of heuristic optimization approaches is

365

that they often work in practice even when their performance cannot be guaranteed to be

366

optimal33–35.

367

Nevertheless, we can push our theory further than the hope expressed in the last para-

368

graph, by considering superiorized versions of algorithms that satisfy some condition. In

369

this paper, the condition that we discuss is strong perturbation resilience. We show below

370

that if P is strongly perturbation resilient, then, for any problem T ∈T, a sequence RT

371

produced by its superiorized version has the following desirable property: For all ε ∈R+, if



is defined for every x ∈Ω, then O (T, ε′, RT) is also defined for every

373

ε′ > ε; in other words, the Superiorized Version of Algorithm P provides an ε′-compatible

374

output. As stated above, the advantage of the superiorized version is that its output is

14

likely to be superior to the output of the original unperturbed algorithm. We point out that

376

strong perturbation resilience is a sufficient, but not necessary, condition for guaranteeing

377

such desirable behavior of the superiorized version, finding additional sufficient conditions

378

and proving that algorithms that we wish to superiorize satisfy such conditions is part of

379

our ongoing research.

380

The superiorized version assumes that we have available a summable sequence (γℓ)∞

381

positive real numbers (for example, γℓ= aℓ, where 0 < a < 1) and it generates, simultane-

384

depends on a specified initial point ¯x ∈Ωand on a positive integer N. It makes use of a

385

logical variable called loop.

(Viii)

set vk,n to be a nonascending vector for φ at xk,n

406

Next we analyze the behavior of the Superiorized Version of Algorithm P.

407

The iteration number k is set to 0 in (i) and xk = x0 is set to its initial value ¯x in (ii). The

408

integer index ℓfor picking the next element from the sequence (γℓ)∞

409

line (iii), it is repeatedly increased by line (xi). The lines (v) - (xix) that follow the repeat

410

in (iv) perform a complete iterative step from xk to xk+1, infinite repetitions of such steps

Xk
∞

k=0. During one iterative step, there is one application of

412

the operator PT, in line (xviii), but there are N steering steps aimed at reducing the value of

413

φ; the latter are done by lines (v) - (xvii). These lines produce a sequence of points

Xk

.

415

We prove the truth of the last sentence by induction on the nonnegative integers. For

416

n = 0, we have by lines (v) and (vi) that xk,0 = xk. But xk ∈Ω, since it is either ¯x that is

417

assumed to be in Ωdue to lines (i) and (ii) or it is in the range Ωof PT due to lines (xviii)

418

and (xix). Now we assume, for any 0 ≤n < N, that xk,n ∈∆and φ

419

show that lines (viii) - (xvii) perform a computation that leads from xk,n to an xk,n+1 ∈∆

Xk

. To see this, observe that line (viii) sets vk,n to be a

421

nonascending vector for φ at xk,n, which implies that (7) is satisfied with x = xk,n and

422

d = vk,n. Line (ix) sets loop to true, and it remains true while searching for the desired

423

xk,n+1, by repeatedly executing the loop sequence that follows line (x). In this sequence,

424

line (xi) increases ℓby 1 and line (xii) sets βk,n to γℓ. Thus for the vector z defined by line

Xk,N

, provided that βk,n is not greater than the δ in (7). Since

(Γℓ)∞

ℓ=0 is a summable sequence of positive real numbers, there must be a positive integer L

427

such that γℓ≤δ, for all ℓ≥L. This implies that if we applied lines (xi) - (xiii) often enough,

428

we would reach a vector z that satisfies z ∈∆and φ (z) ≤φ

429

(xiv) is not satisfied when the process gets to it, then lines (xi) - (xiii) are again executed

430

and eventually we get a vector z for which the condition in line (xiv) is satisfied due to the

432

time xk,n+1 is set to z and so we obtain that xk,n+1 ∈∆and φ

Xk

, as desired.

433

Line (xvii) sets loop to false and so control is returned to line (vii). When this happens for

434

the Nth time, it will be the case that n = N and therefore line (xviii) is used to produce

435

xk+1 ∈Ωand the increasing of k by line (xix) allows us then to move on to the next iterative

436

step. Infinite repetition of such steps produces the sequence RT =

Xk
∞

k=0 of points in Ω.

438

ε′ > ε, the Superiorized Version of Algorithm P produces an ε′-compatible output. Since P

439

is assumed to be strongly perturbation resilient, this desired result follows if we can show

Vk
∞

k=0 of vectors in RJ such that (6) is satisfied for all k ≥0. In view of line

442

(xviii), this is achieved if we can define the βk and the vk so that xk,N = xk + βkvk. This

Βk

vk,n.

(9)

That these assignments result in xk,N = xk + βkvk follows from lines (v) - (xvii). From line

448

condition expressed in (6) is that, for all k ≥0, xk + βkvk ∈∆. This follows from the fact

449

that xk,N = xk + βkvk is assigned its value by line (xvi), but only if the condition expressed

450

in line (xiv) is satisfied.

451

In conclusion, we have shown that the superiorized version of a strongly perturbation

452

resilient algorithm produces outputs that are essentially as constraints-compatible as those

453

produced by the original version of the algorithm. However, due to the repeated steering of

454

the process by lines (vii) - (xvii) toward reducing the value of the optimization criterion φ,

455

we can expect that the output of the superiorized version will be superior (from the point

456

of view of φ) to the output of the original algorithm.

F.

Information on performance comparison with MAP methods

458

Using our notation, the constrained minimization formulation that we are considering is:

460

minimize φ(x), subject to PrT (x) ≤ε.

(10)

The aim of superiorization is not identical with the aim of constrained minimization in (10).

461

One difference is that ε is not “given” in the superiorization context. The superiorization

462

of an algorithm produces a sequence and, for any ε, the associated output of the algorithm

463

is considered to be the first x in the sequence for which PrT (x) ≤ε. The other difference

464

is that we do not claim that this output is a minimizer of φ among all points that satisfy

465

the constraint, but hope only that it is usually an x for which φ(x) is at the small end

466

of its range of values over the set of constraint-satisfying points. This latter difference is

467

generally shared by comparisons of a heuristic approach with an exact approach to solving

468

a constrained minimization problem.

469

The MAP (or regularized) formulation of a physical problem that leads to the constrained

470

minimization problem (10) is the unconstrained minimization problem of the form: Given

472

minimize [φ(x) + βPrT (x)] .

(11)

Formulations of both kinds (i.e, the ones of (10) and of (11)) are widely used for solving

473

medical physics problems and the question “Which of these two formulations leads to faster or

474

better solutions of the underlying physical problem?” is open. Examples of both formulations

475

with various choices for PrT and φ are listed in the beginning parts of the paper of Goldstein

476

and Osher47.

477

We now return to the question raised near the end of Section I: Will superiorization pro-

478

duce superior results to those produced by contemporary MAP methods or is it faster than

479

the better of such methods? As yet, there is very little information available regarding this

480

general question; in fact, we are aware of only one published study45. That study compared

481

a superiorization algorithm with the algorithm of Goldstein and Osher that they refer to

482

as TwIST46 with split Bregman47 as the substep, which is indeed a contemporary method

483

that uses the MAP formulation. (For example, see the discussion of the split Bregman

484

method in56.) The problem S to which the two algorithms were applied was one from the

485

tomographic problem set S defined in (1). ResS as defined in (2) was used as the proximity

18

function and total variation, TV as defined below in (12), was the choice for φ. It is reported

487

in45 that for the outputs of the two algorithms that were being compared, the values of ResS

488

and TV were very similar, but the superiorization algorithm produced its output four times

489

faster than the MAP method.

492

We use tomography to refer to the process of reconstructing a function over a Euclidean

493

space from estimated values of its integrals along lines (that are usually, but not necessarily,

494

straight). The particular reconstruction processes to which our discussion applies are the

495

series expansion methods, see Section 6.3 of55, in which it is assumed that the function to

496

be reconstructed can be approximated by a linear combination of a finite number (say J)

497

of basis functions and the reconstruction task becomes one of estimating the coefficients of

498

the basis functions in the expansion. Sometimes, prior knowledge about the nature of the

499

function to be reconstructed allows us to confine the sought-after vector x of coefficients to

500

a subset Ωof RJ (such as the nonnegative orthant RJ

501

which we integrate, ai ∈RJ to denote the vector whose jth component is the integral of the

502

jth basis function along the ith line, and bi to denote the measured integral of the function

503

to be reconstructed along the ith line. Under these circumstances the constraints come from

504

the desire that, for each of the lines, ⟨ai, x⟩should be close (in some sense) to bi.

505

To make this concrete, consider (1).

506

tomography by grouping the lines of integration into W blocks, with ℓw lines in the wth block.

507

Such groupings often (but not always) are done according to some geometrical condition on

508

the lines (for example, in case of straight lines, we may decide that all the lines that are

509

parallel to each other form one block). In this framework the proximity function Res defined

510

by (2) provides a reasonable measure of the incompatibility of a vector x with the constraints.

511

The algorithm R described by (3) - (5) is applicable to this concrete formulation.

512

There are many optimization criteria that have been used in tomography, see Section

513

6.4 of55, here we discuss the one called total variation (TV ), whose use has been popular

514

in medical physics recently, see as examples20,22,23,41–44. The definition of TV that we use

Related Journal Articles & DOI Links

Selected peer-reviewed publications relevant to 12 Lead ECG Acquisition. Click the DOI to access the full paper (may require institutional access).

Why Choose Us?

Bangalore guidance for robotics, Spectre and autonomous systems projects.