ABSTRACT. We present the construction of the multiscale two-grid preconditioner for diffusion processes. The method is based on the careful design of accurate coarse space approximation by solving local spectral problems for normalized graph Laplacians in overlapped local domains. The choice of smoothing iterations is investigated. We apply the proposed algorithm for solving several applied problems: (a) mixed-dimensional flow model in fractured porous media; (b) anisotropic heat flow with a high level of anisotropy; (c) transport problems in the least square formulation; and (d) discrete graph models in high-contrast networks.

ABSTRACT. We introduce one Gauss-Seidel sweep as an efficient solver for the Gram matrix systems arising in the formulation of s-step Conjugate Gradient methods (CG). Gauss-Seidel in this context is equivalent to a modified Gram-Schmidt projection in the A-inner product x^TAx and is based on a matrix splitting $P^TAP = I + L + L^T$ after diagonal scaling. This approach reduces computational and communication costs while maintaining numerical stability, making it well suited for modern high-performance computing (HPC) architectures b GPU accelerators and APUs such as the AMD MI-300.

The Gauss-Seidel solver effectively handles ill-conditioned matrices, achieving robust convergence with minimal synchronization overhead as block sizes grow. Rapid convergence is explained by way of a backward error analysis based upon Neuman series. Extending this approach to AMG, we tackle the challenge of approximate solves for the coarse-level problem.

The method integrates seamlessly into mixed-precision multi-grid methods, enhancing scalability and robustness. Numerical experiments demonstrate that stagnation is avoided while maintaining low Krylov solver iteration counts with 4X faster execution on the MI-300 GPU, and superior parallel scalability in multilevel solvers for large-scale simulations. These findings underscore the potential of Gauss-Seidel for overcoming computational bottlenecks in s-step CG.

ABSTRACT. In this work, we introduce a new approach for efficiently solving large, sparse linear systems arising from the discretization of elliptic partial differential equations using isogeometric analysis. Our method integrates vector extrapolation techniques with geometric multigrid solvers, significantly enhancing the convergence of standard multigrid iterations with classical smoothers. We further extend this approach to nonlinear problems by employing polynomial-type extrapolation methods to improve the convergence of fixed-point methods, such as the Picard iterative method. This work provides quadratic convergence results for polynomial extrapolation methods. Specifically, a new theoretical result on the correlation between the residual norm and the error norm, as well as a new estimation for the generalized residual norm of some extrapolation methods, are given. Through various numerical experiments, we demonstrate that the proposed combination of multigrid solvers with polynomial extrapolation methods considerably accelerates convergence and provides a robust alternative to existing techniques, such as Anderson acceleration, for solving nonlinear problems efficiently.

ABSTRACT. The Discrete Exterior Calculus (DEC) defines a family of discretized differential operators which preserves certain desirable properties from the Exterior Calculus. CombinatorialSpaces.jl is an open-source Julia library which implements the DEC on a wide range of simplicial complexes, and now offers a geometric multigrid solver over maps between subdivided simplicial complexes. We demonstrate numerical results of multigrid solvers for the Laplacian operator, both as a standalone solver and as a preconditioner for open-source Julia iterative methods libraries.

ABSTRACT. Modelling turbulence around wind turbines can be done with the Reynolds-averaged Navier-Stokes equations. The challenges faced when solving such problems arise from two main sources. The first main challenge is the need to model the Reynolds-stress tensor. There are a wide range of approaches to this problem, often involving eddy-viscosity or closure models. Another source of difficulty when solving such problems is due to the very high Reynolds numbers of such systems, often in the thousands.

In this talk we focus on the aspect of high Reynolds numbers ignoring the Reynolds-stress tensor. We study an attempt at solving such a system with very high Reynolds numbers by combining approaches from genetic programming and continuation methods, with a particular emphasis on using (AMG +) RAS/ILU as part of the solver.

ABSTRACT. Recently we developed the damped block Newton (dBN) method, along with adaptive dBN (AdBN), for the shallow Ritz approximation. The dBN method alternates between updates of the linear (weights) and non-linear parameters (biases). Per each iteration, the linear and the non-linear parameters are updated by exact inversion and one step of a modified, damped Newton method applied to a reduced non-linear system, respectively. The geometric understanding of the parameters has proven to be crucial in implementing an effective solver, particularly for the non-linear parameters. Furthermore, the geometric meaning of the parameters gives insight into constructing cheap error indicators, further improving the already improved efficiency of AdBN.

ABSTRACT. This talk presents a local convergence analysis for block Newton (BN) methods, which are used to find the best neural network approximation for one-dimensional least-squares and elliptic PDE problems. These methods use the block Gauss-Seidel method as an outer iteration between the linear and nonlinear parameters. In each outer iteration, both the linear (weights) and nonlinear parameters (biases) are updated with one step of Newton's method. This analysis frames BN methods as fixed-point iteration methods. The conditions for convergence are derived for various applications (least-squares, diffusion, and diffusion-reaction problems).

ABSTRACT. Machine learning (ML) has been demonstrated as a powerful tool across numerous fields of research. Despite its efficacy, training ML models, particularly large ones, remains challenging and requires significant computing resources and long training times. In numerical PDEs, multilevel methods were one of the most transformative algorithms developed, solving problems of n unknowns that would cost O(n^3) FLOPs with direct solvers or O(n^2) with, e.g., Krylov methods, to O(n) or O(n\log n) using multilevel methods. Recent years have seen increasing interest in applying multilevel concepts to the training of ML models, but, despite some progress, robust and significant speedups have remained elusive. In this talk, I will review the role of preconditioning in the optimization of ML models, then proceed to discuss the challenges of incorporating multilevel ideas, and directions that may have potential.

ABSTRACT. Convolutional Neural Networks (CNNs) are the backbone of many deep learning methods, but optimizing them remains computationally expensive. To address this, we explore multiscale training frameworks and mathematically identify key challenges, particularly when dealing with noisy inputs. Our analysis reveals that in the presence of noise, the gradient of standard CNNs in multiscale training may fail to converge as the mesh-size approaches to , undermining the optimization process. This insight drives the development of Mesh-Free Convolutions (MFCs), which are independent of input scale and avoid the pitfalls of traditional convolution kernels. We demonstrate that MFCs, with their robust gradient behavior, ensure convergence even with noisy inputs, enabling more efficient neural network optimization in multiscale settings. To validate the generality and effectiveness of our multiscale training approach, we show that (i) MFCs can theoretically deliver substantial computational speedups without sacrificing performance in practice, and (ii) standard convolutions benefit from our multiscale training framework in practice.

ABSTRACT. Graph Neural Networks (GNNs) have emerged as a powerful tool for learning and inference on graph-structured data, and are wildly used in a variety of applications, where this usage often takes into considerations a large amount of data. However, training on large graphs is an inefficient and expensive process, both in terms of memory and computations. In this paper, we introduce a novel framework for multiscale training of GNNs, designed to efficiently integrate information across multiple structural levels of a graph. Our approach leverages hierarchical graph representation, taking a smaller amount of nodes into account and use the advantages of coarse graph scales in the training process. We propose several multiscale method for training GNNs and demonstrate our approach on various datasets and learning problems.

ABSTRACT. We present a novel block-preconditioner for the elastic Helmholtz equation, based on a reduction to the acoustic Helmholtz equation. Utilizing approximate commutativity of the underlying differential operators, we suggest a block-triangular preconditioner whose diagonal blocks are acoustic Helmholtz operators. Thus, we enable the solution of the elastic version using virtually any existing solver for the acoustic version as a black box. We prove a sufficient condition for the convergence of our method, that sheds light on the long-questioned role of the commutator in the convergence of approximate commutator preconditioners. We show scalability properties of our preconditioner. Our approach, combined with multigrid solve of each block, performs better than an existing monolithic multigrid approach in terms of time and memory, enabling the solution of wave propagation problems in challenging heterogeneous media in 3D on a laptop.

ABSTRACT. As the discretization error for the solution of a partial differential equation (PDE) decreases, the precision required to store the corresponding coefficients naturally increases. The storage complexity of state-of-the-art solvers therefore grows as O(n log n) where n is the number of degrees of freedom (DoFs). This paper presents a full multigrid method to compute the solution in a compressed format reducing the storage complexity of the solution and intermediate vectors to O(n) bits. This allows a matrix-free implementation to solve elliptic PDEs with an overall linear space complexity. For problems limited by the memory capacity of current supercomputers, we expect a memory footprint reduction of about an order of magnitude compared to state-of-the-art mixed-precision methods. Using linear elements to solve Poisson's equation, irrespective of the problem size, as few as 4 bits for the solution and 7 bits for each of three intermediate vectors are required per DoF on the finest grid, in addition to quantities on the coarser grids and lower-order terms.

ABSTRACT. For domains with sufficient structure and regularity, geometric multigrid solvers are among the fastest for computing the numerical solution to elliptic PDEs; however, for complicated domains, constructing a suitable hierarchy of meshes becomes challenging. We propose a framework for mapping computations from such complex domains to a regular computational domain via diffeomorphism, enabling the use of black-box multigrid. This mapping facilitates regular memory accesses during solves, improving efficiency and scalability, especially on massively parallel processors such as GPUs. Moreover, we show that the diffeomorphic mapping itself may be approximately learned using an invertible neural network, enabling automated application to geometries where no analytic mapping exists.