ABSTRACT. Gradient methods are widely used in optimization problems. In practice, while the smoothness parameter can be estimated utilizing techniques such as backtracking, estimating the strong convexity parameter remains a challenge; moreover, even with the optimal parameter choice, convergence can be slow. In this work, we propose a framework for dynamically adapting the step size and momentum parameters in first-order gradient methods for the optimization problem, without prior knowledge of the strong convexity parameter. The main idea is to use the geometric average of the ratios of successive residual norms as an empirical estimate of the upper bound on the convergence rate, which in turn allows us to adaptively update the algorithm parameters. The resulting algorithms are simple to implement, yet efficient in practice, requiring only a few additional computations on existing information. The proposed adaptive gradient methods are shown to converge at least as fast as gradient descent for quadratic optimization problems. Numerical experiments on both quadratic and nonlinear problems validate the effectiveness of the proposed adaptive algorithms. The results show that the adaptive algorithms are comparable to their counterparts using optimal parameters, and in some cases, they capture local information and exhibit improved performance.

ABSTRACT. This presentation discusses interpolation-based derivative-free trust-region methods for the important black-box optimization, with a particular focus on a new hybrid neural-quadratic model-based optimization approach. We explain why interpolation models combining quadratic models with neural network surrogates are effective for optimization when true derivative information is unavailable, theoretically and numerically. The presentation also reviews the historical development of derivative-free optimization, highlights recent advances in hybrid model construction, and illustrates several motivating applications. Finally, potential directions for future research, including applications in quantum computing and artificial intelligence, will be discussed. This is joint work with Stefan M. Wild at the LBNL.

ABSTRACT. Many optimization algorithms can be viewed as evolution equations, and this perspective suggests importing ideas from numerical PDEs into optimization. In this talk, I discuss three such ideas: preconditioning transformations for adjoint systems of evolution equations, nonlinear splitting methods for adjoint-based optimization, and localized regularization of Hamilton--Jacobi and Hamilton--Jacobi--Bellman equations motivated by shock and caustic formation. The common theme is to modify the underlying dynamics, rather than only the outer optimization iteration, in order to improve robustness and efficiency. Preconditioning can reshape adjoint backpropagation in multiscale PDE-constrained problems, while nonlinear splitting leads to semi-implicit optimization methods with improved stability and computational efficiency. On the Hamilton--Jacobi side, information-geometric regularization provides a localized smoothing of kinks and caustics in value functions, avoiding the global oversmoothing associated with artificial viscosity. Together, these examples illustrate how splitting, preconditioning, and regularization ideas from numerical PDEs can lead to more effective optimization methods.

ABSTRACT. Data science methods (e.g., regression, deep learning) rely heavily on solving optimization problems with a large number of parameters and observations, which becomes a major hurdle when optimization solvers are computationally inefficient or do not scale. The efficiency and scalability of solvers can be understood using worst-case theoretical complexity analysis. These theoretical complexity analyses guide the selection of an appropriate algorithm to solve the optimization problem, which can often impact the quality of the data science application.

Typically, worst-case complexity is analyzed given a fixed error bound. However, we claim that analyzing algorithms for a fixed error bound is incomplete, and not an indicator of the limiting behavior of algorithms as we continue to decrease error bounds. We introduce a new definition for determining complexity rates in an asymptotic fashion, and show that the worst-case complexity for deterministic first-order methods cannot be achieved asymptotically. Moreover, we will construct functions over unbounded domains that can achieve slightly better complexity rates asymptotically.

Finally, we use different techniques for complexity analysis, namely average-case and smoothed-case analysis, to better indicate the performance of algorithms over practical, large-scale problems in data science in order to assist users with algorithm selection. We also validate the efficacy of our theoretical results using software packages developed in Julia.

ABSTRACT. Polytopal methods are a class of finite element methods that can be applied on polytopal meshes, i.e, meshes that are composed of general polygons or polyhedra. Polytopal methods have become increasingly popular since they can deal with complex geometries without the need for complex meshing tools. These often involve hybrid spaces, where degrees of freedom are attached to mesh entities of different dimensions.

The design of efficient, robust, scalable solvers for linear systems arising from these kind of discretizations is important to make them competitive with traditional methods on real world applications. One family of such scalable preconditioners are geometric multigrid methods.

Multigrid methods for nonconforming hybrid discretisations have been studied in the literature, mainly for HDG and HHO. These earlier approaches usually restricted to nested multigrid mesh hierarchies with planar faces at all levels, drastically reducing their applicability. Furthermore, the analysis of these methods is often restricted to conforming, nested, simplicial mesh hierarchies and relies on conforming piecewise linear spaces.

Efforts to circumvent this issue turn to non-nested hierarchies, where the coarser levels are defined on geometrically coarsened meshes with planar faces. However, these algorithms rely on complex of the geometrical algorithms, which have only been developed in 2D and never extended to 3D.

In this talk, we will present the first optimal geometric multigrid solver for hybrid high-order discretizations that can handle arbitrary polytopal agglomeration hierarchies in both two and three dimensions. The key ingredient is the use of modified skeleton spaces, which naturally accommodate non-planar interfaces arising during coarsening while reducing the number of degrees of freedom. We prove robust convergence with respect to the mesh size and the number of levels, and we validate our results numerically on a range of agglomeration-based mesh hierarchies. The approach extends naturally to other hybrid discretizations such as hybridizable discontinuous Galerkin and Weak Galerkin methods. This talk is based on `Geometric Multigrid solvers for Hybrid High-Order methods on polytopal meshes` (https://arxiv.org/abs/2603.00860).

ABSTRACT. Note: this abstract is based on a paper submitted (currently in review) to the SIAM Journal on Scientific Computing (SISC). A preprint can be found at the following link: https://arxiv.org/abs/2512.08555

1. Introduction

The efficient numerical solution of the Poisson equation remains a central challenge across computational physics, with applications spanning electromagnetism, gravitational problems, and, most notably, incompressible fluid dynamics. In many time‑dependent simulations, the Poisson solve constitutes the dominant computational cost, particularly in large‑scale problems where the physical phenomena exhibit widely disparate spatial scales. This contribution presents a novel hybrid multigrid-FFT Poisson solver designed to operate on fully adaptive multiresolution grids, while supporting arbitrary combinations of periodic, bounded, and unbounded boundary conditions. Particular emphasis is placed on unbounded domains, a regime of great interest in external aerodynamics and vortex‑dominated flows, where the solution must decay naturally at infinity.

The proposed solver is implemented within murphy [1], a massively parallel, block‑structured, collocated multiresolution framework. To address the intrinsic challenges posed by unbounded directions, the scheme relies on the FLUPS library [2], an FFT‑based Poisson solver.

2. Methodology

The methodology introduces an adaptive multigrid strategy designed to solve the Poisson equation on three‑dimensional multiresolution grids while supporting high‑order accuracy and unbounded boundary conditions. The approach replaces the classical hierarchy of global uniform meshes with a composite grid built from locally refined regions generated by the murphy multiresolution framework. This structure enables the use of the Full Approximation Scheme (FAS), in which both solution and residual fields are transferred between levels to maintain consistency across resolution jumps and to avoid nonlinear artefacts originating from interpolation‑based ghost reconstructions.

A standard adaptive V‑cycle with FAS is used, combining Gauss-Seidel smoothing, full‑weighting restriction, and trilinear prolongation. At the coarsest level (kept intentionally large to preserve parallel load balancing) the Poisson equation is solved directly using a Fourier‑based solver. The spatial discretization employs compact Mehrstellen stencils of second, fourth, or sixth order. These stencils maintain a radius‑one structure, thus avoiding extra communication overheads, while higher accuracy is enforced through a modified right-hand side constructed via high‑order finite-difference operators. This compact formulation provides high-order convergence without increasing halo sizes.

A key methodological contribution concerns the imposition of unbounded boundary conditions. The coarse‑grid direct solver applies lattice Green’s functions (LGFs) corresponding exactly to the finite-difference operator used throughout the multigrid hierarchy. Ensuring this operator–kernel compatibility is essential, as mismatched LGFs lead to non‑convergent multigrid cycles. For configurations with two unbounded directions, an LGF construction based on combining a known two‑dimensional unbounded kernel with the discrete Laplacian’s Fourier symbol provides a compatible solution. Finally, the methodology addresses the non‑conservative nature of interpolation across resolution jumps, which causes weak violations of global flux conservation and complicates convergence monitoring in singular (Neumann or periodic) Poisson problems. Instead of relying on residual norms, the solver employs a robust Cauchy‑type convergence criterion between successive solution iterates.

3. Results

The proposed multigrid-FFT Poisson solver is validated across periodic, semi‑unbounded, and fully unbounded domains using the method of manufactured solutions. In all configurations, the solver exhibits the expected convergence behavior: the maximum pointwise error decreases proportionally to the refinement tolerance, following the predicted scaling (see [1]), and achieving second-, fourth-, and sixth‑order accuracy depending on the compact stencil employed. For unbounded domains with one, two, or three unbounded directions, convergence is preserved thanks to the consistent pairing of finite‑difference stencils with their corresponding lattice Green’s functions on the coarse grid.

The adaptive multiresolution discretization reveals two distinct error regimes. At coarser resolutions, the method displays spectral‑like convergence due to localized refinement near sharp features. As refinement intensifies, the 2:1 multiresolution constraint leads to domain‑wide refinement and a transition to the polynomial convergence expected of uniform high‑order schemes.

Parallel performance is assessed through a large‑scale Biot-Savart test case involving mixed periodic–unbounded boundary conditions. Weak scaling experiments performed on the LUMI supercomputer demonstrate excellent scalability of all multigrid components, with only the coarse‑grid FFT solver showing reduced efficiency, an anticipated effect given the small problem size per rank at that level. Despite this, the overall solver maintains near‑perfect weak scaling up to 128 nodes (16,384 cores).

4. Conclusion

The presented solver combines an adaptive multiresolution multigrid method with a high‑performance FFT‑based coarse‑grid solver, delivering an efficient and flexible strategy for Poisson problems involving periodic, bounded, or fully unbounded directions. The approach enforces consistency between finite‑difference stencils and lattice Green’s functions, enabling accurate high‑order solutions even in unbounded configurations. Thanks to compact high‑order stencils and a sufficiently large coarse grid, the method avoids the typical communication and load‑balancing bottlenecks of multigrid schemes.

The solver exhibits the expected accuracy and convergence behavior, with error trends showing both spectral‑like and polynomial regimes depending on the refinement level and the structure of fine‑scale features. Large‑scale performance tests demonstrate excellent weak scalability, with only a limited overhead from the FFT coarse solve, confirming that the combined method is well‑suited for massively parallel adaptive simulations. The solver thus provides a robust and broadly applicable tool for computational physics, with future work aiming to relax refinement constraints and extend the implementation to heterogeneous CPU–GPU architectures.

[1] T. Gillis and W. M. van Rees, MURPHY—A Scalable Multiresolution Framework for Scientific Computing on 3D Block-Structured Collocated Grids, SIAM Journal on Scientific Computing, 44 (2022) [2] P. Balty, P. Chatelain, and T. Gillis, FLUPS - A Flexible and Performant Massively Parallel Fourier Transform Library, IEEE Transactions on Parallel and Distributed Systems, 34 (2023)

ABSTRACT. In this talk innovative enhancements to symbol-based multigrid methods for solving large block-structured linear systems are presented. Block-structure linear systems naturally arise, e.g., in the discretization of PDEs using higher order finite elements.

Our work focuses on an aggregation-based grid transfer operators that map the symbol of a block Toeplitz matrix from matrix-valued to scalar-valued at coarser levels. Based on [1] through a convergence analysis of the Two-Grid Method (TGM), in [2] we establish a direct relationship between the properties of the scalar-valued symbol at the coarser level and those of the original matrix-valued symbol. This analysis enables us to demonstrate the convergence of a V-cycle multigrid approach using standard grid transfer operators designed for scalar Toeplitz systems at the coarse levels.

Building on these insights, we broaden the range of effective smoothers for block Toeplitz matrices, emphasizing the performance of block-oriented techniques such as the relaxed block Jacobi method. We derive general conditions for smoothing parameters, highlighting scenarios where these parameters can be computed with minimal computational overhead. Naturally, this works also for other approaches for block Toeplitz systems like the ones proposed in [3].

The proposed methods are validated on linear systems arising from the discretization of differential equations using Qd Lagrangian finite elements or B-splines with non-maximal regularity. Numerical experiments demonstrate significant computational advantages over existing techniques for block-structured linear systems.

References:

[1] C. An and Y. Su. An aggregation-based two-grid method for multilevel block Toeplitz linear systems. J. Sci. Comput., 98(3):47, 2024. Id/No 54. [2] M. Bolten, M. Donatelli, P. Ferrari, and I. Furci. Analysis on aggregation and block smoothers in multigrid methods for block Toeplitz linear systems. SIAM J. Matrix Anal. Appl., 46(2):1416–1443, 2025. [3] M. Bolten, M. Donatelli, P. Ferrari, and I. Furci. A symbol based analysis for multigrid methods for block-circulant and block-Toeplitz systems. SIAM J. Matrix Anal. Appl., 43(1):405–438, 2022.

ABSTRACT. The HYPRE library is a suite of parallel linear solvers and preconditioners featuring multigrid methods. One of the hallmarks of the library is its conceptual linear system interfaces, which allows users to describe their linear systems in the traditional row-column based matrix format, but also in terms of (logically) structured grids and stencils. This semi-structured interface provides information that HYPRE can then exploit to gain performance. Although the semi-structured interface has been available for many years, users only had access to fully structured and fully unstructured multigrid methods, while practical applications often fall somewhere in between. Recently, the semi-structured component of HYPRE was overhauled to extend functionality from square matrices to rectangular matrices with full support for matrix-vector and matrix-matrix operations, enabling the development of semi-structured solvers such as the new semi-structured algebraic multigrid (SSAMG) method. In this talk, we discuss the new developments, especially the stencil algebra underlying them.