High-Order Retractions on Matrix Manifolds using Projected Polynomials

TL;DR

This paper introduces high-order, structure-preserving approximations of the Riemannian exponential map on matrix manifolds using polynomial-based methods inspired by Bessel polynomials, improving accuracy for manifold computations.

Contribution

It develops a novel family of high-order approximation techniques for the exponential map on matrix manifolds, leveraging polynomial approximations and polar decompositions.

Findings

Achieves error of order O(t^{2n+1}) in approximations

Provides high-order methods for Grassmannian and Stiefel manifolds

Establishes supercloseness of geometric and arithmetic means of unitary matrices

Abstract

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 $\Omega$ is a skew-Hermitian matrix and $t$ is a sufficiently small scalar, then there exists a polynomial of degree $n$ in $t\Omega$ (namely, a Bessel polynomial) whose polar decomposition delivers an approximation of $e^{t\Omega}$ 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.

Tables

Table 1. Table 1: Errors in approximating the exponential of a skew-Hermitian matrix Ω Ω \Omega with the projected polynomials of Theorem  2.1 . Shown above are the errors ‖ 𝒫 ​ ( Θ n ​ ( t ​ Ω ) ) − e t ​ Ω ‖ norm 𝒫 subscript Θ 𝑛 𝑡 Ω superscript 𝑒 𝑡 Ω \|\mathcal{P}(\Theta_{n}(t\Omega))-e^{t\Omega}\| versus t 𝑡 t for n = 1 , 2 , 3 𝑛 1 2 3 n=1,2,3 , where t 0 = 0.01 subscript 𝑡 0 0.01 t_{0}=0.01 and Ω Ω \Omega is a random 1000 × 1000 1000 1000 1000\times 1000 skew-Hermitian matrix.

	$n=1$			$n=2$			$n=3$
$t_0/t$	Error		Order	Error		Order	Error		Order
$1$	$1$	$.607\cdot {10}^0$		$6$	$.945\cdot {10}^{-2}$		$1$	$.312\cdot {10}^{-3}$
$2$	$2$	$.433\cdot {10}^{-1}$	$2.724$	$2$	$.444\cdot {10}^{-3}$	$4.828$	$1$	$.109\cdot {10}^{-5}$	$6.887$
$4$	$3$	$.223\cdot {10}^{-2}$	$2.916$	$7$	$.860\cdot {10}^{-5}$	$4.959$	$8$	$.830\cdot {10}^{-8}$	$6.972$
$8$	$4$	$.091\cdot {10}^{-3}$	$2.978$	$2$	$.474\cdot {10}^{-6}$	$4.990$	$6$	$.932\cdot {10}^{-10}$	$6.993$

Table 2. Table 2: Errors in approximating the Riemannian exponential map on the Grassmannian manifold with the projected polynomials of Theorem  2.3 . Shown above are the errors ‖ 𝒫 ​ ( Y ​ α n ​ ( t 2 ​ H ∗ ​ H ) + t ​ H ​ β n ​ ( t 2 ​ H ∗ ​ H ) ) − Exp Y G ​ r ​ ( t ​ H ) ‖ norm 𝒫 𝑌 subscript 𝛼 𝑛 superscript 𝑡 2 superscript 𝐻 𝐻 𝑡 𝐻 subscript 𝛽 𝑛 superscript 𝑡 2 superscript 𝐻 𝐻 subscript superscript Exp 𝐺 𝑟 𝑌 𝑡 𝐻 \|\mathcal{P}(Y\alpha_{n}(t^{2}H^{*}H)+tH\beta_{n}(t^{2}H^{*}H))-\mathrm{Exp}^{Gr}_{Y}(tH)\| versus t 𝑡 t for n = 1 , 2 , 3 𝑛 1 2 3 n=1,2,3 , where t 0 = 0.01 subscript 𝑡 0 0.01 t_{0}=0.01 , Y 𝑌 Y is a random 2000 × 400 2000 400 2000\times 400 matrix with orthonormal columns, and H 𝐻 H is a random 2000 × 400 2000 400 2000\times 400 matrix satisfying Y ∗ ​ H = 0 superscript 𝑌 𝐻 0 Y^{*}H=0 .

	$n=1$			$n=2$			$n=3$
$t_0/t$	Error		Order	Error		Order	Error		Order
$1$	$5$	$.030\cdot {10}^{-1}$		$9$	$.383\cdot {10}^{-3}$		$7$	$.849\cdot {10}^{-5}$
$2$	$6$	$.951\cdot {10}^{-2}$	$2.855$	$3$	$.090\cdot {10}^{-4}$	$4.924$	$6$	$.352\cdot {10}^{-7}$	$6.949$
$4$	$8$	$.933\cdot {10}^{-3}$	$2.960$	$9$	$.781\cdot {10}^{-6}$	$4.982$	$5$	$.006\cdot {10}^{-9}$	$6.987$
$8$	$1$	$.125\cdot {10}^{-3}$	$2.990$	$3$	$.066\cdot {10}^{-7}$	$4.995$	$3$	$.919\cdot {10}^{-11}$	$6.997$

Table 3. Table 3: Errors in approximating the Riemannian exponential map on the Stiefel manifold with the projected polynomials of Theorem  2.5 . Shown above are the errors ‖ 𝒫 ​ ( Y ​ γ n ​ ( t 2 ​ H ∗ ​ H , t ​ Y ∗ ​ H ) + t ​ H ​ δ n ​ ( t 2 ​ H ∗ ​ H , t ​ Y ∗ ​ H ) ) − Exp Y S ​ t ​ ( t ​ H ) ‖ norm 𝒫 𝑌 subscript 𝛾 𝑛 superscript 𝑡 2 superscript 𝐻 𝐻 𝑡 superscript 𝑌 𝐻 𝑡 𝐻 subscript 𝛿 𝑛 superscript 𝑡 2 superscript 𝐻 𝐻 𝑡 superscript 𝑌 𝐻 subscript superscript Exp 𝑆 𝑡 𝑌 𝑡 𝐻 \|\mathcal{P}(Y\gamma_{n}(t^{2}H^{*}H,tY^{*}H)+tH\delta_{n}(t^{2}H^{*}H,tY^{*}H))-\mathrm{Exp}^{St}_{Y}(tH)\| versus t 𝑡 t for n = 1 , 2 , 3 𝑛 1 2 3 n=1,2,3 , where t 0 = 0.01 subscript 𝑡 0 0.01 t_{0}=0.01 , Y 𝑌 Y is a random 2000 × 400 2000 400 2000\times 400 matrix with orthonormal columns, and H 𝐻 H is a random 2000 × 400 2000 400 2000\times 400 matrix satisfying Y ∗ ​ H = − H ∗ ​ Y superscript 𝑌 𝐻 superscript 𝐻 𝑌 Y^{*}H=-H^{*}Y .

	$n=1$			$n=2$			$n=3$
$t_0/t$	Error		Order	Error		Order	Error		Order
$1$	$1$	$.336\cdot {10}^0$		$1$	$.416\cdot {10}^{-1}$		$1$	$.654\cdot {10}^{-2}$
$2$	$3$	$.013\cdot {10}^{-1}$	$2.149$	$1$	$.931\cdot {10}^{-2}$	$2.875$	$9$	$.724\cdot {10}^{-4}$	$4.089$
$4$	$7$	$.195\cdot {10}^{-2}$	$2.066$	$2$	$.479\cdot {10}^{-3}$	$2.962$	$5$	$.930\cdot {10}^{-5}$	$4.035$
$8$	$1$	$.774\cdot {10}^{-2}$	$2.020$	$3$	$.120\cdot {10}^{-4}$	$2.990$	$3$	$.681\cdot {10}^{-6}$	$4.010$

Equations

$$q_n(t\Omega )=\sum _{k=0}^n\frac{(t\Omega {)}^k}{k!}.$$

$$\Vert U-q_n(t\Omega )\Vert \le \Vert V-q_n(t\Omega )\Vert$$

$$\Vert U-e^{t\Omega }\Vert$$

$$\le \Vert e^{t\Omega }-q_n(t\Omega )\Vert +\Vert q_n(t\Omega )-e^{t\Omega }\Vert$$

$$=O(t^{n+1}).$$

$$q_n(t\Omega )=\sum _{k=0}^n{a}_k{t}^k{\Omega }^k$$

$${\Theta }_n(z)=\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)\frac{(2n-k)!}{(2n)!}(2z{)}^k.$$

$$\mathcal{P}({\Theta }_n(t\Omega ))=e^{t\Omega }+O(t^{2n+1}).$$

$$\mathcal{P}(r_n(t\Omega ))=e^{t\Omega }+O(t^k)$$

$${\Theta }_n(z)=\frac{2^{n}n!}{(2n)!}{\theta }_n(z),$$

$${\theta }_n(z)=\sum _{k=0}^n\frac{(n+k)!}{(n-k)!k!}\frac{z^{n-k}}{2^k}$$

$${\Theta }_0(z)$$

$${\Theta }_1(z)$$

$${\Theta }_2(z)$$

$${\Theta }_3(z)$$

$${\Theta }_4(z)$$

$$T_{Y}Gr(p,m)=\{Y_{\perp }K\mid K\in {\mathbb{C}}^{(m-p)\times p}\},$$

$${\mathrm{Exp}}_Y^{Gr}H=\left(\begin{array}{cc}Y& Y_{\perp }\end{array}\right)\exp \left(\begin{array}{cc}0& -K^{\ast }\\ K& 0\end{array}\right)\left(\begin{array}{c}I\\ 0\end{array}\right).$$

$${\mathrm{Exp}}_Y^{Gr}H=YV\cos (\Sigma )V^{\ast }+U\sin (\Sigma )V^{\ast }.$$

$${\mathrm{Exp}}_Y^{Gr}(tH)=\left(\begin{array}{cc}Y& Y_{\perp }\end{array}\right)\mathcal{P}({\Theta }_n(tZ))\left(\begin{array}{c}I\\ 0\end{array}\right)+O(t^{2n+1}),$$

$$\mathcal{P}({\Theta }_n(tZ))\left(\begin{array}{c}I\\ 0\end{array}\right)=\mathcal{P}\left({\Theta }_n(tZ)\left(\begin{array}{c}I\\ 0\end{array}\right)\right)$$

$${\mathrm{Exp}}_Y^{Gr}(tH)=\mathcal{P}\left(Y{\alpha }_n(t^2{H}^{\ast }H)+tH{\beta }_n(t^2{H}^{\ast }H)\right)+O\left(t^{2n+1}\right),$$

$${\alpha }_n(z)$$

$${\beta }_n(z)$$

$$a_k=\left(\genfrac{}{}{0ex}{}{n}{k}\right)\frac{(2n-k)!}{(2n)!}{2}^k,k=0,1,\dots ,n.$$

$${\mathrm{Exp}}_Y^{Gr}(tH)$$

$$\begin{array}{rl}{\mathrm{Exp}}_Y^{Gr}(tH)& =\mathcal{P}(Y\left(I-\frac{3}{7}{t}^2{H}^{\ast }H+\frac{1}{105}{t}^4(H^{\ast }H{)}^2\right)\\ & +tH\left(I-\frac{2}{21}{t}^2{H}^{\ast }H\right))+O(t^9).\end{array}$$

$$\mathcal{P}\left(\tilde{Y}{\alpha }_n(t^2{\tilde{H}}^{\ast }\tilde{H})+t\tilde{H}{\beta }_n(t^2{\tilde{H}}^{\ast }\tilde{H})\right)$$

$$=V\mathcal{P}\left(Y{\alpha }_n(t^2{H}^{\ast }H)+tH{\beta }_n(t^2{H}^{\ast }H)\right),$$

$$\mathcal{P}(VA)=V\mathcal{P}(A)$$

$${\mathrm{dist}}^{Gr}(X,Y)=\underset{\begin{array}{c}V,W\in {\mathbb{C}}^{p\times p}\\ V^{\ast }V=W^{\ast }W=I\end{array}}{\min }\Vert XV-YW\Vert ,$$

$${\mathrm{dist}}^{Gr}(\mathcal{Q}\left(Y{\alpha }_n(t^2{H}^{\ast }H)+tH{\beta }_n(t^2{H}^{\ast }H)\right),{\mathrm{Exp}}_Y(tH))=O(t^{2n+1}),$$

$$St(p,m)=\{Y\in {\mathbb{C}}^{m\times p}\mid Y^{\ast }Y=I\}.$$

$$T_{Y}St(p,m)=\{Y\Omega +Y_{\perp }K\mid \Omega =-{\Omega }^{\ast }\in {\mathbb{C}}^{p\times p},K\in {\mathbb{C}}^{(m-p)\times p}\},$$

Full text

1 Introduction

Approximating the Riemannian or Lie-theoretic exponential map on a matrix manifold is a task of importance in a variety of applications, including numerical integration on Lie groups [17, 23, 18, 5], optimization on manifolds [2, 8, 4, 30], interpolation of manifold-valued data [33, 34, 16, 13], rigid body simulation [5, 26], fluid simulation [15], and computer vision [36, 29, 11]. Often, special attention is paid to preserving the structure of the exponential map [6], which, for instance, should return a unitary matrix when the input $\Omega$ is skew-Hermitian. In this paper, we construct structure-preserving approximations to the Riemannian exponential map on matrix manifolds using projected polynomials – polynomial functions of matrices which, when projected onto a suitable set, deliver approximations to the Riemannian exponential with a desired order of accuracy. These projected polynomials can be thought of as high-order generalizations of the “projection-like retractions” considered in [3]. The matrix manifolds we consider are:

The group of unitary $m\times m$ matrices. 2. 2.

The Grassmannian manifold $Gr(p,m)$, which consists of all $p$-dimensional linear subspaces of $\mathbb{C}^{m}$, where $m\geq p$. 3. 3.

The Stiefel manifold $St(p,m)=\{Y\in\mathbb{C}^{m\times p}\mid Y^{*}Y=I\}$, where $m\geq p$.

The projector we use is to accomplish this task is the map which sends a full-rank matrix $A\in\mathbb{C}^{m\times p}$ ($m\geq p$) to the nearest matrix with orthonormal columns. The latter matrix is precisely the factor $U$ in the polar decomposition $A=UH$, where $U\in\mathbb{C}^{m\times p}$ has orthonormal columns and $H\in\mathbb{C}^{p\times p}$ is Hermitian positive-definite [9, Theorem 1]. In the case of the Grassmannian manifold, the $QR$ decomposition can be used in place of the polar decomposition, leading to methods with very low computational cost.

Interestingly, in the case of the unitary group and the Grassmannian manifold, superconvergent approximations of the exponential are constructible with this approach. By this we mean that it is possible to construct polynomials of degree $n$ that, upon projection, deliver approximations to the exponential with error of order $n+2$ or higher. The appropriate choices of polynomials turn out to be intimately related to the Bessel polynomials, a well-known orthogonal sequence of polynomials [27], and the resulting approximations have error of order $2n+1$; see Theorems 2.1 and 2.3 and Corollary 2.4.

One of the major advantages of this approach is that it delivers approximations to the exponential on the unitary group, the Grassmannian manifold, and the Steifel manifold that, to machine precision, have orthonormal columns. This is of obvious importance for the unitary group and the Stiefel manifold, and it is even desirable on the Grassmannian manifold, where it is common in computations to represent elements of the Grassmannian – $p$-dimensional subspaces of $\mathbb{C}^{m}$ – as $m\times p$ matrices whose columns form orthonormal bases for those subspaces [8].

Furthermore, when the polar decomposition is adopted as the projector, projected polynomials have the advantage that they can be computed using only rudimentary operations on matrices: matrix addition, multiplication, and inversion. This follows from the fact that the polar decomposition can be computed iteratively [20, Chapter 8]. Rudimentary algorithms for calculating the exponential on the Grassmannian and Stiefel manifolds are particularly desirable, since the most competitive existing algorithms for accomplishing this task involve singular value and/or eigenvalue decompositions [8, Theorem 2.1, Corollary 2.2, and Theorem 2.3], a feature that renders existing algorithms less ideally suited for parallel computation than projected polynomials.

In spite of these advantages, it is worth noting that not all of the constructions in this paper lead to algorithms that outshine their competitors. Diagonal Padé approximations of the exponential deliver, to machine precision, unitary approximations of $e^{\Omega}$ when $\Omega$ is skew-Hermitian [23, p. 97]. It is clear that the projected polynomials we present below (in Theorem 2.1) for approximating $e^{\Omega}$ are more expensive to compute, at least when the comparison is restricted to approximations of $e^{\Omega}$ with equal order of accuracy. On the Stiefel manifold, Padé approximation is not an option, rendering projected polynomials more attractive. However, they are not superconvergent on the Stiefel manifold; see Theorem 2.5. The setting in which projected polynomials appear to shine the brightest is the Grassmannian manifold $Gr(p,m)$, where they provide superconvergent, orthonormal approximations to the Riemannian exponential with algorithmic complexity $O(mp^{2})$; see Theorem 2.3 and Corollary 2.4. To our knowledge, these are the first such approximations (other than the lowest-order versions) to appear in the literature on the Grassmannian manifold.

Structure-preserving approximations of the exponential map on matrix manifolds have a long history, particularly for matrix manifolds that form Lie groups. On Lie groups, techniques involving rational approximation [23, p. 97], splitting [6, 37], canonical coordinates of the second kind [7], and the generalized polar decomposition [24] have been studied, and many of these strategies lead to high-order approximations. For more general matrix manifolds like the Grassmannian and Stiefel manifolds, attention has been primarily restricted to methods for calculating the exponential exactly [8, 1, 2] or approximating it to low order [2, 3, 25, 10]. In this context, structure-preserving approximations of the exponential are commonly referred to as retractions [2, Definition 4.1.1]. High-order retractions on the Grassmannian and Stiefel manifolds have received very little attention, but there are good reasons to pursue them. For instance, in optimization, exactly evaluating the Riemannian Hessian of a function defined on a matrix manifold requires the use of a retraction with second-order accuracy or higher, at least if one is interested in its value away from critical points [2, p. 107]. In addition, existing algorithms for calculating the exponential on the Grassmannian and Stiefel manifolds exactly [8, Theorem 2.1, Corollary 2.2, and Theorem 2.3] are relatively expensive, which raises the question of whether more efficient options, perhaps with nonzero but controllable error, are available. On the Grassmannian manifold, the answer seems to be yes, at least for small-normed input; see Corollary 2.4.

Organization

This paper is organized as follows. In Section 2, we give statements of our results, deferring their proofs to Section 3. Our main results are Theorems 2.1, 2.3, and 2.5, which detail families of approximants to the exponential on the unitary group, the Grassmannian manifold, and the Steifel manifold, respectively. A fourth noteworthy result is Corollary 2.4, which provides a computationally inexpensive variant of the approximants in Theorem 2.3. We also detail two related results, Proposition 2.6 and Theorem 2.7, that concern the supercloseness of the geometric and arithmetic means of unitary matrices. In Section 3, we prove each of the results just mentioned. In Section 4, we describe algorithms for calculating our proposed approximations, with an emphasis on iterative methods for computing the polar decomposition. We conclude that section with numerical examples.

2 Statement of Results

In this section, we give statements of our results. Proofs are detailed in Section 3.

2.1 Exponentiation on the Unitary Group

Our first result deals with the approximation of the exponential of a skew-Hermitian matrix $\Omega\in\mathbb{C}^{m\times m}$ with projected polynomials. To motivate the forthcoming theorem, consider the Taylor polynomial $q_{n}(t\Omega)$ of degree $n$ for $e^{t\Omega}$:

[TABLE]

This quantity, in general, is not an unitary matrix, even though the matrix it aims to approximate, $e^{t\Omega}$, is unitary. If $t$ is sufficiently small (small enough so that $q_{n}(t\Omega)$ is nonsingular), then $q_{n}(t\Omega)$ can be made unitary by computing the polar decomposition $q_{n}(t\Omega)=UH$, where $U\in\mathbb{C}^{m\times m}$ is unitary and $H\in\mathbb{C}^{m\times m}$ is Hermitian positive-definite. The matrix $U$ is easily seen to be an unitary approximation to $e^{t\Omega}$ with error at worst $O(t^{n+1})$, owing to the fact that

[TABLE]

for every unitary matrix $V\in\mathbb{C}^{m\times m}$, where $\|\cdot\|$ denotes the Frobenius norm [9, Theorem 1]. Indeed,

[TABLE]

Below we address the question of whether a better approximation to $e^{t\Omega}$ can be constructed by computing the unitary factor in the polar decomposition of

[TABLE]

for suitably chosen coefficients $a_{k}$. We show that if the coefficients $a_{k}$ are chosen carefully, then an approximation with error of order $t^{2n+1}$ can be constructed. The choice of coefficients $a_{k}$ corresponds to the selection of a Bessel polynomial of degree $n$ in $t\Omega$. In what follows, we use $\mathcal{P}$ to denote the map which sends a full-rank matrix $A\in\mathbb{C}^{m\times p}$ ($m\geq p$) to the factor $\mathcal{P}(A)=U$ in its polar decomposition $A=UH$, where $U\in\mathbb{C}^{m\times p}$ has orthonormal columns and $H\in\mathbb{C}^{p\times p}$ is Hermitian positive-definite [20, Theorem 8.1].

Theorem 2.1**.**

Let $\Omega\in\mathbb{C}^{m\times m}$ be skew-Hermitian, and let $n\geq 0$ be an integer. Define

[TABLE]

Then

[TABLE]

In fact, we will show that the polar decomposition of $\Theta_{n}(t\Omega)$ delivers the highest order approximation of $e^{t\Omega}$ among all polynomials in $t\Omega$ of degree $n$, up to rescaling. That is, if $r_{n}(t\Omega)$ is any other polynomial in $t\Omega$ of degree $n$ satisfying $r_{n}(0)=I$, then

[TABLE]

for some $0\leq k\leq 2n$.

The polynomials (1) are scaled versions of Bessel polynomials [27]. More precisely, we have

[TABLE]

where

[TABLE]

denotes the reverse Bessel polynomial of degree $n$. The first few polynomials $\Theta_{n}(z)$ are given by

[TABLE]

Note that, rather surprisingly, $\Theta_{n}(z)$ agrees with $e^{z}$ only to first order for every $n\geq 1$.

2.2 Exponentiation on the Grassmannian

We now consider the task of approximating the Riemannian exponential map on the Grassmannian manifold $Gr(p,m)$, which consists of all $p$-dimensional subspaces of $\mathbb{C}^{m}$, where $m\geq p$. We begin by reviewing the geometry of $Gr(p,m)$, with an emphasis on computational aspects.

In computations, it is convenient to represent each subspace $\mathcal{V}\in Gr(p,m)$ with a matrix $Y\in\mathbb{C}^{m\times p}$ having orthonormal columns that span $\mathcal{V}$. The choice of $Y$ is not unique, so we are of course thinking of $Y$ as a representative of an equivalence class of $m\times p$ matrices sharing the same column space. With this identification, the tangent space to $Gr(p,m)$ at $Y$ is given by

[TABLE]

where $Y_{\perp}\in\mathbb{C}^{m\times(m-p)}$ is any matrix such that $\begin{pmatrix}Y&Y_{\perp}\end{pmatrix}$ is unitary, and $Y_{\perp}K$ is regarded as a representative of an equivalence class of matrices sharing the same column space [8, p. 15]. With respect to the canonical metric on $Gr(p,m)$, the Riemannian exponential $\mathrm{Exp}^{Gr}_{Y}:T_{Y}Gr(p,m)\rightarrow Gr(p,m)$ at $Y\in\mathbb{C}^{m\times p}$ in the direction $H=Y_{\perp}K\in\mathbb{C}^{m\times p}$ is given by [8, p. 10]

[TABLE]

The goal of this subsection is to construct computationally inexpensive approximations of $\text{Exp}^{Gr}_{Y}H$. In order to be competitive with existing methods, such approximations must have computational complexity $O(mp^{2})$ or better, owing to the following well-known result [8, Theorem 2.3].

Theorem 2.2**.**

[8, Theorem 2.3]** Let $H=U\Sigma V^{*}$ be the thin singular value decomposition of $H$, i.e. $U\in\mathbb{C}^{m\times p}$ has orthonormal columns, $\Sigma\in\mathbb{C}^{p\times p}$ is diagonal with nonnegative entries, and $V\in\mathbb{C}^{p\times p}$ is unitary. Then

[TABLE]

The preceding theorem reveals that $\mathrm{Exp}^{Gr}_{Y}H$ can be computed exactly with $O(mp^{2})$ operations, since this is the cost of computing the thin singular value decomposition of $H\in\mathbb{C}^{m\times p}$. With this in mind, we aim to derive approximations of $\mathrm{Exp}^{Gr}_{Y}H$ with smaller or comparable computational complexity.

Since the matrix $Z:=\begin{pmatrix}0&-K^{*}\\ K&0\end{pmatrix}$ appearing in (2) is skew-Hermitian, an obvious option is to approximate $\exp Z$ in (2) with a projected polynomial $\mathcal{P}(\Theta_{n}(Z))$ in accordance with Section 2.1. This leads to approximants of the form

[TABLE]

which, unfortunately, have computational complexity $O(m^{3})$. Remarkably, we show in Lemma 3.14 below that

[TABLE]

if $t$ is sufficiently small (small enough so that $\Theta_{n}(tZ)$ is nonsingular). This is significant, since the right-hand side of this equality involves the polar decomposition of an $m\times p$ matrix $\Theta_{n}(tZ)\begin{pmatrix}I\\ 0\end{pmatrix}$, which can be computed in $O(mp^{2})$ operations. A few more algebraic manipulations (detailed in Section 3.4) lead to the following scheme for approximating the exponential on the Grassmannian in $O(mp^{2})$ operations.

Theorem 2.3**.**

Let $Y\in\mathbb{C}^{m\times p}$ have orthonormal columns, and let $H\in T_{Y}Gr(p,m)$. Then, for any $n\geq 0$,

[TABLE]

where

[TABLE]

and

[TABLE]

The first few nontrivial approximants provided by Theorem 2.3 read

[TABLE]

Note that, rather interestingly, the commonly used retraction $\mathcal{P}(Y+tH)$ (see, for instance, [1]) is in fact an approximation of $\mathrm{Exp}^{Gr}_{Y}(tH)$ with error $O(t^{3})$, despite its appearance.

Note also that the approximants provided by Theorem 2.3 are rotationally equivariant. That is, if $V\in\mathbb{C}^{m\times m}$ is a unitary matrix, $\widetilde{Y}=VY$, and $\widetilde{H}=VH$, then $\widetilde{H}^{*}\widetilde{H}=H^{*}V^{*}VH=H^{*}H$ and hence

[TABLE]

where the last line follows from the fact that

[TABLE]

for every full-rank $A\in\mathbb{C}^{m\times p}$ ($m\geq p$) and every unitary $V\in\mathbb{C}^{m\times m}$.

Replacing the Polar Decomposition with the QR Decomposition

An extraordinary feature of Theorem 2.3 is that it applies, with slight modification, even if the map $\mathcal{P}$ is replaced by the map $\mathcal{Q}$ which sends a full-rank matrix $A\in\mathbb{C}^{m\times p}$ ($m\geq p$) to the factor $Q$ in the QR decomposition $A=QR$, where $Q\in\mathbb{C}^{m\times p}$ has orthonormal columns and $R\in\mathbb{C}^{p\times p}$ is upper triangular. (The map $\mathcal{Q}$ is denoted $\mathrm{qf}$ in [2].) This follows from the fact that the columns of $A$, $\mathcal{P}(A)$, and $\mathcal{Q}(A)$ span the same space, so $\mathcal{P}(A)$ and $\mathcal{Q}(A)$ represent the same element of the Grassmannian manifold.

The only modification of Theorem 2.3 needed to make this idea precise is to use a genuine distance on the Grassmannian, such as

[TABLE]

to measure the distance between subspaces [8, p. 30]. We then have the following corollary.

Corollary 2.4**.**

Let $Y\in\mathbb{C}^{m\times p}$ have orthonormal columns, and let $H\in T_{Y}Gr(p,m)$. Then, for any $n\geq 0$,

[TABLE]

*where $\alpha_{n}$ and $\beta_{n}$ are given in the statement of Theorem 2.3. *

This corollary is quite powerful, since the factor $Q$ in the QR decomposition of an $m\times p$ matrix can be computed in merely $2mp^{2}-\frac{2}{3}p^{3}$ operations [20, Appendix C], rendering the approximations $\mathcal{Q}\left(Y\alpha_{n}(t^{2}H^{*}H)+tH\beta_{n}(t^{2}H^{*}H)\right)$ very cheap to compute. Note also that these approximations are rotationally equivariant, by an argument similar to the one leading up to (9).

2.3 Exponentiation on the Stiefel Manifold

The next manifold we consider is the Stiefel manifold

[TABLE]

Unlike the Grassmannian, here we do not regard matrices $Y\in\mathbb{C}^{m\times p}$ as representatives of equivalence classes; each $Y\in\mathbb{C}^{m\times p}$ corresponds to a distinct element of $St(p,m)$. The tangent space to $St(p,m)$ at $Y\in St(p,m)$ is given by

[TABLE]

where $Y_{\perp}\in\mathbb{C}^{m\times(m-p)}$ is any matrix such that $\begin{pmatrix}Y&Y_{\perp}\end{pmatrix}$ is unitary [8, Equation 2.5]. With respect to the canonical metric on $St(p,m)$ [8, Section 2.4], the Riemannian exponential $\mathrm{Exp}^{St}_{Y}:T_{Y}St(p,m)\rightarrow St(p,m)$ at $Y\in\mathbb{C}^{m\times p}$ in the direction $H=Y\Omega+Y_{\perp}K\in\mathbb{C}^{m\times p}$ is given by [8, Equation 2.42]

[TABLE]

As an aside, we remark that a different exponential map is obtained if one endows $St(p,m)$ with the metric inherited from the embedding of $St(p,m)$ in Euclidean space [8, Section 2.2]. We do not consider the latter exponential map in this paper.

There exist algorithms for calculating (11) in $O(mp^{2})$ operations, the simplest of which involves calculating the $QR$ decomposition of a certain $m\times p$ matrix and then exponentiating a $2p\times 2p$ skew symmetric matrix [8, Corollary 2.2]. Our aim below is to derive a competitive algorithm for approximating (11) to high order using projected polynomials.

Before doing so, it is important to note that the right-hand side of (11) reduces to more familiar expressions in two special cases. First, when $Y^{*}H=\Omega=0$, the right-hand side of (11) coincides with the right-hand side of (2), the Riemmannian exponential on the Grassmannian. On the other hand, if $m=p$ and $Y=I$, then $K$ and $Y_{\perp}$ are empty matrices and the right-hand side of (11) reduces to $e^{\Omega}$, the exponential of a skew-Hermitian matrix. Thus, in an effort to generalize Theorems 2.1 and 2.3, we seek to approximate (11) with projected polynomials that reduce to the ones appearing in Theorems 2.1 and 2.3 in those special cases.

In view of Theorem 2.3 and the identities $Y^{*}H=\Omega$ and $H^{*}H=-\Omega^{2}+K^{*}K$, it is natural to consider approximations of (11) of the form

[TABLE]

where $q(x,y)$ and $r(x,y)$ are polynomials in the (non-commuting) variables $x$ and $y$. In order to ensure that these approximations recover those appearing in Theorems 2.1 and 2.3, we insist that:

1. (LABEL:enumisec:stiefelexp.i)

$q(x,0)=\alpha_{n}(x)$. 2. (LABEL:enumisec:stiefelexp.ii)

$r(x,0)=\beta_{n}(x)$. 3. (LABEL:enumisec:stiefelexp.iii)

$q(-x^{2},x)$ and $r(-x^{2},x)$ are polynomials of degree at most $n$ and $n-1$, respectively, satisfying

[TABLE]

It turns out that such approximations can be constructed, but they lack the superconvergence enjoyed by the approximations in Theorems 2.1 and 2.3 (unless $\Omega=0$ or $m=p$). The difficulty becomes apparent if one compares $\mathcal{P}(Y+tH)$ with $\mathrm{Exp}^{St}_{Y}(tH)$ for generic $Y\in St(p,m)$ and $H=Y\Omega+Y_{\perp}K\in T_{Y}St(p,m)$. As $t\rightarrow 0$, one observes numerically that $\mathcal{P}(Y+tH)=\mathrm{Exp}^{St}_{Y}(tH)+O(t^{2})$ (unless $\Omega=0$ or $m=p$).111The astute reader may notice that this appears to contradict Theorem 4.9 of [3], which states, among other things, that projective retractions (see [3, Example 4.5]) are automatically second-order. However, it is not the exponential map (11) that $\mathcal{P}(Y+tH)$ approximates to second order. Rather, it is exponential map associated with the metric inherited from the embedding of $St(p,m)$ in $\mathbb{C}^{m\times p}$. This contrasts starkly with the situation in Theorems 2.1 and 2.3, where the polar decomposition of the first-order Taylor approximant of the exponential had superconvergent error $O(t^{3})$. The following theorem confirms this observation and provides a couple of higher-order approximations of (11). In it, we use $\|\cdot\|$ to denote the Frobenius norm.

Theorem 2.5**.**

Let $Y\in St(p,m)$ and $H\in T_{Y}St(p,m)$. Define

[TABLE]

For $n=1,2,3$, we have

[TABLE]

where

[TABLE]

In addition, for every polynomial $q(x,y)$ and $r(x,y)$ satisfying (sec:stiefelexp.i-sec:stiefelexp.iii) ($1\leq n\leq 3$), there exists $Y\in St(p,m)$, $H\in T_{Y}St(p,m)$, $C>0$, and $t_{0}>0$ such that

[TABLE]

*for every $t\leq t_{0}$. *

Written more explicitly, the approximants provided by Theorem 2.5 read

[TABLE]

All of these are rotationally equivariant by an argument similar to the one leading up to (9).

Observe that when $Y^{*}H=\Omega=0$, the right-hand sides of (16-18) reduce to those in (6-7), respectively. Likewise, when $m=p$ (so that $H=Y\Omega$ and $YY^{*}=I$), they reduce to $Y\mathcal{P}(\Theta_{1}(tH))$, $Y\mathcal{P}(\Theta_{2}(tH))$, and $Y\mathcal{P}(\Theta_{3}(tH))$, respectively, which are precisely the approximations of $Ye^{\Omega}$ provided by Theorem 2.1. See Section 3.5 for details.

In view of the complexity of (17) and (18), we do not believe that a general formula (valid for all $n$) can be conveniently written down for polynomials $q(x,y)$ and $r(x,y)$ satisfying (sec:stiefelexp.i-sec:stiefelexp.iii) that deliver approximations of (11) of the form (12) with error of optimal order. (As a matter of fact, the polynomials $\gamma_{3}(x,y)$ and $\delta_{3}(x,y)$ are not even uniquely determined by these conditions.) However, the proofs presented in Section 3.5 demonstrate how one can construct such polynomials.

Note that if the conditions (sec:stiefelexp.i-sec:stiefelexp.iii) are relaxed, then it is straightforward to construct approximations of (11) with error $O(t^{n+1})$: Simply truncate the Taylor series for $\exp\begin{pmatrix}\Omega&-K^{*}\\ K&0\end{pmatrix}$, insert the result into (11), and express the result in terms of $Y$ and $H$ using the identities $H=Y\Omega+Y_{\perp}K$, $Y^{*}H=\Omega$, and $H^{*}H=-\Omega^{2}+K^{*}K$. If desired, the result can be orthonormalized with the map $\mathcal{P}$, retaining the order of accuracy of the approximation. For this reason, Theorem 2.5 is less powerful than Theorems 2.1 and 2.3, and it underscores the complexity of the Stiefel manifold relative to the Grassmannian and the unitary group. For more evidence of the computational difficulties inherent to the Stiefel manifold, we refer the reader to [8].

2.4 Geometric and Arithmetic Means of Unitary Matrices

We conclude this section by stating two results that concern the supercloseness of certain means of unitary matrices. At the surface, these results might not appear to be closely related to Theorems 2.1, 2.3, 2.5, but in fact they follow from the same general theory.

Our first result reveals that the polar decomposition of the (componentwise) linear interpolant of two unitary matrices $U_{1}$ and $U_{2}$ is superclose to the geodesic joining $U_{1}$ and $U_{2}$.

Proposition 2.6**.**

Let $U_{1}\in\mathbb{C}^{m\times m}$ be unitary, let $\Omega\in\mathbb{C}^{m\times m}$ be skew-Hermitian, and let $U_{2}=U_{1}e^{t\Omega}$. Assume that $(1-s)U_{1}+sU_{2}$ is nonsingular for each $s\in[0,1]$. Then

[TABLE]

*for every $s\in[0,1/2)\cup(1/2,1]$. When $s=1/2$, the equality $\mathcal{P}((1-s)U_{1}+sU_{2})=U_{1}e^{st\Omega}$ holds exactly. *

The case $s=1/2$ in the preceding lemma recovers the well-known observation (see [21, Theorem 4.7] and [31, Equation (3.14)]) that the unitary factor $U=\mathcal{P}(\frac{1}{2}(U_{1}+U_{2}))$ in the polar decomposition of $\frac{1}{2}(U_{1}+U_{2})$ is given by

[TABLE]

whenever $U_{1}$ and $U_{2}$ are (non-antipodal) members of the unitary group.

Our second result of this subsection generalizes the preceding proposition in the following way. Let

[TABLE]

denote the weighted arithmetic mean [31, Definition 5.1] of unitary matrices $U_{1},\dots,U_{n}\in\mathbb{C}^{m\times m}$, where $w\in\mathbb{R}^{n}$ is a vector of weights summing to 1. Let

[TABLE]

denote their weighted geometric mean [31, Definition 5.2], where

[TABLE]

denotes the geodesic distance on the unitary group [31, Equation (2.6)]. It can be shown [12, Proposition 4] that $\mathbb{A}(U_{1},\dots,U_{n};w)$ exists whenever $\sum_{i=1}^{n}w_{i}U_{i}$ is nonsingular, and is given explicitly by

[TABLE]

On the other hand, $\mathbb{G}(U_{1},\dots,U_{n};w)$ is characterized implicitly by the condition [31, p. 14]

[TABLE]

The following theorem reveals that if the data $U_{1},\dots,U_{n}$ are nearby, then their weighted arithmetic and geometric means are superclose.

Theorem 2.7**.**

Let $U_{i}:[0,T]\rightarrow\mathbb{C}^{m\times m}$, $i=1,2,\dots,n$, be continuous functions on $[0,T]$ such that $U_{i}(t)$ is unitary for each $t\in[0,T]$. Suppose that there exists $C>0$ such that

[TABLE]

for every $i,j=1,2,\dots,n$ and every $t\in[0,T]$. Then, for any $w\in\mathbb{R}^{n}$ with entries summing to 1,

[TABLE]

3 Proofs

In this section, we prove Theorems 2.1, 2.3, 2.5, and 2.7 and Proposition 2.6. Our proofs are structured as follows. In Section 3.1, we consider the general problem of estimating

[TABLE]

where $A\in\mathbb{C}^{m\times p}$ ($m\geq p$) is a full-rank matrix and $\widetilde{U}\in\mathbb{C}^{m\times p}$ has orthonormal columns. We show in Lemma 3.2 that this quantity can be estimated by measuring (1) the extent to which $\widetilde{U}^{*}A$ fails to be Hermitian and (2) the discrepancy between the range of $A$ and the range of $\widetilde{U}$. We then leverage this lemma to prove Theorem 2.1 in Section 3.2, Theorem 2.3 in Section 3.4, Theorem 2.5 in Section 3.5, and Theorem 2.6 and Proposition 2.7 in Section 3.6.

It turns out that one of the theorems proved below, Theorem 2.1, admits an alternative proof that does not rely on Lemma 3.2. This alternative proof, which relies instead on a relationship between projected polynomials and Padé approximation, is shorter than the one we present in Section 3.2, so we detail it in Section 3.3 for completeness. We have chosen to retain both proofs in this paper for several reasons. The proof in Section 3.2, despite being longer, highlights the versatility of Lemma 3.2, a lemma whose wide-ranging applicability is, in our opinion, one of the key contributions of this paper. The proof in Section 3.2 is also more elementary, in a certain sense, than that in Section 3.3, since the former relies merely on well-known perturbation estimates for the polar decomposition, whereas the latter relies on Padé approximation theory and certain results concerning the commutativity of functions of matrices.

3.1 Perturbations of the Polar Decomposition

We begin our examination of (22) by studying the sensitivity of the polar decomposition to perturbations. In what follows, we continue to use $\|\cdot\|$ to denote the Frobenius norm. We denote the $i^{th}$ largest singular value of a matrix $A\in\mathbb{C}^{m\times p}$ by $\sigma_{i}(A)$. If $m=p$, we use $\rho(A)$ to denote the spectral radius of $A$. If furthermore $A$ has real eigenvalues, we denote them by $\lambda_{1}(A)\geq\lambda_{2}(A)\geq\dots\lambda_{p}(A)$. Note that with this convention, it need not be true that $|\lambda_{1}(A)|\geq|\lambda_{2}(A)|\geq\dots|\lambda_{p}(A)|$.

We denote by

[TABLE]

and

[TABLE]

the Hermitian and skew-Hermitian parts of a square matrix $A$, respectively. Note that since

[TABLE]

we have

[TABLE]

and

[TABLE]

for any square matrix $A$.

We will make use of the following additional properties of the Frobenius norm. For any $A\in\mathbb{C}^{m\times p}$, any $B\in\mathbb{C}^{p\times q}$, any $C\in\mathbb{C}^{p\times p}$, and any $U\in\mathbb{C}^{m\times p}$ with orthonormal columns:

1. (LABEL:enumisec:polarpert.i)

$\|A^{*}B\|\leq\|A^{*}\|\sigma_{1}(B)=\|A\|\sigma_{1}(B)$ [20, Equation (B.7)]. 2. (LABEL:enumisec:polarpert.ii)

$\|U^{*}A\|\leq\|A\|$ (This follows from (sec:polarpert.i) and (23)). 3. (LABEL:enumisec:polarpert.iii)

$\|UC\|=\|C\|$ [20, Problem B.7]. 4. (LABEL:enumisec:polarpert.iv)

$\rho(C)\leq\|C\|$ [20, Equation (B.8)].

Note that (sec:polarpert.i) is sharper than the estimate $\|A^{*}B\|\leq\|A^{*}\|\|B\|$, so we will frequently use (sec:polarpert.i) instead of the latter estimate.

We first recall a result concerning the stability of the polar decomposition under perturbations. A proof is given in [28].

Lemma 3.1**.**

[28, Theorem 2.4]** Let $A,\widetilde{A}\in\mathbb{C}^{m\times p}$ ($m\geq p$) be full-rank matrices with polar decompositions $A=UH$ and $\widetilde{A}=\widetilde{U}\widetilde{H}$, where $U,\widetilde{U}\in\mathbb{C}^{m\times p}$ have orthonormal columns and $H,\widetilde{H}\in\mathbb{C}^{p\times p}$ are Hermitian positive-definite. Then

[TABLE]

Next, we consider a full-rank matrix $A\in\mathbb{C}^{m\times p}$ with polar decomposition $A=UH$, and we use the preceding lemma to show that the distance from $U$ to any other matrix $\widetilde{U}\in\mathbb{C}^{m\times p}$ (sufficiently close to $A$) with orthonormal columns is controlled by two properties: (1) the extent to which $\widetilde{U}^{*}A$ fails to be Hermitian, and (2) the discrepancy between the range of $A$ and the range of $\widetilde{U}$.

Lemma 3.2**.**

Let $A\in\mathbb{C}^{m\times p}$ ($m\geq p$) be a full-rank matrix with polar decomposition $A=UH$, where $U=\mathcal{P}(A)\in\mathbb{C}^{m\times p}$ has orthonormal columns and $H\in\mathbb{C}^{p\times p}$ is Hermitian positive-definite. Then, for any matrix $\widetilde{U}\in\mathbb{C}^{m\times p}$ with orthonormal columns satisfying $\|A-\widetilde{U}\|<1$, we have

[TABLE]

Proof 3.3**.**

Define $\widetilde{H}=\mathrm{sym}(\widetilde{U}^{*}A)$. This matrix is positive-definite, since it is a small-normed perturbation of the identity matrix. Indeed, the relation

[TABLE]

implies $\|\widetilde{H}-I\|\leq\|\widetilde{U}^{*}(A-\widetilde{U})\|$. Thus, using (sec:polarpert.ii) and (sec:polarpert.iv), the smallest eigenvalue of $\widetilde{H}$ satisfies

[TABLE]

Now define $\widetilde{A}=\widetilde{U}\widetilde{H}$. Observe that

[TABLE]

so

[TABLE]

by (sec:polarpert.iii). The right-hand inequality in (26) then follows from Lemma 3.1 upon noting that $\widetilde{A}$ and $\mathrm{sym}(\widetilde{U}^{*}A)$ have the same singular values.

To prove the left-hand inequality in (26), observe that since $H=U^{*}A$ is Hermitian,

[TABLE]

Thus, using (sec:polarpert.i) and (25),

[TABLE]

On the other hand, since $UU^{*}A=UH=A$, we have

[TABLE]

Thus, using (sec:polarpert.i), (sec:polarpert.iii), and the fact that $\sigma_{1}(U^{*}A)=\sigma_{1}(A)$, it follows that

[TABLE]

*Combining (27) and (28) proves the left-hand inequality in (26). *

The following less sharp version of Lemma 3.2, applicable in the square case ($m=p$), will be useful in the upcoming sections.

Lemma 3.4**.**

Let $A$, $U$, and $\widetilde{U}$ be as in Lemma 3.2. If $A$ is square (i.e. $m=p$), then

[TABLE]

Proof 3.5**.**

*Use (26) together with the fact that $\sigma_{1}(A)\leq\|A\|$, $\sigma_{p}(A)^{-1}\leq\|A^{-1}\|$, and $\widetilde{U}\widetilde{U}^{*}=I$ when $\widetilde{U}$ is square. *

3.2 Exponentiation on the Unitary Group

We now prove Theorem 2.1. Fix an integer $n\geq 0$ and consider a polynomial of the form

[TABLE]

with $a_{0}=1$ and $a_{k}\in\mathbb{C}$, $k=1,2,\dots,n$. We aim to find coefficients $a_{k}$ making $\mathcal{P}(q_{n}(t\Omega))-e^{t\Omega}$ small for any skew-Hermitian matrix $\Omega$. Applying Lemma 3.4 with $A=q_{n}(t\Omega)$ and $\widetilde{U}=e^{t\Omega}$ gives

[TABLE]

provided that $t$ is sufficiently small (small enough that $q_{n}(t\Omega)$ has full rank and $\|q_{n}(t\Omega)-e^{t\Omega}\|<1$).

This inequality is of great utility, since $\|q_{n}(t\Omega)\|^{-1}$ and $\|q_{n}(t\Omega)^{-1}\|$ are each $O(1)$ as $t\rightarrow 0$, and

[TABLE]

can be expanded in powers of $t$. Namely,

[TABLE]

where

[TABLE]

The quantity $e^{-t\Omega}q_{n}(t\Omega)-q_{n}(-t\Omega)e^{t\Omega}$ is thus of the highest order in $t$ when the $n$ coefficients $a_{k}$, $k=1,2,\dots,n$, are chosen to make $b_{l}=0$ for $l=1,3,5,\dots,2n-1$. This is achieved when

[TABLE]

as the following lemma shows.

Lemma 3.6**.**

With $a_{0}=1$ and $a_{k}$ given by (30) for $1\leq k\leq n$, we have

[TABLE]

Proof 3.7**.**

Substitution gives

[TABLE]

so it suffices to show that

[TABLE]

for each $l=1,3,5,\dots,2n-1$. To prove this, consider the polynomial

[TABLE]

The coefficient of $z^{n}$ in this polynomial is precisely

[TABLE]

Indeed,

[TABLE]

and taking $j=n$ in the inner summation above gives (32). Now observe that $r(z)$ satisfies the symmetry

[TABLE]

*From this it follows that the coefficients $r_{k}$ satisfy $r_{k}=(-1)^{l}r_{2n-k}$, $k=0,1,\dots,2n$. In particular, $r_{n}=(-1)^{l}r_{n}$, so $r_{n}=0$ when $l$ is odd. *

This completes the proof of Theorem 2.1, since by the inequality (29), the unitary factor in the polar decomposition of the polynomial $\sum_{k=0}^{n}a_{k}t^{k}\Omega^{k}$ with coefficients given by (30) delivers an approximation of $e^{t\Omega}$ with error of order $t^{2n+1}$.

Uniqueness of the solution (30) to (31) is a consequence of the following lemma, which proves that the linear system (31) is nonsingular.

Lemma 3.8**.**

*Fixing $a_{0}=1$, the linear system (31) in the $n$ unknowns $a_{1},a_{2},\dots,a_{n}$ is nonsingular. *

Proof 3.9**.**

Upon rearrangement, (31) reads

[TABLE]

where $x,y\in\mathbb{C}^{n}$ and $M\in\mathbb{C}^{n\times n}$ have entries given by

[TABLE]

An inductive argument shows that the determinant of $M$ is equal to

[TABLE]

*where $l!!=\prod_{j=0}^{\lceil l/2\rceil-1}(l-2j)$ denotes the double factorial. In particular, $\det M\neq 0$, showing that $M$ is nonsingular. *

3.3 Connections with Padé Approximation

We now present an alternative proof of Theorem 2.1 that relies not on Lemma 3.2, but rather on a connection between $\mathcal{P}(\Theta_{n}(t\Omega))$ and the diagonal Padé approximant of $e^{2t\Omega}$.

Our alternative proof will make use of the fact that if $A\in\mathbb{C}^{m\times p}$ ($m\geq p$) has full rank, then

[TABLE]

where $C^{-1/2}$ denotes inverse of the principal square root of a square matrix $C$ with no nonpositive real eigenvalues [20, Theorem 8.1].

It will also make use of the following facts: If $f$ and $g$ are two scalar-valued functions defined on the spectrum of a square matrix $A$, then $f(A)$ and $g(A)$ are well-defined [20, Section 1.2], $f(A)$ commutes with $g(A)$ [20, Theorem 1.13(e)], and the spectrum of $f(A)$ is the image of the spectrum of $A$ under $f$ [20, Theorem 1.13(d)].

Lemma 3.10**.**

Let $q_{n}(z)$ be a polynomial of degree $n\geq 0$ with $q_{n}(0)\neq 0$, and let $\Omega\in\mathbb{C}^{m\times m}$ be skew-Hermitian. For each $t$ sufficiently small (small enough so that $q_{n}(t\Omega)$ is nonsingular), we have

[TABLE]

*Furthermore, $\mathcal{P}(q_{n}(t\Omega))$ commutes with $e^{t\Omega}$. *

Proof 3.11**.**

Since $\Omega$ is skew-Hermitian, $q_{n}(t\Omega)^{*}=q_{n}(-t\Omega)$. Hence, by (33),

[TABLE]

This shows that $\mathcal{P}(q_{n}(t\Omega))=f(t\Omega)$, where $f(z)=q_{n}(z)(q_{n}(-z)q_{n}(z))^{-1/2}$. If $q_{n}(t\Omega)$ is nonsingular, then $f$ is defined on the spectrum of $t\Omega$. Indeed, if $\lambda$ is an eigenvalue of $t\Omega$, then $q_{n}(\lambda)$, being an eigenvalue of $q_{n}(t\Omega)$, is nonzero, and $q_{n}(-\lambda)$, being an eigenvalue of $q_{n}(-t\Omega)=q_{n}(t\Omega)^{*}$, is nonzero. Thus, $\mathcal{P}(q_{n}(t\Omega))$ commutes with any function of $t\Omega$ defined on the spectrum of $t\Omega$, including $e^{t\Omega}$. By similar reasoning, $\left[q_{n}(-t\Omega)q_{n}(t\Omega)\right]^{-1/2}$ commutes with $q_{n}(t\Omega)$, so

[TABLE]

The preceding lemma implies that if $t$ is sufficiently small, then

[TABLE]

Using this identity, it is not hard to see that the polynomial $q_{n}$ for which $\mathcal{P}(q_{n}(t\Omega))-e^{t\Omega}$ is of the highest order in $t$ is precisely that for which $q_{n}(t\Omega)q_{n}(-t\Omega)^{-1}-e^{2t\Omega}$ is of the highest order in $t$. That polynomial is none other than the numerator in the diagonal Padé approximant of $e^{2t\Omega}$, which is precisely $\Theta_{n}(t\Omega)$ [23, p. 97].

3.4 Exponentiation on the Grassmannian

We now prove Theorem 2.3. The proof will consist of two parts. First, we prove the identity (4) by exploiting the block structure of the matrix $Z=\begin{pmatrix}0&-K^{*}\\ K&0\end{pmatrix}$. Then, we insert the right-hand side of (4) into (3) and expand the result to obtain Theorem 2.3.

Throughout this subsection, we make use of the identities

[TABLE]

and

[TABLE]

which hold for for every nonnegative integer $j$.

Lemma 3.12**.**

Let $r(z)=c_{0}+c_{1}z+c_{2}z^{2}+\dots+c_{n}z^{n}$ be a polynomial, let $K\in\mathbb{C}^{(m-p)\times p}$, and let $Z=\begin{pmatrix}0&-K^{*}\\ K&0\end{pmatrix}$. Then

[TABLE]

*where $B=r(Z)\begin{pmatrix}I\\ 0\end{pmatrix}$ and $C=r(Z)\begin{pmatrix}0\\ I\end{pmatrix}$. *

Proof 3.13**.**

*The diagonal blocks of $r(Z)^{*}r(Z)$ are automatically given by $B^{*}B$ and $C^{*}C$, so it suffices to show that the off-diagonal blocks of $r(Z)^{*}r(Z)$ vanish. To this end, observe that the skew-Hermiticity of $Z$ implies $r(Z)^{*}r(Z)=r(-Z)r(Z)$. But $r(Z)^{*}r(Z)$ is Hermitian, so taking the Hermitian part of both sides gives $r(Z)^{*}r(Z)=\mathrm{sym}\left(r(-Z)r(Z)\right)$. Since $\mathrm{sym}(Z^{j})=0$ for odd $j$, it follows that $r(Z)^{*}r(Z)$ is a linear combination of even powers of $Z$, all of which are block diagonal by (34). *

Lemma 3.14**.**

Let $r(z)=c_{0}+c_{1}z+c_{2}z^{2}+\dots+c_{n}z^{n}$ be a polynomial, let $K\in\mathbb{C}^{(m-p)\times p}$, and define $Z=\begin{pmatrix}0&-K^{*}\\ K&0\end{pmatrix}$. If $r(Z)$ has full rank, then

[TABLE]

Proof 3.15**.**

In the notation of Lemma 3.12,

[TABLE]

so

[TABLE]

On the other hand,

[TABLE]

*as well. *

The preceding lemma establishes the identity (4). We now study the quantity $\Theta_{n}(tZ)\begin{pmatrix}I\\ 0\end{pmatrix}$ in more detail.

Lemma 3.16**.**

Let $\alpha_{n}(z)$ and $\beta_{n}(z)$ be as in Theorem 2.3, let $K\in\mathbb{C}^{(m-p)\times p}$, and let $Z=\begin{pmatrix}0&-K^{*}\\ K&0\end{pmatrix}$. Then

[TABLE]

Proof 3.17**.**

Using (34-35), we have

[TABLE]

We are now in a position to prove Theorem 2.3 by substituting the preceding results into (3). Combining Lemmas 3.14 and 3.16, we have

[TABLE]

where the second line follows from (10), and the last line follows from the fact that $H=Y_{\perp}K$, and $Y_{\perp}$ has orthonormal columns. This, together with (3), completes the proof of Theorem 2.3.

3.5 Exponentiation on the Stiefel Manifold

We now turn to the proof of Theorem 2.5. Let $q(x,y)$ and $r(x,y)$ be polynomials in non-commutating variables $x$ and $y$, and define

[TABLE]

For the moment we assume only that $q(0,0)=1$, but later we will make the additional assumptions (sec:stiefelexp.i-sec:stiefelexp.iii) (the first of which implies $q(0,0)=1$). Using the identities $Y=\begin{pmatrix}Y&Y_{\perp}\end{pmatrix}\begin{pmatrix}I\\ 0\end{pmatrix}$, $H=Y\Omega+Y_{\perp}K=\begin{pmatrix}Y&Y_{\perp}\end{pmatrix}\begin{pmatrix}\Omega\\ K\end{pmatrix}$, $H^{*}H=K^{*}K-\Omega^{2}$, and $Y^{*}H=\Omega$, we can write

[TABLE]

where we have suppressed the arguments to $q$ and $r$ in the last line to reduce clutter.

Now let $Z=\begin{pmatrix}\Omega&-K^{*}\\ K&0\end{pmatrix}$ and define

[TABLE]

We aim to bound

[TABLE]

using Lemma 3.2. Since A\big{|}_{t=0}=\begin{pmatrix}Y&Y_{\perp}\end{pmatrix} is unitary, it follows that $\sigma_{i}(A)=O(1)$ as $t\rightarrow 0$ for each $i=1,2,\dots,p$. Thus, it is enough to bound $\|\mathrm{skew}(\widetilde{U}^{*}A)\|$ and $\|(I-\widetilde{U}\widetilde{U}^{*})A\|$. We begin with a lemma.

Lemma 3.18**.**

We have

[TABLE]

Proof 3.19**.**

The first equality follows from a direct calculation, using the fact that $Z$ is skew-Hermitian and $\begin{pmatrix}Y&Y_{\perp}\end{pmatrix}$ is unitary. For the second, observe that

[TABLE]

*The result follows from the fact that the Frobenius norm is unitarily invariant, and $\begin{pmatrix}Y&Y_{\perp}\end{pmatrix}$ and $e^{tZ}$ are unitary. *

The preceding lemma reveals that the order of accuracy of the approximation (12) can be determined by studying the quantity

[TABLE]

Let us carry out this task in order to determine, as an illustration, the highest order approximation of the form (12) that can be achieved using polynomials $q(x,y)$ and $r(x,y)$ satisfying (sec:stiefelexp.i-sec:stiefelexp.iii) with $n=2$. The cases $n=1$ and $n=3$ are handled similarly; we leave those details to the reader. Collectively, these arguments will prove Theorem 2.5.

It is a simple exercise to show that when $n=2$, the only polynomials $q(x,y)$ and $r(x,y)$ satisfying (sec:stiefelexp.i-sec:stiefelexp.ii) are of the form

[TABLE]

where $c$ is a constant. Substituting into (36), writing $e^{-tZ}=\sum_{k=0}^{\infty}\frac{(-tZ)^{k}}{k!}$, and multiplying, one finds after a tedious calculation that

[TABLE]

Hence, by Lemma 3.18 and the symmetry of $K^{*}K$ and $\Omega^{2}$, we have

[TABLE]

The optimal choice of $c$ is $c=-\frac{1}{2}$, giving

[TABLE]

and

[TABLE]

It follows that

[TABLE]

Clearly, no other choice of $c$ will improve this approximant’s order of accuracy, proving (15) for $n=2$.

If it happens that $Y^{*}H=\Omega=0$, then (2) and (11) coincide, and (sec:stiefelexp.i-sec:stiefelexp.ii) and Theorem 2.3 imply

[TABLE]

Likewise, if $m=p$, so that $H=Y\Omega$, $YY^{*}=I$, and $\mathrm{Exp}^{St}_{Y}(tH)=Ye^{t\Omega}$, then (sec:stiefelexp.iii), (10), and Theorem 2.1 imply

[TABLE]

These observations prove (13-14) for the case $n=2$. The proof of Theorem 2.5 is completed by performing analogous arguments for the cases $n=1$ and $n=3$.

3.6 Geometric and Arithmetic Means of Unitary Matrices

We now prove Proposition 2.6 and Theorem 2.7.

Proof of Proposition 2.6

Without loss of generality, consider the case in which $U_{1}=I$, so that

[TABLE]

By Lemma 3.4, it suffices to examine the norm of

[TABLE]

The series expansion of $e^{-st\Omega}\left((1-s)I+se^{t\Omega}\right)$ reads

[TABLE]

Since $\Omega$ is skew-Hermitian, $\mathrm{skew}(\Omega^{k})=0$ for every even $k$, showing that

[TABLE]

When $s=1/2$, each term in the series vanishes, giving

[TABLE]

When $s\neq 1/2$, the first non-vanishing term is of order $t^{3}$, showing that

[TABLE]

The result follows by applying Lemma 3.4.

Proof of Theorem 2.7

Let

[TABLE]

and

[TABLE]

Observe that if $A(t)=U(t)H(t)$ is the polar decomposition of $A(t)$, then, by (19),

[TABLE]

Moreover, using the fact that $\sum_{i=1}^{n}w_{i}=1$, we have

[TABLE]

This shows, by (21), that

[TABLE]

The latter matrix is unitary, so $\|A(t)\|=O(1)$. This is independent of the $U_{i}(t)$ we chose to pull out of the sum, since (21) implies that $U_{1}(0)=U_{2}(0)=\dots=U_{n}(0)$.

Now let $\Omega_{i}(t)=\frac{1}{t}\log(\widetilde{U}(t)^{*}U_{i}(t))$ for each $i$, so that

[TABLE]

Note that $\Omega_{i}(t)=O(1)$ by (21). In addition, by (20),

[TABLE]

Suppressing the dependencies on $t$ for ease of reading, it follows that

[TABLE]

The skew-Hermitian part of $\widetilde{U}^{*}A$ is thus given by

[TABLE]

Since

[TABLE]

and $\|A\|=O(1)$, it follows from Lemma 3.4 that

[TABLE]

i.e.,

[TABLE]

4 Numerical Examples

In this section, we discuss how the projected polynomials proposed in Theorems 2.1, 2.3, and 2.5 can be efficiently computed, focusing on iterative methods for computing the polar decomposition. We then present numerical examples that illustrate their order of accuracy.

4.1 Iterative Methods for Computing the Polar Decomposition

The cost of computing a projected polynomial is largely dominated by the cost of evaluating the map $\mathcal{P}$. This map can be computed efficiently via a number of different iterative methods. The most widely known, applicable when $m=p$, is the Newton iteration

[TABLE]

The iterates $X_{k}$ so defined converge quadratically to the unitary factor $\mathcal{P}(A)=U$ in the polar decomposition $A=UH$ for any nonsingular square matrix $A$ [20, Theorem 8.12]. A closely related iteration, applicable when $m\geq p$, is given by

[TABLE]

These iterates converge quadratically to $\mathcal{P}(A)$ for any full-rank $A\in\mathbb{C}^{m\times p}$ ($m\geq p$) [20, Corollary 8.14(b)]. Finally, the Newton-Schulz iteration

[TABLE]

provides an inverse-free iteration whose iterates $X_{k}$ converge quadratically to $\mathcal{P}(A)$ for any $A\in\mathbb{C}^{m\times p}$ ($m\geq p$) whose singular values all lie in the interval $(0,\sqrt{3})$ [20, Problem 8.20]. For further information, including other iterations for computing $\mathcal{P}(A)$, see [20, Chapter 8].

4.2 Numerical Convergence

We tested the accuracy of the projected polynomials detailed in Theorems 2.1, 2.3, and 2.5 by applying them to randomly generated inputs. To calculate $\mathcal{P}$, we used (37) for square matrices and (39) for rectangular matrices. The results of the tests, detailed in Tables 1-3, corroborate the convergence rates predicted by the theory.

5 Conclusion

This paper has presented a family of high-order retractions on the unitary group, the Grassmannian manifold, and the Stiefel manifold. All of these retractions were constructed by projecting certain matrix polynomials onto the set of matrices with orthonormal columns using the polar decomposition, or, in the case of the Grassmannian, using either the polar decomposition or the QR decomposition. There are several interesting applications and extensions of this strategy that seem worthwhile to pursue. On quadratic Lie groups other than the unitary group, one might consider adopting the same strategy, replacing the polar decomposition with the generalized polar decomposition [32, 22]. It might also be worthwhile to consider projecting rational functions, rather than polynomials, to achieve higher accuracy for comparable cost. It may also be possible to leverage these retractions, together with methods for computing their derivatives [14], to construct high-order approximations of parallel transport operators on matrix manifolds; see [2, Section 8.1.2].

It is worth noting that many of the constructions in this paper might generalize nicely to infinite dimensions. For instance, replacing the matrices $Y$ and $H$ in Theorem 2.3 with quasi-matrices in the sense of [35], one obtains a method for approximating geodesics between finite-dimensional function spaces, with $H^{*}H$ playing the role of a Gramian, and with $\mathcal{P}$ interpreted as the map sending an ordered basis of functions to the nearest ordered, orthonormal basis of functions.

6 Acknowledgements

We wish to thank the developers of the non-commutative algebra package NCAlgebra [19], which we used to carry out some of the calculations that appeared/were mentioned in Section 3.5. The first author is supported by NSF grants CMMI-1334759, DMS-1345013, and DMS-1703719. The second author is supported by NSF grants CMMI-1334759, DMS-1345013, DMS-1411792.