ABSTRACT. A polynomial of A can give a very good approximation of the inverse of A. For sparse matrices, this in some sense gives a sparse version of the inverse. For ill-conditioned matrices, a high degree polynomial may be needed. The GMRES algorithm can be used to develop this polynomial.

Application is given to linear equations with multiple right-hand sides. Once a polynomial is developed, it can be applied to each right-hand side. Cost can be reduced by computing some eigenvalues and then generating a deflated polynomial of lower degree.

To reduce orthogonalization costs, we will try both symmetric and nonsymmetric Lanczos for finding the polynomial and try double (composite) polynomials. Other possible topics are the challenging case of indefinite matrices and trace estimation of the inverse.

ABSTRACT. Indefinite spectra occur for linear equations in many applications. Examples include Helmholtz equations such as the wave equation and in quantum chromodynamics. Indefinite problems can be very difficult for Krylov iterative methods. We investigate adding polynomial preconditioning for such problems and show it can give tremendous improvement. Several difficulties can arise in polynomial preconditioning for indefinite matrices. These are mentioned along with some algorithmic solutions. Eigenvalue problems may also be mentioned.

ABSTRACT. Obtaining a suitable, effective preconditioner is often a tall task -- both theoretically and computationally. In many applications, e.g., subsequent linearizations, topology optimization, optimization and others, we face a sequence of linear systems and obtaining a solid preconditioner for each of them separately is computationally prohibitive. Unfortunately, in many of these cases the naive approach of calculating a powerful preconditioner $P_0$ for the first problem $A_0 x_0 = b_0$ leaves a lot to be desired and hence the natural next step is to try to recycle some information from $P_0$. Among the possible approaches is the idea of the so-called preconditioner map, which suggests to adapt the original efficient right-preconditioner $P_0$ to the new system $A_1 x_1 = b_1$, obtaining the new right-preconditioner in the form $P_1P_0$, where $P_1$ (approximately) satisfies the equation $P_1 \approx A_1^{-1}A_0 (*)$ so that we recover $A_1P_1 \approx A_0P_0$ and hence also the good convergence properties of the initial preconditioned system. The key problem of (approximately) solving the equation (*) has been approached by Carr et al. in 2021, using sparse matrix formulations, obtaining the so-called Sparse Approximate Maps for ``updating'' the preconditioner when the system changes.

We propose, analyze and numerically investigate this framework but using hierarchical matrix formulations instead and obtain Hierarchical Approximate Maps for preconditioner updating. The hierarchical formats such as HODLR and HSS matrices have shown to be efficient when used in the right set-up and for some model problems we show that in this context they can be a powerful tool for re-using the computational resources expanded for obtaining the initial preconditioner $P_0$.

ABSTRACT. The recent proliferation of smart devices is driving interest in including such devices in distributing computing. These “edge” devices present issues for standard distributed linear solver algorithms: (i) synchronization points (e.g. reduce operations) in the solvers mean the entire computation is throttled by the slowest edge device and (ii) the edge devices might not maintain the same precision and integrity of the data assumed by the solvers. This work explores how asynchronous versions of standard linear solvers can mitigate those issues, especially when combined with fault tolerance strategies for detecting and rejecting data corruption. A recently developed approach that applies conjugate gradient locally and an s-step approximate conjugate descent algorithm globally seeks to bring some of the beneficial convergence properties of Krylov-based methods to an asynchronous solver. Recently developed rejection criteria based on convergence theory or empirical statistics are observed to be successful in restoring convergence in the presence of at least some data corruption. Ideally, the general approaches presented for asynchronous algorithms and fault tolerance will motivate further innovation in the form of not only additional linear solvers but eigensolvers as well.

This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. Funded by LLNL LDRD projects 21-FS-007 and 22-ERD-045. LLNL-ABS-859954.

ABSTRACT. We consider solving many block linear systems of the form (A + z_i I)X = B_i, where A is a fixed Hermitian matrix and each z_i is a complex constant. GMRES can solve each individual linear system but lacks a short-term recurrence. We derive a short-term recurrence block Krylov subspace method which efficiently solves each of these systems in sequence. During each linear solve, we converge certain eigenvectors and eigenvalues of the matrix A. These eigenvectors are reused to construct increasingly better initial guesses for subsequent systems. One specific formulation of the random-phase approximation in density functional theory (DFT) necessitates the solution of many of these block systems. Our method achieves fast solution times compared to GMRES and is well-suited for large-scale DFT calculations.

ABSTRACT. Inverse problems can be found in a variety of scientific and engineering applications and involve the reconstruction of unobserved objects from possibly noisy data. This kind of problems are challenging to solve since they are typically ill-posed: the reconstruction is very sensitive to perturbations in the measurements; and real-world applications are often large-scale: resulting in computationally demanding tasks. In the cases where the matrix is not explicitly stored, and where the matrix structure does not allow for closed form approximations of the problem solution, efficient and functional algorithms need to carefully leverage the use of iterative solvers and suitable regularization.

In this talk I will focus on inner-product free Krylov solvers, offering efficient and highly parallelizable alternatives to traditional Krylov methods based on partial re-orthogonalizations, which are particularly suited for implementations that exploit randomization and mixed precision arithmetic. In particular, I will explain how two solvers based on the Hessenberg method perform when solving large-scale ill-posed problems, paying special attention to their iterative regularizing properties as well as presenting a new hybrid variant that incorporates Tikhonov regularization on the projected problems. The two presented methods are called changing minimal residual Hessenberg method (CMRH) and hybrid changing minimal residual Hessenberg method (H-CMRH).

To illustrate the theoretical findings of this work, I will present extensive numerical experiments that highlight the regularizing properties of the discussed methods.

ABSTRACT. l1 regularization is used to preserve edges or enforce sparsity in a solution to an inverse problem. We investigate the Split Bregman (SB) and the Majorization-Minimization (MM) iterative methods that turn this non-smooth minimization problem into a sequence of steps that include solving an l2-regularized minimization problem. To select the regularization parameter at each inner iteration, we extend the generalized cross validation (GCV) selection method to the corresponding inner problem. We present numerical results on signal reconstruction and image deblurring problems, which show that SB-GCV performs well while MM-GCV performs somewhat poorly in comparison.

ABSTRACT. Uncertainty quantification for large-scale inverse problems remains a challenging task. For linear inverse problems with additive Gaussian noise and Gaussian priors, the posterior is Gaussian but sampling can be challenging, especially for problems with a very large number of unknown parameters (e.g., dynamic inverse problems) and for problems where computation of the square root and inverse of the prior covariance matrix are not feasible. Moreover, for hierarchical problems where several hyperparameters that define the prior and the noise model must be estimated from the data, the posterior distribution may no longer be Gaussian, even if the forward operator is linear. Performing large-scale uncertainty quantification for these hierarchical settings requires new computational techniques.

In this work, we consider generalized Golub-Kahan based methods for large-scale, hierarchical Bayesian inverse problems. We consider a hierarchical Bayesian framework and exploit generalized Golub-Kahan based methods to efficiently sample from the posterior distribution. Numerical examples from dynamic photoacoustic tomography and atmospheric inverse modeling demonstrate the effectiveness of the described approaches.

ABSTRACT. X-ray crystallography has been the workhorse technique for solving the structure of macromolecules, and relies on the analysis of X-ray diffraction data from crystalline samples. An ideal, well-ordered crystal consists of a periodic repeat of exactly identical copies of a molecule, resulting in a diffraction pattern that consists of a set of discrete and sharp Bragg peaks. Biological molecules, however, are not static entities, but exist in an ensemble of conformations (shapes), even when pinned to a crystalline lattice. X-ray diffraction from a crystal with copies of the molecule in different conformations consists of Bragg peaks as well as a smoothly varying signal called diffuse scattering. However, traditional crystallographic structure solution methods are focused on handling Bragg data only yielding an ensemble averaged structure, and are unable to reconstruct biologically relevant motions from the diffuse signal.

Here we present a new algorithmic approach based on iterative projection algorithms (IPAs) to extract functional dynamics of macromolecules by reconstructing multiple macromolecular conformations from Bragg and diffuse signal within a self-consistent framework. IPAs aim to solve an inverse problem by posing it as a constraint satisfaction problem by utilizing norm-minimizing projection operators to find the intersection of the constraints, and have been successfully employed to recover phase information in coherent diffractive imaging, electron microscopy, and astronomical imaging. We propose to approach diffuse X-ray scattering as an inverse problem of deriving the models of correlated motion by iterative application of constraints to data consisting of both Bragg peaks and the diffuse signal. We formulate a diffraction model that accounts for the presence of multiple conformations and defects in terms of an occupancy matrix and exploit the properties of the occupancy transform to develop projection operators that satisfy constraints of periodicity, symmetry, and support in the real and reciprocal space to recover the density of different protein conformations that are consistent with the Bragg and diffuse signal. We will show results that demonstrate that our approach works on crystals with conformational disorder, vacancies and stacking faults, and edge effects.

ABSTRACT. The Anderson Acceleration (AA) is a well-established method that allows to speed up or encourage convergence of fixed-point iterations. It has been successfully used in a variety of applications, in particular within the Self-Consistent Field (SCF) iteration method for quantum chemistry and physics computations. In recent years, the Conjugate Residual with OPtimal trial vectors (CROP) algorithm was introduced and shown to have a better performance than the classical Anderson Acceleration with less storage needed. This paper aims to delve into the intricate connections between the classical Anderson Acceleration method and the CROP algorithm. Our objectives include a comprehensive study of their convergence properties, explaining the underlying relationships, and substantiating our findings through some numerical examples. Through this exploration, we contribute valuable insights that can enhance the understanding and application of the mixing acceleration methods in practical computations, as well as the developments of new and more efficient acceleration schemes.

ABSTRACT. Anderson acceleration (AA) is an extrapolation technique used to accelerate the convergence of fixed-point iterations $x_{k+1}=q(x_k)$, $x_k\in\R^n$, but little is known concerning the convergence improvement provided by AA. In this work, an asymptotic convergence analysis is presented for two classes of linear fixed-point maps $q(x) = Mx+b$, where $M\in\R^{2\times2}$ and $x,b\in\R^2$. The extrapolation coefficients that arise from applying AA to these maps are shown to follow a periodic orbit, with no acceleration provided in every other iteration. A connection is established between AA applied to these maps and a restarted variant of AA with an accompanying example comparing the two methods. Additionally, we show that the convergence of AA for symmetric $M$ can be simplified to that of diagonal matrices $M$. The asymptotic convergence results shown in this paper provide a partial understanding of the acceleration property of AA with some numerical tests that indicate the performance of AA for $2\times2$ linear problems.

ABSTRACT. In this talk, we consider Anderson acceleration with window size $m$ (AA(m)) applied to fixed-point iteration for linear systems. We explore some conditions on the $m+1$ initial guesses of AA(m), aiming for the residuals $r_{m+1}=0$. We propose the sufficient and necessary conditions on the $m+1$ initial guesses for $r_{m+1}=0$. These findings can help us better understand the performance between original fixed-point iteration and Anderson acceleration. Meanwhile, it may give us some guidance on the choice of good initial guesses. Moreover, we give examples to show the temporary stalling behavior of Anderson acceleration applied to solving linear systems.

ABSTRACT. Anderson Acceleration (AA) has been widely used to solve nonlinear fixed-point problems because of its rapid convergence characteristics. This work considers a variant where Anderson Acceleration and Picard iteration are alternated, known as the Alternating Anderson-Picard (AAP) method. This method capitalizes on the efficiency of AA in speeding up convergence and the simplicity of Picard iterations to generate historical points at each AA step. We present the theoretical analysis of the convergence properties of AAP. Theoretical analysis indicates that the AA-Picard method is equivalent to a Jacobian-free Newton–Krylov method, leading to a local linear or quadratic convergence rate towards the solution. This equivalence also helps to bound the optimization gain, a concept crucial in quantifying AA's convergence rate, yet its implications have not been fully understood.

ABSTRACT. This talk describes recently-developed algorithms for computing data-sparse representations of rank-structured matrices, specifically H2 matrices. The algorithms are black-box, meaning that they only interact with the matrix to be compressed through its action on vectors, making them useful for tasks like forming Schur complements or matrix-matrix multiplication.

Of particular interest is the algorithm for "Randomized Strong Recursive Skeletonization", or RSRS, that simultaneously both compresses and factorizes an H2-matrix with a strong admissibility criterion in the black-box setting. The work has particular relevance for 3D solvers and preconditioners.

ABSTRACT. In this talk, methods are developed to iteratively embed aggregate pieces of information from a data matrix to generate low-rank approximations. The algorithms are built by sampling the row and column spaces of the matrix simultaneously. The developed two-sided randomized algorithms converge with provable error bounds. A heuristic accelerated scheme is also developed, motivated by the two-sided sampling analysis. We compare our algorithms to one-sided sampled algorithms on a selection of matrices from a standard test-suite and verify that the theoretical convergence rates are numerically realized.

ABSTRACT. In this talk we present an adaptive, partially matrix-free, Hierarchically Semi-Separable matrix construction algorithm that leverages the sparse Johnson-Lindenstrauss transform (SJLT). This algorithm is an extension of Gorman et al. [SIAM J. Sci. Comput. 41(5), 2019] which leveraged Gaussian sketching operators. We first present theoretical work which justifies this extension via concentration bounds. Then, we discuss the implementation details of applying SJLT efficiently. Finally, we demonstrate experimentally that using SJLT instead of Gaussian sketching operators leads to 1.5-2.5x speedups of the HSS construction implementation in the STRUMPACK C++ library while maintaining similar relative error.

ABSTRACT. Randomized or deterministic sketching methods have been effective for the solution of consistent rectangular or square linear systems. In particular, when the matrix is large and the data contains redundancies or errors, highly accurate solutions might not always be desired. Instead of approximating a solution using all the system information, sketching methods use smaller, projected systems. When randomized sketches are used to generate the smaller subproblems, the convergence of the resulting methods is determined by rates that depend on random matrices. When the rate is close to unity, the observed convergence can slow on practical problems. In this work we develop structured sketching methods, that contain an implicit preconditioning mechanism based on known properties of the system. For instance, when the matrix is symmetric indefinite the new methods improve in the number of iterations when compared to state-of-the-art random or deterministic sketching methods.

ABSTRACT. In this talk we present new multigrid reduction-in-time (MGRIT) parallel-in-time methods that solve nonlinear hyperbolic PDEs and hyperbolic systems in a small number of iterations with convergence rate independent of mesh resolution.

In the first part of the talk, we consider scalar nonlinear conservation laws in one spatial dimension. The equations are discretized in space with a conservative finite-volume method using weighted essentially non-oscillatory (WENO) reconstructions, and in time with high-order explicit Runge-Kutta methods. The solution of the global, discretized space-time problem is sought via a nonlinear iteration that uses a novel linearization strategy in the case of non-differentiable discrete equations. At each nonlinear iteration, the linearized problem takes the form of a certain discretization of a linear conservation law over the space-time domain in question. An approximate parallel-in-time solution of the linearized problem is computed with a single MGRIT iteration. The MGRIT iteration employs a novel coarse-grid operator that is a modified conservative semi-Lagrangian discretization and generalizes those we have developed previously for non-conservative scalar linear hyperbolic problems. Numerical tests are performed for the inviscid Burgers and Buckley-Leverett equations. For many test problems, the solver converges in just a handful of iterations with convergence rate independent of mesh resolution, including problems with (interacting) shocks and rarefactions.

In the second part of the talk, we extend the approach to systems of hyperbolic PDEs. We first consider the case of the linear acoustic equations in heterogeneous media. We obtain fast convergence in a small number of iterations with convergence rate independent of mesh resolution, including for problems with discontinuous material parameters. Finally, we show some preliminary results extending the approach to nonlinear hyperbolic systems such as the shallow water equations.

ABSTRACT. In this talk we discuss the application of waveform relaxation (via geometric multigrid) to (2+1)D space-time discretisations of Navier-Stokes. We compare our method against direct linear solvers in addition to investigating the effectivity of the time parallelism of the proposed algorithm.

ABSTRACT. Multigrid Reduction in Time (MGRIT) uses multigrid reduction techniques to enable temporal parallelism for solving initial value problems. It is known that the convergence rate of MGRIT depends in part on the choice of time-stepping operators on the fine- and coarse-grid, which are termed the fine-grid operator and coarse-grid operator respectively. This work explores algebraic approximations to the ideal coarse-grid operator for the linear advection equation, improving the convergence when standard re-discretized coarse-grid operators are used.

ABSTRACT. Despite the fact that Parallel-in-Time (PinT) methods are predicted to become necessary to fully utilize next-generation exascale machines, there are currently no known practical methods which scale well with the length of the time-domain for chaotic problems, due to exponential dependence of the condition number on the fastest chaotic time-scale. While most prior works applying multigrid-in-time to chaotic systems focus on the coarse grid equation, we investigate the relaxation techniques commonly used and demonstrate that they in fact diverge for chaotic systems. Here the novel Least Squares Relaxation (LSR) is presented and proven to be a convergent, PinT relaxation for chaotic PDE systems. Promising preliminary analytical results and numerical experiments with the Lorenz system indicate that LSR may solve the scaling problem for chaotic systems, potentially allowing massive space-time parallelization of turbulent computational fluid dynamics.

ABSTRACT. Traditionally, algebraic multigrid (AMG) research assumed a fixed pointwise smoother. For some problems such as the Maxwell equations, this is not sufficient. Block smoothers are a popular choice in geometric multigrid methods. They are sometimes required to achieve O(n) multigrid convergence. However, the research on using block smoothers on AMG is still in the early stages. We would like to combine the best of both worlds and build a generalized AMG theory with block smoothers. We recently developed an automatic approach for building smoothers in AMG based on the ideal interpolation that is designed to eliminate near-null space components that have local support.

ABSTRACT. Partitioned methods for coupled problems rely on data transfers between subdomains to synchronize the sub- domain equations and enable their independent solution. By treating each subproblem as a separate entity, these methods enable code reuse, increase concurrency and provide a convenient framework for plug-and-play multiphysics simulations. However, accuracy and stability of partitioned methods depends critically on the type of information exchanged between the subproblems. The exchange mechanisms can vary from minimally intrusive remap across interfaces to more accurate but also more intrusive and expensive estimates of the necessary information based on monolithic formulations of the coupled system. These transfer mechanisms are separated by accuracy, performance and intrusiveness gaps that tend to limit the scope of the resulting partitioned methods to specific simulation scenarios. Data-driven system identification techniques provide an opportunity to close these gaps by enabling the construction of accurate, computationally efficient and minimally intrusive data transfer surrogates. This approach shifts the principal computational burden to an offline phase, leaving the application of the surrogate as the sole additional cost during the online simulation phase. In this paper we formulate and demonstrate such a surrogate-based partitioned method for a model advection-diffusion transmission problem by using Dynamic Mode Decomposition (DMD) to learn the dynamics of the interface flux from data. The accuracy of the resulting DMD flux surrogate is comparable to that of a dual Schur complement reconstruction, yet its application cost is significantly lower. Numerical results confirm the attractive properties of the new partitioned approach.

ABSTRACT. The atmosphere and ocean interact through the exchange of surface fluxes which, in global climate models such as E3SM, are derived from a system of nonlinear equations. Currently, these surface fluxes are obtained via a simple fixed point iteration. In this talk, we analyze the existence and uniqueness of solutions of these equations from a dynamical systems perspective. In particular, an analysis of the bifurcation points of the system demonstrates that there are conditions under which these equations have no solutions and conditions under which several solutions exist. We discuss developments towards recent improvements in the fixed point iteration such as regularization techniques to ensure existence and uniqueness of the surface fluxes.

This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. LLNL-ABS-859925-DRAFT.

ABSTRACT. Creating a digital twin for a physical system often requires multi-scale, multi-physics, and multi-fidelity modeling. However, simulation of individual components of a system are often done using specialized software designed with a particular set of physics in mind and need to be coupled together for multiple components. Additionally, either through the identification of candidates for reduced order modeling (ROM) or machine learned models (ML), or in the absence of a physical model such as a data-driven surrogate, it then becomes necessary to couple non-traditional numerical methods together or with other traditional numerical methods such as the finite element method (FEM).

In this talk we present work on a new method for coupling FEM-to-FEM, FEM-to-ROM, and ROM-to-ROM methods using a Lagrange Multiplier (LM) along with a dual-Schur complement. The discussion will cover formulation of the coupling scheme as well as details regarding the construction of LM spaces that lead to a non-singular Schur complement system. Numerical results will be presented demonstrating the accuracy, stability, and conditioning of the approach.

ABSTRACT. Problems involving coupling over a nonoverlapping interface abound, with common examples including virtual interfaces introduced for the purpose of domain decomposition and physical interfaces such as are found in, e.g., multimaterial solid interaction and fluid-structure interaction. Solution of coupled problems can be particularly difficult, as a monolithic fully-coupled system of equations is challenging to develop software for, requiring the combination of multiple areas of domain expertise. They can be even more challenging to solve once posed, due to difficulty in preconditioning the resulting large linear systems. Partitioned approaches for solution of coupled problems enable the use of subdomain solvers developed by domain experts and are more amenable to a ``plug-and-play'' style of framework. However, with partitioned approaches, considerable attention must be given to enforcing coupling conditions. As model reduction techniques continue to mature, it is natural that there is a desire among analysts to use cheaper, more efficient models where possible. This fits closely with the partitioned approach to solving coupled problems. It is towards this goal of solving multifidelity coupled problems involving reduced order models (ROMs) and full order models (FOMs) that we extend a partial differential equation (PDE) constrained optimization technique, developed for FOMs by M. Gunzberger and H. K. Lee, to ROM-ROM coupling. The use of a ROM in one or more subdomains allows for harnessing their efficiency, but also presents a challenge with respect to reduced adjoint systems required to be solved as part of the technique. We explore the ROM-ROM coupling technique, applying it to coupled advection-diffusion equations. We present numerical studies demonstrating the accuracy of the approach along with an investigation related to selecting a reduced order basis for the adjoint systems.

ABSTRACT. Coupled Earth System Model (ESM) simulations involve computationally complex models of the atmosphere, ocean, land surface sea ice, and other components that must act as one system, increasing the complexity. Rigorous spatial coupling between components in ESMs involves field transformations, and communication of data across multiresolution grids while preserving key attributes of interest such as global integrals and local features. The Energy Exascale Earth System Model (E3SM) is one such solver with a strong need to balance numerics and computational performance for simulating climate response on heterogeneous computing architectures.

Significant effort has been invested in moving away from the current Model Coupling Toolkit (MCT) hub-and-spoke coupler in E3SM, which is unaware of the underlying mesh topology but relies on trivial decompositions of linear maps loaded from disk. Unlike MCT, the Mesh Oriented datABase (MOAB) allows a complete description of the topology used by each submodel in E3SM, along with serialization of the discrete coupled scalar and flux solution fields participating in the nonlinear coupling. Such an approach provides the flexibility to use arbitrary grids and adaptive decompositions while requiring mesh intersections, conservative projection operator computation, and remapping weight application to be seamlessly available to transfer coupled field data between component models.

With constraints to recover full reproducibility and verifiability compared to the MCT implementation, and to maximize the computational efficiency without sacrificing discretization accuracy of field data being transferred, we present comparative accuracy and performance results between the two approaches for various high-fidelity cases on large-scale machines. We also discuss future extensions to the coupler to ensure property preservation with high-order projection approaches and novel distributed coupling interfaces to depart from the traditional hub-and-spoke coupler model.