Implicit-Explicit Trust Region Method for Computing Second-Order Stationary Points of A Class of Landau Models

Abstract: We propose an implicit-explicit trust region method for computing second-order stationary points of a class of Landau-type free energy functionals, which correspond to physically (meta-)stable phases. The proposed method is demonstrated through the Landau-Brazovskii (LB) model in this work, while broader applicability to more Landau models of the similar type is straightforwardly extended. The LB energy functional is discretized via the Fourier pseudospectral method, which yields a finite-dimensional nonconvex optimization problem. By exploiting the Hessian structure, specifically, that the interaction potential is diagonal in reciprocal space whereas the bulk energy is diagonal in physical space, we design an adaptive implicit-explicit solver for the trust region subproblem. This solver utilizes the fast Fourier transform to perform efficient matrix-vector products, significantly reducing computational complexity while ensuring provable convergence to the global minimizer of the subproblem. In contrast to existing algorithms that target first-order stationary points, our proposed method can converge to a second-order stationary state, corresponding to a local minimum with theoretical convergence guarantees. Numerical experiments on the LB model demonstrate that the proposed approach efficiently escapes saddle points and significantly outperforms existing first-order schemes. Furthermore, we successfully identify the stable region of the FDDD phase, a structure previously unreported in the LB phase diagram.

• The paper introduces an implicit-explicit trust region algorithm whose FFT-efficient inner solver globally minimizes nonconvex quadratic subproblems and supports convergence to second-order stationary points.
• The method costs O(N log N) per Hessian-vector or inner iteration, escapes saddle points missed by first-order schemes, and identifies HEX, BCC, and lower-energy FDDD phases in numerical tests.
• The study maps a previously unreported stable FDDD region in the Landau–Brazovskii phase diagram, while noting open issues involving degenerate gradients, step-size certification, and broader model validation.

Problem setting and motivation

Landau-type free energy functionals, of the form $E(\psi) = G(\psi) + F(\psi)$ with a nonlocal interaction term $G$ and a nonlinear bulk potential $F$, underpin the modeling of phases and phase transitions in systems ranging from block copolymers to microemulsions. The physically meaningful states—stable and metastable phases—correspond to second-order stationary points (SP-IIs): points where the gradient vanishes and the Hessian is positive semi-definite. First-order stationary points (SP-Is), by contrast, include saddle points that are energetically unstable. Existing computational approaches for such models—direct solution of Euler–Lagrange equations, energy-stable gradient flow integrators (semi-implicit schemes, SAV, IEQ, ETD), and first-order optimization methods such as Bregman proximal gradient—guarantee convergence only to SP-Is, and can therefore terminate at saddle points rather than genuine local minima.

The paper by Bao, Deng, Jiang, and Zhang addresses this gap with an implicit-explicit trust region (IMEX-TR) method that carries provable convergence guarantees to SP-IIs. The method is developed and validated on the Landau–Brazovskii (LB) model,

$$E(\psi) = \frac{1}{|\Omega|}\int_\Omega \frac{1}{2}[(\Delta+1)\psi]^2 + \frac{\tau}{2}\psi^2 - \frac{\gamma}{6}\psi^3 + \frac{1}{24}\psi^4\, d\mathbf{x},$$

subject to mass conservation $\int_\Omega \psi\, d\mathbf{x} = 0$, but the authors state that extension to other Landau models of similar structure is straightforward.

Discretization and Hessian structure

The LB functional is discretized on a Bravais-periodic torus via the Fourier pseudospectral method. Writing the discrete order parameter as $\hat{\bm{\Psi}} = \mathscr{F}_N \bm{\Psi}$, the interaction term contributes a Hessian component $\nabla^2 G = \bm{\Lambda}$, diagonal in reciprocal space with nonnegative entries $[1-\mathbf{k}^\top\mathbf{k}]^2$, while the bulk term contributes $$\nabla^2 F = \mathscr{F}_N \bm{\Gamma}' \mathscr{F}_N^{-1}$$, diagonal in physical space. This complementary diagonal structure is the key structural observation exploited throughout: it permits matrix-vector products with the full Hessian at FFT cost, $\mathcal{O}(N\log N)$ per application, avoiding the $G$0 cost of dense operations. Mass conservation is eliminated exactly through the projection $G$1, reducing the problem to unconstrained minimization over $G$2.

The IMEX-TR outer iteration

At each outer iteration, the method builds the exact quadratic model $G$3 and solves the trust region subproblem $G$4 to global optimality. The outer loop follows the standard acceptance/rejection logic with an adaptive radius update, augmented with a CONTRACT procedure that adjusts $G$5 based on the multiplier $G$6 when steps are rejected. Under Lipschitz continuity of the gradient and Hessian—which the authors establish from coercivity of $G$7 and boundedness of the third derivative on the sublevel set—the iterate sequence satisfies

$G$8

i.e., every accumulation point is an SP-II. This guarantee is inherited from classical trust region theory given globally optimal subproblem solutions; the substantive contribution is therefore the subproblem solver itself.

Global solver for the nonconvex trust region subproblem

The subproblem is nonconvex whenever the Hessian has a negative eigenvalue, and computing its global minimizer is generally hard. The proposed inner iteration treats the reciprocal-space diagonal part $G$9 implicitly and the physical-space diagonal part $F$0 explicitly:

$F$1

Because $F$2 is diagonal and positive definite, each update reduces to element-wise division in reciprocal space, and the multiplier $F$3 is chosen adaptively—zero when the unconstrained step is feasible, otherwise the root of $F$4 found by Newton's method—to enforce feasibility $F$5 with complementarity. Each inner iteration costs one FFT pair plus pointwise operations.

The convergence analysis proceeds in three stages. First, for step sizes satisfying $F$6, the scheme is a descent method for both $F$7 and the Lagrangian $F$8, with a per-step decrease of at least $F$9. Second, using the Kurdyka–Łojasiewicz property of $$E(\psi) = \frac{1}{|\Omega|}\int_\Omega \frac{1}{2}[(\Delta+1)\psi]^2 + \frac{\tau}{2}\psi^2 - \frac{\gamma}{6}\psi^3 + \frac{1}{24}\psi^4\, d\mathbf{x},$$0 together with a bound relating the Lagrangian gradient norm to the step length, the sequence $$E(\psi) = \frac{1}{|\Omega|}\int_\Omega \frac{1}{2}[(\Delta+1)\psi]^2 + \frac{\tau}{2}\psi^2 - \frac{\gamma}{6}\psi^3 + \frac{1}{24}\psi^4\, d\mathbf{x},$$1 is shown to be Cauchy, and every limit point $$E(\psi) = \frac{1}{|\Omega|}\int_\Omega \frac{1}{2}[(\Delta+1)\psi]^2 + \frac{\tau}{2}\psi^2 - \frac{\gamma}{6}\psi^3 + \frac{1}{24}\psi^4\, d\mathbf{x},$$2 satisfies the KKT conditions. Third—and this is the strongest claim in the paper—if $$E(\psi) = \frac{1}{|\Omega|}\int_\Omega \frac{1}{2}[(\Delta+1)\psi]^2 + \frac{\tau}{2}\psi^2 - \frac{\gamma}{6}\psi^3 + \frac{1}{24}\psi^4\, d\mathbf{x},$$3 along the eigenvector of the most negative eigenvalue $$E(\psi) = \frac{1}{|\Omega|}\int_\Omega \frac{1}{2}[(\Delta+1)\psi]^2 + \frac{\tau}{2}\psi^2 - \frac{\gamma}{6}\psi^3 + \frac{1}{24}\psi^4\, d\mathbf{x},$$4 and the initial step is taken as $$E(\psi) = \frac{1}{|\Omega|}\int_\Omega \frac{1}{2}[(\Delta+1)\psi]^2 + \frac{\tau}{2}\psi^2 - \frac{\gamma}{6}\psi^3 + \frac{1}{24}\psi^4\, d\mathbf{x},$$5, then an induction argument shows $$E(\psi) = \frac{1}{|\Omega|}\int_\Omega \frac{1}{2}[(\Delta+1)\psi]^2 + \frac{\tau}{2}\psi^2 - \frac{\gamma}{6}\psi^3 + \frac{1}{24}\psi^4\, d\mathbf{x},$$6 for all $$E(\psi) = \frac{1}{|\Omega|}\int_\Omega \frac{1}{2}[(\Delta+1)\psi]^2 + \frac{\tau}{2}\psi^2 - \frac{\gamma}{6}\psi^3 + \frac{1}{24}\psi^4\, d\mathbf{x},$$7. Combined with the characterization that the global minimizer is the unique stationary point satisfying $$E(\psi) = \frac{1}{|\Omega|}\int_\Omega \frac{1}{2}[(\Delta+1)\psi]^2 + \frac{\tau}{2}\psi^2 - \frac{\gamma}{6}\psi^3 + \frac{1}{24}\psi^4\, d\mathbf{x},$$8, this establishes that $$E(\psi) = \frac{1}{|\Omega|}\int_\Omega \frac{1}{2}[(\Delta+1)\psi]^2 + \frac{\tau}{2}\psi^2 - \frac{\gamma}{6}\psi^3 + \frac{1}{24}\psi^4\, d\mathbf{x},$$9 is the global minimizer of the subproblem for any step size in an interval $\int_\Omega \psi\, d\mathbf{x} = 0$0. Consequently, the outer IMEX-TR iteration inherits SP-II convergence with certainty, not merely asymptotically.

Two assumptions deserve note here. The global optimality argument requires $\int_\Omega \psi\, d\mathbf{x} = 0$1; the degenerate case where the gradient has no component along the most-negative-curvature direction is not covered by the theorem. Additionally, the admissible step size window depends on problem-dependent quantities ($\int_\Omega \psi\, d\mathbf{x} = 0$2, $\int_\Omega \psi\, d\mathbf{x} = 0$3, $\int_\Omega \psi\, d\mathbf{x} = 0$4), and the practical choice $\int_\Omega \psi\, d\mathbf{x} = 0$5 used in experiments is empirical rather than derived from these bounds.

Numerical results

Three benchmark configurations of the LB model compare IMEX-TR against AA-BPG and several gradient flow schemes (SSIS1, SSIS2, SAV, ETD). The findings are consistent across examples:

• In Example 1 ($\int_\Omega \psi\, d\mathbf{x} = 0$6, $\int_\Omega \psi\, d\mathbf{x} = 0$7), first-order methods stagnate at the initial lamellar (LAM) configuration, which eigenvalue analysis confirms is a saddle point with persistent negative Hessian eigenvalues; IMEX-TR converges to the hexagonal (HEX) phase as an SP-II.
• In Example 2 ($\int_\Omega \psi\, d\mathbf{x} = 0$8, $\int_\Omega \psi\, d\mathbf{x} = 0$9), first-order methods collapse to the disordered state, again a saddle point, while IMEX-TR reaches the body-centered cubic (BCC) phase.
• In Example 3 ($\hat{\bm{\Psi}} = \mathscr{F}_N \bm{\Psi}$0, $\hat{\bm{\Psi}} = \mathscr{F}_N \bm{\Psi}$1), all methods reach HEX, but IMEX-TR escapes the HEX basin and identifies a cubic FDDD pattern with energy $\hat{\bm{\Psi}} = \mathscr{F}_N \bm{\Psi}$2, lower than the HEX value $\hat{\bm{\Psi}} = \mathscr{F}_N \bm{\Psi}$3.

Robustness experiments activate IMEX-TR at three stages of a first-order trajectory—initial condition, intermediate state, and converged saddle point—and in all cases the method converges to the same SP-II, demonstrating that it can be used either standalone or as a post-processing escape mechanism for solvers trapped at saddles.

Discovery of the FDDD stable region

The FDDD result is extended into a systematic search over the $\hat{\bm{\Psi}} = \mathscr{F}_N \bm{\Psi}$4 parameter plane. The resulting phase diagram confirms a finite stable region for the cubic FDDD phase in the LB model—a phase previously documented for AB-type diblock copolymers but unreported in existing LB phase diagrams. Importantly, the identified FDDD region does not perturb the stable regions of the face-centered cubic spheres, A15, and $\hat{\bm{\Psi}} = \mathscr{F}_N \bm{\Psi}$5 phases reported in prior work, so the new diagram refines rather than contradicts the established picture. This constitutes a concrete physical contribution attributable directly to the algorithm's ability to traverse negative curvature directions.

Limitations and open questions

Several constraints qualify the results. The convergence theory relies on the KL property of the Lagrangian and on the nondegeneracy condition $\hat{\bm{\Psi}} = \mathscr{F}_N \bm{\Psi}$6; behavior when the gradient is orthogonal to the minimal eigenvector is unresolved within the presented analysis. The theoretical step size bounds involve constants that are not computed explicitly in practice, so the empirical choice $\hat{\bm{\Psi}} = \mathscr{F}_N \bm{\Psi}$7 is not certified a priori. The numerical validation is confined to the single-component LB model in three dimensions with periodic boundaries; while the authors assert straightforward extensibility to other Landau-type functionals and discretizations, no quantitative evidence for such extensions is provided. Finally, the paper does not report wall-clock or iteration-count comparisons quantifying the overhead of the inner loop relative to first-order methods, leaving the computational trade-off between guaranteed second-order convergence and per-iteration cost only implicitly characterized.

Conclusion

This paper couples a trust region framework, whose outer iterates provably converge to SP-IIs, with a structure-exploiting implicit-explicit inner solver that computes the global minimizer of each nonconvex subproblem at $\hat{\bm{\Psi}} = \mathscr{F}_N \bm{\Psi}$8 cost per iteration. The combination yields certified convergence to physically (meta-)stable states of the Landau–Brazovskii model, demonstrably escaping the saddle points at which first-order schemes stagnate, and leads to the identification of a previously unreported stable FDDD region in the LB phase diagram. The main open issues are the degenerate-gradient case in the global optimality proof, a priori step size selection, and empirical validation of the claimed generality beyond the LB benchmark.