GAMM ANLA Book of Abstracts
Henrik Eisenmann: Local definite Multiparameter eigenvalue problems
Multiparameter eigenvalue problems (MEP) combine aspects from eigenvalue problems and linear systems of equations and generalize both. Locally definite MEPs are a generalization of symmetric eigenvalue problems to problems with multiple eigenvalue parameters. We present an approach to compute eigenvalues and eigenvectors of locally definite multiparameter eigenvalue problems by its signed multiindex. The definiteness assumptions are naturally satisfied by multiparameter Sturm-Lioville problems that arise when variables of multidimensional boundary eigenvalue problems can be separated.Tom Werner: On a matrix-Newton-type framework for solving NEPv
In this talk, we want to introduce a novel framework for solving nonlinear matrix equations by inexact matrix-valued Newton-Krylov methods. In this framework, we consider global Krylov subspace methods for solving the correction equation and discuss the benefits of this approach. Finally, we show how to use the proposed algorithm to solve nonlinear eigenvector dependent eigenvalue problems (NEPv) and compare with the well-established self-consistent field iteration on different applications.Tom-Christian Riemer: Low-Rank Methods for IGA with Multiple Patches
Isogeometric analysis (IGA) has become one of the most popular methods for the discretization of partial differential equations (PDEs) on complex domains, which arise in various applications of science and engineering. In many cases, it is necessary to describe the computational domain of a complex geometry using multiple patches, with each patch represented by a tensor-product B-spline parameterization, and adjacent patches connected by interfaces. A key challenge in this setup is solving the discretized equations resulting from a PDE, as the discretization typically produces large mass and stiffness tensors that are costly to assemble and may be infeasible to handle with conventional methods due to their size. Additionally, continuity of the approximation across patch interfaces must be ensured. In this talk, we illustrate how low-rank methods based on the tensor train format, combined with the Block-AMEn (alternating minimal energy) iterative solver, can be extended to multi-patch geometries by integrating an isogeometric tearing and interconnecting (IETI) approach. In this approach, the discretization of the PDE is expressed in a low-rank structure using Kronecker products across spatial dimensions. Continuity on the respective interfaces is enforced by extending the linear system with constraint tensors and Lagrange multipliers. Numerical tests implemented in MATLAB using the GeoPDEs toolbox confirm this theory.Eda Oktay: Reorthogonalized Pythagorean variants of block classical Gram Schmidt
The computation of Krylov subspace methods (KSMs) and orthogonalization procedures requires sparse matrix-vector products (SpMV) and inner products in each iteration. These operations are also called synchronization points, due to the data movement between levels of cache on a node or the message-passing operations necessary in distributed computing, and they cause KSMs to be communication bound. To reduce communication costs, communication-avoiding Krylov methods can be employed, which make use of the sparse matrix powers kernel and a batched or “block"" orthogonalization of basis vectors. Block Classical Gram-Schmidt (BCGS) is a common choice for the block orthogonalization process, as it requires only two synchronization points in addition to an SpMV per iteration. However, BCGS is known to have a high loss of orthogonality, which can contribute to instability and inaccurate solutions in the overall KSM. To balance performance and stability, a so-called “Pythagorean"" variant of BCGS was introduced in [E. Carson, K. Lund, & M. Rozložník. The stability of block variants of classical GramSchmidt. SIAM J. Matrix Anal. Appl. 42(3), pp. 1365–1380, 2021].We introduce two reorthogonalized Pythagorean BCGS variants, each with the same number of synchronization points as BCGS, to improve stability. We show how reorthogonalization can restore a bound on the loss of orthogonality of order of the working precision both in theory and in practice.
Yuxin Ma: On the backward stability of s-step GMRES
Communication, i.e., data movement, is a critical bottleneck for the performance of classical Krylov subspace method solvers on modern computer architectures. Variants of these methods which avoid communication have been introduced, which, while equivalent in exact arithmetic, can be unstable in finite precision. In this work, we present a framework to analyze s-step GMRES with popular block orthogonalization methods and further discuss the instability of s-step GMRES.Xinye Chen: Computing k-means in mixed precision
The k-means algorithm is one of the most popular and critical techniques in data mining and machine learning, and it has achieved significant success in numerous science and engineering domains. Computing k-means to a global optimum is NP-hard in Euclidean space, yet there are a variety of efficient heuristic algorithms, such as Lloyd's algorithm, that converge to a local optimum with superpolynomial complexity in the worst case.In this talk we will discuss the numerical stability of Lloyd's k-means algorithm, and, in particular, we confirm the stability of the widely used distance computation formula. We propose a mixed-precision framework for k-means computation and investigate the effects of low-precision distance computation within the framework. Through extensive simulations on various data clustering and image segmentation tasks, we verify the applicability and robustness of the mixed precision k-means method. We find that, in k-means computation, normalized data is more tolerant to the reduction of precision in the distance computation, while for nonnormalized data more care is needed in the use of reduced precision, mainly to avoid overflow. Our study demonstrates the potential for the use of mixed precision to accelerate the k-means computation and offers some insights into other distance-based machine learning methods.
Robert Altmann: A New Energy-based Modeling Framework for DAEs
TBAJohannes Maly: Algorithmic approaches to recovering sparse and low-rank matrices
In this talk, we consider the problem of recovering an unknown sparse and low-rank matrix from observations gathered in a linear measurement process. We discuss the challenges that come with leveraging several structures simultaneously and present two algorithmic strategies to efficiently approach the problem. Both strategies come with local convergence guarantees.Mehdi Babamehdi: Using nonlinear domain decomposition as smoother in nonlinear multigrid
Nonlinear Partial Differential Equations (PDEs) frequently arise in different fields of science, e.g. material science, fluid dynamics etc. Finite element discretization of such problems usually leads to large nonlinear systems, either due to the need of higher accuracy, having to deal with large physical sizes or both. Solution of such big discretized nonlinear problems needs fast, highly scalable, and parallelize solvers.Nonlinear multigrid is a well-known method for efficiently solving nonlinear boundary value problems. The full approximation scheme (FAS) solves nonlinear problems on fine and coarse grids. To smooth the nonlinear problem a nonlinear smoother must be utilized and since implementation is supposed to be in matrix-free form for an efficient finite element code, this form of implementation of the smoother should be plausible. For this purpose, Nonlinear Restricted Additive Schwarz Method (NRASM) seems to be a very interesting choice. The method converges with the same rate as linear iterations applied to the linearised equation. In addition, it is inherently parallel and proper to be implemented in matrix-free form.
We combine FAS and a nonlinear restricted additive Schwarz method to obtain hybrid NRASM/FAS. The full approximation scheme is used to solve the nonlinear problem and the NRASM is used to smooth the nonlinear boundary value problem in local subdomains on each level of the multigrid method. Within the NRASM, we use the Jacobian-Free Newton Krylov method as a solver on each subdomain. We consider different nonlinear equations in three dimensional space as test problems. Several parameters of the methods were investigated in this study to have a better understanding of influence of the parameters on the efficiency of the method, the behavior of the methods and its convergence rate.
Jacob B. Schroder: Generalized Optimal AMG Convergence Theory for Nonsymmetric and Indefinite Problems
Algebraic multigrid (AMG) is a popular and effective solver for sparse linear systems arising from discretized PDEs, and especially for thesymmetric positive definite (SPD) case. In this SPD setting, the convergence and theory of AMG is well-developed and understood in terms of the A-norm.
However for nonsymmetric systems, such an energy norm is non-existent. For this reason, convergence of AMG for nonsymmetric systems remains an open area of research. In the SPD setting, the classical form of optimal AMG interpolation provides a useful insight in determining the best possible
two-grid convergence rate for a given symmetrized relaxation method. In this work, we discuss a generalization of the optimal AMG convergence theory for
nonsymmetric problems. The theory relies on a certain matrix-induced orthogonality of the left and right eigenvectors of a generalized eigenvalue
problem based on the system matrix and relaxation operator. The result is that we can measure the spectral radius of the two-grid error operator, in a way
that is equivalent to the derivation in the SPD setting, where instead norms are used. We also provide some supporting numerical evidence of the
theoretical estimate, as well as possible practical applications of the theory.
Thomas Mach: Solution of linear ill-posed operator equations by modified truncated singular value expansion
Often the solution of linear ill-posed operator equations, the operator equation is discretized and regularization methods are developed for the discretized problem so obtained, without discussing the ramification of these methods for the infinite-dimensional problem. In particular, these regularization methods may only be applicable to certain linear ill-posed operator equations. We discuss how regularization by a modified truncated singular value decomposition introduced for finite-dimensional problems can be extended to operator equations. In finite dimensions, this regularization method yields approximate solutions of higher quality than standard truncated singular value decomposition. Our analysis in a Hilbert space setting is of practical interest, because the solution method presented avoids introduction of discretization errors during the solution process. We demonstrate the results using Chebfun.This is joint work with Laura Dykes, Mykhailo Kuian, Silvia Noschese, and Lothar Reichel.
Wim Vanroose: Residual quadratic programming active set subspace method
We present the Residual Quadratic Programming Active-Set Subspace (ResQPASS) method that solves large-scale linear least-squares problems with bound constraints on the variables. The problem is solved by creating a series of small problems of increasing size by projecting on the basis of residuals. Each projected problem is solved by the active-set method for convex quadratic programming, warm-started with a working set and solution of the previous problem. The method coincides with conjugate gradients (CG) or, equivalently, LSQR applied to the normal equations when none of the constraints is active. When only a few constraints are active the method converges, after a few initial iterations, like the CG and LSQR. An analysis links the convergence to an asymptotic Krylov subspace. We also present an efficient implementation where QR factorizations of the projected are updated over the inner iterations and Cholesky over the outer iterations.Sascha Portaro: Row-aware Randomized RSVD
We introduce a novel procedure for computing an SVD-type approximation of a tall matrix A. Specifically, we propose a randomization-based algorithm that improves the standard Randomized Singular Value Decomposition (RSVD). Most significantly, our approach, the Row-aware RSVD (R-RSVD), explicitly constructs information from the row space of A. This leads to better approximations to Range(A) while maintaining the same computational cost. The efficacy of the R-RSVD is supported by both robust theoretical results and extensive numerical experiments. Furthermore, we present an alternative algorithm inspired by the R-RSVD, capable of achieving comparable accuracy despite utilizing only a subsample of the rows of A, resulting in a significantly reduced computational cost. This method, that we name the Subsample Row-aware RSVD (Rsub- SVD), is supported by a weaker error bound compared to the ones we derived for the R-RSVD, but still meaningful as it ensures that the error remains under control. Additionally, numerous experiments demonstrate that the Rsub-RSVD trend is akin to the one attained by the R-RSVD when the subsampling parameter is on the order of n, for a mxn real matrix A, with m ≫ n. Finally, we consider the application of our schemes in two very diverse settings which share the need for the computation of singular vectors as an intermediate step: the computation of CUR decompositions by the discrete empirical interpolation method (DEIM) and the construction of reduced-order models in the Loewner framework, a data-driven technique for model reduction of dynamical systems.Marcel Schweitzer: Polynomial preconditioning for the action of the matrix square root and inverse square root
While preconditioning is a long-standing concept to accelerate iterative methods for linear systems, generalizations to matrix functions are still in their infancy. We go a further step in this direction, introducing polynomial preconditioning for Krylov subspace methods that approximate the action of the matrix square root and inverse square root on a vector. Preconditioning reduces the subspace size and therefore avoids the storage problem together with—for non-Hermitian matrices—the increased computational cost per iteration that arises in the unpreconditioned case. Polynomial preconditioning is an attractive alternative to current restarting or sketching approaches since it is simpler and computationally more efficient. We demonstrate this for several numerical examples.Anne Wald: Identification of forces in biophysical experiments
Many processes in cells are driven by the interaction of multiple proteins, for example cell contraction, division or migration. For a better understanding of the underlying mechanics, one is often interested in identifying the forces or stresses that are responsible for such processes. In this talk, we address two relevant biophysical experiments that involve the numerical solution of the respective inverse problems. First, we consider droplets containing an actomyosin network, which consists mainly of two important types of proteins: myosin and actin filaments. Their interaction leads to an evolution of the network which creates a flow inside the droplet by exerting forces on the fluid. These can be determined by solving a parameter identification problem for the Stokes equation. In addition, on a somewhat larger scale, we will discuss traction force microscopy. The goal is to identify forces generated by the cytoskeleton during cell contraction. To this end, the cell is adhered to a soft elastic substrate, whose deformation can be observed under a microscope. From measurements of the deformation we calculate the respective forces, again by means of parameter identification, but for suitable elasticity equations. In both cases we present some numerical results for measurement data.Lukas Klingbiel: A Riemannian Douglas-Rachford Algorithm for Low-Rank and Row-Sparse Matrix Recovery
In this work, we consider a matrix recovery problem under low-rank and row-sparse constraints using a Riemannian Douglas-Rachford algorithm. We introduce a retraction-based Riemannian Douglas-Rachford based on the Riemannian Douglas-Rachford on symmetric Hadamard manifolds. We show local convergence of this method for nonexpansive reflections and manifolds with locally bounded sectional curvature. We give an explicit form of the algorithm on the fixed-rank manifold to solve the matrix recovery problem. In particular, numerical experiments suggest that we require the minimal number of measurements for the recovery of a low-rank and row-sparse matrix.Jakob Stoye: Injectivity Radius and Curvature Bounds on the Stiefel Manifold
This presentation explores key geometric properties of the Stiefel manifold, focusing on its injectivity radius and sectional curvature under both the Euclidean and canonical metric. It is established that the injectivity radius of the Stiefel manifold is precisely π under the Euclidean metric, a result derived from the observation that geodesics on this manifold are space curves of constant Frenet curvatures. Additionally, the sectional curvature of the Stiefel manifold is examined, confirming that the global curvature bounds of 5/4 under the canonical metric and 1 under the Euclidean metric hold true across all admissible dimensions. Notably, the maximum curvature is attained by rank-two matrices in the canonical case and rank-one matrices in the Euclidean case. These findings provide a comprehensive understanding of curvature behavior, offering valuable insights for Riemannian computing methods.Ion Victor Gosea: Extending balanced truncation for model reduction of bilinear systems with quadratic outputs
Bilinear systems with quadratic outputs (BQOs) are a special class of bilinear systems for which the output expression contains a quadratic form of the state variable. Such dynamical systems are characterized by differential and standard equations with bilinear and quadratic nonlinearities and are of use, e.g., when one is interested in observing energies or other quadratic quantities of a bilinear system as an extension of the case of linear systems [1]. However, as far as the authors are aware, model order reduction (MOR) techniques for such systems have not been studied so far. In this contribution, we propose a balanced truncation (BT) approach for the MOR of BQO systems. We use the Volterra series representation of a BQO system and the state-space representation of its dual system [4] to derive the time-domain generalized kernels of these systems and define infinite Gramians (P and Q) for such systems. We prove various relations between the newly-defined Gramians and the energy functionals of BQO systems and show that the Gramians satisfy generalized Lyapunov matrix equations (which are still linear in both P and Q variables). We put together the proposed BT algorithm for BQO systems. We also show that the reduced-order model is balanced in the generalized sense. It is to be noted that solving the Lyapunov equations can be performed by employing iterative numerical methods such as in [3]. Instead of solving these generalized Lyapunov equations, we also propose using truncated Gramians similar to the approach in [2]. We then present an efficient BT algorithm based on such Gramians and also investigate the possibility of deriving error bounds to quantify the approximation quality. Finally, various numerical experiments are reported, such as a semi-discretized heat transfer model with Robin boundary conditions, for which the observed output is given by the energy of the temperature values instead of the average thereof.[1] P. Benner, P. Goyal, and I. Pontes Duff. Gramians, energy functionals, and balanced truncation for linear dynamical systems with quadratic outputs. IEEE Transactions on Automatic Control, 67(2):886–893, 2021.
[2] P. Benner, P. Goyal, and M. Redmann. Truncated Gramians for bilinear systems and their advantages in model order reduction. In P. Benner and et al., editors, Model Reduction of Parametrized Systems, volume 17 of MS&A - Modeling, Simulation and Applications, pages 285–300. Springer, Cham, 2017.
[3] P. Benner, M. Köhler, and J. Saak. Matrix equations, sparse solvers: M-M.E.S.S.-2.0.1 – philosophy, features and application for (parametric) model order reduction. In P. Benner and et al., editors, Model Reduction of Complex Dynamical Systems, volume 171, pages 369–392. Birkhäuser, Cham, 2021.
[4] K. Fujimoto, J. M. Scherpen, and W. Steven Gray. Hamiltonian realizations of nonlinear adjoint operators. IFAC Proceedings Volumes, 33(2):39–44, 2000. IFAC Workshop on Lagrangian and Hamiltonian Methods for Nonlinear Control, Princeton, NJ, USA, 16-18 March 2000.
Jonas Püschel: An Energy-Adaptive Riemannian Conjugate Gradient Method for Kohn--Sham Type Eigenvector Problems
Kohn-Sham type eigenvector problems often arise in computational physics and chemistry and can be formulated as minimization of an energy functional E on the Stiefel manifold St(n,p). We propose to use the Hamiltonian, a space-dependent symmetric linear map which connects to E via the gradient, to construct a Riemannian metric on St(n,p). This allows us to formulate the Riemannian energy-adaptive conjugate gradient method. Numerical experiments illustrate that our method outperforms established Riemannian methods and is competitive with state-of-the-art self consistent field (SCF) based methods.Jan Martin Nicolaus: Nonintrusive model order reduction for stochastic differential equations
This work introduces a non-intrusive model order reduction method provided by extending the operator inference approach to linear and controlled stochastic differential equations with additive noise. The drift and diffusion coefficients of the reduced order model (ROM) are inferred from state observations by solving appropriate least-squares problems. The closeness of the ROM obtained by the presented approach to the intrusive ROM obtained by the proper orthogonal decomposition (POD) method is investigated. Two possible generalisations of the snapshot-based dominant subspace construction to the stochastic case are presented. Numerical experiments are provided to compare the developed approach to POD.This is joint work with Melina Freitag and Martin Redmann.