Nonlinear matrix recovery using optimization on the Grassmann manifold

Florentin Goyens[This author’s work was supported by The Alan Turing Institute under the EPSRC Grant No. EP/N510129/1 and under the Turing Project Scheme.]

Coralia Cartis[This author’s work was supported by The Alan Turing Institute under the Turing Project Scheme.]

Armin Eftekhari

We investigate the problem of recovering a partially observed high-rank matrix whose columns obey a nonlinear structure such as a union of subspaces, an algebraic variety or grouped in clusters. The recovery problem is formulated as the rank minimization of a nonlinear feature map applied to the original matrix, which is then further approximated by a constrained non-convex optimization problem involving the Grassmann manifold.
We propose two sets of algorithms, one arising from Riemannian optimization and the other as an alternating minimization scheme, both of which include first- and second-order variants.
Both sets of algorithms have theoretical guarantees. In particular, for the alternating minimization, we establish global convergence and worst-case complexity bounds. Additionally, using the Kurdyka-Lojasiewicz property, we show that the alternating minimization converges to a unique limit point.
We provide extensive numerical results for the recovery of union of subspaces and clustering under entry sampling and dense Gaussian sampling.  Our methods are competitive with existing approaches and, in particular, high accuracy is achieved in the recovery using Riemannian second-order methods.Keywords: nonlinear matrix recovery, nonconvex optimization, Riemannian optimization, second-order methods.

§ INTRODUCTION

In the matrix recovery problem, one tries to estimate a matrix $M\in \mathbb{R}^{n\times s}$ from partial information.  The low-rank matrix recovery problem deals with instances where the matrix $M$ is low-rank.  This problem has received great attention in the literature, as applications abound in recommender systems and engineering (see <cit.> and the references therein for an overview). It was shown in <cit.> that solving a convex program allows to recover the original matrix $M$ with very high probability, provided enough samples are available. However, solving this convex semi-definite problem for large instances is very costly in time and memory allocation. This has sparked the search for alternative nonconvex formulations of the problem <cit.>. Riemannian optimization methods are used in some of the most efficient algorithm known to date for low-rank matrix completion. These methods solve optimization problems defined on smooth Riemannian manifolds, such as the manifold of fixed-rank matrices <cit.> or the Grassmann manifold <cit.>.

All traditional approaches to matrix completion fail if the matrix $M$ is high-rank. Our work is based on the recent discovery that an adaptation of traditional methods allows to recover specific classes of high-rank matrices <cit.>. This problem is known as nonlinear matrix recovery (or high-rank matrix recovery). Recovering high-rank matrices requires to make assumptions on the structure of $M$. Let $m_1, \dots, m_s$ denote the columns of $M$. When the points $m_i \in \R^n$ belong to a low-dimensional subspace in $\R^n$, low-rank matrix recovery methods can be applied. Nonlinear matrix recovery attempts to recover $M$ when the points $m_i$ are related in a nonlinear way.

Classically, for some integer $ m <ns$, the matrix $M$ satisfies $m$ linear equations of the type $\langle A_i,M\rangle =b_i$ for given matrices $A_i\in  \mathbb{R}^{n\times s} $ and a given vector $b\in \mathbb{R}^m$, where we use the usual inner product $\langle A_i,M\rangle = \trace(A_i^\top M)$. The matrices $A_i$ are assumed to be randomly drawn from a known distribution. One defines the linear operator

\begin{equation}
\mathcal{A}: \R^{n\times s} \longrightarrow \R^m \txt{ where }\A(M)_i = \langle A_i,M\rangle
\label{eq:measurements}
\end{equation}

so as to have the compact notation $\mathcal{A}(M)=b$ for the measurements. When each matrix $A_i$ has exactly one non-zero entry which is equal to $1$, this is known as a matrix completion problem.  The matrix $M$ is then known on a subset $\Omega$ of the complete set of entries $\{1, \dots, n\} \times \{1, \dots, s\}$. Without loss of generality we assume $n\leq s$.

Problem description
Nonlinear matrix recovery methods use features that map the columns of $M$ to a space of higher dimension.
The feature map is defined as

\begin{equation}
\varphi : \mathbb{R}^n \longrightarrow \mathcal{F} : v \longmapsto \varphi (v),
\end{equation}

where $\mathcal{F}$ is a Hilbert space. If $\mathcal{F}$ is finite dimensional, we write $\mathcal{F}= \R^N$ where  $N$ is the dimension of the lifted space, with $N\geq n$.  We obtain the feature matrix $\Phi(M)$ by applying $\varphi$ to each column of $M$,

\begin{equation}
\Phi(M) = \begin{bmatrix}
\varphi(m_1) & \dots & \varphi(m_s)
\end{bmatrix}\in \R^{N\times s}.
\label{eq:features_matrix}
\end{equation}

The map $\varphi$ is chosen using a priori knowledge of the data so that the features of the data points $\varphi(m_i)$ for $i=1, \dots, s$, all belong to the same subspace in $\R^N$. The nonlinear structure in $M$ will cause a rank deficiency in the feature matrix $\Phi(M)$, even though $M$ may be full-rank.  This is illustrated in Figure <ref>.

[scale = 0.5]
[hide axis,colormap/violet]
domain y=-1.3:1.3,

[fill=black] (1,2.5) circle (.1);
[fill=black] (5,3.6) circle (.1);
[fill=black] (5.4,2.9) circle (.1);
[fill=black] (6,3.2) circle (.1);
[fill=black] (6,3.2) circle (.1);
[fill=black] (4.2,2) circle (.1);
[fill=black] (4.5,2.5) circle (.1);
[fill=black] (4,3) circle (.1);
[fill=black] (3,3) circle (.1);
[fill=black] (2,2.3) circle (.1);
[fill=black] (1.7,2.8) circle (.1);
[fill=black] (2.5,3.5) circle (.1);
[right] at (6.5,4) $\mathbb{R}^n$;

[scale = 0.4]

[above] at (-6.6,0.5) $\varphi$;
[->,thick] (-7.6,-0.3) to [out = 0, in = 180] (-5.6,-0.3);
[left color=blue!40,right color=blue!30] (4,0) – (3.6,-2.3) – (-4.6,-0.3) – (-2.7,1.5) – (4,0);
[above] at (4,1) $ \mathcal{F}$;
[fill=black] (3,-1) circle (.1);
[fill=black] (2.5,-0.8) circle (.1);
[fill=black] (2,-1.2) circle (.1);
[fill=black] (1.5,0.3) circle (.1);
[fill=black] (1,-0.5) circle (.1);
[fill=black] (3,-1) circle (.1);

[fill=black] (-3,0.3) circle (.1);
[fill=black] (-2.5,1.2) circle (.1);
[fill=black] (-2,0.7) circle (.1);
[fill=black] (-1.5,0.7) circle (.1);
[fill=black] (-1,-0.5) circle (.1);
[fill=black] (-3.3,-0.5) circle (.1);

The feature map $\varphi$ is chosen to exploit the nonlinear structure.

If the features are infinite dimensional or that $N$ is very large, the feature map should be represented using a kernel, which is known as the kernel trick. The set $\mathcal{F}$ is then called a reproducing kernel Hilbert space.  The kernel map represents the inner product between elements in the Hilbert space of features,

\begin{equation}
k: \R^n \times \R^n \longrightarrow \R:	k(x,y) = \langle \varphi(x),\varphi(y)\rangle_\mathcal{F}.
\label{eq:def_kernel}
\end{equation}

This allows to define the kernel matrix of the data $K(M,M)\in \R^{s \times s}$, with $K_{ij}(M,M) = k(m_i,m_j)$. Throughout, we assume that $r = \rank(\Phi(M))$ is known and smaller than $\min(N,s)$. When $\mathcal{F}= \R^N$, we note that $K(M,M) = \Phi(M)^\top \Phi(M) $ and therefore $\rank \big(K(M,M)\big) = \rank \big(\Phi(M)\big)$. We use the term lifting to denote a kernel or a feature map instinctively.
In <cit.>, Ongie et al. use the monomial kernel for the completion of matrices whose columns belong to an algebraic variety (a set defined by a finite number of polynomial equations). This can notably be applied to a union of subspaces. In <cit.>, Fan et al. use the monomial kernel and the Gaussian kernel on image inpainting problems. In Section <ref>, we detail why polynomial and Gaussian kernels may be used to model data which belongs to, respectively, an algebraic variety or a set of clusters.

*Problem formulation
For an appropriately chosen feature map, the nonlinear matrix recovery problem can be formulated as the rank minimization of the features matrix under the measurements constraint

\begin{equation}
\left\lbrace
\begin{aligned}
& \underset{X\in \R^{n\times s}}{\min}
& & \rank(\Phi(X)) \\
&&&  \Acal(X) = b.
\end{aligned}
\right.
\label{eq:min_rank_phi}
\end{equation}

This seeks to find the matrix which fits the observations using a minimum number of independent features. As is the case for low-rank matrix recovery, minimizing the rank directly is NP-hard and should be avoided <cit.>. It is necessary to find a suitable approximation to the rank function.

Related work
In essence, <cit.> and <cit.> apply different minimization algorithms to the Schatten p-norm of the features which is defined by

\begin{equation}
\left(\sum_{i=1}^{\min(N,s)} \sigma_i (\Phi(X))^p \right)^{1/p}~~\txt{ for } 0<  p \leq 1.
\end{equation}

When $p=1$, the sum of the singular values is the nuclear norm. Both <cit.> and <cit.> use a kernel representation of the features, so that the features are never computed explicitely. In <cit.>, the authors introduce the algorithm NLMC, which applies a quasi-Newton method to minimize the Schatten p-norm. The Schatten p-norm for $0<p\leq 1$ is nonsmooth. This has the benefit of encouraging sparsity in the singular values, but it might prevent fast convergence near a minimizer. The algorithm VMC was introduced in <cit.>. It applies an iterative reweighted least squares (IRLS) which was originally proposed in <cit.> and <cit.> for low-rank matrix recovery and rank minimization. The IRLS framework has the advantage that it generalizes seemlessly to the kernel version.
In <cit.> a truncated version of the Schatten norm is proposed, where only the smallest singular values are minimized,

\begin{equation}
\left(\sum_{i=r+1}^{\min(N,s)} \sigma_i (\Phi(X))^p \right)^{1/p}~~\txt{ for } 0<  p \leq 1.
\end{equation}

where $r = \rank(\Phi(M))$. They use the kernel trick and propose an algorithm which alternates between truncated singular value decompositions and a step of the Adam method with an additional tuning of the stepsizes.

In <cit.>, the authors propose an extension to handle outliers in the data. This is achieved by introducing a sparse matrix in the model, which absorbs the outliers.
In <cit.>, the authors build a tensor representation of the data and apply known matrix completion techniques in the tensor space. Their algorihm, LADMC, is a simple and efficient approach for which they are able to show that the sampling requirements nearly match the information theoretic lower bounds for recovery under a union of subspace model. This is remarkable as the sampling pattern in the tensor space is not random, and low-rank recovery results do not apply directly. Note that the approach in <cit.> is only applicable to matrix completion problems, not matrix sensing.

In <cit.>, a new algorithm KFMC  (kernelized factorization matrix completion) is proposed, which lends itself to online completion. In this setting, the columns of the matrix $M$ are accessible as a stream and the matrix $M$ is never stored in its entirety. They also develop a variant algorithm to deal with out of samples extensions. That is, how to complete a new column without recomputing the model. The offline formulation applies the kernel trick to

\begin{equation}
\left\lbrace
\begin{aligned}
& \min_{X,D,Z}
& &\norm{\Phi(X) - \Phi(D)Z}^2_F + \alpha \norm{\Phi(D)}^2_F + \beta \norm{Z}^2_F \\
& %\text{subject to}
& & X_{ij} = M_{ij}, \; (i,j) \in \Omega, \\
&&&D \in \R^{n\times r}, Z\in \R^{r\times s}, X\in \R^{n\times s}.
\end{aligned}
\right.
\end{equation}

The variable $D\in \R^{n\times r}$ aims to find $r$ points in $\R^n$ so that their features will form a basis for $\Phi(M)$ in the Hilbert space. The last two terms in the objective are added as regularizers to improve the practical performances of the algorithm, as is often done in low-rank matrix completion <cit.>.
In the method DMF, proposed in <cit.>, the lifting is replaced by a deep-structure neural network who is trained to minimize the reconstruction error for the observable entries of $M$. This work showcases the applicability and performance of nonlinear matrix completion with numerous examples including image inpainting and collaborative filtering problems. For data drawn from multiple subspaces, <cit.> proposes a sparse factorization where each subspace is represented in a low rank decomposition. They solve this model with an algorithm in the spirit of PALM <cit.>. They are able to show sampling bounds to guarantee recovery.

Contribution and outline of the paper

In Section <ref> we describe the approach taken to recover high-rank matrices. It consists in using a feature map (or kernel) that exploits the nonlinear structure present in the matrix.
This is applied to data which follows algebraic variety models or grouped in clusters. For these, we respectively use the monomial kernel and the Gaussian kernel. We demonstrate that the Gaussian kernel can be used to perform clustering with missing data, which expands the use cases of nonlinear matrix recovery.

In Section <ref>,  we propose to use a new formulation for nonlinear matrix recovery. We use the feature map to write the recovery problem as a constrained nonconvex optimization problem on the Grassmann manifold. This extends the residual proposed in <cit.> in the context of low-rank matrix completion to the nonlinear case.

We propose to use Riemannian optimization methods to solve the recovery problem, which is new in the context of nonlinear matrix recovery. Riemannian optimization, as described in Section <ref>, provides a framework to design algorithms for problems with smooth constraints. This allows to seamlessly choose between standardized first- and second-order methods. The use of second-order methods allows to recover high-rank matrices up to high accuracy if desired.

Section <ref> presents an alternating minimization algorithm to solve the recovery problem. First- and second-order variants of the alternating minimization are discussed. We prove global convergence of the algorithm to first-order stationary points in Section <ref> and give a global complexity rate to achieve an arbitrary accuracy on the gradient norm from an arbitrary initial guess. In Section <ref>, we also show convergence of the sequence of iterates to a unique limit point using the Kurdyka-Lojasiewicz property.
Our alternating minimization method is a similar approach to the method proposed in <cit.>. We provide extensive convergence analysis, which was not done in <cit.>.

Section <ref> summarizes the applications and algorithms covered in this paper with a framework to solve nonlinear matrix completion.

We conclude with an extensive set of numerical experiments that compare the performances of the optimization approaches and the quality of the solutions that can be obtained (Section <ref>). We discuss the influence of the complexity of the data and the role of model parameters on the recovery. Moreover, we showcase that our approach is very efficient at clustering data with missing information.

Throughout the paper we use a notation consistent with <cit.> for the derivative of a function $f$ defined on a Riemannian manifold.  The unconstrained gradient of a function $f$ is written $\nabla f$, when the domain of $f$ is extended to an embedding Euclidean space. Conversely, we use $\grad f$ for the Riemannian gradient of $f$ defined over a Riemannian manifold. For matrices $A,B \in \Rns$, $\inner{A}{B} = \trace(A\transpose B)$ is the canonical inner product, $\range(A)$ is the column space of $A$, $\Null(A)$ is the null space of $A$. The identity matrix of size $n$ is denoted by $\I_n$ and $\Id$ is the identity operator.

§ THE FEATURE MAP

As mentioned, our approach uses a lifting to a space of features, in the spirit of <cit.><cit.>. In this section, we give more details and examples on the role of this lifting.

The lifting $\varphi$ is chosen such that $\Phi(M)$ is low rank. In addition, $\Phi(X)$ should be high rank if $X$ does not exhibit the same geometrical structure as $M$.

The goal is to find a lifting that reveals the nonlinear relation between the points $m_i$, the columns of $M$. In <cit.> the authors use the polynomial features for data sets living on algebraic varieties.

[(Algebraic varieties <cit.>)]

Let $\R[x]$ be the set of real valued polynomials over $\R^n$. A real (affine) algebraic variety is defined as the zero set of a system of polynomials $P \subset \R[x]$:

\begin{equation}
V(P) = \lbrace x\in \R^{n}: p(x) = 0 \txt{ for all } p \in P\rbrace.
\end{equation}

We say that the matrix $M $ follows an algebraic variety model if every column of $M$ belongs to the same algebraic variety.

\begin{equation}
N(n,d) = \begin{pmatrix}
\end{pmatrix},
\label{eq:N_features}
\end{equation}

be the number of monomials of degree less or equal to $d$ that can be formed with $n$ variables.
The monomial features $\varphi_d$ for some degree $d$ are defined as

\begin{equation}
%\tag{Monomial features}
\varphi_d:\R^n \longrightarrow \R^{N(n,d)} : \varphi_d(x) = \begin{pmatrix}
\vdots \\
\end{pmatrix}
\label{eq:monomial-features}
\end{equation}

where $\bm{\alpha}^i$ for $i = 1, 2,\dots,N(n,d)$ is a multi-index of non-negative integers $(\alpha^i_1,\alpha^i_2,\dots, \alpha^i_n)$ so that
$x^{\bm \alpha^i} := x_1^{\alpha^i_1} x_2^{\alpha^i_2}\dots x_n^{\alpha^i_n}$ and  $\alpha^i_1 + \alpha^i_2 + \dots + \alpha^i_n \leq d$.
The dimension of the lifted space $N(n,d)$ increases exponentially in $n$ and $d$. Therefore, a kernel implementation is usually used in practice for moderate and large dimensions, or more precisely, whenever $s\leq N(n,d)$. The monomial kernel of degree $d$ is defined for any $X,Y\in \R^{n\times s}$ as

\begin{equation}
%\tag{Monomial kernel}
K_d(X,Y) = (X^\top Y + c  \mathbf{1}_{s\times s} )^{ \odot d},
\label{eq:mono_kernel}
\end{equation}

where the value $c\in \R$ is a parameter of the kernel, $\mathbf{1}_{s\times s} $ is a square matrix of size $s$ full of ones and $\odot $ is an entry-wise exponent. If the equations describing the variety are known to be homogeneous, one can set $c=0$. Note that the monomial kernel in (<ref>) is not exactly the kernel associated with the monomial features in (<ref>). Instead, $K_d(X,X) = \tilde{\Phi}_d(X)^\top \tilde{\Phi}_d(X)$ for a similar map of monomials $\tilde{\Phi}_d$ that has different coefficients than $\Phi_d$. Therefore, $\Phi_d(X)$ and $K_d(X,X)$ have the same rank.

In <cit.>, the authors argue why using this lifting is appropriate when the points $m_i$ belong to an algebraic variety. Suppose the variety $V(P)\subset \R^n$ is defined by the set of polynomials $P=\lbrace p_1,\dots,p_k\rbrace$ where each $p_i$ is at most of degree $d$. Then the columns of $X$ belong to the variety $V(P)$ if and only if there exists $C\in \R^{N\times k}$ such that $\Phi_d(X)^\top C=0$, where the columns of $C$ define the coefficients of the polynomials $p_i$ in the monomial basis. This implies that $\rank(\Phi_d(X)) \leq \min(N-k,s)$. This justifies that $\Phi_d(X)$ is rank deficient when there are sufficiently many data points such that $s\geq N-k$, and $X$ follows an algebraic variety model. The second case study below presents a union of subspaces as a particular type of algebraic variety.

[(Union of subspaces)]

Given two affine subspaces $\mathcal{S}_1$, $\mathcal{S}_2\subset \R^n$ of dimension $r_1$ and $r_2$ respectively,  we can write $\mathcal{S}_1 =\lbrace x:q_i(x)=0 \txt{ for } i=1,\dots,n-r_1\rbrace$
and $\mathcal{S}_2 =\lbrace x:p_j(x)=0 \txt{ for } j=1,\dots, n-r_2\rbrace$
where the $q_i$ and $p_j$ are affine functions. The union $ \mathcal{S}_1 \cup \mathcal{S}_2$ can be expressed as the set where all possible products $q_i(x)p_j(x)$ vanish.  Therefore, $ \mathcal{S}_1 \cup \mathcal{S}_2$ is the solution of a system of $(n-r_1)(n-r_2)$ quadratic polynomial equations. Similarly, a union of $k$ affine subspaces of dimensions $r_1 , \dots, r_k$ is a variety described by a system of $\Pi_{i=1}^k (n-r_i)$ polynomial equations of degree $k$.

If the columns of a matrix $X\in \R^{n\times s}$ belong to a union of $p$ affine subspaces of dimension at most $\tilde r$, then for any degree $d\geq 1$, the matrix $\Phi_d(X) \in \R^{N(n,d)\times s}$ of monomial features, with $N(n,d)$ the dimension of the features space defined in equation (<ref>), satisfies

\begin{equation}
\rank \Phi_d(X) \leq p\begin{pmatrix}
\tilde r+d\\
\end{pmatrix}.
\end{equation}

In practice, choosing the degree $d$ of the monomial kernel is a tricky task. In Section <ref>, we discuss the practical choice of this degree and how it impacts the rank of the feature matrix and the possibility to recover $M$. Previous works using the monomials kernel to recover high-rank matrices all restricted themselves to degrees two or three <cit.>. Using a polynomial lift of large degree would seem helpful to capture all the nonlinearity in some data sets. Unfortunately, increasing the degree will grow the dimensions of the optimization problem exponentially. Indeed, the lifted dimension $N(n,d)$ blows up with $d$ for even moderate values of $n$ and the number of data points required is at least $N(n,d)-k$ where $k$ is the number of polynomial equations that define the variety.

We now define the Gaussian kernel which will be used when the columns of the matrix $M$ are grouped in several clusters.

For $X, Y\in \R^{n\times s}$, the entry $(i,j)$ of the Gaussian kernel $K^G: \R^{n\times s} \times \R^{n\times s} \longrightarrow \R^{s\times s}$ is defined as

\begin{equation}
K_{ij}^G(X,Y) = e^{ -\dfrac{\|x_i - y_j\|^2_2}{2\sigma^2}},
\label{eq:gaussian_kernel}
\end{equation}

where $\sigma>0$ is the width of the kernel. The Gaussian kernel acts as a proximity measure. For two columns of $X$, labelled $x_i$ and $x_j$, we observe that $x_i$ being close to $x_j$ gives $K^G_{ij}(X,X) \approx 1$ and if $x_i$ is far from $x_j$ then $K^G_{ij}(X,X) \approx 0$. Therefore, the rank of the Gaussian kernel approximately coincides with the number of clusters in $X$. More precisly, one can show that the singular values whose index exceeds the number of clusters decay rapidly <cit.>. The value of $\sigma$ should be chosen appropriately depending on the size of the clusters.

In Figure <ref>, we present a small data set of 100 data points divided in two clusters in $\R^2$ with the singular values of the Gaussian kernel (in log-scale). We see that the two largest singular values are much greater than the third one, and that the following singular values decrease at an approximately exponential rate. Therefore, the Gaussian kernel is near a low rank matrix for clustered data, which will allow us to complete such data sets from partial measurements.
In <cit.>, the Gaussian kernel was also used effectively on image inpainting and denoising problems.

Clustered data and the singular values of the Gaussian kernel in log-scale.

§ NONLINEAR MATRIX RECOVERY AS AN OPTIMIZATION PROBLEM

Noiseless measurements case
Considering a noiseless measurements case, we would like to minimize the rank of the lifted matrix, as in Equation (<ref>). This is unfortunately known to be intractable, even in the case where $M$ is low-rank <cit.>. We have to resort to approximations of this problem. The second difficulty is the nonlinearity of $\Phi$.

As a nonconvex approximation to (<ref>), we consider the formulation in <cit.> for low-rank problems, and extend it to the nonlinear case. Assuming that $\Phi(M)$ has rank $r$ leads to the following formulation

\begin{equation}
\left\lbrace
\begin{aligned}
& \underset{X,\U}{\min}
& & f(X,\U) := \fronorm{ \Phi(X) - P_{\U}\Phi(X)}^2 \\
& %\text{subject to}
& & \U\in \Grass(N,r) \\
&&& \mathcal{A}(X)=b,
\end{aligned}
\right.
\label{eq:p1}
\end{equation}

where  $\Grass(N,r)$ is the Grassmann manifold, the set of all subspaces of dimension $r$ in $\R^{N}$, $P_{\U}$ is the orthogonal projection on the subspace $\U$ and $\fronorm{.}$ denotes the Frobenius norm. The linear measurements, $\A(X)=b$, are defined in equation (<ref>).  Given $U\in \R^{N\times r}$ such that $\range(U) = \U$ and $U^\top U=I_r$, the projection is given by $\Prm_\U = UU^\top$. In (<ref>) ,the objective function is expected to be nonconvex but smooth for practical choices of $\varphi$, such as case study <ref>, <ref>.
If the variable $\U$ is additionally constrained to be the range of the $r$ leading singular vectors of $\Phi(X)$, the cost function becomes
$\sum_{i = r+1}^{\min(N,s)} \sigma_i(\Phi(X))^2$. The advantage of the formulation in (<ref>) is that is it straightforward to express it as a finite sum of $s$ terms over all the data points. This allows to use stochastic sub-sampling algorithms that scale better to matrices with many columns (large $s$).

The advantage of using the Grassmann manifold, which is a quotient space, instead of the Stiefel manifold of orthogonal matrices $\St(N,r) := \lbrace U\in \R^{N\times r}: U^\top U = \I_r\rbrace$ is that due to the invariance of the cost function with respect to the matrix that represents the subspace $\U$, local optimizers cannot possibly be isolated in a formulation over $\St(N,r)$. Therefore, the fast local convergence rates of some second-order algorithms might not apply on $\St(N,r)$, while they would apply on the quotient manifold.

Consider $U_\perp$ a basis of $\U_\perp$, the orthogonal complement of $\U$. The variable $\U_\perp\in \Grass(N,N-r)$ has a nice interpretation since, if $f(X,\U) =0$, then $U_\perp$ spans $\Null(\Phi(X)^\top)$. In the case of algebraic varieties (case study <ref>), $U_\perp$ gives the coefficients of the polynomials defining the variety in the basis given by $\Phi$. Recovering the equations of the variety is of interest in some applications <cit.>.

Noisy measurements case
When the measurements are known to be noisy, i.e.  $\A(M)=b + \eta$ for some noise $\eta\in \R^m$, the measurement constraint can be lifted into the cost function as a penalty.  This gives

\begin{equation}
\left\lbrace
\begin{aligned}
& \underset{X,\U}{\min}
& &  f_\lambda(X,\U):=  \| \Phi(X) - P_{\U}\Phi(X)\|^2_F  + \lambda \norm{ \mathcal{A}(X) -b}^2_2\\
& %\text{subject to}
& & \U\in \Grass(N,r), \\
\end{aligned}
\right.
\label{eq:p1_noise}
\end{equation}

where the parameter $\lambda>0$ has to be adjusted.  This allows to satisfy the noisy measurements approximately.

§.§ Kernel representation of the features

When working with a kernel instead of a feature map, as shown in (<ref>), we want to find a cost function equivalent to that of (<ref>) which uses the kernel instead of the feature map. We find that they are related in the following way.

Given a feature map $\Phi$ and the associated kernel $K : X\times X \longmapsto \Phi(X)^\top \Phi(X)$, for $\W \in \Grass(s,r)$ we have

\begin{equation}
\fronorm{\Phi (X)^\top - P_{\W}\Phi (X)^\top}^2 = \trace\Big( K(X,X) - P_{\W}K(X,X) \Big).
\label{eq:trace-kernel}
\end{equation}

We write $P_{\W^{ \bot}} =  \I_{s\times s} - P_{\W}$ and find

\begin{align}
\nonumber
\trace\Big( P_{\W^{ \bot}}K(X,X) \Big) &=  \trace\Big( P_{\W^{ \bot}}\Phi(X)\transpose \Phi(X) \Big)\\\nonumber
& = \trace\Big(\Phi (X)  P_{\W^{ \bot}}\Phi(X)^\top \Big) \\
&=   \trace\Big(\Phi (X)  P^\top_{\W^{ \bot}} P_{\W^{ \bot}}\Phi(X)\transpose \Big)\\\nonumber
&  =  \trace\Big((P_{\W^{ \bot}}\Phi(X)^\top)^\top   P_{\W^{ \bot}}\Phi(X)\transpose \Big)\\
%&=   \langle  P_{\W^{ \bot}}\Phi(X)\transpose,  P_{\W^{ \bot}}\Phi(X)\transpose\rangle
& = \fronorm{P_{\W^{ \bot}}\Phi(X)^\top}^2.\nonumber
\end{align}

Using the kernel formula (<ref>) corresponds to finding a subspace $\W$ of dimension $r$ in the row space of $\Phi(X)$.
When a kernel is used to lift the data, the following optimization problem should be solved,

\begin{equation}
\left\lbrace
\begin{aligned}
& \underset{\W,X}{\min}
& & f(X,\W) := \trace\big(P_{\W_\perp} K(X,X)\big) \\
& %\text{subject to}
& & \W\in \Grass(s,r) \\
&&& \mathcal{A}(X)=b.
\end{aligned}
\right.
\label{eq:p1kernel}
\end{equation}

Replacing the features $\Phi$ by the corresponding kernel becomes beneficial when the dimension $N$ of the features is larger than the number of points $s$ and when a convenient formula is available to compute the kernel and its derivatives. For example, in the case of clusters (case study <ref>), the features exist implicitly in an infinite dimensional space and we use the Gaussian kernel to represent them.

In the upcoming sections, we usually describe the algorithms and their properties using the notation of problem (<ref>) with a feature map $\Phi$ and cost function $f(X,\U)$. Unless specified otherwise, the developments also apply to problem (<ref>) and the use of a kernel.

§ RIEMANNIAN OPTIMIZATION ALGORITHMS

Riemannian optimization methods provide an elegant and efficient way to solve optimization problems with smooth nonlinear constraints. The field of Riemannian optimization has  rapidly developed over the past two decades. In particular,
Riemannian optimization methods have proved very efficient in low-rank matrix completion <cit.>. In this section, we investiate the use of Riemannian optimization methods to solve (<ref>) or (<ref>). For an overview of optimization algorithms on Riemannian manifolds, see <cit.>. In order to formally express (<ref>) as a Riemannian optimization problem, we define a notation for the affine subspace that represents the measurements on the matrix $M$,

\begin{equation}
\LAb = \{X\in \R^{n\times s}: \A(X)=b\}.
\label{eq:Lab}
\end{equation}

We form the product manifold
\begin{equation}
\M = \LAb \times \Grass(N,r),
\label{eq:M}
\end{equation}

so that Problem (<ref>) can be viewed as  the unconstrained minimisation of a smooth cost function defined on the manifold $\M$,

\begin{equation}
\left\lbrace
\begin{aligned}
& \underset{(X,\U)}{\min}
&  &  \fronorm{ \Phi(X) - P_{\U}\Phi(X)}^2 \\
& & &(X,\U) \in \M.
\end{aligned}
\right.
\label{eq:p1riem}
\end{equation}

We introduce the notation $z = (X,\U) \in \M$ to denote both variables that appear in the optimization problem. Riemannian optimization algorithms are feasible methods that iteratively exploit the local geometry of the feasible set. Analogously to unconstrained optimization in Euclidean spaces, each iteration of a Riemannian optimization algorithm uses derivatives to build a model that locally approximates the cost function. This model is then fully or approximately minimized. Most commonly, the model uses the gradient and possibly the Hessian or an approximation of it. This yields respectively a first- or second-order method. Riemannian optimization follows these principles, only the model is defined on a local linearisation of the manifold, namely, the tangent space. The Riemannian gradient, written $\grad f$ is a vector belonging to the tangent space of the manifold. The Riemannian Hessian, written $\Hess f$ is a symmetric operator on that tangent space.
In Appendix <ref>, we show how to compute the Euclidean gradient and Hessian for the cost function of (<ref>) in the case of the kernels presented in Section <ref> (monomial kernel, Gaussian kernel). Then, to find their Riemannian counterparts,
$\nabla f$ and $\nabla^2 f$ are projected onto the tangent space of $\M$ using the tools defined later on in this section. When exploring the tangent space, it is necessary to have a tool that allows to travel on the manifold in a direction prescribed by a tangent vector. This operation is called a retraction <cit.>.

A retraction on a manifold $\M$ is a smooth mapping $Retr$ from the tangent bundle $\Trm\M$ to $\M$ with the following properties. Let $\Retr_z: \Trm_z\M \longrightarrow \M$ denote the restriction of $\Retr$ to $\Trm_z \M$.

(i) $\Retr_z(0_z) = z$, where $0_z$ is the zero vector in $\Trm_z\M$;

(ii) the differential of $\Retr_z$ at $0_z$, $\D \Retr_z(0_z)$, is the identity map.

The retraction curves $t\longmapsto \Retr_z(t\eta)$ agree up to first order with geodesics passing through $z$ with velocity $\eta$, around $t=0$.

Let us detail the tools necessary to use Riemannian optimization methods on the two manifolds that compose our search space $\M$. Note that the cartesian product of two Riemannian manifolds is a Riemannian manifold. The geometry of $\LAb$ is rather trivial because the manifold is affine. It must nonetheless be implemented so we will describe how to handle this constraint in a Riemannian way. More generally, this gives a straightforward way do deal with affine equality constraints.

Measurement subspace
At any point $X\in \LAb$, the tangent space is the null space of $\A$,

\begin{equation}
\Trm_X \LAb = \Null(\A) = \{ \Delta \in \R^{n\times s}: \A(\Delta)=0\}.
\end{equation}

Since the tangent space does not depend on $X$, we write $\Trm \LAb$. This tangent space inherits an inner product from the embedding space $\R^{n\times s}$,

\begin{equation}
\langle \Delta_1 , \Delta_2 \rangle = \trace(\Delta_1^\top \Delta_2)~~\txt{for all }\Delta_1,\Delta_2 \in \Trm\LAb.
\end{equation}

The Riemannian gradient is the orthogonal projection of the Euclidean gradient onto  the tangent space $\Trm\LAb$.

From the fundamental theorem of algebra, $\Null(\A) = \range(\A^\top)^\perp$. Therefore we can express $\Prm_{\Null(\A)} = \Id - \Prm_{\range(\A^\top)}$.

The application $\A $ is represented by a flat matrix $A \in \R^{m \times ns}$ such that $\A(X)= AX(:)\in \R^m$, where $X(:)$ is a vector of length $ns$ made of the columns of $X$ taken from left to right and stacked on top of each other. Visually this gives,

\begin{equation}
\A(X)=\begin{pmatrix}
\langle A_1,X\rangle\\
\langle A_2,X\rangle\\
\vdots \\
\langle A_m,X\rangle\\
\end{pmatrix} = \begin{pmatrix}
\vdots\\
\end{pmatrix}
X(:) =:  AX(:).
\label{eq:Asensing}
\end{equation}

The tall matrix $A^\top$ represents the linear application $\Acal^\top$. Viewing $A^\top$ as the matrix of an overdetermined system of linear equations convinces us that $\Prm_{\range(\A^\top)} = A^\top (AA^\top)^{-1}A$. Equivalently, if $Q \in \R^{m\times m}$ is an orthogonal basis for $\range(\A^\top)$ (which can be obtained by a reduced QR factorization of $A^\top$), we can apply $\Prm_{\range(\A^\top)} = Q Q^\top.$  The projection is given by

\begin{equation}
\Prm_{\Null(\A)} = \I_{ns} - A^\top(AA^\top)^{-1}A =  \I_{ns} - Q Q^\top .
\label{eq:PnullA}
\end{equation}

\begin{equation}
\Prm_{\Trm \LAb}: \R^{n\times s} \longrightarrow  \Trm \LAb: \Prm_{\Trm \LAb}(\Delta)  = (\I_{ns} - Q Q^\top) \Delta
\label{eq:PTLAb}
\end{equation}

where $Q \in \R^{m\times m}$ is an orthogonal basis for $\range(\A^\top)$.
In the case of matrix completion, the operator $\A$ selects the known entries of $M$. The description of the feasible subspace is simplified. We write $\calL_{\Omega,b} = \{ X: X_{ij} = M_{ij}, ij\in \Omega\}$ to make explicit that the measurements correspond to matrix completion. The tangent space is  $\Trm \mathrm{L}_{\Omega,b} = \{ \Delta: \Delta_{ij} =0 \txt{ for } ij \in \Omega\}$ and the projection onto $\Trm \mathrm{L}_{\Omega,b}$ simply amounts to setting the entries outside of $\Omega$ to zero,

\[
\Prm_{\Trm \calL_{\Omega,b}} (\Delta) = \left\lbrace
\begin{aligned}
&0 &~\txt{for } ij \in  \Omega\\
&\Delta_{ij} &~\txt{for } ij \notin \Omega.
\end{aligned} \right.
\]

The natural retraction on $\LAb$ for $X\in \LAb$ and $\Delta \in \Trm \LAb$ is given by

\begin{equation}
\Retr_X: \Trm \LAb \longrightarrow  \LAb: \Retr_X(\Delta) = X+\Delta,
\end{equation}

because the manifold is flat. These tools are also needed for the Grassmann manifold and we follow the representation given in <cit.>.

Grassmann  manifold
The Grassmann manifold, written $\Grass(N,r)$, is the set of all linear subspaces of dimension $r$ in $\R^N$. A point $\U$ on the Grassmann is represented by a full-rank matrix $U\in \R^{N\times r}$ such that $\range(U) = \U$. For any orthogonal matrix $Y\in \R^{r\times r}$, the matrix $UY$ is also a valid representation of $\U$ since $\range(UY) = \range(U)$. Denote by $\St(N,r) = \{ U \in \R^{N\times r}: U^\top U = I_r\}$ the Stiefel manifold, the set of matrices with orthonormal columns and the orthogonal group $O(r) = \{ Y\in \R^{r\times r}: Y^\top Y = I_r\}$. The orthogonal group introduces an equivalence relation on the Stiefel manifold, where any two matrices with the same column space are equivalent.
In this regard, the Grassmann is a quotient space,

\begin{equation}
\Grass(N,r) = \St(N,r)/\mathrm{O}(r).
\end{equation}

As a manifold, the Grassmann admits a tangent space at every point $\U$. A tangent vector $\mathcal{H} \in \Trm_\Ucal \Grass(N,r)$ is represented by a matrix $H\in \R^{N\times r}$. Once the representation $U$ has been chosen for $\U$, the tangent space $\Trm_\Ucal \Grass(N,r)$ can be represented by its horizontal lift

\begin{equation}
\Trm_U \Grass(N,r) = \{ H \in \R^{N \times r  }: U^\top H =0\} = \range(U_\perp)
\end{equation}

and the projection onto it is given by

\begin{equation}
\mathrm{P}_{\Trm_U \Grass(N,r)}(H) = \Prm_{U_\perp} H = U_\perp U_\perp^\top H,
\label{eq:PTGr}
\end{equation}

where $U_\perp \in \St(N,N-r)$ is any orthonormal matrix whose column space is the orthogonal complement of $\U$. The representation of the tangent space can be equipped with the usual inner product,

\begin{equation}
\langle H_1, H_2 \rangle_U = \trace(H_1^\top H_2),~~~ \forall H_1, H_2 \in \Trm_U  \Grass(N,r).
\end{equation}

From this Riemannian metric, one can define a norm on the tangent space of the Grassmann, which we label

\begin{equation}
\norm{H} := \sqrt{\langle H, H \rangle_U},~~~ \forall H \in \Trm_U  \Grass(N,r).
\end{equation}

Let us call qf, the mapping that sends a matrix to the $Q$ factor of its (reduced) $QR$ decomposition with $Q\in \St(N,r)$ and $R$ upper triangular with positive diagonal entries. A popular retraction on $\Grass(N,r)$ for $\Ucal\in \Grass(N,r)$ and $\Delta \in \Trm_\Ucal \Grass(N,r)$ is given by

\begin{equation}
\Retr_\U : \Trm_\U \Grass(N,r) \longrightarrow \Grass(N,r); 	\Retr_\Ucal(H) = \range\left(\txt{qf}(U + H)\right).
\end{equation}

We are now in a position to present the Riemannian trust-region algorithm <cit.>. This is an extension of the classical trust-region methods <cit.> to Riemannian manifolds.

Riemannian trust-region (RTR)
At each iterate, RTR builds a local model of the function.  The method sequentially minimizes this model under a ball constraint that prevents undesirably large steps where the model does not accurately represent the function. The trust-region subproblem takes the following form around $z_k \in \M$

\begin{equation}
\min_{\substack{\eta\in \Trm_{z_k}\M\\
\norm{\eta}_{z_k} \leq \Delta_k}} \Hat m_{z_k}(\eta):= f(z_k) + \langle \eta , \grad f(z_k)\rangle_{z_k} + \dfrac{1}{2}\langle \eta, H_k [\eta]\rangle_{z_k},
\label{eq:TRsub}
\end{equation}

where $H_k: \Trm_{z_k}\M \longrightarrow \Trm_{z_k}\M$ is a symmetric operator on $\Trm_{z_k}\M$, $\Delta_k$ is the trust-region radius and the model $\hat m_{z_k}: \Trm_{z_k} \M \longrightarrow \R$ is a quadratic approximation of the pullback $\hat f_{z_k}  = f\circ \Retr_{z_k} $, defined on the tangent space at $z_k\in \M$.

First-order Riemannian trust-region
When the Hessian of the cost function is not available, or to obtain a cheap first order trust region method, one should use the identity in the quadratic model, $H_k = \Id$. This method will be very similar to gradient descent. But the trust-region is used to ensure global convergence, as opposed to a line-search.

Second-order Riemannian trust-region
When the true Hessian of the cost function is available, the classical second-order trust-region method is obtained with $H_k = \Hess f(z_k)$. It is also possible to use an approximation of the true Hessian for $H_k$.

Riemannian trust-region (RTR) <cit.>
[1]
Given:  $z_0 \in \M$ and $0< \Delta_0<\bar{\Delta}$, $\varepsilon_g>0$, $\varepsilon_H>0$ and $0<\rho'< 1/4$
Init:  $k=0$
$\norm{\grad f(z_k)}>\varepsilon_g $
Obtain $\eta_k\in \Trm_{z_k}\M$ satisfying assu:1st-order-decrease-rtr
$\varepsilon_H <\infty$
$\lambdamin(H_k) <-\varepsilon_H$
Obtain $\eta_k\in \Trm_{z_k}\M$ satisfying assu:2nd-order-decrease-rtr
return  $z_k$
return $z_k$
$z_k^+ = \Retr_{z_k}(\eta_k)$
$\rho = f(z_k) - f(z_k^+) / (\hat m_{z_k}(0) - \hat m_{z_k}(\eta_k))$
$\rho <1/4$
$\Delta_{k+1} = \Delta_k/4$
$\rho >3/4$  and  $\norm{\eta_k}= \Delta_k$
$\Delta_{k+1} = \min(2\Delta_k, \bar{\Delta})$
$\Delta_{k+1} = \Delta_k $
$\rho >\rho'$
$z_{k+1} = z^+_k$
$k =k+1$

We apply RTR, as described in Algorithm <ref>, to problem (<ref>) or (<ref>). If a first-order critical point is sought, set $\varepsilon_H = \infty$. The second-order version of RTR provably converges to second-order critical points for any initialization under a weak decrease condition in the subproblems and satisfies global worst-case complexity bounds matching their unconstrained counterparts, as was shown in <cit.>. The local convergence rate is quadratic for an appropriate choice of parameters <cit.>. We note that there is no guarantee on the quality of the stationary point, due to the nonconvexity. Nonetheless, we see in Section <ref> that the method performs very well in practice for nonlinear matrix recovery. We introduce the following assumptions.

There exists $f^*>-\infty$ such that $f(x)\geq f^*$ for all $x\in\M$.

The cost functions of problems (<ref>) and (<ref>) are nonnegative, so assumption:bounded_below is satisfied throughout this paper.
We also state a regularity assumption on the gradient and Hessian of the pullback which was introduced in <cit.>. In Appendix <ref>, we detail how these conditions relate to the smoothness of the Riemannian derivatives and discuss the practicality of these assumptions for problem (<ref>).

[(Lipschitz gradient of the pullback)]
There exists $L_g\geq 0$ such that for all $z = (X,\U)\in \M$, the pullback $\hat f_z = f\circ \Retr_{z}$ has Lipschitz continuous gradient with constant $L_g$, that is, for all $\eta \in \Trm_{z}\M$, it holds that

\begin{equation}
\left| \hat f_z(\eta) - [f(z) + \inner{\eta}{\grad f(z)}]\right| \leq \dfrac{L_g}{2}\norm{\eta}^2.
\label{eq:grad-pullback-lipschitz}
\end{equation}

[(Lipschitz Hessian of the pullback)]
There exists $L_H\geq 0$ such that, for all $z = (X,\U)\in \M$, the pullback $\hat f_z = f\circ \Retr_{z}$ has Lipschitz continuous Hessian with constant $L_H$, that is, for all $\eta \in \Trm_{z}\M$, it holds that

\begin{equation}
\left| \hat f_z(\eta) - \left[f(z) + \inner{\eta}{\grad f(z)} + \dfrac{1}{2}\inner{\eta}{\nabla^2 \hat f_z(0_{z})[\eta]} \right]\right| \leq \dfrac{L_H}{6}\norm{\eta}^3.
\label{eq:hess-pullback-lipschitz}
\end{equation}

Algorithm <ref> is flexible in that it does not specify how the subproblems are solved. We discuss the implementation of RTR in Section <ref>. For the complexity results, the following decreases in the model for first- and second-order steps are required.

There exists $c_2>0$ such that all first-order steps $\eta_k$ satisfy

\begin{equation}
\hat m_k(0_{z_k}) - \hat m_k(\eta_k) \geq c_2 \min\left(\Delta_k, \dfrac{\varepsilon_g}{c_0}\right)\varepsilon_g.
\end{equation}

There exists $c_3>0$ such that all second-order steps $\eta_k$ satisfy

\begin{equation}
\hat m_k(0_{z_k}) - \hat m_k(\eta_k) \geq c_3 \Delta_k^2\varepsilon_H.
\end{equation}

There exists $c_0 \geq 0$ such that, for all first-order steps, $\norm{H_k}\leq c_0$ and $H_k$ is radially linear, that is, $\forall \alpha \geq 0$ and $\forall \eta \in \Trm_{z_k} \M$, $H_k[\alpha \eta] = \alpha H_k[\eta]$.

There exists $c_1\geq 0$ such that, for all second-order steps,

$$\left| \inner{\eta_k}{ \left(\nabla^2 \hat{f}_{z_k} (0_{z_k}) - H_k\right) [\eta_k] } \right| \leq \dfrac{c_1 \Delta_k}{3} \norm{\eta_k}^2.$$

In addition, for all second-order steps, $H_k$ is linear and symmetric.

Define the following constants

\begin{align}
\lambda_g &= \dfrac{1}{4}\min\left(\dfrac{1}{c_0}, \dfrac{c_2}{L_g + c_0}\right) &\text{ and } && \lambda_H &= \dfrac{3}{4}\dfrac{c_3}{L_H + c_1}.
\end{align}

The following (sharp) worst-case bound for Riemannian trust-region was recently established.

Under assu:lip-gradient-tr, assu:1st-order-decrease-rtr, assu:rtr-H-first-order and assuming $\varepsilon_g \leq \dfrac{\Delta_0}{\lambda_g}$, Algorithm <ref> produces an iterate $z_{N_1}$ satisfying $\norm{\grad f(z_{N_1})} \leq \varepsilon_g$ with

\begin{equation}
N_1 \leq \bigO(1/\varepsilon^2_g).
\end{equation}

Furthermore, if $\varepsilon_H <\infty$ then under additionally assu:lip-hessian-tr, assu:2nd-order-decrease-rtr, assu:rtr-H-second-order and assuming $\varepsilon_g \leq \dfrac{c_2}{c_3}\dfrac{\lambda_H}{\lambda_g^2}$ and $\varepsilon_H \leq \dfrac{c_2}{c_3\lambda_g}$, Algorithm <ref> also produces an iterate $z_{N_2}$ satisfying $\grad f(z_{N_2})\leq \varepsilon_g$ and $\lambdamin(H_{N_2}) \geq - \varepsilon_H$ with

\begin{equation}
N_1\leq N_2 \leq \bigO\left( \dfrac{1}{\varepsilon^2\varepsilon_H}\right).
\end{equation}

§ AN ALTERNATING MINIMIZATION ALGORITHM

In this section, we propose an alternating minimization algorithm to solve (<ref>) or (<ref>) (Algorithm <ref>). This comes from the natural separation of the variables into two blocks $X$ and $\U$, yielding two distinct minimization subproblems. Alternating minimization type methods have been very popular in recent years to solve large-scale nonconvex problems <cit.>. This is due to their good practical performances and ease of implementation, as often one or both of the subproblems have a closed form solution. Strictly speaking, this is still a Riemannian optimization approach, as all iterates will be feasible for the constraints. But this section describes a two-block coordinate minimization, whereas the previous section was considering full block variants.

We set the initial guess $X_0$ as any solution of the underdetermined linear system $\mathcal{A}(X)=b$ and $\U_0$ as the span of the $r$ leading singular vectors of $\Phi(X_0)$. The framework is as follows, for $k \geq 0$:With $\U_{k}$ fixed, solve

\begin{equation}
X_{k+1} = \left\lbrace
\begin{aligned}
& \underset{X\in \R^{n \times s}}{\argmin}
& & \fronorm{\Phi (X) - P_{\U_{k}}\Phi (X)}^2 \\
& %\text{subject to}
%& & P_{\Omega}(X)=Y.
& & \mathcal{A}(X)=b.
\end{aligned}
\right.
\label{eq:s1}
\end{equation}

With $X_{k+1}$ fixed, solve

\begin{equation}
\U_{k+1} = \left\lbrace
\begin{aligned}
& \underset{\U}{\argmin}
& & \fronorm{\Phi (X_{k+1}) - P_{\U}\Phi (X_{k+1})}^2 \\
& %\text{subject to}
& &\U\in \Grass(N,r).
\end{aligned}
\right.
\label{eq:s2}
\end{equation}

This separation of the variables takes advantage of the fact that problem (<ref>), even though nonconvex, is solved to global optimality by computing the $r$ leading left singular vectors of the matrix $\Phi(X_{k+1})$. The result is a consequence of the celebrated Eckart-Young-Mirsky theorem, which gives the best rank $r$ approximation in Frobenius norm of a matrix by the $r$ leading terms of the singular value decomposition <cit.>. In particular, let $\Phi(X_k) = \sum_{i=1}^{\min(N,s)} \sigma_i u_i v_i^\top$, then

\begin{equation}
\U_{k+1} = \txt{span}(u_1, \dots, u_r)
\label{eq:truncate_svd}
\end{equation}

is a global minimizer of (<ref>)[Note that the solution need not be unique, in the case where $\sigma_r = \sigma_{r+1}$.]. This truncated singular value decomposition (SVD) is denoted by  in Algorithm <ref>. The singular vectors can be computed to high accuracy in polynomial time <cit.>.

Problem (<ref>) is in general hard to solve to global optimality.  The difficulty comes from the nonconvexity of the cost function, which is due to $\Phi$. One can choose from a variety of first-or second-order methods to find an approximate first-or second-order critical point.  We will present the merits of both possibilities in Section <ref>.

i) First-order version of alternating minimization When only gradient information is available, a first-order method will be used to minimize subproblem (<ref>). For the sake of illustration, in Algorithm <ref> we present a projected gradient descent with line search for (<ref>). The gradient of the cost function is projected onto the null space of $\A$. This ensures that the iterates remain in the feasible set $\LAb$. The line search is a classical backtracking with an Armijo condition for sufficient decrease. Variants in the line search or even constant step sizes are possible.

ii) Second-order version of alternating minimization
In the subproblem (<ref>), it is possible to use second-order methods to speed up the convergence and reach a higher accuracy. We apply RTR on the affine manifold $\LAb$ with the Hessian of the cost function $\nabla_{XX}^2 f(X_k,\U_k)$ in the model.

*	Algorithm <ref> details the alternating minimization where gradient descent with an Armijo line search, a standard inexact procedure in nonconvex optimization, is applied in the subproblem (<ref>). The  Armijo line search is described in Algorithm <ref>.

Alternating minimization scheme for Problem (<ref>) or (<ref>)
[1]
Given: The sensing matrix $A\in \R^{m\times ns}$, measurements $b\in \R^m$, tolerances $\varepsilon_u \geq 0, \varepsilon_x\geq 0$, an estimation of the lifted rank $r$.
Set $k=0$
Find $X_0$ that satisfies $AX_0 = b$
$U_{0} = \verb= truncate_svd= (\Phi(X_0))$  Equation (<ref>)
$ \fronorm{\grad_X f(X_{k},\U_k) } > \varepsilon_x$ or $ \|\grad_\U f(X_{k},\U_k) \| > \varepsilon_u$
Set $X_k^{(0)} = X_k, i=0$
Choose $\varepsilon_{x,k}$ using Equation (<ref>) or (<ref>)
$ \|\grad_X f(X_{k}^{(i)},\U_{k}) \| > \varepsilon_{x,k}$
$ \grad_X f(X_{k}^{(i)},\U_{k}) =  \Prm_{\Trm \LAb}\Big(\nabla_X f(X_{k}^{(i)},\U_{k}) \Big)$ Equation (<ref>)
$ \alpha_k^{(i)} = $ Armijo$\left((X_{k}^{(i)},\U_{k}), - \grad_X f(X_{k}^{(i)},\U_{k})\right)$ Algorithm <ref>
$X_{k}^{(i+1)} = X_k^{(i)} - \alpha_k^{(i)} \grad_X f(X_{k}^{(i)},\U_{k})$
$i =i+1$
$X_{k+1} = X_k^{(i)}$
$ \|\grad_\U f(X_{k+1},\U_k) \| \leq \varepsilon_{u}$
$U_{k+1} = U_k$
$U_{k+1} = \verb= truncate_svd=(\Phi(X_{k+1}))$
$k = k +1$
Output: $(X_k,\U_k)$ such that $\|\grad f(X_k,\U_k)\| \leq \varepsilon_u + \varepsilon_x$.

Armijo($z_k$, $d_k$): Line search with Armijo condition
INPUT: Function $f$ and gradient $\grad_X f$, current iterate $(X_{k}^{(i)},\U_{k})$ and a descent direction $d_k$ such that $\langle \grad_X  f(X_{k}^{(i)},\U_{k}),d_k \rangle < 0$, a sufficient decrease coefficient $\beta\in ]0,1[$, initial step $\alpha_0>0$ and $\tau \in ]0,1[$.
OUTPUT: Step size $\alpha_k^{(i)}$.
[1]
Set $\alpha=\alpha_0$.
$  f(X_k^{(i)} + \alpha d_k,\U_k)  > f(X_{k}^{(i)},\U_{k}) + \beta \alpha \langle  \grad_X  f(X_{k}^{(i)},\U_{k}),d_k \rangle $
$\alpha = \tau \alpha$.
Set $\alpha_k^{(i)} = \alpha$.

*iii) Accuracy of the subproblems solution
Algorithm (<ref>) alternatively solves the subproblems (<ref>) and (<ref>). For the solution of (<ref>), there is no incentive to solve it to very high accuracy early on in the run of the algorithm, as we could still be far from convergence and the variable $\U$ might still change a lot. At iterations $k$, we use the following stopping criterion for some $\varepsilon_{x,k}>0$,

\begin{equation}
\fronorm{\grad_X f(X_{k+1},\U_{k})} \leq \varepsilon_{x,k}.
\end{equation}

We propose the two following strategies for the choice of $\varepsilon_{x,k}$,

\begin{equation}
% \tag{Greedy Decrease condition}
\varepsilon_{x,k} = \varepsilon_x \txt{ for all } k, \label{eq:greedy_decrease}
\end{equation}

\begin{equation}
% \tag{Adaptive Decrease condition}
\varepsilon_{x,k} = \max\left(\varepsilon_x,  \theta \fronorm{\grad_X f(X_k, \U_{k})} \right) \txt{ for some user-chosen } 0<\theta<1.\label{eq:adaptive_decrease}
\end{equation}

To solve (<ref>), it is possible to use a randomized SVD procedure. The randomized SVD is a stochastic algorithm that approximately computes the singular value decomposition of a matrix that exibits a low-rank pattern <cit.>. The matrix must be of low rank or have a fast decay in its singular values for the random SVD to be accurate. As the iterates $X_k$ converge towards the solution $M$, the matrix $\Phi(X_k)$, for which we have to compute  an SVD, becomes low-rank and therefore it is natural to use a randomized SVD in Algorithm <ref>. In the first iterations, for a random starting point of the algorithm, the feature matrix is not expected to be low-rank. In those case, the random SVD should not be used. When the matrix $\Phi(X_k)$ is only approximately low-rank, we can apply power iterations to make the singular values decrease faster. This will make the randomized decomposition more costly, but will improve the accuracy of the computed SVD.

Our strategy is as follows, choose two parameters $0<\tau_1\ll\tau_2<1$. As long as $f(X_{k+1},\U_k) > \tau_2$, use an exact SVD algorithm, without randomization. When $\tau_1 < f(X_{k+1},\U_k)\leq  \tau_2$ we know that the energy of $\Phi(X_{k+1})$ is mostly contained in the  span of $\U_k$ which has dimension $r$.  So we are justified in using a randomized $\SVD$, which we start up with a step of the power method to improve the accuracy. When $f(X_{k+1},\U_k)\leq \tau_1$, the matrix $\Phi(X_{k+1})$ is even closer to low-rank and we no longer need to apply a power iteration before computing the randomized SVD.

§ CONVERGENCE OF THE ALTERNATING MINIMIZATION ALGORITHM

In this section we present convergence results for the alternating minimization Algorithm <ref>. We consider a first-order version where subproblem (<ref>) is minimized with the gradient descent method and the Armijo backtracking line-search (Algorithm <ref>).
We will first show asymptotic convergence of the gradient norms to zero. We also give a worst-case global complexity bound on the number of iterations necessary to achieve a small gradient from an arbitrary initial starting point. Note that we chose the Armijo linesearch for the sake of example. Minor adjustments of the proof below allow to prove similar results for other minimization methods in subproblem (<ref>).

§.§ Assumptions

There exists constants $L_x$ and $L_u$ (which are both independent of $X$ and $\U$) such that for all $z=(X, \U) \in \M$, the pullback $\hat f_z = f\circ \Retr_z$ has a Lipschitz continuous gradient in $X$ and $\U$, with constants $L_x$ and $L_u$ respectively. That is, for all $(\eta_x,\eta_u) \in  \Trm_{(X,\U)}\M$,

\begin{equation}
\left| f\circ \Retr(\eta_x,0) - [f(X,\U) + \langle \grad_X f(X,\U), \eta_x \rangle] \right| \leq \dfrac{L_x}{2} \norm{\eta_x}^2
\label{eq:lipschitz_x_am}
\end{equation}

\begin{equation}
\left| f\circ \Retr(0,\eta_u) - [f(X,\U) + \langle \grad_\U f(X,\U), \eta_u \rangle] \right| \leq \dfrac{L_u}{2} \norm{\eta_u}^2.
\label{eq:lipschitz_u_am}
\end{equation}

In words, this means that, in each variable, the pullback is well approximated by its first-order Taylor approximation.

Note that if assu:lip-gradient-tr holds, then assumption:pullback_lipschitz_altmin holds with $L_x = L_u = L_g$.

Let us discuss, under which conditions on the kernel one can ensure that assu:lip-gradient-tr or assumption:pullback_lipschitz_altmin are satisfied for the cost function of (<ref>). The following discussion requires to use the exponential map as the retraction. The exponential map follows geodesics along the manifold in directions prescribed by tangent vectors. Using the exponential map on $\M$ is not a restriction, as the exponential map on the Grassmann manifold is computable <cit.> and the exponential map on $\LAb$ is the identity.

Consider the cost function of (<ref>) and assume that the retraction being used is the exponential map. If $\D K(X,X)$ is Lipschitz continuous over $\LAb$, then (<ref>) holds where $L_x$ is the Lipschitz constant of $\D K(X,X)$.  If $\fronorm{K(X,X)}\leq M$ for all $X\in \LAb$, condition (<ref>) holds with $L_u = 2M$.

See appendix <ref>

The conditions listed in Proposition <ref> on the kernel and its derivatives are not always satisfied or can be difficult to verify. For instance, the Gaussian kernel $K_G$ is bounded above on $\LAb$, but the monomial kernel $K_d$ is not for any degree $d\geq 1$.  For the Gaussian kernel, the map $\D K_G(X)$ is always Lipschitz continuous on $\LAb$. Whereas for the monomial kernel,  the map $\D K_d(X)$ is Lipschitz continuous for $d\leq 2$, and only locally Lipschitz continuous for $d\geq 3$.

Fortunately, the picture is much simpler if the sequence of iterates $(X_k)_{k\in \mathbb{N}}$ generated by Algorithm <ref> or <ref> is contained in a bounded set. This ensures that we can find Lipschitz constants such that the bounds of assu:lip-gradient-tr, assu:lip-hessian-tr and assumption:pullback_lipschitz_altmin holds at every iterate of the algorithm (and trial points if any), which is all that is needed in the convergence analysis.

Consider the cost function of either (<ref>) or (<ref>) and apply Algorithm <ref> or Algorithm <ref>  with the exponential map as the retraction. If the convex hull of the sequence of iterates $(X_k)_{k\in \mathbb{N}}$ and the trial points of the algorithm is a bounded set, then (<ref>),
(<ref>) and (<ref>)-(<ref>) hold at every iterate $\left(X_k,\U_k\right)_{k\in\mathbb{N}}$ and trial points of the algorithm.

See appendix <ref>.

§.§ Global convergence results

We now carry on with the convergence analysis of the alternating minimization algorithm. The next lemma is an adaptation of the classical descent lemma for the SVD step.
(Descent lemma based on <cit.>)
Let assumption:bounded_below and assumption:pullback_lipschitz_altmin hold for $f:\M \longrightarrow \R$.
Then, for any $k\geq 0$,

\begin{equation}
f(X_{k+1},\U_k) - f^* \geq \dfrac{1}{2L_u}\|\grad_\mathcal{U} f(X_{k+1},\U_k)\|^2.
\end{equation}

where $L_u$ is the Lipschitz constant of the gradient of the pullback (assumption:pullback_lipschitz_altmin) and $f^*$ is the lower bound of $f$ (assumption:bounded_below).

See Appendix <ref>

Throughout this section we use the following notation. Let the number of gradient steps between $X_k$ and $X_{k+1}$ be $n_k \geq 0$ and the intermediate iterates,

$$X_k=X_k^{(0)},~X_k^{(1)},~X_k^{(2)},~\dots,~ X_k^{(n_k)}=X_{k+1} .$$

The next lemma gives upper and lower bounds on the step sizes given by the Armijo linesearch. This is an adaptation of a standard argument for linesearch methods <cit.> where the constraint $\A(X)=b$ is added.

Under assumption:pullback_lipschitz_altmin, for the direction
$-\grad_X f(X_k^{(i)},\U_k) \in T_{X_k}\LAb$,
the linesearch Algorithm <ref> produces a step size $\alpha_k^{(i)}$ that satisfies

\begin{equation}
\underline{\alpha} := \min \left\lbrace \alpha_0,   \dfrac{ 2\tau(1  - \beta)  }{L_x } \right\rbrace\leq \alpha_k^{(i)} \leq \alpha_{0}
\end{equation}

and produces the following decrease

\begin{equation}
f(X_k^{(i)},\U_k) -  f(X_k^{(i+1)},\U_k) \geq   \beta \alpha \fronorm{\grad_X  f(X_k^{(i)},\U_k)}^2,
\label{eq:armijo_decrease_grad}
\end{equation}

where $X_k^{(i+1)} = X_k^{(i)} - \alpha_k^{(i)} \grad_X f(X_k^{(i)},\U_k)$.

See Appendix <ref>

We are now ready to prove global convergence of the alternating minimization algorithm.

Let assumption:pullback_lipschitz_altmin hold for $f:\M \longrightarrow \R$ from (<ref>) or (<ref>).
Let $\varepsilon_x=0$, $\varepsilon_u=0$ and use Equation (<ref>) to set $\varepsilon_{x,k}$, for any starting point $(X_0,\U_0)\in \M$,  Algorithm <ref> produces a sequence $\Big(X_{k}, \U_{k}\Big)_{k\in \mathbb{N}}$ such that

\begin{equation}
\lim_{\substack{k\longrightarrow\infty }}\norm{ \grad f(X_{k}, \U_{k})}=0.
\label{eq:grad_to_zero}
\end{equation}

First note that $f$ is bounded below by $f_* = 0$. For any $k \geq 0$,

\begin{align}
\norm{ \big( \grad_X f(X_{k}, \U_{k}),  \grad_\U f(X_{k}, \U_{k}) \big)} &\leq \fronorm{ \grad_X f(X_{k}, \U_{k})}+ \norm{ \grad_\U f(X_{k}, \U_{k})}\\
&\leq  \fronorm{ \grad_X f(X_{k}, \U_{k})}
\label{eq:boundgradUX}
\end{align}

since $\varepsilon_u =0$.
Given that each step is non-increasing,

\begin{align*}
f(X_{k},\U_{k})- f(X_{k+1},\U_{k+1}) &\geq f(X_{k},\U_{k})- f(X_{k+1},\U_{k})  \\
&\geq f(X_{k},\U_{k}) - f(X_{k}^{(1)},\U_{k})  \\
&\geq \beta \alpha_k^{(0)} \fronorm{ \grad_X f(X_k, \U_k)}^2\\
&\geq \beta \underline{\alpha} \fronorm{ \grad_X f(X_k, \U_k)}^2,
%  \label{eq:1_armijo_decrease}
\end{align*}

where we used Lemma <ref> about Armijo steps.
Summing over all iterations gives a telescopic sum on the left-hand side. For any $\bar{k}\geq 0$,

\begin{equation}
f(X_0,\U_0)- f^* \geq f(X_0,\U_0)- f(X_{\bar{k}},\U_{\bar{k}})
\geq \beta  \underline{\alpha} \sum_{k=0}^{\bar{k}} \fronorm{ \grad_X f(X_k, \U_k)}^2.\\
\end{equation}

The series is convergent since it is bounded independently of $\bar{k}$. Letting $\bar{k}\longrightarrow \infty$ gives $\sum_{k=0}^{\infty} \fronorm{ \grad_X f(X_k, \U_k)}^2<\infty$ and therefore

\begin{equation}
\lim_{k\longrightarrow\infty}\fronorm{\grad_X f(X_k,\U_k)}=0.
\end{equation}

We have $\norm{\grad_\U f(X_k,\U_k)}=0$ for all $k\geq 0$ since $\varepsilon_u=0$. This corresponds to taking exact SVDs. Taking $k\longrightarrow \infty$ in equation (<ref>) gives convergence of the gradient norms to zero (<ref>).

Let assumption:pullback_lipschitz_altmin hold for $f:\M \longrightarrow \R$ from (<ref>) or (<ref>).
Let $\varepsilon_x>0$, $\varepsilon_u>0$, for any starting point $z_0=(X_0,\U_0)\in \M$, Algorithm <ref> produces a sequence $\Big(X_{k}, \U_{k}\Big)_{k\in \mathbb{N}}$ such that

\begin{equation}
\norm{\Big(\grad_X f(X_{k}, \U_{k}), \grad_\U f(X_{k}, \U_{k}) \Big)}\leq  \varepsilon_x + \varepsilon_u
\end{equation}

is achieved using at most $N_{grad}$ gradient steps and $N_{svd}$ singular value decompositions with

\begin{align}
N_{grad}  &\leq \dfrac{L_x(f(z_0) -f_*)}{C \varepsilon^2_x} & \text{ and } && N_{svd} &\leq \dfrac{2L_u(f(z_0) -f_*)}{\varepsilon^2_u},
\end{align}

where $C= \beta \tau (1-\beta) $ is a constant depending on parameters of the line search (<ref>).

Note that $f$ is bounded below by $f_*=0$. Define $N_{iter}$ as the number of iterations performed by Algorithm <ref>, the smallest $k$ such that $\norm{\grad_X f(X_k,\U_k) }\leq \varepsilon_x$ and $\norm{\grad_\U f(X_k,\U_k) }\leq\varepsilon_u$.
Let $N_{svd}$ be the number of SVDs that have to be performed to get $\| \grad_{\U}f(X_{k+1},\U_{k})\|\leq \varepsilon_u$, at which point the algorithm would return without performing another computation.
For any $k\leq N_{svd}$, from Lemma <ref> we have

\begin{equation}
f(X_{k+1},\U_k)-f_* \geq \frac{1}{2L_u}\|\grad_\U f(X_{k+1},\U_k)\|^2 \geq \frac{1}{2L_u}\varepsilon_u^2.
\end{equation}

Summing from $k=0$ to $N_{svd}$ gives,

\begin{equation}
f(z_0) - f_* \geq f(z_0) - f(z_{N_{svd}}) \geq \sum_{k=0}^{N_{svd}}\frac{\varepsilon_u^2}{2L_u} = \frac{\varepsilon_u^2 N_{svd}}{2L_u}.
\end{equation}

Hence, this bounds the number of SVDs to ensure  $\| \grad_{\U}f(X_{k+1},\U_{k})\|\leq \varepsilon_u$, as

\begin{equation}
N_{svd} \leq 2L_u \frac{(f(z_0) - f_*)}{\varepsilon_u^2}.
\end{equation}

For $0\leq i\leq n_k -1$, we have $ \fronorm{ \grad_X f(X_{k}^{(i)},\U_k)}^2 \geq  \varepsilon_{x,k}^2 $ by definition since the stopping criterion is $ \fronorm{\grad_X f(X_{k}^{(n_k)},\U_k)}\leq \varepsilon_{x,k}$.
Combined with the Armijo decrease this gives

\begin{equation}
f(X_{k}^{(i)},\U_k)- f(X_{k}^{(i+1)},\U_k) \geq  \alpha^{(i)}_k\beta \| \grad_X f(X_{k}^{(i)},\U_k)\|^2 \geq \alpha^{(i)}_k\beta \varepsilon_{x,k}^2.
\label{eq:armijo_decreaseUX}
\end{equation}

We sum these bounds for the $n_k$ gradient steps from  $X_k$ to $X_{k+1}$,

\begin{align}
\sum_{i=0}^{n_k-1}\left[ f(X_{k}^{(i)},\U_k)- f(X_{k}^{(i+1)},\U_k) \right] \geq \sum_{i=0}^{n_k-1} \alpha^{(i)}_k\beta \varepsilon_{x,k}^2 . \\
\intertext{Using that the step sizes $\alpha^{(i)}_k$ are bounded below by $\underline{\alpha} = \tau (1-\beta)/L_x$ (Lemma~\ref{lemma:step_bound}),}
f(X_{k},\U_k)- f(X_{k+1},\U_k) \geq n_k \underline{\alpha}\beta \varepsilon_{x,k}^2 ~~ \forall k.
\end{align}

The SVD is nonincreasing, meaning $f(X_{k},\U_k) - f(X_{k+1},\U_{k+1}) \geq  f(X_{k},\U_k)- f(X_{k+1},\U_k)$. This yields,

\begin{equation}
\begin{aligned}
f(X_{k},\U_{k})- f(X_{k+1},\U_{k+1}) \geq n_k \underline{\alpha}\beta \varepsilon_{x,k}^2 \geq n_k \underline{\alpha}\beta \varepsilon_{x}^2 ~~ \forall k \leq N_{iter}, \\
\end{aligned}
\end{equation}

as both (<ref>) and (<ref>) satisfy $\varepsilon_{x,k} \geq \varepsilon_x$. We sum once again over the iterations,

\begin{equation}
\begin{aligned}
f(X_{0},\U_0)- f_* \geq f(X_{0},\U_{0})- f(X_{N_{iter}+1},\U_{N_{iter}+1}) \geq \sum^{N_{iter}}_{k=0} n_k\underline{\alpha}\beta \varepsilon_{x}^2  \\
\end{aligned}
\end{equation}

We conclude that

\begin{equation}
\dfrac{(f(z_{0})- f_*) }{ \underline{\alpha}\beta \varepsilon_x^2} \geq \sum^{N_{iter}}_{k=0} n_k =: N_{grad} .
\end{equation}

A similar algorithm using fixed step sizes will also converge, provided the step sizes are small enough.

If the Armijo linesearch in Algorithm <ref> is replaced by a gradient descent with constant stepsizes $\alpha$ satisfying $\alpha < \dfrac{2}{L_x}$, Algorithm <ref> also converges

\begin{equation}
\dis \lim_{k\longrightarrow\infty}\norm{ \Big(\grad_X f( X_{k}, \U_{k}),   \grad_\U f( X_{k},\U_{k}) \Big)}=0.
\label{eq:grad_to_zero2}
\end{equation}

We also have the worst case bound

\begin{equation}
N_{grad}  \leq \dfrac{L_x(f_0 -f_*)}{\alpha \varepsilon^2_x}.
\label{eq:N_grad}
\end{equation}

We derive the usual descent lemma from Lipschitz continuity of the gradient. This gives

\begin{equation}
f\big(X_k - \alpha \grad_X f(X_k,\U_k), \U_k\big) \leq f(X_k,\U_k) - \alpha \norm{\grad_X f(X_k,\U_k)}^2 + \alpha^2L_x/2\norm{\grad_X f(X_k,\U_k)}^2
\end{equation}

which simplifies to

\begin{equation}
f(X_k,\U_k) - f(X_k^{(1)},\U_k) \geq (\alpha - \alpha^2L_x/2)\norm{\grad_X f(X_k,\U_k)}^2.
\end{equation}

This bound replaces the Armijo decrease of Equation (<ref>).
The rest of the proofs from Theorems <ref> and <ref> holds verbatim with stepsize $\alpha$ for every iteration. Note that for $\alpha>0$, the factor $(\alpha - \alpha^2L_x/2)$ is positive for $\alpha<2/L_x$ and is maximized at $\alpha= 1/L_x$.

§ CONVERGENCE OF THE ITERATES USING THE KURDYKA-LOJASIEWICZ PROPERTY

This section proves convergence of the sequence of iterates to a unique stationary point for a simplified version of the alternating minimization scheme. This section considers an algorithm where only one gradient step is performed in between the truncated SVDs (Algorithm <ref>). This is analogue to the algorithm described in <cit.> which does not provide theoretical convergence guarantees. Our observations indicate that Algorithm <ref> is expected to behave similarly to Algorithm <ref> in the limit. Asymptotically, there is usually only one gradient step needed between two truncated SVDs. It is only in the early iterations that Algorithm <ref> differs by making several gradient steps in between SVDs. For the purpose of this theoretical section, we will assume that the singular value decompositions in Algorithm <ref> are exact and not approximated or randomized. This corresponds to setting $\varepsilon_u=0$ in Algorithm <ref>. This section is written using the notation of a feature matrix $\Phi$ as in problem (<ref>), but the results apply similarly to problem (<ref>) if one assumes that the Lipschitz condition assumption:lipschitz applies to a kernel $K$ instead of $\Phi$.

A simple alternating minimization scheme for Problem (<ref>) or (<ref>)
[1]
Given: The sensing matrix $A\in \R^{m\times ns}$, measurements $b\in \R^m$, a tolerance $\varepsilon _x>0$, estimation of the lifted rank $r$.
Set $k=0$
Find $X_0$ that satisfies $AX_0 =b$.
$U_{0} = \verb= truncate_svd=(\Phi(X_{0}))$ Equation (<ref>)
$ \|\grad_X f(X_{k},\U_k) \| > \varepsilon_x$
$ \grad_X f(X_{k},\U_{k}) = \mathrm{P}_{\mathrm{T}\LAb} (\nabla_X f(X_{k},\U_{k}))$ Equation (<ref>)
$ \alpha_k = \text{Armijo}\left( (X_k,\U_k), - \grad_X f(X_k,\U_k)\right)$ Algorithm <ref>
$X_{k+1} = X_k - \alpha_k \grad_X f(X_{k},\U_{k})$
$U_{k+1} = \verb= truncate_svd=(\Phi(X_{k+1}))$  exact SVD, not randomized
Output: $(X_k,\U_k)$ such that $\|\grad f(X_k,\U_k)\| \leq \varepsilon_x$.

Let us define a distance on the manifold $\M$ (Equation (<ref>)).

Given two subspaces $\U_1,\U_2 \in \Grass(N,r)$, the canonical angles $\theta_i $ for $i = 1,\dots,r$ are defined as  $\theta_i= cos^{-1}(\sigma_i)$ where $\sigma_i$ are the $r$ singular values of $U_1^\top U_2 $, with $\range(U_1) = \calU_1$ and $\range(U_2) = \calU_2$. For all \(
\U_1, \U_2 \in \Grass\) define \( \dist( \U_1, \U_2) := \sqrt{ \sum_{i=1}^r \sin^2 \theta_i}.\)
This gives a distance on $\M$,

\begin{equation}
\forall (X_1,\U_1), (X_2,\U_2) \in \M,~~~~\dist\Big( (X_1,\U_1), (X_2,\U_2) \Big) := \sqrt{ \fronorm{X_1 - X_2 }^2 + \sum_{i=1}^r \sin^2 \theta_i}.
\label{eq:dist-sin-theta}
\end{equation}

We will prove finite length of the sequence of iterates in this metric on $\M$. A more mainstream approach to define the distance between $\calU_1$ and $\calU_2$ on the Grassmann would be to use $\fronorm{\Theta}$ instead of $\fronorm{\sin \Theta}$, where $\Theta = \txt{diag}(\theta_i)$ is the diagonal matrix containing the principal angles. We do so because the distance (<ref>) makes it easier to derive perturbation bounds for the SVD and is equivalent to the usual distance.

The following assumption ensures a useful non-degeneracy of the spectrum of the feature matrix.
[(Gap between the singular values)]
There exists $\delta > 0$ such that, the accumulation points of the sequence generated by Algorithm <ref> satisfy

\begin{equation}
\sigma_{r}(\Phi(X)) - \sigma_{r+1}(\Phi(X)) \geq \delta>0.
\end{equation}

This property ensures that the minimizer of the function $f( X,\cdot) : \mathrm{Grass}(N,r) \to \R$ is well defined, i.e., that its truncated SVD is unique.  As $\sigma_{r+1}(\Phi(X)) \geq 0$, this assumption also implies that

\[ \sigma_{r}(\Phi(X)) \geq \delta >0 \]

In particular it means that we cannot overestimate the rank. If the true rank is $r-1$, then $\sigma_r =0$ and the assumption does not hold. Let us stress that this is an artefact of the convergence proof and does not imply poor practical performance of the algorithm when the rank is overestimated. We investigate this in the numerics Section <ref>. We will need Assumption <ref> to derive a Lipschitz continuity result on the SVD.

We now show the two main lemmas, inspired by <cit.>.

Assume that Algorithm <ref> generates a bounded sequence of iterates. Then, there exists $\rho_2 >0$ such that, for all $k\in \mathbb{N}$,

\begin{equation}
\norm{\grad f(X_{k+1},\U_{k+1})} \leq \rho_2 \dist\Big(  (X_{k+1},\U_{k+1}) , (X_{k},\U_{k}) \Big)
\label{eq:rho2-1}
\end{equation}

with $\rho_2 :=2(L_g + 1/\alpha_0)$.

The expression

\begin{equation}
X_{k+1} = X_{k} -\alpha_k \grad_X f(X_{k},\U_k)
\end{equation}

\begin{equation}
\grad_X f(X_{k+1},\U_{k+1}) = (  X_{k} - X_{k+1})/\alpha_k + \grad_X f(X_{k+1},\U_{k+1}) - \grad_X f(X_{k},\U_{k})
\end{equation}

We now show that $\grad f$ is $L_g$-Lipschitz continuous in the sense of Definition <ref> for some $L_g\geq 0$. Since the sequence of iterates is bounded, the Hessian, which is continuous, is bounded over the closure of the convex hull of the iterates. By Proposition <ref>, $\grad f$ is Lipschitz continuous.
Using the triangular inequality gives

\begin{align}
\norm{\grad_X f(X_{k+1}, \U_{k+1})} &\leq \norm{  X_{k+1} - X_k}/\alpha_0+ \norm{ \grad_X f(X_{k+1},\U_{k+1}) - \grad_X f(X_{k},\U_{k})}\\\nonumber
&\leq    \dist\Big(  (X_{k+1},\U_{k+1}) , (X_{k},\U_{k}) \Big)/\alpha_0 + L_g   \dist\Big(  (X_{k+1},\U_{k+1}) , (X_{k},\U_{k}) \Big)\\ \nonumber
&\leq (1/\alpha_0+L_g)  \dist\Big(  (X_{k+1},\U_{k+1}) , (X_{k},\U_{k}) \Big).
\end{align}

This gives (<ref>) recalling that, since $\grad_\calU f(X_{k+1}, \calU_{k+1}) = 0$,

\begin{align}
\norm{\grad f(X_{k+1},\U_{k+1})} &= \norm{\grad_X f(X_{k+1},\U_{k+1})}
\end{align}

Further auxiliary results are needed.

Let $C, B \in \mathbb{R}^{N \times s}$, such that

\[ C = \begin{bmatrix}
U_1^C & U_2^C & U_3^C\\
\end{bmatrix}
\begin{bmatrix}
\Sigma_1^C & 0 \\
0 & \Sigma_2^C  \\
0 &  0   \\
\end{bmatrix}
\begin{bmatrix}
V_1^{C \top}  \\
V_2^{C \top}
\end{bmatrix},  \qquad B = \begin{bmatrix}
U_1^B & U_2^B & U_3^B\\
\end{bmatrix}
\begin{bmatrix}
\Sigma_1^B & 0 \\
0 & \Sigma_2^B \\
0 & 0 \\
\end{bmatrix}
\begin{bmatrix}
V_1^{B \top}  \\
V_2^{B \top}
\end{bmatrix}, \]

with $U_1^{C,B} \in \R^{N \times r}$, $U_2^{C,B} \in \R^{N \times (s-r)}$, $U_3^{C,B} \in \R^{N \times (N-s)}$, $\Sigma_1^{C,B} \in \R^{r \times r}$, $\Sigma_2^{C,B} \in \R^{(s-r) \times (s-r)}$, $V_1^{C,B} \in \R^{s \times r}$, $V_1^{C,B} \in \R^{s \times (s-r)}$. We assume $\Sigma_1^{C}(1,1) \geq \Sigma_1^{C}(2,2) \geq \dots \geq \Sigma_1^{C}(r,r) \geq \Sigma_2^{C}(1,1) \geq \dots \geq \Sigma_2^{C}(s-r,s-r)$, and similarly for $B$. If there exists $\delta > 0$ such that

\begin{equation}
\min_{\substack{ 1\leq i \leq r\\
r+1\leq j \leq s-r}}  \rvert \Sigma_1^B(i,i) - \Sigma_{2}^C(j,j) \rvert \geq \delta
\label{eq:wedin_assumption1}
\end{equation}

\begin{equation} \Sigma_1^B(r,r) \geq \delta,  \end{equation}

\begin{equation}
\fronorm{\sin \Theta}^2  \leq \dfrac{2 \fronorm{B-C}^2}{\delta^2}
\label{eq:wedin2}
\end{equation}

with $\Theta$ the matrix of the principal angles between $\mathrm{range}(U_1^C)$ and $\mathrm{range}(U_1^B)$.

The following lemma is a direct consequence of Wedin's theorem.

Let $Y, \check Y \in \R^{N \times s}$. Consider the singular value decomposition of $Y =  U \Sigma V^\top$, with $\sigma_1 \geq \sigma_2 \geq \dots \geq \sigma_s$, $ U \in \R^{N \times s}$, $\Sigma = \diag(\sigma_i) \in \R^{s \times s}$ and $V \in \R^{s \times s}$. Let us also write $U_r :=  U(:,1:r)$ a matrix whose columns span the left principal subspace associated to the $r$ largest singular values. Similarly, $\check Y =  \check U \check \Sigma \check V^\top$ with $\check U_r := \check U(:,1:r)$.
If there exists $\delta > 0$ such that $\sigma_r - \sigma_{r+1} \geq \delta$ and  $\check \sigma_r - \check\sigma_{r+1} \geq \delta$, then

\begin{equation}
\dist( \check \U,  \U)^2 \leq \frac{2}{\delta^2} \fronorm{\check Y - Y}^2,
\end{equation}

where $\dist(\U, \check \U) = \sqrt{\sum_{i=1}^r \sin(\theta_i)^2}$ (with $\theta_i$ the principal angles between $\U$ and $ \check \U$) is the distance between the subspaces $\U$ and $\check \U$. Note that the condition $\sigma_r - \sigma_{r+1} \geq \delta$ can be replaced by $\check \sigma_r - \check \sigma_{r+1} \geq \delta$, by switching the definition of $Y$ and $\check Y$.

The result follows from the $\sin \Theta$ bound (<ref>) in Wedin's theorem. Let us verify the assumptions. From the assumptions we know  that $\sigma_{r} > 0$ and $\check \sigma_{r} > 0$. If Wedin's theorem does not apply, assumption (<ref>) does not apply with $B =Y$ and $ C = \check Y $, and neither does it apply with the roles reversed ($B =\check Y$ and $C = Y$). In that case, we must have $\sigma_i = \check \sigma_{j}$, for some $i \leq r$, $j \geq r+1$, and $\check \sigma_l = \sigma_{m}$, for $l \leq r, m \geq r+1$.
However, since the singular values are ordered decreasingly, this gives:

\[\sigma_{m} \leq \sigma_i = \check \sigma_{j} \leq \check \sigma_l = \sigma_{m}, \]

which implies

\[\sigma_{m} = \sigma_i = \check \sigma_{j} = \check \sigma_l. \]

This is a contradiction with $\sigma_r - \sigma_{r+1} \geq \delta$ and $\check \sigma_r - \check \sigma_{r+1} \geq \delta$.
Therefore, these conditions guarantee that Wedin's theorem applies.

In the next lemma we combine the previous bound with the Lipschitz continuity of $\Phi$.
[(Lipschitz continuity of the features)]
For Problem (<ref>), there exists $L_{\Phi} \geq 0$ such that for any $X_{k}, X_{k+1}$ produced by Algorithm <ref>,
$\fronorm{\Phi(X_{k+1}) - \Phi(X_k)} \leq L_{\Phi} \fronorm{ X_{k+1} - X_k}$. For Problem (<ref>), there exists $L_{K}\geq 0$ such that $\fronorm{K(X_{k+1},X_{k+1}) - K(X_k,X_k)} \leq L_{\Phi} \fronorm{ X_{k+1} - X_k}$.

If we assume that the sequence $\left(X_k\right)_{k\in \mathbb{N}}$ is bounded, which we do in the main result of this section (Theorem <ref>), then it is sufficient for the features and kernel to be locally Lipschitz in order for assumption:lipschitz to hold. We also note that if the level set $\{	(X,\U) \in \M: f(X,\U) \leq f(X_0,\U_0)\}$ is bounded, then the iterates are contained in a bounded set since Algorithm <ref> is a descent method.

Let assumption:sigma and assumption:lipschitz hold. It follows that

\[\dist( \U_{k},\U_{k+1})^2 \leq \frac{2 L_\Phi^2}{\delta^2}  \fronorm{X_{k+1} - X_k}^2.  \]

By definition of $\U_{k}$, Lemma <ref> ensures that

\begin{equation}
\dist( \U_{k},\U_{k+1})^2 \leq \frac{2}{\delta^2} \fronorm{\Phi(X_{k+1}) - \Phi(X_k)}^2.
\end{equation}

Indeed, \(
\U_k = \texttt{truncate-svd}(\Phi(X_{k}))\) is composed of the $r$ first left singular vectors of $\Phi(X_{k})$. They are uniquely defined due to assumption:sigma.
The result then follows from the Lipschitz continuity of $\Phi$.

This lemma allows us to show the following crucial result.

Assume that assumption:sigma and assumption:lipschitz hold at $X_k$ and $X_{k+1}$. Then, there exists $\rho_1 > 0$, independent of $k$, such that the iterates of Algorithm <ref> satisfy

\begin{equation}
f(X_k,U_k) - f(X_{k+1},U_{k+1}) \geq \rho_1\dist\Big( (X_k,\U_k),(X_{k+1},\U_{k+1})\Big)^2.
\label{eq:suff_decrease_rho1}
\end{equation}

From the Armijo decrease of Lemma <ref>,

\begin{equation}
f(X_k,\U_k) - f(X_{k+1},\U_k) \geq   \dfrac{\beta}{\alpha_0}\fronorm{X_k  - X_{k+1}}^2,
\end{equation}

where $\alpha_0$ is the largest step allowed by the backtracking. Set $M^2 = 2L_\Phi^2/\delta^2$.
Using that $f(X_{k+1},\U_{k+1}) \leq f(X_{k+1}, \U_k)$, we get

\begin{align}
f(\U_k, X_k) - f(\U_{k+1}, X_{k+1}) &\geq f(\U_k, X_k) - f(\U_{k}, X_{k+1}) \\
&\geq \dfrac{\beta}{\alpha_0} \fronorm{X_{k+1} - X_k}^2   \\
&= \dfrac{\beta}{\alpha_0(1 + M^2)} (1 + M^2) \fronorm{X_{k+1} - X_k}^2  \\
&\geq  \dfrac{\beta}{\alpha_0(1 + M^2)}\left( \fronorm{X_{k+1} - X_k}^2 + \dist^2( \U_k,\U_{k+1}) \right)
\end{align}

where the last inequality comes from Lemma <ref>. This establishes (<ref>)
with $\rho_1 := \dfrac{\beta}{\alpha_0(1 + M^2)}$.

We now show convergence of the gradient norms to zero for Algorithm <ref>.

Set $\varepsilon_x = 0$,	for any starting point $z_0 = (X_0,\U_0)\in \M$, Algorithm <ref> applied to (<ref>) or (<ref>) produces a sequence $\Big(X_{k}, \U_{k}\Big)_{k\in \mathbb{N}}$ such that

\begin{equation}
\lim_{\substack{k\longrightarrow\infty }}\norm{ \grad f(X_{k}, \U_{k})}=0.
\label{eq:grad_to_zeroKL}
\end{equation}

Using Lemmas <ref> and <ref> gives,

\begin{align}
f(X_k,\U_k) - f(X_{k+1},\U_{k+1}) &\geq \rho_1\dist\Big( (X_k,\U_k),(X_{k+1},\U_{k+1})\Big)^2\\
&\geq \rho_1/\rho_2 \norm{ \grad f(X_{k+1},\U_{k+1})}^2.
\end{align}

The telescopic sum is bounded above by $f(z_0)$ independently of $k$ since $f$ is nonnegative, which ensures $\lim_{k\longrightarrow\infty}  \norm{ \grad f(X_{k+1},\U_{k+1})} = 0$.

The bounds in Lemmas <ref> and <ref> are standard and hold for most descent methods. The values $\rho_1, \rho_2$ depend on the specifics of the algorithm used <cit.>.

We now define the Kurdyka-Lojasiewicz inequality on Riemannian manifolds, which was already introduced in <cit.>.
A locally Lipschitz function $f: \M \longrightarrow \R$ satisfies the Kurdyka-Lajasiewicz inequality at $x \in \M$ iff there exist $\eta \in (0,\infty)$, a neighborhood $V \subset \M$ of $x$, and a continuous concave function $\kappa : [0,\eta] \longrightarrow [0,\infty [ $ such that

• $ \kappa(0) = 0$,

• $\kappa$ is continuously differentiable on $(0,\eta)$,

• $\kappa' > 0$ on $(0,\eta)$,

• For every $y \in V$ with $f(x) < f(y) < f(x) + \eta$, we have

\kappa'(f(y) - f(x)) \norm{\grad f(y)} \geq 1.

If $f$ satisfies the KL inequality at every point $x\in \M$ we call $f$ a KL function.

Let $\{ a_k \}_{k\in \mathbb{N}}$ be a sequence of nonnegative numbers. If $\dis \sum_{k = 1}^\infty \dfrac{a_k^2}{a_{k-1}}$ converges, then $\dis \sum_{k = 1}^\infty a_k$ converges as well.

This is a standard result. A proof can be found in <cit.>.

Assume that Algorithm <ref> is applied to problem (<ref>) or (<ref>), for case study <ref> or <ref>, generates a bounded sequence $(X_k,\U_k)_{k\in \mathbb{N}}$. If assumption:sigma and assumption:lipschitz hold, then, the sequence has finite length, that is

\begin{equation}
\sum_{k=1}^\infty \dist\Big((X_k,\U_k) ,(X_{k+1},\U_{k+1}) \Big) <\infty.
\label{eq:finite_length}
\end{equation}

Therefore $\left(X_k,\U_k\right)_{k\in \mathbb{N}}$ converges to a unique point $(X_*,\U_*)$, which is a critical point of $f$ on $\Mcal$.

For case-studies <ref> and <ref>, the lifting is an algebraic or exponential function. These functions are known to be the KL functions <cit.>, so the cost function $f$ is a KL function (Definition <ref>).
For convenience, we write $z_k = (X_k,\U_k)$. Since the sequence $(z_k)_{k\in \mathbb{N}}$ is bounded, there is a subsequence $(z_{k_q})_{q\in \mathbb{N}}$ which converges to some $\bar{z}\in\M$. Let $\omega(z_0)$ denote the set of limit points for some starting point $z_0$. The set $\omega(z_0)$ is bounded by assumption and clearly closed, therefore it is compact.  We want to show that $\omega(z_0)$ is a singleton, i.e. $\omega(z_0)= \{\bar{z}\}$. The function $f$ is continuous, which implies  $\lim_{q\longrightarrow \infty} f(X_{k_q},\U_{k_q}) = f(\bar{z})$. Since $f(z_{k_q})_{q\in \mathbb{N}}$ is non-increasing, the function $f$ is also constant on $\omega(z_0)$. Since $f$ is a KL function, for every point $z \in \omega(z_0)$, there exists a neighbourhood $V_z$ of $z$ and a continuous concave function $\kappa_z: [0,\eta_z]\rightarrow [0,\infty [$ of class $C^1$ on $]0,\eta_z[$ with $\kappa_z(0)= 0$, $\kappa'_z >0$ on $]0,\eta_z[$ such that, for all $y \in V(z)$ with $f(z) < f(y) < f(z) + \eta_z$, we have

\begin{equation}
\kappa_z'(f(y) - f(z)) \norm{\grad f(y)} \geq 1.
\end{equation}

By compactness of $\omega(z_0)$, we find a finite number of points $\bar{z}_1, \dots, \bar{z}_p \in \omega(z_0)$ such that $\dis\cup_{i=1}^p V_{\bar{z}_i}$ covers $\omega(z_0)$. We choose $\varepsilon >0$, such that $V:= \{z \in \M: \dist\big(z,\omega(z_0)\big) < \varepsilon\}$ is contained in $ \cup_{i = 1}^p V_{\bar{z}_i}$. Then, we set $\eta = \min_{i=1,\dots,p}\eta_{\bar{z}_i}$, $\kappa'(t) = \max_{i = 1,\dots, p} \kappa'_{\bar{z}_i}(t) $ and $\kappa(t) = \int_0^t \kappa'(\tau)d\tau$. We claim that for every $z \in \omega(z_0)$, and $y \in V$, with $ f(z) < f(y) < f(z) + \eta$, we have

\begin{equation}
\kappa'(f(y) - f(z)) \norm{\grad f(y)} \geq 1.
\end{equation}

Indeed, there exists some $\bar{z}_i$ such that $y \in V_{\bar{z}_i}$. Then, from the definition of $\kappa'$ and the fact that $f$ is constant on $\omega(z_0)$,

\begin{equation}
\kappa'(f(y) - f(z)) \norm{\grad f(y)} \geq \kappa'_{\bar{z}_i}(f(y) - f({\bar{z}_i})) \norm{\grad f(y)} \geq 1.
\end{equation}

For $\eta >0$ given above,  there exists $k_0 $ such that for all $ k>k_0$,

\begin{equation}
f(z_k) < f(\bar{z}) +\eta.
\end{equation}

By definition of the accumulation points,  there exists $ k_1 $ such that  for all $ k>k_1$,

\begin{equation}
\dist(z_k,\omega(z_0)) <\varepsilon.
\end{equation}

Since $\sigma_r (\Phi(X))> \sigma_{r+1}(\Phi(X))$ for any $X$ such that $(X,\U)\in \omega(z_0)$ (assumption:sigma), by continuity of the singular values, there exists $\bar{\delta}>0$ such that for all points $z_k$ satisfying $\dist(z_k,\omega(z_0)) <\bar{\delta}$, we have $ \sigma_r (\Phi(X_k))> \sigma_{r+1}(\Phi(X)_k)$. Again by definition, there exists $k_2$ such that for all $ k>k_2$,

\begin{equation}
\dist(z_k,\omega(z_0)) <\bar{\delta}.
\end{equation}

For $k>l= \max\{k_0, k_1,k_2\}$, we have

\begin{equation}
\kappa'(f(z_k) - f(\bar{z})) \norm{\grad f(z_k) } \geq 1.
\end{equation}

Using  $\norm{\grad f(z_{k+1})} \leq \rho_2 \dist\Big(z_k ,z_{k+1} \Big)$ (Equation (<ref>)), gives

\begin{equation}
\kappa'(f(z_k) - f(\bar{z})) \geq \dfrac{1}{\rho_2 \dist\Big(z_{k-1},z_k \Big)}.
\end{equation}

Concavity of $\kappa$ gives

\begin{equation}
\kappa\Big(f(z_k) - f(\bar{z})\Big) - \kappa\Big(f(z_{k+1}) - f(\bar{z})\Big) \geq \kappa'\Big(f(z_k) - f(\bar{z})\Big) \Big(f(z_k) - f(z_{k+1})\Big).
\end{equation}

Since $k>l\geq k_2$, we have that $\rho_1 \dist^2\Big(z_k,z_{k+1} \Big) \leq f(z_k) - f(z_{k+1})$ by Equation (<ref>) ,

\begin{equation}
\kappa\Big(f(z_k) - f(\bar{z})\Big) - \kappa\Big(f(z_{k+1}) - f(\bar{z})\Big) \geq  \dfrac{1}{\rho_2  \dist\Big(z_{k-1},z_k \Big)} \rho_1  \dist^2\Big( z_k,z_{k+1} \Big),
\end{equation}

and so

\begin{equation}
\dfrac{ \dist^2\Big(z_k,z_{k+1} \Big)}{ \dist\Big(z_{k-1},z_k\Big)} \leq \dfrac{\rho_2}{\rho_1}\kappa\Big(f(z_k) - f(\bar{z})\Big) - \kappa\Big(f(z_{k+1}) - f(\bar{z})\Big)\\
\label{eq:736}
\end{equation}

For any $N>l$, we sum (<ref>) for all $l\leq k \leq N$, using that the right hand side is a telescopic sum,

\begin{align}
\nonumber
\sum_{k \geq l}^N \dfrac{ \dist^2\Big(z_k ,z_{k+1} \Big)}{\dist\Big(z_{k-1} ,z_k \Big)} &\leq \sum_{k \geq l}^N \dfrac{\rho_2}{\rho_1}  \Big[ \kappa\Big(f(z_k) - f(\bar{z})\Big) - \kappa\Big(f(z_{k+1}) - f(\bar{z})\Big)\Big] \\\nonumber
&\leq  \dfrac{\rho_2}{\rho_1}\Big[   \kappa\Big(f(z_l) - f(\bar{z})\Big) - \kappa\Big(f(z_{N}) - f(\bar{z})\Big)\Big] \\\nonumber
&\leq  \dfrac{\rho_2}{\rho_1}\Big[   \kappa\Big(f(z_l) - f(\bar{z})\Big) - \kappa\Big(f(\bar{z}) - f(\bar{z})\Big)\Big] \\
&=   \dfrac{\rho_2}{\rho_1}   \kappa\Big(f(z_l) - f(\bar{z})\Big).
\label{eq:740}
\end{align}

where we used that $f(\bar{z}) \leq f(z_N)$ , $\kappa$ is increasing and $\kappa(0)=0$. Letting $N\longrightarrow \infty$ in (<ref>), we deduce that the left-hand side of (<ref>) converges.
By Lemma <ref>, $\sum_{k \geq l}^\infty \dist\Big(z_k,z_{k+1} \Big) $ also converges and therefore

\begin{equation}
\sum_{k=1}^\infty  \dist\Big(z_k ,z_{k+1} \Big) < \infty.
\end{equation}

This concludes the proof and shows finite length of the sequence of iterates, which implies convergence of the Cauchy sequence $\big(X_k, \U_k\big)_{k\in \mathbb{N}}$ to a unique point $(X_*, \U_*)$.

§ NUMERICAL EXPERIMENTS

In this section we validate our approach with numerical results on randomly generated test problems. We also compare the performances of the different algorithms we propose.

§.§ Implementation of the algorithms

Let us describe the implementation of the different algorithms  and variants that are considered.  is a first-order version of alternating minimization (Algorithm <ref>) which uses gradient descent with Armijo linesearch to solve subproblem (<ref>). It uses the monomial kernel (Equation (<ref>)). The degree of the kernel that gives the best results is almost always $d=2$. We set the constant $c=1$ in the monomial kernel.
In Altmin2, a second-order trust region method using the exact Hessian is applied to the minimization of (<ref>). This is the only difference with Altmin1. The default values for the parameters of Algorithm <ref> and the Gaussian and monomial kernels are presented in the table below.

Parameter     Default value     Parameter     Default value

$\varepsilon_x$, $\varepsilon_u$      $10^{-6}$     $c$ in Equation (<ref>)     1

$\varepsilon_{x,k}$     Equation (<ref>)      $\alpha_0$ in Algorithm <ref>     2

$\sigma$ in Equation (<ref>)     2.5     $\tau$ in Algorithm <ref>     0.5

$d$ in Equation (<ref>)     2     $\beta$ in Algorithm <ref>     $10^{-4}$

Our code is available at <https://github.com/flgoyens/nonlinear-matrix-recovery> in both Matlab and Python. We use the Manopt <cit.> and Pymanopt <cit.> libraries for optimization on manifold solvers. The Riemannian trust-region RTR2, which implements Algorithm <ref>, is the corresponding Manopt solver for optimization on manifolds. We used a second-order version with the Hessian in the model and the default parameters of the solver. The maximum number of iterations is set at $500$ for RTR2. In Manopt, the subproblems are solved with a truncated conjugate gradient method and the final termination criterion is only a first-order condition (the norm of the gradient) which we set at $10^{-6}$ for RTR2. Pymanopt uses automatic differentiation and does not require to compute the derivatives by hand, while the Manopt uses finite differences if the Hessian is not given as an input.

§.§ Test problems

We describe the set of parameters that we want to vary and test the dependence of each algorithm with respect to these parameters.

Union of subspaces
Case study <ref> depends on the following parameters: ambient dimension $n$,
number of subspaces, dimension of each subspace, number of points on each subspace. To generate a random union of subspaces, we place the same number of points on each subspace and take subspaces of the same dimension. We calculate a basis for a random subspace and generate each point on that subspace by taking a random combination of the columns of that basis.

For case study <ref>, the parameters defining a point cloud divided in clusters in $\R^n$ are: the number of clusters, the number of points in each cluster and the standard deviation $\sigma_c$ of each cluster. We first generate random centers in $\R^n$.  We then add to each center a cluster of points with multivariate Gaussian distribution with zero mean and covariance $\sigma_c^2 \Id$ with $\sigma_c = 0.5$.

§.§ Testing methodology

Throughout we say that an algorithm successfully recovers the matrix $M\in \R^{n\times s}$ if it returns a matrix $X^{output}$ such that the root mean square error (RMSE) is below $10^{-3}$,

\begin{equation}
\mathrm{RMSE}(M,X^{output}) :=\norm{X^{output} - M}_F/\sqrt{ns} \leq 10^{-3}.
\end{equation}

Our goal is to test the ability of our methods to recover the original matrix $M$. We measure the performance against an increase in difficulty of the problem for several parameters. Parameters that increase the difficulty of the recovery include:

* Reducing the number of measurements $m$;

* Increasing the rank in the feature space.

In the case of unions of subspaces, for a fixed number of points, the rank of $\Phi_d(M)$ depends on the number and the dimension of the subspaces, as indicated by Proposition <ref>. For clusters, the rank increases with the number of clusters.
The undersampling ratio is defined as $\delta = \dfrac{m}{ns}$, it is the number of measurements over the number of entries in $M$.
We present phase transition results to numerically show which geometries can be recovered and which undersampling ratios are needed.
Typical phase transition plots for matrix completion vary the undersampling ratio and the rank of the matrix <cit.>. For union of subspaces, the rank of the feature space is difficult to control, therefore we vary the number and dimension of the subspaces. For each value of the varying parameter, we generate $10$ random matrices $M$ that follow the desired structure. We try to recover each with varying $\delta$ from $0.1$ to $0.9$ for a random initial guess. If the RMSE is below $10^{-3}$ in the maximum number of iterations allowed by the algorithm, we consider the recovery to be successful. The phase transition plots record which of the 10 random problems is solved for each configuration. In Figures <ref> through <ref> the grayscale indicates the proportion of problems solved, with white $= 100\%$ of instances solved and black $=0\%$.

§.§ Numerical results

§.§.§ Comparing the performance of RTR and Alternating minimization

Figure <ref> compares the performance of RTR2 (Algorithm <ref> using a second-order Taylor model), Altmin1 and Altmin2 which are first- and second-order alternating minimization (Algorithm <ref>). We chose a problem of matrix completion over a union of subspaces. We find that RTR2  has a local quadratic rate of convergence, which makes it the method of choice if we want to recover $M$ to high accuracy. Both Altmin1 and Altmin2 make faster progress during the early iterations; thus these methods should be considered if the required accuracy is low. We also observed that, in general, the distance to the solution $M$ is of the same order of magnitude than the gradient norm. That is, using an algorithm, such as RTR2, which terminates with a smaller gradient norm yields a greater accuracy for the recovery. We noticed that the first-order methods, as well as the second-order Altmin2, typically stall numerically when the gradient norm is around $10^{-7}$, but that is not the case for RTR2.

Comparing alternating minimization (first-order Altmin1 and second-order Altmin2) with the Riemannian trust-region algorithm (RTR2) for a union of subspaces recovery.

§.§.§ Recovery of unions of subspaces

We now how the parameters at play affect the recovery for data that follows a union of subspaces model.

Degree of the polynomial features
Deciding which degree $d$ to use in practice requires a careful choice. Previous works limit themselves to $d=2$ and $d=3$.  This is understandable because the lifted dimension $N(n,d)$ increases exponentially with $d$, and so the dimension of the Grassmannian variable in (<ref>) becomes too large and the problem becomes practically intractable for even moderate values of $d$ and $n$. For example, $N(n=20, d=5) = 53130$ while $N(n=20, d=2) = 231$.

Therefore, the other natural option is to solve the kernel-based problem (<ref>), where the lifted dimension $N$ does not appear explicitly and the Grassmann has dimension $s\times r$. This is attractive because, a priori, $s$ may not be as large as $N$. However, there are still some requirements on the number of samples needed to allow recovery. We need to ensure that $s \geq N- q$ where $q$ is the number of linearly independent vectors $v$ such that  $v^\top \Phi_d(M) = 0$. That is, $s$ needs to be large enough so that $\rank(\Phi_d(M))$ is not limited by $s$ but the dimension of the variety. Some analysis in <cit.> shows that the number of points $s$ needed to allow recovery increases exponentially with $d$. For that reason, if $d$ is not small, it is not realistic to solve problems where $s$ is large enough to enable recovery. The monomial basis is also known to be ill-conditioned for large degrees. Therefore, the performances of these algorithms for recovery are poor beyond very small values of $d$. In general, $d=2$ seems to give the best results for the majority of data sets.

In Figure <ref>, we solve the recovery problems using RTR2 for an increasing number of data points, and we do that for a monomial kernel of degree one, two and three to compare the recovery that is possible for each degree.
In Figure <ref>, the degree used is $d=1$. For $n=15$, the dimension of the lifted space is $N(15,1) = 16$. For a large number of data points spread over 4 subspaces of dimension 2, the rank of the monomial kernel is 9. This explains why recovery is impossible when $s\leq 9$, since the kernel is not rank deficient at the solution $M$.
In Figure <ref>, the degree used is $d=2$. For $n=15$, the dimension of the lifted space $N(15,2) = 136$. For a large number of data points spread over  4 subspaces of dimension 2, the rank of the monomial kernel is 21. This explains why recovery is impossible when $s\leq 21$, since the kernel is not rank deficient at the solution $M$.
In Figure <ref>, the degree used is $d=3$. For $n=15$, the dimension of the lifted space is $N(15,3) = 816$. For a large number of data points spread over 4 subspaces of dimension 2, the rank of the monomial kernel is 37. This  explains why recovery is impossible when $s\leq 37$, since the kernel is not rank deficient at the solution $M$. We notice that the recovery is still poor for $s\geq  37$.

Phrase transition for data belonging to a union of $4$ subspaces of dimension $2$ in $\mathbb{R}^{15}$ for an increasing number of data points spread across the subspaces. Each square gives the proportion of problems solved over $50$ randomly generated problems, $\text{white}=100\%$ of instance recovered and black $=0\%$,

The dimension of the subspaces that we aim to recover plays a role in the possibility to recover. In Figure <ref>, we increase  the dimension of the subspaces while the other parameters of the data remain fixed. For a fixed number of data points $s$, increasing the dimension of the subspaces increases the rank of the monomial features (see Proposition <ref>), and therefore, if the dimension of the subspaces becomes too large, the recovery is compromised. For subspaces of dimension smaller than 4 in $\R^{10}$, we can observe good recovery depending on the undersampling ratio.

Phrase transition for data belonging to a union of 2 subspaces of increasing dimension in $\mathbb{R}^{10}$ with 20 points on each subspace. Each square gives the proportion of problems solved over $50$  randomly generated problems, $\text{white}=100\%$ of instance recovered and black $=0\%$, using the polynomial kernel of degree $d=2$ as lifting.

The same phenomenon is observed when the number of subspaces is increased for a fixed number of data points, see Figure <ref>.

[Monomial kernel degree $d=1$]

[Monomial kernel degree $d = 2$]

[Monomial kernel degree $d = 3$]

Phrase transition for data belonging to a union of an increasing number of subspaces  of dimension $2$ in $\mathbb{R}^{15}$ with $150$ data points spread across the subspaces. Each square gives the proportion of problems solved over $10$  randomly generated problems, $\text{white}=100\%$ of instance recovered and black $=0\%$.

§.§.§ Clustering with missing data

In the case of clusters (case study <ref>), there is a noise inherent to the model because the kernel at the solution is only approximately low-rank. In fact, the matrix $K^G(M,M)$ has numerical full rank, but there is a big gap in the singular values. These matrices are notoriously difficult to recover in low rank matrix completion. For this reason the recovery error $\fronorm{M-X^*}^2$ rarely converge to high accuracy and recovering the matrix up to 2 digits of accuracy is typical. This completed matrix $X^*$ allows to do a clustering of the data points starting from missing entries. We are interested in determining when the matrix $X^*$ has the same clustering as $M$.
We use the Rand index to measure the compatibility of two different clusterings of the same set <cit.>. We can see in Figure <ref> that for 5 clusters or less, the original clustering can be recovered even though up to 40% of the entries in the original matrix are missing.

Percentage of problems correctly clustered for different sampling rates and increasing number of clusters. 50 random instances generated in each case; $\text{white}=100\%$ of instance correctly clustered and black $=0\%$. Clusters belong to $\R^5$ with 20 points in each cluster.

§.§.§ Robustness to measurement noise

In applications, it is common to assume some noise on the measurements, namely, $\langle A_i, M\rangle = b_i + \xi_i =:  \tilde b_i$ where $ \xi_i\in \R$ is some noise. In the following numerical test we generate white Gaussian noise, i.e. $\xi_i \sim \mathcal{N}(0,\sigma^2)$ for some variance $\sigma^2$. The problem formulation then becomes Problem (<ref>) with $\lambda>0$ as the penalty parameter that should be tuned based on the noise level.

To recover a matrix without any estimate of the noise level $\sigma$, we use a scheme which increases $\lambda$ over successive calls to the solver while warm starting each solve with the previous solution to (<ref>). We have found that starting with the value $\lambda=10^{-6}$ is satisfactory and we multiply $\lambda$ by a factor 10 at each iteration. Figure <ref> shows, for three different noise levels, the evolution of the solution of (<ref>), labelled $X^*$, as the penalty parameter $\lambda$ increases. We see that $\norm{ \A(X^*)-\tilde b }$ and $\norm{ \A(X^*)-b }$ are roughly equal for a range of values of $\lambda$. As $\lambda$ continues to increase, $\norm{ \A(X^*)-b }$ further decreases which is expected but useless as these solutions give to much weight to the noisy measurements. We find the best value of $\lambda$ to be where the approximated rank of the features ($\norm{\Phi(X^*) - P_\U \Phi(X^*)}$) is minimized in the range of $\lambda$ where $\norm{ \A(X^*)-\tilde b }$ is on a plateau before decreasing further. Table <ref> shows the accuracy of the solution $X^*$ for that choice of $\lambda$. We see that both the infeasibility ($\norm{ \A(X^*)-b }$) and the distance to the solution ($\norm{X^*-M}$) are proportional to the noise level and decreases with the later. This shows that the warm start scheme to find a good value for $\lambda$ in conjunction with Problem (<ref>) handles the presence of noise in the measurements very well.

[$\sigma = 10^{-2}$]

[$\sigma = 10^{-3}$]

[$\sigma = 10^{-4}$]

Solutions for noisy problems as a function of the parameter $\lambda$ on the horizontal axis.

Standard deviation     $\norm{ \A(X^*)-\tilde b }$     $\norm{ \A(X^*)-b }$     $\norm{X^*-M}$     $\norm{\A(M)-b}$

$\sigma = 10^{-2}$     $2*10^{-1}$     $7*10^{-2}$     $8*10^{-2}$      $0.2$

$\sigma = 10^{-3}$     $2*10^{-2} $     $8*10^{-3}$     $8*10^{-3}$     $0.02$

$\sigma = 10^{-4}$     $2*10^{-3} $     $8*10^{-4}$     $9*10^{-4}$     $0.002$

Quality of the solution $X^*$ for different levels of noise $\sigma$ in the measurements.

§.§.§ Robustness to a bad estimate of the rank

In the case of a union of subspaces, the polynomial features are exactly low rank. Then, it is important to have an accurate upper bound on the rank. Recovery is sometimes possible if the upper bound is close to the correct value. If the estimated rank is less than the exact rank or much too large, recovery will normally fail. This intuition is guided by the cost function that we use. If the variable $\U$ is artificially constrained to be the leading singular vectors of $\Phi(X)$, that is $\U= \ttt{truncate-svd}(\Phi(X))$, then it can be substituted and the cost function in (<ref>) simplifies to

\begin{equation}
\min_X \sum_{i = r+1}^{\min(N,s)} \sigma_i^2(\Phi(X)).
\end{equation}

The cost function represents the energy in the tail of the singular value decomposition, where $r$ is the estimation of the rank. In Figure <ref>, the data belongs to a union of 2 subspaces of dimension $2$ in $\R^{15}$, with a total of $s= 150$ data points. With a kernel of degree $2$, the lifted dimension is $N(15,2) = 136$.

[Phase transition]

[Singular values of $k_2(M)$]

Impact of an incorrect estimate of the rank for the completion of a union of subspaces. Each square gives the proportion of problems solved over $50$  randomly generated problems, $\text{white}=100\%$ of instance recovered and black $=0\%$, using the polynomial kernel of degree $d=2$ as lifting.

§ CONCLUSION

In this work,  we study in depth the problem of nonlinear matrix completion where one tries to recover a high rank matrix that exhibit low rank structure in some feature space.
In therms of the use cases considered, in addition to the union of subspaces and algebraic varieties, we propose the use of the Gaussian kernel for clustering problems with missing data, which we believe is novel in the context of nonlinear matrix completion.

We propose a novel formulation for the nonlinear matrix completion problem which is inspired from low-rank matrix completion techniques using the Grassmann manifold.
We then show how Riemannian optimization and alternating minimization methods can be applied effectively to solve this optimization problem.
The algorithms we propose, come with strong global convergence results to critical points and worst-case complexity guarantees. In addition, we show that the alternating minimization algorithm converges to a unique limit point using the Kurdyka-Lojasiewicz property.

We provide extensive numerical results that attest to the efficiency of the approach to recover high-rank matrices drawn from union of subspaces or clustered data. We note that the second-order Riemannian trust-region method allows to recover with high accuracy.  We expose the difficulty of using polynomials of high degree in the monomial kernel, as they require an exponentially increasing number of sample points to allow recovery. Our approach proves to be efficient at clustering a data set despite the presence of missing entries and our approach also shows great robustness against the presence of noise in the measurements.

§.§ Acknowledgement

The authors would like to thank Estelle Massart and Greg Ongie for interesting discussions and their helpful ideas.

§ DERIVATIVES OF COST FUNCTIONS

In this section we compute by hand the first- and second-order derivatives of the cost function that appears in the optimization problem (<ref>). To compute the derivative of a matrix valued function, we write a Taylor expansion and identify the gradient by looking at first order terms. Before computing derivatives for a specific kernel, we look at the Lipschitz continuity properties of the gradient problem, for an arbitrary kernel.

§.§ Lipschitz properties

Let us consider the cost function of (<ref>), defined over $\M = \LAb \times \Grass(s,r)$. We work out conditions on the problem which ensure that the Riemannian gradient is Lipschitz continuous. Lipschitz continuity of a vector field on a smooth manifold is defined as follows using a vector transport.
A vector field $V$ on a connected manifold $\M$ is L-Lipschitz continuous if, for all $x,y\in \M$ with $\dist(x,y) <\textrm{inj}(x)$,

$$\norm{\mathrm{PT}^\gamma_{0\leftarrow 1} V(y) - V(x)} \leq L \dist(x,y),$$

where $\gamma: [0,1] \rightarrow \M$ is the unique minimizing geodesic connecting $x$ to $y$ and $\textrm{PT}^\gamma_{0\leftarrow 1}$ denotes the parallel transport along $\gamma$.

Functions with Lipschitz continuous gradient exhibit the following regularity condition for the pullback $\hat{f} = f\circ \Retr$, provided the retraction used is the exponential map, $\Retr = \Exp$.
If $f:\M \rightarrow \R$ has L-Lipschitz continuous gradient, then

\left| f(\Exp_x(s)) - f(x) - \inner{s}{\grad f(x)} \right| \leq \dfrac{L}{2}\norm{s}^2

for all $(x,s)$ in the domain of the exponential map.

In order to show Lipschitz continuity of the gradient, we use both the definition, and the following proposition which uses an upper bound on the derivative of the gradient.
If $f:\M \rightarrow \R$ is twice continuously differentiable on a manifold $\M$, then $\grad f$ is L-Lipschitz continuous if and only if $\Hess f(x)$ has operator norm bounded by L for all $x$, that is, if for all x we have

\norm{\Hess f(x)} = \max_{\mathclap{\substack{x\in T_x\M \\
\norm{s}=1}}} ~\norm{\Hess f(x)[s]}\leq L.

We first compute the Euclidean gradient with respect to each variable,

\begin{equation}
\nabla_X f(X,\W) = \D K(X)^* \Prm_{\W^\perp}
\end{equation}

\begin{equation}
\nabla_\W f(X,\W) = -2 K(X)\W.
\end{equation}

This naturally gives

\begin{equation}
\nabla^2_{\W\W} f(X,\W)[\Delta] = - 2 K(X)\Delta.
\end{equation}

Consider the cost function of (<ref>) and assume that the retraction being used is the exponential map. If $\D K(X)$ is Lipschitz continuous over $\LAb$, then (<ref>) holds where $L_x$ is the Lipschitz constant of $\D K(X)$.  If $\fronorm{K(X)}\leq M$ for all $X\in \LAb$, condition (<ref>) holds with $L_u = 2M$.

For a given $X\in \LAb$, consider the function $f_\W(X,\cdot): \Grass(s,r) \rightarrow \R$. Its Riemannian Hessian is such that for $\Delta \in \Trm_\W \Grass(s,r)$, $\Hess f_\W(\W,X)[\Delta] = 	-2 \Prm_{\W_\perp} K(x)\Delta$. Hence $\norm{\Hess f_\W(\W,X)}\leq \fronorm{k(X)}$. Using <cit.> Corollary 10.45, the vector field $\grad f_\W (X,\cdot)$ is Lipschitz continuous with constant $L_\W = 2\fronorm{k(x)}$. If the kernel is upper-bounded for all $X\in \LAb$, the constant $L_\W$ is independent of $X$. This implies that (<ref>) holds.

We  also analyze Lipschitz continuity of the vector field $\grad f_X$.

\begin{align}
\fronorm{\grad_X f(X_1,\W) - \grad_X f(X_2, \W) } &= \fronorm{\Prm_{T \LAb} \left( \nabla_X f(X_1,\W) - \nabla_X f(X_2,\W) \right)}\\
&\leq \fronorm{ \nabla_X f(X_1,\W) - \nabla_X f(X_2,\W) }\\
&\leq \fronorm{\left(\D K(X_1)^* - \D K(X_2)^*\right) \Prm_{\W_\perp}}\\
&\leq \norm{\D K(X_1)^* - \D K(X_2)^*}_2 \fronorm{\Prm_{\W_\perp}} \\
&\leq \norm{\D K(X_1) - \D K(X_2)}_2.
\end{align}

If $\D K(X)$ is $L_x$-Lipschitz over $\LAb$, we can write

\begin{equation}
\fronorm{\grad_X f(X_1,\W) - \grad_X f(X_2, \W) } \leq L_x \fronorm{X_1-X_2}
\end{equation}

and $\grad_X f(.,\W)$ is also $L_x$-Lipschitz, where the constant $L_x$ is independent of $\W\in \Grass(s,r)$.  This implies that (<ref>) holds.

Consider the cost function of either (<ref>) or (<ref>) and apply Algorithm <ref> or Algorithm <ref>  with the exponential map as the retraction. If the convex hull of the sequence of iterates $(X_k)_{k\in \mathbb{N}}$ and the trial points is a bounded set, then (<ref>),
(<ref>) and (<ref>)-(<ref>) hold at every iterate $\left(X_k,\U_k\right)_{k\in\mathbb{N}}$ and trial points of the algorithm.

Provided the kernel is a smooth function, the derivatives of the cost function are continuous. As a consequence of the Weierstrass theorem, the derivatives are bounded on the closure of the convex hull of the iterates, which is compact ($\Grass(s,r)$ is compact). If the Hessian is bounded on the closure of convex hull of the iterates, the gradient is Lipschitz continuous on that set (Proposition <ref>) and therefore assu:lip-gradient-tr and assumption:pullback_lipschitz_altmin hold with the exponential map as the retraction (Propositon <ref>). The continuity of the third-order derivatives implies assu:lip-hessian-tr in a similar way.

§.§ Monomial kernel

We wish to find the Euclidean derivatives of the cost function in (<ref>) for the monomial kernel defined in (<ref>). First, we make the following developments. Up to first order in $\Delta_X$,

\begin{equation}
\begin{aligned}
K_d(X+\Delta_X) &= K_d(X) + d K_{d-1}(X)\odot(X^\top \Delta_X + \Delta_X^\top X) + \mathcal{O}(\Delta_X^2).
\end{aligned}
\end{equation}

For $\Delta_W = W_\perp B$, we have

\begin{equation}
\begin{aligned}
P_{W+\Delta_W} &= WW\transpose + W\Delta_W\transpose + \Delta_W W\transpose,\\
&= P_W + W\Delta_W\transpose + \Delta_W W\transpose +  \mathcal{O}(\Delta_W^2).\\
\end{aligned}
\end{equation}

Let us write
$$f(X,W) = \tr\left(K_d(X) - P_{W}K_d(X)\right). $$

We find the gradient in $X$ using direct computation,

\begin{equation}
\begin{aligned}
\nabla_X f(X,W) &= \nabla_X  \tr\left(P_{W_\perp }K_d(X)\right),\\
&= 2 d X\left(K_{d-1}(X) \odot P_{W_\perp }\right).
\end{aligned}
\end{equation}

since $P_{W_\perp }$ is symmetric.
Quite naturally we find $ \nabla_W f(X,W)$ with the expansion,

\begin{equation}
\begin{aligned}
f(X,W+ \Delta_W ) &= \tr\left(K_d(X) - P_{W+\Delta_W}K_d(X)\right) ,\\
&=  \tr\left(K_d(X) - (P_W + W\Delta_W^T + \Delta_W W^T ) K_d(X)\right) ,\\
&=  \tr\left(  P_{W_\perp }K_d(X) - ( W\Delta_W^T + \Delta_W W^T ) K(_dX)\right) ,\\
&=  \tr\left(  P_{W_\perp }K_d(X)\right)  -\tr  \left( (W\Delta_W^T + \Delta_W W^T ) K_d(X)\right) ,\\
&=  f(X,W)  -\tr  \left( W\Delta_W^T K_d(X)\right) -\tr \left( \Delta_W W^T K(X) )\right) ,\\
&=  f(X,W)  +\langle \Delta_W,-2  K_d(X)W\rangle .\\
\end{aligned}
\end{equation}

Ans so we observe $ \nabla_W f(X,W) = -2  K_d(X)W$.
Quickly we have $ \nabla^2_W f(X,W)[E] = -2  K_d(X)E$. Now we want to find the second derivative in $X$

\begin{equation}
\begin{aligned}
\nabla_X f(X+\Delta_X,W ) &=  2 d (X+\Delta_X \left(K_{d-1}(X+\Delta_X) \odot P_{U_\perp }\right)\\
&=   2 d (X+\Delta_X) \left(  \left[K_{d-1}(X) + (d-1) K_{d-2}(X)\odot(X^\top\Delta_X + \Delta_X^\top X)\right]  \odot P_{W_\perp }\right)\\
&=   2 d X \left( K_{d-1}(X) \odot P_{W_\perp }\right) +  2 d(d-1) X \left(  K_{d-2}(X)\odot(X^\top \Delta_X + \Delta_X^\top X) \odot P_{W_\perp }\right)  \\  				&~~~+ 2 d \Delta_X \left( K_{d-1}(X) \odot P_{W_\perp }\right) + \mathcal{O}(\Delta_X^2)\\
&=   \nabla_X f(X,W )  +  2 d(d-1) X \left(  K_{d-2}(X)\odot(X^\top \Delta_X + \Delta_X^\top X) \odot P_{W_\perp }\right)  \\
&~~~+ 2 d \Delta_X \left( K_{d-1}(X) \odot P_{W_\perp }\right) + \mathcal{O}(\Delta_X^2).\\
\end{aligned}
\end{equation}

And so we identify

\begin{equation}
\nabla_X^2 f(X,W)[\Delta_X] =2 d(d-1) X \left(  K_{d-2}(X)\odot(X^\top \Delta_X + \Delta_X^\top X) \odot P_{W_\perp }\right) + 2 d \Delta_X \left( K_{d-1}(X) \odot P_{W_\perp }\right).
\end{equation}

Now we need the cross derivatives $ \nabla_W \nabla_X f(X,W)[\Delta_W]$ and $ \nabla_X \nabla_W f(X,W)[\Delta_W]$.

\begin{equation}
\begin{aligned}
\nabla_W f(X+\Delta_X,W ) 	&= -2  K(X+\Delta_X)W 	\\
&= 	  -2  \left( K_d(X) + d K_{d-1}(X)\odot(X^\top \Delta_X + \Delta_X^\top X) \right) W \\
&= 	  -2   K_d(X)W -2d\left(  K_{d-1}(X)\odot(X^\top \Delta_X + \Delta_X^\top X) \right) W \\
&= 	 \nabla_U f(X,W )  -2d\left(  K_{d-1}(X)\odot(X^\top \Delta_X + \Delta_X^\top X) \right) W. \\
\end{aligned}
\end{equation}

$$ \nabla_X \nabla_W f(X,W)[\Delta_X] = -2d\left(  K_{d-1}(X)\odot(X^\top \Delta_X + \Delta_X^\top X) \right) W \in s \times r.$$

Similarly we find

\begin{equation}
\begin{aligned}
\nabla_X f(X,W+\Delta_W ) 	&= 2 d X\left(K_{d-1}(X) \odot  P_{(W+\Delta_W)_\perp } \right)	\\
&=  2 d X\left(K_{d-1}(X) \odot ( P_{W_\perp }-W\Delta_W^T - \Delta_W W^T  ) \right)		   \\
&= 2 d X\left(K_{d-1}(X) \odot  P_{W_\perp } \right)	 - 2 d X\left(K_{d-1}(X) \odot  ( W\Delta_W^T - \Delta_W W^T  )\right) \\
&=  \nabla_X f(X,W ) - 2 d X\left(K_{d-1}(X) \odot  ( W\Delta_W^T - \Delta_W W^T  )\right). \\
\end{aligned}
\end{equation}

And so

$$\nabla_W \nabla_X f(X,W)[\nabla_W] = - 2 d X\left(K_{d-1}(X) \odot  ( W\Delta_W^T - \Delta_W W^T  )\right) \in n \times s .$$

In the end

\begin{equation}
\nabla^2 f(X,W) \begin{pmatrix}
\Delta_X\\
\Delta_W \\
\end{pmatrix} = \begin{pmatrix}
\nabla_X^2 f(X,W)[\Delta_X] + \nabla_W \nabla_X f(X,W)[\nabla_W]\\
\nabla_X \nabla_W f(X,W)[\Delta_X] + \nabla_W^2 f(X,W)[\Delta_W].\\
\end{pmatrix}
\end{equation}

§.§ Gaussian kernel

Consider the Gaussian kernel defined in (<ref>). A direct computation gives

\begin{equation}
\nabla_X f(X,W)=  -\dfrac{2}{\sigma^2}X \left( \mathrm{diag} \left(\mathrm{sum}(K^{G}\odot P_{W^\perp},1)\right)-K^{G}\odot P_{W^\perp}\right)
\end{equation}

for $K^{G}$ the Gaussian kernel and $\mathrm{sum}(K^{G}\odot P_{W^\perp},1)$ is the vector with the sum of each column of the matrix $K^{G}\odot P_{W^\perp}$.
And similarly to the monomial kernel above, we have

\begin{equation}
\nabla_W f(X,W) = -2K^{G}(X)W.
\end{equation}

We do not compute the Hessian by hand for the Gaussian kernel. We either use automatic differentiation or finite differences of the gradient in the algorithm.

§ PROOFS FOR SECTION 6 (CONVERGENCE OF THE ALTERNATING MINIMIZATION ALGORITHM)

(Descent lemma based on <cit.>)
Let assumption:bounded_below and assumption:pullback_lipschitz_altmin hold for $f:\M \longrightarrow \R$.
Then, for any $k\geq 0$,

\begin{equation}
f(X_{k+1},\U_k) - f^* \geq \dfrac{1}{2L_u}\|\grad_\mathcal{U} f(X_{k+1},\U_k)\|^2.
\end{equation}

where $L_u$ is the Lipschitz constant of the gradient of the pullback (assumption:pullback_lipschitz_altmin) and $f^*$ is the lower bound of $f$ (assumption:bounded_below).

We follow the development of <cit.>. By Lipschitz continuity of the gradient we have,

\begin{equation}
\left| f\left(X_{k+1}, R_{\calU_k}(\eta)\right) - [f(X_{k+1},\calU_k) + \langle \grad_\calU f(X_{k+1},\calU_k), \eta \rangle] \right| \leq \dfrac{L_u}{2} \norm{\eta}^2~~~~ \forall \eta \in T_\calU\Grass.
\end{equation}

Let $\eta = -\dfrac{1}{L_u} \grad_\calU f(X_{k+1},\calU_k)$ and define $\calU^+ = \Retr_{\calU_k}\left(\dfrac{-1}{L_u} \grad_\calU f(X_{k+1},\calU_k)\right)$ , which gives

\begin{align}
f(X_{k+1},\U^+) &\leq f(X_{k+1},\U_k) + \langle \grad_\U f(X_{k+1},\U_k), \dfrac{-1}{L_u} \grad_U f(X_{k+1},\U_k)\rangle \\\nonumber
&~~~~~+\dfrac{L_u}{2}\|\dfrac{-1}{L_u} \grad_\U f(X_{k+1},\U_k)\|^2\\
&\leq  f(X_{k+1},\U_k) -\dfrac{1}{2L_u}\| \grad_\U f(X_{k+1},\U_k)\|^2.
\end{align}

A fortiori we can conclude

\begin{equation}
f(X_{k+1},\U_k) - f^* \geq    f(X_{k+1},\calU_k) -f(X_{k+1},\U^+) \geq \dfrac{1}{2L_u}\|\grad_\U f(X_{k+1},\U_k)\|^2.
\end{equation}

Under assumption:pullback_lipschitz_altmin, for the direction
$-\grad_X f(X_k^{(i)},\U_k) \in T_{X_k}\LAb$,
the linesearch Algorithm <ref> produces a step size $\alpha_k^{(i)}$ that satisfies

\begin{equation}
\underline{\alpha} := \min \left\lbrace \alpha_0,   \dfrac{ 2\tau(1  - \beta)  }{L_x } \right\rbrace\leq \alpha_k^{(i)} \leq \alpha_{0}
\end{equation}

and produces the following decrease

\begin{equation}
f(X_k^{(i)},\U_k) -  f(X_k^{(i+1)},\U_k) \geq   \beta \alpha \fronorm{\grad_X  f(X_k^{(i)},\U_k)}^2,
\label{eq:armijo_decrease_grad_apx}
\end{equation}

where $X_k^{(i+1)} = X_k^{(i)} - \alpha_k^{(i)} \grad_X f(X_k^{(i)},\U_k)$.

It is clear from the algorithm that $\alpha_k^{(i)} \leq \alpha_{0}$. The Armijo condition also ensures (<ref>).
For any $\alpha>0$, Lipschitz continuity of the gradient gives

\begin{align}
f\Big(X_k^{(i)} - \alpha \grad_X f(X_k^{(i)},\U_k) ,\U_k\Big) &\leq f(X_k^{(i)},\U_k) - \alpha \fronorm{ \grad_X f(X_k^{(i)},\U_k)}^2 \\
&~~~~~~~~+ \alpha^2 \dfrac{L_x}{2} \fronorm{\grad_X f(X_k^{(i)},\U_k) }^2.
\end{align}

Hence the Armijo condition (<ref>) is satisfied whenever

\begin{align}
-\alpha \norm{ \grad_X f(X_k^{(i)},\U_k) }^2 + \alpha^2 \dfrac{L_x}{2} \fronorm{ \grad_X f(X_k^{(i)},\U_k)}^2  \leq -\alpha\beta \fronorm{ \grad_X f(X_k^{(i)},\U_k) }^2,
\end{align}

which simplifies to

\begin{equation}
\alpha \leq\dfrac{ 2(1  - \beta)  }{L_x }=:\alpha_{max}.
\end{equation}

If $\alpha_0$ satisfies Armijo, then $\alpha_k^{(i)} = \alpha_0$. Otherwise, we have $\alpha_k^{(i)} = \tau \alpha_l$ where $\alpha_l> \alpha_{max}$ is the last $\alpha$ that does not satisfy Armijo and $\alpha_{l+1} = \tau \alpha_l$ satisfies Armijo. In this case we have $\alpha_k^{(i)} \geq \tau \alpha_{max}=   \dfrac{ 2\tau(1  - \beta)  }{L_x }$.