Parallel Optimization: Theory, Algorithms, and Applications,

Organization of the Book

• The book has three parts—Theory, Algorithms, and Applications—with Chapter 1 introducing parallel computing.Source: Organization of the Book · p. 12
• The theory part develops foundations for iterative projection, model decomposition, and interior point algorithms.These foundations support the algorithms developed in Chapters 5–8.Source: Organization of the Book · pp. 12–13
• The applications part models problems across matrix estimation, image reconstruction, radiation therapy, transportation, and planning under uncertainty.Source: Organization of the Book · p. 13
• Chapters 14 and 15 cover parallel implementation and numerical experiments, with implementation discussion aimed at classes of machines rather than one hardware model.Source: Organization of the Book · p. 13
• The chapters can serve as a sequential text or reference; readers can study applications selectively, while understanding convergence proofs requires the relevant theory.Source: Organization of the Book · p. 13

Suggested course outlines

• Course paths can emphasize constrained-optimization theory and algorithms, structured numerical methods, or applications, with different chapter selections for each.An additional parallel-computing course can use algorithms and applications for implementation exercises, then Chapter 15 to test performance.Source: xiv · pp. 15–16

• The contents list the Algebraic Reconstruction Technique for systems of equations.Source: xxii · p. 24
• The chapter includes iterative data refinement, covering its fundamentals and applications in medical imaging.Source: xxii · p. 24
• A separate section addresses selective use of iterative algorithms for inversion problems.Source: xxii · p. 24

• The chapter begins with problem definition and continuous forward and inverse models.Source: §11 · p. 24
• Its listed computational topics include discretizing the feasibility problem and computational inversion of the data, alongside consequences, limitations, and experimental results.Source: §11 · p. 24
• The outline treats treatment-plan combination separately for feasible and infeasible cases.Source: §11 · p. 24

12.2 Problem Formulations

• The formulation material covers transportation problems and multicommodity network flow problems.Source: §12 · p. 24

12.3 Sample Applications

• The sample applications include covering positions in stock options, air-traffic control, and traffic routing.Source: §12 · p. 24

12.4 Iterative Algorithms for Multicommodity Network Flow Problems

• The iterative-algorithm material includes a row-action algorithm for quadratic transportation problems and extensions to generalized networks.Source: §12 · p. 24

1.1.3 Control and data parallelism

• Data parallelism scales with available data, whereas control parallelism is limited by serial work; with infinitely many processors, speedup remains bounded by 1/s.For control parallelism, s is the serial fraction of the code.Source: §1.1.3 · p. 40, Data Parallelism · p. 41

Data Parallelism

• Data-parallel execution partitions problem data into tasks, with virtual-to-physical processor mapping guided by communication needs and the target network.A physical processor may emulate multiple virtual processors, and modern data-parallel languages can provide mapping primitives.Source: Data Parallelism · p. 41

1.2 How Does Parallelism Affect Computing?

• Algorithms should be flexible enough to run across a range of parallel architectures rather than being tied to one machine’s specifics.The book emphasizes examining iterative algorithms’ general structure to assess their parallel suitability and adaptability.Source: §1.2 · pp. 43–44

1.3 A Classification of Parallel Algorithms

• Parallelism may arise from algorithm structure or problem structure; iterative projection methods reorganize operations, while model decomposition and interior-point methods exploit model structure.The text identifies these as three broad algorithmic classes with distinct ways of exploiting parallelism.Source: §1.3 · p. 44

1.3.1 Parallelism due to algorithm structure: Iterative projection algorithms

• Iterative projection algorithms are classified by row- or block-action updates and by whether iterations proceed sequentially or simultaneously.Block-iterative schemes may use fixed blocks or changing blocks, subject to convergence restrictions.Source: §1.3.1 · pp. 45–46

Parallel Computing with Iterative Algorithms

• Sequential row-action updates can run concurrently when they change different vector components, while simultaneous updates require synchronization to combine intermediate results.Block choices can balance task difficulty across processors; fine-grain machines favor many small blocks, while coarse-grain machines favor fewer, larger blocks.Source: Parallel Computing with Iterative… · pp. 47–48

1.3.2 Parallelism due to problem structure: Model decomposition and interior point algorithms

• A block-separable objective is a sum of functions on variable subsets; the section introduces a constrained optimization model with objective sum_k f_k(x_k).The displayed formulation places the concatenated variable vector x in a constraint set D.Source: §1.3.2 · p. 48, K Optimize F(x) ~ L fk(xk) · p. 48

K Optimize F(x) ~ L fk(xk)

• Model decomposition removes complicating constraints in modified problems, yielding independent subproblems that can be solved in parallel and repeated until an acceptable original-problem solution is reached.For multicommodity flow, the block-angular constraints can be represented using per-commodity flow constraints and shared capacity constraints.Source: K Optimize F(x) ~ L fk(xk) · p. 49

Interior Point Algorithms

• Interior-point methods exploit block-angular structure by solving independent smaller linear-algebra subsystems in parallel, although the algorithms themselves are not inherently parallel.The main computational effort is factorizing and solving a system of equations.Source: Interior Point Algorithms · p. 51

1.4 Measuring the Performance of Parallel Algorithms

• Efficiency is relative speedup divided by processor count, while scalability measures maintaining constant efficiency as processors and problem size increase.Amdahl’s law bounds relative speedup by the serial fraction; performance measures are partial, since high speedup or FLOPS may not imply faster solution overall.Source: §1.4 · pp. 52–53

1.5 Notes and References

• The notes point readers to work on parallel optimization, block-separability, block-angular structure, and parallel-algorithm performance measures.They cite Kindervater, Lenstra, and Rinnooy Kan (1989), Zenios (1994b), and Barr and Hickman (1993), among others.Source: §1.5 · pp. 53, 55

Projections

• Generalized distances may be asymmetric and need not satisfy the triangle inequality; under suitable conditions, they define projections onto convex sets by minimizing distance.The chapter presents generalized distances and projections as tools for subsequent parallel optimization algorithms.Source: Projections · p. 58

2.1 Bregman Functions and Generalized Projections

• A Bregman generalized distance is the gap between a function and its tangent hyperplane at the reference point.The supporting passage interprets this gap geometrically for a convex function.Source: §2.1 · p. 59
• For every Bregman function, its generalized distance is nonnegative and equals zero exactly when its two arguments coincide.Source: §2.1 · pp. 60–61
• Every nonempty intersection of a closed convex set with the zone has a unique generalized projection of any point in the zone.Source: §2.1 · p. 61
• The negative of Shannon’s entropy is a Bregman function on the strictly positive orthant, with its associated distance given by the Kullback–Leibler cross-entropy function.Source: §2.1 · p. 62

2.2 Generalized Projections onto Hyperplanes

• For a zone-consistent Bregman function, a system of equations uniquely determines both the generalized projection onto a hyperplane and its parameter.Source: §2.2 · p. 64
• For points outside a hyperplane, the projection parameter has the sign of the hyperplane offset from the point; parameters for points in opposing half-spaces have inverse signs.Source: §2.2 · pp. 65–66
• Successive generalized projections onto two parallel hyperplanes yield the same final point as direct projection onto the second, with the successive parameters summing to the direct parameter.Source: §2.2 · p. 66

2.3 Bregman Functions on the Whole Space

• The section turns to whole-space Bregman functions because the six defining conditions can be difficult to verify in specific instances.Source: §2.3 · p. 67
• A C² function with positive definite Hessian and f(x) → +∞ as ||x|| → ∞ belongs to the whole-space Bregman class under sufficient, not necessary, conditions.The coercivity condition is not necessary: the text gives a twice continuously differentiable Bregman function for which f(x)/x → 0 as x → +∞.Source: §2.3 · pp. 68, 70

2.4 Characterization of Generalized Projections

• Under Theorem 2.4.1’s assumptions, a generalized projection is characterized by (u − x*, ∇f(x) − ∇f(x*)) ≤ 0 for every u in the feasible intersection.The theorem also gives the inequality D_f(P(y), y) ≤ D_f(z, y) − D_f(z, P(y)) for feasible z.Source: §2.4 · pp. 71–72

2.5 Csiszar ~.p-divergences

• Csiszar φ-divergences are nonnegative, vanish only when their arguments agree, have bounded level sets, and are jointly convex and strictly convex in the first argument.For fixed y*, dφ(y*, yν) → 0 if and only if yν → y*.Source: §2.5 · p. 74
• The divergence d'P is strictly convex in x and jointly convex in (x,y); its Kullback-Leibler case is the only known overlap with Bregman distances, up to linear terms and multiplicative constants.The Kullback-Leibler cross-entropy arises from φ1(t) = t log t - t + 1.Source: §2.5 · pp. 74–75

2.6 Notes and References

• The notes trace Bregman-function foundations to Bregman (1967b) and Censor and Lent (1981), and correct De Pierro and Iusem's strict-convexity assumption to positive definiteness.The notes also point to Csiszar (1991) and related axiomatic work on the families of generalized distances and divergences.Source: §2.6 · pp. 76–77

Proximal Minimization with D-Functions

• The proximal scheme replaces a quadratic regularizer with any Bregman distance and uses parallel algorithms on the resulting subproblems to solve problems that may lack parallel algorithms directly.The chapter establishes the mathematical foundation; applications to real-world problems are deferred to Chapter 15.Source: Proximal Minimization with D-Functions · p. 79

3.1 The Proximal Minimization Algorithm

• Section 3.1 introduces the quadratic proximal minimization method using a positive parameter sequence and an arbitrary initial vector y0.Source: §3.1 · p. 79, Algorithm 3.1.1 The Quadratic Proximal… · p. 80

Algorithm 3.1.1 The Quadratic Proximal Minimization Algorithm

• Algorithm 3.1.1 defines PMD using a Bregman function f and a positive sequence satisfying condition (3.5).Source: Algorithm 3.1.1 The Quadratic Proximal… · p. 80

Algorithm 3.1.2 The Proximal Minimization Algorithm with D-functions (PMD)

• The algorithm requires every iterate to remain in S because D_f is defined on S × S, restricting admissible choices of f and S.The text says this condition holds trivially when X ⊆ S and also holds for the entropy case.Source: Algorithm 3.1.2 The Proximal… · p. 80

3.2 Convergence Analysis of the PMD Algorithm

• For f ∈ B(S), c > 0, and F bounded below over X, the regularized objective has a unique minimizer over X ∩ S for every y ∈ S.The minimizer is denoted X_J(y,c).Source: §3.2 · p. 80
• If D_J is jointly convex, the regularized value function is convex; if f is also twice continuously differentiable, that function is continuously differentiable.The differentiability result also assumes joint convexity of D_f.Source: §3.2 · pp. 81–83
• When the Bregman function’s zone is R^n, the minimizers of the regularized value function coincide with the solution set X*.The section defines S* as the set of minimizers of that value function.Source: §3.2 · p. 83
• Under the stated assumptions, PMD iterates converge to an element of X* when X* ∩ S is nonempty and c_v > 0 with liminf c_v > 0.The assumptions include f ∈ B(S) and that every generated iterate remains in S.Source: §3.2 · p. 83

3.3 Special Cases: Quadratic and Entropic PMD

• Quadratic f recovers the original quadratic proximal algorithm, while negative Shannon entropy yields an entropy-type PMD whose practical advantage depends on problem structure and efficient step solvers.For linear programming, the entropy iteration becomes a linearly constrained entropy optimization problem, to which parallel Chapter 6 algorithms apply.Source: §3.3 · p. 86

3.4 Notes and References

• Censor and Zenios (1992) proposed PMD and established its convergence; later work extended it to maximal monotone operators and approximate iterations.The notes also cite alternative convergence proofs and related entropic proximal-map approaches.Source: §3.4 · p. 87
• Iusem and Teboulle (1993, 1995) studied convergence rates for PMD and EPMC, while Auslender and Haddou (1995) and Burachik (1995) extended them to variational inequalities.Other cited work addressed Bregman functions for variational inequalities, generalized Bregman functions, and relationships between PMD, central paths, and Cauchy trajectories.Source: §3.4 · p. 88

Augmented Lagrangians

• Augmented Lagrangian methods are closely related to proximal minimization and form part of a framework connecting penalty and barrier methods with primal-dual algorithms.The text also identifies augmented Lagrangians as suitable problem modifiers for model decomposition and barrier methods as useful for primal-dual path following.Source: Augmented Lagrangians · p. 89

4.1 Penalty Methods

• Penalty methods penalize constraint violations, and the constrained optimum is obtained as the limit of solutions to unconstrained penalized problems.A penalty function is continuous, nonnegative, and zero exactly on the feasible set; it is exact if a finite penalty parameter yields the original problem’s solution.Source: §4.1 · pp. 89–90

Algorithm 4.1.1 The Penalty Algorithmic Scheme

• Under the stated continuity, closedness, and growth or boundedness conditions, every accumulation point of the penalty iterates is optimal and p(x^v) tends to zero.The penalty parameters increase monotonically to infinity while each iteration solves an unconstrained penalized problem.Source: Algorithm 4.1.1 The Penalty… · pp. 90–91

4.2 Barrier Methods

• A barrier function is continuous and nonnegative in the interior of a robust feasible set and diverges toward every boundary point, keeping interior searches away from the boundary.The barrier method augments the objective with a positive multiple of the barrier and decreases that parameter toward zero.Source: §4.2 · pp. 92–93

Algorithm 4.2.1 The Barrier Algorithmic Scheme

• With a nonempty, closed, robust feasible set, continuous objective, and coercivity or boundedness, every accumulation point of barrier iterates is optimal and c_vq(x^v) tends to zero.The barrier parameters are positive and decrease monotonically to zero; the objective values also decrease monotonically across iterations.Source: Algorithm 4.2.1 The Barrier… · pp. 93–94

4.3 The Primal-Dual Algorithmic Scheme

• For equality-constrained problems, regularity and a positive-definite Lagrangian Hessian make the multiplier a local dual maximizer with the primal optimum as its value.The primal-dual approach uses this duality result to motivate alternating Lagrangian minimization and steepest-ascent updates in the dual variables.Source: §4.3 · p. 96

Algorithm 4.3.1 The Primal-Dual Algorithmic Scheme

• The primal-dual algorithm alternates Lagrangian minimization with steepest-ascent multiplier updates, but requires local convexity and can converge slowly; separability can make its minimizations efficient.Its positive, nondecreasing stepsizes govern the multiplier updates, while separable problems decompose the minimization into smaller subproblems.Source: Algorithm 4.3.1 The Primal-Dual… · pp. 96–97

4.4 Augmented lagrangian Methods

• Augmented Lagrangian methods combine penalty and primal-dual ideas to address their drawbacks; the described equivalent reformulation assumes convex F, a nonempty polyhedral feasible set, f(0)=0, and c>0.The text says these methods can perform better than penalty and classical primal-dual methods when properly employed under mild technical conditions.Source: §4.4 · pp. 97–98
• The augmented Lagrangian scheme alternates minimizing the augmented Lagrangian and ascending in the dual, with the multiplier update using a deflected-gradient direction and positive step parameters.For penalty functions in PE, the ascent nature of the update can be guaranteed.Source: §4.4 · p. 98, Algorithm 4.4.1 The Augmented… · p. 100

Algorithm 4.4.1 The Augmented Lagrangian Algorithmic Scheme

• For linear equality constraints, the augmented Lagrangian algorithm is equivalent to proximal minimization on the dual; with a quadratic penalty, it is precisely the quadratic proximal method, while general Bregman equivalence remains open.The text also identifies a similar equivalence between the exponential method of multipliers and entropic PMD.Source: Algorithm 4.4.1 The Augmented… · pp. 101, 103
• The augmented Lagrangian algorithm with linear equality constraints is closely related to proximal minimization, assuming each minimization is uniquely attained.For constant penalty parameters, the matrices B tend to the identity as the primal iterates approach a local solution, yielding a fixed-stepsize steepest-ascent iteration near a dual optimum when the stated derivative condition holds.Source: Algorithm 4.4.1 The Augmented… · p. 101

4.5 Notes and References

• Augmented Lagrangian methods originated with Hestenes (1969) and Powell (1969); the exponential multiplier method was proposed by Kort and Bertsekas (1972), with its entropic PMD equivalence analyzed by Tseng and Bertsekas (1993).Bertsekas (1982) provides an extensive treatment of nonquadratic penalties and references.Source: §4.5 · pp. 103–104

Problems

• The convex feasibility problem seeks a point in the nonempty intersection of finitely many closed convex sets, and the chapter surveys projection, subgradient, block-iterative, and simultaneous methods.It also describes replacing projections with operations such as separating-hyperplane projections and using nonsequential updates.Source: Problems · pp. 108–109

5.1 Preliminaries: Control Sequences and Relaxation Parameters

• Control sequences select sets or blocks for iterative steps; almost cyclic control requires every C consecutive indices to cover all sets and allows more flexibility than cyclic control, while repetitive control revisits every index.Remotest-set and approximately remotest-set controls are also described, with cyclic control included among the latter.Source: §5.1 · pp. 109–110
• Relaxation parameters can over- or undershoot an iterative step, adding implementation flexibility and affecting practical performance; convergence depends on technical conditions specific to each method.The text says their effects and convergence conditions are discussed in the sections that follow.Source: §5.1 · p. 110

Algorithm 5.2.1 The Method of Successive Orthogonal Projections (SOP)

• SOP successively projects from an arbitrary starting point onto individual convex sets, but generally requires solving a distance-minimization problem unless those projections are readily computed.Hyperplanes and half-spaces are examples where projections are easily calculated.Source: §5.2 · p. 110, Algorithm 5.2.1 The Method of… · p. 111, §5.4.1 · p. 116

5.4.1 The method of successive orthogonal projections

• CSP replaces SOP’s unknown projection-point gradient and auxiliary distance minimization with a gradient or subgradient at the known current iterate, though initial convergence may be slower.This trades potentially less efficient initial progress for avoiding distance-minimization problems.Source: §5.3 · p. 112, §5.4.1 · pp. 116–117

Algorithm 5.3.1 The Method of Cyclic Subgradient Projections (CSP)

• Under continuous convex constraint functions, a nonempty feasible intersection, and uniformly bounded subgradients near a solution, CSP iterates converge to a feasible solution.The theorem also requires an almost-cyclic control sequence and relaxation parameters bounded away from 0 and 2.Source: Algorithm 5.3.1 The Method of Cyclic… · pp. 113–114

5.4.2 A remotest set controlled subgradient projections method

• Cyclic control avoids identifying the remotest set, which would require projections onto all sets at each step and can be costly when a problem has many functions.The text identifies this flexibility as practically significant for real-world problems with many functions.Source: §5.4.2 · p. 117

Algorithm 5.4.1 Oettli's Scheme

• Oettli’s scheme applies subgradients of a monotonic norm to positive constraint violations; the l∞ norm yields remotest-set control, but no known norm yields cyclic control.The scheme’s control is determined by the choice of monotonic norm.Source: Algorithm 5.4.1 Oettli's Scheme · p. 118

Algorithm 5.4.2 The Relaxation Method of Agmon, Motzkin, and Schoenberg (AMS)

• For linear inequalities, cyclically controlled SOP and CSP coincide with the cyclic Agmon–Motzkin–Schoenberg relaxation method, whose relaxed update lies between the iterate and its reflection.With unity relaxation, the update is the orthogonal projection onto the violated half-space; a satisfied inequality leaves the iterate unchanged.Source: §5.4.4 · p. 118, Algorithm 5.4.2 The Relaxation Method… · pp. 118–119

5.4.5 Kaczmarz's algorithm for systems of linear equations and its nonlinear extension

• Kaczmarz's algorithm for linear equations can be derived by replacing each equality with opposite inequalities; it was independently discovered for image reconstruction as ART.Source: §5.4.5 · p. 119

Algorithm 5.4.3 Kaczmarz's Algorithm, Algebraic Reconstruction Technique (ART)

• For all v, relaxation parameters Av are confined to ε1 ≤ Av < 2 − ε2 for arbitrarily small positive ε1 and ε2, and the control sequence is almost cyclic on I.Source: Algorithm 5.4.3 Kaczmarz's Algorithm,… · p. 120

Algorithm 5.4.4 Nonlinear Extension of Kaczmarz's Algorithm

• The nonlinear Kaczmarz extension converges locally under conditions including second-order differentiability, without convexity assumptions, and incorporates no relaxation parameters.The stated version uses cyclic control and unity relaxation.Source: Algorithm 5.4.4 Nonlinear Extension of… · p. 121

5.5 The (6, 71)-Aigorithm

• For 0 < δ ≤ 1 and 0 < η ≤ 1, the (δ, η)-Algorithm converges to a point in the nonempty intersection of the convex sets.Its generated sequence is Fejér-monotone with respect to the solution set.Source: §5.5 · pp. 124–125
• Algorithm 5.5.2 requires an interior point in each set, avoids direct projections onto the sets, and converges to a solution of the convex feasibility problem.It constructs a separating supporting hyperplane and projects onto that hyperplane.Source: §5.5 · pp. 126–127
• The scheme is shown to be a (δ, γ)-algorithm with positive parameter d/T, using Fejér monotonicity and boundedness.Source: §5.5 · pp. 127–128
• For convex inequalities with cyclic control, the special scheme uses subgradients to construct supporting hyperplanes and a user-chosen relaxation parameter for its update.It begins with points strictly satisfying each inequality; the cited passage also establishes that the subgradient used is nonzero.Source: §5.5 · pp. 128–129

5.6 The Block-Iterative Projections Algorithm

• BIP encompasses sequential and fully simultaneous projections while allowing block count, size, and set assignments to change dynamically, subject to a technical weight condition.Source: §5.6 · p. 129

Algorithm 5.6.1 The Block-Iterative Projections Algorithm (BIP)

• Fair weight sequences and relaxation parameters define BIP; one-set weights recover SOP, weights on all sets give Cimmino-type updates, and fixed blocks yield an AMS relaxation variant.Source: Algorithm 5.6.1 The Block-Iterative… · pp. 129–130

Algorithm 5.6.2 The Block-AMS Algorithm

• Block-AMS permits changing block sizes and constraint assignments; Figure 5.9 illustrates a unity-relaxation step for three inequalities.The relaxation parameter is confined to the interval stated for the algorithm.Source: Algorithm 5.6.2 The Block-AMS Algorithm · p. 130, Figure 5.9 · p. 131

5.6.1 Convergence of the BIP algorithm

• With nonempty Q and fair weights, BIP converges to a point in the intersection of sets whose weights sum to infinity, with relaxation parameters bounded away from 0 and 2.For fixed weights, the text says this intersection is Q.Source: §5.6.1 · p. 133

5.7 The Block-Iterative (6, 17}-Aigorithm

• Under nonempty Q, fair weight vectors, and parameters in the stated positive unit range, the block-iterative algorithm converges to a feasible point.It combines approximate per-set projections within weighted blocks.Source: §5.7 · pp. 135–136

5.8 The Method of Successive Generalized Projections

• SGP replaces SOP’s orthogonal projections with successive generalized projections; under Assumption 5.8.1, its iterates converge into the nonempty intersection of finite closed convex sets.The algorithm starts from an arbitrary point.Source: §5.8 · p. 136

Algorithm 5.8.1 The Method of Successive Generalized Projections (SGP)

• Under the stated feasibility conditions and zone consistency, every SGP sequence converges to a point in Q.The SGP iteration uses cyclic generalized projections and allows no relaxation parameters.Source: Algorithm 5.8.1 The Method of… · pp. 136–137
• SGP replaces orthogonal projections with generalized projections, but only with unity relaxation; its sequence is DrFejer-monotone with respect to Q.The method leads to special-purpose optimization algorithms and allows different types of generalized projections within one algorithm.Source: Algorithm 5.8.1 The Method of… · pp. 137–138

5.9.1 The product space setup

• The convex feasibility problem is represented in product space as finding a point in the intersection of the product set and the diagonal subspace.Component Bregman functions are combined into an induced function using positive weights.Source: §5.9.1 · pp. 139–141

5.9.2 Generalized projections in the product space

• In product space, projection onto the product set is computed by independently projecting each component onto its corresponding set.The positive weighted sum of component D-functions is minimized when each component D-function reaches its minimum.Source: §5.9.2 · p. 141

5.9.3 The simultaneous multiprojections algorithm and the split feasibility problem

• Applying SGP to the product set and diagonal subspace with cyclic control alternates between those two sets.The induced Bregman function is assumed zone consistent with respect to both sets.Source: §5.9.3 · pp. 142–143

Algorithm 5.9.1 Successive Generalized Projections (SGP) for Two Sets in the Product Space

• The two-set product-space SGP construction is translated into an algorithm in R^n, initialized from a point in the diagonal set intersected with S.The translation uses the product-space membership condition and the component projection results.Source: Algorithm 5.9.1 Successive Generalized… · p. 143

Algorithm 5.9.2 The Simultaneous Multiprojections Algorithm

• Under the stated assumptions, the simultaneous multiprojections algorithm converges to a point in Q intersected with S.With squared Euclidean Bregman functions it is Cimmino’s fully simultaneous projections method; other choices yield new simultaneous algorithms.Source: Algorithm 5.9.2 The Simultaneous… · pp. 143–144
• The simultaneous multiprojections algorithm also addresses split feasibility: find x in C such that Ax is in Q, where A is a linear bijection.The stated applications include phase retrieval and image restoration or recovery.Source: Algorithm 5.9.2 The Simultaneous… · p. 144
• The simultaneous multiprojections algorithm accommodates two different projections within one iterative process and has been used to solve problems with these projection difficulties.The SOP-based formulation is guaranteed to converge, but practical applications may make one set difficult to characterize and one nonorthogonal projection difficult to perform.Source: Algorithm 5.9.2 The Simultaneous… · p. 145

5.10 Automatic Relaxation for Linear Interval Feasibility Problems

• The interval method handles each hyperslab as a pair of inequalities and incorporates automatic relaxation, with external relaxation parameters also available.It exploits the interval structure rather than treating the system as 2m one-sided inequalities for the AMS method.Source: §5.10 · p. 145

Algorithm 5.10.1 Automatic Relaxation Method (ARM)

• ARM is a row-action method that uses unchanged input data, accesses one system row per step, and needs only the preceding iterate, supporting use on very large, sparse systems.It works with slabs rather than their individual bounding hyperplanes and requires no matrix operations.Source: Algorithm 5.10.1 Automatic Relaxation… · p. 147
• ARM automatically varies its step within a slab: far from the slab, its move-back correction is small; near the slab, the iterate moves closer to a bounding hyperplane.Its relaxation strategy starts near one and increases toward two as the sequence approaches feasibility; external parameters can further control implementation.Source: Algorithm 5.10.1 Automatic Relaxation… · p. 147

Algorithm 5.10.2 Algebraic Reconstruction Technique 3 (ART3)

• ART3 addresses interval inequalities motivated by measurement or modeling errors, projecting to the median hyperplane when far from a slab and reflecting at a bounding hyperplane when nearer.It uses an enveloping hyperslab to distinguish distance regimes; the provided passage also describes the case when the iterate is inside the basic slab.Source: Algorithm 5.10.1 Automatic Relaxation… · p. 148, Algorithm 5.10.2 Algebraic… · p. 148

I I I I ' I ' : : : : ',, ------~-----. ------·------.. ------ '

• ART3 projects onto a bounding hyperplane outside the enveloping hyperplane, reflects within it, while ARM varies relaxation continuously outside the hyperslab.The figure compares Kaczmarz's method (ART), ART3, AMS, and ARM.Source: I I I I ' I ' : : : : ',,… · p. 150, Figure 5.11 · p. 150
• Every sequence produced by ARM converges to a solution of the interval linear inequalities problem, provided it is feasible.The proof reformulates the problem as convex feasibility and applies the CSP algorithm for convex inequalities.Source: I I I I ' I ' : : : : ',,… · p. 151

5.11 Notes and References

• The notes trace origins and extensions of feasibility methods, including ARM's proposal by Censor (1985), and point to reviews of the broader literature.They also identify Goffin's relaxation principle as the basis for this section's algorithms.Source: §5.11 · pp. 151, 152, 154, 155

Constrained Optimization Problems

• The chapter presents row-action methods for very large sparse systems and treats linear and nonlinear equalities and inequalities.Its unified treatment of linearly constrained optimization uses generalized Bregman distances and projections, and its roadmap includes norm and entropy optimization algorithms, underrelaxation, and hybrid methods.Source: Constrained Optimization Problems · pp. 156–157

6.1.1 The problem

• Real-world equation systems may be underdetermined, inconsistent, ill-conditioned, or corrupted by measurement, noise, or discretization error, making an exact solution less desirable.These conditions motivate choosing a different solution concept when appropriate.Source: §6.1.1 · p. 157

6.1.2 Approaches and solution concepts

• The user chooses a solution concept for the particular application, using the system's information in a way that reflects limited confidence in its equations.Different approaches yield different solution concepts.Source: §6.1.2 · p. 157

The Feasibility Approach

• The feasibility approach seeks a point within prescribed positive tolerances of every equation hyperplane, with tolerances chosen to avoid an empty or excessively large feasible region.The resulting interval feasibility problem can be treated as paired inequalities or as one-sided inequalities.Source: The Feasibility Approach · p. 158

The Optimization Approach

• The optimization approach selects a feasible point by an objective, favors methods that exploit slab constraints, and may use a secondary objective when solutions are nonunique.The regularization approach replaces the constrained system with minimizing f(x) + p·g(Ax − b), where p weights the two convex-function terms; least-squares and entropy regularizations are noted.Source: The Optimization Approach · p. 158, The Regularization Approach · p. 159

The Distance Reduction Approach

• Distance reduction minimizes a generalized Bregman distance from b to Ax, yielding an exact solution when feasible and otherwise a point in Q closest to R(A).The squared Euclidean choice gives least squares, while negative Shannon entropy yields a Kullback-Leibler problem addressable by EM.Source: The Distance Reduction Approach · pp. 159–160

6.1.3 The special computational environment

• The computational environment may involve huge, sparse, unstructured systems, urgent solution deadlines, and limited computing resources.Source: §6.1.3 · p. 160

6.2 Row-Action Methods

• Row-action methods update sequentially using one matrix row at a time, without altering or operating on the whole matrix, and can support row generation.The framework emphasizes modest arithmetic demands; cyclic and almost-cyclic controls are favored for practicality.Source: §6.2 · p. 161

6.3 Bregman's Algorithm for Inequality Constrained Problems

• Under strong zone consistency, nonempty Z0, and almost-cyclic control, Bregman’s algorithm converges to the solution; each iterate also satisfies the currently selected inequality.The convergence theorem additionally assumes f ∈ B(S); the iterate property states xv+1 ∈ Qi(v).Source: §6.3 · pp. 165–166
• Every sequence generated by Algorithm 6.3.1 converges to the unique solution of (6.12)–(6.14).Source: §6.3 · p. 169

6.4 Algorithm for Interval-Constrained Problems

• The interval method uses one dual variable per interval constraint and half as many iterative steps as applying Algorithm 6.3.1 to the doubled one-sided formulation.Source: §6.4 · p. 171

Algorithm 6.4.1 Row-Action Algorithm for Interval Convex Programming

• Algorithm 6.4.1 processes each interval’s paired inequalities in almost-cyclic order, choosing which side first from the current dual-variable sign; almost-cyclic control enables this flexibility.Source: Algorithm 6.4.1 Row-Action Algorithm… · p. 173
• Under strong zone consistency, Q ∩ S nonempty, nonzero interval widths, almost-cyclic control, and U ≠ 0, every generated sequence converges to a solution.Source: Algorithm 6.4.1 Row-Action Algorithm… · p. 173

6.5 Row-Action Algorithms for Norm Minimization

• The squared-norm specializations yield Kaczmarz for linear equalities, Hildreth for linear inequalities, and ART4 for linear interval constraints.These are special cases of the general iterative algorithms, Algorithm 6.3.1 and Algorithm 6.4.1.Source: §6.5 · p. 176

Algorithm 6.5.1 Kaczmarz's Algorithm for Norm Minimization

• Kaczmarz becomes a norm-minimization method when initialized in the range of A^T, but its general relaxation-parameter convergence requires separate treatment.The cited general iterative schemes cover the unrelaxed and underrelaxed cases, not the full range of relaxation parameters.Source: §6.5.1 · p. 176, Algorithm 6.5.1 Kaczmarz's Algorithm… · p. 177

Algorithm 6.5.2 Hildreth's Algorithm

• Hildreth’s algorithm updates primal and dual iterates while changing only one dual component per step, with the change chosen to guarantee dual ascent.For unity relaxation convergence follows from Algorithm 6.3.1; underrelaxed and general relaxation cases are addressed separately.Source: §6.5.2 · p. 178, Algorithm 6.5.2 Hildreth's Algorithm · pp. 178–179

6.5.3 ART4-An algorithm for norm minimization over linear intervals

• ART4 extends Hildreth’s algorithm to interval constraints, motivated by replacing inconsistent equalities with feasible intervals in image reconstruction.At each step it uses the median operation to select an update that places the iterate inside the slab or on one of its bounding hyperplanes.Source: §6.5.3 · pp. 179–180
• Every sequence generated by ART2 converges to a solution of the stated problem with interval and box constraints.ART2 handles box constraints by resetting any variable that violates a bound to its nearest bound.Source: §6.5.3 · p. 181

Algorithm 6.6.1 Bregman's Row-Action Method for Linear Equality Constrained Entropy Maximization

• Bregman’s row-action method for equality-constrained entropy maximization requires solving a nonlinear system at each iteration, creating a gap between the theoretical method and practical implementation.MART is offered as an alternative with closed-form iterative updates rather than an inner iterative solve.Source: Algorithm 6.6.1 Bregman's Row-Action… · p. 182

Algorithm 6.6.2 Multiplicative Algebraic Reconstruction Technique (MART)

• MART’s convergence analysis assumes feasibility, nonnegative nonzero constraint coefficients with positive right-hand sides, and normalization of A so its entries are at most 1.Its stated control sequence is almost cyclic, and its relaxation parameters are underrelaxed and bounded away from zero.Source: Algorithm 6.6.2 Multiplicative… · p. 183
• Under Assumptions 6.6.1–6.6.3, every MART sequence converges to the unique minimizer of the specified entropy optimization problem.The minimizer is for problem (6.106)–(6.107) with Q = Q1 as in (6.108).Source: Algorithm 6.6.2 Multiplicative… · p. 183

Algorithm 6.7.1 Block-Iterative MART

• With positive initialization, Algorithm 6.7.1 keeps every iterate positive and uses almost-cyclic block selection.The control sequence is almost cyclic over the set of blocks.Source: Algorithm 6.7.1 Block-Iterative MART · pp. 184–185

6. 7 Block-Iterative MART Algorithm

• The block-iterative MART partitions equality constraints into weighted blocks, and under the stated assumptions its sequence L increases monotonically.The blocks are selected by an almost-cyclic control sequence.Source: §6 · pp. 183–184, Algorithm 6.7.1 Block-Iterative MART · p. 185

Algorithm 6.7.1 Block-Iterative MART

• Under Assumptions 6.6.1–6.6.3, the block algorithm’s sequence L converges.Its monotonic increase is established separately; convergence follows from an upper bound.Source: Algorithm 6.7.1 Block-Iterative MART · p. 186
• Under the stated assumptions, block-algorithm iterates are bounded, and every accumulation point belongs to the feasible set Q1 ∩ R_+^n.The feasibility result applies to any accumulation point of a generated sequence.Source: Algorithm 6.7.1 Block-Iterative MART · p. 187

S I (ar,xs)- br I+ (ar,xs) · max I hj(v(C),s) -1 j, (6.146) lSJSn

• Under Assumptions 6.6.1–6.6.3, every sequence generated by Algorithm 6.7.1 converges to the unique solution of problem (6.106)–(6.107) with Q = Q1.The proof establishes optimality over Q1 ∩ R^n_+ and uses strict convexity of f(x) = −ent x to establish uniqueness.Source: S I (ar,xs)- br I+ (ar,xs) · max I… · pp. 188–189

6.8 Underrelaxation Parameters and Extension of the Family of Bregman Functions

• The extended Bregman-function framework permits essential boundary discontinuities by requiring partial level sets to be bounded and bounded away from the singularity set.The singularity set may be all, part, or none of the zone boundary; ordinary Bregman functions are included as the case S = S.Source: §6.8 · pp. 190–191

Definition 6.8.4 (Zone Consistency of Extended Bregman Functions)

• Strong zone consistency requires extended projections onto a hyperplane and the parallel hyperplanes between it and the point to remain in the zone.For a given hyperplane, the first-order conditions identify a unique projection and associated projection parameter.Source: Definition 6.8.4 (Zone Consistency of… · p. 192

Algorithm 6.8.1 Underrelaxed Algorithm for Extended Bregman Functions Over Linear Inequality Constraints

• Underrelaxation uses an unrelaxed extended projection onto a parallel hyperplane lying between the current point and the original hyperplane.The algorithm selects control indices almost cyclically; strong zone consistency keeps its iterates in S.Source: Definition 6.8.4 (Zone Consistency of… · p. 194, Algorithm 6.8.1 Underrelaxed Algorithm… · p. 194

Assumption 6.8.1

• Algorithm 6.8.1 maintains nonnegative dual iterates satisfying ∇f(xv) = −A^Tzv, and its Lagrangian values are nondecreasing.The stated assumptions include strong zone consistency for each constraint hyperplane and an almost-cyclic control sequence.Source: Assumption 6.8.1 · pp. 195–196
• The iterates are bounded when relaxation parameters satisfy only 0<alpha_v<=1, but the later convergence proof requires them to be bounded away from zero.The boundedness result for the iterates holds under the weaker relaxation condition.Source: Assumption 6.8.1 · p. 197
• Under Assumption 6.8.1, all accumulation points lie in Q∩S, and Algorithm 6.8.1 converges to a solution of (6.158).Strict convexity of problem (6.158)–(6.160) ensures a unique accumulation point.Source: Assumption 6.8.1 · pp. 197–198, Given p E h, define · p. 200

Given p E h, define

• For linear equalities, modified Algorithm 6.8.2 eliminates dual iterates and extends Bregman’s algorithm to underrelaxation and extended functions; its iterates generally differ unless relaxation is unity.Theorem 6.8.1 also holds for the equality-constraint case.Source: Given p E h, define · p. 200

6.9 The Hybrid Algorithm: A Computational Simplification

• The hybrid algorithm retains Bregman’s structure but replaces costly inner-loop projection calculations with a closed-form MART step, facilitating parallel computation across constraints.The passage characterizes MART as an approximation to Bregman’s algorithm; projection systems can burden computation and numerical errors can cause practical deviations.Source: §6.9 · pp. 201–203
• When x can be expressed as a solution function h(y,a,s), projection reduces to a scalar root; a suitable underrelaxation strategy may yield a closed-form update.The scalar-root formulation follows by substituting x=h(y,a,s) into the remaining projection equation.Source: §6.9 · p. 203
• When condition (6.215) holds, a closed-form formula replaces projection calculations and yields a computationally simpler algorithm.The hybrid construction is motivated by entropy maximization over linear constraints.Source: §6.9 · p. 204

Algorithm 6.9.1 The Hybrid Algorithm

• The hybrid algorithm retains Algorithm 6.8.1’s structure, uses a closed-form coefficient reminiscent of MART, and inherits convergence guarantees under specified underrelaxation choices.The algorithm also requires condition (6.215); additional matrix and vector conditions may be needed to obtain the closed-form formula.Source: Algorithm 6.9.1 The Hybrid Algorithm · p. 205

6.9.1 Hybrid algorithms for Shannon's entropy

• For Shannon’s entropy, conditions on A and b yield a closed-form coefficient for the hybrid method, with relaxation parameters allowed to vary across iterations.Rows with zero right-hand side can be removed beforehand because nonnegativity forces variables with nonzero coefficients in those rows to be zero.Source: §6.9.1 · p. 206

Algorithm 6.9.2 A Hybrid Algorithm for Shannon's Entropy

• With linear equality constraints and the stated assumptions, Algorithm 6.9.2 reduces to MART, which preceded and inspired the hybrid algorithms.The text identifies MART as the forerunner of these entropy-maximization hybrid algorithms.Source: Algorithm 6.9.2 A Hybrid Algorithm for… · p. 207

6.9.2 Algorithms for the Burg entropy function

• Burg’s entropy is an extended Bregman function; the section gives a necessary-and-sufficient solvability criterion and proves existence of a unique positive hyperplane projection under stated conditions.The solvability criterion applies when the feasible set and matrix are nonempty and nonzero; the projection result assumes a nonempty positive hyperplane intersection and y > 0.Source: §6.9.2 · pp. 207, 209

Algorithm 6.9.4 Hybrid Algorithm for Burg's Entropy Maximization Over Inequalities

• Under condition (6.254), the Burg-entropy hybrid algorithm uses a closed-form update coefficient, with relaxation parameters satisfying (6.227).The coefficient uses the stated endpoint formulas for the positive and negative right-hand-side cases.Source: §6.9.2 · p. 211, Algorithm 6.9.4 Hybrid Algorithm for… · p. 211
• Under condition (6.254), the hybrid Burg algorithm has 0 < α_v ≤ 1 for all v, and its iterates remain bounded and bounded away from the singularity set.The notes attribute the hybrid algorithm to Censor, De Pierro, and Iusem (1991).Source: Algorithm 6.9.4 Hybrid Algorithm for… · p. 212, §6.10 · p. 218

6.9.3 Renyi's entropy function

• Rényi entropy is an extended Bregman function with a singularity at zero; replacing the objective with the reciprocal-power alternative removes that singularity and makes Algorithms 6.8.1 and 6.9.1 applicable.The alternative is described as preferable for practical linearly constrained Rényi entropy optimization.Source: §6.9.3 · pp. 213–214

Model Decomposition Algorithms

• Model decomposition handles large sparse or linked problems by optimizing components separately while coordinating their solutions through a master problem.Examples include time periods linked by inventory flows and geographic regions linked by trade or traffic.Source: Model Decomposition Algorithms · p. 219
• Parallel Dantzig-Wolfe and Benders methods often achieve modest speedups because a serial master program becomes an increasingly restrictive bottleneck as iterations proceed.The cited parallel implementations generally solve subproblems in parallel but leave the master program on one processor.Source: Model Decomposition Algorithms · p. 220

7.1 General Framework of Model Decompositions

• A model decomposition removes complicating constraints from a modified problem, then iteratively updates that problem until its solution is sufficiently close to the original problem’s solution.The modifier and algorithm for solving the modified problem together determine whether the decomposition is suitable for parallel computation.Source: §7.1 · p. 221

7.1.1 Problem modifiers

• The first problem modifier eliminates complicating constraints by adding a penalty or barrier function to the objective.The modified objective is F(x) + cp(x).Source: §7.1.1 · p. 222
• Exact penalties can make solutions of the modified problem solutions of the original, while barrier problems with decreasing parameters converge to an original-problem solution.The penalty removes complicating constraints but may remain nonseparable; barrier iterates remain in the interior when initialized there, and the barrier function may also be nonseparable.Source: §7.1.1 · pp. 222–223
• Variable splitting replicates each block and uses an augmented Lagrangian to eliminate the linking constraints, yielding a dual problem solvable by the method of multipliers.The multiplier update is simple, but the quadratic cross-products keep the augmented Lagrangian nonseparable; its update operations can be executed efficiently in parallel.Source: §7.1.1 · pp. 223–224

7.1.2 Solution algorithms

• Linearization and diagonal approximation turn nonseparable objectives with block-decomposable constraints into K subproblems that can be solved in parallel.Linearization is used for nonlinear functions; diagonal approximation is used for nonseparable quadratic functions such as the augmented-Lagrangian term.Source: §7.1.2 · pp. 224–225

Algorithm 7.1.2 The Frank-Wolfe Algorithm

• Frank-Wolfe linearizes the objective at the current iterate, solves a linear program for a descent direction, then uses a one-dimensional line search.The linearized subproblem decomposes into K independent linear programs over the block constraint sets, enabling parallel solution.Source: Algorithm 7.1.2 The Frank-Wolfe… · p. 226
• Frank-Wolfe can zigzag and cannot reach an optimum on a constraint-set facet by a line search between a vertex and an interior point.Simplicial decomposition addresses this limitation by using multiple generated vertices in a nonlinear master program.Source: Algorithm 7.1.2 The Frank-Wolfe… · pp. 226–227

Algorithm 7.1.3 The Simplicial Decomposition Algorithm

• Simplicial decomposition avoids Frank-Wolfe zigzagging by optimizing over a master program on generated vertices’ convex hull, which is usually much smaller than the original problem.Vertices are generated through linear-programming subproblems; vertices with zero weight can be removed from the master program.Source: Algorithm 7.1.2 The Frank-Wolfe… · p. 227, Algorithm 7.1.3 The Simplicial… · p. 228

Solution Algorithms Based on Diagonalization

• Diagonalization approximates the nonseparable quadratic terms ||xk - zk||^2 in the augmented Lagrangian with a separable quadratic function.It approximates the Hessian by a diagonal matrix, using first-order approximations of the cross-products around the current iterate.Source: Solution Algorithms Based on… · p. 229
• Diagonalization replaces augmented-Lagrangian minimization with a diagonal quadratic approximation that separates into K independent block problems solvable in parallel.The algorithm is described as using the method of multipliers with this separable approximation in Step 1.Source: Solution Algorithms Based on… · pp. 229–230

7.2 The Linear-Quadratic Penalty (LQP) Algorithm

• The LQP algorithm smooths the nondifferentiable exact penalty, enabling gradient-based methods and offering a natural way to handle soft constraints that need not be satisfied exactly.The smoothing is uniform as its parameter tends to zero; it can introduce approximation error, for which the section says an a priori bound can be computed.Source: §7.2 · p. 231

Algorithm 7.2.1 The linear-Quadratic Penalty (LQP) Algorithm

• The LQP algorithm repeatedly minimizes the smoothed penalty problem, stopping when the iterate meets the current feasibility tolerance and that tolerance is no greater than the user-set final threshold.If the iterate is feasible but the tolerance remains too large, the algorithm reduces the tolerance; if it is infeasible, it increases the penalty parameter.Source: Algorithm 7.2.1 The linear-Quadratic… · pp. 232–233

7.2.1 Analysis of the e-smoothed linear-quadratic penalty function

• The exact and smoothed penalty objectives differ by at most Lμε^2 for any x, μ > 0, and ε ≥ 0.The bound is Proposition 7.2.1’s a priori bound on smoothing error.Source: §7.2.1 · p. 233

L L 0 :S LP(9t(x))- LPE(gt(x))::; LJ.L~, for all x E IRnK,

• For an ε-feasible minimizer of the smoothed problem, the original objective gap is bounded above by Lμε.The proposition compares an optimal solution of the original constrained problem with an optimal solution of the smoothed penalty problem.Source: L L 0 :S LP(9t(x))- LPE(gt(x))::;… · p. 234

L L 0 ~ Uo(x*) + LP(YI(x*)))- Uo(x) + Lp,(gz(x))) ~ LJ.Li· (7.52)

• When μ is at least each complicating-constraint multiplier z_l, the smoothed penalty value at the constrained optimum is within Lμε^2 of its value at any x in X.This bounds the constrained optimum’s suboptimality for the smoothed penalty problem even though the two optima need not coincide.Source: L L 0 ~ Uo(x*) + LP(YI(x*)))- Uo(x) +… · p. 235

L L - L z[(gt(x)- gt(x*)) + J.L L g((x)

• At the original optimum x*, the smoothed penalty objective exceeds the true optimal value by at most Lμε/2.As ε vanishes, the classical exact-penalty result is recovered: an original optimum also optimizes the penalty problem for some μ.Source: L L - L z[(gt(x)- gt(x*)) + J.L L g((x) · p. 237

7.2.2 e-exactness properties of the LQP function

• The smooth LQP penalty inherits approximate exactness from its nonsmooth counterpart, attaining ε-feasibility for a finite penalty value.The section relates optimality conditions for the penalized and original problems, with feasibility and complementary slackness as exceptions.Source: §7.2.2 · p. 238, Optimality Conditions for (7.45) · p. 239

Optimality Conditions for (7.45)

• A penalty parameter μ ≥ μ*/τ_ε ensures an ε-feasible solution; for a given problem, the threshold is independent of ε.The threshold μ* exceeds each Lagrange multiplier z_l.Source: Optimality Conditions for (7.45) · p. 239
• The algorithm increases the penalty parameter when a solution is not ε-feasible, decreases smoothing after successful minimizations, and stops when ε is below the accuracy tolerance.At termination, Proposition 7.2.8 supplies an upper bound on the difference between the original and smoothed-penalty optimal values.Source: Optimality Conditions for (7.45) · p. 242

7.3 Notes and References

• Notes trace decomposition methods from Dantzig (1963), Geoffrion (1970), and Lasdon (1970) through parallel Dantzig-Wolfe and Benders methods to Ferris and Mangasarian’s (1991, 1994) constraint-distribution methods.The notes also cite parallelizations of Dantzig-Wolfe (Dantzig and Wolfe (1960)) and Benders decomposition (Benders (1962)).Source: §7.3 · pp. 242–243

Algorithms

• Parallel matrix factorization can solve the systems required by primal-dual path-following algorithms for structured problems, improving the performance of an already efficient method.The chapter treats linear and quadratic programs and develops parallel factorization procedures for their required systems.Source: Algorithms · pp. 246–247

8.1 The Primal-Dual Path Following Algorithm for Linear Programming

• As μ tends to zero, the central path converges to primal and dual optima, and both objective values converge to their common optimal value.For each μ > 0, the path is defined by the simultaneous primal-dual barrier solution.Source: §8.1 · p. 249

Algorithm 8.1.1 The Primal-Dual Path Following Algorithm for Linear Programs

• Starting from positive primal and slack variables, the method tracks the central path with Newton steps, linearizes complementarity, and solves a symmetric positive definite system for the dual direction.It then recovers the primal and slack directions from the dual direction.Source: §8.1 · pp. 249–250

8.1.1 Choosing the step lengths

• 0.9995 is a typical fraction-to-boundary factor for step lengths chosen to keep primal and slack variables positive.The primal and dual step lengths are fractions of the respective maximum lengths that preserve positivity.Source: §8.1.1 · p. 251

8.1.2 Choosing the barrier parameter

• A duality-gap-based practical update reduces μ rapidly near the optimum, whereas the theoretical choice that assures polynomial convergence is very slow in practice.The practical update keeps μ large when the duality gap indicates the iterate is far from optimal.Source: §8.1.2 · p. 252

8.2 The Primal-Dual Path Following Algorithm for Quadratic Programming

• For quadratic programs, logarithmic barriers remove inequality constraints and Newton's method solves the barrier optimality conditions; as μ tends to zero, feasible sequences converge to primal and dual solutions.For separable quadratic functions, the linear- and quadratic-program algorithms are virtually identical.Source: §8.2 · pp. 252–253
• Under boundedness and nonempty-interior feasibility, the primal-dual system has a unique positive solution for each μ > 0, and the algorithm tracks its central path as μ decreases.Source: §8.2 · p. 254

Algorithm 8.2.1 The Primal-Dual Path Following Algorithm for Quadratic Programs

• Algorithm 8.2.1 starts with positive x and z, decreases the barrier parameter each iteration, and uses the minimum of primal and dual step lengths to preserve positivity.Source: Algorithm 8.2.1 The Primal-Dual Path… · pp. 255–256

8.3 Parallel Matrix Factorization Procedures for the Interior Point Algorithm

• For block-angular constraints and diagonal scaling, the factorization exploits structure arising in stochastic programs and multicommodity network flows, where systems can be very large.The source notes that the product matrix AΘA^T may nevertheless be dense.Source: §8.3 · pp. 256–257

8.3.1 The matrix factorization procedure for the dual step direction calculation

• A generalized Sherman-Morrison-Woodbury formula underpins the factorization, expressing the dual-step solution through block solves and a smaller structured system.Source: §8.3.1 · pp. 257, 259

Procedure 8.3.1 Matrix factorization for dual step direction calculation

• Procedure 8.3.1 solves independent block systems, forms G1 and G2, and then solves the resulting dense systems to calculate the dual step direction.Source: Procedure 8.3.1 Matrix factorization… · pp. 259–260

Decompositions for Parallel Computing

• Parallel decomposition assigns local block computations to processors, with communication to assemble and distribute shared quantities; alternatives replicate dense calculations or parallelize dense linear algebra.Source: Decompositions for Parallel Computing · pp. 260–261

calculation

• The calculation routine performs block solves in parallel, handles dense G1 and G2 systems on the master processor, broadcasts q, and gathers the final dual-step vector.Source: calculation · p. 261

8.4 Notes and References

• The notes credit Karmarkar (1984) with proposing projective scaling and establishing polynomial complexity, while identifying affine scaling as earlier work by Dikin (1967, 1974).Source: §8.4 · pp. 261–262
• The notes cite work on primal-dual path-following algorithms, implementation, parallel matrix factorization, stochastic programming, and optimization software.For example, they credit Monteiro and Adler (1989a) and Kojima, Mizuno, and Yoshise (1989) with independently developing the LP algorithm, and describe OB1 in Lustig, Marsten, and Shanno (1991).Source: §8.4 · pp. 262–263

Matrix Estimation Problems

• Matrix estimation adjusts a given matrix to satisfy linear restrictions, with matrix balancing as a special case across diverse application areas.The chapter’s applications include economics, transportation, telecommunications, statistics, demography, and stochastic modeling.Source: Matrix Estimation Problems · p. 266, §9.1 · p. 267, §9.2.1 · p. 273

9.1.1 Economics: social accounting matrices (SAMs)

• Raw social accounting matrix estimates are generally unbalanced, so they must be adjusted to satisfy accounting identities equating each agent’s total expenditures and receipts.Balancing also updates an older matrix to match new row and column totals when current transaction-level data are unavailable.Source: §9.1.1 · pp. 267–268

9.1.2 Transportation: estimating Origin-Destination Matrices

• Origin-destination matrices estimated from observed arc flows are chosen to match those flows while minimizing distance from a reference matrix.When each origin-destination pair has a single route, estimation from origin outflows and destination inflows can be formulated as matrix balancing.Source: §9.1.2 · p. 269

9.1.3 Statistics: estimating contingency tables

• Contingency-table cell probabilities are adjusted to agree with known marginal distributions when direct population sampling is prohibitively expensive.The passage also notes extensions to four-way tables with available one-, two-, or three-way margins.Source: §9.1.3 · p. 270

9.1.4 Demography: modeling interregional migration

• Migration matrices reconcile outdated interregional flow patterns with newer regional net in- and out-migration figures.The doubly constrained minimum-information model yields a balancing problem, sometimes with an additional constraint on total migrant travel distance.Source: §9.1.4 · pp. 270–271

9.1.5 Stochastic modeling: estimating transition probabilities

• Transition probabilities can be estimated from observed state proportions by minimizing quadratic discrepancy across periods, subject to nonnegative entries and row sums of one.When observations outnumber states, the passage says an exact solution generally does not exist; the resulting estimation problem uses a discrepancy objective rather than a prior transition matrix.Source: §9.1.5 · pp. 271–272

9.2 Mathematical Models for Matrix Balancing

• The chapter’s formulations specify constraints on matrix integrals and how the estimated matrix relates to the given matrix, distinguishing matrix balancing from broader matrix estimation.Considered constraints include prescribed row and column sums, entry bounds, total-sum conditions, and equality between row and column sums; relationship types include form relations.Source: §9.2 · pp. 272–273, §9.2.1 · p. 273

9.2.1 Matrix estimation formulations

• Matrix estimation is framed either axiomatically or by minimizing a distance from the original matrix, including Frobenius-norm and negative-entropy objectives.The Frobenius choice yields a linearly constrained quadratic optimization problem.Source: §9.2.1 · p. 274

L Xij (log (:'1 ) - 1) . (i,j)EA ZJ

• The two common balancing tasks seek a nearby matrix with prescribed row and column sums or, for a square matrix, equal row and column sums.Social accounting matrices are given as examples of both tasks; the consistency restrictions admit infinitely many matrices, so a notion of nearness is needed.Source: L Xij (log (:'1 ) - 1) . (i,j)EA ZJ · pp. 274–275

9.2.2 Network structure of matrix balancing problems

• Matrix balancing maps to transportation on a bipartite graph for the row-and-column-sum task and to transshipment on a square-matrix graph for the flow-conservation task.Nonzero matrix entries correspond to arcs, and the transshipment flow-conservation conditions are the consistency requirements for Problem 9.2.2.Source: §9.2.2 · pp. 275–277

9.2.3 Entropy optimization models for matrix balancing

• Entropy optimization treats a matrix as nearby when its entropic distance from the original is minimized, a formulation justified as least biased and widely used with network-flow constraints.Variables for zero entries of the original matrix are omitted, and the corresponding estimated entries remain zero.Source: §9.2.3 · p. 277
• Interval-constrained entropy models accommodate uncertain row and column totals or relaxed row-column balance by bounding sums or their differences.The models impose lower and upper bounds on row and column sums, or a positive tolerance on each row-column difference.Source: §9.2.3 · pp. 278–279

9.3 Iterative Algorithms for Matrix Balancing

• Because entropy is a Bregman function, the chapter specializes a general row-action algorithm to interval constraints, then obtains the equality-constrained algorithms as special cases.The interval formulations are Problems 9.2.5 and 9.2.6; the equality formulations are Problems 9.2.3 and 9.2.4.Source: §9.3 · p. 279

9.3.1 The range-RAS algorithm (RRAS)

• Range-RAS solves Problem 9.2.5 by adapting row-action interval convex programming; its dual ascent performs exact minimization because the constraint matrix is 0-1.Source: §9.3.1 · p. 279

Algorithm 9.3.1 The range-RAS (RRAS) Algorithm

• The RRAS procedure initializes the entries from A, then repeatedly applies row scaling and column scaling.The algorithm specifies initialization in Step 0, row scaling in Step 1, column scaling in Step 2, and returns to Step 1.Source: Algorithm 9.3.1 The range-RAS (RRAS)… · pp. 282–283

Decompositions for Parallel Computing

• RRAS row scaling and flow updates can run concurrently across rows and incident arcs, and column scaling parallelizes similarly.The scaling parameters for different rows are independent because their corresponding nodes share no arcs; row scaling can use up to m processors.Source: Decompositions for Parallel Computing · p. 283

9.3.2 The RAS scaling algorithm

• For fixed row and column sums, RAS alternates row and column scaling to force each node’s incident-flow total to its desired value.With the stated input, the range-RAS formulation coincides with classic RAS; the algorithm initializes from a_ij, scales rows, then columns, and repeats.Source: §9.3.2 · p. 283, Algorithm 9.3.2 The RAS Algorithm · pp. 283–284

9.3.3 The range-DSS algorithm (ROSS)

• Range-DSS derives from Algorithm 6.4.1 for interval-constrained entropy optimization and updates flows using a cyclic sequence of row or column indices.At each iteration, it obtains scaling parameters from nonnegative roots of quadratic equations before applying the update.Source: §9.3.3 · p. 284

Algorithm 9.3.3 The range-DSS (ROSS) Algorithm

• The RDSS algorithm selects a cyclic index, computes sums and nonnegative scaling roots, updates incident flows, and repeats.The update applies the computed parameters to flows incident to the selected row or column’s graph node.Source: Algorithm 9.3.3 The range-DSS (ROSS)… · p. 286

Decompositions for Parallel Computing

• In RDSS, scaling-parameter computation is serial across selected nodes, but sums and updates for the selected node’s incident arcs can be parallelized.Computing the sums requires communication among processors.Source: Decompositions for Parallel Computing · p. 286
• Graph coloring enables concurrent RDSS steps on nonadjacent nodes, and a color-based almost cyclic control sequence is permitted because convergence does not require cyclic control.Nodes with no shared arcs can have their scaling parameters computed and incident flows updated concurrently.Source: Decompositions for Parallel Computing · pp. 286–287

9.3.4 The diagonal similarity scaling (DSS) algorithm

• With f = 0, RDSS reduces to DSS for Problem 9.2.4; DSS adjusts incident flows to satisfy node conservation and can be parallelized similarly.The DSS algorithm computes scaling parameters and updates iteratively under a cyclic control sequence.Source: §9.3.4 · p. 287, Algorithm 9.3.4 The DSS Algorithm · pp. 287–288

9.4 Notes and References

• Censor and Zenios (1991) developed interval-constrained entropy algorithms and a massively parallel range-RAS implementation; the notes trace RAS and DSS antecedents and parallelization approaches.The notes identify Kruithof (1937) and Deming and Stephan (1940) among independent developers of RAS, and describe DSS extensions by Schneider (1989, 1990).Source: §9.4 · pp. 289–290

Image Reconstruction from Projections

• In practice, reconstruction uses discrete, finite-ray estimates of projection data, while the linear integral model is an idealized but tractable starting point despite measurement error and possible nonlinearity.The idealized continuous problem is to recover f from p; Radon's inversion formula expresses f in terms of p.Source: Image Reconstruction from Projections · p. 293

10.1 Transform Methods and the Fully Discretized Model

• The fully discretized transmission CT model represents the image with pixel values and each measured ray as a weighted sum, yielding the matrix system y = Ax.The pixel weights are the lengths of ray intersections with pixels, and the model uses an m × n projection matrix.Source: §10.1 · pp. 294–295
• The series-expansion model is flexible: functionals can differ from line integrals and basis functions can vary, with attenuated line integrals in ECT producing nonlinear equations for simultaneous activity and attenuation estimation.Transform methods keep a continuous formulation until implementation, whereas the series-expansion approach discretizes at the outset.Source: §10.1 · pp. 296, 301
• Given system features such as dimensions, sparsity, inconsistency, and ill-conditioning, users choose among feasibility, optimization, distance minimization, and regularization; this choice determines which mathematical solutions count for image reconstruction.Source: §10.1 · pp. 298–299

10.2 A Fully Discretized Model for Positron Emission Tomography

• The chapter next examines a fully discretized model for positron emission tomography (PET).The chapter outline identifies PET as Section 10.2.Source: Image Reconstruction from Projections · p. 293, §10.2 · p. 301
• The PET model treats voxel concentrations as constant Poisson emission rates and uses detector-tube counts and geometry-based detection probabilities to form a likelihood.It approximates the physics: positron travel, non-180° photon angles, and tissue absorption are omitted and require model extensions.Source: §10.2 · pp. 302, 304

10.2.1 The expectation-maximization algorithm

• The relaxed EM algorithm converges to a solution for α in (0, 2); convergence holds from any positive start for α in (0, 1], but is only local for α in (1, 2).The unrelaxed iteration converges from any positive starting vector.Source: §10.2.1 · p. 305

10.3 A Justification for Entropy Maximization in Image Reconstruction

• For noisy or underdetermined tomography measurements, the reconstruction is framed as maximizing the posterior probability P(x|y), equivalently the likelihood times the prior.The measurement system may have no solutions or infinitely many solutions.Source: §10.3 · p. 306

Assumption 10.3.1

• Under the stated assumptions, maximizing the posterior reduces to a strictly convex, coercive entropy-plus-error objective with a unique minimizer.The model assumes equally likely fine configurations, relates pixel values to particle fractions, and takes independent zero-mean Gaussian measurement errors.Source: Assumption 10.3.1 · pp. 307–308

10.4 Algebraic Reconstruction Technique (ART) for Systems of Equations

• Variable-block ART unifies row-action ART, block-Kaczmarz, and Cimmino SIRT while allowing blocks to change dynamically in number, size, and structure.The weights attached to equations must not vanish.Source: §10.4 · p. 310

Algorithm 10.4.1 Variable-block ART-The General Scheme

• For a consistent system, every sequence generated by variable-block ART converges to a solution when the weight vectors satisfy property (10.60).Source: Algorithm 10.4.1 Variable-block… · p. 311
• Variable-block ART can change the number, sizes, and composition of equation blocks at each step, enabling dynamic block strategies and multilevel reconstruction.The block strategy may vary during iteration to accelerate initial convergence toward an acceptable reconstructed image; the method also includes row-action, Cimmino-type, and fixed-block cases.Source: Algorithm 10.4.1 Variable-block… · pp. 311–312

10.5 Iterative Data Refinement in Image Reconstruction

• Iterative data refinement repeatedly estimates the discrepancy between actual and ideal measurements, with a few iterations potentially reducing the initial discrepancy.Source: §10.5 · p. 312

10.5.1 The fundamentals of iterative data refinement

• If U is linear with bounded norm ||U||<1, IDR converges regardless of initialization; iterative result refinement can converge even when the measurement sequence does not.The latter distinction is possible when the linear recovery operator S has a nontrivial null space.Source: §10.5.1 · pp. 314–317
• When successive iterates remain in a targeted-contraction set, each is closer to the ideal measurement than its predecessor, though choosing a suitable operator is application-dependent.For ill-posed problems, the amount of constraining or prior knowledge needed for a reasonable output is also application-dependent.Source: §10.5.1 · p. 316
• With invertible linear A, IDR corresponds to the parallel-chord iteration, and its convergence is guaranteed when B is linear and ||I − μBA^-1||<1.When the ideal and actual measurement operators act as intended on the ideal result, IDR converges to the idealized measurement.Source: §10.5.1 · p. 318
• Starting from y0∈Gβ for 0<β<1, IDR either converges to ideal measurement f or moves closer to f than any element of G~.A stronger sufficient condition involving y0∈Gβ for some β<1 ensures convergence, but the authors say it is unlikely to hold in practice because U only estimates f−μy.Source: §10.5.1 · pp. 320–321

10.5.2 Applications in medical imaging

• The medical-imaging framework covers beam-hardening correction, SPECT attenuation correction, and reconstruction with incomplete data; forcing closed-form reconstructions into a pixel basis adds approximation error.Pixel bases are necessary for display, but the authors identify this added error as avoidable when using the natural basis.Source: §10.5.2 · p. 322
• ||fj−y0||=2.950 initially; in the illustrative beam-hardening run, improvement is rapid at first, slows, then the error slightly worsens by iteration 3.The authors note that this simple example’s behavior matches beam-hardening correction methods reported in the literature.Source: §10.5.2 · pp. 325–326, Table 10.2 · p. 326
• In the beam-hardening example, the actual device operator differs from the modeled operator even at the ideal image, so IDR can promise initial improvement but not convergence to y.The example’s tissue values cannot all be recovered from a pixel value shared by tissue types a and b.Source: §10.5.2 · pp. 325–326

Attenuation Correction in Single Photon Emission Computed

• For SPECT, body attenuation prevents direct measurement of activity line integrals, while the reported method produced an excellent approximation to x at the fifth iterate.The method can treat BS as a single operator, avoiding explicit production of a pixel-based version of each Syv.Source: Attenuation Correction in Single… · pp. 326–327

Image Reconstruction with Incomplete Data

• In realistic simulations, IDR reconstructions of x were very good after about 5 to 10 iterations, with incomplete CT data handled by estimating missing components of Ax.The known components are retained, while missing components are replaced using the current reconstruction; incorporating prior knowledge such as positivity can improve practical usefulness.Source: Image Reconstruction with Incomplete… · pp. 327–328

10.6 On the Selective Use of Iterative Algorithms for Inversion Problems

• Iterative algorithms are potentially advantageous when analytic inversion is unavailable or transform methods are not efficiently implementable, including sparse or limited-view CT data.The section distinguishes whether an operator or its inverse has a closed-form representation from whether the inversion can be implemented efficiently; iterative methods may be the only viable approach in some cases.Source: §10.6 · pp. 328, 330
• When both transform and iterative methods are implementable, their relative performance depends on the object and the application; an iterative method may still be pursued if its recovery is no worse but takes more computation.For example, the Radon transform in x-ray CT has an analytic inverse, whereas the attenuated Radon transform in SPECT has a known analytic inverse only for constant attenuation.Source: §10.6 · pp. 330–331

10.7 Notes and References

• Radiation-therapy planning motivates full discretization and row-action or block-iterative methods, with the aim of automating acceptable treatment setups from physician specifications.The planning problem can be difficult when tumors lie near critical organs, and the acceptable dose range between tumor control and normal-tissue injury is narrow.Source: §10.7 · pp. 335, 337

11.1 Problem Definition and the Continuous Model

• The treatment-planning model begins by describing the gantry-motion plane in polar coordinates and defining dose distribution D(r, θ) as absorbed dose at a patient-section point.The angle is measured clockwise from the positive y-axis.Source: §11.1 · p. 337

11.1.1 The continuous forward problem

• The forward problem maps a specified radiation-intensity distribution to patient dose, but its dose operator lacks a closed-form analytic representation for accurate practical calculation.The approach uses a computer program, incorporating empirical data such as lookup tables, rather than simplifying the physics to obtain a closed-form formula.Source: §11.1.1 · pp. 338–340

11.1.2 The continuous inverse problem

• Treatment planning seeks a radiation field that meets prescribed target doses while respecting maximum doses to critical organs and other tissues.Because some prescribed dose distributions may be unattainable, the inverse problem is formulated as satisfying upper and lower dose bounds rather than an exact equality.Source: §11.1.2 · pp. 340–341

11.2 Discretization of the Feasibility Problem

• The model discretizes the patient into pixels and represents the radiation field as nonnegative combinations of unit-intensity rays, yielding linear interval dose constraints.The discretization neglects scatter and apportions primary dose from each source among rays that pass through or immediately straddle each pixel.Source: §11.2 · pp. 341, 342, 343, 345

11.3 Computational Inversion of the Data

• The computational procedure precalculates and apportions dose coefficients, then iteratively adjusts ray intensities to satisfy physician constraints; rays missing tumor regions remain at zero.The stated procedure initializes all ray intensities to zero and omits scatter in its forward dose calculation.Source: §11.3 · p. 347
• The basic ray-intensity solution serves as input to beam reduction, followed by correction for scattered radiation after the beams are reduced.The original dose calculations account for primary dose only.Source: §11.3 · p. 348

11.4 Consequences and Limitations

• Present finite-width machines cannot accurately deliver finely discretized pencil-ray fields; implementable plans seek at most about 10 source columns with linear or piecewise-linear intensity changes.The described current approach selects source positions with the greatest ray-intensity contributions and repeats the feasibility calculation for a limited set of positions.Source: §11.4 · p. 349

11.5 Experimental Results

• The experimental model applies fully simultaneous Cimmino to pixel dose constraints, using 24 beam positions and 50 rays per beam with fixed regional weights.Weights are 10 for the target and eyes, 3 for the brain stem, and 1 for unlabelled tissue.Source: §11.5 · pp. 349–350

Right eye Left eye

• Repeated Cimmino screening reduced the 24-beam plan to beams at 60° and 120°, matching the dosimetrist’s two angles; dose-prescription errors appeared as beams were eliminated.The six-beam and two-beam configurations are shown with dose distributions and violations, but the supplied captions do not give numerical values.Source: Right eye Left eye · pp. 352–353
• Beam-reduction screening used ray-intensity sums first and adjacent-intensity smoothness to distinguish beams with similar sums.The initial 24 beams were narrowed to 11, then six, four, three, and finally two by repeated Cimmino applications.Source: Right eye Left eye · pp. 352–353

11.6 Combination of Plans in Radiotherapy

• Plan combination addresses cases where existing approximate plans each violate a significant dose region, seeking a jointly acceptable and technically realizable plan.If no such combination exists, the treatment requirements can be considered for systematic relaxation.Source: §11.6 · p. 357

11.6.2 The feasible case

• An acceptable plan combination exists if and only if the dual linear feasibility problem (11.21) is feasible, which simplex Phase 1 can test.The passage also identifies simplex or an interior-point algorithm as options for solving the linear programs.Source: §11.6.2 · pp. 360–361
• The organ-dose-minimizing combination had objective value 185.31 Gy; the target-dose-maximizing alternative delivered 54.4 Gy to the target.Reducing the target lower bound from 54 Gy to 53 Gy reduced total radiation to non-target organs by 4.12 Gy.Source: §11.6.2 · p. 362

11.6.3 The infeasible case

• An infeasible plan can be made feasible by relaxing dose bounds, including lowering the target minimum from 54 Gy to 53.7 Gy or raising the fourth constraint’s upper bound from 44 Gy to 44.2 Gy.The latter change yields y = (14.73, 24.55, 39.27, 44.18, 54)T.Source: §11.6.3 · p. 363

Algorithm 11.6.1 Simultaneous ART3

• With suitable relaxation parameters, simultaneous ART3 converges globally when feasible and locally, when infeasible, to a priority-weighted least-squares point.With equal weights, one resulting plan is y = (14.75, 24.54, 39.25, 44.13, 44.15, 53.9)T; weights (0.05, 0.05, 0.05, 0.05, 0.05, 0.75) yield y = (14.8, 24.62, 39.35, 44.23, 44.27, 53.99)T.Source: Algorithm 11.6.1 Simultaneous ART3 · pp. 365–367

11.7 Notes and References

• The text distinguishes its fully discretized inverse radiation-treatment approach from fixed-plan comparisons and sequentially varying a few parameters, using discretization to handle continuous-modeling difficulties.The patient section and radiation sources are both fully discretized.Source: §11.7 · p. 369

Multicommodity Network Flow Problems

• The chapter’s network-flow methods target parallel computation, emphasizing quadratic or entropic objectives while allowing linear-cost formulations through the PMD framework.The linear-cost combination is investigated empirically in Section 15.6.Source: Multicommodity Network Flow Problems · p. 371

12.1 Preliminaries

• A network’s node-arc incidence matrix encodes its structure; networks are capacitated when arcs have flow upper bounds and are pure when their arc coefficients equal 1.Networks with other coefficients are called generalized.Source: §12.1 · pp. 371–372

12.2 Problem Formulations

• The general flow formulation uses a convex, continuously differentiable objective, flow-conservation constraints, variable bounds, and a network-structured matrix that may be block diagonal.It also permits complicating nonnetwork constraints, and specializations include transportation and multicommodity network-flow problems.Source: §12.2 · p. 373

12.2.1 Transportation problems

• A single-commodity transportation problem is a bipartite network-flow special case with origin supplies, destination demands, and bounded arc flows.Its objective is typically separable, with convex, continuously differentiable costs for individual arcs.Source: §12.2.1 · pp. 373–375

12.2.2 Multicommodity network flow problems

• Multicommodity network flows conserve each commodity separately while coupling them through shared arc capacities, represented by generalized upper-bounding constraints.The model minimizes the sum of commodity costs subject to each commodity’s flow-conservation and individual bounds, plus joint capacity limits.Source: §12.2.2 · pp. 375–378

12.3.1 Example one: Covering positions in stock options

• Minimum stock-option margin can be formulated as a transportation problem that chooses which securities to pair, accounting for margin reductions from hedging options.The model pairs long and short positions using pairwise margin requirements to determine the minimum margin.Source: §12.3.1 · pp. 378–379

12.3.2 Example two: Air-traffic control

• Air-traffic network optimization balances cost-effective route allocation against congestion reduction, subject to route limits that maintain safety.Congestion depends nonlinearly on the aircraft using a route, and delays may be needed to avoid congestion, especially at peak times.Source: §12.3.2 · p. 379

12.3.3 Example three: Routing of traffic

• Communication traffic must conserve flows by origin-destination pair, while shared links create joint capacity limits and congestion-dependent, nonseparable delays.These interactions link the commodities through the cost function, producing traffic-assignment problems.Source: §12.3.3 · p. 380

Flow Problems

• The row-action algorithms are tailored to multicommodity transportation, extending from single-commodity quadratic transportation to entropic network objectives and multicommodity network flows.The general row-action method can be simplified for pure networks by setting arc multipliers to 1; linear objectives can be handled by embedding the algorithms in the PMD framework.Source: Flow Problems · p. 381

12.4.1 Row-action algorithm for quadratic transportation problems

• For quadratic transportation, the row-action method projects an iterate onto a node’s flow-conservation constraint by adjusting flows on incident arcs and updating the associated dual variable.If the current iterate already satisfies the constraint, the projection leaves it unchanged.Source: §12.4.1 · pp. 381–382, S;- L Yij · p. 382

S;- L Yij

• The quadratic row-action method also projects flows onto nonnegativity and upper-bound constraints, with corresponding primal and dual updates.The bound projections are triggered when a current arc flow is below zero or above its upper limit.Source: S;- L Yij · pp. 382–383

Algorithm 12.4.1 Row-Action Algorithm for Quadratic Transportation Problems

• For quadratic transportation, the row-action algorithm iterates over origin nodes, destination nodes, and flow bounds, with relaxation parameters 0 < E ≤ A ≤ 2 − E ensuring asymptotic convergence.The algorithm initializes primal and dual vectors, then repeats the three update steps.Source: S;- L Yij · p. 383, Algorithm 12.4.1 Row-Action Algorithm… · pp. 383–384

12.4.2 Extensions to generalized networks

• For generalized networks, quadratic objectives yield closed-form projection parameters, whereas generalized nonlinear systems may require iterative solution, such as Newton’s method for twice continuously differentiable objectives.Closed-form solutions are also possible for some nonlinear objectives on pure networks because network constraint coefficients are 0 or 1.Source: §12.4.2 · p. 385
• A one-secant-step approximation gives a closed-form entropy projection parameter and preserves asymptotic convergence, avoiding an exact nonlinear solve at the cost of potentially more iterations.The passage notes that the efficient iterations can offset the increased iteration count in total solution time and enable SIMD parallelism.Source: §12.4.2 · pp. 385, 387, 388

Algorithm 12.4.2 Row-action Algorithm for Entropy Optimization Generalized Transportation Problems

• The generalized transportation entropy algorithm uses exponentiated dual prices and relaxation parameters 0 < E ≤ λ ≤ 1; its positive multiplicative updates need no nonnegative-orthant projection.Setting multipliers to one gives the pure entropy version, which on pure networks is identical to RAS except for the flow-bound step.Source: §12.4.2 · p. 388, Algorithm 12.4.2 Row-action Algorithm… · p. 390

Decompositions for Parallel Computing

• The generalized-network algorithm calculations retain the pure-network algorithm’s structure, enabling parallelization across network updates and arcs.The text states that the algorithm can use concurrent calculations and updates, and refers to Chapter 14 for implementation data structures.Source: Decompositions for Parallel Computing · pp. 384, 385, 390

12.4.3 Row-action algorithm for quadratic multicommodity transportation problems

• For quadratic multicommodity transportation, the algorithm repeats the single-commodity transportation updates and projects commodity flows onto the shared GUB constraints on each arc.The projection handles aggregate flow below zero or above the arc’s upper bound.Source: §12.4.3 · pp. 390–392

Algorithm 12.4.3 Row-action Algorithm for Quadratic Multicommodity Transportation Problems

• Algorithm 12.4.3 initializes the origin, destination, and flow-bound dual prices at zero, then begins iterative solution of each commodity’s problem.Source: Algorithm 12.4.3 Row-action Algorithm… · p. 392

Decompositions for Parallel Computing

• The row-action algorithm parallelizes updates across origin nodes, destination nodes, arcs, and commodities, with processor counts depending on the update step.The passage specifies K m_o processors for Step 1.1 projection parameters and dual variables, and K n processors for flow updates.Source: Decompositions for Parallel Computing · pp. 393–394

12.5 A Model Decomposition Algorithm for Multicommodity Network Flow Problems

• The model decomposition removes complicating GUB constraints from the feasible set and penalizes their violations using a linear-quadratic penalty function.It separates solutions satisfying flow conservation and bounds from those satisfying the complicating constraints.Source: §12.5 · p. 394, §12.5.1 · p. 394

Algorithm 12.5.1 The LQP Algorithm for Multicommodity Network Flow Problems

• The LQP algorithm repeatedly solves a smoothed penalty problem, increasing the penalty parameter and decreasing approximation error until termination criteria are met.Independent network-flow solutions for each commodity can provide a starting point; if they satisfy the GUB constraints, the original problem is already solved.Source: §12.5.1 · pp. 394–395, Algorithm 12.5.1 The LQP Algorithm for… · p. 395

Decompositions for Parallel Computing

• Combining the LQP algorithm with simplicial decomposition yields parallel linear network subproblems and a small, parallelizable coordinating master program.The procedure parallelizes well when there are many commodity blocks or when each block is large.Source: Decompositions for Parallel Computing · p. 396

12.6 Notes and References

• The notes attribute the LQP method to Zenios, Pinar, and Dembo (1994), its multicommodity application to Pinar and Zenios (1992), and a comparative study to Pinar and Zenios (1993).They also note massively parallel data structures developed by Pinar and Zenios (1994a).Source: §12.6 · p. 398

Planning Under Uncertainty

• Because uncertain, noisy, incomplete, or erroneous data can materially affect real systems, the chapter develops stochastic linear programming and robust optimization models for planning under uncertainty.It cautions that mean-value models can yield impractical solutions when the mean is unrealizable, while worst-case models can be overly conservative.Source: Planning Under Uncertainty · p. 400

13.1 Preliminaries

• The preliminaries define probability spaces, discrete scenarios, expectation, variance, and conditional expectation given a finite partition.The conditional expectation is defined on each partition event with positive probability.Source: §13.1 · pp. 401–402, E[Q I A]~ P(~i) L, Q(w)dP(w) · p. 402

13.2 The Newsboy Problem

• The newsboy example chooses newspaper purchases before random demand to maximize expected profit, assuming fixed recourse and ignoring higher moments and alternative reactions.Unsold papers are sold as scrap at a loss, while demand determines the number sold at a profit.Source: §13.2 · pp. 402–403

13.3 Stochastic Programming Problems

• Stochastic programs minimize an expected objective subject to expected constraints, but multidimensional integrals can be difficult to evaluate, motivating approximation and structured decomposition methods.The integrands may be nondifferentiable, nonconvex, or discontinuous; many models can be formulated as large-scale programs with specially structured constraint matrices.Source: §13.3 · pp. 403–404

13.3.1 Anticipative models

• Anticipative models choose decisions before future observations and can impose chance constraints requiring a specified reliability level for uncertain constraints.The threshold depends on the application and the cost of constraint violations.Source: §13.3.1 · p. 404

13.3.2 Adaptive models

• Adaptive models let decisions depend on partial information through conditional expectations, with no information reducing to an anticipative model and complete information yielding a distribution model.The most interesting cases involve partial information becoming available after some decisions have been made.Source: §13.3.2 · pp. 404–405

13.3.3 Recourse models

• Two-stage recourse combines first-stage decisions made before uncertainty is observed with adaptive second-stage decisions, minimizing first-stage cost plus expected recourse cost.Two-stage formulations can yield large-scale deterministic nonlinear programs with special constraint-matrix structure suitable for decomposition and parallel computation.Source: §13.3.3 · pp. 406–407

Deterministic Equivalent Formulation

• With a finite discrete scenario set, each realization receives its own second-stage decision, yielding a large-scale deterministic equivalent program.The resulting constraint matrix has block-angular structure.Source: Deterministic Equivalent Formulation · pp. 407–408, N Minimize f(x) + I> 8 q(y 8 ,w8 ) · p. 408

N Minimize f(x) + I> 8 q(y 8 ,w8 )

• Replicating first-stage variables by scenario makes the system block-diagonal, while equality restrictions enforce nonanticipativity; multistage recourse nests decisions across information stages.For discrete finite distributions, the multistage model can also be formulated as a large-scale deterministic equivalent program.Source: N Minimize f(x) + I> 8 q(y 8 ,w8 ) · pp. 408–410

13.4 Robust Optimization Problems

• Robust optimization separates noise-free structural decisions from controls adjusted to noisy data; in some cases it matches two-stage recourse but permits additional flexibility.The structural variables depend only on noise-free constraints, whereas control variables can be adjusted after noisy inputs are known.Source: §13.4 · p. 410
• Control variables adjust production in response to disruptions or changes in demand or yield, while design variables set system structure and production-module size.The passage says optimal values depend on both the realization of uncertain parameters and the optimal design-variable values.Source: §13.4 · p. 410
• The stated model minimizes (c, x) + (d, y).The passage introduces this as the structure of the optimization model of interest.Source: §13.4 · p. 410

X E JR~O, (13.25) y E lR~', (13.26)

• The robust model balances solution robustness with feasibility robustness through a scenario-dependent objective and penalties for control-constraint violations.The objective can reflect preferences over mean value and variability, while the feasibility penalty accounts for violations that may remain unavoidable across scenarios.Source: X E JR~O, (13.25) y E lR~', (13.26) · pp. 411–412

13.5 Applications

• The applications overview covers robust diet optimization, planning under uncertainty in production and inventory, and matrix balancing, with financial planning treated in Section 13.6.Source: §13.5 · p. 413

13.5.1 Robust optimization for the diet problem

• With an acceptable ±2E vitamin C error, robust diets need no remediation across apple scenarios, while the linear-programming diet exceeds that margin in seven scenarios.The robust diet is somewhat more expensive, and its cost rises as nutritional robustness increases.Source: §13.5.1 · p. 414, Figure 13.2 · p. 416, Figure 13.1

13.5.2 Robust optimization for planning capacity expansion

• Uncertain demand and other operating disruptions make capacity expansion difficult, while worst-case capacity planning can be expensive and leave facilities underutilized.Managers must weigh the marginal profit from peak-period sales against the cost of maintaining excess capacity.Source: §13.5.2 · pp. 414–415

The Multiperiod Stochastic Program

• The multiperiod capacity model makes plant configuration decisions first and production decisions later for each demand scenario and time period.It models one product across multiple sites, allows extension to multiple products, and maximizes direct and diverted-demand profits less retooling costs, subject to demand and capacity constraints.Source: The Multiperiod Stochastic Program · pp. 416, 418
• Variance-based robustness trades expected profit against risk, while scenario profit floors limit downside risk and yield a linear program.The variance objective becomes a nonlinear program with continuous and integer variables and is unsuitable when profit distributions are asymmetric; raising the profit floor guarantees higher scenario profits at the cost of expected profit.Source: The Multiperiod Stochastic Program · pp. 419–420

13.5.3 Robust optimization for matrix balancing

• The matrix-balancing model addresses inconsistent observations by penalizing a norm of constraint infeasibilities in the objective.It is presented as an alternative to scaling the observations or using interval constraints.Source: §13.5.3 · pp. 420–421

Row-action Algorithm for Robust Optimization of Matrix Balancing

• Because the model’s quadratic and entropy objective terms are Bregman functions with strong zone consistency, a row-action algorithm can solve it.A single Newton step gives a closed-form approximation to the row update and preserves asymptotic convergence.Source: Row-action Algorithm for Robust… · pp. 421–422

Algorithm 13.5.1 Row-Action Algorithm for Robust Matrix Balancing Optimization Model

• Algorithm 13.5.1 initializes the variables, then iteratively updates rows and columns of the matrix constraints.After each pass over both sets of constraints, it increments the iteration count and repeats.Source: Algorithm 13.5.1 Row-Action Algorithm… · p. 423

13.6 Stochastic Programming for Portfolio Management

• A three-stage stochastic program models fixed-income portfolio decisions at t0, t1, and t2, with extension to more stages described as straightforward.Decisions anticipate future scenarios and adapt as market information becomes available; the model maximizes expected utility of terminal wealth.Source: §13.6 · pp. 424–425

13.6.1 Notation

• The portfolio notation defines scenario-dependent prices, cashflow factors, and nonnegative buy, sell, hold, and riskless-asset investment decisions across stages.Amortization factors account for changes such as calls and prepayments, while cash accrual factors capture payments and accrued interest.Source: §13.6.1 · p. 426

13.6.2 Model formulation

• The portfolio model tracks securities sold and held and riskless-asset investment, enforces staged cashflow and inventory balances, and maximizes expected utility of marked-to-market terminal wealth.The model’s objective depends on terminal wealth across scenario paths.Source: §13.6.1 · p. 427, §13.6.2 · pp. 427–429

13.7 Stochastic Network Models

• In the United States, ground delays averaged 2,000 hours per day in 1986, and a West German study estimated avoidable 1990 ground-holding delay costs at $1.5 billion.The passage equates the daily delay total to grounding 250 airplanes, a carrier the size of Delta Airlines.Source: §13.7 · pp. 430–431
• At Logan airport, stochastic-programming dynamic models achieved substantial total-delay reductions compared with commonly used static models.The application used data from Logan airport in Boston.Source: §13.7 · p. 431
• The stochastic network problem can be split into scenario-specific deterministic network-flow problems, coupled by non-anticipativity constraints on first-stage decisions.Without the coupling constraints, the scenario problems decompose independently.Source: §13.7 · pp. 431–432

13.7.2 Component-wise representation of the stochastic network problem

• The component-wise formulation represents first-stage and second-stage decisions as network-arc flows, with scenario-specific flow balances, bounds, and equality constraints linking replicated first-stage flows.The network partitions nodes and arcs into first-stage and second-stage sets.Source: §13.7.2 · p. 433, L Psfij(xfj) + L Psqfj(Yfj) (13.80)… · p. 434

13.8 Iterative Algorithm for Stochastic Network Optimization

• The proposed algorithm specializes a general row-action method to stochastic network problems with quadratic objectives and iterates over flow-conservation, bound, and non-anticipativity constraints.Its projection steps address flow-conservation constraints at network nodes.Source: §13.8 · pp. 434–435
• The algorithm updates primal variables and node dual variables for network-flow constraints, then applies steps for simple bounds and non-anticipativity constraints.Its scenario-subproblem calculations are repeated in Step 1, followed by the non-anticipativity step and another iteration.Source: §13.8 · p. 436, Algorithm 13.8.1 Row-Action Algorithm… · p. 441

Projection on Simple Bound Constraints

• For simple bounds, the projection sets a variable below zero to zero or a variable above its upper bound to that bound, and updates the corresponding dual price.When the variable is within its bounds, the dual price is set to zero.Source: Projection on Simple Bound Constraints · pp. 436–437

Projections on Non-anticipativity Constraints

• Repeated almost-cyclic projections onto non-anticipativity constraints converge within tolerance, and their limit can be calculated analytically rather than iteratively.For two replicated variables, an individual projection gives their weighted average; the derivation for the limit uses a system with a unique solution.Source: Projections on Non-anticipativity… · pp. 437–439

The Row-action Algorithm for Quadratic Stochastic Networks

• The complete quadratic stochastic network algorithm combines scenario-subproblem updates, simple-bound projections, and projections on first-stage arcs to enforce non-anticipativity.It returns to the scenario-subproblem step after each iteration.Source: The Row-action Algorithm for Quadratic… · p. 440, Algorithm 13.8.1 Row-Action Algorithm… · pp. 441–442

Decompositions for Parallel Computing

• The calculations for separate scenarios can run concurrently, while graph coloring permits simultaneous node updates without shared arcs; a Jacobi variant is another option.Bound projections can run concurrently across arcs and scenarios, and non-anticipativity calculations across first-stage arcs can also run concurrently.Source: Decompositions for Parallel Computing · p. 442

13.9 Notes and References

• The notes trace the network row-action method to Nielsen and Zenios (1993a, 1996a) and describe graph coloring and Jacobi approaches to parallelism.They report that graph coloring converges faster in practice, while Jacobi permits more processors and may be substantially faster on large-scale problems with enough processors (Zenios and Lasken, 1988a, 1988b).Source: §13.9 · p. 445

Decompositions for Parallel Computing

• Parallel implementation partitions independent algorithm operations or data into tasks; scheduling assigns tasks to processors, while synchronization orders execution and information exchange.The chapter focuses on task partitioning, assuming enough processors to assign one per task; scheduling and synchronization depend on the architecture and are not addressed.Source: Decompositions for Parallel Computing · pp. 446–447

14.1 Vector-Random Access Machine (V-RAM)

• The V-RAM extends a RAM with vector memory and a parallel vector processor, and its parallelism scales with problem data or the machine’s vector length.Vector lengths may differ, and when physical processors are insufficient, one processor can emulate multiple virtual processors.Source: §14.1 · pp. 447–448

14.1.1 Parallel prefix operations

• V-RAM scans compute associative prefixes, spreads replicate an associative aggregate across a vector, and segmented scans restart the operation at segment boundaries.Segment boundaries are encoded by bits set to 1 at each segment’s start; the chapter uses scans and spreads as vector instructions.Source: §14.1.1 · pp. 448–449

14.2.1 Mapping a dense matrix

• Dense matrices can be mapped by checkerboard blocks across processors or by striped blocks of complete rows or columns, with striped mapping using at most n virtual processors for an n x n matrix.With p processors, checkerboard blocks have dimensions n/√p x n/√p.Source: §14.2.1 · p. 449

The V-RAM Model of a Dense Matrix

• A dense-matrix V-RAM layout depends on the operation: row and column sums use separate row-wise or column-wise vectors, mapped via checkerboard or corresponding striped partitions.The row-sum model uses n vectors of length n, each holding one row; the column-sum model similarly stores one column per vector.Source: The V-RAM Model of a Dense Matrix · p. 450

14.2.2 Mapping a sparse matrix

• Sparse matrices use compressed row-wise storage of nonzeros and can be striped by rows, but row blocks should balance nonzero counts to avoid processor load imbalance.Randomly assigning rows to blocks is offered as a simple way to mix dense and sparse rows and yield approximately equal nonzero counts.Source: §14.2.2 · pp. 450–452
• A sparse-matrix V-RAM stores row entries in segmented vectors, enabling row sums with segmented add-scan; simultaneous row and column sums use concatenated row- and column-format segments.The simultaneous-sum model uses 2q virtual processors in 2n segments, with one segment for each row and column.Source: §14.2.2 · pp. 452–453

14.3.1 Data parallel computing with RAS

• Data-parallel RAS computes dense row and column scaling in parallel, while sparse RAS uses compressed row and column copies and communicates updated entries between formats before column processing.The dense implementation uses add-spread operations for row and column partial sums; the sparse implementation uses segmented operations on row data and transfers updated values to the column-wise representation.Source: Solving Dense Problems · p. 454, Solving Sparse Problems · p. 455

Solving Dense Problems

• For dense problems, Step 2 sums each column and distributes its total, then uses the locally available row total to estimate a scaling factor and scale entries.Each VP divides the row total by the computed partial sum before applying the factor to the matrix entries.Source: Solving Sparse Problems · p. 455

14.3.2 Control parallel computing with RAS

• RAS scales rows before columns; independent row updates can run in parallel, with task scheduling trading load balance against task-initiation overhead.Alternatively, processors can receive row blocks sized to have approximately equal numbers of nonzeros.Source: §14.3.2 · p. 456

Algorithm 14.3.1 Simultaneous RAS Algorithm

• Simultaneous RAS updates row and column scaling factors and matrix entries in parallel; as a block-iterative MART special case, it needs no separate convergence analysis.Both row and column sums can be computed simultaneously, allowing simultaneous RAS to use more processors than RAS.Source: §14.3.2 · p. 456, Algorithm 14.3.1 Simultaneous RAS… · pp. 456–457

14.4 Parallel Computing for Image Reconstruction

• Block-MART supports three control-parallel image-reconstruction schemes, each a special case of the algorithm, so separate mathematical convergence analysis is unnecessary.The schemes illustrate the flexibility of block-iterative algorithms for parallel computation.Source: §14.4 · pp. 457–458

J J ~ t t)' iEft(v)

• The basic image update accumulates ray contributions in a block correction array, then updates affected pixels with vector operations.BLAS routines support implementation on vector and parallel architectures.Source: J J ~ t t)' iEft(v) · pp. 458–459

14.4.1 Parallelism within a block

• Within a block, ray geometry and measurements and correction terms can be computed concurrently, but pixel updates need synchronization when rays share pixels.Selecting rays for a block also requires synchronization, though the passage describes that step as not particularly time-consuming.Source: §14.4.1 · p. 460

14.4.2 Parallelism with independent blocks

• Independent blocks group rays from the same view that intersect no common pixels, avoiding update conflicts within that view; synchronization is needed before processing a new view.Blocks can therefore be processed in parallel without synchronizing pixel updates until the algorithm changes views.Source: §14.4.2 · pp. 460–461

14.4.3 Parallelism between views

• The between-view simultaneous block-MART algorithm is mathematically equivalent to a fully simultaneous algorithm, while offering a different parallel implementation based on per-view corrections.Because rays from each view update pixel intensities concurrently, temporary correction storage is aggregated before the image vector is updated.Source: §14.4.3 · p. 462

14.5 Parallel Computing for Network-structured Problems; 14.5.4 Solving network-structured problems

• V-RAM models for network-structured problems generalize the sparse transshipment model, with scenario-specific vectors for stochastic programs and extra vectors for shared constraints.Scenarios can be processed concurrently, but first-stage split variables must converge to the same value; multicommodity flows must respect shared arc bounds.Source: §14.5 · p. 462, §14.5.4 · p. 467

14.5.1 Solving dense transportation problems

• The dense transportation V-RAM model uses at most m^2 K virtual processors by operating on origin-node vectors and then destination-node vectors.The arc-related vectors are symmetric between origin and destination nodes, allowing those vectors to be shared.Source: §14.5.1 · p. 464

14.5.2 Solving sparse transportation problems; 14.5.3 Solving sparse transshipment graphs

• Sparse transportation uses degree-length vectors and segmented operations, while sparse transshipment stores each arc twice and needs 2n virtual processors.In sparse transportation, origin- and destination-side copies of arc variables must be reconciled when switching between node iterations.Source: §14.5.2 · pp. 464, 466, §14.5.3 · p. 467

14.6 Parallel Computing with Interior Point Algorithms; 14.6.1 The communication schemes on a hypercube

• The stochastic-programming interior-point procedure parallelizes independent submatrix operations, with communication for forming shared quantities, broadcasting, and gathering partial results.The target is a distributed-memory MIMD machine whose hypercube has 2^d processor nodes for dimension d.Source: §14.6 · pp. 467–468, §14.6.1 · p. 468

14.6.1 The communication schemes on a hypercube

• STB broadcasts, STG gathers, and ADE accumulates distributed data using nearest-neighbor exchanges; each takes the minimum d communication steps on a d-cube.The text notes that minimizing communication steps matters because communication costs are typically high relative to computation.Source: §14.6.1 · pp. 468–470

14.6.2 The parallel implementation on a hypercube

• The hypercube implementation distributes scenarios across processors, forms global sums at node 0, performs small dense-matrix work there, and executes the remaining steps largely in parallel.The final accumulation of the full correction vector at the master processor can be omitted if the next iteration needs a distributed vector.Source: §14.6.2 · pp. 470–471

14.6.3 An alternative parallel implementation

• On the CM-5, add-spread forms global sums on every processor, enabling redundant parallel execution of Steps 2c–2e and eliminating the broadcast of q1 before Step 3.The hypercube implementation instead forms G1 and the related sum at node 0, performs those steps serially there, then broadcasts q1.Source: §14.6.3 · p. 471

Procedure 14.6.1 Parallel matrix factorization for dual step direction calculation on the Connection Machine CM-5

• Procedure 14.6.1 uses add-spread to form G1 and p1 on all processors, then redundantly solves the small systems and computes direction components across processors.Processors hold scenario-specific data and solve the corresponding systems in parallel.Source: Procedure 14.6.1 Parallel matrix… · pp. 471–472

14.7 Notes and References

• The notes cite prior work on parallel interior-point implementations, hypercube communication, performance modeling, and CM-5 fat-tree communication (Jessup, Yang, and Zenios (1994a, 1994b); Saad and Schultz (1989); Leiserson et al. (1992)).Source: §14.7 · p. 473

Numerical Investigations

• The chapter’s experiments assess parallel-algorithm performance by comparing parallel implementations with the same or competing algorithms on contemporary workstations or mainframes.The experiments are described as typical and illustrative; the discussion does not evaluate machine cost or ease of parallel programming.Source: Numerical Investigations · pp. 474–475

15.1 Reporting Computational Experiments on Parallel Machines

• Reporting guidance distinguishes machine-efficiency tests, where time may suffice, from code-efficiency tests, where work includes execution, input/output, and memory-access time.The chapter follows the guidance only generally, aiming to provide an overview of parallel-algorithm performance.Source: §15.1 · p. 476

on Parallel Machines

• The experiments recommend documenting algorithm and code, computing and communication environments, testing methodology, and performance metrics, while reporting results for difficult, typical problems with relevant factors and variability.The section also reports that a linear-quadratic penalty model-decomposition algorithm was implemented with control parallelism and tested on large-scale military-logistics problems.Source: on Parallel Machines · pp. 476–477, §15.4.3 · p. 488

15.2 Matrix Balancing

• Large-scale problems took several hours on current mainframes and several minutes on vector supercomputers; on the same equality-constrained test problems, RRAS was roughly twice as slow as RAS, but both finished within a few seconds.Source: §15.2 · p. 478, §15.2.1 · p. 479

15.2.2 Control parallel implementations

• Control-parallel RAS on an Alliant FX/8 achieved significant speedups that increased almost linearly with processor count, while simultaneous RAS appeared superlinear because it was more efficient even on one processor.Comparing parallel simultaneous RAS with its own uniprocessor implementation would yield the same speedup curve as RAS.Source: §15.2.2 · p. 480, Figure 15.1 · p. 481

15.3 Image Reconstruction

• A vectorized CRAY implementation improved block-MART performance by about tenfold, while the CRAY achieved a 100-fold solution-time reduction over a SUN-3 workstation.The parallel schemes within a block and across independent blocks each achieved more than 75 percent efficiency; the comparison used 45 blocks per view.Source: §15.3 · pp. 481–482, Table 15.7 · p. 482, Table 15.8 · p. 483
• For SHEPP512, 30 blocks per view was optimal: intermediate block sizes balanced large tasks and low synchronization against enough tasks for processor load balancing.The optimal block size depends on the problem and parallel machine; neither one equation per block nor all equations of a view in one block was best in this experiment.Source: §15.3 · p. 483, Figure 15.2 · p. 484

15.4.1 Row-action algorithm for transportation problems

• On the CM-2, dense row-action code ran at 1.5 GFLOPS and sparse code at 275 MFLOPS, yet sparse code outperformed dense code on the very sparse test problems.The dense code ran at a higher FLOPS rate, but the sparse implementation solved those problems faster.Source: §15.4.1 · pp. 484–485, Table 1 · p. 485
• The CM-2 row-action algorithm was up to two orders of magnitude faster than state-of-the-art serial algorithms on an IBM 3090 for extremely large, dense quadratic problems.The serial methods were an equilibration algorithm and an exact Newton method.Source: Table 15.12 · p. 486, §15.4.1 · p. 486

15.4.2 Row-action algorithm for multicommodity transportation problems

• The CM-2 data-parallel multicommodity row-action algorithm solved problems with millions of variables within 4–10 minutes; on an Alliant FX/8, efficiency exceeded 80 percent for sufficiently many commodities.On the Alliant, speedup did not rise monotonically with commodity count, and good load balancing required commodity counts to be multiples of processor count.Source: Table 15.13 · p. 487, §15.4.2 · p. 487, Figure 15.3 · p. 488

15.4.3 Linear-quadratic penalty (LQP) algorithm for multicommodity network flow problems

• Using eight CRAY processors, LQP achieved speedups of up to four, but serial master work and uneven subproblem difficulty limited efficiency to below 50%.The variation in network-simplex pivots across commodity subproblems contributed to load imbalance.Source: §15.4.3 · p. 489, Table 15.15 · p. 490, Figure 15.4 · p. 491
• LQP outperformed the general-purpose linear programming code OB1 on a single processor, and its advantage increased with parallelism.Both codes ran on a CRAY Y-MP.Source: §15.4.3 · p. 490, Figure 15.5 · p. 491

15.5.2 Row-action algorithm for nonlinear stochastic networks

• The row-action algorithm solved large stochastic network problems in under half a minute, with solution time not directly dependent on scenario count or deterministic-equivalent size.Larger problems could use more processors, while iteration counts did not necessarily increase with problem size.Source: Table 15.20 · p. 496, §15.5.2 · p. 497
• The row-action implementation was scalable and substantially outperformed the general-purpose optimizer OBN, with only marginal time increases as problem size grew.On a CM-2 with 32K processing elements, it solved four times as many scenarios in almost the time the 8K machine needed for 512 scenarios: 108.2 seconds.Source: §15.5.2 · pp. 497–498, Figure 15.7 · p. 497

15.6.1 Solving linear network problems

• On extremely large linear network problems, the parallel PMD algorithms ran successfully when the serial network-simplex code could not solve them.For smaller problems, the serial code on a CRAY Y-MP was faster; the performance gap narrowed as problem size increased.Source: §15.6.1 · p. 501, Figure 15.10 · p. 502

15.6.2 Solving linear stochastic network problems

• For smaller stochastic network tests, OSL was 1.3 to 3.7 times faster than QPP, but QPP’s relative performance improved with problem size and its Random3 objective was within 0.02% of OSL’s.For a similarly generated three-stage problem with four times as many scenarios, QPP took 4.2 times longer, whereas OSL took 19.3 times longer.Source: §15.6.2 · pp. 501, 503, Table 15.26 · p. 504
• QPP achieved high accuracy on extremely large stochastic network problems, although execution times were substantial.The passage refers to the large test problems summarized in Tables 15.27 and 15.28.Source: §15.6.2 · p. 503, Table 15.27 · p. 505

15.7.2 Connection Machine CM-2

• 3 GFLOPS was the peak optimization performance reported for a fully configured 64K-PE CM-2; the QPP test problems used 32K PEs and reported solution times in CM hours.The CM-2 is a fine-grain SIMD system with 16 PEs per chip, a 12-dimensional hypercube interconnect, and local memory; floating-point accelerators serve clusters of 32 PEs.Source: Table 15.28 · p. 505, §15.7.2 · p. 506

15.7 .3 Connection Machine CM-5

• The largest CM-5 system used in the experiments had 64 processors, while the largest configuration installed at the time had 1,000.The CM-5 is a distributed-memory MIMD machine with SPARC-based nodes, vector functional units, and a fat-tree communication network.Source: §15.7 · p. 506

15.8 Notes and References

• Jessup, Yang, and Zenios report parallel factorization results for structured stochastic-programming matrices on the Intel iPSC/860 and CM-5.The notes also cite a hypercube communication-model analysis and later work on parallel matrix factorization with interior-point methods.Source: §15.8 · p. 508
• Nielsen and Zenios report PMD results for linear network problems and stochastic network problems, extending their experiments to multistage cases.The notes cite Nielsen and Zenios (1993b, 1993c) for linear problems, (1993d) for two-stage problems, and (1996a) for multistage problems.Source: §15.8 · pp. 508–509