Low-distortion Subspace Embeddings in Input-sparsity Time and Applications to Robust Linear Regression

Introduction

Background

Given an $n\times d$ matrix $A$ and $p\in[1,\infty]$, let

For simplicity, we will use $\kappa_{p}$, $\sigma_{p}^{\min}$, and $\sigma_{p}^{\max}$ when the underlying matrix is clear.

Given an $n\times d$ matrix $A$ and $p\in[1,\infty]$, we always have

For a discussion of ellipsoidal rounding, we refer readers to Clarkson et al. . In this paper, we simply cite the following lemma, which is based on ellipsoidal rounding.

Stable distributions.

where $X_{i}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\operatorname{\mathcal{D}}$ and $X\sim\operatorname{\mathcal{D}}$. By “$X\simeq Y$”, we mean $X$ and $Y$ have the same distribution.

By a result due to Lévy , it is known that $p$-stable distributions exist for $p\in(0,2]$; and from Chambers et al. , it is known that $p$-stable random variables can be generated efficiently, thus allowing their practical use. Let us use $\operatorname{\mathcal{D}}_{p}$ to denote the “standard” $p$-stable distribution, for $p\in$, specified by its characteristic function $\psi(t)=e^{-|t|^{p}}$. It is known that $\operatorname{\mathcal{D}}_{1}$ is the standard Cauchy distribution, and that $\operatorname{\mathcal{D}}_{2}$ is the Gaussian distribution with mean and variance $2$.

Tail inequalities.

We note two inequalities from Clarkson et al. regarding the tails of the Cauchy distribution.

For $i=1,\ldots,m$, let $C_{i}$ be $m$ (not necessarily independent) standard Cauchy variables, and $\gamma_{i}>0$ with $\gamma=\sum_{i}\gamma_{i}$. Let $X=\sum_{i}\gamma_{i}|C_{i}|$. For any $t>1$,

For simplicity, we assume that $m\geq 3$ and $t\geq 1$, and then we have $\Pr[X>t\gamma]\leq 2\log(mt)/t$.

For $i=1,\ldots,m$, let $C_{i}$ be independent standard Cauchy random variables, and $\gamma_{i}\geq 0$ with $\gamma=\sum_{i}\gamma_{i}$. Let $X=\sum_{i}\gamma_{i}|C_{i}|$. Then, for any $t>0$,

We also note the following result about Gaussian variables. This is a direct consequence of Maurer’s inequality (), and we will use it to derive lower tail inequalities for $p$-stable distributions.

For $i=1,\ldots,m$, let $G_{i}$ be independent standard Gaussian random variables, and $\gamma_{i}\geq 0$ with $\gamma=\sum_{i}\gamma_{i}$. Let $X=\sum_{i}\gamma_{i}|G_{i}|^{2}$. Then, for any $t>0$,

In addition, $\Pi A$ can be computed in $\operatorname{\mathcal{O}}(\operatorname{nnz}(A))$ time.

The construction of $\Pi$ in this theorem is the same as the construction in Clarkson and Woodruff . For them, $s=\operatorname{\mathcal{O}}((d/\epsilon)^{4}\log^{2}(d/\epsilon))$ in order to achieve $(1\pm\epsilon)$ distortion with a constant probability. Theorem 1 shows that it actually suffices to set $s=\operatorname{\mathcal{O}}((d^{2}+d)/\epsilon^{2})$. Surprisingly, the proof is rather simple. Let $X=U^{T}\Pi^{T}\Pi U$, where $U$ is an orthonormal basis for $\operatorname{\mathcal{A}}_{2}$. Compute $\mathbf{E}[\|X-I\|_{F}^{2}]$ and apply Markov’s inequality to $\|X-I\|_{F}^{2}\leq\epsilon^{2}$, which implies $\|X-I\|_{2}\leq\epsilon$ and hence the embedding result. See Appendix A.1 for a complete proof.

Remark. The $\operatorname{\mathcal{O}}(\operatorname{nnz}(A))$ running time is indeed optimal, up to constant factors, for general inputs. Consider the case when $A$ has an important row $a_{j}$ such that $A$ becomes rank-deficient without it. Thus, we have to observe $a_{j}$ in order to compute a low-distortion embedding. However, without any prior knowledge, we have to scan at least a constant portion of the input to guarantee that $a_{j}$ is observed with a constant probability, which takes $\operatorname{\mathcal{O}}(\operatorname{nnz}(A))$ time. Note that this optimality result applies to general $p$.

In addition, $\Pi A$ can be computed in $\operatorname{\mathcal{O}}(\operatorname{nnz}(A))$ time.

Instead of analyzing $\operatorname{\mathcal{D}}_{p}$ directly, for any $p\in(1,2)$, we establish an order among the Cauchy distribution, the $p$-stable distribution, and the Gaussian distribution, and then we derive upper and lower tail inequalities for the $p$-stable distribution similar to the ones we used to prove Theorem 2. We state these technical results here since they are of independent interest. We start with the following lemma, which is proved in Appendix A.4 and which establishes this order.

For any $p\in(1,2)$, there exist constants $\alpha_{p}>0$ and $\beta_{p}>0$ such that

Our numerical results suggest that the constants $\alpha_{p}$ and $\beta_{p}$ are not too far away from $1$. See Figure 1, which plots of the CDFs of $|X_{p}/2|^{p}$ for $p=1,0,1.1,\ldots,2.0$, based on which we conjecture $|X_{p_{1}}/2|^{p_{1}}\succeq|X_{p_{2}}/2|^{p_{2}}$, for all $1\leq p_{1}\leq p_{2}\leq 2$. This implies that $2^{p-1}|C|\succeq|X_{p}|^{p}$ and $|X_{p}|^{p}\succeq 2^{p-2}|X_{2}|^{2}\simeq 2^{p-1}|G|^{2}$, which therefore provides a value for the constants $\alpha_{p}$ and $\beta_{p}$.

Lemma 8 suggests that we can use Lemma 5 (regarding Cauchy random variables) to derive upper tail inequalities for general $p$-stable distributions and that we can use Lemma 7 (regarding Gaussian variables) to derive lower tail inequalities for general $p$-stable distributions. The following two lemmas establish these results; the proofs of these lemmas are provided in Appendix A.5 and Appendix A.6, respectively.

Given $p\in(1,2)$, for $i=1,\ldots,m$, let $X_{i}$ be $m$ (not necessarily independent) random variables sampled from $\operatorname{\mathcal{D}}_{p}$, and $\gamma_{i}>0$ with $\gamma=\sum_{i}\gamma_{i}$. Let $X=\sum_{i}\gamma_{i}|X_{i}|^{p}$. Assume that $m\geq 3$. Then for any $t\geq 1$,

For $i=1,\ldots,m$, let $X_{i}$ be independent random variables sampled from $\operatorname{\mathcal{D}}_{p}$, and $\gamma_{i}\geq 0$ with $\gamma=\sum_{i}\gamma_{i}$. Let $X=\sum_{i}\gamma_{i}|c_{i}|$. Then,

In addition, $\Pi A$ can be computed in $\operatorname{\mathcal{O}}(\operatorname{nnz}(A))$ time.

In addition, $\Pi A$ can be computed in $\operatorname{\mathcal{O}}(\operatorname{nnz}(A)\cdot d\log d)$ time.

Improving the Embedding Dimension

time. The second term comes from solving a subsampled problem of size $\operatorname{\mathcal{O}}(d^{3+p/2}\log(1/\epsilon)/\epsilon^{2})\times d$.

Remark. We have stated our results in the previous sections as $\operatorname{poly}(d)$ without stating the value of the polynomial because there are numerous trade-offs between the conditioning quality and the running time. For example, let $p=1$. We can use a rounding algorithm instead of QR to compute the $R$ matrix. If we use the input-sparsity time embedding with the $\operatorname{\mathcal{O}}(d)$-rounding algorithm of , then the running time to compute the $(1\pm\epsilon)$-distortion embedding is $\operatorname{\mathcal{O}}(\operatorname{nnz}(A)\cdot\log n+d^{8}/\epsilon^{2})$ and the embedding dimension is $\operatorname{\mathcal{O}}(d^{6.5}/\epsilon^{2})$ (ignoring $\log$ factors). If, on the other hand, we use QR to compute $R$, then the running time is $\operatorname{\mathcal{O}}(\operatorname{nnz}(A)\cdot\log n+d^{7}/\epsilon^{2})$ and the embedding dimension is $\operatorname{\mathcal{O}}(d^{8}/\epsilon^{2})$. However, with the result from this section, the running time is simply $\operatorname{\mathcal{O}}(\operatorname{nnz}(A)\cdot\log n+\operatorname{poly}(d)+\operatorname{\mathcal{T}}_{p}(\epsilon;d^{3+p/2}/\epsilon^{2},d))$ and the $\operatorname{poly}(d)$ term can be absorbed by the $\operatorname{nnz}(A)$ term.

Acknowledgments

References

Appendix A Appendix

Let the $n\times d$ matrix $U$ be an orthonormal basis for the range of the $n\times d$ matrix $A$. Rather than proving the theorem by establishing that

where $s_{ij}$ is the $(i,j)$-th element of $S$, $d_{j}$ is the $j$-th diagonal element of $D$, and $u_{jk}$ is the $(j,k)$-th element of $U$. We will use the following facts in the proof:

Given these results, it is easy to obtain that

For any $\delta\in(0,1)$, set $s=(d^{2}+d)/(\epsilon^{2}\delta)$. Then, by Markov’s inequality,

Therefore, with probability at least $1-\delta$, we have $\|X-I\|_{2}\leq\|X-I\|_{F}\leq\epsilon$, which implies

We start with the following result, which establishes the existence of the so-called Auerbach’s basis of a $d$-dimensional normed vector space. For our proof, we will only need its existence and not an algorithm to construct it.

(Auerbach ) Let $(\operatorname{\mathcal{A}},\|\cdot\|)$ be a $d$-dimensional normed vector space. There exists a basis $\{e_{1},\ldots,e_{d}\}$ of $\operatorname{\mathcal{A}}$, called Auerbach basis, such that $\|e_{k}\|=1$ and $\|e^{k}\|^{*}=1$ for $k=1,\ldots,d$, where $\{e^{1},\ldots,e^{n}\}$ is a basis of $\operatorname{\mathcal{A}}^{*}$ dual to $\{e_{1},\ldots,e_{n}\}$.

For any $y=Ux\in Y$, we have $\|y\|_{1}=\|Ux\|_{1}\geq\|x\|_{\infty}=1$,

and thus $\|y\|_{1}\leq\|v\|_{1}=d$. Define $Y^{L}=\{y\in Y\,|\,\|y^{L}\|_{1}\geq\frac{1}{2}\|y\|_{1}\}$ and $Y^{H}=Y\backslash Y^{L}$. Given $S$, define a mapping $\phi:\{1,\ldots,n\}\to\{1,\ldots,s\}$ such that $s_{\phi(j),j}=1$, $j=1,\ldots,n$, and split $L$ into two subsets: $\hat{L}=\{j\in L\,|\,\phi(j)\in\phi(H)\}$ and $\bar{L}=L\backslash\hat{L}$. Consider these events:

$\operatorname{\mathcal{E}}_{U}$: $|\Pi U|_{1}\leq\omega_{1}d\log d$ for some $\omega_{1}>0$.

$\operatorname{\mathcal{E}}_{L}$: $\|Sv^{L}\|_{\infty}\leq\omega_{2}/(d\log d)$ for some $\omega_{2}>0$.

$\operatorname{\mathcal{E}}_{H}$: $\phi(j_{1})\neq\phi(j_{2}),\ \forall\,j_{1}\neq j_{2},\ j_{1},j_{2}\in H$.

$\operatorname{\mathcal{E}}_{C}$: $\min_{j\in|H|}|c_{j}|\geq\omega_{3}/(d^{2}\log^{2}d)$ for some $\omega_{3}>0$.

$\operatorname{\mathcal{E}}_{\hat{L}}$: $|\Pi U^{\hat{L}}|_{1}\leq\omega_{4}/(d^{2}\log^{2}d)$ for some $\omega_{4}>0$.

Recall that we set $s=\omega d^{5}\log^{5}d$ in Theorem 2. We will show that, with $\omega$ sufficiently large and proper choices of $\omega_{1}$, $\omega_{2}$, $\omega_{3}$, and $\omega_{4}$, the event $\operatorname{\mathcal{E}}_{U}$ leads to an upper bound of $\|\Pi y\|_{1}$ for all $y\in\text{range}(A)$, $\operatorname{\mathcal{E}}_{U}$ and $\operatorname{\mathcal{E}}_{L}$ lead to a lower bound of $\|\Pi y\|_{1}$ for all $y\in Y^{L}$ with probability at least $0.9$, and $\operatorname{\mathcal{E}}_{H}$, $\operatorname{\mathcal{E}}_{\hat{L}}$, and $\operatorname{\mathcal{E}}_{C}$ together imply an lower bound of $\|\Pi y\|_{1}$ for all $y\in Y^{H}$.

Provided $\operatorname{\mathcal{E}}_{U}$, we have

For any $y\in\text{range}(A)$, we can find an $x$ such that $y=Ux$. Then,

Provided $\operatorname{\mathcal{E}}_{L}$, for any fixed $y\in Y^{L}$, we have

By assumption $\operatorname{\mathcal{E}}_{L}$ and $\|y^{L}\|_{1}\geq\frac{1}{2}\|y\|_{1}\geq\frac{1}{2}$, we obtain the result. ∎

Assume both $\operatorname{\mathcal{E}}_{U}$ and $\operatorname{\mathcal{E}}_{L}$. If $\omega_{1}$ and $\omega_{2}$ satisfy

for some $\delta\in(0,1)$ regardless of $d$, then, with probability at least $1-\delta$, we have

Set $\epsilon=1/(2+8\omega_{1}d\log d)$ and create an $\epsilon$-net $Y^{L}_{\epsilon}\subseteq Y^{L}$ such that for any $y\in Y^{L}$, we can find a $y_{\epsilon}\in Y^{L}_{\epsilon}$ such that $\|y-y_{\epsilon}\|_{1}\leq\epsilon$. Since $\|y\|_{1}\leq d$ for all $y\in Y^{L}$, there exist such an $\epsilon$-net with at most $(3d/\epsilon)^{d}$ elements (Bourgain et al. ). By Lemma 13, we can apply a union bound for all the elements in $Y^{L}_{\epsilon}$:

For any $y\in Y^{L}$, we have, noting that $y-y_{\epsilon}\in\text{range}(A)$,

So we establish a lower bound for all $y\in Y^{L}$. ∎

Provided $\operatorname{\mathcal{E}}_{H}$ and $\operatorname{\mathcal{E}}_{\hat{L}}$, if $\omega_{3}>4\omega_{4}$, we have

which creates a lower bound for all $y\in Y^{H}$. ∎

We continue to show that, with $\omega$ sufficiently large, by setting $\tau=\omega^{1/4}/(d\log^{2}d)$ and choosing $\omega_{1}$, $\omega_{2}$, $\omega_{3}$, and $\omega_{4}$ properly, we have each event with probability at least $1-0.08=0.92$ and thus

Moreover, the condition in Lemma 14 holds with $\delta=0.1$, and the condition in Lemma 15 holds. Therefore, $\Pi=SC$ has the desired property with probability at least $0.5$, which would conclude the proof of Theorem 2.

With probability at least $0.92$, $\operatorname{\mathcal{E}}_{U}$ holds with $\omega_{1}=500(1+\log\omega)$.

Setting $\omega_{1}=500(1+\log\omega)$ and $t=\omega_{1}\log d$, we have

We assume that $\log d\geq 1$ and $\log\omega\geq 1$. ∎

For any $\delta\in(0,0.1)$, if $s\geq d/\tau$, we have,

Let $X_{ij}=s_{ij}v^{L}_{j}$. We have $\mathbf{E}[X_{ij}]=v^{L}_{j}/s$, $\mathbf{E}[X_{ij}^{2}]=(v^{L}_{j})^{2}/s$, and $0\leq X_{ij}\leq v^{L}_{j}\leq\tau$. Fixed $i$, $X_{ij}$ are independent, $j=1,\ldots,n$. By Bernstein’s inequality,

where we use Holder’s inequality: $\|v^{L}\|_{2}^{2}\leq\|v^{L}\|_{1}\|v^{L}\|_{\infty}\leq d\tau$. To obtain a union bound for all $i$ with probability $1-\delta$, we need

Given $\delta<0.1$, it suffices to choose $s=d/\tau$ and $t=2\log(d/(\delta\tau))\tau$. Note that $\|v^{L}\|_{1}/s\leq\|v\|_{1}/s=\tau$. We have

Increasing $s$ will decrease the failure rate, so it holds for all $s\geq d/\tau$. ∎

With probability at least $0.92$, $\operatorname{\mathcal{E}}_{L}$ holds with $\omega_{2}=(15+\log\omega)/\omega^{1/4}$.

By Lemma 17, with probability at least $0.92$, $\operatorname{\mathcal{E}}_{L}$ holds with

With the above choices of $\omega_{1}$ and $\omega_{2}$, the condition in Lemma 13 holds with $\delta=0.1$ for sufficiently large $\omega$.

With $\omega_{1}=500(1+\log\omega)$, and $\omega_{2}=(15+\log\omega)/\omega^{1/4}$, the first term in

increases much slower than the second term as $\omega$ increases, while both are at the order of $d\log d$. Therefore, if $\omega$ is sufficiently large, the condition hold with $\delta=0.1$. ∎

If $\omega\geq 160$, event $\operatorname{\mathcal{E}}_{H}$ holds with probability at least $0.92$.

Given $j_{1},j_{2}\in H$ and $j_{1}\neq j_{2}$, let $X_{j_{1}j_{2}}=1$ if $\phi(j_{1})=\phi(j_{2})$ and $X_{j_{1}j_{2}}=0$ otherwise. It is easy to see that $\Pr[X_{j_{1}j_{2}}=1]=\frac{1}{s}$. Therefore,

With probability at least $0.92$, event $\operatorname{\mathcal{E}}_{C}$ holds with $\omega_{3}=1/(8\omega^{1/4})$.

$|H|$ is at most $d/\tau=\omega^{1/4}d^{2}\log^{2}d$. Then

Therefore, $\omega_{3}=1/(8\omega^{1/4})$ would suffice. ∎

With probability at least $0.92$, event $\operatorname{\mathcal{E}}_{\hat{L}}$ holds with $\omega_{4}=25000(1+\log\omega)/\omega^{3/4}$. Thus with $\omega$ sufficiently large and the above choice of $\omega_{3}$, the condition in Lemma 15 $\omega_{3}>4\omega_{4}$ holds.

Assume that $|U^{\hat{L}}|_{1}\leq\frac{25}{\omega^{3/4}d^{2}\log^{3}d}$. Similar to the proof of Lemma 16, we have

It suffices to choose $t=1000(1+\log\omega)\log d$ to make the RHS less than $0.04$. So with probability at least $0.92$, we have $\operatorname{\mathcal{E}}_{\hat{L}}$ holds with $\omega_{4}=25000(1+\log\omega)/\omega^{3/4}$. ∎

By Theorem 2 and Lemma 3, we know that Steps 2 and 4 of Algorithm 1 succeed with a constant probability. Conditioning on this event, we have

where the last inequality is due to $\epsilon<1/2$. By Theorem 2, Step 2 takes $\operatorname{\mathcal{O}}(\operatorname{nnz}(A))$ time, and Step 3 takes $\operatorname{\mathcal{O}}(\operatorname{poly}(d))$ time because $\Pi A$ has $\operatorname{\mathcal{O}}(\operatorname{poly}(d)$ rows. Then, by Lemma 3, Step 4 takes $\operatorname{\mathcal{O}}(\operatorname{nnz}(A)\cdot\log n)$ time, and Step 5 takes $\operatorname{\mathcal{T}}_{1}(\epsilon/4;\operatorname{\mathcal{O}}(\operatorname{poly}(d)\log(1/\epsilon)/\epsilon^{2}),d)$ time. Therefore, the total running time of Algorithm 1 is as stated.

A.4 Proof of Lemma 8

Next, we state the following lemma, which is due to Nolan .

(Nolan [27, Thm. 1.12]) Let $X\sim\operatorname{\mathcal{D}}_{p}$ with $p\in[1,2)$. Then as $x\to\infty$,

where $c_{p}=\sin\frac{\pi p}{2}\cdot\Gamma(p)/\pi$.

By Lemma 23, it follows that, as $t\to\infty$,

Hence, there exist $\alpha_{p}^{\prime}>0$ and $t_{1}>0$ such that for all $t>t_{1}$,

Note that all the $p$-stable distributions with $p\in$ have finite and positive density at $x=0$. Therefore, there exists $\alpha_{p}^{\prime\prime}>0$ such that for all $0\leq t\leq t_{1}$,

Let $\alpha_{p}=\max\{\alpha_{p}^{\prime},\alpha_{p}^{\prime\prime}\}$. We get $\alpha_{p}|C|\succeq|X_{p}|^{p}$. For the Gaussian distribution, we have, as $t\to\infty$,

which converges to zero much faster than $t^{-1}$, so we can apply similar arguments to obtain $\beta_{p}$.

A.5 Proof of Lemma 9 (Upper Tail Inequality for p𝑝p-stable Distributions)

Let $C_{i}=F_{c}^{-1}(F_{p}(X_{i}))$, $i=1,\ldots,m$, where $F_{c}$ is the CDF of the standard Cauchy distribution and $F_{p}$ is the CDF of $\operatorname{\mathcal{D}}_{p}$. $C_{i}$ follows the standard Cauchy distribution, and, by Lemma 8, we have $\alpha_{p}|C_{i}|\geq|X_{i}|^{p}$. Therefore, for any $t\geq 1$,

A.6 Proof of Lemma 10 (Lower Tail Inequality for p𝑝p-stable Distributions)

Let $G_{i}$ be independent random variables sampled from the standard Gaussian distribution, $i=1,\ldots,m$. By Lemma 8, we have

The lower tail inequality from Lemma 7 concludes the proof.

A.8 Proof of Theorem 7 (Improving the Embedding Dimension)

Each of Steps 1, 3, and 5 of Algorithm 2 succeeds with a constant probability. We can control the success rate of each by adjusting the constant factor in the embedding dimension, such that all steps succeed with a constant probability. Conditioning on this event, we have $\kappa_{p}(AR^{-1})=6d$ because

By Lemma 1, $\bar{\kappa}_{p}(AR^{-1})\leq 6d^{1/p+1}$, and then by Lemma 3, the embedding dimension of $S$ is $\operatorname{\mathcal{O}}(\bar{\kappa}_{p}^{p}(AR^{-1})d^{|p/2-1|}d\log(1/\epsilon)/\epsilon^{2})=\operatorname{\mathcal{O}}(d^{3+p/2}\log(1/\epsilon)/\epsilon^{2})$.