--- title: Iterated Penalty Methods in Optimization url: https://www.emergentmind.com/topics/iterated-penalty-method type: topic --- # Iterated Penalty Methods in Optimization 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 [1008.0401] [1210.1301] [1712.04767] [2201.02425] [2509.17899] [1401.4541]. 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 | [1210.1301] | | Double-loop penalty method | Inner approximate augmented solve, outer dual/penalty update | [1712.04767] | | Adaptive penalty / active-set limit | Penalty recomputed from current residual and driven to zero | [2201.02425] | | Penalty-based but not classical outer-loop | Single penalized reformulation with fixed design parameters | [2211.04088] | | Single-loop penalty dual-primal ALM | Fixed \(\beta\), dual-first then primal update | [2305.03922] | 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\) [1210.1301]. 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 [1712.04767]. 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\) [2211.04088]. 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 [2305.03922]. ## 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 \[ \min\{A_sx-b_s:s\in\mathbb S\}=0, \] obtained from monotone, stable finite-difference discretization, typically with fully implicit time stepping [1008.0401]. 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 \(A_s\) are \(M\)-matrices, hence invertible with \(A_s^{-1}\ge 0\). The discrete HJB system then has a unique solution. The penalized approximation fixes a distinguished control \(s_0\in\mathbb S\) and \(\rho>0\), and solves \[ (A_{s_0}x_\rho-b_{s_0})-\rho\sum_{s\in\mathbb S\setminus\{s_0\}}\max\{b_s-A_sx_\rho,0\}=0. \] The penalty term measures violations of the inequalities \(A_sx-b_s\ge 0\). The paper proves existence and uniqueness of the penalized solution, uniform boundedness independent of \(\rho\), the residual estimate \[ \left\|\min\{A_sx_\rho-b_s:s\in\mathbb S\}\right\|_\infty \le \frac{C}{\rho}, \] and the first-order penalty accuracy \[ \|x^*-x_\rho\|_\infty\le \frac{C}{\rho}, \] where \(x^*\) is the exact discrete HJB solution [1008.0401]. The inner solver is a Newton-like active-set iteration \[ J_G(x^n)(x^{n+1}-x^n)=-G(x^n), \] equivalently \[ \left(A_{s_0}+\rho\sum_{s\in\mathbb S\setminus\{s_0\}}A_s^{x^n}\right)x^{n+1} = b_{s_0}+\rho\sum_{s\in\mathbb S\setminus\{s_0\}}b_s^{x^n}. \] Each linear system involves an \(M\)-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 \(1\)–\(2\) iterations, the penalty method usually takes \(3\)–\(4\) iterations, and penalty runtime is roughly \(3\)–\(4\) times longer [1008.0401]. A closely related development extends the penalty framework to **continuous control** and **HJB obstacle problems**. Here the discrete HJB problem is \[ \min_{u\in\mathbf U}\{A_ux-b_u\}=0, \] and the penalized equation is \[ (A_{u_0}x_{\rho} - b_{u_0}) -\rho \max_{u\in\mathbf{U}}\Pi(b_u-A_u\,x_{\rho}) =0, \] with diagonal penalty operator \[ \Pi(y)=\big(\pi_1(y_1),\dots,\pi_N(y_N)\big)^{tr}, \] where each \(\pi_i\) is non-decreasing, vanishes on \((-\infty,0]\), and is positive on \((0,\infty)\) [1105.5954]. For the canonical choice \(\Pi(y)=\max\{y,0\}\), the residual satisfies \[ \big\|\min_{u\in\mathbf U}\{A_u x_\rho-b_u\}\big\|_\infty\le C/\rho, \] and under the lower derivative bound \[ \inf\left\{\frac{\partial \pi_i}{\partial y}(y):y\in(0,\infty),\ i\in\mathcal N\right\}\ge c_{\min}>0, \] the solution error satisfies \(\|x^*-x_\rho\|_\infty\le C/\rho\) [1105.5954]. The corresponding obstacle formulation penalizes \(\tilde A z-\tilde b\ge 0\) in \[ \max_{u\in\mathbf{U}}\{A_u\,z_{\rho}-b_u\}-\rho\Pi\big(\tilde{b}-\tilde{A}z_{\rho}\big) =0. \] 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 \[ \big(A_{u^{\max}(z^n)} + \rho A^{+}(z^n)\big)z^{n+1} =b_{u^{\max}(z^n)} + \rho b^{+}(z^n), \] and the numerical tests report first-order slopes \(0.992\) and \(0.910\) for the HJB and obstacle examples, respectively [1105.5954]. ## 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 \(\beta_k\to\infty\). For unilateral constraints \(Gx-g\le 0\), Ito and Takeuchi consider the exact penalty functional \[ J(x)=F(x)+\beta\,\psi(Gx-g),\qquad \psi(y)=\sum_i\max(0,y_i), \] equivalently \[ J(x)=F(x)+\beta\sum_{i=1}^m \max(0,(Gx-g)_i) \] [1210.1301]. Exactness is linked to the KKT system \[ F^\prime(x)+G^t\mu=0,\qquad \mu=\max(0,\mu+(Gx-g)), \] and if \((\bar x,\bar\mu)\) solves this system, then \(\bar x\) solves the exact-penalty optimality condition provided \[ \beta\ge \max_i |\bar \mu_i|=\|\bar\mu\|_\infty. \] Because \(J\) is nonsmooth, the paper introduces a regularized penalty \[ J_\epsilon(x)=F(x)+\beta \psi_{\epsilon}(Gx-g), \] with the piecewise-defined smoothing \(\phi_\epsilon\) of \(s\mapsto \max(0,s)\), and derives the implicit iteration \[ (\alpha\,P + \beta\,G^t \chi_\epsilon(x^k) G )d^k = -J^\prime_\epsilon(x^k),\qquad x^{k+1}=x^k+d^k. \] The method is a single-loop primal iteration with fixed \(\beta\) and fixed \(\epsilon\) during the solve. It yields descent, global convergence under the stated assumptions, and local superlinear convergence once \(x^k\) is close to the solution and \(Gx^k-g\le \epsilon\), because then \(\chi_k=N_k\) and the method inherits the semismooth Newton rate [1210.1301]. A matrix-free convex variant arises in exact-penalty subproblems of the form \[ \min_{x\in X}\; g^Tx+\frac12 x^THx+\operatorname{dist}(Ax+b\mid C), \] with \(H\succeq 0\) and product-structured \(C\) [1402.1917]. The new Iterative Re-Weighting Algorithm minimizes at each step the quadratic model \[ \hat G_{(x^k,\epsilon^k)}(x) = g^Tx+\frac12 \sum_{i\in \mathcal I_0} w_i(x^k,\epsilon^k)\, \big\|A_ix+b_i-P_{C_i}(A_ix^k+b_i)\big\|_2^2, \] where \[ w_i(x,\epsilon)= \left(\operatorname{dist}_2^2(A_ix+b_i\mid C_i)+\epsilon_i^2\right)^{-1/2}, \] and decreases the relaxation vector \(\epsilon^k\) adaptively. The paper proves global convergence under loose assumptions and at most \(O(1/\varepsilon^2)\) iterations to reach \(\varepsilon\)-optimality of the objective function for both IRWA and its ADAL alternative [1402.1917]. For equality-constrained nonlinear optimization, Fletcher’s smooth exact penalty gives a further reinterpretation. Estrin, Friedlander, Orban, and Saunders implement \[ \phi_\sigma(x):=f(x)-c(x)^T y_\sigma(x), \] where \[ y_\sigma(x):= \arg\min_{y} \left\{ \|A(x)y-g(x)\|_2^2+\sigma\, c(x)^T y \right\}, \] and show that the computational kernel is a structured linear system that can be reused for the penalty value, gradient, and Hessian-vector products [1910.04300]. The penalty is locally exact for finite \(\sigma\) above the threshold \[ \sigma^* := \lambda_{\max}^+\!\left( P(x^*)\, H(x^*,y^*)\, P(x^*) \right), \] and the regularized continuation \[ \delta_k \gets \max\left\{ \min\left\{\|\nabla \phi_\sigma(x_k;\delta_{k-1})\|,\ \delta_{k-1}\right\}, \ \nu(\delta_{k-1}) \right\} \] yields a genuine outer-inner iterative penalty/regularization method with local \(R\)-quadratic convergence [1910.04300]. ## 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 \[ \min_{\mathbf{x}} \quad g(\mathbf{x}) + \sum_{i=1}^K h_i(\mathbf{x}_i) \quad\text{s.t.}\quad \sum_{i=1}^K \mathbf{A}_i \mathbf{x}_i = \mathbf{b},\quad \mathbf{x}_i\in\mathcal X_i. \] The augmented penalty function is \[ L_\rho(\mathbf{x},\mathbf{y}) = g(\mathbf{x}) + \sum_{i=1}^K h_i(\mathbf{x}_i) + \left\langle \mathbf{y}, \sum_{i=1}^K \mathbf{A}_i \mathbf{x}_i - \mathbf{b} \right\rangle + \frac{1}{2\rho} \left\| \sum_{i=1}^K \mathbf{A}_i \mathbf{x}_i - \mathbf{b} \right\|^2. \] At outer iteration \(r\), the method approximately solves the augmented subproblem and then either updates the multiplier, \[ \mathbf{y}^{r+1} = \mathbf{y}^r + \frac{1}{\rho^r}\mathbf{c}(\mathbf{x}^{r+1}), \] if feasibility is adequate, or decreases \(\rho\) via \[ \rho^{r+1} = \theta \rho^r,\qquad 0<\theta<1, \] if feasibility is not [1712.04767]. 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 \[ \min_{v \in V} \frac12 \langle Av, v\rangle - \langle f, v\rangle \quad \text{subject to } v \le g. \] The complementarity operator \[ M(\phi,\varphi) := \phi - [\phi + c\varphi]_+ \] satisfies \[ M(\phi,\varphi)=0 \Longleftrightarrow \phi \ge 0,\ \varphi \le 0,\ \langle \phi,\varphi\rangle = 0, \] so the primal equation is \(M(f-Au,\ u-g)=0\) [2201.02425]. The smoothed penalty/Newton step uses \[ \alpha_\rho^k = [\,f-Au^k + c(u^k-g)\,]_\rho' \] and \[ \bigl((I-\alpha_\rho^k)A + \alpha_\rho^k c\bigr)\delta u = M(f-Au^k,\ u^k-g), \qquad u^{k+1}=u^k+\delta u. \] The penalty field \(\rho(x)\) is recomputed from the current residual by solving \[ \rho - \epsilon \Delta \rho = \gamma \, |M(f-Au,\ u-g)| \quad \text{in } \Omega, \qquad \mathbf n \cdot \nabla \rho = 0 \quad \text{on } \partial\Omega. \] As the residual decreases, \(\rho\to 0\), and the derivative weight converges to the active-set indicator: \[ \alpha_\rho^k \to \mathds{1}_{\mathcal A^k}. \] This establishes that the primal-dual active set method is exactly the \(\rho\downarrow 0\) limit of the regularized penalty/Newton method [2201.02425]. For large collections of linear inequalities, \[ \langle a_i,x\rangle - b_i \le 0,\qquad i\in[m], \] a Huber-loss-based single-loop incremental penalty method uses the time-varying penalized objective \[ F_k(x)=f(x)+\gamma_k H_k(x),\qquad H_k(x)=\frac{1}{m}\sum_{i=1}^m h_{\delta_k}(x;a_i,b_i), \] with one-sided Huber penalty \[ h_\delta(x;a,b)= \begin{cases} \frac{\langle a,x\rangle-b}{\|a\|}, & \langle a,x\rangle-b>\delta,\\[1mm] \frac{(\langle a,x\rangle-b+\delta)^2}{4\delta\|a\|}, & -\delta\le \langle a,x\rangle-b\le \delta,\\[1mm] 0, & \langle a,x\rangle-b<-\delta . \end{cases} \] The stochastic update is \[ x_{k+1} = x_k - s_k\Bigl[\tilde\nabla f(x_k)+\gamma_k\nabla h_{\delta_k}(x_k;a_{i_k},b_{i_k})\Bigr], \] with conditions \[ \gamma_k\to\infty,\qquad \sum_{k=1}^\infty s_k=\infty,\qquad \sum_{k=1}^\infty s_k\gamma_k\delta_k<\infty,\qquad \sum_{k=1}^\infty s_k^2\gamma_k^2<\infty. \] Under convexity assumptions, the iterates converge to a solution almost surely and in expectation, with weighted-average rates \[ O\!\left(\frac{\ln^{1/2+\epsilon} k}{\sqrt k}\right) \quad\text{for convex }f, \qquad O\!\left(\frac{\ln^{\epsilon} k}{k}\right) \quad\text{for strongly convex }f \] [2311.00874]. These examples also delimit the notion. A penalty-function-based decentralized bilevel algorithm replaces consensus constraints by graph quadratic penalties, for example \[ \frac{1}{2\alpha}x^\top(I_{nd_1}-\acute W)x, \] and proves \(\|x_i^*-\dot x^*\|=\mathcal O(\alpha)\) as \(\alpha\to 0\), but it does not run a classical outer penalty schedule; \(\alpha\) and \(\gamma\) are fixed design parameters in a single alternating decentralized algorithm [2211.04088]. ## 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 \[ Q_h=\operatorname{div}X_{h,\sim}=\operatorname{div}V_h, \] the pressure can be represented as the discrete divergence of an auxiliary velocity-like variable [2509.17899]. The inhomogeneous Dirichlet problem uses compatible boundary data \(g_h\in X_{h,\sim}\), and the iteration is initialized with \[ \phi_{h,0}=0. \] At step \(i\), \[ \nu (\nabla u_{h,i},\nabla v_h)+\rho (\operatorname{div}u_{h,i},\operatorname{div}v_h) = \langle f,v_h\rangle + (\operatorname{div}\phi_{h,i-1},\operatorname{div}v_h) \qquad \forall v_h\in V_h, \] with \(u_{h,i}=g_h\) on \(\partial\Omega\), followed by \[ \phi_{h,i} := \phi_{h,i-1} - \rho\, u_{h,i}, \qquad p_{h,i}:=\operatorname{div}\phi_{h,i}. \] Equivalently, \[ p_{h,i}=p_{h,i-1}-\rho\,\operatorname{div}u_{h,i}. \] The convergence factor is \[ \theta := \frac{\nu}{\nu+\rho \overline\beta^2}<1, \] and the paper proves \[ \|\nabla(u_h-u_{h,i})\|_{L^2(\Omega)} \le \theta^{\,i-1}\, \|\nabla(u_h-u_{h,1})\|_{L^2(\Omega)}, \] together with pressure convergence and asymptotic pressure robustness as \(i\to\infty\) [2509.17899]. A notable feature specific to this exactly divergence-free setting is monotonic decrease of the divergence norm: \[ \|\operatorname{div}u_{h,i}\|_{L^2(\Omega)} \le \|\operatorname{div}u_{h,i-1}\|_{L^2(\Omega)}. \] 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 \(\mathcal O(p^d)\) velocity unknowns per element, whereas SCIP reduces the iterative global system to \(\mathcal O(p^{d-1})\) boundary unknowns by using specially constructed spaces \(\tilde{\mathbf X}_B\) and \(\tilde{\mathbf X}_B^\dagger\) [2301.01818]. The reduced iteration is \[ a_\lambda(\tilde{\mathbf u}^n,\mathbf v) = L(\mathbf v)+(\operatorname{div}\tilde{\mathbf w}^n,\operatorname{div}\mathbf v) \qquad \forall \mathbf v\in \tilde{\mathbf X}_B^\dagger, \] with \[ \tilde{\mathbf w}^{n+1} = \tilde{\mathbf w}^n - \lambda \tilde{\mathbf u}^n. \] After convergence, the interior Stokes corrections are solved once, element by element. The resulting error bound is still geometric: \[ \|\operatorname{div}\tilde{\mathbf u}^n\| \le \sqrt d \left[ \frac{M(M+\alpha)^4}{\alpha^4\beta_X^2\lambda} \right]^n \|\tilde{\mathbf u}-\tilde{\mathbf u}^0\|_1, \] although the theoretical lower bound on \(\lambda\) is larger than for the standard method [2301.01818]. For steady incompressible Navier–Stokes, the classical iterated penalty Picard method uses \[ u_k\cdot \nabla u_{k+1} + \nabla p_{k+1} - \nu \Delta u_{k+1} = f, \qquad \varepsilon p_{k+1} + \nabla\cdot u_{k+1} = \varepsilon p_k. \] The associated fixed-point map \(G\) is Lipschitz continuous and Lipschitz continuously Fréchet differentiable in the norm \[ \|(v,q)\|_X := \sqrt{\nu \|\nabla v\|^2 + \varepsilon \|q\|^2}, \] which permits the application of Anderson acceleration theory [2105.09339]. The accelerated residual estimate has leading factor \[ \theta_k(1-\beta_k+\beta_k C_L), \] 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 \(1\), unaccelerated IPP can fail on challenging cases, whereas AAIPP with moderate or large memory converges effectively [2105.09339]. ## 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 \[ Ax=y \] in Hilbert spaces, Jin and Zou propose the two-step scheme \[ \xi_{n+1}=\xi_n - t_n A^*(\alpha_n I +AA^*)^{-1} (A x_n-y), \qquad x_{n+1} = \arg \min_{x\in \mathcal X} \left\{ \Theta(x) - \langle \xi_{n+1}, x \rangle \right\}, \] where \(\Theta\) is a proper, lower semicontinuous, convex penalty [1402.6544]. The first step involves only the operator \(A\); the second only the penalty. When \[ \Theta(x)=\frac12\|x\|^2, \] the method reduces to classical nonstationary iterated Tikhonov regularization. With nonsmooth penalties such as \[ \Theta(x)=\frac{1}{2\beta}\|x\|^2+\|x\|_1 \quad\text{or}\quad \Theta(x)=\frac{1}{2\beta}\|x\|^2+\mathrm{TV}(x), \] the scheme promotes sparsity or piecewise constancy while remaining computationally split [1402.6544]. The Banach-space counterpart uses the Bregman distance \[ D_\xi \Theta(\bar x,x):=\Theta(\bar x)-\Theta(x)- \langle \xi, \bar x-x\rangle, \] and defines the noisy-data iterate by \[ x_n^\delta \in \arg \min_{x\in D(F)} \left\{\frac{1}{r} \|F(x)-y^\delta\|^r + \alpha_n D_{\xi_{n-1}^\delta} \Theta(x, x_{n-1}^\delta)\right\}, \] \[ \xi_n^\delta = \xi_{n-1}^\delta-\frac{1}{\alpha_n} F'(x_n^\delta)^* J_r(F(x_n^\delta)-y^\delta). \] Under the stated assumptions—uniform convexity of \(\Theta\), weak closedness and differentiability of \(F\), and \(\sum_{n=1}^\infty \alpha_n^{-1}=\infty\)—the exact-data iterates converge strongly, with \[ \lim_{n\rightarrow \infty} \|x_n-x_*\|=0, \qquad \lim_{n\rightarrow \infty} D_{\xi_n} \Theta(x_*, x_n)=0, \] and the discrepancy principle \[ \|F(x_{n_\delta}^\delta)-y^\delta\|\leq \tau\delta<\|F(x_n^\delta)-y^\delta\|, \qquad 0\leq n0 \text{ very small}. \] Instead of directly minimizing \(\Phi\), the modified augmented Lagrangian method introduces \(\xi\) and rewrites the problem as \[ \min_{x,\xi}\quad f(x) + \frac{\omega}{2}\,\|\xi\|_2^2 \quad\text{s.t.}\quad c(x) + \omega \xi = 0. \] Eliminating \(\xi\) after applying ALM leads to the modified root function \[ F(x,\lambda)= \begin{pmatrix} \nabla f(x)-\nabla c(x)(\lambda_k+\lambda)\\ c(x)+\omega\lambda_k+(\omega+\rho)\lambda \end{pmatrix} =0, \] which is designed for problems where the direct quadratic-penalty objective is badly scaled [1804.08072]. 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. Source: https://www.emergentmind.com/topics/iterated-penalty-method