The Fast Cauchy Transform and Faster Robust Linear Regression

Introduction

Preliminaries

The following two Bernstein-type tail inequalities are useful because they give tail bounds without reference to the number of i.i.d. trials. The first bound is due to Maurer , and the second is an immediate application of the first.

Let $X_{i}\geq 0$ be independent random variables with $\sum_{i}\hbox{\bf{E}}[X_{i}^{2}]<\infty$, and define $X=\sum_{i}X_{i}$. Then, for any $t>0$,

Let $x_{i}$ be i.i.d. Bernoulli random variables with probability $p$, and let $X=\sum_{i\in[n]}\xi_{i}x_{i}$, where $\xi_{i}\geq 0$, with $\sum_{i\in[n]}\xi_{i}=\xi$ and $\sum_{i\in[n]}\xi_{i}^{2}\leq\xi^{2}/\beta^{2}$. Then, for any $t>0$,

The proof is a straightforward application of Lemma 1 to $Z=\sum_{i\in[n]}\xi_{i}(1-x_{i})$. ∎

The Cauchy distribution, having density $p(x)={1\over\pi}{1\over 1+x^{2}}$, is the unique $1$-stable distribution. If $C_{1},\ldots,C_{M}$ are independent Cauchys, then $\sum_{i\in[M]}\gamma_{i}C_{i}$ is distributed as a Cauchy scaled by $\gamma=\sum_{i\in[M]}|\gamma_{i}|$. The Cauchy distribution will factor heavily in our discussion, and bounds for sums of Cauchy random variables will be used throughout. We note that the Cauchy distribution has undefined expectation and infinite variance.

The following upper and lower tail inequalities for sums of Cauchy random variables are proved in Appendix A. The proof of Lemma 3 is similar to an argument of Indyk , though in that paper the Cauchy random variables are independent. As in that paper, our argument follows by a Markov bound after conditioning on the magnitudes of the Cauchy random variable summands not being too large, so that their conditional expectations are defined. However, in this paper, the Cauchy random variables are dependent, and so after conditioning on a global event, the expectations of the magnitudes need not be the same as after this conditioning in the independent case.

Lemma 4 is a simple application of Lemma 1, while Lemma 5 was shown in ; we include the proofs for completeness.

For $i\in[m]$, let $C_{i}$ be $m$ (not necessarily independent) Cauchy random variables, and $\gamma_{i}>0$ with $\gamma=\sum_{i\in[m]}\gamma_{i}$. Let $X=\sum_{i\in[m]}\gamma_{i}|C_{i}|$. Then, for any $t\geq 1$,

Remark. The bound has only logarithmic dependence on the number of Cauchy random variables and does not rely on any independence assumption among the random variables. Even if the Cauchys are independent, one cannot substantially improve on this bound due to the nature of the Cauchy distribution. This is because, for independent Cauchys, $\sum_{i}\gamma_{i}|C_{i}|\geq|\sum_{i}\gamma_{i}C_{i}|$, and the latter sum is itself distributed as a Cauchy scaled by $\gamma$. Hence for independent Cauchys, $\hbox{\bf{Pr}}[X\geq\gamma t]={2\over\pi}\tan^{-1}t=\Omega({1\over t})$.

For $i\in[r]$, let $C_{i}$ be independent Cauchy random variables, and $\gamma_{i}\geq 0$ with $\gamma=\sum_{i\in[r]}\gamma_{i}$ and $\sum_{i\in[r]}{\gamma_{i}^{2}}\leq\gamma^{2}/\beta^{2}$. Let $X=\sum_{i\in[r]}\gamma_{i}|C_{i}|$. Then, for any $t\geq 0$,

where $\displaystyle\delta\leq 2\operatorname{exp}\left({-as\varepsilon^{2}{\|Zx\|}_{1}\over(2+{2\over 3}\varepsilon){\|Z\|}_{1}{\|x\|}_{\infty}}\right).$

Main Technical Result: the Fast Cauchy Transform

This FCT construction first preprocesses by a deterministic low-coherence “spreading matrix,” then rescales by Cauchy random variables, and finally samples linear combinations of the rows. Let $\delta\in(0,1]$ be a parameter governing the failure probability of our algorithm. Then, we construct $\Pi_{1}$ as

$B\in\R^{r_{1}\times 2n}$ has each column chosen independently and uniformly from the $r_{1}$ standard basis vectors for $\R^{r_{1}}$; we will set the parameter $r_{1}=\alpha d\log{d\over\delta}$, where $\delta$ controls the probability that our algorithms fail and $\alpha$ is a suitably large constant;

$C\in\R^{2n\times 2n}$ is a diagonal matrix with diagonal entries chosen independently from a Cauchy distribution; and

(For completeness, we remind the reader that the (non-normalized) $n\times n$ matrix of the Hadamard transform $H_{n}$ may be defined recursively as follows:

There is a distribution (given by the above construction) over matrices $\Pi_{1}\in\R^{r_{1}\times n}$, with $r_{1}=O(d\log d+d\log{1\over\delta})$, such that for an arbitrary (but fixed) $A\in\R^{n\times d}$, and for all $x\in\R^{d}$, the inequalities

Further, for any $y\in\R^{n}$, the product $\Pi_{1}y$ can be computed in $O(n\log r_{1})$ time.

Setting $\delta$ to a small constant, since $\sqrt{s}=r_{1}^{3}$ and $r_{1}=O(d\log d)$, it follows that $\kappa=O(d^{4}\log^{4}d)$ in the above theorem.

The existence of such a $\Pi_{1}$ satisfying bounds of the form (2) was established by Sohler and Woodruff . Here, our contribution is to show that $\Pi_{1}$ can be factored into structured matrices so that the product $\Pi_{1}A$ can be computed in $O(nd\log d)$ time. We also remark that, in additional theoretical bounds provided by the FJLT, high-quality numerical implementations of variants of the Hadamard transform exist, which is an additional plus for our empirical evaluations of Theorem 1 and Theorem 2.

Our proof of this theorem uses a tail bound for ${\|Bg\|}_{2}^{2}$ in terms of ${\|g\|}_{2}$ and ${\|g\|}_{1}$, where $g$ is any positive vector in $\R^{n}$, and $B$ is the matrix used in our FCT construction. ${\|Bg\|}_{2}^{2}=\sum_{j}{\gamma_{j}^{2}}$ where $\gamma_{j}=\sum_{i}B_{ji}g_{i}$ are anti-correlated random variables. To get concentration, we independently bounded $\gamma_{j}^{2}$ in our proof which required $s=r_{1}^{6}$ to obtain the high probability result; this resulted in the bound $\kappa=O(d^{4}\log^{4}d)$.

2 FCT2 Construction: via a Fast Johnson-Lindenstrauss Transform

This FCT construction first preprocesses by a FJLT and then rescales by Cauchy random variables. Recall that $\delta\in(0,1]$ is a parameter governing the failure probability of our algorithm; and let $\eta>0$ be a generic arbitrarily small positive constant (whose value may change from one formula to another). Let $r_{1}=c\cdot d\log{d\over\delta}$, $s=c^{\prime}\cdot(d+\log{n\over\delta})$, and $t=s^{2+\eta}$, where the parameters $c,c^{\prime}>0$ are appropriately large constants. Then, we construct $\Pi_{1}\in\R^{r_{1}\times n}$ as

$C\in\R^{r_{1}\times ns/t}$ is a matrix of independent Cauchy random variables; and

There is a distribution (given by the above construction) over matrices $\Pi_{1}\in\R^{r_{1}\times n}$, with $r_{1}=O(d\log{d\over\delta})$, such that for arbitrary (but fixed) $A\in\R^{n\times d}$, and for all $x\in\R^{d}$, the inequalities

hold with probability $1-\delta$, where $\kappa=O({d\over\delta}(d+\log{n\over\delta})^{1+\eta}\log d)$. Further, for any $y\in\R^{n}$, the product $\Pi_{1}y$ can be computed in $O(n\log{d\over\delta})$ time.

Remark. For $\log n<d$, FCT2 gives a better dependence of the distortion on $d$, but more generally FCT2 has a dependence on $\log n$. This dependence arises because the random FJLT matrix does not give a deterministic guarantee for spreading out a vector whereas the low coherence matrix used in FCT1 does give a deterministic guarantee. This means that in using the union bound, we need to overcome a factor of $n$.

Remark. The requirement $t\geq s^{2+\eta}$ is set by the restriction in Lemma 9 in the proof of Theorem 2. In the bound of Theorem 2, $\kappa=\kappa^{\prime}\sqrt{t}$, where $\kappa^{\prime}=O(d\log(r_{1}d))$ arises from Theorem 12, which originally appeared in . If a stronger version of Lemma 9 can be proved that relaxes the restriction $t\geq s^{2+\eta}$, then correspondingly the bound of Theorem 2 will improve.

A basis $U$ for the range of $A$ is $(\alpha,\beta)$-conditioned if ${\|U\|}_{1}\leq\alpha$ and for all $x\in\R^{d}$, ${\|x\|}_{\infty}\leq\beta{\|Ux\|}_{1}$. We will say that $U$ is well-conditioned if $\alpha$ and $\beta$ are low-degree polynomials in $d$, independent of $n$.

Given an $n\times d$ matrix $A$, let $\Pi_{1}\in\R^{r_{1}\times n}$ be any projection matrix such that for any $x\in\R^{d}$,

For example, it could be constructed with either of the FCT constructions described in Section 3, or with the “slow” Cauchy Transform of , or via some other means. After computing the matrix $\Pi_{1}$, the FastL1Basis algorithm of Figure 1 consists of the following steps: construct $\Pi_{1}A$ and an $R$ such that $\Pi_{1}A=QR$, where $Q$ has orthonormal columns (for example using a QR-factorization of $\Pi_{1}A$); and then return $U=AR^{-1}=A(Q^{T}\Pi_{1}A)^{-1}$.

The next theorem and its corollary are our main results for the FastL1Basis algorithm; and this theorem follows by combining our Theorem 2 with Theorems 9 and 10 of . The proof of this theorem may be found in Appendix D.

For any $A\in\R^{n\times d}$, the basis $U=A{R^{-1}}$ constructed by FastL1Basis$(A)$ of Figure 1 using any $\Pi_{1}$ satisfying (3) is a $(d\sqrt{r_{1}},\kappa)$-conditioned basis for the range of $A$.

If $\Pi_{1}$ is obtained from the FCT2 construction of Theorem 2, then the resulting $U$ is an $(\alpha,\beta)$-conditioned basis for $A$, with $\alpha=O(d^{3/2}\log^{1/2}d)$ and $\beta=O(d^{2+\eta}\log d)$, with probability $1-\delta$. The time to compute the change of basis matrix ${R^{-1}}$ is $O(nd\log d+d^{3}\log d)$, assuming $\log n=O(d)$ and $\delta>0$ is a fixed constant.

Remark. Our constructions that result in $\Pi_{1}$ satisfying (3) do not require that $A\in\R^{n\times d}$; they only require that $A$ have rank $d$, and so can be applied to any $A\in\R^{n\times m}$ having rank $d$. In this case, a small modification is needed in the construction of $U$, because $R\in\R^{d\times m}$, and so we need to use $R^{\dagger}$ instead of $R^{-1}$. The running time will involve terms with $m$. This can be improved by processing $A$ quickly into a smaller matrix by sampling columns so that the range is preserved (as in ), which we do not discuss further.

and an $x^{*}$ achieving this minimum. We start with our main algorithm and theorem for this problem; and we then describe how a somewhat more sophisticated version of the algorithm yields improved running time bounds.

Prior work has shown that there is a diagonal sampling matrix $D$ with a small number of nonzero entries so that $\hat{x}=\operatorname{argmin}\nolimits_{x\in\R^{d}}\lVert D(Ax-b)\rVert_{1}$ satisfies

The next theorem summarizes our main quality-of-approximation results for the FastCauchyRegression algorithm of Figure 2. It improves the $O(nd^{2}+\operatorname{poly}(d\varepsilon^{-1}\log n))$ algorithm of , which in turn improved the result in . (Technically, the running time of is $O(nd^{\omega^{+}-1}+\operatorname{poly}(d\varepsilon^{-1}\log n))$, where $\omega^{+}$ is any constant larger than the exponent for matrix multiplication; for practical purposes, we can set $\omega^{+}=3$.) Our improved running time comes from using the FCT and a simple row-norm estimator for the row-norms of a well-conditioned basis. The proof of this theorem may be found in Appendix E.

Given are $\varepsilon\in(0,1)$, $\rho>0$, $A\in\R^{n\times d}$ and $b\in\R^{n}$. FastCauchyRegression$(A,b)$ constructs a coreset specified by the diagonal sampling matrix $D$ and a solution vector $\hat{x}\in\R^{d}$ that minimizes the weighted regression objective ${\|D(Ax-b)\|}_{1}$. The solution $\hat{x}$ satisfies, with probability at least $1-{1\over d^{\rho}}$ ($\rho>0$ is a constant),

Further, with probability $1-o(1)$, the entire algorithm to construct $\hat{x}$ runs in time

Given are $\varepsilon\in(0,1)$, $\rho>0$, $A\in\R^{n\times d}$ and $b\in\R^{n}$. OptimizedFastCauchyRegression$(A,b)$ constructs a coreset specified by the diagonal sampling matrix $D$ and a solution vector $\hat{x}\in\R^{d}$ that minimizes the weighted regression objective ${\|D(Ax-b)\|}_{1}$. The solution $\hat{x}$ satisfies, with probability at least $1-{1\over d^{\rho}}-{1\over\log^{\rho}n}$,

Further, with probability $1-o(1)$, the entire algorithm to construct $\hat{x}$, runs in time

Note that our algorithms and results also extend to multiple regression with $b\in\R^{n\times k}$, a fact that will be exploited in the next section.

A detailed inspection of the proof of Theorem 4 in Section 4.2 (see Appendix E for the proof) reveals that nowhere is it necessary that $x$ be a vector, i.e., the whole proof generalizes to a matrix $Z$. In particular, the inequalities in (9) continue to hold, since if they hold for every vector $x$, then it must hold for a matrix $Z$ because ${\|XZ\|}_{1}=\sum_{j\in[k]}{\|XZ^{(j)}\|}_{1}$. Similarly, if Lemma 13 continues to hold for vectors then it will imply the desired result for matrices, and so the only change in all the algorithms and results is that the short dimension of $X$ changes from $d+1$ to $d+k$. Thus, by shrinking $\delta$ by an additional factor of $k$, and taking a union bound we get a relative error approximation for each individual regression. We refer to this modified algorithm, where a matrix $B$ is input and the optimization problem in the last step is modified appropriately, as FastCauchyRegression$(A,B)$, overloading notation in the obvious way. This discussion is summarized in the following theorem.

Given $\varepsilon\in(0,1)$, $\rho>0$, a matrix $A\in\R^{n\times d}$ and $B\in\R^{n\times k}$, FastCauchyRegression$(A,B)$ constructs a coreset specified by the diagonal sampling matrix $D$ and a solution ${\hat{W}}\in\R^{d\times k}$ that minimizes the weighted multiple regression objective ${\|D(AW-B)\|}_{1}$. The solution $\hat{W}$ satisfies, with probability at least $1-{1\over(d+k)^{\rho}}$,

Further, with probability $1-o(1)$, the entire algorithm to construct $\hat{W}$, runs in time

First, we can save an extra factor of $(d+k)$ in $s$ in the above theorem if all we want is a relative error approximation to the entire multiple regression and we do not need relative error approximations to each individual regression.

where the constraint set is ${\cal C}=\{W\in\R^{d\times d}:W_{ii}=-1\}$. Since the constraint set effectively places an independent constraint on each column of $W$, after some elementary manipulation, it is easy to see that this regression is equivalent to the $d$ individual regressions to obtain $w_{j}^{*}$. Indeed, for an optimal solution $W^{*}$, we can set $w_{j}^{*}=W^{*(j)}$.

where (a) is from the $(1+\varepsilon)$-optimality of the constrained multiple regression as analyzed in Appendix E and (b) is because $j^{*}$ attained minimum error among all $j\in[d]$. This discussion is summarized in the following theorem.

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

The following lemma characterizes the relationship between these two quantities.

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

Thus, $\kappa_{p}(A)\leq d^{|1/p-1/2|}\alpha\beta$. ∎

Although it is easier to describe sampling algorithms in terms of $\bar{\kappa}_{p}$, after we show the equivalence between $\kappa_{p}$ and $\bar{\kappa}_{p}$, it will be easier for us to discuss conditioning algorithms in terms of $\kappa_{p}$, which naturally connects to ellipsoidal rounding algorithms.

Therefore, we have $\kappa_{p}(A{R^{-1}})\leq\kappa$. So a $\kappa$-rounding of $\mathcal{C}$ leads to a $\kappa$-conditioning of $A$.

2 Fast Ellipsoidal Rounding

Applying Theorem 8 to the convex set $\mathcal{C}=\{x\,|\,\|Ax\|_{p}\leq 1\}$, with the separation oracle described via a subgradient of $\|Ax\|_{p}$ and the initial rounding provided by the “$R$” matrix from the QR decomposition of $A$, we improve the running time of the algorithm used by Clarkson and by Dasgupta et al. from $\mathcal{O}(nd^{5}\log n)$ to $\mathcal{O}(nd^{3}\log n)$ while maintaining an $\mathcal{O}(d)$-conditioning. The proof of this theorem may be found in Appendix G.2.

When $d>\log n$, $\kappa_{p}(AR^{-1})=\mathcal{O}(d^{1+3\cdot|1/p-1/2|})$ and hence $\bar{\kappa}_{p}(AR^{-1})=\mathcal{O}(d^{1+3\cdot|1/p-1/2|+\max\{1/p,1/2\}})$ by Lemma 6. Note that, even for the case when $p=1$, we have $\bar{\kappa}_{p}(AR^{-1})=\mathcal{O}(d^{7/2})$, which is slightly better than FCT2 (see Corollary 1). However, we have to solve a rounding problem of size $ns/t\times d$ in the step 2 of FastLpBasis, which requires storage and work depending on $n$.

Numerical Implementation and Empirical Evaluation

$A_{3}$, the leading submatrix of the SNP matrix used by Paschou et al. . The SNP matrix is of size $492516\times 2250$, from the Human Genome Diversity Panel and the HapMap Phase 3 dataset. See for more descriptions of the data.

$A_{4}$, the leading submatrix of the TinyImages matrix created by Torralba et al. . The original images are in RGB format. We convert them to grayscale intensity images, resulting a matrix of size $8e7\times 1024$.

To implement FCT1 and FCT2 for our empirical evaluations, we have to fix several parameters in Theorems 2 and 1, finding a compromise between theory and practice. We choose $r_{1}=\lceil 2d\log d\rceil$ except $r_{1}=2d$ for GT. We choose $s=\lceil 2d\log d\rceil$ and $t=2d^{2}$ for FCT1, and $s=2^{\lceil 2\log_{2}(2d\log d)\rceil}$ for FCT2. Although those settings don’t follow Theorems 2 and 1 very closely, they seem to be good for practical use. Since all the transforms are randomized algorithms that may fail with certain probabilities, for each test matrix and each transform, we take $50$ independent runs and show the first and the third quartiles of $\bar{\kappa}_{1}$ in Tables 2 and 3.

The true signal $x^{*}$ is a standard Gaussian vector.

Each row of the design matrix $A$ is a canonical vector, which means that we only estimate a single entry of $x^{*}$ in each measurement. The number of measurements on the $i$-th entry of $x^{*}$ is twice as large as that on the $(i+1)$-th entry, $i=1,\ldots,14$. We have $2.62$ billion measurements on the first entry while only $0.16$ million measurements on the last. Imbalanced measurements apparently create difficulties for sampling-based algorithms.

Since the problem is separable, we know that an optimal solution is simply given by the median of responses corresponding to each entry.

We first check the overall performance of these sampling algorithms, measured by relative errors in $1$-, $2$-, and $\infty$-norms. The results are shown in Table 4.

Since the algorithms are all randomized, we show the first and the third quartiles of the relative errors in $100$ independent runs. We see that CT clearly performs the best, followed by GT. UNIF works but it is about a magnitude worse than CT. NOCD is close to UNIF at the first quartile, but makes very large errors at the third. Without conditioning, NOCD is more likely to sample outliers because the response from a corrupted measurement is much larger than that from a normal measurement. However, those corrupted measurements contain no information about $x^{*}$, which leads to NOCD’s poor performance. UNIF treats all the measurements the same, but the measurements are imbalanced. Although we sample $100000$ measurements, the expected number of measurements on the last entry is only $3.05$, which downgrades UNIF’s overall performance.

Conclusion

References

Appendix A Proofs of Technical Cauchy Lemmas

The proof uses similar techniques to the bounds due to Indyk for sums of independent clipped half-Cauchy random variables. Fix $M>0$ (we will choose $M$ later) and define the events

and $F=\cap_{i\in[m]}F_{i}$. Note that $F\cap F_{i}=F$. Using the pdf of a Cauchy and because $\tan^{-1}x\leq x$, we have that:

By a union bound, $\hbox{\bf{Pr}}[F]\geq 1-{2m\over\pi M}$. Further, $\hbox{\bf{Pr}}[F|F_{i}]\hbox{\bf{Pr}}[F_{i}]=\hbox{\bf{Pr}}[F\cap F_{i}]=\hbox{\bf{Pr}}[F]$, hence $\hbox{\bf{Pr}}[F|F_{i}]=\hbox{\bf{Pr}}[F]/\hbox{\bf{Pr}}[F_{i}]$. We now bound \hbox{\bf{E}}\left[|C_{i}|\ \bigl{|}\ F\right]. First, observe that

Next, since $F_{i}\cap F=F$, we have that

Finally, by using the pdf of a Cauchy, \hbox{\bf{E}}\left[|C_{i}|\bigl{|}F_{i}\right]={{1\over\pi}\log(1+M^{2})}/{\hbox{\bf{Pr}}[F_{i}]}, and so

By Markov’s inequality and because $\hbox{\bf{Pr}}[X\geq\gamma t|\bar{F}]\leq 1$, we have:

A.2 Proof of Lemma 4 (Cauchy Lower Tail Inequality)

To bound the lower tail, we will use Lemma 1. By homogeneity, it suffices to prove the result for $\gamma=1$. Let $Z_{i}=\gamma_{i}\min(|C_{i}|,M)$. Clearly $Z_{i}\leq\gamma_{i}|C_{i}|$ and so defining $Z=\sum_{i}Z_{i}$, we have that $Z\leq X$ and $\hbox{\bf{Pr}}[X\leq 1-t]\leq\hbox{\bf{Pr}}[Z\leq 1-t]$. Thus, we have that

where the last step holds by Lemma 1 for $1-t<\hbox{\bf{E}}[Z]$. Using the distribution of the half-Cauchy, one can verify using standard techniques that by choosing $M\approx 1.6768$, $\hbox{\bf{E}}[Z_{i}]=\gamma_{i}$ and $\hbox{\bf{E}}[Z_{i}^{2}]\leq{3\over 2}\gamma_{i}^{2}$, so $\sum_{i}\hbox{\bf{E}}[Z_{i}]=1$ and $\sum_{i}\hbox{\bf{E}}[Z_{i}^{2}]\leq{3\over 2}\sum_{i}\gamma_{i}^{2}\leq{3\over 2\beta^{2}}$. It follows that $\hbox{\bf{Pr}}[Z\leq 1-t]\leq\operatorname{exp}\left(-t^{2}/{\textstyle{3\over\beta^{2}}}\right),$ and the result follows.

First, observe that ${\|DZx\|}_{1}=\sum_{i\in[n]}D_{ii}|Z_{(i)}x|$, and since $\hbox{\bf{E}}[D_{ii}]=1$, $\hbox{\bf{E}}[{\|DZx\|}_{1}]=\sum_{i\in[n]}|Z_{(i)}x|={\|Zx\|}_{1}$. Next, observe that

because when $\hat{p}_{i}=1$, that row must be sampled, and so does not contribute to the deviation. So, we only need to analyze the RHS of the above equation. From now on, we only consider those $i$ with $\hat{p}_{i}<1$, in which case $\hat{p}_{i}=s\cdot t_{i}$, where $t_{i}\geq a{\|Z_{(i)}\|}_{1}/{\|Z\|}_{1}$. Let $Q_{i}$ be the (positive) random variable $D_{ii}|Z_{(i)}x|$; either $Q_{i}=0$ or

where we defined $\gamma={1\over a}{\|Z\|}_{1}{\|x\|}_{\infty}$. We can also obtain a bound for $\sum_{\hat{p}_{i}<1}\hbox{\bf{Var}}[Q_{i}]$:

where, in the last inequality, we used the upper bound for $Q_{i}$ and we further upper bounded by summing over all $i\in[n]$. Let $Q=\sum_{i}Q_{i}$ with $Q_{i}\leq\gamma$; the standard Bernstein bound states that

Plugging in our bounds for $\sum_{i}\hbox{\bf{Var}}[Q_{i}]$ and $\gamma$, we deduce that

The lemma follows after some simple algebraic manipulations.

Appendix B Proof of Theorem 1 (Fast Cauchy Transform (FCT1))

Before presenting the proof, we describe the main idea. It follows a similar line of reasoning to , and it uses an “uncertainty principle” (which we state as Lemma 7 below).

The uncertainty principle we prove follows from the fact that the concatenation of the Hadamard matrix with the identity matrix is a dictionary of low coherence. For background, and similar arguments to those we use in Lemma 7 below, see Section 4 of . In particular, see Claim 4.1 and Lemma 4.2 of that section.

To prove the upper bound, we use the existence of a $(d,1)$-conditioned basis $U$ and apply $\Pi_{1}$ to this basis to show that ${\|\Pi_{1}Ux\|}_{1}$ cannot expand too much, which in turn means that ${\|\Pi_{1}Ax\|}_{1}$ cannot expand too much (for any $x$). To prove the lower bound, we show that the inequality holds with exponentially high probability, for a particular $y$; and we then use a suitable $\gamma$-net to obtain the result for all $y$.

We now proceed with the proof of Theorem 1. We will first prove an upper bound (Proposition 1) and then a lower bound (Proposition 2); the theorem follows by combining Propositions 1 and 2.

With probability at least $1-\delta$, for all $x\in\R^{d}$, $\lVert\Pi_{1}Ax\rVert_{1}\leq\kappa\lVert Ax\rVert_{1}$, where $\kappa=O({d\sqrt{s}\over\delta}\log(r_{1}d))$.

Let $U\in\R^{n\times d}$ be a $(d,1)$-conditioned basis (see Definition 2 below) for the column space of $A$, which implies that for some $Z\in\R^{d\times d}$ we can write $A=UZ$. Since $\lVert\Pi_{1}Ax\rVert_{1}\leq\kappa\lVert Ax\rVert_{1}$ if and only if $\lVert\Pi_{1}UZx\rVert_{1}\leq\kappa\lVert UZx\rVert_{1}$, it suffices to prove the proposition for $U$. By construction of $U$, for any $x\in\R^{d}$, $\lVert x\rVert_{\infty}\leq\lVert Ux\rVert_{1}$, and so

Thus it is enough to show that $\lVert\Pi_{1}U\rVert_{1}\leq\kappa$. We have

Since $\lVert G_{s}z_{i}\rVert_{1}\leq\sqrt{2s}\lVert G_{s}z_{i}\rVert_{2}\leq\sqrt{4s}\lVert z_{i}\rVert_{2}\leq\sqrt{4s}\lVert z_{i}\rVert_{1}$, it follows that

Applying this to $y=U^{(j)}$ for $j\in[d]$ yields

since $\lVert U\rVert_{1}\leq d$ because $U$ is $(d,1)$-conditioned.

The $(i,j)$ entry of $BC\hat{U}$ is $\sum_{k}B_{ik}C_{kk}{\hat{U}}_{kj},$ which is a Cauchy scaled by $\gamma_{ij}=\sum_{k}|B_{ik}\hat{U}_{kj}|$. So,

Hence, we can apply Lemma 3 with $\gamma={\|\hat{U}\|}_{1}$ and $m=r_{1}d$ to obtain

Setting the RHS to $\delta$, it suffices that $t=O({1\over\delta}\log(dr_{1}))$. Thus, with probability at least $1-\delta$,

Before we prove the lower bound, we need the following lemma which is derived using a sparsity result for matrices with unit norm rows and low “coherence,” as measured by the maximum magnitude of the inner product between distinct rows ($G_{s}$ is a matrix with low coherence). This result mimics results in .

For $G=\left[\begin{smallmatrix}H_{s}\\ I_{s}\end{smallmatrix}\right]$ and any $z\in\R^{s}$, $\lVert Gz\rVert_{1}\geq{1\over 2}s^{1/4}\lVert z\rVert_{2}$.

We can assume $\lVert z\rVert_{2}=1$, and so $\lVert Gz\rVert_{2}^{2}=2$. Let $G_{(S^{\prime})}$ be $k$ rows of $G$, with $\kappa$ of them coming from $H_{s}$ and $k-\kappa$ from $I_{s}$. $G_{(S^{\prime})}G_{(S^{\prime})}^{T}=I+\Lambda$ where $\Lambda$ is a symmetric $2\times 2$ block matrix $\left[\begin{smallmatrix}\bm{0}&{1\over\sqrt{s}}Q\\ {1\over\sqrt{s}}Q^{T}&\bm{0}\\ \end{smallmatrix}\right]$ where the entries in $Q\in\R^{\kappa\times(k-\kappa)}$ are $\pm 1$, and so $\lVert Q\rVert_{2}\leq\sqrt{\kappa(k-\kappa)}\leq{1\over 2}k$.

Now, given any $z$, we set $k=2\beta\sqrt{s}$ with $\beta={2\over 5}$, and choose $G_{(S^{\prime})}$ to be the rows corresponding to the $k$ components of $Gz$ having largest magnitude, with $G_{(S)}$ being the rows with indices in $[s]\setminus S^{\prime}$. Then $\lVert G_{(S^{\prime})}z\rVert_{2}^{2}\leq 1+\beta$, and so the entry in $G_{(S^{\prime})}z$ with smallest magnitude has magnitude at most $a=\sqrt{(1+\beta)/k}=\sqrt{(1+\beta)/2\beta}s^{-1/4}$. We now consider $G_{(S)}z$. Since $\lVert Gz\rVert_{2}^{2}=2$, $\lVert G_{(S)}z\rVert_{2}^{2}\geq 1-\beta$; further, all components have magnitude at most $a$ (as all the components of $G_{(S)}z$ have smaller magnitude than those of $G_{(S^{\prime})}z$). $\lVert G_{(S)}z\rVert_{1}$ is minimized by concentrating all the entries into as few components as possible. Since the number of non-zero components is at least $(1-\beta)/a^{2}=2\beta s^{1/2}(1-\beta)/(1+\beta)$, giving these entries the maximum possible magnitude results in

(where we used $\beta={2\over 5}$). We are done because $\lVert Gz\rVert_{1}\geq\lVert G_{(S)}z\rVert_{1}$ ∎

We now prove the lower bound. We assume that Proposition 1 holds for $\Pi_{1}$, which is true with probability at least $1-\delta$ for $\kappa$ as defined in Proposition 1. Then, by a union bound, both Propositions 1 and 2 hold with probability at least

($\delta$ and $\kappa$ are from Proposition 1). Since $s^{1/2}=r_{1}^{3}$, by choosing $r_{1}=\alpha d\log{d\over\delta}$ for large enough $\alpha$, the final probability of failure is at most $2\delta$, because $\kappa=O({d\sqrt{s}\over\delta}\log(r_{1}d))=O(\operatorname{poly}(d))$.

Assume Proposition 1 holds. Then, for all $x\in\R^{d}$, $\lVert\Pi_{1}Ax\rVert_{1}\geq\lVert Ax\rVert_{1}$ holds with probability at least

First we will show a result for fixed $y\in\R^{n}$, summarized in the next lemma.

$\displaystyle\hbox{\bf{Pr}}\left[{\|\Pi_{1}y\|}_{1}<2{\|y\|}_{1}\right]\leq\operatorname{exp}\left(-{{r_{1}\over 48}}\right)+\operatorname{exp}\left(-{s^{1/2}\over 8r_{1}^{2}}+\log r_{1}\right)$

Given this lemma, the proposition follows by putting a $\gamma$-net $\Gamma$ on the range of $A$ (observe that the range of $A$ has dimension at most $d$). This argument follows the same line as in Sections 3 and 4 of . Specifically, let $L$ be any fixed (at most) $d$ dimensional subspace of $\R^{n}$ (in our case, $L$ is the range of $A$). Consider the $\gamma$-net on $L$ with cubes of side $\gamma/d$. There are $(2d/\gamma)^{d}$ such cubes required to cover the hyper-cube ${\|y\|}_{\infty}\leq 1$; and, for any two points $y_{1},y_{2}$ inside the same $\gamma/d$-cube, $\lVert y_{1}-y_{2}\rVert_{1}\leq\gamma$. From each of the $\gamma/d$-cubes, select a fixed representative point which we will generically refer to as $y^{*}$; select the representative to have $\lVert y^{*}\rVert_{1}=1$ if possible. By a union bound and Lemma 8,

We will thus condition on the high probability event that ${\|\Pi_{1}y^{*}\|}_{1}\geq 2{\|y^{*}\|}$ for all $y^{*}$. For any $y\in L$ with ${\|y\|}_{1}=1$, let $y^{*}$ denote the representative point for the cube in which $y$ resides (${\|y^{*}\|}_{1}=1$ as well). Then ${\|y-y^{*}\|}\leq\gamma$.

where the last inequality holds using Proposition 1. By choosing $\gamma=1/\kappa$ and recalling that ${\|y^{*}\|}_{1}=1$, we have that ${\|\Pi_{1}y\|}_{1}\geq 1$, with probability at least

We have that ${\|g\|}_{2}^{2}=\sum_{i}{\|G_{s}z_{i}\|}_{2}^{2}=2\sum_{i}{\|z_{i}\|}_{2}^{2}=2{\|y\|}_{2}^{2}$, and

where the last inequality is because $B^{(i)}$ is a standard basis vector. To bound $\sum_{j}\gamma_{j}^{2}$, we will show that $\gamma_{j}$ is nearly uniform. Since $\gamma_{j}$ is a weighted sum of independent Bernoulli random variables (because $B_{ji}$ and $B_{jk}$ are independent for $i\not=k$), we can use Lemma 2 with $\xi_{i}=|g|_{i}$ and $1-p=1-1/r_{1}\leq 1$, and so $\sum_{i}\xi_{i}={\|g\|}_{1}$ and $\sum_{i}\xi_{i}^{2}={\|g\|}^{2}_{2}$; setting $t=1/r_{1}$ in Lemma 2:

By a union bound, none of the $\gamma_{j}$ exceed $2{\|g\|}_{1}/r_{1}$ with probability at most $r_{1}\operatorname{exp}\left(-{s^{1/2}}/{8r_{1}^{2}}\right)$. We assume this high probability event, in which case $\sum_{j}\gamma_{j}^{2}\leq 4{\|g_{1}\|}^{2}_{1}/r_{1}$. We can now apply Lemma 4 with $\beta^{2}=r_{1}/4$ and $t={1\over 2}$ to obtain

By a union bound, ${\|BCg\|}_{1}\geq{1\over 2}{\|g\|}_{1}$ with probability at least $1-\operatorname{exp}\left(-{{r_{1}\over 48}}\right)-\operatorname{exp}\left(-{s^{1/2}\over 8r_{1}^{2}}+\log r_{1}\right)$. Scaling both sides by $4$ gives the lemma. ∎

Appendix C Proof of Theorem 2 (Fast Cauchy Transform (FCT2))

We will need results from prior work, which we paraphrase in our notation.

holds with probability at least $1-c_{1}e^{-c_{2}k\varepsilon^{2}}$ (w.r.t. $G$), for global constants $c_{1},c_{2},c_{3}>0$.

Remark. This is the standard Johnson-Lindenstrauss property with the additional requirement on ${\|Gx\|}_{1}$. Essentially it says that $Gx$ is a nearly uniform, so that ${\|Gx\|}_{1}\approx\sqrt{s}{\|Gx\|}_{2}\approx\sqrt{s}{\|x\|}_{2}$.

Let $L$ be any (fixed) $d$ dimensional subspace of $\R^{t}$, and $G$ an $s\times t$ matrix sampled from a distribution having the MJLP property. Given $\varepsilon\in(0,{1\over 3}]$, let $s=36(k+{8d\over c_{3}\varepsilon}+\log(2c_{1}))/c_{2}\varepsilon^{2}=O({k\over\varepsilon^{2}}+{d\over\varepsilon^{3}})$. Then, with probability at least $1-e^{-k}$, for every $x\in\R^{t}$,

We also need a result on how the matrix of Cauchy random variables $C$ behaves when it hits a vector $y$. The next theorem is Theorem 5 of . For completeness and also to fix some minor errors in the proof of , we give a proof of Theorem 12 in Appendix I.

where $\kappa^{\prime}=O({d\over\delta}\log(r_{1}d))$.

Note that for $\delta$ fixed to some small error probability, $r_{1}=O(d\log d)$, and the product $Cy$ in the theorem above can be computed in time $O(r_{1}n)=O(nd\log d).$

The vector $z_{i}\in\operatorname{colsp}A_{(\{i\})}$, noting that $\operatorname{colsp}A_{(\{i\})}$ is a subspace of $\R^{t}$ of dimension at most $d$. Let $U_{i}\in\R^{t\times d}$ be an orthonormal basis for $\operatorname{colsp}A_{(\{i\})}$, and let $z_{i}=U_{i}w_{i}$. Setting $\varepsilon={1\over 2}$ in Lemma 10, and recalling that $G$ is $s\times t$, $k$ in Lemma 10 can be expressed as $k={c_{2}s\over 144}-{16d\over c_{3}}-\log(2c_{2})$. Applying a union bound, we have that for all $i\in[n/t]$ with probability at least $1-2c_{1}\cdot{n\over t}\cdot\operatorname{exp}(-{c_{2}s\over 144}+{16d\over c_{3}})$, that for all $y=Ax$ (and corresponding $z_{i}\in\R^{t}$), it holds that

holds with probability at least $1-\delta-2c_{1}\cdot{n\over t}\operatorname{exp}(-{c_{2}s\over 144}+{16d\over c_{3}})\geq 1-2\delta$ (by choosing $s\geq{144\over c_{2}}({16d\over c_{3}}+\log{2c_{1}n\over\delta t})$). The theorem follows because $\log n\leq d$ and hence $\kappa^{\prime}=O({d\over\delta}\log d)$, $s=O(d+\log{1\over\delta})$ and $t=O(s^{2+\eta})$.

Clearly $U=AR^{-1}$ is in the range of $A$ and has the same null-space otherwise $\Pi_{1}A$ would not preserve lengths to relative error. Therefore $U$ is a basis for the range of $A$. Consider any $x\in\R^{d}$. The first claim of the theorem follows from the following derivations:

(a) follows from the lower bound in (3), because it holds for every column of $AR^{-1}$; (b) follows because by the construction of $R$, $\Pi_{1}AR^{-1}$ has $d$ orthonormal columns; finally, (c) follows from the upper bound in (3).

Finally, to obtain the Corollary, if $\Pi_{1}$ satisfying (3) is constructed using Theorem 2 with small fixed probability of failure $\delta$, then $r_{1}=O(d\log d)$ and $\kappa=O(d^{2+\eta}\log d)$. The running time to compute ${R^{-1}}$ is obtained by summing $O(nd\log d)$ (to compute $\Pi_{1}A$) and $O(r_{1}d^{2})=O(d^{3}\log d)$ (to obtain $R^{-1}\in\R^{d\times d}$).

Let ${\cal C}^{\prime}=\{y={R^{-1}}x:x\in{\cal C}\}$ be a linear transform of the constraint set. We start with a basic lemma that allow us to use $U$ instead of $X$. This lemma says that if we can construct a sampling matrix $D$ for $U$ under the constraint ${\cal C}^{\prime}$ such that solving the down-sampled problem for $U$ gives a $(1+\varepsilon)$-approximation, then that same sampling matrix works for $X$ under the constraint ${\cal C}$.

Let $U=X{R^{-1}}$ and $D$ any diagonal sampling matrix as in Figure 2. Suppose that for any $\hat{y}\in{\cal C}^{\prime}$ that minimizes ${\|DUy\|}$, $\hat{y}$ is a $(1+\varepsilon)$-approximation for the problem $\min_{y\in{\cal C}^{\prime}}{\|Uy\|}$. Let $\hat{x}$ be any solution to $\min_{x\in{\cal C}}{\|DXx\|}$. Then $\hat{x}$ is a $(1+\varepsilon)$-approximation for the problem $\min_{x\in{\cal C}}{\|Xx\|}$.

Select $\hat{y}={R^{-1}}\hat{x}\in{\cal C}^{\prime}$. For any $y\in{\cal C}^{\prime}$, there is some $x\in{\cal C}$ with $y={R^{-1}}x$, and we have:

where (a) is by the optimality of $\hat{x}$. So, $\hat{y}$ minimizes ${\|DUy\|}_{1}$, hence for all $y\in{\cal C}^{\prime}$, ${\|U\hat{y}\|}_{1}\leq(1+\varepsilon){\|Uy\|}_{1}$. Now consider any $x\in{\cal C}$ and let $y={R^{-1}}x\in{\cal C}^{\prime}$. Then,

By Theorem 3, $U=X{R^{-1}}$ is an $(\alpha,\beta)$-conditioned basis for the range of $X$, where

So, ${\|U\|}_{1}\leq q\sqrt{r_{1}}$ and for all $y\in\R^{q}$, ${\|y\|}_{\infty}\leq\kappa{\|Uy\|}_{1}$. We next show that $\lambda_{i}$ estimates ${\|U_{(i)}\|}_{1}$. The following lemma is a straightforward application of a Chernoff bound to independent half-Cauchys (see also Claims 1 and 2 and Lemmas 1 and 2 in ).