1 Introduction↩︎

We aim to minimize a continuously differentiable function \(f\colon {\mathbb{R}}^{m \times n} \to {\mathbb{R}}\) over the set of bounded-rank matrices: \[\begin{align} \label{eq:problem} \min_{X \in {\mathbb{R}}^{m \times n}} f(X) && \text{subject to} && \mathop{\mathrm{rank}}X \leq r, \end{align}\tag{1}\] for some \(r < \min(m, n)\). This is a classical problem: see Section 1.2 for a literature review.

The feasible set is a nonconvex, nonsmooth algebraic variety: \[\begin{align} \mathbb{R}_{\leq r}^{m \times n}= \{X \in {\mathbb{R}}^{m \times n} : \mathop{\mathrm{rank}}(X) \leq r \}. \end{align}\] Such problems are computationally hard in general [2], but we can still aim for local solutions. There exist at least three paradigms to do so. The first consists of methods that optimize directly over \(\mathbb{R}_{\leq r}^{m \times n}\). This notably includes the projected gradient descent (PGD) algorithm [3]–[5] and variants [6]. These are first-order methods with limited local convergence rates (at most linear).

Another approach is to optimize only over the maximal rank stratum \[\begin{align} \mathbb{R}_r^{m \times n}= \{X \in {\mathbb{R}}^{m \times n} : \mathop{\mathrm{rank}}(X) = r\}, \end{align}\] thereby ignoring matrices of rank strictly less than \(r\). That stratum accounts for almost all of \(\mathbb{R}_{\leq r}^{m \times n}\) and it is a smooth manifold. Thus, it is possible to endow it with a Riemannian structure and roll out off-the-shelf second-order optimization algorithms [7]–[9]. Those algorithms generate sequences of matrices of rank \(r\). These may converge to matrices of rank strictly less than \(r\), in which case they “fall off” the manifold: theory for optimization on manifolds breaks down in that event, which makes it difficult to provide a priori guarantees. This is likely to happen when the parameter \(r\) overestimates the actual rank of the solution. Additionally, the Hessian can grow unbounded near lesser rank matrices [9].

The third, and perhaps most popular, class of techniques is to smoothly parameterize the feasible set \(\mathbb{R}_{\leq r}^{m \times n}\). Specifically, let \(\mathcal{M}\) be a smooth manifold and let \(\varphi\colon \mathcal{M} \to {\mathbb{R}}^{m \times n}\) be a smooth map such that \(\varphi(\mathcal{M}) = \mathbb{R}_{\leq r}^{m \times n}\). Then minimizing \(f\) over \(\mathbb{R}_{\leq r}^{m \times n}\) is equivalent to minimizing the lifted function \(g= f\circ \varphi\) over \(\mathcal{M}\), for which we can apply standard (Riemannian) optimization algorithms: this last paradigm is what we explore here.

Unfortunately, even when \(f\) has favorable properties, the lifted cost function \(g= f\circ \varphi\) may not. As an example, consider the prominent parameterization \(\varphi(L, R) = L R^\top\) with domain \(\mathcal{M}= {\mathbb{R}}^{m \times r} \times {\mathbb{R}}^{n \times r}\). We refer to it as the LR parameterization (see also Appendix 7). It reduces the feasible set to a simple Euclidean space, yet the lifted cost function \(g\) has several drawbacks for theoretical analyses:²

1. Unbounded sublevel sets for \(g\).
   This is because the fibers \(\varphi^{-1}(X) = \{(L, R) : L R^\top = X\}\) are unbounded. An indirect effect of this is that \(g\) may not have a Lipschitz continuous gradient even if \(f\) does. This hinders global convergence analyses.

2. No local Polyak–Łojasiewicz (\(\text{P\L}\)) condition around minimizers of \(g\) in general.
   Indeed, say 1 has a unique local minimizer \(X\). Its fiber \(\varphi^{-1}(X)\) is a set of local minimizers for \(g\). If \(g\) is twice continuously differentiable (\({\mathrm{C}^{2}}\)) and satisfies local \(\text{P\L}\) around minimizers, then those form a locally smooth set [10]. Yet \(\varphi^{-1}(X)\) is a nonsmooth set when \(\mathop{\mathrm{rank}}(X) < r\). Thus, \(g\) lacks local \(\text{P\L}\) even in this simple scenario. That hinders fast local convergence guarantees.

Accordingly, a promising direction is to explore other smooth parameterizations of \(\mathbb{R}_{\leq r}^{m \times n}\) whose fibers would be smooth and bounded. As it happens, [1] proposed one such parameterization. They call it a desingularization (a term that comes from algebraic geometry). Their construction relies on the following simple fact: a matrix \(X\) of size \(m \times n\) has rank at most \(r\) if and only if its kernel has dimension at least \(n - r\), that is, its kernel contains a subspace of dimension \(n - r\). The set of subspaces of dimension \(s\) in \({\mathbb{R}}^n\) is the Grassmannian, which we define in terms of orthogonal projectors³ as: \[\begin{align} \mathrm{Gr}(n, s) = \big\{P \in \mathrm{Sym}(n) : P^2 = P \text{ and } \mathop{\mathrm{rank}}(P) = s\big\}, \end{align}\] where \(\mathrm{Sym}(n)\) is the set of real symmetric matrices of size \(n\). In this way, \(\mathrm{Gr}(n, s)\) is an embedded submanifold of \(\mathrm{Sym}(n)\). The desingularization of bounded-rank matrices is the parameterization \[\begin{align} \label{eq:desing} \mathcal{M}= \big\{(X, P) \in {\mathbb{R}}^{m \times n} \times \mathrm{Gr}(n, n - r) : X P = 0\big\} && \text{and} && \varphi(X, P) = X. \end{align}\tag{2}\] We confirm in Proposition 1 of Section 2 that \(\mathcal{M}\) is indeed a smooth manifold of dimension \((m + n - r)r\), which is the same as those of \(\mathbb{R}_{\leq r}^{m \times n}\) and \(\mathbb{R}_r^{m \times n}\). Given a pair \((X, P) \in \mathcal{M}\), the equality \(X P = 0\) implies that the kernel of \(X\) contains a subspace of dimension \(n - r\). As such, the projector \(P\) acts as a certificate that \(\mathop{\mathrm{rank}}X \leq r\).

Figure 1 illustrates the parameterization of \(\mathbb{R}_{\leq r}^{m \times n}\) with \(\mathcal{M}\).

[1] further describe how to compute efficiently with points and tangent vectors of \(\mathcal{M}\), and they build a second-order optimization algorithm that borrows ideas from optimization on manifolds (e.g., retractions) and from constrained optimization (e.g., KKT conditions). This hybrid method does not come with convergence guarantees.

To exploit their parameterization further, we endow \(\mathcal{M}\) with a Riemannian metric and with (a few) retractions. This gives access to a large collection of general-purpose optimization algorithms, including second-order ones [7]. Then, building on work by [12], we exploit the favorable geometry of the fibers of the desingularization to establish both global and local convergence guarantees for optimization on \(\mathbb{R}_{\leq r}^{m \times n}\). It appears that such guarantees had not been jointly available for a single algorithm and for a wide class of cost functions. We conclude with numerical experiments.

1.1 Contributions↩︎

We develop the geometry of \(\mathcal{M}\) beyond that of [1], derive optimization tools with a focus on numerical efficiency, secure both global and local convergence guarantees, and run numerical experiments. More precisely:

• In Section 2, we review the geometry of \(\mathcal{M}\), specify parsimonious numerical representations and provide several retractions. Anticipating our needs for convergence analysis, we study geometric properties of preimages of \(\varphi\).

• In Section 3 we equip \(\mathcal{M}\) with a family of Riemannian structures. We deduce formulas to compute the Riemannian gradient and Hessian of \(g= f\circ \varphi\) in terms of \(f\), efficiently. We also build second-order retractions.

• Table ¿tbl:tab:complexities? summarizes the computational complexity of basic operations.

• We exploit the favorable properties of \(\mathcal{M}\) to obtain convergence guarantees:

  • In Section 4.2 we show that \(\varphi\) preserves compactness of sublevel sets and also certain Lipschitz properties for the cost function. We recover global convergence guarantees (accumulation at stationary points of 1 ) for a wide range of algorithms, in a similar fashion to [12].

  • In Section 4.4 we prove that \(\varphi\) can preserve a type of Polyak–Łojasiewicz condition around minimizers. Assuming convergence to one of those minimizers, we deduce fast local rates for several algorithms, including superlinear rates for second-order algorithms (going beyond [12]).

Theorems 17 and 18 outline these guarantees. Together, they yield algorithms with both global convergence and fast local rates under some assumptions on the cost function.

Finally, in Section 5 we experimentally evaluate algorithms for matrix completion tasks. Our implementation is available in the open-source toolbox Manopt [13] and code for the numerical experiments is in a public GitHub repository.⁴

1.2 Related work↩︎

1.2.0.1 Context and applications.

Low-rank optimization is a classical problem with applications in recommender systems [14], Euclidean embedding [15], neuro-imaging [16] and image denoising [17]. Matrix completion [18], [19] is arguably the most prevalent low-rank optimization problem. There is extensive literature exploring under which conditions a low-rank matrix can be recovered from some of its entries [20]–[23].

There are mainly two approaches to formulate low-rank optimization problems. The first one is to regularize an unconstrained problem (over \(m \times n\) matrices) with a nuclear norm penalty to promote low-rank solutions [24]–[26]. The second one imposes a strict limit on the rank of the matrix, resulting in a constrained problem. This is the direction we follow. (To solve the unconstrained problem with nuclear norm penalty, some algorithms exploit the low-rank structure of the solution by also explicitly imposing a rank constraint [27], [28].)

1.2.0.2 Algorithms.

Many existing algorithms for 1 cater to a specific cost function for a given application: see for example [29]–[38] for weighted low-rank approximation, matrix completion and matrix sensing. Some algorithms are instead general-purpose, in the sense that they can deal with a wide range of cost functions. We compete with those, some of which we list now.

Projected gradient descent techniques [5], [6] optimize directly over the singular variety \(\mathbb{R}_{\leq r}^{m \times n}\): see below.

[8] and [9] propose to optimize over the fixed-rank manifold \({\mathbb{R}}_r^{m \times n}\) to take advantage of Riemannian optimization tools. These techniques are efficient in practice but theory breaks when iterates approach lower-rank strata.

[39] parameterize \(\mathbb{R}_r^{m \times n}\) with quotient geometries. [12] consider several parameterizations of \(\mathbb{R}_{\leq r}^{m \times n}\) (including LR and the desingularization of [1]) and propose a special trust-region method: we build on this extensively in later sections.

1.2.0.3 Global convergence guarantees.

A method is globally convergent if its iterates accumulate only at stationary points. Surprisingly few algorithms are known to satisfy this guarantee on \(\mathbb{R}_{\leq r}^{m \times n}\). We review them here.

[40], [41], and [42] proved that PGD is globally convergent. In contrast, the P\(^2\)GD variant [6] can generate a sequence converging to a non-stationary point [12], [43], [44]. [45], [46] refine P\(^2\)GD to ensure global convergence.

There are also guarantees for optimization through a parameterization \(\varphi \colon \mathcal{M}\to {\mathbb{R}}^{m \times n}\) such that \(\varphi(\mathcal{M}) = \mathbb{R}_{\leq r}^{m \times n}\). In this case, one obstacle is that a critical point of the function \(g= f \circ \varphi\) may not map to a stationary point of 1 . But other correspondences are known (see also Section 4.1). Notably, [47] prove that for the LR parameterization, all second-order critical points of \(g\) map to stationary points of 1 . [12], [43] extend this to the desingularization (and more) to obtain a globally convergent trust-region type algorithm. We build on that to add local guarantees.

1.2.0.4 Local convergence guarantees.

Under non-degeneracy assumptions, several algorithms enjoy local convergence guarantees to minimizers of maximal rank. For example, we can instantiate general Riemannian optimization theorems [7] for the particular case of \({\mathbb{R}}_r^{m \times n}\) [9]. This readily gives local convergence guarantees for standard algorithms (including gradient descent and Newton’s method). Separately, [6] obtain capture-type results and convergence rates for P\(^2\)GD based on Łojasiewicz-type arguments. [48] prove local linear convergence to maximal rank points for an alternating descent algorithm with the LR parameterization. [49] derive similar local guarantees for different alternating rules and step sizes (including line search).

Ensuring fast local convergence to matrices of non-maximal rank is more delicate. There exist only a couple of results in the literature. [6] obtain convergence rates (at most linear) for a specific retraction-free gradient method, and these are effective even at non-maximal rank points. Also, [50] study a general tensor version of the matrix sensing problem. Their setting particularizes to matrix sensing over the constant rank stratum \({\mathbb{R}}_r^{m \times n}\). This yields local linear convergence for Riemannian gradient descent. Remarkably, it holds even when the sequence converges to a point of rank less than \(r\), though it is for a specific cost function.

In this paper, we systematize such results to a broader range of cost functions and algorithms with the desingularization geometry. It also allows for superlinear rates of convergence.

1.2.0.5 Positive semidefinite case.

Another related problem is the minimization of a function over the set of symmetric positive semidefinite matrices with bounded rank. One popular approach is the Burer–Monteiro parameterization \(\varphi \colon Y \mapsto YY^\top\) of the feasible set. Unlike the LR parameterization, the fibers of \(\varphi\) are compact, which is valuable for theoretical analyses. This makes the positive semidefinite case markedly different [51]–[58].

1.3 Conventions↩︎

Given a smooth manifold \(\mathcal{N}\), we let \(\mathrm{T}\mathcal{N}\) denote its tangent bundle, and \(\mathrm{T}_x \mathcal{N}\) the tangent space at a point \(x \in \mathcal{N}\).

The identity matrix of size \(r\) is \(I_r\). The Stiefel manifold is \(\mathrm{St}(n, r) = \{V \in {\mathbb{R}}^{n \times r} : V^\top\mkern-1muV = I_r\}\). The orthogonal group is \(O(n) = \{Q \in {\mathbb{R}}^{n \times n} : Q^\top\mkern-1muQ = I_n\}\). Given a matrix \(V \in \mathrm{St}(n, r)\), we often let \(V_\perp \in \mathrm{St}(n, n - r)\) denote any matrix such that \(\big[\begin{matrix} V & V_\perp \end{matrix}\big]\) is in \(O(n)\). (We also say that \(V_\perp\) is an orthonormal completion of \(V\).)

For a linear operator \(\mathcal{L}\) we let \(\|\mathcal{L}\|_\mathrm{op}\) denote its operator norm. For a matrix \(A\) we let \(\|A\|_2\) denote its spectral norm and \(\|A\|_\mathrm{F}\) its Frobenius norm. The singular values of \(A\) are ordered as \(\sigma_1(A) \geq \sigma_2(A) \geq \sigma_3(A) \geq \cdots\) so that \(\|A\|_2 = \sigma_1(A)\).

For a real-valued function \(f\) on a smooth manifold, we say \(x\) is a critical point if the gradient of \(f\) at \(x\) is zero. When minimizing \(f\) constrained to a set \(S\), we say \(x\) is a stationary point if \(\mathrm{D}f(x)[v] \geq 0\) for all \(v\) in the (Bouligand) tangent cone to \(S\) at \(x\).

2 Geometry of the desingularization↩︎

This section explores the geometry of \(\mathcal{M}\): we rederive its tangent spaces, consider some retractions, and analyze the structure of the fibers of \(\varphi\).

The dimension of \(\mathcal{M}\) is \((m + n - r)r\), which is much smaller than \(m n\) when \(r \ll \min(m, n)\). For efficiency, we should store points and tangent vectors using a number of floating point numbers close to \(\dim \mathcal{M}\). We follow an approach based on the singular value decomposition (SVD) as in [9], with some (but limited) resemblance to the specifics in [1]. Here \(\mathrm{St}\) denotes the Stiefel manifold (details are in Section 1.3).

Note that this representation is not unique (for example when \(X = 0\) any \(U\) applies). However, a point \((X, P)\) is stored numerically as a triplet \((U, \Sigma, V)\) and therefore all computations are expressed in terms of that representation (see also Table ¿tbl:tab:complexities?). If \(V_\perp \in \mathrm{St}(n, n - r)\) is an orthonormal completion of \(V\) then \(P = V_\perp^{} V_\perp^\top\). We often use this in theoretical analyses, but algorithms never have to actually build this prohibitively large matrix.

2.1 Smooth structure and tangent spaces↩︎

The set \(\mathcal{M}\) is an embedded submanifold of \(\mathcal{E}= {\mathbb{R}}^{m \times n} \times \mathrm{Sym}(n)\), as proven by [1]. We provide another proof below, after a technical lemma.

It is easy to see that \(\mathcal{M}\) is closed and connected since \(\mathrm{Gr}(n, n - r)\) is so. These properties are desirable for algorithmic convergence. In contrast, \({\mathbb{R}}_r^{m \times n}\) is not closed in \({\mathbb{R}}^{m \times n}\).

We now derive the tangent spaces of \(\mathcal{M}\). Given a point \((X, P)\), we find a parameterization of the tangent space at \((X, P)\) that relies on its representation (see Definition 1).

The expression given in Proposition 3 reveals a parsimonious representation for tangent vectors. We can represent each with a pair \((K, V_p) \in {\mathbb{R}}^{m \times r} \times {\mathbb{R}}^{n \times r}\), requiring only \((m + n)r\) real numbers.

Given \((U, \Sigma, V)\), this representation of \((\dot{X}, \dot{P})\) is unique. We express (most) tangent vectors with this representation for numerical efficiency. For theoretical analyses we sometimes decompose \(V_p = V_\perp L\) for some matrix \(L \in {\mathbb{R}}^{(n - r) \times r}\).

2.2 Retractions↩︎

Let \(\mathrm{T}\mathcal{M}\) denote the tangent bundle of \(\mathcal{M}\). A retraction is a smooth map \(\mathrm{R}\colon \mathrm{T}\mathcal{M}\to \mathcal{M}\) to navigate on \(\mathcal{M}\). Given a point on \(\mathcal{M}\) and a tangent direction at that point, it provides a new point on \(\mathcal{M}\). We let \(\mathrm{R}_{(X, P)} \colon \mathrm{T}_{(X, P)}\mathcal{M}\to \mathcal{M}\) denote the restriction of a retraction \(\mathrm{R}\) at \((X, P) \in \mathcal{M}\). See [7] and [60] for background and references.

Many retractions are already known for \({\mathbb{R}}_r^{m \times n}\) [61]. [1] proposed a retraction for \(\mathcal{M}\) based on the Q-factor retraction on \(\mathrm{Gr}(n, n - r)\). (This is one popular retraction for the Grassmannian.) We show below that any retraction on \(\mathrm{Gr}(n, n - r)\) yields a retraction on \(\mathcal{M}\) in the same way. (Later, in Section 3.3, we also consider second-order retractions.) We express all computations in terms of the parsimonious representations given in Definitions 1 and 2.

We now specialize Proposition 4 to the Q-factor retraction and briefly demonstrate how to compute it (this can be compared to [1]). Let \((X, P) \in \mathcal{M}\) and \((\dot{X}, \dot{P}) \in \mathrm{T}_{(X, P)}\mathcal{M}\) have representations \((U, \Sigma, V)\) and \((K, V_p)\). Let \(V + V_p = QR\) be a thin QR decomposition, which can be computed in \(O(nr^2)\) operations. Then the retracted point is \((\bar X, \bar P)\) where we define \(\bar X = (X + \dot{X}) Q Q^\top\) and \(\bar P = I - Q Q^\top\). We can efficiently obtain a representation for it. First, form the matrix \(W = (X + \dot{X}) Q = (U \Sigma + K) V^\top\mkern-1muQ + (U \Sigma) V_p^\top\mkern-1muQ\). This can be done in \(O((m + n)r^2)\) operations using the representations of \(X\) and \(\dot{X}\). Now let \(W = \bar U \bar \Sigma H^\top \in {\mathbb{R}}^{m \times r}\) be a thin SVD, which can be computed in \(O(mr^2)\) operations. Then \((\bar X, \bar P)\) has representation \((\bar U, \bar \Sigma, \bar V)\) where \(\bar V = Q H\).

2.3 Geometry of fibers and preimages↩︎

The fiber of a point \(X \in \mathbb{R}_{\leq r}^{m \times n}\) is the set \(\varphi^{-1}(X)\). In this section, we study more generally the preimages \(\varphi^{-1}(\mathcal{S})\) of subsets \(\mathcal{S}\subseteq \mathbb{R}_{\leq r}^{m \times n}\). Their geometries affect both global and local convergence properties of optimization algorithms, as we explain in Section 4. A first notable property is that the preimages of compact sets are compact.

We now seek to understand when \(\varphi^{-1}(\mathcal{S})\) is an embedded submanifold of \(\mathcal{M}\). This will be important to secure the Polyak–Łojasiewicz condition in Section 4.4. We start with the preimages of constant rank strata.

We now generalize this result to embedded submanifolds of \({\mathbb{R}}^{m \times n}\) included in a constant rank stratum, and also give an expression for the tangent spaces.

3 Riemannian geometry↩︎

We now endow the tangent spaces of \(\mathcal{M}\) with inner products to turn it into a Riemannian manifold. This provides a notion of gradients, Hessians, and second-order retractions.

3.1 A family of metrics with parameter \(\boldsymbol{\alpha}\)↩︎

We equip the embedding space \(\mathcal{E}= {\mathbb{R}}^{m \times n} \times \mathrm{Sym}(n)\) with the inner product \[\begin{align} \label{eq:inner-product} \inner{\cdot}{\cdot}_\alpha \colon \mathcal{E}\times \mathcal{E}\to {\mathbb{R}}&& \text{defined as} && \sinner{(Y_1, Z_1)}{(Y_2, Z_2)}{\big}_\alpha = \inner{Y_1}{Y_2} + \alpha \inner{Z_1}{Z_2}, \end{align}\tag{5}\] for some constant parameter \(\alpha > 0\), where the inner products on the right-hand side are the standard Frobenius ones. This generalizes the canonical choice \(\alpha = 1\). We usually omit the subscript \(\alpha\) because there is no ambiguity. The inner product on \(\mathcal{E}\) induces inner products on the tangent spaces of \(\mathcal{M}\), turning it into a Riemannian submanifold of \(\mathcal{E}\). The set \(\mathcal{M}\) is closed in \(\mathcal{E}\) (Proposition 2) so we deduce that \(\mathcal{M}\) is a complete manifold [60].

Such inner products can be computed in \(O((m + n)r)\) operations using the representation of tangent vectors in Definition 2. For this, we define the diagonal matrix \[\begin{align} \label{eq:sfactor} S(\alpha) = 2 \alpha I_r + \Sigma^2, \end{align}\tag{6}\] which appears in multiple expressions below.

At a point \((X, P) \in \mathcal{M}\), the inner product 5 induces a normal space \(\mathrm{N}_{(X, P)}\mathcal{M}\). It is the orthogonal complement of \(\mathrm{T}_{(X, P)}\mathcal{M}\) in \(\mathcal{E}\).

From the normal spaces we can deduce an expression for the orthogonal projectors onto tangent spaces. We use them in Section 3.2 to compute Riemannian gradients and Hessians. The orthogonal projection of \((Y, Z) \in \mathcal{E}\) onto \(\mathrm{T}_{(X, P)}\mathcal{M}\) is \[\begin{align} \mathop{\mathrm{proj}}_{(X, P)}(Y, Z) = \mathop{\mathrm{arg\,min}}_{(\dot{X}, \dot{P}) \in \mathrm{T}_{(X, P)} \mathcal{M}} \big\|(\dot{X} - Y, \dot{P} - Z)\big\|^2, \end{align}\] where the norm is induced by the inner product 5 . The (unique) solution \((\dot{X}, \dot{P})\) is a tangent vector at \((X, P)\). We seek a representation (Definition 2) for it. (See also Lemma 7.)

3.2 Gradients and Hessians↩︎

Consider a function \(f\colon {\mathbb{R}}^{m \times n} \to {\mathbb{R}}\) and let \(g= f\circ \varphi\) be the lifted function over \(\mathcal{M}\). Let \(\nabla f\) and \(\nabla^2f\) denote the Euclidean gradient and Hessian of \(f\), and \(\nabla g\) and \(\nabla^2g\) denote the Riemannian gradient and Hessian of \(g\). In this section we express \(\nabla g\) and \(\nabla^2g\) in terms of \(\nabla f\) and \(\nabla^2f\). Importantly, \(g\) depends only on the first component of the pair \((X, P)\). This considerably simplifies expressions for the derivatives.

3.2.0.1 First-order derivatives.

We first find a representation (Definition 2) for \(\nabla g\).

The Riemannian gradient can be computed efficiently if \(\nabla f(X)^\top\mkern-1muU\) and \(\nabla f(X)V\) can be computed efficiently—see Table ¿tbl:tab:complexities?. This can be the case when \(\nabla f(X)\) is appropriately structured (for example if it is sparse). From the expression of the gradient we deduce the first-order optimality conditions.

3.2.0.2 Second-order derivatives.

In this section we assume that \(f\) (and hence \(g\)) is \({\mathrm{C}^{2}}\). We derive the Hessian of \(g\) as a function of the Hessian of \(f\). The proof is computational and can be found in Appendix 8. Notice again that \(\nabla^2g\) can be computed efficiently if \(\nabla f\) and \(\nabla^2f\) are appropriately structured—see Table ¿tbl:tab:complexities?.

From this, we deduce the second-order optimality conditions. Notice that these are independent of \(\alpha\), which is expected since first- and second-order criticality conditions do not depend on the Riemannian metric.

3.3 Second-order retractions↩︎

Second-order retractions are retractions (Definition 3) that satisfy an additional regularity requirement [7]. As they provide better approximations to geodesics, they can be useful both for numerical and theoretical reasons. In this section we list a couple that are efficiently computable.

To define second-order retractions, we need the following notion. A curve \(c \colon {\mathbb{R}}\to \mathcal{M}\) is also a curve in \(\mathcal{E}\). Let \(\ddot c\) denote its (extrinsic) acceleration in \(\mathcal{E}\). Further let \(c''\) denote its (intrinsic, Riemannian) acceleration as a curve on \(\mathcal{M}\). It holds that \(c''(t) = \mathop{\mathrm{proj}}_{\mathrm{T}_{c(t)} \mathcal{M}} \ddot c(t)\).

3.3.0.1 A second-order retraction: metric projection.

For Riemannian submanifolds of a Euclidean space, a natural way to define the retraction consists in taking a step in the embedding space and to project back onto the manifold. This is the metric projection retraction, and it is second order [63], [60]. More precisely, given a point \((X, P) \in \mathcal{M}\) and a tangent vector \((\dot{X}, \dot{P}) \in \mathrm{T}_{(X, P)}\mathcal{M}\), metric projection is \[\begin{align} \label{eq:proj-retr} \mathrm{R}_{(X, P)}(\dot{X}, \dot{P}) &\in \mathop{\mathrm{arg\,min}}_{(\bar X, \bar P) \in \mathcal{M}} \big\|\big(X + \dot{X} - \bar X, P + \dot{P} - \bar P\big)\big\|^2, \end{align}\tag{7}\] where the norm is induced by 5 . (The solution is unique for small tangent vectors but not in general: more on this below.) Algorithm 2 describes an efficient procedure with complexity \(O((m + n)r^2)\) to compute a representation \((\bar U, \bar \Sigma, \bar V)\) for a solution to 7 . We justify this with a lemma then a proposition.

The metric projection retraction is not uniquely defined everywhere (and hence is not smooth everywhere). To see this, consider a point with representation \((U, \Sigma, V)\) and a tangent vector with representation \((K, V_\perp L)\). Choose \(\Sigma = I_r\), \(K = -(2 \alpha + 1) U\) and \(L = (4 \alpha^2 + 2 \alpha)^{\frac{1}{2}} J\), where \(J \in {\mathbb{R}}^{(n - r) \times r}\) contains the first \(r\) columns of \(I_{n - r}\) (assuming \(r \leq n - r\)). Then the matrix \(C\) defined in the proof of Proposition 15 is \[\begin{align} C = 2\alpha(2 \alpha + 1) \begin{bmatrix} V & V_\perp \end{bmatrix} \begin{bmatrix} I_{2r} & 0\\ 0& 0\end{bmatrix} \begin{bmatrix} V & V_\perp \end{bmatrix}^\top. \end{align}\] The top \(2r\) eigenvalues of \(C\) are the same, so the metric projection is not unique (see Lemma 4). This motivates the next part.

3.3.0.2 A globally smooth second-order retraction.

We proceed to derive a computable second-order retraction that is smoothly defined on the whole tangent bundle. It is based on the polar retraction for the Grassmann manifold [60] (as opposed to the Q-factor example in Section 2.2). The proof is provided in Appendix 9.

In Proposition 16, the matrix \(Z \in {\mathbb{R}}^{n \times r}\) always has rank \(r\), so the retraction is well-defined on all of \(\mathrm{T}\mathcal{M}\) (unlike the metric projection retraction). It is also possible to compute a representation of the retracted point in \(O((m + n)r^2)\) operations. Indeed, building the matrices \(Z\) and \(Q\) takes \(O((m + n)r^2)\) operations. Then building the matrix \(W = (X + \dot{X}) Q = (U \Sigma + K) V^\top\mkern-1muQ + (U \Sigma) V_p^\top\mkern-1muQ\) takes \(O((m + n)r^2)\) operations. Now let \(W = \bar U \bar \Sigma H^\top \in {\mathbb{R}}^{m \times r}\) be a thin SVD, which can be computed in \(O(mr^2)\) operations. Then \((\bar X, \bar P) = \mathrm{R}_{(X, P)}(\dot{X}, \dot{P})\) has representation \((\bar U, \bar \Sigma, \bar V)\), where \(\bar V = Q H\).

4 Convergence guarantees↩︎

In this section we use the desingularization to obtain convergence guarantees for low-rank optimization. Again, \(f\colon {\mathbb{R}}^{m \times n} \to {\mathbb{R}}\) is the function we want to minimize over \(\mathbb{R}_{\leq r}^{m \times n}\)—that is problem 1 —and \(g= f\circ \varphi\) is the lifted function over \(\mathcal{M}\). The reparameterized optimization problem is hence \[\begin{align} \label{eq:lifted-problem} \min_{(X, P) \in \mathcal{M}}\, g(X, P) && \textrm{ with } && g(X, P) = f(X). \end{align}\tag{8}\]

We first establish a global convergence theorem. Specifically, we derive algorithms whose iterates accumulate only at stationary points of 1 . This is not straightforward because the feasible set \(\mathbb{R}_{\leq r}^{m \times n}\) is nonsmooth. Consequently, a sequence of points in \(\mathbb{R}_{\leq r}^{m \times n}\) may converge to a non-stationary point, even when the stationarity measure of the sequence converges to zero. [12] exhibited such sequences generated by the algorithm P\(^2\)GD [6]. They also designed second-order algorithms that provably accumulate only at stationary points of 1 using reparameterizations through various \(\varphi\) (including the desingularization). [40], [41], [45], [46], and [42] obtained guarantees for first-order algorithms, including projected gradient descent.

We propose returning to the reparameterization approach, here with the desingularization, because it allows for second-order algorithms. As we detail later, these methods further enjoy fast local convergence.

The global convergence theorem below is a variant of [12]. The argument is the same, but here it is applied to any descent algorithm, that is, one that generates sequences of iterates with non-increasing function values.

Many standard algorithms satisfy the descent condition, such as gradient descent (with line search or appropriate constant step size) and trust-region methods. Except in some pathological cases, they also accumulate only at second-order critical points.

We additionally aim for fast local convergence. In general, the minimizers of 8 are not isolated. This is because the fiber \(\varphi^{-1}(X)\) of a local minimizer \(X\) of 1 is a continuum of local minimizers of 8 whenever \(\mathop{\mathrm{rank}}X < r\) (see Proposition 6). Consequently, we cannot assume local strong convexity. In the literature, it is typical to resort instead to the Polyak–Łojasiewicz (\(\text{P\L}\)) condition (see Definition 5 below), which requires the squared gradient norm to locally dominate the relative function value around local minimizers. In fact, \(\text{P\L}\) is necessary for fast local convergence of gradient descent [64].

Accordingly, we want to obtain \(\text{P\L}\) for \(g\) based on hypotheses on \(f\). This is nontrivial because the set \(\mathbb{R}_{\leq r}^{m \times n}\) is nonsmooth. The \(\text{P\L}\) condition is not even well defined for 1 and we must assume something different.

For this reason, we rely on the Morse–Bott (\(\text{MB}\)) property. It is locally equivalent to \(\text{P\L}\) on manifolds and requires positive definiteness of the Hessian along the normal space to the set of minimizers. In Section 4.3 we define the variant 12 , a type of \(\text{MB}\) condition for \(f\) that takes into account the nonsmoothness of \(\mathbb{R}_{\leq r}^{m \times n}\), allowing us to deduce \(\text{P\L}\) for \(g\).

Our results encompass existing local convergence theory for minimizers of maximal rank. They also enable fast convergence to minimizers of non-maximal rank. As an illustration, if \(\nabla^2f(X) \succ 0\) at a minimizer \(X \in \mathbb{R}_{\leq r}^{m \times n}\) of 1 then \(f\) satisfies 12 at \(X\). In this simple case, the minimizer is (locally) unique, but the definition of 12 allows for more complicated sets of minimizers. See also Remark 25.

Many standard algorithms have fast local rates when \(\text{P\L}\) holds. This includes gradient descent (linear), regularized Newton (quadratic) and trust-regions (superlinear). See [65], [66] for gradient descent and regularized Newton in the Euclidean case, and [10], [67] for trust-regions and extensions to Riemannian manifolds as is needed here.

4.1 Points of interest of the two problems↩︎

To study convergence guarantees we must understand how the points of interest of 1 and 8 relate. The parameterization \(\varphi\) may indeed induce spurious critical points, that is, points \((X, P) \in \mathcal{M}\) that are critical for 8 but such that \(\varphi(X, P) = X\) is not stationary for 1 . However, there exist some (already known) relations that we briefly mention here.

If \(X \in \mathbb{R}_{\leq r}^{m \times n}\) is a stationary point (respectively a local minimizer) for 1 then the fiber \(\varphi^{-1}(X)\) is composed of critical points (respectively local minimizers) for 8 : this is so because \(\varphi\) is smooth.

Conversely, a critical point (respectively a local minimizer) for 8 may not be a stationary point (respectively a local minimizer) for 1 . (This is partly because \(\varphi\) is not open [43].) However, the following known facts provide welcome guarantees:

• Let \((X, P) \in \mathcal{M}\) with \(\mathop{\mathrm{rank}}X = r\). If \((X, P)\) is a critical point (respectively a second-order critical point, a local minimizer) for 8 then \(X\) is a stationary point (respectively a second-order stationary point, a local minimizer) for 1 .

• If \((X, P)\) is second-order critical for 8 then \(X\) is stationary for 1 .

See [43] and [44] for proofs. These facts are analogous to the ones established by [47] for the LR parameterization. In some cases, these properties can be used to deduce the strict saddle property for the LR parameterization, and therefore we expect the same should be true for the desingularization.

4.2 Global convergence↩︎

This section provides the ingredients at the basis of Theorem 17.

4.2.0.1 Sublevel sets.

A notable property of the parameterization \(\varphi\) is that it preserves compactness of sublevel sets. This is not the case for the parameterization \((L, R) \mapsto LR^\top\) because its fibers are unbounded.

We deduce Theorem 17 from this and standard compactness arguments. Indeed, iterates of descent algorithms stay in the initial sublevel set. If it is compact, they also accumulate to at least one point. Section 4.1 further provides the stationarity guarantees. This type of argument appears in [12] for trust-region methods in the context of low-rank optimization. We cannot say the same for the LR parameterization because it is a priori not possible to ensure that the iterates remain in a bounded domain (unless we modify the algorithm, e.g., with rebalancing).

4.2.0.2 Lipschitz continuity.

Lipschitz properties of the cost function are often instrumental to obtain descent algorithms and to ensure that iterates accumulate at first- or second-order critical points. When \(g\) is \({\mathrm{C}^{2}}\), standard continuity arguments readily give that \(\nabla g\) is Lipschitz continuous in any compact domain. The same is true for the Hessian \(\nabla^2g\) when \(g\) is \({\mathrm{C}^{3}}\). It is possible to bound the Lipschitz constants and gain some control. We now study \(\nabla g\) and \(\nabla^2g\) for this purpose.

In particular, if \(f\) is \(L\)-Lipschitz continuous on \({\mathbb{R}}^{m \times n}\) then \(g\) is \(L\)-Lipschitz continuous on \(\mathcal{M}\). We further obtain a bound for the Hessian.

The bound is still finite when \(\sigma_r(X) = 0\) thanks to the term in \(\alpha\). In contrast, for a \({\mathrm{C}^{2}}\) function on \({\mathbb{R}}_r^{m \times n}\), the Riemannian Hessian can blow up as \(\sigma_r(X)\) approaches zero.

4.3 Local convergence with a \(\text{P\L}\)-type condition↩︎

The Polyak–Łojasiewicz (\(\text{P\L}\)) condition is a standard assumption to derive local convergence guarantees for optimization algorithms [65], [68]–[71].

We seek to obtain the \(\text{P\L}\) condition for \(g\) via appropriate assumptions on \(f\). However, it is a priori not clear what to assume on \(f\) because \(\mathbb{R}_{\leq r}^{m \times n}\) is nonsmooth. In particular, \(\text{P\L}\) is not well defined for \(f\) around points of non-maximal rank.

To deal with these singular points, we resort to the local Morse–Bott property (which is equivalent to \(\text{P\L}\) for \({\mathrm{C}^{2}}\) functions on manifolds), and we propose a version for \(f\) that accommodates the nonsmoothness of \(\mathbb{R}_{\leq r}^{m \times n}\). As a preliminary, we need the notion of tangent cone.

4.3.0.1 Tangent cones.

At a point \(X \in \mathbb{R}_{\leq r}^{m \times n}\), the tangent cone to \(\mathbb{R}_{\leq r}^{m \times n}\) is defined as \[\begin{align} \mathrm{T}_X \mathbb{R}_{\leq r}^{m \times n}= \bigg\{ \lim_{k \to \infty} \frac{X_k - X}{t_k} : X_k \in \mathbb{R}_{\leq r}^{m \times n}, t_k > 0 \text{ for all k},\, t_k \to 0 \text{ and the limit exists}\bigg\}, \end{align}\] which is a closed (not necessarily convex) cone. A simple expression is known for this cone [6]. At a point \(X\) of rank \(s\), we have \[\label{eq:tgt-bounded} \begin{align} \mathrm{T}_X {\mathbb{R}}_{\leq r}^{m \times n} = \bigg\{& \begin{bmatrix} U_s & U_{s, \perp} \end{bmatrix} \begin{bmatrix} A & B\\ C & D \end{bmatrix} \begin{bmatrix} V_s & V_{s, \perp} \end{bmatrix}^\top :\\ &\qquad \text{A, B, C, D are of appropriate sizes and \mathop{\mathrm{rank}}(D) \leq r - s} \bigg\}, \end{align}\tag{10}\] where \(U_s \in \mathrm{St}(m, s)\) and \(V_s \in \mathrm{St}(n, s)\) span the column and row spaces of \(X\) respectively, and \(U_{s, \perp}\) and \(V_{s, \perp}\) are orthonormal completions. The tangent space \(\mathrm{T}_X {\mathbb{R}}_s^{m \times n}\) has the same set description as 10 with the additional constraint that \(D = 0\).

4.3.0.2 The Morse–Bott property.

We now define the Morse–Bott property for \(g\), and a variant for \(f\) tailored to \(\mathbb{R}_{\leq r}^{m \times n}\). It requires \(f\) and \(g\) to be \({\mathrm{C}^{2}}\). We let \(\mathcal{S}_f\) and \(\mathcal{S}_g\) denote the sets that contain all the local minimizers of 1 and 8 respectively.

This condition is equivalent to \(\text{P\L}\) for \({\mathrm{C}^{2}}\) functions, up to arbitrarily small loss in the constant \(\mu\) [10].

It is now natural to ask: under what conditions on \(f\) is it the case that \(g= f\circ \varphi\) satisfies 11 ? We suggest one answer: a Morse–Bott type condition with an additional requirement on \(\mathcal{S}_f\). In Section 4.4, we show the implication \(\text{MB}{}_f\Rightarrow \text{MB}{}_g\) (which could fail without the extra condition on \(\mathcal{S}_f\)).

Note the two differences between 11 and 12 . For the latter, (i) we require the solution set \(\mathcal{S}_f\) to be included in a constant rank stratum \({\mathbb{R}}_s^{m \times n}\), and (ii) the inner product lower bound needs to hold only for matrices \(\dot{X}\) in the tangent cone at \(X\). This second fact makes the condition depend only on the values that \(f\) takes in \(\mathbb{R}_{\leq r}^{m \times n}\) (as opposed to depending also on the values of \(f\) in a neighborhood of \(\mathbb{R}_{\leq r}^{m \times n}\) in \({\mathbb{R}}^{m \times n}\)). More explicitly, if two functions \(f_1\) and \(f_2\) coincide on \(\mathbb{R}_{\leq r}^{m \times n}\) then the quantities \(\inner{\dot{X}}{\nabla^2f_i(X)[\dot{X}]}\) are the same when \(\dot{X} \in \mathrm{T}_X \mathbb{R}_{\leq r}^{m \times n}\). (This can be seen using that \(\mathbb{R}_{\leq r}^{m \times n}\) is geometrically derivable [72].) This is arguably a desirable property, as it makes our assumption intrinsic.

4.4 Transitivity of a \(\text{P\L}\)-type condition↩︎

In this section we prove that 12 implies 11 . Theorem 18 is then a direct consequence. We study the maximal rank stratum and lower-rank strata separately, with different arguments.

4.4.0.1 Maximal rank stratum.

We first examine \(\text{P\L}\) on constant rank strata and deduce \(\text{P\L}\) transitivity for the maximal rank stratum. In this part we only require \(f\) (and hence \(g\)) to be continuously differentiable (\({\mathrm{C}^{1}}\)). We derive below an expression for the tangent spaces of constant rank strata, which we need for the analysis.

It is informative to particularize Proposition 22 to \(s = r\). Indeed, the maximal-rank stratum of \(\mathcal{M}\) is open (see Remark 8), so the restriction \(g_r\) coincides with \(g\) on \(\varphi^{-1}({\mathbb{R}}_r^{m \times n})\); in particular, their gradients (hence gradient norms) agree there. Assuming that \(f_r\) satisfies 9 around a point \(X \in {\mathbb{R}}_r^{m \times n}\), Proposition 22 implies that \(g_r\) satisfies 9 around \(\varphi^{-1}(X)\), and therefore so does \(g\) on the same neighborhood. However, this reasoning fails for lower-rank strata because they are not open, so we handle those next by a different argument.

4.4.0.2 Non-maximal rank strata.

For lower-rank strata, we work with the Morse–Bott property, as defined in Section 4.3. We first prove a technical lemma and then state the main result of this section.

Theorem 18 is then a direct consequence.

When \(g\) satisfies \(\text{MB}\) (and equivalently \(\text{P\L}\)), we can deduce local convergence results for standard algorithms, as described below Theorem 18.

5 Numerical experiments↩︎

In this section we run matrix completion experiments to compare the LR parameterization, fixed-rank optimization, and the desingularization (with three metric parameters: \(\alpha = 1/20\), \(\alpha = 1/2\) and \(\alpha = 5\)). Matrix completion does not satisfy the assumptions of Section 4—namely, compact sublevel sets and the Morse–Bott condition of Theorems 17 and 18 (see Remark 25). Nevertheless, it remains an important problem in low-rank optimization, which motivates its use here for illustration.

Given \(A \in {\mathbb{R}}^{m \times n}\) and a mask \(\Omega \in \{0, 1\}^{m \times n}\), the cost function \(f\colon {\mathbb{R}}^{m \times n} \to {\mathbb{R}}\) is defined as \(f(X) = \frac{1}{2}\|(X - A) \odot \Omega\|_{\mathrm{F}}^2\), where \(\odot\) denotes the Hadamard product. The (Euclidean) gradient and Hessian are \(\nabla f(X) = (X - A) \odot \Omega\) and \(\nabla^2f(X)[\dot{X}] = \dot{X} \odot \Omega\). Both are sparse when \(\Omega\) is sparse: this is exploited in the code.

We set the dimensions as \(m = n = 5000\) in all the experiments. The rank \(r^*\) of \(A\) is comparatively small (see details below). We construct \(A = U_A^{} \Sigma_A^{} V_A^\top\) for some matrix of singular values \(\Sigma_A \in {\mathbb{R}}^{r^* \times r^*}\), and where \(U_A \in \mathrm{St}(m, r^*)\) and \(V_A \in \mathrm{St}(n, r^*)\) are sampled uniformly at random. The mask \(\Omega\) is sparse with an oversampling factor of \(5\), meaning that it contains \(5(m + n - r) r\) non-zero entries. They are chosen uniformly at random.

We interface the implementation with the Matlab toolbox Manopt [13]. For all geometries, we run a (Riemannian) trust-region algorithm [73] with the truncated conjugate gradient subproblem solver.⁵ Evaluating the cost function, the gradient, and Hessian vector products takes similar times for the three geometries. For a given experiment, all the algorithms start from the same initial random point in \(\mathbb{R}_{\leq r}^{m \times n}\) (with balanced factors for the LR parameterization). We pick the initialization close to zero (compared to \(A\)) because we found that it significantly improves performance for all. Specifically, we sample \(U \in \mathrm{St}(m, r), V \in \mathrm{St}(n, r), \sigma_1, \dots, \sigma_r \in \interval{0}{10^{-3}}\) all independently and uniformly at random, and then set the initial point as \(U \mathop{\mathrm{diag}}(\sigma_1, \dots, \sigma_r) V^\top\).

For each experiment we report two figures. Their vertical axes display the cost function value. The horizontal axis of the first figure is the number of outer iterations of trust-regions. The horizontal axis of the second figure is the time in seconds.

5.0.0.1 Overestimation of the rank.

We let the target matrix \(A\) have rank \(r^* = 10\) and optimize with \(r = 20\). The 10 non-zero singular values of \(A\) are sampled uniformly at random in the interval \(\interval{\frac{1}{2}}{1}\). Figure 3 reports the results.

In this scenario the convergence is particularly fast for the desingularization.

5.0.0.2 Exponential decay.

We let \(A\) have rank \(r^* = 20\) with singular values that follow an exponential decay as \(\sigma_i(A) = 0.9^{i - 1}\) for \(i = 1, \dots, r^*\). Figure 4 displays the results with an optimization rank \(r = r^* = 20\).

In this setting both LR and fixed-rank optimization are faster than the desingularization. Figure 5 reports the results for the exact same experiment but with an optimization rank of \(r = 30 > r^*\).

This is detrimental for LR and fixed-rank optimization but the desingularization is robust to this rank overestimation.

5.0.0.3 Other solvers for matrix completion.

This paper deals with general algorithms for bounded-rank optimization. The numerical experiments compare general algorithms on a specific task, namely, matrix completion. There exist many solvers tailored for the latter. For example, the algorithms proposed in [35], [36], [38] exhibit state-of-the-art performances in practice. We experimented with GNMR [38] and empirically observed the following. This algorithm is particularly efficient to recover the matrix \(A\), even in hard settings where the oversampling factor is close to \(1\), or when the condition number is large. However, iterations are expensive and we found that GNMR is slower than the trust-region method. We do not include these results for this reason.

6 Perspectives↩︎

We close with a few directions for future research.

• The desingularization is particularly relevant when the true (or numerical) rank \(r^*\) of the target matrix is strictly less than the optimization rank \(r\). This is especially convenient when \(r^*\) is not known because it allows the user to set \(r\) as a conservative upper bound of \(r^*\). In that scenario, we may want to nudge the algorithm towards solutions of rank less than \(r\). This could be achieved with a nuclear norm penalty [27], thus mixing both a hard and a soft low-rank prior.

• Riemannian optimization theory breaks for algorithms running over \({\mathbb{R}}_r^{m \times n}\) when the sequence converges to a point of rank strictly less than \(r\) (because it lies outside of the manifold). Notwithstanding, even in that setting, we sometimes observe fast local convergence. [50] explain that behavior for gradient descent on specific cost functions. In order to generalize such results, it may be possible to use the desingularization with the family of Riemannian metrics we propose, and to study algorithms in the limit where the metric parameter \(\alpha\) goes to zero. Indeed, \({\mathbb{R}}_r^{m \times n}\) (with its usual metric) is isometric to the maximal rank stratum of \(\mathcal{M}\) when \(\alpha = 0\) (see Remark 12).

• Our local convergence results focus on the Polyak–Łojasiewicz condition. They may generalize to the broader class of Łojasiewicz inequalities. Especially, one might ask: does \(g\) satisfy a Łojasiewicz inequality if \(f\) satisfies an adequately modified version of that inequality for \(\mathbb{R}_{\leq r}^{m \times n}\)? A starting point could be [6].

Acknowledgments↩︎

We thank Yuetian Luo and Anru Zhang for offering detailed pointers to relevant literature, as well as Eitan Levin, Joe Kileel, Bart Vandereycken and Guillaume Olikier for numerous discussions.

7 A note about regularizing the LR parameterization↩︎

It is common to regularize the LR parameterization. The cost function becomes \(g(L, R) = f(LR^\top) + \lambda \rho(L, R)\), where \(f\) is the function to minimize over \(\mathbb{R}_{\leq r}^{m \times n}\), \(\rho \colon {\mathbb{R}}^{m \times r} \times {\mathbb{R}}^{n \times r} \to {\mathbb{R}}\) is the regularizer, and \(\lambda > 0\) is the regularization intensity.

One possible choice is \(\rho(L, R) = \|L\|_\mathrm{F}^2 + \|R\|_\mathrm{F}^2\). This promotes solutions with small nuclear norm [74]. While this is useful in several contexts, we do not entertain it here because it changes the minimizers of the cost function.

As introduced in [75], the regularizer \(\rho(L, R) = \|L^\top L - R^\top R\|_\mathrm{F}^2\) promotes balanced factors. It does not change the global minimizers of the problem because the LR factors of a solution can always be balanced. This regularizer has other nice properties, as described in [76]. However, in general the regularized cost \(g\) does not satisfy \(\text{P\L}\) at points of non-maximal rank. This can happen even for simple costs such as \[\begin{align} (L, R) \mapsto \frac{1}{2}\|LR^\top - A\|_\mathrm{F}^2 + \frac{\lambda}{4}\|L^\top L - R^\top R\|_\mathrm{F}^2 \end{align}\] with \(\mathop{\mathrm{rank}}A < r\). For example, when \(m = n = r = 1\) and \(A = 0\), the cost becomes \((x, y) \mapsto \frac{1}{2}x^2 y^2 + \frac{\lambda}{4}(x^2 - y^2)^2\), and the Hessian is zero at \((0, 0)\).

8 Second-order derivatives↩︎

The expression of the Hessian involves the derivative of the projector onto tangent spaces (see [77], [60]). To find it we define four auxiliary functions: \[\label{eq:qs} \begin{align} &\xi_0(X) = (2\alpha I + XX^\top)^{-1}, && && \xi_1(X) = XX^\top\mkern-1mu\xi_0(X) = I - 2 \alpha \xi_0(X),\\ &\xi_2(X) = \xi_0(X)X, && && \xi_3(X) = (2\alpha I + X^\top\mkern-1muX)^{-1}. \end{align}\tag{19}\] With \(X = U \Sigma V^\top\) an SVD, the definition of \(S(\alpha)\) 6 and direct computations give \[\label{eq:qs-usv} \begin{align} &\xi_0(X) = U S(\alpha)^{-1} U^\top + \frac{1}{2\alpha} U_\perp^{} U_\perp^\top, && && \xi_1(X) = U\Sigma^2 S(\alpha)^{-1} U^\top,\\ &\xi_2(X) = U S(\alpha)^{-1} \Sigma V^\top, && && \xi_3(X) = V S(\alpha)^{-1} V^\top + \frac{1}{2\alpha} P. \end{align}\tag{20}\] Let \(\mathop{\mathrm{sym}}(A) = (A + A^\top)/2\) denote the symmetric part of a matrix \(A\). Using the functions defined in 19 , we first express the orthogonal projector (Proposition 11) as a function of the point \((X, P)\) rather than its representation \((U, \Sigma, V)\).

We can now scrutinize the differential of the orthogonal projection. Given \((X, P) \in \mathcal{M}\) and \((\dot{X}, \dot{P}) \in \mathrm{T}_{(X, P)}\mathcal{M}\), we define \(\mathcal{P}\colon \mathcal{E}\to \mathcal{E}\) as \[\begin{align} \label{eq:def-projdiff} \mathcal{P}(Y, Z) = \mathrm{D}\!\left( (X, P) \mapsto \mathop{\mathrm{proj}}_{(X, P)}(Y, Z)\right)\!(X, P)[\dot{X}, \dot{P}]. \end{align}\tag{22}\] We explicitly compute this object because it plays a role in the Hessian of \(g\).

We are now ready to express the Riemannian Hessian \(\nabla^2g\) as a function of \(\nabla f\) and \(\nabla^2f\).

9 Second-order retractions↩︎

The lemma below characterizes when the retractions given in Proposition 4 are second-order (see Definition 4).

References↩︎

---

1. Ecole Polytechnique Fédérale de Lausanne (EPFL), Institute of Mathematics. This work was supported by the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract number MB22.00027.↩︎

2. One might also worry about \(g\) suffering from spurious critical points introduced by the parameterization: see Section 4.1 for reassurances.↩︎

3. Equivalently, we could define the Grassmannian as a quotient (Stiefel manifold by orthogonal group): see [11] for a review. That is the approach of [1].↩︎

4. https://github.com/qrebjock/desingularization-experiments↩︎

5. We set the power for the residual stopping criterion [73] as \(\theta = \sqrt{2} - 1\) (instead of the usual default \(\theta = 1\)). Indeed, it appears that setting \(\theta < 1\) helps convergence for all the methods.↩︎