Solving linear systems of the form \((A + \gamma UU^T) x = b\)

I will discuss the iterative solution of large linear systems of equations in which the coefficient matrix is the sum of two terms, a sparse matrix \(A\) and a possibly dense, rank deficient matrix of the form \(\gamma UU^T\), where \(\gamma > 0\) is a parameter which in some applications may be taken to be 1. The matrix \(A\) itself can be singular, but I assume that the symmetric part of \(A\) is positive semidefinite and that \(A+\gamma UU^T\) is nonsingular. Linear systems of this form arise frequently in fields like optimization, fluid mechanics, computational statistics, finance, and others. I will describe a promising iterative strategy for the solution of such linear systems. The performance of the proposed approach is demonstrated by means of numerical experiments on linear systems from different application areas.

This is joint work with Chiara Faccio (Scuola Normale Superiore).

A New Stabilization Algorithm for Cut-cell Meshes

Embedded boundary meshes greatly ease the mesh generation problem for complex geometries.  But when used with explicit finite volume methods, they lead to the ' small cell' problem: cut cells are not stable when used with the time step needed for the regular uncut cells.  We present a new algorithm called State Redistribution that is easy to implement as a post-processing step using mesh information already available in many codes. We show examples of smooth and shocked flows demonstrating its performance.

Krylov Subspace Recycling for Evolving Structures

Krylov subspace recycling is a powerful tool when solving a long series of large, sparse linear systems that change only slowly over time. In PDE-constrained shape optimization, these series appear naturally, as typically hundreds or thousands of optimization steps are needed with only small changes in the geometry. In this setting, however, applying Krylov subspace recycling can be a difficult task. As the geometry evolves, in general, so does the finite element mesh defined on or representing this geometry, including the numbers of nodes and elements and element connectivity.

We developed a method to map the approximate invariant subspaces to be recycled from one mesh to another. The method is very effective for shape optimization and should be applicable in rather general settings. The method will be explained and numerical results will be presented.

Finite Element Methods for Elliptic Problems with Rough Coefficients

We will present finite element methods for elliptic problems with rough coefficients, where the only assumption is on the lower and upper bounds of the eigenvalues of the diffusion matrix. These finite element methods are variants of the localized orthogonal decomposition method in numerical homogenization. Their construction and analysis only involve basic results in finite element methods, domain decomposition methods, and numerical linear algebra.

This is joint work with Jose Garay and Li-yeng Sung

Spatially distributed models of brain metabolism: where biochemistry meets applied mathematics

Most computational challenges in modeling the complex interplay between electrophysiology, hemodynamics and metabolism in human brain come from their multiscale nature, both in time and space. Diffusion of electrolytes and metabolites through extracellular space and glial syncytia, and the neurotransmitter dynamics play a very important role in understanding the coordination of these three different aspects of brain. In this talk we highlight how different mathematical ideas were combine to design a spatially distributed computational model of brain multiphysiology.

Asynchronous Iterative Methods

Daniel B. Szyld, with his collaborators, started working on asynchronous iterative methods in the 1990s, and his work includes an influential survey paper with Andreas Frommer. This talk will describe some recent research on asynchronous iterative methods done in collaboration with Daniel, including asynchronous domain decomposition methods. The talk will be introductory in nature, emphasizing some of the unintuitive behavior of asynchronous iterative methods that currently lack mathematical understanding.

Flexible infinite GMRES for parameterized linear systems

We consider parameterized linear systems \(A(\mu) x(\mu) = b\) for many different \(\mu\), where \(A(\mu)\) is large, sparse and nonsingular, with a nonlinear analytic dependence on \(\mu\). In this work we propose to compute a partial parameterization \(\tilde{x} \approx x(\mu)\), where \(\tilde{x}(\mu)\) is cheap to evaluate for many values of \(\mu\). Our method is based on a companion linearization without truncation, where \(\mu\) appears only linearly. This allows us to combine the well-established Krylov subspace method for linear systems, preconditioned GMRES, with algorithms for nonlinear eigenvalue problems (NEPs) to generate a basis for the Krylov subspace. Moreover, in this work we consider a flexible framework for GMRES, allowing for an inexact iteration dependent application of the preconditioner. This novel approach leads to a significant improvement in complexity over the method infinite GMRES, as proposed in [Jarlebring, Correnty 2021]. The method is analyzed analogously to the standard convergence theory for the method GMRES for linear systems. More specifically, the error is estimated based on the magnitude of the parameter \(\mu\), the inexactness of the preconditioning and the spectrum of the linear companion matrix, which corresponds to the reciprocal solutions to the corresponding NEP. The competitiveness of the method is illustrated with large-scale problems arising from a finite element discretization of a Helmholtz equation with parameterized material coefficient.

Diagonal scalings for improving the accuracy of computed eigenvalues of arbitrary pencils

We present algorithms to construct diagonal scalings for arbitrary matrix pencils \(\lambda B-A\), in which both \(A\) and \(B\) are complex matrices (square or nonsquare). The goal of such diagonal scalings is to "balance'' in some sense the row and column norms of the pencil without changing its eigenvalues in such a way that they are computed more accurately. We will see that the problem of scaling a matrix pencil is equivalent to the problem of scaling the row and column sums of a particular nonnegative matrix. However, it is known that there exist square and nonsquare nonnegative matrices that cannot be scaled arbitrarily. To address this issue, we consider an approximate embedded problem, in which the corresponding nonnegative matrix is square and can always be scaled. The new scaling methods are then based on the Sinkhorn-Knopp algorithm for scaling a square nonnegative matrix with total support to be doubly stochastic or on a variant of it. We illustrate numerically that the new scaling techniques for pencils improve the accuracy of the computation of their eigenvalues. In contrast with previous approaches, which are only valid for regular pencils, the new one works for arbitrary pencils.

Optimal order multigrid preconditioners for the distributed control of parabolic equations with coarsening in space and time

We devise multigrid preconditioners for linear-quadratic space-time distributed parabolic optimal control problems. While our method is rooted in earlier work on elliptic control, the temporal dimension presents new challenges in terms of algorithm design and quality. Our primary focus is on the \(cG(s)dG(r)\) discretizations which are based on functions that are continuous in space and discontinuous in time, but our technique is applicable to various other space-time finite element discretizations. We construct and analyse two kinds of multigrid preconditioners: the first is based on full coarsening in space and time, while the second is based on semi-coarsening in space only. Our analysis, in conjunction with numerical experiments, shows that both preconditioners are of optimal order with respect to the discretization in case of \(cG(1)dG(r)\) for \(r = 0, 1\) and exhibit a suboptimal behaviour in time for Crank–Nicolson. We also show that, under certain conditions, the preconditioner using full space-time coarsening is more efficient than the one involving semi-coarsening in space, a phenomenon that has not been observed previously. Our numerical results confirm the theoretical findings.

Using GMRES to Precondition Arnoldi

Daniel Szyld has been a great champion of Krylov subspace methods, enriching our understanding of their behavior and discovering ways to accelerate their performance. Inspired by his spirit, we will describe how the GMRES residual polynomial can be deployed as a preconditioner for Arnoldi eigenvalue computations. This approach, which is rich in matrix-vector products, can provide an appealing alternative to shift-invert operations, provided one attends to some important stability considerations. (This talk describes joint work with Jennifer Loe and Ron Morgan.)

Adventures with helices: swimming and dissolving at the microscale

Microorganisms often navigate a complex environment composed of a viscous fluid with suspended microstructures such as elastic polymers and filamentous networks. These microstructures can have similar length scales to the microorganisms, leading to complex swimming dynamics. Some microorganisms secrete enzymes that dynamically change the elastic properties of the viscoelastic networks they move through. In addition to nature's organisms, helical microrobots have been engineered with the goals of mucin gel penetration or dissolving blood clots. Here we discuss a model based upon a regularized Stokeslet boundary element method that captures the motion of a helical microswimmer as it penetrates, remodels, and moves through a viscoelastic network. We know that these models would be greatly enhanced by more efficient and robust methods in numerical linear algebra, and hope to spark the interests of this audience full of experts!

Stochastic trace estimation

Computing the trace of a matrix function \(f(A)\), in particular of the matrix inverse, is computationally very demanding when \(A\) is large and sparse. We review stochastic and probing methods for the computation of the trace and present new results based on multilevel Monte Carlo approaches relying on a multigrid hierarchy.

Multiscale Finite Element Methods for an Elliptic Optimal Control Problem with Rough Coefficients.

The solution of multiscale elliptic problems with non-separable scales and high contrast in the coefficients by standard Finite Element Methods (FEM) is typically prohibitively expensive since it requires the resolution of all characteristic lengths to produce an accurate solution. Numerical homogenization methods such as Localized Orthogonal Decomposition (LOD) methods provide access to feasible and reliable simulations of such multiscale problems. These methods are based on the idea of a generalized finite element space where the generalized basis functions are obtained by modifying standard coarse FEM basis functions to incorporate relevant microscopic information in a computationally feasible procedure. Using this enhanced basis one can solve a much smaller problem to produce an approximate solution whose accuracy is comparable to the solution obtained by the expensive standard FEM. Based on the LOD methodology, we investigate multiscale finite element methods for an elliptic distributed optimal control problem with rough coefficients. Numerical results obtained with our parallel implementation of the method illustrate our theoretical results.

Optimal Krylov Space Approximation of Rational Functions

We show how to efficiently generate the optimal (in a certain norm) Krylov space approximation to a function \(D(A)^{-1} N(A) b\) and we give bounds on the accuracy of such approximations.

Hierarchical matrix recovery using matrix-vector products

Can one recover a hierarchical matrix efficiently from only matrix-vector products? And, if so, then how many are needed? In this talk, we describe a stable algorithm for recovering an \(N \times N\) unknown hierarchical matrix using \(\mathcal{O}(k\log_2(N))\) input-output pairs of matrix-vector products, where \(k\) is the rank of the matrix's off-diagonal blocks. We do this by carefully constructing randomized vectors that use the hierarchical structure of the matrix. While there are existing algorithms that use a recursive ``peeling" procedure based on elimination, this approach can be numerically unstable in floating point arithmetic. Instead, we will use a recursive projection procedure to derive a more stable algorithm. We extend our work to the continuous setting by applying the algorithm to learn Green's functions associated with elliptic ordinary differential equations.

Optimal Size of the Block in Block GMRES on GPUs: Computational Model and Experiments

The block version of GMRES (BGMRES) is most advantageous over the single right hand side (RHS) counterpart when the cost of communication is high while the cost of floating point operations is not. This is the case on modern Graphics Processing Units(GPUs), while it is generally not the case on traditional Central Processing Units (CPUs). In this talk, experiments on both GPUs and CPUs are shown that compare the performance of BGMRES against GMRES as the number of RHS increases. The experiments indicatethat there are many cases in which BGMRES is slower than GMRES on CPUs, but faster on GPUs. Furthermore, when varying the number of RHS on the GPU, there is an optimal number of RHS where BGMRES is clearly most advantageous over GMRES. A computational modelis developed using hardware specific parameters, showing qualitatively where this optimal number of RHS is, and this model also helps explain the phenomena observed in the experiments.

Running on water, dry sand, and quicksand

The natural world is filled with an enormous diversity of substrates, yet, until recently, much of our understanding of how animals move was based upon studies of animals running along homogeneously smooth, solid trackways or treadmills. Running on yielding surfaces such as water (Newtonian), dry sand (non-Newtonian), and quicksand (non-Newtonian, shear-thickening) is fundamentally a different dynamic problem than when running on solid ground. Behaviors of all three materials change linearly or nonlinearly with how a force is applied. For example, dry sand makes transitions during a step among solid-like, fluid-like, and gas-like states. In contrast, quicksand solidifies during high-shear, high-velocity movements and is fluid-like, otherwise. Furthermore, whereas running on solid ground involves relatively simpler kinetics of exerted force translating into an equal and opposite reaction force, running across fluids and granular media is determined by a combination of added mass and dissipative effects. While designing a robot that can effectively traverse surfaces with such rapidly changing, complex properties is an engineer’s nightmare, animals do this daily and seemingly effortlessly. A key feature that distinguishes animals from engineered machines is their ability to use the same morphological appendage in multifunctional ways. In this talk, I will discuss our findings on bulk material behavior under different biologically-relevant intrusion conditions and connect this back to what the animals are doing above, and below, the surface to achieve the high levels of running performance we observe in nature.

Digital twins: Modeling, simulation and control of modern energy systems (Public Lecture)

The next level of digitization will create digital twins of every technological product or process. An important application for such digital twins is the modeling, simulation and operation of energy networks (e.g. power, gas, or district heating networks).

To model, simulate and control such networks in a sound and risk and error controlled way, a new modeling, simulation and optimization paradigm is needed. While automated modularized modeling is common in some domains like circuit design or multibody dynamics, it becomes increasingly challenging when systems or numerical solvers from different physical domains are coupled, due to largely different scales or modeling accuracy, and very different software technologies.

A recent system theoretic approach to address these challenges is the use of network based energy based modeling via constrained port-Hamiltonian (pH) systems, where the coupling is done in a physically meaningful way via energy variables. Furthermore, for each subsystem a whole model hierarchy can be employed ranging from very fine grain models to highly reduced surrogate models arising from model reduction or data based modeling. The model hierarchy allows adaptivity not only in the discretization but also in the model selection.

This talk will demonstrate the success of the new paradigm at the hand of several real world examples of energy networks.

Super Polynomials that Approximate the Inverse of a Matrix

We look at approximating the inverse of a matrix with a polynomial of the matrix. For a sparse matrix, this in some sense gives a sparse version of the inverse. It will be debated whether it is fair to call such polynomials “Superhero Polynomial”. And if so, are they more like Superman or Batman? A polynomial related to the GMRES polynomial is used, and it is implement using the GMRES polynomial roots. We show that this polynomial can often give a remarkably accurate approximation to the inverse. For ill-conditioned matrices, a high degree polynomial is needed and a composite of two polynomials may be most efficient and stable. Applications include large systems of linear equations with multiple right-hand sides and some deflation may be mentioned, hopefully with no deflation of the audience.

Error Estimates for Golub-Kahan Bidiagonalization with Tikhonov Regularization for Ill-posed Operator Equations

Linear ill-posed operator equations arise in various areas of science and engineering. The presence of errors in the operator and the data often makes the computation of an accurate approximate solution difficult. We compute an approximate solution of an ill-posed operator equation by first determining an approximation of the operators of generally fairly small dimension by carrying out a few steps of a continuous version of the Golub-Kahan bidiagonalization (GKB) process to the noisy operator. Then Tikhonov regularization is applied to the low-dimensional problem so obtained and the regularization parameter is determined by solving a low-dimensional nonlinear equation. The effect of replacing the original operator by the low-dimensional operator obtained by the GKB process on the accuracy of the solution is analyzed, as is the effect of errors in the operator and data. Computed examples that illustrate the theory are presented.

This talk presents joint work with A. Alqahtani, T. Mach, and R. Ramlau.

Characterization and convergence improvement of some Krylov subspace methodes for solving linear systems

NOSAS: Nonoverlapping Spectral Additive Schwarz Methods

We introduce novel non-overlapping additive Schwarz methods and focus on elliptic problems with highly heterogeneous coefficients. One of the NOSAS preconditioners is based on the approximation of the assembled Schur complement matrix by its diagonal part plus a low-rank perturbation, so, the Sherman–Morrison–Woodbury formula can be applied.   The dimension of this perturbation is the sum of the number of small eigenvalues of each local Schur complement. Another NOSAS preconditioner can be interpreted as glueing together local local modes subdomain eigenfunctions across the subdomain interfaces by the lowest energy. We also characterize these low eigenvalues with respect to the geometry of the coefficients (number of connected island/channels with large coefficients that touch the boundary of the subdomains. If time permits, we consider economical and three-level versions. The NOSAS methods are robust with respect to any diffusion coefficient and have good parallelization properties. Numerical results are presented.

Joint work with Yu Yi (U. of Hong Kong) and Maksymilian Dryja (U. of Warsaw)

Low Tensor-Train Rank Methods to Solve Sylvester Tensor Equations

In this talk, we focus on designing algorithms to solve Sylvester tensor equations in tensor-train (TT) format. The cornerstone of our solvers is the factored alternating direction implicit (fADI) method for Sylvester matrix equations. For different rank structures of the displacement tensor, we combine fADI with TTSVD or parallel-TTSVD to obtain the TT cores of the solution tensor directly. The solvers have optimal running complexity, and are practical in applications such as solving 3D Poisson equations on the cube. This talk is based on joint work with Alex Townsend and Maximilian Ruth.

Adaptive meshing for inverse problems

The estimation of distributed parameters from indirect observations through discretization is often complicated by the fact that the problemsare ill-posed, and discretization errors act like highly correlated noise. In this talk, some novel ideas of using adaptive meshing in computational inverse problems are discussed.

Golub-Kahan bidiagonalization with hybrid projection and recycling for streaming and very large problems

We discuss Golub-Kahan-based methods that exploit compression and recycling techniques to solve a broad class of inverse problems where memory requirements or high computational cost may otherwise be prohibitive. For problems with many unknown parameters that require many iterations, these methods compress and recycle the set of solution basis vectors to reduce the number of vectors that must be stored, while obtaining a solution accuracy comparable to that of standard methods. In other scenarios, such as streaming data problems or inverse problems with multiple datasets, hybrid projection methods with recycling can be used to efficiently integrate previously computed information for faster and better reconstruction. The algorithm allows accessing the matrix one block at a time, this can be a block of rows or a block of columns. The possibility of solving a linear system with only a subset of columns of the matrix available at a time also allows a range of interesting solution strategies for general matrices.

Building a bridge between numerical linear algebra and theoretical computer science

The numerical linear algebra (NLA) and theoretical computer science (TCS) research communities work on similar linear algebra problems from strikingly different viewpoints. During the pandemic, we organized an online fortnightly seminar to explore the connections between these two communities. In this talk, I discuss some observations in how the two communities formulate problems, design and analyse algorithms, and publicise their findings. The aim of the seminar was to foster future collaborations between NLA and TCS and generally bring about a greater appreciation for each other's work. The organizing committee for the seminar was Ilse Ipsen (NCSU), Mike Mahoney (Berkeley), Yuji Nakatsukasa (University of Oxford), Nikhil Srivastava (Berkeley), Alex Townsend (Cornell), and Joel Tropp (Caltech), along with founding members Daniel Kressner (EPFL) and Cleve Moler (MathWorks).

On the Development of Overlapping Schwarz Methods

This talk will try to trace how the overlapping Schwarz methods developed from a means to prove the existence of solution of elliptic problems to a major class of iterative algorithms for very large finite element problems. The core of this story is how the coarse global problems have been modified to handle irregular subdomains and coefficients. A discussion of recent results using adaptive algorithms to select the basis for efficient coarse spaces will also be included.

Register to attend at THIS FORM so that we know you are coming.

• Plane: There are more many daily nonstop flights to Philadelphia International Airport (PHL). From the airport, one can travel via SEPTA mass transit (25 minutes) or a taxi (15-25 minutes) to cover the eight miles from PHL to Center City.
• Train: Amtrak Acela and commuter trains arrive at 30th Street Station frequently, which is a short walk/ride to Center City. Philadelphia is only 1 hour and 20 minutes from New York City and 1 hour and 45 minutes from Washington, D.C., via Acela Express.
• Bus: For those who prefer bus travel, Greyhound, Megabus, and BoltBus offer convenient and economical travel to Philadelphia. These offer daily arrivals from all parts of the country.
• Car: Interstates 95 and 76, and the Pennsylvania and New Jersey Turnpikes, provide access from all directions. Parking in garages is plentiful Center City. Pay with cash or credit or purchase a reloadable Smart Card. See philapark.org for more information.

Temple University has its roots in the education of non-traditional working students who were not able to attend university courses but wanted to advance their educations. It began as a tutoring program delivered by a well-known Philadelphia minister, Russell Conwell, for young working people. More information can be found here.

This conference will take place on the Main Campus of Temple University, located in North Philadelphia on Broad Street. The conference hotel is also on Broad Street. This campus can be reached by car. Thus a group of attendees staying at the same hotel could also share an Uber/Taxi to get up to Temple. However, the north-south running Orange Subway Line has a stop at Temple's main campus on Broad Street and may be more convenient for many attendees. Some directions can be found at this web page.

The subway stop is on the corner of the hotel Broad and Locust st. Take the Orange Broad Street line Northbound LOCAL. Get off at Cecil B Moore (Temple University), go North one block, cross Montogomery Ave. and one more block, pass the Performing Arts Center (formerly the Baptist Temple), enter campus through the gates, first building on the left is Mitten Hall. Attendees will need to SIGN IN when entering Mitten Hall. The conference room is 250 Mitten Hall. You can walk to the end and take the elevator to the 2nd floor. Or go one flight up the stairs, turn right as if you were exiting the building and go up one more flight.

The conference will primarily take place at Temple University's Mitten Hall 250. The exception is the Public Interest Lecture, which will take place in SERC 116.

A map of Temple's campus has also been prepared.

Temple University has an extensive list of dining options near campus. The majority of the dining options are found on Liacouras Walk, the Student Center (13th St & Montgomery Ave), Morgan Hall (Broad St & Cecil B Moore Ave), and adjacent to View at Montgomery (11th St & Montgomery Ave). Below is a list of a few recommendations along with their locations.

If you want to have a casual lunch with a larger group

• Student Center
• Morgan Hall

Salads and healthy-ish options

• Playa Bowl (Liacouras Walk)
• Saladworks (Student Center)
• Honeygrow (below Morgan Hall, southeast corner of Broad St & Cecil B Moore Ave)
• Chipotle Mexican Grill (11th St & Montgomery Ave)

If you want a beer with your lunch

• Maxi's Pizza (Liacouras Walk)

Sandwiches and Quick Meals

• Potbelly Sandwiches (11th St & Montgomery Ave)
• Old Nelson (11th St Montgomery Ave)
• Subway (Liacouras Walk)

Non healthy options

• Chick-Fil-A (Student Center)
• Burger Fi (Student Center)
• Blaze Pizza (11th St & Montgomery Ave)
• Maxi’s Pizza (Liacouras Walk)

Coffee

• Richie's Cafe (11th & Berks St.; near SERC)
• Saxby’s (Liacouras Walk, Alter Hall)
• Cloud (food truck at 13th & Norris St)
• Starbucks (Student Center, Morgan Hall)

A map of Temple's campus may also be helpful.

The organizers have arranged for a block of rooms to be made available for booking at the DoubleTree Hotel in Center City Philadelphia. Booking can be arranged at this webpage. The advantages of the DoubleTree is that it is on Broad Street, under which runs the Orange Subway Line. Temple University is also on the Orange Line, north of Center City.

In addition, Philadelphia boasts many other hotels and other accomodations available through sites such as Airnb, etc.

As of this writing, there is a requirement at Temple University for everyone to use either (1) N95 or KN95 masks or (2) double mask. Temple University has a webpage describing the covid situation and protocols put into place.