The generalized orthogonal Procrustes problem in the high noise regime

TL;DR

This paper addresses the challenge of estimating point clouds from noisy observations with unknown rotations in high noise regimes, proposing an invariant polynomial-based method that is statistically optimal and analyzing its variance.

Contribution

It introduces a simple, efficient invariant polynomial approach for the generalized orthogonal Procrustes problem, effective across different noise levels, and provides new bounds on Cholesky factorization sensitivity.

Findings

Method is statistically optimal in high noise regimes.

Proposed approach adapts to noise levels for accurate recovery.

Improved bounds on Cholesky factorization sensitivity are established.

Abstract

We consider the problem of estimating a cloud of points from numerous noisy observations of that cloud after unknown rotations, and possibly reflections. This is an instance of the general problem of estimation under group action, originally inspired by applications in 3-D imaging and computer vision. We focus on a regime where the noise level is larger than the magnitude of the signal, so much so that the rotations cannot be estimated reliably. We propose a simple and efficient procedure based on invariant polynomials (effectively: the Gram matrices) to recover the signal, and we assess it against fundamental limits of the problem that we derive. We show our approach adapts to the noise level and is statistically optimal (up to constants) for both the low and high noise regimes. In studying the variance of our estimator, we encounter the question of the sensivity of a type of thin Cholesky factorization, for which we provide an improved bound which may be of independent interest.

Equations

$$Y_i=Q_{i}X+\sigma E_i,$$

$$\mathcal{O}(d)$$

$$X_1\sim X_2⟺X_1=Q{X}_2\text{ for some }Q\in \mathcal{O}(d).$$

$$[X]$$

$$\underset{\hat{X}\in {\mathbb{R}}^{d\times k},{\hat{Q}}_1,\dots ,{\hat{Q}}_N\in \mathcal{O}(d)}{\min }\sum _{i=1}^N\Vert {\hat{Q}}_i^T{Y}_i-\hat{X}{\Vert }_{\mathrm{F}}^2.$$

$$\hat{X}=\frac{1}{N}\sum _{i=1}^N{\hat{Q}}_i^T{Y}_i.$$

$$\rho ([X_1],[X_2])=\underset{Q\in \mathcal{O}(d)}{\min }\Vert X_1-Q{X}_2{\Vert }_{\mathrm{F}},$$

$$\underset{{\hat{Q}}_1,\dots ,{\hat{Q}}_N\in \mathcal{O}(d)}{\min }\sum _{i\ne j}\Vert {\hat{Q}}_i{Y}_i-{\hat{Q}}_j{Y}_j{\Vert }_{\mathrm{F}}^2,$$

$$\underset{{\hat{Q}}_1,\dots ,{\hat{Q}}_N\in \mathcal{O}(d)}{\max }\sum _{i\ne j}⟨{\hat{Q}}_i{Y}_i,{\hat{Q}}_j{Y}_j⟩,$$

$$p(QX)$$

$$p(X)$$

$$p_0+p_1(X)+\cdots +p_r(X)=p(X)=p(QX)=p_0+p_1(QX)+\cdots +p_r(QX).$$

$$p(X)=\sum _{i,j}{m}_{ij}⟨x_i,x_j⟩=⟨M,X^{T}X⟩.$$

$$p(X)$$

$$p(QX)$$

$$=⟨\mathrm{vec}(X)\mathrm{vec}(X{)}^T,(I_k\otimes Q{)}^{T}A(I_k\otimes Q)⟩.$$

$$A$$

$$A_{ij}$$

$$p(X)$$

$$=⟨\mathrm{vec}(X),\mathrm{vec}(X{M}^T)⟩$$

$$=⟨X,X{M}^T⟩$$

$$=⟨X^{T}X,M⟩.$$

$$Y_i^T{Y}_i=X^{T}X+\sigma \left(X^T{Q}_i^T{E}_i+E_i^T{Q}_{i}X\right)+{\sigma }^2{E}_i^T{E}_i.$$

$${\hat{M}}_N$$

$${\hat{G}}_N$$

$$\underset{\hat{X}\in {\mathbb{R}}^{d\times k}}{\min }\Vert {\hat{G}}_N-{\hat{X}}^T\hat{X}{\Vert }_{\mathrm{F}}.$$

$$d{\hat{\sigma }}_N^2$$

$$\mathcal{M}={\mathbb{R}}_{\ast }^{d\times k}/\sim ,$$

$$\mathcal{N}=\{G\in {\mathrm{Sym}}_k:\mathrm{rank}(G)=d,G⪰0\},$$

$$\rho ([X],[\tilde{X}])$$

$$\displaystyle\begin{array}[]{ccccc}\mathcal{L}_{X}&:&\mathbb{R}^{d\times k}&\to&\mathrm{Sym}_{k}\\ &&\dot{X}&\mapsto&X^{T}\dot{X}+\dot{X}^{T}X.\end{array}$$

$${\mathrm{V}}_X$$

Full text

2000 Math Subject Classification: 34K30, 35K57, 35Q80, 92D25

1 Introduction

We consider the problem of estimating $k$ labeled points in $\mathbb{R}^{d}$, with $k\geq d$. This cloud of points, which we call the parameter, is represented as a matrix $X$ of size $d\times k$. We restrict ourselves to the case where the smallest singular value of $X$ is bounded away from zero, that is, the cloud spans all $d$ dimensions. We observe $N$ independent measurements $Y_{1},\ldots,Y_{N}$ of $X$, following the model

[TABLE]

where $\sigma>0$ is the standard deviation of the noise, $E_{1},\ldots,E_{N}$ are independent noise matrices in $\mathbb{R}^{d\times k}$ with independent, standard Gaussian entries and $Q_{1},\ldots,Q_{N}$ are drawn uniformly and independently at random from the orthogonal group,

[TABLE]

where $I_{d}$ is the identity matrix of size $d$. In what follows, we refer to orthogonal matrices as rotations, bearing in mind that (in our meaning) they may also include a reflection. The method we propose in Section 2.2 applies under relaxed assumptions on the distributions of $Q_{i}$ and $E_{i}$: our assumptions here serve to streamline exposition. This problem we investigate belongs to a larger class of estimation problems under group actions (Bandeira et al., 2017b).

Notice that the distribution of the observations $Y_{i}$ is unchanged if $X$ is replaced by $QX$, for any orthogonal $Q$. As a result, we can only hope to recover the cloud $X$ up to a global rotation. Accordingly, we define an equivalence relation $\sim$ over $\mathbb{R}^{d\times k}$:

[TABLE]

This equivalence relation partitions the parameter space into equivalence classes

[TABLE]

The set of equivalence classes is the quotient space $\mathbb{R}^{d\times k}/\!\sim$. The distribution of the measurements $Y_{i}$ is parameterized by $[X]$, which we aim to estimate.

A natural approach to estimate $[X]$ would be to estimate the rotations $Q_{i}$ first (seen as latent or nuisance variables), align the observations $Y_{i}$ using the estimated rotations, and average. Typical of those approaches, the maximum likelihood estimator (MLE) is a solution of the following non-convex optimization problem:

[TABLE]

For any fixed choice of estimators $\hat{Q}_{1},\ldots,\hat{Q}_{N}$, the optimal estimator for $\hat{X}$ according to the above is

[TABLE]

This estimator can be plugged into the cost function of (4), reducing the problem to that of estimating only the rotations, upon which $\hat{X}$ can be deduced from (5). A number of papers focus on the resulting problem, called synchronization of rotations (Singer, 2011).

Such approaches, however, necessarily fail at low signal-to-noise ratio (SNR). Specifically, we argue that if the noise level is too large, no procedure can reliably determine the latent rotations $Q_{i}$. Essentially, this is because, when $\sigma$ is too large, the distribution of $QX+\sigma E$ is indistinguishable from that of $Q^{\prime}X+\sigma E$, where $Q$ and $Q^{\prime}$ are two rotations.

To see this, consider the following strictly simpler problem: we observe $Y=sQX+\sigma E$, where $Q\in\mathcal{O}(d)$, $X\in\mathbb{R}^{d\times k}$ and $\sigma>0$ are known, while $E\in\mathbb{R}^{d\times k}$ has independent standard Gaussian entries and $s$ is uniformly sampled from $\{+1,-1\}$, both unknown. An estimator $\psi$ (deterministic) assigns an estimate of $s$ to an observed $Y$. Estimating $s$ in this context is strictly simpler than estimating $Q_{i}$’s to any reasonable accuracy in our model (where furthermore $X$ is unknown), as we are only asked to determine whether $Q_{i}$ is some given rotation, or its opposite. Yet, even this simpler problem is hopeless at low SNR, directly implying the impossibility of estimating the rotations in problem (4):

Proposition 1.1**.**

*For any tolerance $\tau\in(0,1/2)$, there exists a critical noise level $\sigma_{0}$ such that, for any estimator $\psi$, if $\sigma>\sigma_{0}$, then the probability of error $\mathbb{P}(\psi(Y)\neq s)$ exceeds $\tau$. *

(The proof relies on the optimality of the likelihood ratio test for Gaussian distributions, see Appendix B.)

Because of this fundamental obstruction, in this paper, we aim to estimate $[X]$ directly from observations $Y_{i}$, bypassing any estimation (even implicit) of the latent $Q_{i}$’s. To do so, we follow a trend in signal processing that consists in estimating $[X]$ from features of the observations that are invariant under the group action of $\mathcal{O}(d)$  (Tukey, 1984), (Sadler and Giannakis, 1992), (Giannakis, 1989), (Abbe et al., 2018a), (Bandeira et al., 2017b), (Perry et al., 2017), (Boumal et al., 2018).

Specifically, consider the Gram matrix of $X$, that is, $X^{T}X$: it is invariant under orthogonal transformations since $(QX)^{T}(QX)=X^{T}X$ for any $Q\in\mathcal{O}(d)$. Thus, up to noise terms that we will handle, the Gram matrices of the observations $Y_{i}$ reveal information about the Gram matrix of $X$ without the need to estimate the $Q_{i}$’s. In Section 2, we call upon invariant theory to argue that no other polynomial invariant features are necessary, in the sense that (a) they would be redundant with the Gram matrix, and (b) the Gram matrix is sufficient to fully characterize the equivalence class $[X]$.

Based on these observations, we proceed (still in Section 2) to derive an estimator for the Gram matrix of $X$ from the given observations (1)—this also requires estimating the noise level $\sigma$, which we discuss. From the estimated Gram matrix, we construct an estimator for the equivalence class $[X]$. Since we only have access to a finite number $N$ of observations, we can only hope to recover an approximation of the Gram matrix. Accordingly, in Section 3 we study the sensitivity of the mapping from the Gram matrix to the sought equivalence class $[X]$. This reduces to showing stability of the factorization of a rank-$d$ positive semidefinite matrix. To address this question, we propose a new analysis of such matrix factorization, with a geometric proof.

In Section 4, we show that the proposed approach is statistically optimal, that is, it makes the best use of available samples. Moreover, we show that the mean squared error (MSE) of our estimator behaves like $O\left(\sigma^{2}/N+\sigma^{4}/N\right)$, which is shown to be adaptively optimal. Indeed, this highlights the existence of two regimes: at low SNR, our estimator’s MSE behaves like $O(\sigma^{4}/N)$, while at high SNR we obtain an MSE of order $O(\sigma^{2}/N)$ which matches the regime we would get if the rotations were known. In particular, the MSE can be driven to zero at any noise level, provided the number of observations $N$ is sufficiently large. We give an explicit characterization of those two regimes and we support our claim with numerical simulations in Section 5.

To make sense of MSE in estimating $[X]$, we need a distance on the quotient space. A natural choice consists in computing the Frobenius distance between two aligned representatives, that is,

[TABLE]

where $\lVert X\rVert_{\mathrm{F}}=\sqrt{\mathrm{Tr}(X^{T}X)}$. We sometimes write simply $\rho(X_{1},X_{2})$, where it is clear that we mean the distance between $[X_{1}]$ and $[X_{2}]$. Computing this distance is known as the orthogonal Procrustes problem: it can be done efficiently via singular value decomposition (SVD) (Schönemann, 1966) (see also Appendix A).

Related work

Procrustes problems consist in finding correspondences between shapes that have been transformed through translation, rotation or dilation (Schönemann, 1966; Ten Berge, 1977). They notably find application in multivariate analysis (Hurley and Cattell, 1962; Green and Carroll, 1976), multidimensional scaling (Borg and Groenen, 2005), computer vision (Zhang, 2000) and natural language processing (Xing et al., 2015; Smith et al., 2017; Grave et al., 2019). An extensive survey on applications of the Procrustes problem can be found in (Gower and Dijksterhuis, 2005).

A particular line of work on such problems has focused on estimating each rotation before estimating the orbit (that is, equivalence class) of the matrix $X$. In particular, the MLE (4) is cast into the following non-convex optimization problem:

[TABLE]

or equivalently:

[TABLE]

where $\langle A,B\rangle=\mathrm{Tr}(A^{T}B)$ is the inner product we use throughout. Among many, this problem was investigated by Nemirovski (Nemirovski, 2007) and Man-Cho So (Man-Cho So, 2010). A particularization of this problem amounts to the little Grothendieck problem (Khot and Naor, 2012; Naor et al., 2013). A natural way to deal with the non-convexity and ensuing computational complexity is to study a convex relaxation of (8). Following the line initiated by the seminal work of Goemans & Williamson (Goemans and Williamson, 1995), Bandeira et al. (Bandeira et al., 2016) and many others investigated such semidefinite relaxations. In contrast, we focus on a noise regime where $Q_{i}$ s cannot be estimated; hence these approaches cannot succeed.

The problem we consider finds its original motivation in cryo-electron microscopy (cryo-EM), an imaging technique used in structural biology to estimate the three-dimensional shape of a molecule from two-dimensional projections of that molecule under random and unknown orientations—see (Singer, 2018) for a review of mathematical aspects of this task. Our setting and approach are also closely connected to the multi-reference alignment problem and recent literature on the topic, where noisy realizations of a randomly, cyclically shifted version of a vector are observed (Bandeira et al., 2014).

Close to the approach taken in the present work, Giannakis (Giannakis, 1989), Sadler and Giannakis Sadler and Giannakis (1992) and more recently Bendory et al. (Bendory et al., 2018) and Abbe et al. (Abbe et al., 2018a) among others, consider estimating moments of the signal that are invariant under cyclic shifts. The advantage of such formulation is its validity for any SNR regime, provided sufficiently many samples are available. In particular, this approach has been shown to be optimal both in terms of rate of estimation (Bandeira et al., 2017b) and sample complexity (Perry et al., 2017).

2 Invariant features approach

We focus on estimating $[X]$ from samples drawn according to (1), disregarding the variables $Q_{1},\ldots,Q_{N}$. Based on the discussion above, it is apparent that any method for estimating $[X]$ which (explicitly or implicitly) relies on estimating the rotations reliably must fail beyond a certain noise level. Hence, we take a different approach. For each observation $Y_{i}$, we compute polynomial functions of $Y_{i}$ that, aside from the noise $E_{i}$, are invariant under the rotation $Q_{i}$. Such functions are called invariant features. Then, our estimation problem reduces to that of estimating $[X]$ from the estimated invariant features. This approach completely bypasses estimation of the latent variables. We first explain the method with some background below, before showing that it does not break down at high noise levels.

2.1 Invariant polynomials

Invariant theory is concerned with polynomials that are invariant under some group action. In our case, this specializes to the following central definition.

Definition 2.1**.**

A multivariate polynomial $p\colon\mathbb{R}^{d\times k}\to\mathbb{R}$ is said to be invariant under the action of the orthogonal group if, for all $X\in\mathbb{R}^{d\times k}$ and for all $Q\in\mathcal{O}(d)$,

[TABLE]

The goal of this section is to identify all invariant polynomials in our specific setting. A standard observation in invariant theory is that it is sufficient to consider homogeneous invariant polynomials.

Proposition 2.2**.**

Consider a multivariate polynomial $p$ of degree $r$, decomposed into a sum of homogeneous parts $p_{1},\ldots,p_{r}$ where $p_{i}$ has degree $i$:

[TABLE]

If $p$ is invariant under the action of the orthogonal group, that is, if for all $Q\in\mathcal{O}(d)$, $p(QX)\equiv p(X)$, then each homogeneous part $p_{i}$ is itself invariant under that action.

Proof 2.3**.**

By invariance of $p$, we obtain for all $X$ and orthogonal $Q$ that

[TABLE]

By identifying the homogeneous parts of the polynomials, we obtain for all $i$ that $p_{i}(X)\equiv p_{i}(QX)$ for all $Q\in\mathcal{O}(d)$, hence each $p_{i}$ is itself invariant.

Homogeneous invariant polynomials necessarily have even degree. In particular, there are no interesting invariant polynomials of degree zero or one.

Proposition 2.4**.**

No homogeneous polynomial of odd degree is invariant under orthogonal group action.

Proof 2.5**.**

*By contradiction, let $p(X)$ be a homogeneous polynomial of odd degree, invariant under orthogonal group action. Then, $p(X)=p(QX)$ for all orthogonal $Q$. This holds in particular for $Q=-I_{d}$, so that $p(X)=p(-X)=-p(X)$ (we used that the degree is odd in the last equality). Thus, $p(X)=0$ for all $X$, which contradicts the fact that $p$ has odd degree. *

We now give an elementary statement of a key property of degree-two invariant polynomials: they can be expressed as a linear function of the Gram matrix of $X$, namely $X^{T}X$. We use the notation $[k]=\{1,\ldots,k\}$ and $\langle A,B\rangle=\mathrm{Tr}(A^{T}B)$.

Proposition 2.6**.**

Any homogeneous polynomial $p(X)$ of degree two that is invariant under the action of the orthogonal group is a linear combination of scalar products between the column vectors of $X$, denoted by $x_{1},\ldots,x_{k}$. In other words, there exist coefficients $m_{ij}$ for all $i,j\in[k]$ such that

[TABLE]

Proof 2.7**.**

Any homogeneous polynomial of degree two in $X$ can be written as

[TABLE]

where $\mathrm{vec}(X)$ vectorizes a matrix by stacking its columns, and $A$ is a coefficient matrix of size $dk\times dk$. Using the property $\mathrm{vec}(ABC)=(C^{T}\otimes A)\mathrm{vec}(B)$ where $\otimes$ denotes the Kronecker product, we find that

[TABLE]

If $p$ is invariant, then $p(X)=p(QX)$ for all $X$ and $Q\in\mathcal{O}(d)$. Hence, by identification:

[TABLE]

Let $A$ be a block matrix with blocks $A_{ij}\in\mathbb{R}^{d\times d}$ for $i,j$ ranging in $[k]$. The above states:

[TABLE]

for all $i,j\in[k]$ and $Q\in\mathcal{O}(d)$. This has the following consequences:

Considering each $Q$ in the set of diagonal matrices with diagonal entries in $\{-1,+1\}$ shows that all the elements off the diagonal of $A_{ij}$ are equal to their opposite, implying that $A_{ij}$ is diagonal. 2. 2.

Considering each $Q$ in the set of permutation matrices implies that all the diagonal elements are equal.

Thus, $A_{ij}=m_{ij}I_{d}$ and $A=M\otimes I_{d}$, with $M\in\mathbb{R}^{k\times k}$. As a result:

[TABLE]

*In words, $p(X)$ is a linear combination of the Gram matrix entries. *

A question that naturally arises is: what are the other invariant polynomials? The first fundamental theorem of the orthogonal group (see (Kǎc, 1994, Theorem. 14-1.2) for instance) provides an answer to this question.

Theorem 2.8**.**

The functions $g_{ij}\colon(x_{1},\ldots,x_{k})\mapsto\langle x_{i},x_{j}\rangle,1\leq i\leq j\leq k$, generate the ring of invariant polynomials, that is: any invariant polynomial is a polynomial combination of the degree two invariants $g_{ij}$.

Thus, invariant polynomials of degree higher than two do not carry further information about $[X]$. The next natural question is: are the invariants sufficient to fully characterize the equivalence classes? The following classical theorem from invariant theory provides a positive answer to this question in the case of compact groups (see (Kǎc, 1994, Theorem. 6-2.2) for instance).

Theorem 2.9**.**

(Informal) The full invariant ring characterizes the equivalence classes.

Theorems 2.8 and 2.9 combined imply that $X_{1}^{T}X_{1}=X_{2}^{T}X_{2}$ if and only if $[X_{1}]=[X_{2}]$, which is a well-known fact here derived through the prism of invariant features.

2.2 Estimation algorithm

Above, we have shown that the problem of recovering $[X]$ can be reduced (without loss) to that of estimating the Gram matrix $G=X^{T}X$, with the advantage that the latter is invariant under orthogonal transformations. We build on this observation to propose a concrete algorithm. Consider the Gram matrix of an observation as in equation. (1):

[TABLE]

By the strong law of large numbers, their empirical mean converges almost surely to their expectation (a characterization of the fluctuations for finite $N$ follows):

[TABLE]

Here, we used independence of the $Y_{i}$’s, independence of $Q_{i}$ and $E_{i}$ for each $i$, and the fact that individual entries of each $E_{i}$ are independent with mean zero and unit variance.

If the noise level $\sigma$ is known, we can get an unbiased estimator for $X^{T}X$ as

[TABLE]

Since $\hat{G}_{N}$ is expected to be close to $X^{T}X$ for large $N$, it is reasonable to consider an estimator $[\tilde{X}_{N}]$ for the equivalence class $[X]$ where $\tilde{X}_{N}$ is a solution of the optimization problem

[TABLE]

Well-known extremal properties of the eigenvalue decomposition of a symmetric matrix tell us that an optimizer $\tilde{X}_{N}$ can be obtained by computing $d$ dominant, orthonormal eigenvectors of $\hat{G}_{N}$, and scaling them by the square root of their corresponding eigenvalues. Since $\hat{G}_{N}$ and $\hat{M}_{N}$ share the same eigenvectors, with eigenvalues related through $\lambda_{\ell}(\hat{G}_{N})=\lambda_{\ell}(\hat{M}_{N})-d\sigma^{2}$, this computation can equivalently be executed from $\hat{M}_{N}$ directly. This procedure is summarized in Algorithm 1. Provided $\hat{M}_{N}$ has an eigengap separating its $d$th and $(d+1)$st largest eigenvalues, this procedure uniquely defines $[\tilde{X}_{N}]$. In Section 3, we argue that such an eigengap exists with high probability if $N$ is sufficiently large, and we bound the error $\rho([X],[\tilde{X}_{N}])$.

We note that the Gram matrix estimator $\hat{G}_{N}$ could be replaced by a more robust estimator, for example based on median-of-means as favored in (Bandeira et al., 2017a), leveraging work by [Nemirovsky and Yudin] (Nemirovsky and Yudin, 1983) and more recently by [Joly et al.] (Joly et al., 2017). Such refinements are not necessary under our assumption of Gaussian noise but could be useful for heavy-tailed noise.

2.3 Estimation when $\sigma$ is unknown

If $\sigma$ is unknown, it can be estimated from the eigenvalues of $\hat{M}_{N}$ (11). Indeed, in the limit of $N$ going to infinity, the $k-d$ smallest eigenvalues are equal to $d\sigma^{2}$. For finite $N$, they fluctuate around that value. Thus, with $\lambda_{1}\geq\cdots\geq\lambda_{k}$ the eigenvalues of $\hat{M}_{N}$, a possible estimator $\hat{\sigma}_{N}\geq 0$ for $\sigma$ is defined by

[TABLE]

where the second form is computationally favorable. The resulting procedure is summarized as Algorithm 2.

We mention that in the case where $X$ is centered, that is, when $X\bm{1}=0$, $\sigma$ can be reliably estimated by computing the empirical variance of the $dN$ i.i.d. entries of samples $\frac{1}{\sqrt{k}}Y_{i}\bm{1}\sim\mathcal{N}(0,\sigma^{2}I_{d})$ for $i=1,\ldots,N$.

3 Stability of the estimator

The consistency of $\hat{G}_{N}$ (12) as an estimator of the Gram matrix is guaranteed by the law of large numbers. However, for finite $N$, we can only hope to estimate the Gram matrix approximately: we characterize the expected errors here. From this approximate Gram matrix, we obtain our estimator $[\tilde{X}_{N}]$ by solving the optimization problem (13) (which is akin to forming a type of thin Cholesky factorization of $\hat{G}_{N}$ after projecting the latter to the positive semidefinite matrices): we call this step the Gram inversion. To understand the final error on our estimator, we need to study the sensitivity of Gram inversion: we start with this.

3.1 Sensitivity of Gram inversion

Consider the function $f([X])=X^{T}X$. This is a map from the quotient space

[TABLE]

where $\sim$ is the equivalence relation defined in (3) and $\mathbb{R}^{d\times k}_{*}$ is the set of matrices in $\mathbb{R}^{d\times k}$ of full rank $d$, to the set

[TABLE]

where $\mathrm{Sym}_{k}$ is the set of symmetric matrices of size $k$ and $G\succeq 0$ means $G$ is positive semidefinite. In Section 2, we argued that $[X]$ can be recovered uniquely from $X^{T}X$, meaning that $f$ is globally invertible. In this section, we are concerned with the sensitivity of the inverse, $f^{-1}$.

A bound on the sensitivity appears in (Tu et al., 2016, Lem. 5.4). We repeat it here.

Lemma 3.1**.**

For any $[X],[\tilde{X}]\in\mathcal{M}$,

[TABLE]

This result establishes a Lipschitz constant for $f^{-1}$ in the vicinity of $[X]$, with respect to the distance $\rho$ on $\mathcal{M}$ (6) and the Frobenius distance on $\mathcal{N}$. Through a geometric argument, we confirm that the coefficient $\sigma_{d}(X)$ in the denominator cannot be avoided, and we show $L$ must be at least $1/\sqrt{2}$ (that is $0.71..$ compared to $1.10..$ above). Then, we lean on Lemma 3.1 to obtain a new bound with (essentially) that optimal constant. To do so, we use the proof mechanics proposed by [Chang and Stehle] (Chang and Stehlé, 2010) in their study of the stability of the Cholesky decomposition of (strictly) positive definite matrices.

We equip the quotient space $\mathcal{M}$ with a smooth structure as a Riemannian quotient manifold of $\mathbb{R}^{d\times k}_{*}$ with the standard inner product $\langle\cdot,\cdot\rangle$ (Absil et al., 2008, §3.4).111In contrast, the quotient space $\mathbb{R}^{d\times k}/\!\sim$ (without rank restriction) does not admit such a smooth structure, because not all its equivalence classes have the same dimension as submanifolds of $\mathbb{R}^{d\times k}$. Likewise, we endow $\mathcal{N}$ with a smooth structure as a Riemannian submanifold of $\mathrm{Sym}_{k}$ with the standard inner product $\langle\cdot,\cdot\rangle$, as in (Vandereycken et al., 2009). With these smooth structures, $f\colon\mathcal{M}\to\mathcal{N}$ is a smooth function. Detailed background on both geometries and their relations can be found in (Massart and Absil, 2020).

The distance $\rho$ happens to be the geodesic distance on $\mathcal{M}$ (Massart and Absil, 2020). Furthermore, since $\mathcal{N}$ is a Riemannian submanifold of $\mathrm{Sym}_{k}$ equipped with the trace inner product, the Frobenius distance between close-by points of $\mathcal{N}$ is an excellent approximation for the geodesic distance between them. As a result, locally around $X^{T}X$, the operator norm (the largest singular value) of the differential of $f^{-1}$ at $X^{T}X$ reveals the local Lipschitz constant of $f^{-1}$ with respect to these distances.

The inverse function theorem (Lee, 2012, Thm. 4.5) states that the differential of $f^{-1}$ at $X^{T}X$ is the inverse of the differential of $f$ at $[X]$. Accordingly, we first study the singular values of the differential of $f$ at $[X]$.

As a preliminary step, for $X\in\mathbb{R}^{d\times k}_{*}$, consider the differential $\mathcal{L}_{X}$ of the map $X\mapsto X^{T}X$ (its relation to $f$ is elucidated below):

[TABLE]

Clearly, the following subspace is included in the kernel of $\mathcal{L}_{X}$:

[TABLE]

We can thus restrict our attention to the orthogonal complement of $\mathrm{V}_{X}$, which we denote by $\mathrm{H}_{X}=(\mathrm{V}_{X})^{\perp}$. By definition, $\dot{X}\in\mathbb{R}^{d\times k}$ is orthogonal to $\mathrm{V}_{X}$ if and only if $\langle\dot{X},\Omega X\rangle=0$ for all skew-symmetric matrices $\Omega$, hence:

[TABLE]

Restricted to $\mathrm{H}_{X}$, the smallest singular value of $\mathcal{L}_{X}$ is positive.

Proposition 3.2**.**

Given $X\in\mathbb{R}^{d\times k}_{*}$, with notation as above, consider the restriction of $\mathcal{L}_{X}$ as an operator from $\mathrm{H}_{X}$ to $\mathcal{L}_{X}(\mathrm{H}_{X})$. That operator is invertible and its smallest singular value is $\sqrt{2}\sigma_{d}(X)$.

Proof 3.3**.**

See Appendix C.

Using tools from differential geometry (Absil et al., 2008, Section 3.5.8)(Massart and Absil, 2020), it can be shown that $\mathrm{H}_{X}$ (equipped with the standard inner product) is isometric to the tangent space of $\mathcal{M}$ at $[X]$, so that the singular values of $\mathcal{L}_{X}$ on $\mathrm{H}_{X}$ are equal to the singular values of the differential of $f$ at $[X]$. Thus, calling upon the inverse function theorem, we conclude that the largest singular value of the differential of $f^{-1}$ at $X^{T}X$ is $1/\sqrt{2}\sigma_{d}(X)$. In turn, this shows the constant $L$ in Lemma 3.1 must be at least $\frac{1}{\sqrt{2}}$, and the coefficient $\sigma_{d}(X)$ cannot be removed.

Lemma 3.1 and Proposition 3.2 together allow us to show our main result regarding the local sensitivity of Gram inversion, using a technique by Chang & Stehle (Chang and Stehlé, 2010).

Theorem 3.4**.**

Consider two matrices $X,\tilde{X}\in\mathbb{R}^{d\times k}_{*}$. If their Gram matrices $G=X^{T}X$ and $\tilde{G}=\tilde{X}^{T}\tilde{X}$ are close, specifically, if

[TABLE]

then the equivalence classes must be close too:

[TABLE]

where $\rho$ is the distance between equivalence classes defined in (6).

Crucially, notice that for small $\lVert G-\tilde{G}\rVert_{\mathrm{F}}$, we have $1-\sqrt{1-\frac{2\lVert G-\tilde{G}\rVert_{\mathrm{F}}}{\sigma_{d}^{2}(X)}}\approx\frac{1}{\sigma_{d}^{2}(X)}\lVert G-\tilde{G}\rVert_{\mathrm{F}}$. Hence, the right-hand side of (22) behaves like $\frac{1}{\sqrt{2}\sigma_{d}(X)}\lVert G-\tilde{G}\rVert_{\mathrm{F}}$, which by the discussion above cannot be improved.

Proof 3.5**.**

Let $U\Sigma V^{T}$ be the SVD of $X\tilde{X}^{T}$. The orthogonal matrix $Q=UV^{T}$ optimally aligns $X$ and $\tilde{X}$, in the sense that $\rho([X],[\tilde{X}])=\lVert X-Q\tilde{X}\rVert_{\mathrm{F}}$. Then,

[TABLE]

Since the theorem statement depends on $X$ and $\tilde{X}$ only through $[X]$ and $[\tilde{X}]$, without loss of generality, suppose that $X$ and $\tilde{X}$ are already rotationally aligned, that is, $Q=I$. Then, $X\tilde{X}^{T}$ is symmetric, positive semidefinite. Define $\Delta X=\tilde{X}-X$: notice that $\Delta X{X}^{T}$ is also symmetric, and

[TABLE]

Rearranging, we get:

[TABLE]

where $\mathcal{L}_{X}$ is the operator defined in (19). Since $\Delta X$ is in the subspace $\mathrm{H}_{X}$ defined in (21), it is also orthogonal to the null space $\mathrm{V}_{X}$ of $\mathcal{L}_{X}$. As a result, we may write

[TABLE]

where $\mathcal{L}_{X}^{\dagger}$ is the Moore–Penrose pseudo-inverse of $\mathcal{L}_{X}$. Proposition 3.2 then implies:

[TABLE]

Reorganizing, we get the following inequality:

[TABLE]

The right-hand side of this inequality is a quadratic in $\lVert\Delta X\rVert_{\mathrm{F}}$. Under our assumptions, the two roots of this quadratic, $\xi_{-}$ and $\xi_{+}$, are real and nonnegative:

[TABLE]

Since $\lVert\Delta X\rVert_{\mathrm{F}}$ satisfies (23), it lies outside the open interval defined by $(\xi_{-},\xi_{+})$. This means there are two possibilities: either $\lVert\Delta X\rVert_{\mathrm{F}}\in[0,\xi_{-}]$, or $\lVert\Delta X\rVert_{\mathrm{F}}\geq\xi_{+}$. We aim to exclude the latter. To do so, notice that Lemma 3.1 together with our proximity assumption on the Gram matrices implies:

[TABLE]

*This allows to conclude that $\rho([X],[\tilde{X}])=\lVert\Delta X\rVert_{\mathrm{F}}\leq\xi_{-}$. Upon factoring out $\frac{\sigma_{d}(X)}{\sqrt{2}}$ in the expression for $\xi_{-}$, this completes the proof. *

3.2 Upper bounds on cloud estimation error

Theorem 3.4 quantifies how a good estimator for the Gram matrix of $[X]$ can be turned into a good estimator for $[X]$ itself. In this part, we first show quantitatively that, with high probability, we can indeed have a good estimator for the Gram matrix. Afterwards, we connect this result with the above theorem to produce a bound on the estimation error of $[X]$.

The first result relies on standard concentration bounds for quadratic forms of Gaussian random variables. (We remark that it is possible to relax the assumptions to require subgaussian noise rather than Gaussian noise.) We assume $\sigma$ is known, and we use the notation $\lVert X\rVert_{\mathrm{op}}=\sigma_{1}(X)$ for the operator norm. The projection to the set of positive semidefinite matrices of rank at most $d$ is necessary to apply Theorem 3.4, and causes no difficulties in practice.

Theorem 3.6**.**

Let $Y_{1},\ldots,Y_{N}$ be i.i.d. observations drawn from model (1) and let $\hat{G}_{N}$ be the Gram estimator as defined in (12). Since the true Gram matrix $G=X^{T}X$ is positive semidefinite with rank at most $d$, project $\hat{G}_{N}$ to that set; in the notation of Algorithm 1:

[TABLE]

Then, for any $\delta\in(0,1)$, with probability at least $1-\delta$:

[TABLE]

Proof 3.7**.**

For all $u\in\mathbb{S}^{k-1}$ (the unit sphere in $\mathbb{R}^{k}$), consider:

[TABLE]

Notice that $\|(Q_{i}X+\sigma E_{i})u\|^{2}=\|(X+\sigma Q_{i}^{T}E_{i})u\|^{2}$ is equal in distribution to $\|(X+\sigma E_{i})u\|^{2}$ since entries of $E_{i}$ are standard Gaussian. Thus, the first part of (24) is distributed like a sum of squared norms of i.i.d. non-centered Gaussian vectors. Reorganizing standard concentration bounds for noncentral $\chi^{2}$ (see for instance (Birgé, 2001, Lemma 8.1)) gives, for all $u$ in $\mathbb{S}^{k-1}$,

[TABLE]

We then cover the sphere $\mathbb{S}^{k-1}$ with an $\varepsilon$-net. A union bound over all the elements of the net (see for instance (Vershynin, 2012, Section 5.2.2)) yields

[TABLE]

Since $G$ is positive semidefinite and has rank $d$, the projection $\tilde{G}_{N}$ satisfies:

[TABLE]

By the triangle inequality, we obtain:

[TABLE]

Since $\mathrm{rank}(G),\mathrm{rank}(\tilde{G}_{N})\leq d$, it follows that $\mathrm{rank}(G-\tilde{G}_{N})\leq 2d$ and we get:

[TABLE]

As a result,

[TABLE]

Taking $t=\log\left(\frac{2\cdot 5^{k}}{\delta}\right)$ implies the final result:

[TABLE]

with probability at least $1-\delta$.

Combining Theorem 3.6 with a stability result on the thin Cholesky decomposition gives the main result on the stability of the proposed estimator, as measured with the distance $\rho$ (6).

Corollary 3.8**.**

Let $Y_{1},\ldots,Y_{N}$ be $N$ i.i.d. samples drawn according to (1). Let $[\tilde{X}]$ be the estimator returned by Algorithm 1. Then, for large $N$, with probability at least $1-\delta$,

[TABLE]

where $L$ can be taken as $1/\sqrt{2(\sqrt{2}-1)}$ or, for large enough $N$, arbitrarily close to $1/\sqrt{2}$.

Proof 3.9**.**

*The first claim follows from Theorem 3.6 and Lemma 3.1. For $N$ large enough, Theorem 3.6 shows that the assumption of Theorem 3.4, namely, $\lVert\tilde{G}_{N}-G\rVert_{\mathrm{F}}\leq\frac{\sigma_{d}^{2}(X)}{2}$, is satisfied with high probability. In that scenario, combining the two theorems yields the second result. *

4 Statistical optimality of the estimator

Estimating the equivalence class of $X$ from samples of the form (1) is a particular instance of an estimation problem under a group action. Bandeira et al. (Bandeira et al., 2017a) showed that the statistical complexity of estimation problems under a group action is connected to the structure of the group acting on the parameter. In particular, it is shown that the minimum number of samples required to reliably estimate the parameter grows as $O(\sigma^{2p})$ where $p$ is the smallest degree of invariant polynomials required to fully characterize the equivalence classes. We have shown in Section 2 that, in our case, $p=2$. In this section, we build on those results to show minimax lower bounds on the estimation of the equivalence class $[X]$. We also provide matching upper bounds, hence showing the statistical optimality of our estimator in the low SNR regimes.

To fix scale and to avoid pathological cases, throughout this section we assume that $X$ belongs to the space

[TABLE]

where $\sigma_{d}(X)$ is the $d$th (that is, smallest) singular value of $X$, and $\eta>0$ is fixed.

4.1 Lower bound on the estimation error for high noise regimes

In the presence of large noise (that is, for large $\sigma$), the MSE of any estimator of the equivalence class of $X$ scales with $\sigma$ as $\sigma^{4}$: we make this precise in the following theorem.

Theorem 4.1**.**

Suppose we observe $N$ samples $Y_{1},\ldots,Y_{N}$ drawn independently according to (1). Then the so-called minimax risk for estimating the orbit of $X$ satisfies

[TABLE]

for sufficiently large $\sigma$, where $\rho$ is as in (6) and $\mathcal{X}$ is defined by (27), and the infimum is taken over all possible estimators, random or deterministic.

The general study of lower bounds for estimation under a group action has been addressed by Abbe et al (Abbe et al., 2018b), using Chapman–Robbins bounds. Here, we take a different approach, proving minimax rates using two ingredients: a tight bound on the Kullback–Leibler (KL) divergence and a bound on the packing number of a particular metric space. However, since we are not aware of any result on the packing number of our parameter space for the metric $\rho$, we consider a strict subset of this parameter space for which tight bounds on the packing number are known. Specifically, we start by noticing that the Grassmannian $\mathrm{G}(k,d)$—the set of $d$ dimensional subspaces of $\mathbb{R}^{k}$—is in correspondence with a subset of the parameter space. Building on this observation, we restrict ourselves to the strictly simpler problem where the matrix $X$ satisfies the condition $XX^{T}=I_{d}$. In this restricted setting, there is a one-to-one correspondence between an equivalence class and an element of the Grassmannian. We then use a result on the covering number of the Grassmannian to control its local packing number. This result was used by Cai et al. (Cai et al., 2013) in the context of optimal rates of estimation for the principal subspace of a covariance matrix under a sparsity assumption. Then, we show a tight bound of the KL divergence. This bound is a particular case of a more general result on the KL divergence of samples observed under the action of a group.

Packing number of the Grassmannian: we start by defining a metric on $\mathrm{G}(k,d)$. This metric on the Grassmanian is shown to be equivalent to the distance between equivalence classes $\rho$. We then show tight lower and upper bounds on the covering number of $\mathrm{G}(k,d)$ for the aforementioned metric. This result is due to [Szarek] (Szarek, 1982). We then use a result on the local packing number of the Grassmannian: leveraging the previous result, for any $\alpha\in(0,1)$ and for any $\varepsilon$ small enough, we give a lower bound on the $\alpha\varepsilon$-packing number of a ball of radius $\varepsilon$. This follows the technique proposed by Yang and Barron (Yang and Barron, 1999). 2. 2.

Tight control of the KL divergence: we state a lemma giving a bound on the $\mathrm{KL}$ divergence of the distribution $\mathbb{P}_{X}$ of samples drawn according to (1). This lemma follows from a result by Bandeira et al. (Bandeira et al., 2017a).

We start by giving the result on the covering number of the Grassmannian.

Lemma 4.2**.**

[Cai et al. (Cai et al., 2013), Lemma 1] Define the metric on $\mathrm{G}(k,d)$ by $\tilde{\rho}(V,U)=\lVert V^{T}V-U^{T}U\rVert_{\mathrm{F}}$. Then, for any $\varepsilon\in\left(0,\sqrt{2\min(d,k-d)}\right]$, we have

[TABLE]

where $\mathfrak{N}(\mathrm{G}(k,d),\varepsilon)$ is the $\varepsilon$-covering number of $\mathrm{G}(k,d)$ with respect to the metric $\tilde{\rho}$ and $c_{0}$, $c_{1}$ are absolute constants.

Building on the previous result, we can state a result on the local packing of $G(k,d)$.

Lemma 4.3**.**

Let $B(V,\varepsilon)=\{U\in\mathrm{G}(k,d):\tilde{\rho}(U,V)\leq\varepsilon\}$, $\alpha\in(0,1)$ and $\varepsilon\in(0,\varepsilon_{0}]$. Then, there exists $V^{*}\in\mathrm{G}(k,d)$ such that:

[TABLE]

where $\mathfrak{M}(\mathrm{G}(k,d),\varepsilon)$ is the $\varepsilon$-packing number of $\mathrm{G}(k,d)$ with respect to the metric $\tilde{\rho}$ and $c_{0}$, $c_{1}$ are the absolute constants in Lemma 4.2.

Proof 4.4**.**

The original proof can be found in (Yang and Barron, 1999). For convenience, we provide the proof in Appendix D.

Finally, we state a bound on the KL divergence of the distribution $\mathbb{P}_{X}$ of samples drawn according to (1).

Lemma 4.5**.**

For given dimensions $d$ and $k$, there exists a universal constant $C$ such that, for any $X_{1}$, $X_{2}\in\mathbb{R}^{d\times k}$ with $\rho(X_{1},X_{2})\leq\dfrac{\lVert X_{1}\rVert_{\mathrm{F}}}{3}$ and for any $\sigma>1$, we have:

[TABLE]

Proof 4.6**.**

See Appendix E.

Combining these three results, we obtain Theorem 4.1. The proof technique is inspired from (Yang and Barron, 1999) and (Cai et al., 2013). The full proof can be found in Appendix D. The constant in the lower bound exhibits a dependence in the dimension, through the metric entropy of the Grassman manifold. As we focus on the dependence in $\sigma$ and $N$ (noise level and number of samples), we leave the question of whether this dependence in the dimension is tight for future research.

4.2 Lower bound on the sample complexity for high noise regimes

We now show that, still when $\sigma$ is large, the number of samples drawn according to (1) necessary to reliably estimate the equivalence class of $[X]$ grows as $O(\sigma^{4})$. For that matter, we provide matching upper and lower bounds: the upper bound on the KL divergence obtained in Lemma 4.5 in combination with Neyman–Pearson’s lemma (see [ (Rigollet and Hutter, 2017), Lemma 4.3] for instance) yields the lower bound, while properties of the estimator obtained with Algorithm 1 give the upper bound. This result is similar in its message to, but technically different from, Theorem 4.1.

Theorem 4.7**.**

Let $\tau\in(0,1/2)$ be arbitrary. There exists a constant $\tilde{c}$ (possibly function of $\tau$) such that the following holds: let $X_{1},X_{2}\in\mathcal{X}$ belong to two distinct equivalence classes, and let $Y_{1},\ldots,Y_{N}$ be drawn according to (1), where $X$ is either $X_{1}$ or $X_{2}$, with probability 1/2 each. Any test $\psi$ whose task it is to decide whether $X$ is $X_{1}$ or $X_{2}$ has probability of error at least $1/2-\tau$ whenever $N<\tilde{c}\sigma^{4}$.

Moreover, there exists a constant $\tilde{C}$ such that Algorithm 1 outputs $\hat{X}_{N}$ satisfying

[TABLE]

whenever $N\geq\tilde{C}\dfrac{\sigma^{4}}{\delta\cdot\varepsilon^{2}}.$

Proof 4.8**.**

We start by proving that no statistical procedure can reliably distinguish between two equivalence classes whenever the number of samples is less than $O(\sigma^{4})$. In particular, for any test $\psi$ using $N$ samples we have:

[TABLE]

where the first inequality follows from Neyman–Pearson’s Lemma while the second follows from Pinsker’s inequality. Therefore, if $N\leq\dfrac{8}{C\rho^{2}(X_{1},X_{2})}\tau^{2}\sigma^{4}$, we get that:

[TABLE]

yielding the sought lower bound on the probability of error. The constant $\tilde{c}$ can be set uniformly against the choice of $X_{1},X_{2}$ since $\mathcal{X}$ is bounded.

For the upper bound, we simply notice that, for large $\sigma$, by Markov’s inequality,

[TABLE]

where the first inequality follows from Markov’s inequality while the second inequality follows from the computations in Appendix F.2.

4.3 A comment regarding the special orthogonal group

In our model, we consider the case where the orthogonal transformations acting on the cloud are in $\mathcal{O}(d)$. In some applications, it may be more appropriate to assume rotations $Q_{1},\ldots,Q_{N}$ in the special orthogonal group $\mathcal{SO}(d)$, that is, with determinant $+1$ (no reflections). However, the problem in $\mathcal{SO}(d)$ is harder. Indeed, by the first fundamental theorem of invariant theory for the special orthogonal group, the entries of the Gram matrix $X^{T}X$ and the $d\times d$ minors of $X$ generate the ring of invariant polynomials (see (Kraft and Procesi, 1996, Prop. 10.2) for instance). Since our algorithm only recovers the Gram matrix, that is, the orbits of matrix $X$ up to a reflection, it falls short of estimating such reflection. Estimating the reflection can be done by estimating the determinants of the submatrices of $X$ of size $d\times d$, which requires $O(\sigma^{2d})$ samples at high noise level, as further discussed in (Bandeira et al., 2017a).

4.4 Noise regimes

We here summarize the adaptive optimality property of our estimator, characterized by a phase transition around a critical noise level.

Proposition 4.9**.**

The estimator in Algorithm 1 exhibits a phase transition around a critical noise level, namely:

• –

If $\sigma^{2}\gg 1$, the MSE of the estimator behaves like $O\left(\sigma^{4}/N\right)$, which is optimal by Theorem 4.1.

• –

If $\sigma^{2}\ll 1$, the MSE of the estimator behaves like $O\left(\sigma^{2}/N\right)$, which is optimal given the fact that the problem is strictly harder than estimating the matrix $X$ provided the rotations are available, in which case the lower bound on the estimation is $O(\sigma^{2}/N)$.

Proof 4.10**.**

*The proof of the upper bound relies on a computation of bounds on the MSE of the estimator proposed in Algorithm 1: see Appendix F. *

5 Numerical experiments

In this section, we provide numerical support for our theoretical predictions.222Code to generate the figures: https://github.com/thomaspdl/GeneralizedProcrustes. In particular, we highlight the adaptative optimality of Algorithm 1 and validate the heuristic proposed for the estimation of $\sigma^{2}$.

5.1 Noise regimes

Figure 1 illustrates the accuracy of Algorithm 1 over a large range of values of $N$ (number of observations) and $\sigma$ (noise level), with $k=100$ and $d=3$ fixed. A clear phase transition is visible, delineating a regime where our estimator is accurate, and one where it fails. Crucially, the phase transition illustrates the adaptiveness of the estimator. Specifically, as predicted, for small $\sigma$, the number of observations needs to grow as $N\sim\sigma^{2}$ in order to preserve a constant MSE while compensating for noise. For large $\sigma$, this relationship deteriorates and we require $N\sim\sigma^{4}$. As explained in Section 4, this is not merely a requirement of our estimator: any estimator is subject to the same needs, and it is a positive feature of Algorithm 1 that it adapts automatically to both regimes.

To illustrate these two noise regimes, the graph is overlaid with two lines. The red line of slope $1/2$ represents a noise level varying as $O(\sqrt{N})$. Its intercept is chosen as follows: consider an oracle which knows the rotations $Q_{i}$ affecting the measurements $Y_{i}$. This oracle can compute the MLE simply by undoing the rotations, then averaging. It is easy to see that the MSE of that oracle is $\frac{\sigma^{2}dk}{N}$ in Frobenius distance (it may be slightly less in $\rho$ distance). The red line shows the relationship between $N$ and $\sigma$ when that oracle has a relative MSE of $0.95$. It is interesting to see that, for low noise levels, the invariants-based estimator has performance similar to that oracle. The blue line of slope $1/4$ represents a noise level varying as $O(N^{\frac{1}{4}})$: its intercept is chosen manually, for illustration.

5.2 Estimation of $\sigma^{2}$

To illustrate the accuracy of the approximation made in (14), we compute the mean empirical error of our estimator. In particular, for different numbers of observations $N$, we estimate the average relative error of the estimator reported in Fig. 2. We observe a small empirical relative error decreasing with the number of samples.

6 Perspectives

Throughout this paper, we only consider the case where one cloud of points must be estimated. An interesting variant is heterogeneous Procrustes, where each observation consists of a noisy, rotated version of a cloud of points picked at random among $K$ possibilities, $X_{1},\ldots,X_{K}$; for example:

[TABLE]

where $w_{1},\ldots,w_{K}\geq 0$ sum to 1. Several approaches to solve this problem are possible: taking a moments-based method similar to the one proposed by Hsu and Kakade (Hsu and Kakade, 2013) or inverting invariant features using a nonconvex optimization approach as in (Boumal et al., 2018; Ma et al., 2020) could be two interesting directions. Regarding the latter, one idea is to consider these fourth-order $\mathcal{O}(d)$-invariants: $Y^{T}Y\otimes Y^{T}Y$, where $\otimes$ is the Kronecker product.

Another possible extension is to consider observations with unlabeled clouds of points, that is: we observe several copies of $X$ after an unknown rotation, and also after an unknown permutation of the points to model the fact that we do not know which point is which across various observations. This can be modeled as

[TABLE]

where $Q$ is an unknown orthogonal matrix of size $d\times d$ as usual, and $P$ is an unknown permutation matrix of size $k\times k$. Assuming i.i.d. Gaussian noise, in distribution, $Y=Q(X+\sigma E)P$, so that the singular values of $Y$ are equal to those of $X+\sigma E$: up to the noise, the singular values are (non-polynomial) invariants. Studying the distribution of the singular values of $Y$ as a function of those of $X$ may allow one to estimate the general shape of the cloud $X$ (specifically, its singular values) without the need to register clouds (no rotation estimation, and no point correspondence estimation). This connects to principal component analysis.

Yet another extension is to consider projected observations. For example, $X$ is a cloud of points in three dimensions ($d=3$), but observations are of the form $Y=P(QX+\sigma E)$, where $P$ is a projector to a two-dimensional plane (for example, the camera plane). Equivalently,

[TABLE]

where $U$ is an unknown matrix of size $3\times 2$ with orthonormal columns, and $E$ has size $2\times k$. Assuming $U$ is uniformly distributed, $\mathbb{E}[UU^{T}]=\frac{2}{3}I_{3}$, so that $Y^{T}Y$, while not an invariant, does reveal information about $X$ in expectation since $\mathbb{E}[Y^{T}Y]=\frac{2}{3}X^{T}X+d\sigma^{2}I_{k}$. Much of the machinery developed in the present paper applies directly to this extended setting: it would be interesting to study it further.

In all of these, one could also include the possibility that clouds are not centered, that is: they are observed only after an unknown translation, on top of other transformations.

Acknowledgments

We thank Afonso Bandeira, Adrien Bilal, Jianqing Fan, Joe Kileel and Joao Pereira for insightful discussions.

Funding

NB is partially supported by NSF award DMS-1719558.

Data Availability Statement

The data underlying this article are available on GitHub, at https://github.com/thomaspdl/GeneralizedProcrustes.

Appendix A Procrustes distance

Here, we show how to explicitly compute the distance $\rho$ (6). Since:

[TABLE]

finding the optimal $Q$ is equivalent to solving:

[TABLE]

By SVD we can write $X_{1}X_{2}^{T}=U\Sigma V^{T}$ and the solution of the problem is given by the polar factor of $X_{1}X_{2}^{T}$:

[TABLE]

Appendix B Impossibility of estimation for the rotations

Here, we give a proof of Proposition (1.1). We start by reminding the statement of the proposition below.

Proposition B.1**.**

Let us consider the following hypothesis testing problem on the distribution of the unknown distribution $\mathbb{P}$ of the sample:

[TABLE]

Then for any test $\psi$ and for any precision $\delta>0$, it is possible to find a sufficiently large noise level $\sigma_{0}$ such that, for any noise level $\sigma$ larger than $\sigma_{0}$, the sum of type I and type II errors is large:

[TABLE]

Hence, if the two hypotheses are equally likely; the probability of error is at least $1/2-\delta$.

Proof B.2**.**

The likelihood ratio test consists of studying the ratio of the densities $L(Y)=\dfrac{f_{1}(Y)}{f_{2}(Y)}$ given observation $Y$. If $L(Y)>1$, then $H_{1}$ is kept; otherwise, $H_{2}$ is kept. Hence, $\mathbb{P}_{1}(\psi^{*}=2)=\mathbb{P}_{1}(L(Y)<1)$.

[TABLE]

This yields:

[TABLE]

Yet,

[TABLE]

with $\mu=\lVert X\rVert_{\mathrm{F}}^{2}$, $\tilde{\sigma}=\sigma\lVert X\rVert_{\mathrm{F}}$ and $\xi\sim\mathcal{N}(0,1)$. Hence,

[TABLE]

where $F$ is the cumulative density function of a standard Gaussian random variable. Similarly, we get:

[TABLE]

Hence, by Neyman–Pearson’s lemma (Neyman and Pearson, 1933), it is sufficient to take $\sigma$ such that $F\left(-\dfrac{\lVert X\rVert_{\mathrm{F}}}{\sigma}\right)=\dfrac{1-2\delta}{2}$ to get that for any test $\psi$:

[TABLE]

where $F$ is the cumulative density function of a standard normal distribution.

Appendix C Stability of the Cholesky decomposition

We give a proof of Proposition 3.2 which is used in the proof of Theorem (3.4) about the stability of the estimator given by Algorithm 1.

Proof C.1**.**

Let $X=U\Sigma V^{T}$ be the thin SVD of $X$. Here, $V\in\mathbb{R}^{k\times d}$ has orthonormal columns, $U$ is an orthogonal matrix of size $d\times d$ and $\Sigma$ is a diagonal matrix of size $d$ with entries $\sigma_{1}\geq\cdots\geq\sigma_{d}>0$. It is always possible to pick $V_{\perp}$ (a complement of the orthonormal basis $V$) such that $\begin{bmatrix}V&V_{\perp}\end{bmatrix}$ is an orthogonal matrix of size $k\times k$. Hence, for any $\dot{X}\in\mathbb{R}^{d\times k}$, there exist matrices $A,B$ of appropriate size such that

[TABLE]

Using this parameterization,

[TABLE]

We use this to derive an SVD of $\mathcal{L}_{X}$ restricted to $\mathrm{H}_{X}$, that is, for $\dot{X}$ such that $A\Sigma=\Sigma A^{T}$. To this end, consider the orthonormal basis of $\mathrm{H}_{X}$ composed of the following elements, and its image through $\mathcal{L}_{X}$ (below, $e_{i}$ denotes the $i$th column of the identity matrix of appropriate dimension as indicated by context):

*With $A=e_{i}e_{i}^{T},B=0$ for $1\leq i\leq d$, *

$\mathcal{L}_{X}(\dot{X})=2\sigma_{i}\begin{bmatrix}V&V_{\perp}\end{bmatrix}\begin{bmatrix}e_{i}e_{i}^{T}&0\\ 0&0\end{bmatrix}\begin{bmatrix}V&V_{\perp}\end{bmatrix}^{T}$; 2. 2.

*With $A=\frac{\sigma_{i}e_{i}e_{j}^{T}+\sigma_{j}e_{j}e_{i}^{T}}{\sqrt{\sigma_{i}^{2}+\sigma_{j}^{2}}},B=0$ for $1\leq i<j\leq d$, *

$\mathcal{L}_{X}(\dot{X})=\sqrt{2(\sigma_{i}^{2}+\sigma_{j}^{2})}\begin{bmatrix}V&V_{\perp}\end{bmatrix}\begin{bmatrix}\frac{e_{i}e_{j}^{T}+e_{j}e_{i}^{T}}{\sqrt{2}}&0\\ 0&0\end{bmatrix}\begin{bmatrix}V&V_{\perp}\end{bmatrix}^{T}$; 3. 3.

*With $A=0,B=e_{i}e_{j}^{T}$ for $1\leq i\leq d$ and $1\leq j\leq k-d$, *

$\mathcal{L}_{X}(\dot{X})=\sqrt{2}\sigma_{i}\begin{bmatrix}V&V_{\perp}\end{bmatrix}\begin{bmatrix}0&\frac{e_{i}e_{j}^{T}}{\sqrt{2}}\\ \frac{e_{j}e_{i}^{T}}{\sqrt{2}}&0\end{bmatrix}\begin{bmatrix}V&V_{\perp}\end{bmatrix}^{T}$.

To verify that the inputs yield an orthonormal basis of $\mathrm{H}_{X}$ as announced, check that each of these choices yields a matrix $\dot{X}$ in $\mathrm{H}_{X}$; they are indeed orthonormal; and they are in sufficient number to cover $\dim\mathrm{H}_{X}=dk-\dim\mathrm{V}_{X}=dk-\frac{d(d-1)}{2}$. As the outputs of $\mathcal{L}_{X}$ applied to the basis elements are also orthogonal, it is clear that the singular values of $\mathcal{L}_{X}$ restricted to $\mathrm{H}_{X}$ are:

$2\sigma_{1},\ldots,2\sigma_{d}$; 2. 2.

$\sqrt{2(\sigma_{i}^{2}+\sigma_{d}^{2})}$* for $1\leq i<j\leq d$; and* 3. 3.

$\sqrt{2}\sigma_{1},\ldots,\sqrt{2}\sigma_{d}$, each repeated $k-d$ times.

(And of course, the singular values of $\mathcal{L}_{X}$ on $\mathrm{V}_{X}$ are zero, $d(d-1)/2$ times.)

Appendix D Minimax lower bound

D.1 Proof of Lemma 4.3

Proof D.1**.**

Let us consider $G_{\varepsilon}$, a minimal $\varepsilon$-cover of $\mathrm{G}(k,d)$ for the semi metric $\tilde{\rho}$, i.e. such that:

[TABLE]

Then:

[TABLE]

For contradiction, assume that for all $U$ in $G_{\varepsilon}$ we have

[TABLE]

This implies:

[TABLE]

which contradicts (29). Hence, there exists $U^{*}\in\mathrm{G}_{\varepsilon}$ such that:

[TABLE]

where the second inequality follows from Lemma 4.2. The fact that $\mathfrak{M}(E,\varepsilon)\geq\mathfrak{N}(E,\varepsilon)$ for any $E$ allows to conclude the proof.

D.2 Proof of Theorem 4.1

We will use the following lemma, originally given in (Tu et al., 2016).

Lemma D.2**.**

For any $X_{2}\in\mathbb{R}^{d\times k}$ obeying $\rho(X_{2},X_{1})\leq\dfrac{1}{4}\lVert X_{1}\rVert_{\mathrm{op}}$, we have

[TABLE]

Proof D.3**.**

For all $Q$ in $\mathcal{O}(d)$ we have

[TABLE]

Taking the infimum of the right-hand side over $\mathcal{O}(d)$ gives the result.

Proof D.4** (Proof of Theorem 4.1).**

We restrict ourself to equivalence classes $[X]$ such that $XX^{T}=I_{d}$. There is a one-to-one mapping between such equivalence classes and elements of $\mathrm{G}(k,d)$, since $X^{T}X$ is then the orthogonal projector to the space spanned by the columns of $X^{T}$ (see (Absil et al., 2008, Section 3.4.4) for instance). We get that:

[TABLE]

where $\tilde{\rho}$ is the metric defined in Lemma 4.2. Lemma 4.3 allows us to consider a point $U^{*}$ and $m$ points $\{U_{1},\ldots,U_{m}\}\subset B(U^{*},\varepsilon)$, such that $\alpha\varepsilon\leq\tilde{\rho}(U_{i},U_{j})\leq 2\varepsilon$ for all $i\neq j$ and $m\geq\left(\dfrac{c_{0}}{\alpha c_{1}}\right)^{d(k-d)}$. Without loss of generality, we can assume that $\varepsilon\leq\dfrac{1}{2}$. By standard arguments (see for instance chapter 4 of (Rigollet and Hutter, 2017)) we can lower bound the minimax risk by the probability of error for testing a finite number of hypothesis:

[TABLE]

Moreover, by Fano’s inequality (Fano, 1952) (see Theorem 4.19 in (Rigollet and Hutter, 2017)) we can lower bound the probability of error for the multiple hypothesis testing problem $\{\mathbb{P}_{U_{i}}:i\in[m]\}$ by

[TABLE]

By Lemma 4.5, for all pair $i,j$ we have that:

[TABLE]

where by hypothesis $\tilde{\rho}(U_{i},U_{j})\leq 2\varepsilon$. Theorem 3.4 implies that for small enough $\varepsilon$, $\rho(U_{i},U_{j})\leq\sqrt{2}\varepsilon$. Hence for small enough $\varepsilon$ and large enough $\sigma$ we have that

[TABLE]

*Hence, the minimax error admits the following lower bound: *

[TABLE]

for any $\varepsilon\in(0,\varepsilon_{0}]$ and $\alpha\in(0,1)$. By picking

[TABLE]

if $\varepsilon_{0}^{2}\geq\dfrac{\sigma^{4}}{N}\cdot\dfrac{d(k-d)\log(2)}{6C^{{}^{\prime}}}$, and $\varepsilon^{2}=\dfrac{\varepsilon_{0}^{2}}{2}$ otherwise, we get the following inequality:

[TABLE]

which yields the lower bound. The upper bound directly follows from the upper bound on the MSE of the estimator, which can be found in Appendix F.2 .

Appendix E Bound on the KL divergence

In this section, we prove the tight bound on the $\mathrm{KL}$ divergence of the distribution of samples (1) stated in Lemma 4.5. Our result is a direct consequence of Proposition 7.8 in (Bandeira et al., 2017b).

Proof E.1**.**

By the properties of the orthogonal group and the Haar measure, we have:

[TABLE]

We now bound KL divergence between $\mathbb{P}_{X_{1}}$ and $\mathbb{P}_{X_{2}}$ by bounding the $\chi^{2}$-divergence between $\mathbb{P}_{X_{1}}$ and $\mathbb{P}_{X_{2}}$. The density of $\mathbb{P}_{X}$, $f_{X}$ can be written as:

[TABLE]

where $f$ is the density of a standard $dk$-dimensional Gaussian. Since $\lVert X\rVert_{\mathrm{F}}^{2}\leq d$ by hypothesis, we obtain by Jensen’s inequality:

[TABLE]

Hence the $\chi^{2}$ divergence can then be bounded as:

[TABLE]

where $E$ is a $d\times k$ matrix with i.i.d. entries $\sim\mathcal{N}(0,1)$. The first inequality comes from applying inequality (30) to the denominator, and the second equality comes from a change of variables. By Fubini’s theorem we can change the order of expectations. Using the fact that for a standard Gaussian $U$, for any $\lambda$, we have $\mathbb{E}[e^{\lambda U}]=e^{\frac{\lambda^{2}}{2}}$, we obtain that:

[TABLE]

We use this to expand the square in the above expression for $\chi^{2}(\mathbb{P}_{X_{1}},\mathbb{P}_{X_{2}})$:

[TABLE]

Since $Q_{1}$,$\tilde{Q}_{1}$, $Q_{2}$ and $\tilde{Q}_{2}$ are drawn according to the Haar measure on the orthogonal group, according to Fubini’s theorem we can expand each term as a power series and exchange summation and expectation. Using the fact that $\langle u,v\rangle^{l}=\langle u^{\otimes l},v^{\otimes l}\rangle$ and that for all i.i.d. vectors $x$ and $y$ we have that

[TABLE]

we obtain :

[TABLE]

with $\Delta_{l}=\mathbb{E}_{Q}[\mathrm{vec}(QX_{1})^{\otimes l}-\mathrm{vec}(QX_{2})^{\otimes l}]$, where the expectation is taken over the Haar measure. In particular,

[TABLE]

Moreover, by Lemma E.2 below, for all $l\geq 2$ we have:

[TABLE]

Hence, for any $\sigma>1$ we obtain:

[TABLE]

The inequality $\mathrm{KL}(\mathbb{P}_{X_{1}}||\mathbb{P}_{X_{2}})\leq\chi^{2}(\mathbb{P}_{X_{1}},\mathbb{P}_{X_{2}})$ (see Lemma 2.2 in (Tsybakov, 2009) for instance) finishes the proof of the claim.

The general structure of the proof of the lemma below is as in Lemma B.12 of (Bandeira et al., 2017b), with adaptations as needed.

Lemma E.2**.**

Given $X_{1},X_{2}\in\mathbb{R}^{d\times k}$ such that $\rho(X_{1},X_{2})<\dfrac{\lVert X_{1}\rVert_{\mathrm{F}}}{3}$, for all $l\geq 1$,

[TABLE]

Proof E.3**.**

Without loss of generality we assume that $X_{1}$ and $X_{2}$ are rotationally aligned, i.e., $\rho(X_{1},X_{2})=\lVert X_{1}-X_{2}\rVert_{\mathrm{F}}=\varepsilon\lVert X_{1}\rVert_{\mathrm{F}}$, with $\varepsilon<\dfrac{1}{3}$. By Jensen’s inequality,

[TABLE]

Expanding the norm, it is easy to check that

[TABLE]

where $\gamma=\dfrac{\langle X_{1},X_{2}-X_{1}\rangle}{\lVert X_{1}\rVert_{\mathrm{F}}^{2}}$. By Cauchy–Schwarz we have $|\gamma|\leq\varepsilon<\dfrac{1}{3}$. Moreover, $2\gamma+\varepsilon^{2}\leq 3\varepsilon<1$. By the binomial theorem, for all $x$ such that $|x|<1$, there exists an $r_{l}\leq 2^{l}x^{2}$ such that

[TABLE]

Hence,

[TABLE]

and $\lVert\Delta_{l}\rVert_{\mathrm{F}}^{2}\leq\lVert X_{1}\rVert_{\mathrm{F}}^{2l}12(2)^{l}\varepsilon^{2}=\lVert X_{1}\rVert_{\mathrm{F}}^{2l-2}12(2)^{l}(\lVert X_{1}\rVert_{\mathrm{F}}\varepsilon)^{2}\leq 12(2d)^{l}\rho^{2}(X_{1},X_{2}).$

Appendix F Statistical properties of the Estimator

We now compute an upper bound on the MSE of our estimator in order to prove the second part of Proposition 4.9. We start by computing the covariance matrix of the estimator in part F.1, before proving the upper bound in part F.2.

F.1 Covariance matrix of the estimator $\hat{M}_{N}$

It is straightforward to get that

[TABLE]

We now compute the covariance matrix $\Sigma$ of $\hat{M}_{N}$ i.e.

[TABLE]

For that purpose, we define $F_{N}$ as

[TABLE]

Each of the entry of $F_{N}$ can be written as

[TABLE]

where X=\left[\begin{array}[]{cccc}\vrule&\vrule&&\vrule\\ x_{1}&x_{2}&\ldots&x_{k}\\ \vrule&\vrule&&\vrule\end{array}\right], E_{i}=\left[\begin{array}[]{cccc}\vrule&\vrule&&\vrule\\ e_{1}^{(i)}&e_{2}^{(i)}&\ldots&e_{k}^{(i)}\\ \vrule&\vrule&&\vrule\end{array}\right] and

H_{N}=\dfrac{1}{N}\displaystyle\sum_{i=1}^{N}E_{i}=\left[\begin{array}[]{cccc}\vrule&\vrule&&\vrule\\ h_{1}&h_{2}&\ldots&h_{k}\\ \vrule&\vrule&&\vrule\end{array}\right]. Since the entries of $F_{N}$ have zero mean, to compute their variance we start by computing

[TABLE]

Hence, taking the expectation we get

[TABLE]

Equivalently,

[TABLE]

where E=\left[\begin{array}[]{cccc}\vrule&\vrule&&\vrule\\ e_{1}&e_{2}&\ldots&e_{k}\\ \vrule&\vrule&&\vrule\end{array}\right] is a $d\times k$ random matrix with i.i.d. Gaussian coefficients and $A_{s,t,u,v}=\mathbb{E}[\langle e_{s},e_{t}\rangle\langle e_{u},e_{v}\rangle]$. By expanding $A_{s,t,u,v}$, we can write

[TABLE]

In order to compute $\mathbb{E}\left[z_{s}z_{t}z_{u}z_{v}\right]$, we distinguish seven cases below:

[TABLE]

In other words,

[TABLE]

We therefore get that the entry of matrix $\Sigma$ defined in (31) corresponding to the covariance between elements $(s,t)$ and $(u,v)$ of $\hat{M}_{N}$ is

[TABLE]

where $G=X^{T}X$.

F.2 Mean Squared Error of $\hat{G}_{N}$

Using Lemma 3.1 for example (Tu et al., 2016, Lem. 5.4), we find

[TABLE]

where $\tilde{G}_{N}=\underset{H\succeq 0:\mathrm{rank}(H)\leq d}{\mathrm{argmin}}\lVert\hat{G}_{N}-H\rVert_{\mathrm{F}}$ is the projection of $\hat{G}_{N}$ on the cone of positive semidefinite matrices of rank at most $d$. Since $\mathrm{rank}(G)=d$, $\tilde{G}_{N}$ satisfies

[TABLE]

By triangle inequality we have

[TABLE]

Hence

[TABLE]

Building on the computations from subsection F.1, we get that

[TABLE]

Finally, (33) and (34) imply that

[TABLE]

This directly implies that, in the case where $\sigma\ll 1$, we have $\mathbb{E}[\rho^{2}(X,\hat{X})]=O\left(\dfrac{\sigma^{2}}{N}\right)$, while in the case where $\sigma\gg 1$, we have $\mathbb{E}[\rho^{2}(X,\hat{X})]=O\left(\dfrac{\sigma^{4}}{N}\right)$.