Low-Rank Matrix Recovery with Scaled Subgradient Methods: Fast and Robust Convergence Without the Condition Number

Introduction

using first-order methods (e.g. gradient descent). While tremendous progress has been made in recent years [CLC19], applying vanilla gradient descent to the above smooth formulation suffers from two sources of ill-conditioning that preclude a desirable computational efficiency from classical optimization principles:

Due to the heavy-tailed nature of certain measurement operators, such as those encountered in phase retrieval [CLS15] and quadratic sampling [SWW17], the least-squares formulation (2) may suffer from a large smoothness parameter (and hence a large condition number of the loss function) that scales at least linearly with respect to the ambient dimension, leading to a conservative choice of stepsizes and a high iteration complexity when the problem dimension is large.

Due to the composite nature of the formulation (2), the iteration complexity of vanilla gradient descent is further exacerbated by the condition number of the underlying low-rank matrix $\bm{X}_{\star}$, which could be large in many applications of interest.

While there have been encouraging activities [CCD+21, MWCC19, LMCC21, TMC20] that try to alleviate these issues regarding ill-conditioning, none of the existing first-order approaches are able to simultaneously remove both sources of ill-conditioning and achieve fast convergence. Therefore, the goal of the current paper is to develop first-order methods that are guaranteed to converge at a fast rate that is almost dimension-free and independent of the condition number, even in the presence of corruptions.

In this paper, we propose to minimize the following nonsmooth and nonconvex loss function known as the least absolute deviations, which measures the residual sum of absolute errors

Here, $\bm{S}_{t}\in\partial f(\bm{L}_{t}\bm{R}_{t}^{\top})$ is a subgradient of $f(\bm{X})\coloneqq\sum_{i=1}^{m}\left|\mathcal{A}_{i}(\bm{X})-y_{i}\right|$ at $\bm{L}_{t}\bm{R}_{t}^{\top}$, and $\eta_{t}>0$ is a sequence of carefully-chosen stepsizes. Compared with vanilla subgradient methods, our new method (4) scales or preconditions the search directions $\bm{S}_{t}\bm{R}_{t}$ and $\bm{S}_{t}^{\top}\bm{L}_{t}$ by $(\bm{R}_{t}^{\top}\bm{R}_{t})^{-1}$ and $(\bm{L}_{t}^{\top}\bm{L}_{t})^{-1}$, respectively.Under appropriate conditions, the inverse matrices always exist; in practice, one can use the pseudo-inverse matrices to avoid numerical instabilities. As explained in [TMC20] where a similar preconditioning trick was employed for smooth formulations, the scaled subgradient enables better search directions and therefore larger stepsizes. Our main results can be summarized as follows:

Under general geometric assumptions on $f(\cdot)$ such as restricted rank-$r$ Lipschitz continuity and sharpness conditions, we demonstrate that the convergence rate of scaled subgradient methods using both Polyak’s and geometrically decaying stepsizes is independent of the condition number of $\bm{X}_{\star}$.

Instantiating our theory under the mixed-norm restricted isometry property (RIP) of the measurement operator, we demonstrate state-of-the-art computational guarantees for low-rank matrix sensing and quadratic sampling even when the observations are noisy and corrupted by outliers. This leads to improvements over the computational complexity of scaled gradient methods in [TMC20] for heavy-tailed measurement ensembles, as well as of vanilla subgradient methods in [CCD+21]. Table 1 provides a detailed comparison of the local iteration complexities of the proposed scaled subgradient method in comparison with these prior algorithms, highlighting its robustness to heavy-tailed observations, outliers, as well as a large condition number of the true matrix $\bm{X}_{\star}$.

Our work leverages exciting advances in nonsmooth optimization [CCD+21] and scaled first-order methods [TMC20] for low-rank matrix recovery. Our arguments are concise, which avoid the need of sophisticated trajectory-dependent analysis as have been used in [MWCC19, LMCC21] to achieve rapid and robust convergence guarantees.

2 Related work

Low-rank matrix recovery has been a target of intense interest in the last decade; we invite the readers to [DR16, CC18, CLC19] for recent overviews, and limit our discussions to the most relevant literature in the sequel.

Nonsmooth objective functions, such as the least absolute deviations, have been adopted earlier in both convex and nonconvex formulations of low-rank matrix recovery, including phase retrieval [Han17, DDP17, QZEW17, ZZLC17, DR19], blind deconvolution [Día19], quadratic sampling [LSC17, CL16, CCD+21, BL20], low-rank matrix sensing [CCD+21, Li13, WGMM13, LZSV20], robust synchronization [WS13], to name a few. Our work is most closely related to and generalizes the vanilla subgradient method in [CCD+21], by establishing novel performance guarantees of scaled subgradient methods for robust low-rank matrix recovery.

Scaled first-order methods for low-rank matrix recovery.

Variants of the scaled gradient methods are proposed in [MAS12, TW16, TMC20] for minimizing the least-squares formulation (2), where strong statistical and computational complexities are first established in [TMC20]. To the best of our knowledge, the current paper is the first work that provides rigorous statistical and computational guarantees for scaled subgradient methods for addressing nonsmooth formulations. When it comes to problems with heavy-tail observations such as quadratic sampling, while it is possible to establish faster convergence rates of vanilla gradient descent over the smooth least-squares loss function through a tailored analysis [MWCC19, LMCC21] via leave-one-out arguments, it is unclear if similar analyses are viable for scaled gradient methods (ScaledGD) in [TMC20]. Unfortunately, a direct application of the performance guarantee of ScaledGD on minimizing the smooth least-squares loss function leads to a much slower rate in terms of the problem dimension (see Table 1) for quadratic sampling. In contrast, our analysis for scaled subgradient methods yields strong guarantees in a more straightforward manner since the nonsmooth loss function has much better geometric properties [CCD+21].

Robust low-rank matrix recovery via nonconvex optimization.

A pleasant side benefit of nonsmooth formulations is the added robustness to adversarial outliers under a simple algorithm design – the low-rank factors are updated essentially in the same manner regardless of the presence of outliers. In comparison, other nonconvex methods based on smooth formulations often need to introduce some special treatments to mitigate outliers before updating the low-rank factors, e.g. truncation or thresholding [ZCL16, LCZL20, LZSV20], which can be cumbersome to tune properly.

Condition number independent rate of convergence.

It is well-known that first-order methods such as gradient descent exhibit poor scaling with respect to the condition number of the low-rank matrix. Possible remedies include alternating least-squares in the factored space [JNS13, HW14], or spectral methods over the matrix space [JMD10]. However, these approaches either require the inversion of a large matrix or a higher memory footprint, compared with the scaled first-order methods adopted herein.

3 Paper organization and notation

The rest of this paper is organized as follows. Section 2 describes the proposed scaled subgradient method and its connections to existing methods. Section 3 provides the theoretical guarantees for the scaled subgradient method in terms of both statistical and computational complexities, which are then instantiated to robust low-rank matrix sensing and quadratic sampling. Section 4 illustrates the superior empirical performance of the proposed method. Finally, we conclude in Section 5. The proofs are deferred to the appendix.

Problem Formulation and Algorithms

In this section, we formulate the low-rank matrix recovery problem, followed by a detailed description of the proposed scaled subgradient method.

Without loss of generality, we define the ground truth low-rank factors as

so that $\bm{X}_{\star}=\bm{L}_{\star}\bm{R}_{\star}^{\top}$. Moreover, we denote the ground truth stacked factor matrix as

Assume that we have access to a number of observations $\bm{y}=\{y_{i}\}_{i=1}^{m}$ of $\bm{X}_{\star}$, given as

where $\mathcal{A}(\bm{X}_{\star})=\{\mathcal{A}_{i}(\bm{X}_{\star})\}_{i=1}^{m}$ is the measurement ensemble, $\bm{w}=\{w_{i}\}_{i=1}^{m}$ denotes the bounded noise, and $\bm{s}=\{s_{i}\}_{i=1}^{m}$ models arbitrary corruptions. The goal of low-rank matrix recovery is to reconstruct $\bm{X}_{\star}$ from the noisy and corrupted observations $\bm{y}$ in a statistically and computationally efficient manner.

2 Scaled subgradient method

Consider the following nonsmooth and nonconvex optimization problem over the factors

where $f(\cdot)$ is a nonsmooth surrogate of the observation residuals. Of particular interest is the residual sum of absolute errors, defined as

Correspondingly, the minimizer is called the least absolute deviations (LAD) solution.

Let us denote the stacked factor matrix in the $t$-th iterate as $\bm{F}_{t}\coloneqq[\bm{L}_{t}^{\top},\bm{R}_{t}^{\top}]^{\top}$. Given an initialization $\bm{F}_{0}=[\bm{L}_{0}^{\top},\bm{R}_{0}^{\top}]^{\top}$, the proposed scaled subgradient method (ScaledSM) proceeds as

where $\bm{S}_{t}\in\partial f(\bm{L}_{t}\bm{R}_{t}^{\top})$ is a subgradient of $f(\cdot)$ at $\bm{L}_{t}\bm{R}_{t}^{\top}$ (and hence $\bm{S}_{t}\bm{R}_{t}\in\partial_{\bm{L}}\mathcal{L}(\bm{L}_{t},\bm{R}_{t})$ and $\bm{S}_{t}^{\top}\bm{L}_{t}\in\partial_{\bm{R}}\mathcal{L}(\bm{L}_{t},\bm{R}_{t})$), and $\eta_{t}>0$ is some properly selected stepsize, which we discuss next.

We consider the following two choices of stepsize schedules:

If we know the optimal value $f(\bm{X}_{\star})$, we can invoke the following Polyak’s stepsize, given by

where the denominator is the squared norm of the subgradient under a scaled metric concerted with the preconditioners. This schedule is implementable, for example, when the observations are noise-free, leading to $f(\bm{X}_{\star})=0$. However, when the observations are noisy and corrupted, it is not viable to know $f(\bm{X}_{\star})$ beforehand.

In general, we can apply the geometrically decaying stepsize originally introduced in [Gof77], given by

where the denominator is similarly scaled as (14), and $\lambda>0$ and $q\in(0,1)$ are some parameters to be specified. This choice is broadly applicable when dealing with noisy and corrupted observations.

Compared with the vanilla subgradient method, which proceeds according to

the update rule (13) scales the subgradient $\bm{S}_{t}\bm{R}_{t}$ and $\bm{S}_{t}^{\top}\bm{L}_{t}$ by $(\bm{R}_{t}^{\top}\bm{R}_{t})^{-1}$ and $(\bm{L}_{t}^{\top}\bm{L}_{t})^{-1}$, respectively; see [TMC20] for its counterpart in smooth problems. An important highlight of the scaled subgradient method is that the update rule is covariant with respect to the ambiguity of low-rank matrix factorization. To see this, imagine that we modify the $t$-th updates as

both the Polyak’s stepsize (14) and the geometrically decaying stepsize (15) do not change, since

which holds similarly for $\|\bm{S}_{t}^{\top}\bm{L}_{t}(\bm{L}_{t}^{\top}\bm{L}_{t})^{-1/2}\|_{\operatorname{\mathsf{F}}}^{2}$;

The next $(t+1)$-th iterate can be written as

and similarly $\widetilde{\bm{R}}_{t+1}=\bm{R}_{t+1}\bm{Q}^{-\top}$. Therefore, all the iterates are covariant with respect to the invertible transform (17).

where $\mathcal{A}^{*}(\cdot)$ is the adjoint operator of $\mathcal{A}(\cdot)$, and $\bm{r}_{t}\coloneqq\mathcal{A}(\bm{L}_{t}\bm{R}_{t}^{\top})-\bm{y}$ is the residual using the $t$-th iterate. Consequently, the scaled subgradient method follows the update rule

where $\bm{S}_{t}\in\partial f(\bm{L}_{t}\bm{L}_{t}^{\top})$ is a subgradient of $f(\cdot)$ at $\bm{L}_{t}\bm{L}_{t}^{\top}$. Our theory applies to this PSD case in a straightforward manner.

Theoretical Guarantees

In this section, we first provide the theoretical guarantees of the scaled subgradient method under general geometric assumptions on $f(\cdot)$, and then instantiate them to concrete problems including robust low-rank matrix sensing and quadratic sampling.

We start by introducing the following geometric properties of the loss function $f(\cdot)$, which play a key role in the convergence analysis.

The first condition is similar to the usual Lipschitz property of a function.

The second geometric condition is akin to the (one-point) strong convexity of a function, with the key difference that strong convexity adopts the squared Euclidean norm whereas the following one uses the plain Euclidean norm.

For notational simplicity, if a function $f$ is both restricted $L$-Lipschitz continuous and $\mu$-sharp, we denote

In some cases, e.g. in the presence of noise, the loss function $f(\cdot)$ only satisfies an approximate restricted sharpness property, which is detailed below.

As shall be seen in Section 3.3, these conditions can be ensured for proper choices of the loss function as long as the observation operator $\mathcal{A}(\cdot)$ satisfies certain mixed-norm RIP, which holds for a wide number of practical problems.

2 Main results

Motivated by [TMC20], we measure the performance of $\bm{F}=[\bm{L}^{\top},\bm{R}^{\top}]^{\top}$ using the following error metric

which takes into consideration both the representational ambiguity of the factorization up to invertible transforms and the scaling effect of preconditioners. In comparison, the more standard distance metric [MLC21] in the analysis of vanilla gradient methods reads as follows

which is inadequate to delineate the power of preconditioning. See [TMC20] for more discussions.

We start with stating the linear convergence of the scaled subgradient method when $f(\cdot)$ satisfies both the rank-$r$ restricted $L$-Lipschitz continuity and $\mu$-sharpness. The proof is deferred to Appendix B.

and the scaled subgradient method in (13) adopts either Polyak’s stepsizes in (14) or geometrically decaying stepsizes in (15) with $\lambda=\sqrt{\frac{\sqrt{2}-1}{2}}0.02\sigma_{r}(\bm{X}_{\star})/\chi_{f}^{2}$ and $q=\sqrt{1-0.16/\chi_{f}^{2}}$. Then for all $t\geq 0$, the iterates satisfy

Theorem 1 shows that the iterates of the scaled subgradient method converges at a linear rate; to reach $\epsilon$-accuracy, i.e. $\|\bm{L}_{t}\bm{R}_{t}^{\top}-\bm{X}_{\star}\|_{\operatorname{\mathsf{F}}}\leq\epsilon\sigma_{r}(\bm{X}_{\star})$, it takes at most $O(\chi_{f}^{2}\log\frac{1}{\epsilon})$ iterations, which, importantly, is independent of the condition number $\kappa$ of $\bm{X}_{\star}$. In addition, it is still possible to ensure approximate reconstruction when only the approximate restricted sharpness property holds, as shown in the next theorem. Again, we postpone the proof to Appendix C.

Theorem 2 shows that as long as the relaxation parameter $\xi$ is sufficiently small, i.e. $\xi\lesssim\sigma_{r}(\bm{X}_{\star})\mu/\chi_{f}$, then the scaled subgradient method with geometrically decaying stepsizes converges at a linear rate until an error floor is hit. In particular, the iterates satisfy $\|\bm{L}_{t}\bm{R}_{t}^{\top}-\bm{X}_{\star}\|_{\operatorname{\mathsf{F}}}\leq 30\xi/\mu$ after at most $O(\chi_{f}^{2})$ iterations up to logarithmic factors.

For simplicity of exposition, we have fixed the values of $\lambda$ and $q$ for the geometrically decaying stepsizes in the above theorems. It is possible to allow a wider range of $\lambda$ and $q$ by slightly modifying the arguments without sacrificing the linear convergence. In practice, these parameters should be tuned in order to yield optimal performance.

3 A case study: robust low-rank matrix recovery

We now apply the above theorems to robust low-rank matrix recovery, which showcases the superior performance of the scaled subgradient method.

We start with the observation model (10) with clean measurements, i.e. $\bm{w}=\bm{0}$ and $\bm{s}=\bm{0}$. To proceed, we assume that the observation operator $\mathcal{A}(\cdot)$ satisfies the following mixed-norm RIP.

The next proposition verifies that the loss function (12) satisfies restricted Lipschitz continuity and sharpness properties under the mixed-norm RIP.

If $\mathcal{A}(\cdot)$ satisfies rank-$2r$ mixed-norm RIP with constants $(\delta_{1},\delta_{2})$, then $f(\bm{X})=\|\mathcal{A}(\bm{X})-\bm{y}\|_{1}=\|\mathcal{A}(\bm{X}-\bm{X}_{\star})\|_{1}$ in (12) satisfies the rank-$r$ restricted $L$-Lipschitz continuity and $\mu$-sharpness with

With the geometric characterization of $f(\cdot)$ in place, we immediately have the following corollary that captures the performance of the scaled subgradient method when $\mathcal{A}(\cdot)$ satisfies the mixed-norm RIP.

Noisy and corrupted case.

We now consider the observation model (10) where the noise $\bm{w}$ is bounded with $\|\bm{w}\|_{1}\leq\sigma_{w}$ and $\|\bm{s}\|_{0}=p_{s}m$, where $p_{s}\in[0,1/2)$ is the fraction of outliers. Following [CCD+21], we further introduce another important property of $\mathcal{A}(\cdot)$.

where $\mathcal{A}_{\mathcal{S}}(\bm{M})=\{\mathcal{A}_{i}(\bm{M})\}_{i\in\mathcal{S}}$ and $\mathcal{A}_{\mathcal{S}^{c}}(\bm{M})=\{\mathcal{A}_{i}(\bm{M})\}_{i\in\mathcal{S}^{c}}$.

The next proposition verifies that the loss function in (12) satisfies restricted Lipschitz continuity and approximate sharpness properties under the mixed-norm RIP (cf. Definition 4) and the $\mathcal{S}$-outlier bound (cf. Definition 5).

Denote the support of the outlier $\bm{s}$ as $\mathcal{S}$. Suppose that $\mathcal{A}(\cdot)$ satisfies rank-$2r$ mixed-norm RIP with $(\delta_{1},\delta_{2})$ and $\mathcal{S}$-outlier bound with $\delta_{3}$, then $f(\bm{X})$ in (12) satisfies rank-$r$ restricted $L$-Lipschitz continuity and $\xi$-approximate $\mu$-sharpness with

Similar to the previous noise-free case, this immediately leads to performance guarantees of the scaled subgradient method when $\mathcal{A}(\cdot)$ satisfies both the mixed-norm RIP and the $\mathcal{S}$-outlier bound.

We now instantiate the above general guarantee to the following two types of observation operators. For simplicity, we assume there is no dense noise, i.e. $\sigma_{w}=0$; see Table 1 for a summary of the comparisons.

matrix sensing: the measurement operator $\mathcal{A}_{i}(\cdot)$ is defined as $\mathcal{A}_{i}(\bm{X}_{\star})=\frac{1}{m}\langle\bm{A}_{i},\bm{X}_{\star}\rangle$, where the matrix $\bm{A}_{i}$ is composed of i.i.d. Gaussian entries $\mathcal{N}(0,1)$.The same guarantee also holds for sub-Gaussian measurements. It is shown in [CCD+21] (see also [LZSV20]) that $\mathcal{A}(\cdot)$ satisfies the mixed-norm RIP and $\mathcal{S}$-outlier bound with

as long as $m\gtrsim\frac{(n_{1}+n_{2})r}{(1-2p_{s})^{2}}\log(\frac{1}{1-2p_{s}})$. Hence, the scaled subgradient method converges linearly to $\epsilon$-accuracy in $O\left(\frac{1}{(1-2p_{s})^{2}}\log\frac{1}{\epsilon}\right)$ iterations provided that it is initialized properly, making it robust simultaneously to ill-conditioning of the matrix $\bm{X}_{\star}$ and the presence of the outliers.

as long as $m\gtrsim\frac{nr^{2}}{(1-2p_{s})^{2}}\log(\frac{\sqrt{r}}{1-2p_{s}})$. Hence, the scaled subgradient method converges linearly to $\epsilon$-accuracy in $O\left(\frac{r}{(1-2p_{s})^{2}}\log\frac{1}{\epsilon}\right)$ iterations, as long as it is seeded with a good initialization. In comparison, the iteration complexity of the scaled gradient descent method over the least-squares loss function depends polynomially with respect to $n$, due to the heavy-tailed nature of the observation operator, let alone its sensitivity to the outliers.

The above discussions are limited to the local iteration complexity, assuming a good initialization satisfying (20) is available. In the absence of outliers, a standard spectral method can be used, as shown in [TMC20]. In the presence of outliers, a truncated spectral method could be used; see e.g. [ZCL16, LCZL20].

Numerical Experiments

In this section, we conduct numerical experiments to corroborate our theory.

Since the vanilla subgradient method (VanillaSM) has been extensively benchmarked against other methods and established as state-of-the-art in [CCD+21], we focus on comparing the proposed scaled subgradient method (ScaledSM) to VanillaSM in the sequel. In general, the geometrically decaying stepsize (15) is a more practical choice than the Polyak’s stepsize (14), especially in the presence of noise and outliers. Nonetheless, using properly tuned geometrically decaying stepsizes essentially matches the performance of using Polyak’s stepsizes, for both VanillaSM [LZSV20] and ScaledSM, the latter of which we shall illustrate in the ensuing experiments. As such, we adopt Polyak’s stepsizes in the comparisons below, to emulate the scenario where both methods are tuned to operate under its largest allowable stepsizes and achieve the fastest convergence. In addition, both algorithms start from the same initialization.

We consider two low-rank matrix estimation tasks discussed in Section 3.3. Recall the observation model in (10) and its entrywise version in (9), which we repeat below for convenience:

Quadratic sampling. Here, the measurement operator $\mathcal{A}_{i}(\cdot)$ is defined as $\mathcal{A}_{i}(\bm{X}_{\star})=\frac{1}{m}\langle\bm{a}_{i}\bm{a}_{i}^{\top},\bm{X}_{\star}\rangle$, where $\bm{a}_{i}$ is composed of i.i.d. Gaussian entries $\mathcal{N}(0,1)$. The ground truth matrix $\bm{X}_{\star}$ is positive semi-definite, and is generated via its compact SVD $\bm{X}_{\star}=\bm{U}_{\star}\bm{\Sigma}_{\star}\bm{U}_{\star}^{\top}$, where $\bm{U}_{\star}$ and $\bm{\Sigma}_{\star}$ are generated in the same manner described above.

Denote the index set of the remaining measurements after discarding $p_{s}$ fraction with largest amplitudes as $\mathcal{I}=\{i:|y_{i}|\leq|\bm{y}|_{(\lceil p_{s}m\rceil)}\}$, where $|\bm{y}|_{(k)}$ denotes the $k$th largest amplitude of $\bm{y}$. The truncated spectral method in [ZCL16, LCZL20] is used for initialization, where we apply the standard spectral method only on the subset $\mathcal{I}$ of the measurements. For matrix sensing, it follows the prescription in [LCZL20], and for quadratic sampling, it follows [LMCC21].

Fig. 1 shows the relative reconstruction error $\|\bm{X}_{t}-\bm{X}_{\star}\|_{\operatorname{\mathsf{F}}}/\|\bm{X}_{\star}\|_{\operatorname{\mathsf{F}}}$ for matrix sensing without outliers (in (a)) and with $20\%$ outliers (i.e. $p_{s}=0.2$ in (b)) under different condition numbers $\kappa$, where $\bm{X}_{t}$ is the estimated low-rank matrix at the $t$-th iteration. Fig. 2 shows the relative reconstruction error for quadratic sampling under the same setting. It can be seen that ScaledSM is insensitive to $\kappa$ and converges as a fast rate that is independent with $\kappa$, while the convergence of VanillaSM slows down dramatically with the increase of $\kappa$. In addition, both algorithms still converge linearly in the presence of outliers, thanks to the robustness of the least absolute deviations.

Comparisons of stepsize schedules.

We now compare the geometrically decaying stepsize with the Polyak’s stepsize for ScaledSM, which essentially mirrors similar experiments conducted in [LZSV20] for VanillaSM. We run ScaledSM for at most $T=1000$ iterations, and stop early if the relative error achieves $10^{-12}$. Fig. 5 and Fig. 6 show the performance comparisons of ScaledSM under various stepsize schedules for matrix sensing and quadratic sampling, respectively. For both figures, (a) shows the final relative error of ScaledSM using geometrically decaying stepsizes under various $(\lambda,q)$, where we see that ScaledSM converges as long as $\lambda$ is not too large and $q$ is not too small. We further plot the relative error versus the iteration count for ScaledSM using geometrically decaying stepsizes with a fixed $q$ and various $\lambda$ in (b), and with a fixed $\lambda$ and various $q$ in (c), where the performance using Polyak’s stepsizes is plotted for comparison. It can be seen that using Polyak’s stepsizes yields the fastest convergence. Indeed, if properly tuned, geometrically decaying stepsizes match Polyak’s stepsizes, as shown in (d). In general, we find that there is a wide range of parameters for geometrically decaying stepsizes where ScaledSM converges in a fast speed comparable to that of using Polyak’s stepsizes, as long as $\lambda$ is not too large and $q$ is not too small.

Discussions

This paper proposes scaled subgradient methods to minimize a family of nonsmooth and nonconvex formulations for low-rank matrix recovery—in particular, the residual sum of absolute errors—and guarantees its convergence at a rate that is almost dimension-free and independent of the condition number, even in the presence of corruptions. We illustrate the effectiveness of our approach by providing state-of-the-art performance guarantees for robust low-rank matrix sensing and quadratic sampling. In the future, it is of interest to study the performance of scaled subgradient methods for other signal estimation and statistical inference tasks, such as training student-teacher neural networks [DDKL20], as well as using random initializations [CCFM19].

Acknowledgements

The work of T. Tong and Y. Chi is supported in part by ONR under the grants N00014-18-1-2142 and N00014-19-1-2404, by ARO under the grant W911NF-18-1-0303, and by NSF under the grants CAREER ECCS-1818571, CCF-1806154 and CCF-1901199.

References

Appendix A Technical Lemmas

whenever the minimum is achieved.If there exist multiple minimizers, we arbitrarily choose one as $\bm{Q}$.

In particular, taking $\widetilde{\bm{X}}=\bm{X}+\mathcal{P}_{r}(\bm{S})$ arrives at

where the last equality follows from the definition (24). Note that $\mathcal{P}_{r}(\bm{S})$ has rank at most $r$. By the rank-$r$ restricted $L$-Lipschitz continuity of $f(\cdot)$, we have

Combining the above inequality with (27), we conclude $\|\bm{S}\|_{\operatorname{\mathsf{F}},r}\leq L$. ∎

Appendix B Proof of Theorem 1

Suppose that the $t$-th iterate $\bm{F}_{t}$ obeys the condition

Lemma 1 ensures that $\bm{Q}_{t}$, the optimal alignment matrix between $\bm{F}_{t}$ and $\bm{F}_{\star}$ exists. For notational convenience, we denote $\bm{L}\coloneqq\bm{L}_{t}\bm{Q}_{t}$, $\bm{R}\coloneqq\bm{R}_{t}\bm{Q}_{t}^{-\top}$, $\bm{\Delta}_{L}\coloneqq\bm{L}-\bm{L}_{\star}$, $\bm{\Delta}_{R}\coloneqq\bm{R}-\bm{R}_{\star}$, $\bm{S}\coloneqq\bm{S}_{t}$, and $\epsilon\coloneqq 0.02/\chi_{f}$. By the definition

and the relation $\|\bm{A}\bm{B}\|_{\operatorname{\mathsf{F}}}\geq\|\bm{A}\|_{\operatorname{\mathsf{F}}}\sigma_{r}(\bm{B})\geq\|\bm{A}\|\sigma_{r}(\bm{B})$, we have

where in the first line, we used the fact that the update rule (13) is covariant with respect to $\bm{Q}_{t}$, implying that

We proceed to bound $\mathfrak{S}_{1}$ and $\mathfrak{S}_{2}$. The term $\mathfrak{S}_{1}$ can be bounded by

where the second line follows from the condition (30) and Lemma 3 (cf. (23b)):

For the term $\mathfrak{S}_{2}$, note that

has rank at most $r$. Hence we can invoke Lemma 4 (cf. (26)) to obtain

where the second line follows from the triangle inequality, and the third line follows from $\|\bm{S}\|_{\operatorname{\mathsf{F}},r}\leq L$ (cf. Lemma 5), (30), and

Plugging collectively the bounds for $\mathfrak{S}_{1}$ and $\mathfrak{S}_{2}$ into (33) yields

Similarly, we can obtain the control of $\|(\bm{R}_{t+1}\bm{Q}_{t}^{-\top}-\bm{R}_{\star})\bm{\Sigma}_{\star}^{1/2}\|_{\operatorname{\mathsf{F}}}^{2}$. Combine them together to reach

Using the subgradient optimality of $\bm{S}$, we obtain

together with (29), which further implies that

Before proceeding to different cases of stepsize schedules, we record two useful properties. First, by the restricted $\mu$-sharpness of $f(\cdot)$ together with Lemma 2, we have

On the other end, by Lemma 4 (cf. (25)), we have

where the second line follows from $\|\bm{S}\|_{\operatorname{\mathsf{F}},r}\leq L$ (cf. Lemma 5) and

Let $\eta_{t}=\eta_{t}^{\mathsf{P}}$ be the Polyak’s stepsize in (14), which is

where the second line follows since $\bm{L}_{t}\bm{R}_{t}^{\top}=\bm{L}\bm{R}^{\top}$, $\bm{L}_{t}(\bm{L}_{t}^{\top}\bm{L}_{t})^{-1}\bm{L}_{t}^{\top}=\bm{L}(\bm{L}^{\top}\bm{L})^{-1}\bm{L}^{\top}$ and $\bm{R}_{t}(\bm{R}_{t}^{\top}\bm{R}_{t})^{-1}\bm{R}_{t}^{\top}=\bm{R}(\bm{R}^{\top}\bm{R})^{-1}\bm{R}^{\top}$. Plugging (37) into (34), we have

where the second line follows from (35) and $\chi_{f}=L/\mu$.

To continue, combining (35) and (36), we can lower bound the Polyak’s stepsize (37) as

where the contraction rate $\rho(\epsilon,\chi_{f})$ is

Under the condition $\epsilon=0.02/\chi_{f}$, we calculate $(1-\rho(\epsilon,\chi_{f}))\chi_{f}^{2}$ as

thus $\rho(\epsilon,\chi_{f})\leq 1-0.16/\chi_{f}^{2}$. We conclude that

B.2 Convergence with geometrically decaying stepsizes

Let $\eta_{t}=\eta_{t}^{\mathsf{G}}$ be the geometrically decaying stepsize in (15), which is

where the first line follows from (35) and $\chi_{f}=L/\mu$, and the second line follows from $\eta_{t}\geq\frac{\lambda q^{t}}{\sqrt{2}L}$ due to (36). We now aim to show that

in an inductive manner. Assume the above induction hypothesis holds at the $t$-iteration. By the setting of parameters, i.e.

where the contraction rate $\rho(\epsilon,\chi_{f})$ matches exactly (39). Therefore, under the condition $\epsilon=0.02/\chi_{f}$, we have $\rho(\epsilon,\chi_{f})\leq 1-0.16/\chi_{f}^{2}$, thus we conclude that

Appendix C Proof of Theorem 2

We start by introducing the short-hand notation $d_{t}\coloneqq(1-0.13/\chi_{f}^{2})^{t/2}0.02\sigma_{r}(\bm{X}_{\star})/\chi_{f}$. The parameters are set as

Follow the same derivations as the proof of Theorem 1 until (34). Plugging the stepsize (40) into (34), together with the approximate restricted sharpness property

Under the conditions $\chi_{f}\geq 1$ and $\epsilon=0.02/\chi_{f}\leq 0.02$, the above relation can be simplified to

We next prove the theorem by induction, where the base case is established trivially by the initial condition. By the induction hypothesis, the distance at the $t$-th iterate is bounded by

If $d_{t}\geq 20\xi/\mu$, or equivalently, $\xi\leq 0.05\mu d_{t}$, in view of (41), we have

where the third line uses the condition (40), and the last line holds since $d_{t}\geq 0$.

Appendix D Proof of Proposition 1

For $\bm{X}_{1}$ and $\bm{X}_{2}$ where $\bm{X}_{1}-\bm{X}_{2}$ has rank at most $2r$, we have

where the second line follows from the inverse triangle inequality and the assumed rank-$2r$ mixed-norm RIP (cf. Definition 4) of $\mathcal{A}(\cdot)$. As a result, we have $L=\delta_{2}$. On the other end, we note

where the first equality uses $f(\bm{X}_{\star})=0$ and the second inequality follows from the rank-$2r$ mixed-norm RIP; thus $\mu=\delta_{1}$.

Appendix E Proof of Proposition 2

where the second line follows from the inverse triangle inequality and the rank-$2r$ mixed-norm RIP; hence $L=\delta_{2}$. For approximate restricted sharpness, note that

where the second and the fourth lines follow from the triangle inequality, the third line follows from the definition of $\mathcal{S}$, and the last line follows from the definition of the $\mathcal{S}$-outlier bound and the noise upper bound $\|\bm{w}\|_{1}\leq\sigma_{w}$. Therefore, we have $\mu=\delta_{3}$ and $\xi=2\sigma_{w}$.