Iterated Penalty Methods in Optimization

• Iterated Penalty Method is a family of algorithms that reformulate constrained problems into penalized formulations solved via repeated iterative updates.
• These methods vary by using fixed or adaptive penalties, dual updates, and active-set strategies to handle nonsmoothness and improve convergence.
• Applications include discrete Hamilton–Jacobi–Bellman systems, incompressible flow, and inverse problems, offering robust convergence with geometric rates.

Iterated penalty method denotes a family of algorithms in which an original constrained, complementarity, or saddle-point problem is replaced by a penalized formulation and then solved through repeated updates. In the cited literature, the term covers several distinct but related constructions: penalized discrete Hamilton–Jacobi–Bellman systems solved by monotone Newton-like iterations, regularized exact penalties for inequality-constrained optimization, double-loop penalty–multiplier schemes, adaptive penalty limits of active-set methods, Uzawa-type divergence-penalty solvers for incompressible flow, and nonstationary Tikhonov iterations with uniformly convex penalties (Witte et al., 2010, Ito et al., 2012, Shi et al., 2017, Boon et al., 2022, Eickmann et al., 22 Sep 2025, Jin et al., 2014). This diversity is substantive rather than terminological: in some works the iteration acts on a single penalized problem with fixed penalty parameters, whereas in others the penalty parameter itself is updated across outer iterations.

1. Terminological scope and recurring designs

Across the literature, “iterated penalty method” does not designate a single canonical template. The main variants differ in where the iteration occurs: on the state variable for one penalized problem, on the penalty parameter in an outer continuation loop, on multiplier-like corrections, or on active sets induced by the current residual.

Interpretation	Defining feature	Representative source
Regularized exact penalty solve	Fixed $\beta$, fixed $\epsilon$, iterative primal update	(Ito et al., 2012)
Double-loop penalty method	Inner approximate augmented solve, outer dual/penalty update	(Shi et al., 2017)
Adaptive penalty / active-set limit	Penalty recomputed from current residual and driven to zero	(Boon et al., 2022)
Penalty-based but not classical outer-loop	Single penalized reformulation with fixed design parameters	(Nazari et al., 2022)
Single-loop penalty dual-primal ALM	Fixed $\beta$, dual-first then primal update	(Liu et al., 2023)

A central distinction is between classical outer-loop penalty continuation and iterative solution of one penalized problem. Ito and Takeuchi explicitly state that their method is not a classical outer iterated penalty method $(\beta_k\uparrow\infty)$, but rather an iterative solver for a single exact penalty problem after regularization by $\epsilon$ (Ito et al., 2012). By contrast, penalty dual decomposition is a double-loop method whose inner iterations solve a nonconvex nonsmooth augmented Lagrangian approximately and whose outer iteration updates dual variables and/or a penalty parameter (Shi et al., 2017).

This distinction matters algorithmically. In the first case, the main challenge is nonsmoothness and the design of a globally convergent inner iteration. In the second, the main challenge is coordinating feasibility improvement, penalty strengthening, and multiplier correction. A further qualification is provided by decentralized bilevel optimization: the penalty-function-based method of Li et al. is penalty-based, but it does not implement a standard iterated penalty method with a sequence $\mu_t$ or $\alpha_t$; instead it solves a single penalized reformulation with essentially fixed $\alpha$ and $\gamma$ (Nazari et al., 2022). A similar fixed-parameter, single-loop perspective appears in the penalty dual-primal augmented Lagrangian method, where the quadratic penalty/proximal terms are iterated but not scheduled to infinity (Liu et al., 2023).

2. Penalized discrete HJB systems

One of the most explicit algorithmic realizations of an iterated penalty method appears in the numerical solution of discrete Hamilton–Jacobi–Bellman equations. For a finite control set $\mathbb S$, Forsyth and Labahn consider the nonlinear algebraic system

$\epsilon$0

obtained from monotone, stable finite-difference discretization, typically with fully implicit time stepping (Witte et al., 2010). Under matrix assumptions—non-positive off-diagonal entries, positive diagonal entries, non-negative row sums, at least one positive row sum, and common location of the positive row sum—the matrices $\epsilon$1 are $\epsilon$2-matrices, hence invertible with $\epsilon$3. The discrete HJB system then has a unique solution.

The penalized approximation fixes a distinguished control $\epsilon$4 and $\epsilon$5, and solves

$\epsilon$6

The penalty term measures violations of the inequalities $\epsilon$7. The paper proves existence and uniqueness of the penalized solution, uniform boundedness independent of $\epsilon$8, the residual estimate

$\epsilon$9

and the first-order penalty accuracy

$\beta$0

where $\beta$1 is the exact discrete HJB solution (Witte et al., 2010). The inner solver is a Newton-like active-set iteration

$\beta$2

equivalently

$\beta$3

Each linear system involves an $\beta$4-matrix, the iterates are monotone, and finite termination follows because only finitely many active-row compositions are possible. In the reported finance examples, policy iteration often converges in $\beta$5–$\beta$6 iterations, the penalty method usually takes $\beta$7–$\beta$8 iterations, and penalty runtime is roughly $\beta$9–$(\beta_k\uparrow\infty)$0 times longer (Witte et al., 2010).

A closely related development extends the penalty framework to continuous control and HJB obstacle problems. Here the discrete HJB problem is

$(\beta_k\uparrow\infty)$1

and the penalized equation is

$(\beta_k\uparrow\infty)$2

with diagonal penalty operator

$(\beta_k\uparrow\infty)$3

where each $(\beta_k\uparrow\infty)$4 is non-decreasing, vanishes on $(\beta_k\uparrow\infty)$5, and is positive on $(\beta_k\uparrow\infty)$6 (Witte et al., 2011). For the canonical choice $(\beta_k\uparrow\infty)$7, the residual satisfies

$(\beta_k\uparrow\infty)$8

and under the lower derivative bound

$(\beta_k\uparrow\infty)$9

the solution error satisfies $\epsilon$0 (Witte et al., 2011). The corresponding obstacle formulation penalizes $\epsilon$1 in

$\epsilon$2

For smooth penalties, Newton with line search is globally convergent and locally quadratically convergent near the penalized solution. For the nonsmooth max-penalty, Newton-like active-set updates solve linear systems such as

$\epsilon$3

and the numerical tests report first-order slopes $\epsilon$4 and $\epsilon$5 for the HJB and obstacle examples, respectively (Witte et al., 2011).

3. Exact penalties and regularized primal iterations in optimization

In finite-dimensional constrained optimization, the iterated penalty idea frequently appears as the repeated solution of a single exact penalty problem rather than as $\epsilon$6. For unilateral constraints $\epsilon$7, Ito and Takeuchi consider the exact penalty functional

$\epsilon$8

equivalently

$\epsilon$9

(Ito et al., 2012). Exactness is linked to the KKT system

$\mu_t$0

and if $\mu_t$1 solves this system, then $\mu_t$2 solves the exact-penalty optimality condition provided

$\mu_t$3

Because $\mu_t$4 is nonsmooth, the paper introduces a regularized penalty

$\mu_t$5

with the piecewise-defined smoothing $\mu_t$6 of $\mu_t$7, and derives the implicit iteration

$\mu_t$8

The method is a single-loop primal iteration with fixed $\mu_t$9 and fixed $\alpha_t$0 during the solve. It yields descent, global convergence under the stated assumptions, and local superlinear convergence once $\alpha_t$1 is close to the solution and $\alpha_t$2, because then $\alpha_t$3 and the method inherits the semismooth Newton rate (Ito et al., 2012).

A matrix-free convex variant arises in exact-penalty subproblems of the form

$\alpha_t$4

with $\alpha_t$5 and product-structured $\alpha_t$6 (Burke et al., 2014). The new Iterative Re-Weighting Algorithm minimizes at each step the quadratic model

$\alpha_t$7

where

$\alpha_t$8

and decreases the relaxation vector $\alpha_t$9 adaptively. The paper proves global convergence under loose assumptions and at most $\alpha$0 iterations to reach $\alpha$1-optimality of the objective function for both IRWA and its ADAL alternative (Burke et al., 2014).

For equality-constrained nonlinear optimization, Fletcher’s smooth exact penalty gives a further reinterpretation. Estrin, Friedlander, Orban, and Saunders implement

$\alpha$2

where

$\alpha$3

and show that the computational kernel is a structured linear system that can be reused for the penalty value, gradient, and Hessian-vector products (Estrin et al., 2019). The penalty is locally exact for finite $\alpha$4 above the threshold

$\alpha$5

and the regularized continuation

$\alpha$6

yields a genuine outer-inner iterative penalty/regularization method with local $\alpha$7-quadratic convergence (Estrin et al., 2019).

4. Augmented, adaptive, and stochastic penalty updates

The most explicit outer-inner penalty architecture in the supplied literature is penalty dual decomposition for linearly constrained, block-structured, nonsmooth, and nonconvex problems of the form

$\alpha$8

The augmented penalty function is

$\alpha$9

At outer iteration $\gamma$0, the method approximately solves the augmented subproblem and then either updates the multiplier,

$\gamma$1

if feasibility is adequate, or decreases $\gamma$2 via

$\gamma$3

if feasibility is not (Shi et al., 2017). This is a penalty-based outer-inner method with inexact inner solves and convergence to KKT solutions.

An adaptive penalty limit of active-set methods is developed for the obstacle problem

$\gamma$4

The complementarity operator

$\gamma$5

satisfies

$\gamma$6

so the primal equation is $\gamma$7 (Boon et al., 2022). The smoothed penalty/Newton step uses

$\gamma$8

and

$\gamma$9

The penalty field $\mathbb S$0 is recomputed from the current residual by solving

$\mathbb S$1

As the residual decreases, $\mathbb S$2, and the derivative weight converges to the active-set indicator: $\mathbb S$3 This establishes that the primal-dual active set method is exactly the $\mathbb S$4 limit of the regularized penalty/Newton method (Boon et al., 2022).

For large collections of linear inequalities,

$\mathbb S$5

a Huber-loss-based single-loop incremental penalty method uses the time-varying penalized objective

$\mathbb S$6

with one-sided Huber penalty

$\mathbb S$7

The stochastic update is

$\mathbb S$8

with conditions

$\mathbb S$9

Under convexity assumptions, the iterates converge to a solution almost surely and in expectation, with weighted-average rates

$\epsilon$00

(Nedich et al., 2023).

These examples also delimit the notion. A penalty-function-based decentralized bilevel algorithm replaces consensus constraints by graph quadratic penalties, for example

$\epsilon$01

and proves $\epsilon$02 as $\epsilon$03, but it does not run a classical outer penalty schedule; $\epsilon$04 and $\epsilon$05 are fixed design parameters in a single alternating decentralized algorithm (Nazari et al., 2022).

5. Divergence-penalty iterations in incompressible flow

For exactly divergence-free finite element discretizations of Stokes flow, the iterated penalty method takes the form of a velocity-space Uzawa iteration. In the Scott–Vogelius setting with

$\epsilon$06

the pressure can be represented as the discrete divergence of an auxiliary velocity-like variable (Eickmann et al., 22 Sep 2025). The inhomogeneous Dirichlet problem uses compatible boundary data $\epsilon$07, and the iteration is initialized with

$\epsilon$08

At step $\epsilon$09,

$\epsilon$10

with $\epsilon$11 on $\epsilon$12, followed by

$\epsilon$13

Equivalently,

$\epsilon$14

The convergence factor is

$\epsilon$15

and the paper proves

$\epsilon$16

together with pressure convergence and asymptotic pressure robustness as $\epsilon$17 (Eickmann et al., 22 Sep 2025). A notable feature specific to this exactly divergence-free setting is monotonic decrease of the divergence norm: $\epsilon$18

For high-order Scott–Vogelius discretizations, the statically condensed iterated penalty method preserves the same geometric convergence mechanism while reducing the size of the global system solved at each iteration. Standard IP requires solving on $\epsilon$19 velocity unknowns per element, whereas SCIP reduces the iterative global system to $\epsilon$20 boundary unknowns by using specially constructed spaces $\epsilon$21 and $\epsilon$22 (Ainsworth et al., 2023). The reduced iteration is

$\epsilon$23

with

$\epsilon$24

After convergence, the interior Stokes corrections are solved once, element by element. The resulting error bound is still geometric: $\epsilon$25 although the theoretical lower bound on $\epsilon$26 is larger than for the standard method (Ainsworth et al., 2023).

For steady incompressible Navier–Stokes, the classical iterated penalty Picard method uses

$\epsilon$27

The associated fixed-point map $\epsilon$28 is Lipschitz continuous and Lipschitz continuously Fréchet differentiable in the norm

$\epsilon$29

which permits the application of Anderson acceleration theory (Rebholz et al., 2021). The accelerated residual estimate has leading factor

$\epsilon$30

so Anderson acceleration improves the linear convergence rate by the gain factor from the underlying least-squares problem. The numerical experiments show that with penalty parameter $\epsilon$31, unaccelerated IPP can fail on challenging cases, whereas AAIPP with moderate or large memory converges effectively (Rebholz et al., 2021).

6. Regularization and inverse-problem interpretations

In inverse problems, iterated penalty methods are closely tied to nonstationary Tikhonov regularization and Bregman geometry. For the linear equation

$\epsilon$32

in Hilbert spaces, Jin and Zou propose the two-step scheme

$\epsilon$33

where $\epsilon$34 is a proper, lower semicontinuous, convex penalty (Jin et al., 2014). The first step involves only the operator $\epsilon$35; the second only the penalty. When

$\epsilon$36

the method reduces to classical nonstationary iterated Tikhonov regularization. With nonsmooth penalties such as

$\epsilon$37

the scheme promotes sparsity or piecewise constancy while remaining computationally split (Jin et al., 2014).

The Banach-space counterpart uses the Bregman distance

$\epsilon$38

and defines the noisy-data iterate by

$\epsilon$39

$\epsilon$40

Under the stated assumptions—uniform convexity of $\epsilon$41, weak closedness and differentiability of $\epsilon$42, and $\epsilon$43—the exact-data iterates converge strongly, with

$\epsilon$44

and the discrepancy principle

$\epsilon$45

yields the regularization property for noisy data (Jin et al., 2014). The paper emphasizes that this gives strong convergence and convergence in Bregman distance, in contrast with earlier results that provided only weak convergence for a subsequence.

A different reformulation appears when the target objective already contains a large quadratic penalty,

$\epsilon$46

Instead of directly minimizing $\epsilon$47, the modified augmented Lagrangian method introduces $\epsilon$48 and rewrites the problem as

$\epsilon$49

Eliminating $\epsilon$50 after applying ALM leads to the modified root function

$\epsilon$51

which is designed for problems where the direct quadratic-penalty objective is badly scaled (Neuenhofen, 2018). This suggests a broader interpretation of iterated penalty methods: they may target the original constrained problem, a penalized approximation of it, or a reformulation of an already penalized objective, depending on which form yields better-conditioned iterative subproblems.

Taken together, these strands show that the iterated penalty method is best understood as an algorithmic family rather than a single procedure. The unifying pattern is the repeated solution of penalized subproblems whose penalty term is used not merely to encode constraints, but also to shape local geometry, recover hidden dual variables, or eliminate difficult spaces such as pressure bases. The principal axes of variation are the role of the penalty parameter, the presence or absence of multiplier updates, the choice between nonsmooth and smoothed penalties, and whether the iteration is driven by active-set identification, Newton linearization, Bregman proximalization, or Uzawa-type residual correction.