The Method of Moving Asymptotes (MMA)

Table of Contents

• Overview
• Original Optimization Problem
• Convex Approximation
  • Asymptotic Approximation
  • Construction from Function Value and Gradient
  • Resulting Convex Subproblem
  • Key Properties of This Approximation
  • Intuition
• Moving Asymptotes
• Subproblem Solution
• KKT Convergence
• Implementation Details
• Notes
  • Note on the DIP subsolver formulation
  • Key difference from the primal MMA subproblem
  • Practical implication
  • Relation to the topopt_in_petsc implementation
• References

Overview

The Method of Moving Asymptotes (MMA) is a gradient-based optimization algorithm widely used in topology optimization. It is particularly effective for large-scale, constrained, non-linear optimization problems.

The method transforms the original non-convex optimization problem into a sequence of strictly convex subproblems that are easier to solve.

Original Optimization Problem

The MMA implementation in Neko-TOP solves problems of the form

\[ \begin{aligned} \min_{x,z,y} \quad & f_0(x) + a_0 z + \sum_{i=1}^m \left( c_i y_i + \frac{1}{2} d_i y_i^2 \right) \\ \text{s.t.} \quad & f_i(x) - a_i z - y_i \le 0, \quad i = 1, \dots, m, \\ & x_j^{\min} \le x_j \le x_j^{\max}, \quad j = 1, \dots, n, \\ & z \ge 0, \quad y_i \ge 0. \end{aligned} \]

Here:

• \(x\): design variables
• \(f_0\): objective function
• \(f_i\): constraint functions
• \(y_i\), \(z\): auxiliary variables for constraint relaxation

---

Convex Approximation

At each iteration, MMA constructs a separable convex approximation of the original non-linear problem by replacing each function \( f_i(x) \) with an asymptotic approximation built from (note that \(f_0\) is also approximated in the same way):

• the current function value \( f_i(x^k) \)
• the first-order sensitivities \( \nabla f_i(x^k) \)

---

Asymptotic Approximation

Instead of a standard Taylor expansion, MMA uses a asymptotic model of the form:

\[ f_i(x) \;\approx\; \tilde{f}_i(x) = \sum_{j=1}^n \left( \frac{p_{ij}}{u_j - x_j} + \frac{q_{ij}}{x_j - l_j} \right) - b_i \]

where:

• \( l_j \), \( u_j \) are the moving asymptotes
• \( p_{ij}, q_{ij} \) are non-negative coefficients constructed from the gradients
• \( b_i \) is chosen such that the approximation is exact at the current iteration: \( \tilde{f}_i(x^k) = f_i(x^k) \)

---

Construction from Function Value and Gradient

The coefficients \( p_{ij} \), \( q_{ij} \) are derived from the sensitivities \( \frac{\partial f_i}{\partial x_j} \) to ensure:

• First-order consistency:

  \[ \nabla \tilde{f}_i(x^k) \approx \nabla f_i(x^k) \]

• Convexity (by enforcing \( p_{ij}, q_{ij} \ge 0 \))

In practice (as implemented in mma_gensub):

\[ \begin{aligned} p_{ij} &\sim \max\!\left(\frac{\partial f_i}{\partial x_j}, 0\right) \cdot (u_j - x_j)^2 \\ q_{ij} &\sim \max\!\left(-\frac{\partial f_i}{\partial x_j}, 0\right) \cdot (x_j - l_j)^2 \end{aligned} \]

with small regularization terms added for numerical stability. The exact form depends on the subsolver.

For the dip subsolver, on both the CPU and device backends, it follows MMA::GenSub in the topopt_in_petsc code by Aage (with constraintModification = false): only the objective ( \( i = 0 \)) is regularized,

\[ \begin{aligned} p_{0j} &= \left( \max\!\left(\frac{\partial f_0}{\partial x_j}, 0\right) + 0.001 \left|\frac{\partial f_0}{\partial x_j}\right| + \frac{0.5 \cdot 10^{-6}}{u_j - l_j} \right) (u_j - x_j)^2, & p_{ij} &= \max\!\left(\frac{\partial f_i}{\partial x_j}, 0\right) (u_j - x_j)^2, \\[8pt] q_{0j} &= \left( \max\!\left(-\frac{\partial f_0}{\partial x_j}, 0\right) + 0.001 \left|\frac{\partial f_0}{\partial x_j}\right| + \frac{0.5 \cdot 10^{-6}}{u_j - l_j} \right) (x_j - l_j)^2, & q_{ij} &= \max\!\left(-\frac{\partial f_i}{\partial x_j}, 0\right) (x_j - l_j)^2, \end{aligned} \]

for \( i = 1, \dots, m \).

For the pdip subsolver, on both the CPU and device backends, it follows mmasub by Svanberg, where the objective and all constraints ( \( i = 0, \dots, m \)) are regularized:

\[ \begin{aligned} p_{ij} &= \left( 1.001 \max\!\left(\frac{\partial f_i}{\partial x_j}, 0\right) + 0.001 \max\!\left(-\frac{\partial f_i}{\partial x_j}, 0\right) + \frac{10^{-5}}{\max(x_{\text{diff},j}, 10^{-5})} \right) (u_j - x_j)^2 \\[8pt] q_{ij} &= \left( 0.001 \max\!\left(\frac{\partial f_i}{\partial x_j}, 0\right) + 1.001 \max\!\left(-\frac{\partial f_i}{\partial x_j}, 0\right) + \frac{10^{-5}}{\max(x_{\text{diff},j}, 10^{-5})} \right) (x_j - l_j)^2 \end{aligned} \]

Here \( x_{\text{diff},j} = x_j^{\max} - x_j^{\min} \), restricted by the move limit if one is set.

---

Resulting Convex Subproblem

Using these approximations, MMA solves the following convex separable problem:

\[ \min_x \sum_{j=1}^n \left( \frac{p_{0j}}{u_j - x_j} + \frac{q_{0j}}{x_j - l_j} \right) + a_0 z + \sum_{i=1}^m \left(c_i y_i + \frac{1}{2} d_i y_i^2\right) \]

subject to:

\[ \sum_{j=1}^n \left( \frac{p_{ij}}{u_j - x_j} + \frac{q_{ij}}{x_j - l_j} \right) + a_i z + y_i \le b_i \]

---

Key Properties of This Approximation

• Uses only local information:
  • function values \( f_i(x^k) \)
  • gradients \( \nabla f_i(x^k) \)
• Separable in \( x_j \) → efficient for large-scale problems
• Convex by construction → guarantees well-posed subproblems
• Asymptotic behavior near bounds:
  • singularities at \( x_j \to l_j \) and \( x_j \to u_j \).
  • more conservative bounds ("alpha", "beta") for the updated design variables are chosen to prevent them from hitting bounds caused by too aggressive updates.

---

Intuition

Unlike a Taylor expansion, MMA builds a curvature-aware approximation using moving asymptotes where

\[ \frac{1}{u_j - x_j}, \quad \frac{1}{x_j - l_j} \]

act as barrier-like functions that:

• match the local gradient,
• enforce convexity,
• and guide the design smoothly within admissible bounds.

---

Moving Asymptotes

A key feature of MMA is the adaptive update of asymptotes:

• They control step size and stability for each design variable
• They prevent oscillations
• They improve convergence robustness

The update is governed by:

• asyinit (initial spacing)
• asyincr (expansion factor to push the asymptotes apart to help with faster convergence using larger steps)
• asydecr (contraction factor to pull the asymptotes together to take more conservative steps)

---

Subproblem Solution

The convex subproblem is solved using a primal-dual interior point method "pdip" and a pure dual interior point method "dip".

In this implementation both subsolvers support both CPU and device (GPU) execution, and the alignment of the dip subsolver with the topopt_in_petsc code (see the sections below) applies to both: the device backend reproduces the same subproblem, subsolver and KKT residual as the CPU backend, on both the CUDA and HIP backends; see Relation to the topopt_in_petsc implementation below for their verification status.

---

KKT Convergence

Convergence is evaluated using Karush-Kuhn-Tucker (KKT) conditions:

• Infinity norm: residumax
• Euclidean norm: residunorm

The optimizer stops when residumax is below optimization.solver.tolerance. The residual is evaluated at the updated design, with the multipliers of the last subproblem.

For the dip subsolver, on both the CPU and device backends, the residual follows MMA::KKTresidual in the topopt_in_petsc code by Aage. It consists of

• the gradient of the Lagrangian, \( \partial f_0 / \partial x_j + \sum_{i} \lambda_i \partial f_i / \partial x_j - \xi_j + \eta_j \), where the bound multipliers \( \xi_j, \eta_j \) are estimated for design variables within \( 10^{-5} \) of \( x_j^{\min} \) or \( x_j^{\max} \),
• the complementarity of the bounds, \( \xi_j (x_j - x_j^{\min}) \) and \( \eta_j (x_j^{\max} - x_j) \),
• the complementarity of the constraints, \( \sum_{i} \lambda_i (a_i z + y_i - f_i(x)) \).

The true bounds \( x^{\min}, x^{\max} \) are used, not the move limited ones.

For the dip subsolver, the device backend also sets residumax to residunorm whenever the latter is NaN: the maximum reductions can drop a NaN, while the 2-norm keeps it, so a NaN in the design, the sensitivities or the constraint values never reads as convergence. The CPU backend does not have this guard.

Since each subproblem is only solved down to the barrier parameter epsimin, the constraint term does not drop much below \( m \) times the last barrier level, which is the smallest power of ten above epsimin (between epsimin and \( 10 \) epsimin). The tolerance should be set above that.

For the pdip subsolver the residual is the full KKT system of the primal-dual method (Svanberg).

---

Implementation Details

The implementation is encapsulated in the mma_t type and includes the following parameters that can be set in the case file:

Name	Description	Default
mma.max_iter	Max iterations for subproblem	100
mma.epsimin	Smallest barrier parameter of the subsolvers, must be positive	\( 10^{-9} \sqrt{m + n} \)
mma.asyinit	Initial asymptote distance	0.2
mma.asyincr	Asymptote expansion	1.05
mma.asydecr	Asymptote contraction	0.65
mma.move_limit	Move limit for updating the design variables	0.2
mma.a0	MMA param	1.0
mma.a	MMA param	0.0
mma.c	MMA param	1000.0
mma.d	MMA param	0.0
mma.backend	cpu or device	auto based on Neko backend
mma.subsolver	Subsolver type	dip
mma.scale	Scaling factor applied to constraint functions ( f_i ) and their sensitivities. This does not affect the objective function ( f_0 ). It is used to improve numerical conditioning when constraint magnitudes and sensitivities differ significantly from those of the objective function.	1.0
mma.auto_scale	If true, sensitivity and function values for the constraint \( f_i \) are scaled at each iteration with a different value such that we get \( f_1(x_k)= \)mma.scale. This would be an adaptive scaling based on the value of the first constraint.	false

---

Notes

• The constant term in the approximation (see Eq. 3.5 in Svanberg) is omitted, as it does not affect the minimization ( \( b_0 \) in the approximation of \( f_0 \)).

Note on the DIP subsolver formulation

The dip (Dual Interior Point) subsolver solves a slightly different but equivalent reformulation of the MMA subproblem.

Instead of directly minimizing the primal convex approximation in \((x, y, z)\), the method forms the Lagrangian dual problem:

\[ \Psi(\lambda) = \sum_{j=1}^{n} \min_{x_j} \left\{ L_x(x_j, \lambda) \;\middle|\; \alpha_j \le x_j \le \beta_j \right\} + \min_{z \ge 0} L_z(z, \lambda) + \sum_{i=1}^{m} \min_{y_i \ge 0} L_y(y_i, \lambda) \]

and then solves:

\[ \max_{\lambda \ge 0} \; \Psi(\lambda) \]

Key difference from the primal MMA subproblem

The DIP formulation uses the Lagrangian:

\[ \begin{aligned} L(x,y,z,\lambda) = &\sum_{j=1}^{n} \left( \frac{p_{0j} + \sum_{i=1}^{m} \lambda_i p_{ij}}{u_j - x_j} + \frac{q_{0j} + \sum_{i=1}^{m} \lambda_i q_{ij}}{x_j - l_j} \right) - \sum_{i=1}^{m} \lambda_i b_i \\ &+ \sum_{i=1}^{m} \left[ (c_i - \lambda_i) y_i + \frac{1}{2} y_i^2 \right] + \left(a_0 - \sum_{i=1}^{m} \lambda_i a_i\right) z + \frac{1}{20} z^2 \end{aligned} \]

so the quadratic terms for \(y_i\) and \(z\) are enforced to make sure that we can solve the minimization problems, analytically. On both the CPU and device backends the weight \( 1/20 \) of \( z^2 \) follows the topopt_in_petsc code by Aage, which gives \( z = \max\left(0, 10 \left(\sum_{i} \lambda_i a_i - a_0\right)\right) \).

Practical implication

Compared to a standard primal MMA subsolve:

• the quadratic terms for \(y_i\) and \(z\) are considered in the solver regardless of what parameters user set in the case file.
• the primal subproblem is not solved directly
• instead, each iteration computes:
  • analytical minimizers in \(x_j, y_i, z\)
  • and then performs a dual ascent on \(\lambda\)
• DIP is cheaper per iteration

Relation to the topopt_in_petsc implementation

The dip subsolver follows MMA.cc in the topopt_in_petsc code by Aage (subproblem, dual solver and its stopping rules, KKT residual). On the CPU backend it reproduces the iterates of MMA.cc to round-off on the same input. The remaining differences from MMA.cc are

• the default asymptote parameters: Neko-TOP uses 0.2, 1.05, 0.65 for asyinit, asyincr, asydecr, while MMA.cc uses 0.5, 1.2, 0.7,
• the initial multipliers are \( \lambda_i = \max(1, c_i/2) \) instead of \( c_i/2 \), which only differs for \( c_i < 2 \),
• the span \( x_j^{\max} - x_j^{\min} \) is bounded below by \( 10^{-5} \) also for the initial asymptotes, so fixed variables do not break the approximation,
• variables with \( x_j \) within machine precision of \( \alpha_j \) or \( \beta_j \) are treated as active in the dual Hessian, where MMA.cc tests whether the unclamped minimizer lies outside \( [\alpha_j, \beta_j] \),
• the dual Newton system is solved with a pivoted LU factorization (LAPACK).

The device backend implements the same dip subproblem, dual subsolver, stopping rule and KKT residual as the CPU backend, and agrees with it to round-off on the same finite input; the linear system is instead factorized with cuSOLVER (CUDA) or hipSOLVER (HIP). Both the CUDA and HIP backends implement this; it is verified on CUDA hardware, and the HIP kernels mirror the CUDA ones but have not been compiled or run on AMD hardware, so HIP's agreement with the CPU and CUDA backends is unverified.

The topopt_in_petsc driver stops when the largest design change drops below \( 0.01 \) and does not evaluate the KKT residual. In Neko-TOP the design change criterion is optimization.solver.stop_design_change, next to the KKT tolerance.

The topopt_in_petsc driver also normalizes the objective to \( f_0(x^0) = 10 \) before every update. Neko-TOP passes the objective unscaled, while constants such as the regularization \( 0.5 \cdot 10^{-6} \) and the default \( c_i = 1000 \) are absolute. For objectives of a very different magnitude, the weight of the objective can be used to bring it to a similar scale.

---

References

• Svanberg, K. (1987). The method of moving asymptotes—a new method for structural optimization. International Journal for Numerical Methods in Engineering, 24(2), 359-373.
• Svanberg, K. (1993). The Method of Moving Asymptotes (MMA) with Some Extensions. Optimization of Large Structural Systems, 191–207.
• Svanberg, K. (2002). A class of globally convergent optimization methods based on conservative convex separable approximations. SIAM Journal on Optimization, 12(2), 555-573.
• Aage, N., & Lazarov, B. S. (2013). Parallel framework for topology optimization using the method of moving asymptotes. Structural and Multidisciplinary Optimization, 47(4), 493-505.