Evan S. Gawlik

Evan S. Gawlik, Melvin Leok

We derive a family of high-order, structure-preserving approximations of the Riemannian exponential map on several matrix manifolds, including the group of unitary matrices, the Grassmannian manifold, and the Stiefel manifold. Our derivation is inspired by the observation that if $Ω$ is a skew-Hermitian matrix and $t$ is a sufficiently small scalar, then there exists a polynomial of degree $n$ in $tΩ$ (namely, a Bessel polynomial) whose polar decomposition delivers an approximation of $e^{tΩ}$ with error $O(t^{2n+1})$. We prove this fact and then leverage it to derive high-order approximations of the Riemannian exponential map on the Grassmannian and Stiefel manifolds. Along the way, we derive related results concerning the supercloseness of the geometric and arithmetic means of unitary matrices.

Evan S. Gawlik, Yuji Nakatsukasa, Brian D. Sutton

We introduce a backward stable algorithm for computing the CS decomposition of a partitioned $2n \times n$ matrix with orthonormal columns, or a rank-deficient partial isometry. The algorithm computes two $n \times n$ polar decompositions (which can be carried out in parallel) followed by an eigendecomposition of a judiciously crafted $n \times n$ Hermitian matrix. We prove that the algorithm is backward stable whenever the aforementioned decompositions are computed in a backward stable way. Since the polar decomposition and the symmetric eigendecomposition are highly amenable to parallelization, the algorithm inherits this feature. We illustrate this fact by invoking recently developed algorithms for the polar decomposition and symmetric eigendecomposition that leverage Zolotarev's best rational approximations of the sign function. Numerical examples demonstrate that the resulting algorithm for computing the CS decomposition enjoys excellent numerical stability.

Evan S. Gawlik, Adrian J. Lew

We derive upper bounds on the difference between the orthogonal projections of a smooth function $u$ onto two finite element spaces that are nearby, in the sense that the support of every shape function belonging to one but not both of the spaces is contained in a common region whose measure tends to zero under mesh refinement. The bounds apply, in particular, to the setting in which the two finite element spaces consist of continuous functions that are elementwise polynomials over shape-regular, quasi-uniform meshes that coincide except on a region of measure $O(h^γ)$, where $γ$ is a nonnegative scalar and $h$ is the mesh spacing. The projector may be, for example, the orthogonal projector with respect to the $L^2$- or $H^1$-inner product. In these and other circumstances, the bounds are superconvergent under a few mild regularity assumptions. That is, under mesh refinement, the two projections differ in norm by an amount that decays to zero at a faster rate than the amounts by which each projection differs from $u$. We present numerical examples to illustrate these superconvergent estimates and verify the necessity of the regularity assumptions on $u$.

Evan S. Gawlik

We construct a family of iterations for computing the principal square root of a square matrix $A$ using Zolotarev's rational minimax approximants of the square root function. We show that these rational functions obey a recursion, allowing one to iteratively generate optimal rational approximants of $\sqrt{z}$ of high degree using compositions and products of low-degree rational functions. The corresponding iterations for the matrix square root converge to $A^{1/2}$ for any input matrix $A$ having no nonpositive real eigenvalues. In special limiting cases, these iterations reduce to known iterations for the matrix square root: the lowest-order version is an optimally scaled Newton iteration, and for certain parameter choices, the principal family of Padé iterations is recovered. Theoretical results and numerical experiments indicate that the iterations perform especially well on matrices having eigenvalues with widely varying magnitudes.

Evan S. Gawlik

In [E. S. Gawlik, Zolotarev iterations for the matrix square root, arXiv preprint 1804.11000, (2018)], a family of iterations for computing the matrix square root was constructed by exploiting a recursion obeyed by Zolotarev's rational minimax approximants of the function $z^{1/2}$. The present paper generalizes this construction by deriving rational minimax iterations for the matrix $p^{th}$ root, where $p \ge 2$ is an integer. The analysis of these iterations is considerably different from the case $p=2$, owing to the fact that when $p>2$, rational minimax approximants of the function $z^{1/p}$ do not obey a recursion. Nevertheless, we show that several of the salient features of the Zolotarev iterations for the matrix square root, including equioscillatory error, order of convergence, and stability, carry over to case $p>2$. A key role in the analysis is played by the asymptotic behavior of rational minimax approximants on short intervals. Numerical examples are presented to illustrate the predictions of the theory.

Evan S. Gawlik

A widely used approximation of the Gaussian curvature on a triangulated surface is the angle defect, which measures the deviation between $2π$ and the sum of the angles between neighboring edges emanating from a common vertex. We show that the linearization of the angle defect about an arbitrary piecewise constant Regge metric is related to the classical Hellan-Herrmann-Johnson finite element discretization of the div-div operator. Integrating this relation leads to an integral formula for the angle defect which is well-suited for analysis and generalizes naturally to higher order. We prove error estimates for these high-order approximations of the Gaussian curvature in $H^k$-Sobolev norms of integer order $k \ge -1$.

Evan S. Gawlik, Melvin Leok

We study schemes for interpolating functions that take values in the special orthogonal group $SO(n)$. Our focus is on interpolation schemes obtained by embedding $SO(n)$ in a linear space, interpolating in the linear space, and mapping the result onto $SO(n)$ via the closest point projection. The resulting interpolants inherit both the order of accuracy and the regularity of the underlying interpolants on the linear space. The values and derivatives of the interpolants admit efficient evaluation via either explicit formulas or iterative algorithms, which we detail for two choices of embeddings: the embedding of $SO(n)$ in the space of $n \times n$ matrices and, when $n=3$, the identification of $SO(3)$ with the set of unit quaternions. Along the way, we point out a connection between these interpolation schemes and geodesic finite elements. We illustrate the utility of these interpolation schemes by numerically computing minimum acceleration curves on $SO(n)$, a task which is handled naturally with $SO(n)$-valued finite elements having $C^1$-continuity.

Evan S. Gawlik, Melvin Leok

We derive iterative methods for computing the Fréchet derivative of the map which sends a full-rank matrix $A$ to the factor $U$ in its polar decomposition $A=UH$, where $U$ has orthonormal columns and $H$ is Hermitian positive definite. The methods apply to square matrices as well as rectangular matrices having more rows than columns. Our derivation relies on a novel identity that relates the Fréchet derivative of the polar decomposition to the matrix sign function $\mathrm{sign}(X) = X (X^2)^{-1/2}$ applied to a certain block matrix $X$.

Maurizio M. Chiaramonte, Evan S. Gawlik, Hardik Kabaria et al.

Universal meshes have recently appeared in the literature as a compu- tationally efficient and robust paradigm for the generation of conforming simpli- cial meshes for domains with evolving boundaries. The main idea behind a univer- sal mesh is to immerse the moving boundary in a background mesh (the universal mesh), and to produce a mesh that conforms to the moving boundary at any given time by adjusting a few of elements of the background mesh. In this manuscript we present the application of universal meshes to the simulation of brittle fracturing. To this extent, we provide a high level description of a crack propagation algorithm and showcase its capabilities. Alongside universal meshes for the simulation of brit- tle fracture, we provide other examples for which universal meshes prove to be a powerful tool, namely fluid flow past moving obstacles. Lastly, we conclude the manuscript with some remarks on the current state of universal meshes and future directions.