Search

The injectivity radius of a manifold is an important quantity, both from a theoretical point of view and in terms of numerical applications. It is the largest possible radius within which all geodesics are unique and length-minimizing. In consequence, it is the largest possible radius within which calculations in Riemannian normal coordinates are well-defined. A matrix manifold that arises frequently in a wide range of practical applications is the compact Stiefel manifold of orthogonal $p$-frames in $\mathbb{R}^n$. We observe that geodesics on this manifold are space curves of constant Frenet curvatures. Using this fact, we prove that the injectivity radius on the Stiefel manifold under the Euclidean metric is $\pi$., Comment: 10 pages

Riemannian optimization is concerned with problems, where the independent variable lies on a smooth manifold. There is a number of problems from numerical linear algebra that fall into this category, where the manifold is usually specified by special matrix structures, such as orthogonality or definiteness. Following this line of research, we investigate tools for Riemannian optimization on the symplectic Stiefel manifold. We complement the existing set of numerical optimization algorithms with a Riemannian trust region method tailored to the symplectic Stiefel manifold. To this end, we derive a matrix formula for the Riemannian Hessian under a right-invariant metric. Moreover, we propose a novel retraction for approximating the Riemannian geodesics. Finally, we conduct a comparative study in which we juxtapose the performance of the Riemannian variants of the steepest descent, conjugate gradients, and trust region methods on selected matrix optimization problems that feature symplectic constraints., Comment: 24 pages, 2 figures

Since the popularization of the Stiefel manifold for numerical applications in 1998 in a seminal paper from Edelman et al., it has been exhibited to be a key to solve many problems from optimization, statistics and machine learning. In 2021, H\"uper et al. proposed a one-parameter family of Riemannian metrics on the Stiefel manifold, subsuming the well-known Euclidean and canonical metrics. Since then, several methods have been proposed to obtain a candidate for the Riemannian logarithm given any metric from the family. Most of these methods are based on the shooting method or rely on optimization approaches. For the canonical metric, Zimmermann proposed in 2017 a particularly efficient method based on a pure matrix-algebraic approach. In this paper, we derive a generalization of this algorithm that works for the one-parameter family of Riemannian metrics. The algorithm is proposed in two versions, termed backward and forward, for which we prove that it conserves the local linear convergence previously exhibited in Zimmermann's algorithm for the canonical metric., Comment: 20 pages, 6 figures, 3 tables

Arguably, geodesics are the most important geometric objects on a differentiable manifold. They describe candidates for shortest paths and are guaranteed to be unique shortest paths when the starting velocity stays within the so-called injectivity radius of the manifold. In this work, we investigate the injectivity radius of the Stiefel manifold under the canonical metric. The Stiefel manifold $St(n,p)$ is the set of rectangular matrices of dimension $n$-by-$p$ with orthogonal columns, sometimes also called the space of orthogonal $p$-frames in $\mathbb{R}^n$. Using a standard curvature argument, Rentmeesters has shown in 2013 that the injectivity radius of the Stiefel manifold is bounded by $\sqrt{\frac{4}{5}}\pi$. It is an open question, whether this bound is sharp. With the definition of the injectivity radius via cut points of geodesics, we gain access to the information of the injectivity radius by investigating geodesics. More precisely, we consider the behavior of special variations of geodesics, called Jacobi fields. By doing so, we are able to present an explicit example of a cut point. In addition, since the theoretical analysis of geodesics for cut points and especially conjugate points as a type of cut points is difficult, we investigate the question of the sharpness of the bound by means of numerical experiments.

Methods and algorithms that work with data on nonlinear manifolds are collectively summarized under the term `Riemannian computing'. In practice, curvature can be a key limiting factor for the performance of Riemannian computing methods. Yet, curvature can also be a powerful tool in the theoretical analysis of Riemannian algorithms. In this work, we investigate the sectional curvature of the Stiefel and Grassmann manifold. On the Grassmannian, tight curvature bounds are known since the late 1960ies. On the Stiefel manifold under the canonical metric, it was believed that the sectional curvature does not exceed 5/4. Under the Euclidean metric, the maximum was conjectured to be at 1. For both manifolds, the sectional curvature is given by the Frobenius norm of certain structured commutator brackets of skew-symmetric matrices. We provide refined inequalities for such terms and pay special attention to the maximizers of the curvature bounds. In this way, we prove for the Stiefel manifold that the global bounds of 5/4 (canonical metric) and 1 (Euclidean metric) hold indeed. With this addition, a complete account of the curvature bounds in all admissible dimensions is obtained. We observe that `high curvature means low-rank', more precisely, for the Stiefel and Grassmann manifolds under the canonical metric, the global curvature maximum is attained at tangent plane sections that are spanned by rank-two matrices, while the extreme curvature cases of the Euclidean Stiefel manifold occur for rank-one matrices. Numerical examples are included for illustration purposes., Comment: 29 pages, 5 figures

In this paper, we propose two methods for multivariate Hermite interpolation of manifold-valued functions. On the one hand, we approach the problem via computing suitable weighted Riemannian barycenters. To satisfy the conditions for Hermite interpolation, the sampled derivative information is converted into a condition on the derivatives of the associated weight functions. It turns out that this requires the solution of linear systems of equations, but no vector transport is necessary. This approach treats all given sample data points equally and is intrinsic in the sense that it does not depend on local coordinates or embeddings. As an alternative, we consider Hermite interpolation in a tangent space. This is a straightforward approach, where one designated point, for example one of the sample points or (one of) their center(s) of mass, is chosen to act as the base point at which the tangent space is attached. The remaining sampled locations and sampled derivatives are mapped to said tangent space. This requires a vector transport between different tangent spaces. The actual interpolation is then conducted via classical vector space operations. The interpolant depends on the selected base point. The validity and performance of both approaches is illustrated by means of numerical examples.

Gradient-enhanced Kriging (GE-Kriging) is a well-established surrogate modelling technique for approximating expensive computational models. However, it tends to get impractical for high-dimensional problems due to the size of the inherent correlation matrix and the associated high-dimensional hyper-parameter tuning problem. To address these issues, a new method, called sliced GE-Kriging (SGE-Kriging), is developed in this paper for reducing both the size of the correlation matrix and the number of hyper-parameters. We first split the training sample set into multiple slices, and invoke Bayes' theorem to approximate the full likelihood function via a sliced likelihood function, in which multiple small correlation matrices are utilized to describe the correlation of the sample set rather than one large one. Then, we replace the original high-dimensional hyper-parameter tuning problem with a low-dimensional counterpart by learning the relationship between the hyper-parameters and the derivative-based global sensitivity indices. The performance of SGE-Kriging is finally validated by means of numerical experiments with several benchmarks and a high-dimensional aerodynamic modeling problem. The results show that the SGE-Kriging model features an accuracy and robustness that is comparable to the standard one but comes at much less training costs. The benefits are most evident for high-dimensional problems with tens of variables.

The real symplectic Stiefel manifold is the manifold of symplectic bases of symplectic subspaces of a fixed dimension. It features in a large variety of applications in physics and engineering. In this work, we study this manifold with the goal of providing theory and matrix-based numerical tools fit for basic data processing. Geodesics are fundamental for data processing. However, these are so far unknown. Pursuing a Lie group approach, we close this gap and derive efficiently computable formulas for the geodesics both with respect to a natural pseudo-Riemannian metric and a novel Riemannian metric. In addition, we provide efficiently computable and invertible retractions. Moreover, we introduce the real symplectic Grassmann manifold, i.e., the manifold of symplectic subspaces. Again, we derive efficient formulas for pseudo-Riemannian and Riemannian geodesics and invertible retractions. The findings are illustrated by numerical experiments, where we consider optimization via gradient descent on both manifolds and compare the proposed methods with the state of the art. In particular, we treat the 'nearest symplectic matrix' problem and the problem of optimal data representation via a low-rank symplectic subspace. The latter task is associated with the problem of finding a 'proper symplectic decomposition', which is important in structure-preserving model order reduction of Hamiltonian systems., Comment: 32 pages, 8 figures

Solving the so-called geodesic endpoint problem, i.e., finding a geodesic that connects two given points on a manifold, is at the basis of virtually all data processing operations, including averaging, clustering, interpolation and optimization. On the Stiefel manifold of orthonormal frames, this problem is computationally involved. A remedy is to use quasi-geodesics as a replacement for the Riemannian geodesics. Quasi-geodesics feature constant speed and covariant acceleration with constant (but possibly non-zero) norm. For a well-known type of quasi-geodesics, we derive a new representation that is suited for large-scale computations. Moreover, we introduce a new kind of quasi-geodesics that turns out to be much closer to the Riemannian geodesics., Comment: 8 pages, 1 figure, conference GSI2021

We address the problem of computing Riemannian normal coordinates on the real, compact Stiefel manifold of orthogonal frames. The Riemannian normal coordinates are based on the so-called Riemannian exponential and the Riemannian logarithm maps and enable to transfer almost any computational procedure to the realm of the Stiefel manifold. To compute the Riemannian logarithm is to solve the (local) geodesic endpoint problem. Instead of restricting the consideration to geodesics with respect to a single selected metric, we consider a family of Riemannian metrics introduced by H\"uper, Markina and Silva-Leite that includes the Euclidean and the canonical metric as prominent examples. As main contributions, we provide (1) a unified, structured, reduced formula for the Stiefel geodesics for the complete family of metrics, (2) a unified method to tackle the geodesic endpoint problem, (3) an improvement of the existing Riemannian log map under the canonical metric. The findings are illustrated by means of numerical examples, where the novel algorithms prove to be the most efficient methods known to this date., Comment: 32 pages, 10 figures

The Grassmann manifold of linear subspaces is important for the mathematical modelling of a multitude of applications, ranging from problems in machine learning, computer vision and image processing to low-rank matrix optimization problems, dynamic low-rank decompositions and model reduction. With this mostly expository work, we aim to provide a collection of the essential facts and formulae on the geometry of the Grassmann manifold in a fashion that is fit for tackling the aforementioned problems with matrix-based algorithms. Moreover, we expose the Grassmann geometry both from the approach of representing subspaces with orthogonal projectors and when viewed as a quotient space of the orthogonal group, where subspaces are identified as equivalence classes of (orthogonal) bases. This bridges the associated research tracks and allows for an easy transition between these two approaches. Original contributions include a modified algorithm for computing the Riemannian logarithm map on the Grassmannian that is advantageous numerically but also allows for a more elementary, yet more complete description of the cut locus and the conjugate points. We also derive a formula for parallel transport along geodesics in the orthogonal projector perspective, formulae for the derivative of the exponential map, as well as a formula for Jacobi fields vanishing at one point., Comment: 38 pages, 4 figures

The main contribution of this paper is twofold: On the one hand, a general framework for performing Hermite interpolation on Riemannian manifolds is presented. The method is applicable, if algorithms for the associated Riemannian exponential and logarithm mappings are available. This includes many of the matrix manifolds that arise in practical Riemannian computing application such as data analysis and signal processing, computer vision and image processing, structured matrix optimization problems and model reduction. On the other hand, we expose a natural relation between data processing errors and the sectional curvature of the manifold in question. This provides general error bounds for manifold data processing methods that rely on Riemannian normal coordinates. Numerical experiments are conducted for the compact Stiefel manifold of rectangular column-orthogonal matrices. As use cases, we compute Hermite interpolation curves for orthogonal matrix factorizations such as the singular value decomposition and the QR-decomposition., Comment: 27 pages, 13 figures

One approach to parametric and adaptive model reduction is via the interpolation of orthogonal bases, subspaces or positive definite system matrices. In all these cases, the sampled inputs stem from matrix sets that feature a geometric structure and thus form so-called matrix manifolds. This work will be featured as a chapter in the upcoming Handbook on Model Order Reduction (P. Benner, S. Grivet-Talocia, A. Quarteroni, G. Rozza, W.H.A. Schilders, L.M. Silveira, eds, to appear on DE GRUYTER) and reviews the numerical treatment of the most important matrix manifolds that arise in the context of model reduction. Moreover, the principal approaches to data interpolation and Taylor-like extrapolation on matrix manifolds are outlined and complemented by algorithms in pseudo-code., Comment: 37 pages, 4 figures, featured chapter of upcoming "Handbook on Model Order Reduction"

In this work, we consider rank-one adaptations $X_{new} = X+ab^T$ of a given matrix $X\in \mathbb{R}^{n\times p}$ with known matrix factorization $X = UW$, where $U\in\mathbb{R}^{n\times p}$ is column-orthogonal, i.e. $U^TU=I$. Arguably the most important methods that produce such factorizations are the singular value decomposition (SVD), where $X=UW=U\Sigma V^T$, and the QR-decomposition, where $X = UW = QR$. An elementary approach to produce a column-orthogonal matrix $U_{new}$, whose columns span the same subspace as the columns of the rank-one modified $X_{new} = X +ab^T$ is via applying a suitable coordinate change such that in the new coordinates, the update affects a single column and subsequently performing a Gram-Schmidt step for reorthogonalization. This may be interpreted as a rank-one adaptation of the $U$-factor in the SVD or a rank-one adaptation of the $Q$-factor in the QR-decomposition, respectively, and leads to a decomposition for the adapted matrix $X_{new} = U_{new}W_{new}$. By using a geometric approach, we show that this operation is equivalent to traveling from the subspace $\mathcal{S}= \text{ran}(X)$ to the subspace $\mathcal{S}_{new} =\text{ran}(X_{new})$ on a geodesic line on the Grassmann manifold and we derive a closed-form expression for this geodesic. In addition, this allows us to determine the subspace distance between the subspaces $\mathcal{S}$ and $\mathcal{S}_{new}$ without additional computational effort. Both $U_{new}$ and $W_{new}$ are obtained via elementary rank-one matrix updates in $\mathcal{O}(np)$ time for $n\gg p$., Comment: 15 pages, 1 figure, MATLAB code