Problem definition. The goal is to recover an $n\times q$ rank-$r$ matrix ${\bm{X}}^{*}=[\bm{x}^{*}_{1},\bm{x}^{*}_{2},\dots,\bm{x}^{*}_{q}]$ from $m$ linear projections (sketches) of each of its $q$ columns, i.e. from

where each $\bm{y}_{k}$ is an $m$-length vector, $[q]:=\{1,2,\dots,q\}$, and the measurement/sketching matrices $\bm{A}_{k}$ are mutually independent and known. The setting of interest is low-rank (LR), $r\ll\min(n,q)$, and undersampled measurements, $m<n$. Our guarantees assume that each $\bm{A}_{k}$ is random-Gaussian: each entry of it is independent and identically distributed (i.i.d.) standard Gaussian.

We also study the magnitude-only measurements’ setting, LRPR [4, 5, 6]. This involves recovering ${\bm{X}}^{*}$ from

Here $|\bm{z}|$ takes the entry-wise absolute value of entries of the vector $\bm{z}$.

Notation. Everywhere, $\|.\|_{F}$ denotes the Frobenius norm, $\|.\|$ without a subscript denotes the (induced) $l_{2}$ norm (often called the operator norm or spectral norm), $\|\bm{M}\|_{\max}$ is the maximum magnitude entry of the matrix $\bm{M}$, ^⊤ denotes matrix or vector transpose, and $|\bm{z}|$ for a vector $\bm{z}$ denotes element-wise absolute values. $\bm{I}_{n}$ (or sometimes just $\bm{I}$) denotes the $n\times n$ identity matrix. We use $\bm{e}_{k}$ to denote the $k$-th canonical basis vector, i.e., the $k$-th column of $\bm{I}$. For any matrix ${\bm{{Z}}}$, $\bm{z}_{k}$ denotes its $k$-th column.

We say ${\bm{U}}$ is a basis matrix if it contains orthonormal columns. For basis matrices ${\bm{U}}_{1},{\bm{U}}_{2}$, we use

as the Subspace Distance (SD) measure. For two $r$-dimensional subspaces, this is the $l_{2}$ norm of the sines of the $r$ principal angles between $\mathrm{span}({\bm{U}}_{1})$ and $\mathrm{span}({\bm{U}}_{2})$. $\mathrm{SD}({\bm{U}}_{1},{\bm{U}}_{2})$ is symmetric when ${\bm{U}}_{1},{\bm{U}}_{2}$ are both $n\times r$ basis matrices. Notice here we are using the Frobenius SD, unlike many recent works including our older work [5] that use the induced 2-norm one. This is done because it enables us to prove the desired guarantees easily. We reuse the letters $c,C$ to denote different numerical constants in each use with the convention that $c<1$ and $C\geq 1$. The notation $a\in\Omega(b)$ means $a\geq Cb$ while $a\in O(b)$ means $a\leq Cb$. We use $\mathbbm{1}_{\text{statement}}$ to denote an indicator function that takes the value 1 if statement is true and zero otherwise.

For a vector $\bm{w}$, we sometimes use $\bm{w}(k)$ to denote the $k$-th entry of $\bm{w}$. For a vector $\bm{w}$ and a scalar $\alpha$, $\mathbbm{1}(\bm{w}\leq\alpha)$ returns a vector of 1s and 0s of the same length as $\bm{w}$, with 1s where $(\bm{w}(k)\leq\alpha)$ and zero everywhere else. We use $\circ$ to denote the Hadamard product. Thus $\bm{z}:=\bm{w}\circ\mathbbm{1}(\bm{w}\leq\alpha)$ zeroes out entries of $\bm{w}$ larger than $\alpha$, while keeping the smaller ones as is.

For ${\bm{X}}^{*}$ which is a rank-$r$ matrix, we let

denote its reduced (rank $r$) SVD, i.e., ${\bm{U}}^{*}{}$ and ${\bm{V}^{*}}^{\top}$ are matrices with orthonormal columns (basis matrices), ${\bm{U}}^{*}{}$ is $n\times r$ and ${\bm{V}^{*}}$ is $r\times q$, and ${\bm{\Sigma}^{*}}$ is an $r\times r$ diagonal matrix with non-negative entries. We use $\kappa:={\sigma_{\max}^{*}}/{\sigma_{\min}^{*}}$ to denote the condition number of ${\bm{\Sigma}^{*}}$. This is not the condition number of ${\bm{X}}^{*}$ (whose minimum singular value is zero). We let $\bm{B}^{*}:={\bm{\Sigma}^{*}}\bm{V}^{*}{}{}$ and we use $\bm{b}^{*}_{k}$ to denote its $k$-th column.

We use the phrase $\epsilon$-accurate recovery to refer to $\mathrm{SD}({\bm{U}},{\bm{U}}^{*}{})\leq\epsilon$ or $\|{\bm{X}}-{\bm{X}}^{*}\|_{F}\leq\epsilon\|{\bm{X}}^{*}\|_{F}$ or both.

Assumption. Another way to understand (1) is as follows: each scalar measurement $\bm{y}_{ki}$ ($i$-th entry of $\bm{y}_{k}$) satisfies

with $\bm{a}_{ki}{}^{\top}$ being the $i$-th row of $\bm{A}_{k}$. Observe that the measurements are not global, i.e., no $\bm{y}_{ki}$ is a function of the entire matrix ${\bm{X}}^{*}$. They are global for each column ($\bm{y}_{ki}$ is a function of column $\bm{x}^{*}_{k}$) but not across the different columns. We thus need an assumption that enables correct interpolation across the different columns. The following assumption, which is a slightly weaker version of incoherence (w.r.t. the canonical basis) of right singular vectors suffices for this purpose.

Right singular vectors incoherence is the assumption $\max_{k}\|\bm{v}^{*}_{k}\|\leq\mu\sqrt{r/q}$. Since $\bm{b}^{*}_{k}={\bm{\Sigma}^{*}}\bm{v}^{*}_{k}$, this implies that the above holds. Incoherence of both left and right singular vectors was introduced for guaranteeing correct “interpolation” for the LRMC problem [3, 2].

Existing solutions for LRcCS and LRPR. Since it is always possible to obtain magnitude-only measurements ${\bm{y}_{(mag)}}_{k}$ from linear ones $\bm{y}_{k}$ as ${\bm{y}_{(mag)}}_{k}=|\bm{y}_{k}|$, a solution to LRPR also automatically solves LRcCS under the same assumptions. Hence the AltMin algorithm for LRPR from [4, 5] is the first provably correct solution for LRcCS. Of course, since LRcCS is an easier problem than LRPR, we expect a direct solution to LRcCS to need weaker assumptions. As we show in this paper, this is indeed true. A more recent work [7] studied the noisy version of LRcCS and developed a convex relaxation (mixed norm minimization) to provably solve it. Its time complexity is not discussed in the paper, however, it is well known that solvers for convex programs are much slower when compared to direct iterative algorithms: they either require number of iterations proportional to $1/\sqrt{\epsilon}$ or the per-iteration cost has cubic dependence on the problem size (here $((n+q)r)^{3}$) [2]. Thus, if $q\leq n$, its time complexity $O(mqnr\cdot\min(1/\sqrt{\epsilon},n^{3}r^{3}))$. In [6], we provided the best possible guarantee for the AltMin algorithm for solving LRPR, and hence LRcCS. We discuss these results in detail in Sec. II-D and summarize them in Table I.

Other well-studied LR recovery problems. The multivariate regression (MVR) problem, studied in [1], is our problem with $\bm{A}_{k}=\bm{A}$. However this is a very different setting than ours because, with $\bm{A}_{k}=\bm{A}$, the different $\bm{y}_{k}$’s are no longer mutually independent. As a result, one cannot exploit law of large numbers’ arguments over all $mq$ scalar measurements $\bm{y}_{ki}$. Consequently, the required value of $m$ can never be less than $n$. The result of [1] shows that $m$ of order $(n+q)r$ is both necessary and sufficient. LRMS involves recovering ${\bm{X}}^{*}$ from $\bm{y}_{i}=\langle\bm{A}_{i},{\bm{X}}^{*}\rangle,\ i=1,2,\dots,mq$ with $\bm{A}_{i}$ being dense matrices, typically i.i.d. Gaussian [2]. Thus all measurements are i.i.d. and global: each contains information about the entire quantity-of-interest, here ${\bm{X}}^{*}$. Because of this, for LRMS, one can prove a LR Restricted Isometry Property (RIP) that simplifies the rest of the analysis. This is what makes it very different from, and easier than, our problem.

LRMC, which involves recovering ${\bm{X}}^{*}$ from a subset of its observed entries, is the most closely related problem to ours since it also involves recovery from non-global measurements. The typical model assumed is that each matrix entry is observed with probability $p$ independent of others [3, 2]. Setting unobserved entries to zero, this can be written as $\bm{y}_{jk}=\delta_{jk}{\bm{X}}^{*}_{jk}$ with $\delta_{jk}\stackrel{{\scriptstyle\mathrm{iid}}}{{\thicksim}}Bernoulli(p)$. LRMC measurements are both row-wise and column-wise local. To allow correct interpolation across both rows and columns, it needs the incoherence assumption on both its left and right singular vectors. For our problem, the measurements are global for each column, but not across the different columns. For this reason, only right singular vectors’ incoherence is needed. In fact, because of the nature of our measurements, even if left incoherence were assumed, it would not help. This asymmetry in our measurement model and the fact that our measurements are unbounded (each $\bm{y}_{ki}$ is a Gaussian r.v) are two key differences between LRMC and LRcCS that prevent us from borrowing LRMC proof techniques for our work. Here symmetric means: if we replace ${\bm{X}}^{*}$ by its transpose, the probability distribution of the set of measurements does not change. Bounded means that the measurements’ magnitude has a uniform bound. This bound is $\|{\bm{X}}^{*}\|_{\max}$ for LRMC measurements.

Non-convex (iterative, not convex relaxation based) LRMC algorithms with the best sample complexity are GD-based. There are two common approaches for designing GD algorithms in the LR recovery literature, and in particular for LRMC. The first is to use standard projected GD on ${\bm{X}}$ (projGD-X), also referred to as Iterative Hard Thresholding: at each iteration, perform one step of GD for minimizing the squared loss cost function, $\tilde{f}({\bm{X}})$, w.r.t. ${\bm{X}}$, followed by projecting the resulting matrix onto the space of rank $r$ matrices (by SVD). This was studied in [14, 15] for solving LRMC. This is shown to converge geometrically with a constant GD step size, while needing only $\Omega((n+q)r^{2}\log^{2}n\log^{2}(1/\epsilon))$ samples on average.

The second is to let ${\bm{X}}={\bm{U}}\bm{B}$ where ${\bm{U}}$ is $n\times r$ and $\bm{B}$ is $r\times q$ and perform alternating GD for the cost function $f({\bm{U}},\bm{B}):=\tilde{f}({\bm{U}}\bm{B})$, i.e., update $\bm{B}$ with one step of GD for minimizing $f({\bm{U}},\bm{B})$ while keeping ${\bm{U}}$ fixed at its previous value, and then do the same for ${\bm{U}}$ with $\bm{B}$ fixed, and repeat. Since the ${\bm{X}}={\bm{U}}\bm{B}$ factorization is not unique, i.e., ${\bm{X}}={\bm{U}}{\bm{R}}^{-1}{\bm{R}}\bm{B}$ for any invertible $r\times r$ matrix ${\bm{R}}$, this approach can result in the norm of one of ${\bm{U}}$ or $\bm{B}$ growing in an unbounded fashion, while that of the other decreases at the same rate, causing numerical problems. A typical approach to resolve this issue, and one that was used for LRMC [16, 17], is to change the cost function to minimize to $f({\bm{U}},\bm{B})+\lambda f_{2}({\bm{U}},\bm{B})$ where $f_{2}({\bm{U}},\bm{B}):=\|{\bm{U}}^{\top}{\bm{U}}-\bm{B}\bm{B}^{\top}\|_{F}$ is the “norm-balancing term” (helps ensure that norms of ${\bm{U}}$ and $\bm{B}$ remain similar). We henceforth refer to this approach as altGDnormbal. The sample complexity bound for this approach is similar to that for projGD-X. But, it needs a GD step size of order $1/r$ or smaller [16, 17]; making it $r$-times slower than projGD-X.

AltGDmin algorithm. The main idea is as follows. Express ${\bm{X}}$ as ${\bm{X}}={\bm{U}}\bm{B}$ and alternatively update ${\bm{U}}$ and $\bm{B}$ as follows: (a) keeping $\bm{B}$ fixed at its previous value, update ${\bm{U}}$ by a GD step for it for the cost function $f({\bm{U}},\bm{B})$ followed by projecting the output onto the space of matrices with orthonormal columns; and (b) keeping ${\bm{U}}$ fixed at its previous value, update $\bm{B}$ by minimizing $f({\bm{U}},\bm{B})$ over it. Because of the column-wise decoupled form of our measurement model, step (b) is as fast as the GD step and thus the per-iteration time complexity of AltGD-Min is equal to that of any other GD method such as projGD-X or altGDnormbal. This decoupling (which means that, given ${\bm{U}}$, $\bm{b}_{k}$ only depends on $\bm{x}^{*}_{k}$, and not on the other columns of ${\bm{X}}^{*}$) also allows us to get the desired tight-enough bound on $\max_{k}\|\bm{b}_{k}-{\bm{U}}^{\top}\bm{x}^{*}_{k}\|$ and hence on $\max_{k}\|\bm{x}_{k}-\bm{x}^{*}_{k}\|$. This, and the fact that we use the gradient w.r.t. ${\bm{U}}$ in our algorithm, means that the summands in the gradient, and in other error bound terms, are nice-enough sub-exponential random variables (r.v.s): sub-exponential r.v.s whose maximum sub-exponential norm is small enough (is proportional to $(r/q)$), so that the summation can be bounded w.h.p. under the desired sample complexity.

AltGDmin analysis. When we analyzed the AltMin approach for LRPR [5, 6], we could directly modify proof techniques from AltMin for LRMC [2] for getting a bound on $\mathrm{SD}({\bm{U}},{\bm{U}}^{*}{})$ in terms of the bound on this distance from the previous iteration. We cannot do this for AltGD-Min because the algorithm itself is different from the two GD approaches studied for solving LRMC. We instead analyze AltGD-Min by a novel use of the fundamental theorem of calculus [18] that, along with other linear algebra tricks, helps us get a bound on $\mathrm{SD}({\bm{U}},{\bm{U}}^{*}{})$ which has the desired property: the terms in it are sums of nice-enough sub-exponentials. See Lemma 3.4 and its proof. The use of this result is motivated by its use in [19], and many earlier works, where it is used in a standard way: to bound the Euclidean distance, $\|\bm{x}-\bm{x}^{*}\|$, for standard GD to solve the PR problem for recovering a single vector $\bm{x}^{*}$. Thus, at the true solution $\bm{x}=\bm{x}^{*}$, the gradient of the cost function was zero. In our case, there are two differences: (i) we need to bound the subspace distance error, and (ii) our algorithm is not standard GD, and this means that $\nabla_{U}f({\bm{U}}^{*}{}{\bm{U}}^{*}{}^{\top}{\bm{U}},\bm{B})\neq 0$. We explain our approach in Sec. III-B.

AltGDmin initialization. The standard LR spectral initialization approach cannot be used because its summands are sub-exponential r.v.s that are not nice-enough. We give a detailed explanation in Appendix A. We address this issue by borrowing the truncation idea from the PR literature [20, 21, 5]. But, in our case, truncation is applied to a non-symmetric matrix. Thus the sandwiching arguments developed for symmetric matrices in [20], and modified in [21, 5], cannot be borrowed. We need a different argument which is used for proving Lemma B.2 and is briefly explained in Sec. III-D.

The LRcCS and LRPR problems occur in projection imaging applications involving sets of images, e.g., dynamic MRI [8, 9, 10], federated LR sketching [11, 7], and dynamic Fourier ptychography (LRPR) [22]. In MRI, Fourier projections of the region of interest, e.g., a cross-section of the brain or the heart, are acquired one coefficient at a time, making the scanning (data acquisition) quite slow. Hence, reduced sample complexity enables accelerated scanning. Since medical image sequences are usually slow changing, the LR model is a valid assumption for a time sequence [8, 9, 10]. In our notation, $\bm{x}^{*}_{k}$ is the vectorized version of the $k$-th image of the sequence and there are a total of $q$ images. The matrices $\bm{A}_{k}$ are random Fourier, i.e., $\bm{A}_{k}=\bm{H}_{k}{\bm{F}}$ where ${\bm{F}}$ is the $n\times n$ matrix that models computation of the 2D discrete Fourier transform as a matrix-vector operation, and $\bm{H}_{k}$ is an $m\times n$ random sampling “mask” matrix that models the frequency selection. In [13], we have shown the power of AltGD-Min for fast undersampled dynamic MRI of medical image sequences. It is both much faster, and in most cases, also provides better reconstructions, than many existing solutions from the MRI literature.

Large scale usage of smartphones results in large amounts of geographically distributed data, e.g., images. There is a need to compress/sketch this data before storing it. Sketch refers to a compression approach where the compression end is low complexity, usually simple linear projections [11, 7]. Consider the setting where different subsets of columns of ${\bm{X}}^{*}$ (each column corresponds to one vectorized image) are available at each of the $\rho\leq q$ nodes. The goal is to sketch them so that they can be correctly recovered using a federated algorithm. We can store the sketches $\bm{y}_{k}:=\bm{A}_{k}\bm{x}^{*}_{k}$ with $\bm{A}_{k}$’s being i.i.d. Gaussian. This way we store a total of only $mq$ scalars, with $mq$ of order roughly just $(n+q)r^{2}$. Traditional LR sketching approaches, e.g., [23], are designed for centralized settings and will not be efficient in a distributed setting.

In Sec. II, we develop AltGD-Min, give its guarantee for solving LRcCS, and compare it with existing results. We state and prove the two theorems that help prove our main result in Sec. III. This section also contains brief proof outlines before the actual proofs. The lemmas used in these proofs are proved in Sec. IV. The extension for solving LRPR is developed, and its guarantee is stated and proved, in Sec. V. We discuss the limitations of our results in Sec. VI. Simulation experiments are provided in Sec. VII. We conclude in Sec. VIII.

We would like to design a fast GD algorithm to find the matrix ${\bm{X}}$ that minimizes the squared-loss cost function $\tilde{f}({\bm{X}}):=\sum_{k=1}^{q}\|\bm{y}_{k}-\bm{A}_{k}\bm{x}_{k}\|^{2}.$ For reasons described earlier, we decompose ${\bm{X}}={\bm{U}}\bm{B}$ and develop an alternating GD-min (AltGD-Min) approach for the squared loss function,

Starting with a careful initialization for ${\bm{U}}$ explained below, AltGD-Min proceeds as follows. At each new iteration,

• •

  Min-B: update $\bm{B}$ by solving $\bm{B}\leftarrow\arg\min_{\tilde{\bm{B}}}f({\bm{U}},\tilde{\bm{B}})$. Since $\bm{b}_{k}$ only occurs in the $k$-th summand of $f({\bm{U}},\bm{B})$, this decouples to a much simpler column-wise least squares (LS) problem: $\bm{b}_{k}\leftarrow\arg\min_{\tilde{\bm{b}}_{k}}\|\bm{y}_{k}-\bm{A}_{k}{\bm{U}}\tilde{\bm{b}}_{k}\|^{2}$. This is solved in closed form as $\bm{b}_{k}=(\bm{A}_{k}{\bm{U}})^{\dagger}\bm{y}_{k}$ for each $k$; here $\bm{M}^{\dagger}:=(\bm{M}^{\top}\bm{M})^{-1}\bm{M}^{\top}$.

• •

  ProjGD-U: update ${\bm{U}}$ by one GD step for it, $\hat{\bm{U}}^{+}\leftarrow{\bm{U}}-\eta\nabla_{U}f({\bm{U}},\bm{B})$, followed by projecting $\hat{\bm{U}}^{+}$ onto the space of matrices with orthonormal columns to get the updated ${\bm{U}}^{+}$. We get ${\bm{U}}^{+}$ by QR decomposition: $\hat{\bm{U}}^{+}\overset{\mathrm{QR}}{=}{\bm{U}}^{+}{\bm{R}}^{+}$.

Notice that, because of the decoupling for $\bm{B}$, the min step only involves solving $q$ $r$-dimensional Least Squares (LS) problems, in addition to also first computing the matrices, $\bm{A}_{k}{\bm{U}}$. Computing the matrices needs time of order $mnr$, and solving one LS problem needs time of order $mr^{2}$. Thus, the LS step needs time $O(q\max(mnr,mr^{2}))=O(mqnr)$ since $r\leq n$. This is equal to the time needed to compute the gradient w.r.t. ${\bm{U}}$; and thus, the per-iteration cost of AltGD-Min is only $O(mqnr)$. The QR decomposition of an $n\times r$ matrix takes time only $nr^{2}$.

Since $f({\bm{U}},\bm{B})$ is not a convex function of the unknowns $\{{\bm{U}},\bm{B}\}$, a careful initialization is needed. Borrowing the spectral initialization idea from LRMC and LRMS solutions, we should initialize ${\bm{U}}_{0}$ by computing the top $r$ singular vectors of

Clearly the expected value of the $k$-th column of this matrix equals $\bm{x}^{*}_{k}$ and thus $\mathbb{E}[{\bm{X}}_{0,full}]={\bm{X}}^{*}$. But, as we explain next, it is not clear how to prove that this matrix concentrates around ${\bm{X}}^{*}$. Observe that it can also be written as

Its summands are independent sub-exponential r.v.s with maximum sub-exponential norm $\max_{k}\|\bm{x}^{*}_{k}\|\leq\mu\sqrt{r/q}{\sigma_{\max}^{*}}$. This is too large and does not allow us to bound $\|{\bm{X}}_{0,full}-{\bm{X}}^{*}\|$ under the desired sample complexity; see Appendix A. To resolve this issue, we borrow the truncation idea from earlier work on PR [20, 5] and initialize ${\bm{U}}_{0}$ as the top $r$ left singular vectors of

where $\alpha:=\tilde{C}\frac{\sum_{ki}(\bm{y}_{ki})^{2}}{mq}$ and $\bm{y}_{k,trunc}(\alpha):=\bm{y}_{k}\circ\mathbbm{1}(|\bm{y}_{k}|\leq\sqrt{\alpha})$. We set $\tilde{C}$ in our main result. Observe that we are summing over only those $i,k$ for which $\bm{y}_{ki}^{2}$ is not too large (is not much larger than its empirically computed average value). This truncation filters out the too large (outlier-like) measurements and sums over the rest. Theoretically, this converts the summands into sub-Gaussian r.v.s which have lighter tails than the un-truncated ones. This allows us to prove the desired concentration bound. Different from the above setting, in [20, 5], truncation was applied to symmetric positive definite matrices and was used to convert summands that were heavier-tailed than sub-exponential to sub-exponential.

We summarize the complete algorithm in Algorithm 1. This uses sample-splitting which is a commonly used approach in the LR recovery literature [2, 14, 15] as well as in other compressive sensing settings. It helps ensure that the measurement matrices in each iteration for updating ${\bm{U}}$ and $\bm{B}$ are independent of all previous iterates. This allows one to use concentration bounds for sums of independent r.v.s. We provide a detailed discussion in Sec. VI-A.

We can prove the following result.

Sample complexity The sample complexity (total number of samples needed to achieve $\epsilon$-accurate recovery) is $m_{\mathrm{tot}}=\sum_{\tau=0}^{T}m_{\tau}\geq m_{0}+T\min_{t\geq 1}m_{t}$. From the above result, this needs to satisfy $m_{tot}q\geq C\kappa^{6}\mu^{2}(n+q)r^{2}\log(1/\epsilon)\log(\kappa)$ and $m_{tot}>C\kappa^{2}\max(r,\log q,\log n)\log(1/\epsilon)$.

Time complexity Let $m\equiv m_{t}$. The initialization step needs time $mqn$ for computing ${\bm{X}}_{0}$; and time of order $nqr$ times the number of iterations used in the $r$-SVD step. Since we only need a $\delta_{0}$-accurate initial estimate of $\mathrm{span}({\bm{U}}^{*}{})$, with $\delta_{0}=c/\kappa^{2}$, order $\log(\kappa)$ number of iterations suffice for this SVD step. Thus the complexity is $O(nq(m+r)\cdot\log\kappa)=O(mqn\cdot\log\kappa)$ since $m\geq r$. One gradient computation needs time $O(mqnr)$. The QR decomposition needs time of order $nr^{2}$. The update of columns of $\bm{B}$ by LS also needs time $O(mqnr)$ (explained earlier). As we prove above, we need to repeat these steps $T=O(\kappa^{2}\log(1/\epsilon))$ times. Thus the total time complexity is $O(mqn\log\kappa+\max(mqnr,nr^{2},mqnr)\cdot T)=O(\kappa^{2}mqnr\log(1/\epsilon)\log\kappa)$.

Communication complexity The communication complexity per node per iteration for a federated implementation is just order $nr$. Thus, the total is $O(nr\cdot\kappa^{2}\log(1/\epsilon))$.

Thus, we have the following corollary.

Observe that the above results show that after $T=C\kappa^{2}\log(1/\epsilon)$ iterations, $\mathrm{SD}({\bm{U}}_{T},{\bm{U}}^{*}{})\leq\epsilon$, ${\|\bm{x}_{k}-\bm{x}^{*}_{k}\|}\leq\epsilon{\|\bm{x}^{*}_{k}\|}$, and $\|{\bm{X}}_{T}-{\bm{X}}^{*}\|_{F}\leq 1.4\epsilon\|{\bm{X}}^{*}\|$. The RHS in the third bound does indeed contain $\|{\bm{X}}^{*}\|$ (the induced 2-norm). This is correct because, $\mathrm{SD}(.,.)$ is a Frobenius norm subspace distance. We explain this in Sec. III-B.

An algorithm is called linear time if its time complexity is the same order as the time needed to load all input data. In our case, this is $O(mqn)$. Treating $\kappa$ as a constant, the AltGD-Min complexity is worse than linear-time by a factor of only $r\log(1/\epsilon)$. As can be seen from Table I, the same is also true for the fastest LRMC solution, projGD-X [15]. For LRMC, linear time is $O(mq)$. To our best knowledge, this is the case for the fastest algorithms for all LR problems.

Consider the sample complexity. The degrees of freedom (number of unknowns) of a rank-$r$ $n\times q$ matrix are $(n+q)r$. A sample complexity of $\Omega((n+q)r)$ samples (or, sometimes this times log factors) is called “optimal”. Thus, ignoring the log factors, our sample complexity of $m_{tot}q\gtrsim(n+q)r^{2}$ is sub-optimal only by a factor of $r$. As can also be seen from Table I, this suboptimality matches that of the best results for LRMC solutions that are not convex relaxation based [15, 16, 17]. The need for exploiting incoherence while obtaining the high probability bounds on the recovery error terms is what introduces the extra factor of $r$ for both LRMC and LRcCS. LRMC has been extensively studied for over a decade and there does not seem to be a way to obtain an (order-) optimal sample complexity guarantee for it except when studying convex relaxation solutions (which are much slower).

In addition, we also need $m\gtrsim\max(r,\log q,\log n)$. This is redundant except for very large $q,n$. This is needed because, the recovery of each column of $\bm{B}^{*}$ is a decoupled $r$-dimensional LS problem. We analyze this step in Lemma 3.3; notice that the bound on the recovery error of column $k$ holds w.p. at least $1-\exp(r-cm)$. By union bound, it holds for all $q$ columns w.p. at least at least $1-q\exp(r-cm)=1-\exp(\log q+r-cm)$. This probability is at least $1-n^{-10}=1-\exp(-10\log n)$ if $m\gtrsim\max(r,\log q,\log n)$.

There are two existing solutions for LRcCS – AltMin [4, 5, 6] and the convex relaxation (mixed norm minimization) [7]. Mixed norm is defined as $\|{\bm{X}}\|_{mixed}:=\inf_{\{{\bm{U}},\bm{V}:{\bm{U}}\bm{V}={\bm{X}}\}}\|{\bm{U}}\|_{F}\max_{k\in[q]}\|\bm{v}_{k}\|$, where ${\bm{U}}$ is $n\times r$ and $\bm{V}:=[\bm{v}_{1},\bm{v}_{2},\dots\bm{v}_{q}]$ is an $r\times q$ matrix. In our notation, for the noise-free case ($\sigma=0$), their main result states the following.