\mathcal B^{M,\bfw}_\rho \coloneqq \braces*{\bfy \in \mathcal V: \norm*{\bfy}_{M,\bfw} \leq \rho}.
\]

Similarly define the sphere $\mathcal S^{M,\bfw}_\rho$ of radius $\rho>0$ as

\[
\mathcal S^{M,\bfw}_\rho \coloneqq \braces*{\bfy \in \mathcal V: \norm*{\bfy}_{M,\bfw} = \rho}.
\]

If $\bfw = \mathbf{1}_n$, then we simply write $\mathcal B_\rho^{M}$ and $\mathcal S_\rho^{M}$, respectively.

Additional useful properties that we will need are included in Appendix <ref>, including sufficient conditions for triangle inequality and net constructions.

§.§ Sensitivities for M-Estimators

Because $M$-estimators are defined as coordinate-wise sums, one can naturally define analogues of sensitivities, just as was done for $\ell_p$ norms.

Let $\bfA\in\mathbb R^{n\times d}$ and let $\norm*{\cdot}_M$ be an $M$-norm. Then, the $i$th $M$-sensitivity is defined as

\[
\bfs_i^M(\bfA) \coloneqq \sup_{\bfx\in\mathbb R^d, \bfA\bfx\neq 0} \frac{M(\abs*{[\bfA\bfx](i)})}{\norm*{\bfA\bfx}_M^{p_M}}
\]

and the total $M$-sensitivity is defined as

\[
\mathcal T^M(\bfA)\coloneqq \sum_{i=1}^n \bfs_i^M(\bfA).
\]

Let $\bfw\geq\mathbf{1}_n$ be a set of weights. Then, the $i$th weighted $M$-sensitivity is defined as

\[
\bfs_i^{M,\bfw}(\bfA) \coloneqq \sup_{\bfx\in\mathbb R^d, \bfA\bfx\neq 0} \frac{\bfw_i M(\abs*{[\bfA\bfx](i)})}{\norm*{\bfA\bfx}_{M,\bfw}^{p_M}}
\]

and the total weighted $M$-sensitivity is defined as

\[
\mathcal T^{M,\bfw}(\bfA)\coloneqq \sum_{i=1}^n \bfs_i^{M,\bfw}(\bfA).
\]

When $M(x) = \abs*{x}^p$, i.e. for the case of $\ell_p$ norms, it is known that sampling with probabilities proportional to upper bounds on sensitivities yields subspace embeddings [6, 24, 17]. Analogous results are known as well for $M$-estimators [22, 21, 23] and Orlicz norms [61].

Let $\bfA\in\mathbb R^{n\times d}$, let $\norm*{\cdot}_M$ be an $M$-norm, and let $\bfw\geq\mathbf{1}_n$ be a set of weights. Let $m$ be an oversampling parameter. Then, a random set of weights $\bfw'$ is sampled according to sensitivity upper bounds $\tilde\bfs_i^{M,\bfw}(\bfA) \geq \bfs_i^{M,\bfw}(\bfA)$ (see Definition <ref>) if

\[
\bfw_i' \coloneqq \begin{cases}
\bfw_i/\bfp_i & \text{w.p. $\bfp_i$} \\
0 & \text{otherwise}
\end{cases}
\]

where $\bfp_i \coloneqq \min\{1,m\cdot \tilde\bfs_i^{M,\bfw}(\bfA)\}$.

Note that in the case of $M$-estimators, the lack of scale invariance means that we get norm preservation guarantees for spheres rather than for entire subspaces. That is, we can get the following lemma, similar Lemma 43 of [21]:

Let $\bfA\in\mathbb R^{n\times d}$. Let $\eps\in(0,1)$, $\delta>0$, and let $\rho \geq 1$. Let $M:\mathbb R_{\geq0}\to\mathbb R_{\geq0}$ satisfy the conditions of Definition <ref>, and furthermore that

* $M^{1/p_M}$ is subadditive

* $M$ is polynomially bounded below with degree $q_M$ and constant $c_L$ (see Definition <ref>)

Let $\bfw\geq\mathbf{1}_n$ be a set of weights. Let $\tilde\bfs_i^{M,\bfw}(\bfA)\geq\bfs_i^{M,\bfw}(\bfA)$ be sensitivity upper bounds. Let $m\geq m_0$ be an oversampling parameter larger than some

\[
m_0 = O\parens*{\frac{d}{\eps^2}\parens*{\log\frac1\eps}\parens*{\log\frac1\delta}}.
\]

Let $\bfw'\geq \mathbf{1}_n$ be sampled according to Definition <ref>. Then with probability at least $1 - \delta$,

\[
\norm*{\bfy}_{M,\bfw'}^{p_M} = (1\pm\eps)\norm*{\bfy}_{M,\bfw}^{p_M}
\]

for all $\bfy\in\mathcal S_\rho^M$. Furthermore,

\[
\E\nnz(\bfw') \leq m\sum_{i=1}^n\tilde\bfs_i^{M,\bfw} = m\tilde{\mathcal T}^{M,\bfw}(\bfA)
\]

The proof is by a standard Bernstein bound, and is included in Appendix <ref>.

§.§ Efficient Sensitivity Bounds

We first show that algorithmically, one can compute upper bounds to the $M$-estimator sensitivities that sum to at most $O(d^{\max\{1,p_M/2\}}\log^2 n + \tau)$ in time

\[
O\parens*{\nnz(\bfA)\log^3 n + \frac{nT}{\tau}\log n},
\]

where $T = T(n,d)$ is such that constant factor $\ell_{p_M}$ Lewis weight approximation for an $n\times d$ matrix $\bfB$ takes $O(\nnz(\bfB)\log n + T)$ time. For example, it is known that $\ell_p$ Lewis weights for $0 < p < 4$ can be approximated up to constant factors in $O(\nnz(\bfA)\log n + d^{\omega})$ time, so for $\tau = T = d^{\omega}$, we obtain a nearly input sparsity time algorithm that computes upper bounds to $M$-estimator sensitivities that sum to at most $O(d^{\max\{1,p_M/2\}}\log^2 n + d^{\omega})$. Note that in applications (see applications to active regression (Section <ref>), Orlicz norm regression (Section <ref>), and robust subspace approximation (Section <ref>) in this work), this is enough to compute a set of $\poly(d)\log^2 n$ rows that approximates the original matrix well, at which point we can compute sensitivities that sum to only $O(d\log^2 n)$ in an additional $\poly(d\log n)$ time.

The algorithm draws ideas from a theorem of <cit.>, which shows an input sparsity time algorithm for locating “heavy entries” for the Tukey loss, which is equivalent to finding coordinates with high Tukey sensitivity.

Sensitivity upper bounds
input: Matrix $\bfA \in \mathbb R^{n \times d}$, $M$-norm $M$, parameter $1\leq\tau\leq n$.

output: Upper bounds $\tilde\bfs_i^M(\bfA)$ on $\bfs_i^M(\bfA)$.
[1]
Initialize $\tilde\bfs_i^M(\bfA) \leftarrow 2\tau/n$.
$t\in[O(\log n)]$
Hash the rows of $\bfA$ into $B = 10\cdot 2^r$ buckets $S_1, S_2, \dots, S_B$.
Compute $O(1)$-approximate $\ell_p$ Lewis weights of each of the $B$ buckets $\bfA\mid_{S_1}, \bfA\mid_{S_2}, \dots, \bfA\mid_{S_B}$.
For any row $i$ with an $\ell_{p_M}$ Lewis weight of at least $\Omega(1)$, set $\tilde\bfs_i^M(\bfA) \leftarrow \max\{2/2^r, \tilde\bfs_i^M(\bfA)\}$.

[Main Result for Sensitivity Bounds]theoremSensitivityBound
Let $\norm*{\cdot}_M$ be an $M$-norm. Let $1\leq \tau\leq n$ be a parameter. Then, with probability at least $99/100$, Algorithm <ref> computes sensitivity upper bounds $\tilde\bfs_i^M(\bfA) \geq \bfs_i^M(\bfA)$ that sum to at most

\[
\tilde{\mathcal T}^M(\bfA) \coloneqq \sum_{i=1}^n \tilde\bfs_i^M(\bfA) = O(d^{\max\{1,p_M/2\}} \log^2 n + \tau).
\]

If constant factor Lewis weight approximation takes for an $n\times d$ matrix $\bfB$ takes $O(\nnz(\bfB)\log n + T)$ time (see Theorem <ref>), then the total running time is

\[
O\parens*{\nnz(\bfA)\log^3 n + \frac{nT}{\tau}\log n}.
\]

We first show correctness of the algorithm, then show the sensitivity bound, and finally the running time guarantee.

Let $r\in[\ceil{\log_2 (n/\tau)}]$. Consider a coordinate $i\in[n]$ that has $M$-sensitivity between $1/2^r$ and $2/2^r$ and let $\bfy = \bfA\bfx$ be a corresponding vector which satisfies

\[
\frac{M(|\bfy_i|)}{\norm*{\bfy}_M^{p_M}} \in \left[\frac1{2^r},\frac2{2^r}\right].
\]

Note then that there are at most $2^r-1$ entries $j$ of $\bfy$ such that $\abs*{\bfy_j} > \abs*{\bfy_i}$. In our algorithm (line <ref>), we randomly hash the $n$ rows of $\bfA$ into $B = 10\cdot 2^r$ buckets. Then, the probability that any one of these $2^r-1$ entries is hashed to the same bucket as $i$ is $1/B$, so by a union bound, the probability that $\bfy_i$ has the largest absolute value in its hash bucket is at least $1 - 2^r / B \geq 9/10$. Call this event $\mathcal E$.

Now let $S$ be the set of indices which hash to the same bucket as row $i$ and let $S' = S\setminus\{i\}$. By Markov's inequality, with probability at least $9/10$, $i$ is hashed to a bucket such that the $M$-norm of all other entries is most $\norm*{\bfy\mid_{S'}}_M^{p_M} \leq \norm*{\bfy}_M^{p_M}/2^r$, where $\bfy\mid_{S'}$ is the restriction of $\bfy$ to the indices in $S'$. Call this event $\mathcal F$.

Condition on $\mathcal E$ and $\mathcal F$. We have by Lemma <ref> that

\[
\frac{|\bfy_i|^{p_M}}{\norm*{\bfy\mid_{S}}_{p_M}^{p_M}} \geq \frac{M(|\bfy_i|^{p_M})}{\norm*{\bfy\mid_{S}}_M^{p_M}} \geq \frac{M(|\bfy_i|^{p_M})}{M(|\bfy_i|^{p_M}) + \norm*{\bfy\mid_{S'}}_M^{p_M}} \geq \frac{1/2^r}{2/2^r + 1/2^r} = \frac13.
\]

The above holds with probability at least $4/5$. Thus, by repeating the hashing process $O(\log n)$ times, with probability at least $1 - 1/(100 n\log_2 n) = 1 - 1/\poly(n)$, there exists some trial where the $\ell_{p_M}$ sensitivity of the $i$th row in the matrix $\bfA\mid_S$ is at least $1/3$. In this trial, our algorithm will correctly set $\tilde\bfs_i^M(\bfA) \geq 2/2^r$ (line <ref>). By a union bound over $O(\log (n/\tau))$ levels $r$ and the $n$ rows, our algorithm succeeds with probability at least $99/100$.

Sensitivity Bound.

By Lemma <ref>, the $\ell_{p_M}$ sensitivities sum to at most $d^{\max\{1,{p_M}/2\}}$. Thus, each time we compute $O(1)$-approximate $\ell_{p_M}$ Lewis weights (line <ref>), we find at most $O(d^{\max\{1,{p_M}/2\}})$ entries with $M$-sensitivity at least $2/2^r$. Thus, for each $r$ and each iteration, we increase the sum of our upper bounds on $M$-sensitivities by a total of at most

\[
B \cdot O(d^{\max\{1,{p_M}/2\}}) \cdot \frac2{2^r} = O(d^{\max\{1,{p_M}/2\}}).
\]

This occurs at most $O((\log n)(\log(n/\tau)))$ times, and we start at a sensitivity bound of

\[
\frac{2\tau}{n}\cdot n = O(\tau)
\]

so our upper bounds on the sensitivities sum to at most

\[
O(d^{\max\{1,{p_M}/2\}}(\log n)(\log(n/\tau)) + \tau) \leq O(d^{\max\{1,{p_M}/2\}}\log^2 n + \tau).
\]

Running Time.

For a given $r$, the dominating running time cost of the inner-most loop of Algorithm <ref> is the computation of $\ell_{p_M}$ Lewis weights for $O(2^r)$ matrices whose sparsities sum to $\nnz(\bfA)$. Thus, if Lewis weight computation for an $n\times d$ matrix $\bfB$ takes $O(\nnz(\bfB)\log n + T)$ time, then the total running time is

\[
\sum_{r=1}^{\ceil{\log_2(n/\tau)}}O(\log n)\cdot O\parens*{\nnz(\bfA)\log n + 2^r T} = O\parens*{\nnz(\bfA)\log^3 n + \frac{nT}{\tau}\log n}.\qedhere
\]

As noted by [66], if we can control the sensitivities of functions $M_1$ and $M_2$, then it is straightforward to control the sensitivities of the sum of these two functions, i.e., $M = M_1 + M_2$. This applies to our algorithm as well. Suppose that $\bfs_i^M(\bfA) \in [1/2^r, 2/2^r]$ and let $\bfy = \bfA\bfx$ be such that

\[
\frac{M_1(\abs*{\bfy_i}) + M_2(\abs*{\bfy_i})}{\norm*{\bfy}_{M_1}^{p_{M_1}} + \norm*{\bfy}_{M_2}^{p_{M_2}}} \in \bracks*{\frac1{2^r}, \frac2{2^r}}.
\]

\[
\frac{M_1(\abs*{\bfy_i})}{\norm*{\bfy}_{M_1}^{p_{M_1}}} + \frac{M_2(\abs*{\bfy_i})}{\norm*{\bfy}_{M_2}^{p_{M_2}}} \geq \frac1{2^r}
\]

so then there is some $j\in\{1,2\}$ such that

\[
\frac{M_j(\abs*{\bfy_i})}{\norm*{\bfy}_{M_j}^{p_{M_j}}} \geq \frac12\cdot \frac1{2^r},
\]

so the $\ell_{p_j}$ Lewis weight of the $i$th coordinate must be large by using a similar proof as Theorem <ref>. Thus, we can obtain similar sensitivity upper bounds up to a constant factor loss. Similarly, if $M_2$ is a “flat” sensitivity function in the sense of [66], that is, if

\[
\sup_{\bfx\in\mathbb R^d}\frac{M_2(\abs*{[\bfA\bfx](i)})}{\norm*{\bfA\bfx}_{M_2}^{p_{M_2}}} = O\parens*{\frac1n}
\]

for all $i\in[n]$, then this just means that either the $M_1$ sensitivity is large, or the $M$ sensitivity is at most $O(1/n)$, in which case we still get the same bounds.

Although Theorem <ref> only handles unweighted $M$-estimators, this result can be generalized to weighted $M$-estimators by splitting into level sets, similarly to Lemma 39 of [21].

Let $\norm*{\cdot}_M$ be an $M$-norm. Let $\bfw\geq\mathbf{1}_n$ be a set of weights. Let $N\coloneqq \ceil*{\log_2(1+\norm*{\bfw}_\infty)}$. For $j\in[N]$, let

\[
T_j \coloneqq \braces*{i\in[n] : 2^{j-1}\leq \bfw_i < 2^j},
\]

and let $\bfA\mid_{T_j}$ denote the restriction of $\bfA$ to the rows of $T_j$. Then,

\[
\bfs_i^{M,\bfw}(\bfA) \leq 2\cdot \bfs_i^M(\bfA\mid_{T_j})
\]

for $i\in T_j$.

Let $i\in T_j$ for some $j\in[N]$. We have that

\begin{align*}
\bfs_i^{M,\bfw}(\bfA) &= \sup_{\bfx\in\mathbb R^n, \bfA\bfx\neq 0} \frac{\bfw_i M([\bfA\bfx](i))}{\norm*{\bfA\bfx}_{M,\bfw}^{p_M}} \\
&\leq \sup_{\bfx\in\mathbb R^n, \bfA\bfx\neq 0} \frac{\bfw_i M([\bfA\mid_{T_j}\bfx](i))}{\norm*{\bfA\mid_{T_j}\bfx}_{M,\bfw}^{p_M}} \\
&\leq \sup_{\bfx\in\mathbb R^n, \bfA\bfx\neq 0} \frac{2^j M([\bfA\mid_{T_j}\bfx](i))}{2^{j-1}\norm*{\bfA\mid_{T_j}\bfx}_{M}^{p_M}} \\
&= 2\cdot \bfs_i^{M}(\bfA\mid_{T_j})
\end{align*}

as desired.

This leads to an algorithm that achieves guarantees similar to Theorem <ref> for weighted $M$-sensitivities, up to a loss of a factor of $N$ in the running time and sensitivity bound.

Let $\norm*{\cdot}_M$ be an $M$-norm. Let $\bfw\geq\mathbf{1}_n$ be a set of weights. Define $N$ as in Lemma <ref>. There is an algorithm that computes weighted $M$-estimator sensitivities that sum to at most

\[
O(Nd^{\max\{1,p_M/2\}}\log^2 n + N\tau)
\]

in time

\[
O\parens*{\nnz(\bfA)\log^3 n + N\frac{nT}{\tau}\log n},
\]

where $T$ is such that constant factor Lewis weight approximation for an $n\times d$ matrix $\bfB$ takes $O(\nnz(\bfB) + T)$ time (see Theorem <ref>).

This is simply the result of applying Theorem <ref> on the $N$ matrices $\bfA\mid_{T_j}$ as defined in Lemma <ref>. Note that the $\nnz(\bfA\mid_{T_j})$ terms add up to $\nnz(\bfA)$ in the running time.

§.§ Sharper Sensitivity Bounds

We show that we may modify the proof of our input sparsity time algorithm to show that the sum of sensitivities is at most $O(d^{\max\{1,p_M/2\}}\log n)$, if we do not need to efficient algorithms for constructing these sensitivities.

Let $\norm*{\cdot}_M$ be an $M$-norm. Then, the total $M$-sensitivity of $\bfA$ is at most

\[
O(d^{\max\{1,p_M/2\}}\log n).
\]

Our idea is essentially to run Algorithm <ref> with $\tau = d$ without the $O(\log n)$ repetitions of the hashing process.

Let $r\in[\ceil{\log_2 n}]$ and let $I_r$ be the set of coordinates with $M$-sensitivity in $[1/2^r, 2/2^r]$. Suppose we hash the rows of $\bfA$ into $B = 10\cdot 2^r$ buckets. Then, as in the proof of Theorem <ref>, for each $i\in I_r$, there is at least a $9/10$ probability that $i$ has $\ell_{p_M}$ Lewis weight at least $1/3$ in its hash bucket. Thus, the number of such $i$ is $(9/10) \abs*{I_r}$ in expectation, so there exists some hashing such that at least $(9/10) \abs*{I_r}$ of the indices $i\in I_r$ have $\ell_p$ Lewis weight at least $1/3$ in its hash bucket. However, there can be at most $B\cdot d^{\max\{1,p/2\}}$ such indices, so we must have that

\[
\frac{9}{10}\abs*{I_r} \leq B \cdot d^{\max\{1,p/2\}}
\]

\[
\abs*{I_r} = O(B \cdot d^{\max\{1,p/2\}}) = O(2^r\cdot d^{\max\{1,p/2\}}).
\]

By summing over the $r$, we obtain a bound of

\[
\sum_{r=1}^{\ceil{\log_2 n}} \frac2{2^r}\abs*{I_r} \leq \sum_{r=1}^{\ceil{\log_2 n}}\frac2{2^r}O(2^r \cdot d^{\max\{1,p/2\}}) = O(d^{\max\{1,p/2\}}\log n) = O(d^{\max\{1,p/2\}}\log n)
\]

on the total $M$-sensitivity, as claimed.

§.§ Sensitivity Lower Bounds

Finally, we show that our sensitivity upper bounds are tight by showing that the Tukey loss can have a total sensitivity as large as $\Omega(d\log(n/d))$. We also show a weaker lower bound of $\Omega(d\log\log(n/d))$ for the Huber loss. This is in contrast to sensitivities for the $\ell_p$ loss for $0 < p < \infty$, which is always at most $d^{\max\{1,p/2\}}$ due to the existence of Lewis bases [42, 63], and thus has no dependence on $n$. The necessity for a dependence on $n$ can be attributed to the lack of scale invariance for these $M$-estimator losses. A similar observation has been made previously in <cit.>, which shows that the column subset selection problem with the entrywise Huber loss exhibits a lower bound of $\Omega(\sqrt{\log n})$ columns, also attributed to the lack of scale invariance.

We simultaneously handle the Tukey and Huber losses by analyzing the $\ell_2$-$\ell_p$ loss for $p\in[0,1]$, which grows quadratically near the origin and as $\ell_p$ away from the origin, and is polynomially bounded above with degree $2$.

Define the $\ell_2$-$\ell_p$ loss of width $\tau$ to be

\[
M(x) = \begin{cases}
x^2 & \abs*{x} \leq \tau \\
(\tau^2 / \tau^p) \cdot x^p & \abs*{x} > \tau
\end{cases}.
\]

For $d \geq 1$ and $n\geq d$, there exists an $n\times d$ matrix $\bfA$ with total $M$-sensitivity that is at least

\[
\mathcal T^T(\bfA) \geq \begin{cases}
\Omega\parens*{d\log\frac{n}{d}} & \text{if $p \in [0, 1)$} \\
\Omega\parens*{d\log\log\frac{n}{d}} & \text{if $p = 1$}.
\end{cases}
\]

Let $\ell = \floor*{\log_2 n}$ and let $\bfx\in\mathbb R^n$ be a vector with $2^i$ coordinates of value $\tau/2^i$ for $i\in[\ell]$. We will show a sensitivity lower bound of $\Omega(\ell) = \Omega(\log n)$ for the $n\times 1$ matrix formed by the vector $\bfx$. By considering $d$ disjoint copies of this vector, each on $n/d$ coordinates, this implies a lower bound of $\Omega(d\log(n/d))$.

Let $j\in[\ell]$. Then,

\begin{align*}
\norm*{2^j\cdot \bfx}_M^2 &= \sum_{i=1}^\ell 2^i \cdot M\parens*{\tau\frac{2^j}{2^i}} \\
&\leq \frac{\tau^2}{\tau^p}\sum_{i=1}^j 2^i\cdot\parens*{\tau\frac{2^j}{2^i}}^p + \sum_{i=j+1}^\ell 2^i\cdot  \parens*{\tau\frac{2^j}{2^i}}^2 \\
&= \tau^2 2^{pj}\sum_{i=1}^j 2^{(1-p)i} + \tau^2 2^{2j}\sum_{i=j+1}^\ell \frac1{2^i} \\
&= \begin{cases}
O(\tau^2 \cdot2^j) & \text{if $p\in[0,1)$} \\
O(\tau^2 \cdot j2^j) & \text{if $p = 1$} \\
\end{cases}
\end{align*}

so for each $j\in[\ell]$, there are $2^j$ coordinates $i$ such that

\begin{align*}
\frac{M(2^j\cdot\bfx_i)}{\norm*{2^j\cdot \bfx}_M^2} &= \begin{cases}\Omega\parens*{\frac{\tau^2}{\tau^2 \cdot 2^j}} & \text{if $p\in[0,1)$} \\
O\parens*{\frac{\tau^2}{\tau^2 \cdot j2^j}} & \text{if $p = 1$} \\
\end{cases} \\
&= \begin{cases}\Omega\parens*{\frac{1}{2^j}} & \text{if $p\in[0,1)$} \\
O\parens*{\frac{1}{j2^j}} & \text{if $p = 1$} \\
\end{cases}.
\end{align*}

Thus, the sum of sensitivities for the Tukey loss for this matrix is at least

\[
\sum_{j=1}^\ell 2^j \cdot \Omega\parens*{\frac1{2^j}} = \Omega\parens*{\ell} = \Omega\parens*{\log n}
\]

for $p\in[0,1)$ and

\[
\sum_{j=1}^\ell 2^j \cdot \Omega\parens*{\frac1{j2^j}} = \Omega\parens*{\log\ell} = \Omega\parens*{\log\log n}
\]

for $p=1$.

§ APPLICATIONS: ACTIVE REGRESSION FOR M-ESTIMATORS

We now show that our algorithmic ideas for $\ell_p$ active regression may be applied to handle active regression for a general class of $M$-estimators, using our sensitivity results in Section <ref>.

Throughout this section, we assume the following. Let $\bfA\in\mathbb R^{n\times d}$ and let $\mathcal V = \Span(\bfA)$. Let $M:\mathbb R_{\geq0}\to\mathbb R_{\geq0}$ satisfy the conditions of Definition <ref>, and furthermore that

* $M^{1/p_M}$ is subadditive

* $M$ is polynomially bounded below with degree $q_M$ and constant $c_L$ (see Definition <ref>)

In the linear regression problem with the $M$-loss, we are given a matrix $\bfA\in\mathbb R^{n\times d}$ and a vector $\bfb\in\mathbb R^n$ we must solve the following optimization problem:

\[
\min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx - \bfb}_M.
\]

In the active setting, we wish to solve the above problem while querying as few entries of $\bfb$ as possible. We will also consider the weighted version, i.e.

\begin{equation}\label{eqn:m-linear-regression}
\min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx - \bfb}_{M,\bfw}
\end{equation}

for a set of weights $\bfw\geq\mathbf{1}_n$. Throughout this section, let for a given weighted $M$-norm $\norm*{\cdot}_{M,\bfw}$, let

\[
\OPT \coloneqq \min_\bfx\norm*{\bfA\bfx-\bfb}_{M,\bfw}.
\]

§.§ Constant Factor Approximation

We adapt Lemmas 7 and 8 of [24] to first obtain a constant factor solution $\bfx_c$ to the regression problem.

Let $\bfw\geq\mathbf{1}_n$ be a set of weights. Let $\bfw'\geq\mathbf{1}_n$ be a random set of weights that satisfies the following:

* $\E_{\bfw'}\norm*{\bfy}_{M,\bfw'} = \norm*{\bfy}_{M,\bfw}$ for any fixed $\bfy\in\mathbb R^n$

* $\bfw'$ is a $1/8$-subspace embedding for any $\mathcal S_\rho^M$ for any fixed $\rho>0$ with probability at least $9/10$, that is,

\[
\Pr\braces*{\norm*{\bfA\bfx}_{M,\bfw'} = \parens*{1\pm\frac18}\norm*{\bfA\bfx}_{M,\bfw}, \forall \bfA\bfx\in\mathcal S_\rho^M} \geq \frac{9}{10}
\]

Let $\tilde\bfx$ satisfy

\[
\norm*{\bfA\tilde\bfx - \bfb}_{M,\bfw'} \leq \kappa\min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx - \bfb}_{M,\bfw'}.
\]

Then, with constant probability,

\[
\norm*{\bfA\tilde\bfx - \bfb}_{M,\bfw} \leq O(\kappa)\min_{\bfx\in\mathbb R^d} \norm*{\bfA\bfx - \bfb}_{M,\bfw}.
\]

\[
\OPT \coloneqq \min_{\bfx\in\mathbb R^d} \norm*{\bfA\bfx - \bfb}_{M,\bfw}
\]

and let $\bfx^*$ be the minimizer achieving this value. Then by assumption, we have that

\[
\E_{\bfw'} \norm*{\bfA\bfx^* - \bfb}_{M,\bfw'} = \OPT
\]

so by Markov's inequality, with probability at least $9/10$, $\norm*{\bfA\bfx^* - \bfb}_{M,\bfw'}\leq 10\cdot\OPT$. Condition on this event. Note then that

\begin{equation}\label{eqn:m-regression-opt-value}
\norm*{\bfA\tilde\bfx - \bfb}_{M,\bfw'} \leq \kappa\norm*{\bfA\bfx^* - \bfb}_{M,\bfw'} \leq 10\kappa\cdot\OPT.
\end{equation}

Now suppose for contradiction that

\[
\norm*{\bfA\tilde\bfx - \bfb}_{M,\bfw} > \frac{25\kappa}{c_L^{1/p_M}}\cdot\OPT.
\]

\begin{align*}
\norm*{\bfA\tilde\bfx - \bfA\bfx^*}_{M,\bfw} &\geq \norm*{\bfA\tilde\bfx - \bfb}_{M,\bfw} - \norm*{\bfA\bfx^* - \bfb}_{M,\bfw} && \text{triangle inequality} \\
&> \frac{25\kappa}{c_L^{1/p_M}}\cdot\OPT - \OPT = \frac{24\kappa}{c_L^{1/p_M}}\cdot\OPT.
\end{align*}

We now apply the subspace embedding property, with $\rho = (24\kappa/c_L^{1/p_M})\cdot\OPT$ so that

\[
\norm*{\bfA\bfx}_{M,\bfw'} \geq \parens*{1\pm\frac18} \norm*{\bfA\bfx}_{M,\bfw}
\]

for all $\bfA\bfx$ with $\norm*{\bfA\bfx}_{M,\bfw} = (24\kappa/c_L^{1/p_M})\cdot\OPT$. We then apply Lemma <ref> to see that

\[
\norm*{\bfA\bfx}_{M,\bfw'} \geq c_L^{1/p_M} \frac78 \cdot \frac{24\kappa}{c_L^{1/p_M}}\cdot\OPT = 21\kappa\cdot\OPT
\]

for all $\bfA\bfx$ with $\norm*{\bfA\bfx}_{M,\bfw} \geq (24\kappa/c_L^{1/p_M})\cdot\OPT$. Thus,

\begin{align*}
\norm*{\bfA\tilde\bfx - \bfb}_{M,\bfw'} &\geq \norm*{\bfA\tilde\bfx - \bfA\bfx^*}_{M,\bfw'} - \norm*{\bfA\bfx^* - \bfb}_{M,\bfw'} && \text{triangle inequality} \\
&> 21\kappa\cdot\OPT - 10\kappa\cdot\OPT  = 11\kappa\cdot\OPT.
\end{align*}

This contradicts Equation (<ref>), so we conclude that

\[
\norm*{\bfA\tilde\bfx - \bfb}_{M,\bfw} \leq \frac{25\kappa}{c_L^{1/p_M}}\cdot\OPT = O(\kappa)\OPT.\qedhere
\]

§.§ Relative Error Approximation

After an initial constant factor approximation, we show that sensitivity sampling on the residual yields a relative error approximation.

We first reduce to considering only $\bfA\bfx$ with $\norm*{\bfA\bfx}_{M,\bfw} = O(\OPT)$.

Let $\bfA\in\mathbb R^{n\times d}$ and let $\bfb\in\R^n$. Let $\bfw\geq\mathbf{1}_n$ and let $\norm*{\cdot}_M$ be an $M$-norm that satisfies an approximate triangle inequality, that is,

\[
M^{1/p_M}(a+b) \leq O(1)(M^{1/p_M}(a) + M^{1/p_M}(b)).
\]

Suppose that $\bfx_c\in\R^d$ satisfies

\[
\norm*{\bfA\bfx_c - \bfb}_{M,\bfw} \leq O(1)\min_\bfx\norm*{\bfA\bfx - \bfb}_{M,\bfw} = O(\OPT)
\]

and let $\bfb' = \bfb - \bfA\bfx_c$. Let $\bar\bfx$ satisfy

\[
\norm*{\bfA\bar\bfx-\bfb'}_{M,\bfw} \leq O(1)\min_\bfx\norm*{\bfA\bfx-\bfb'}_{M,\bfw} = O(\OPT).
\]

Then, $\norm*{\bfA\bar\bfx}_{M,\bfw} = O(\OPT)$.

We have by approximate triangle inequality that

\[
\norm*{\bfA\bar\bfx}_{M,\bfw} \leq O(1)(\norm*{\bfA\bar\bfx-\bfb'}_{M,\bfw} + \norm*{\bfb'}_{M,\bfw}) \leq O(1)(O(\OPT) + O(\OPT)) = O(\OPT).
\]

Next, we show that we only need to preserve the cost on a certain subset of coordinates where the entries of $\bfb'$ are relatively small.

Consider the setting of Lemma <ref>. Further suppose that either

* for any $a,b\in\R$ with $M(a)\leq\eps^{p_M}M(b)$,

\[
M(a+b) = (1\pm O(\eps))M(b)
\]

* $M^{1/p_M}$ is subadditive

Let $\tilde\bfs_i^{M,\bfw}(\bfA) \geq \bfs_i^{M,\bfw}(\bfA)$ be upper bounds on the weighted $M$-sensitivities of $\bfA$. Let

\[
\mathcal B \coloneqq \braces*{i\in[n] : \bfw_i M(\bfb_i') > \frac{\tilde\bfs_i^{M,\bfw}(\bfA)\OPT^{p_M}}{\eps^{p_M}}}
\]

Let $\bfw'$ be a set of random weights such that for any fixed $\bfy\in\R^n$, $\E_{\bfw'}\norm*{\bfy}_{M,\bfw'}^{p_M} = \norm*{\bfy}_{M,\bfw}^{p_M}$, and such that with probability at least $1-\delta$ satisfies

\begin{equation}\label{eq:good-coordinates-guarantee}
\norm*{(\bfA\bfx-\bfb')\mid_{\overline{\mathcal B}}}_{M,\bfw'}^{p_M} = \norm*{(\bfA\bfx-\bfb')\mid_{\overline{\mathcal B}}}_{M,\bfw}^{p_M} \pm O(\eps)\OPT^{p_M}
\end{equation}

for every $\bfA\bfx$ with $\norm*{\bfA\bfx}_{M,\bfw} = O(\OPT)$. Then with probability at least $1-\delta$, there is a $C$ with $\abs*{C} \leq O(\OPT^{p_M}/\delta)$ such that

\begin{align*}
\abs*{\norm*{\bfy - \bfb'}_{M,\bfw'}^{p_M} - \norm*{\bfy - \bfb'}_{M,\bfw'}^{p_M} - C} = O(\eps/\delta)\OPT^{p_M}.
\end{align*}

First note that assumption (2) implies assumption (1), since that implies that

\[
M^{1/p_M}(a+b) = M^{1/p_M}(b)\pm M^{1/p_M}(a) = M^{1/p_M}(b) \pm \eps M^{1/p_M}(b)
\]

and raising both sides to the $p_M$th power gives assumption (1).

Now note that for any $i\in \mathcal B$ and $\bfy = \bfA\bfx$ with $\norm*{\bfA\bfx}_{M,\bfw} = O(\OPT)$, we have that

\[
\bfw_i M([\bfA\bfx](i)) \leq \tilde\bfs_i^{M,\bfw}(\bfA)\cdot\norm*{\bfA\bfx}_{M,\bfw}^{p_M} = O(1)\tilde\bfs_i^{M,\bfw}(\bfA)\cdot\OPT^{p_M} \leq O(\eps^{p_M})\bfw_i M(\bfb'(i)).
\]

It follows by assumption (1) that

\[
M(\abs*{\bfy_i - \bfb'_i}) = (1\pm O(\eps))M(\abs*{\bfb'_i})
\]

Then, we have that

\begin{align*}
\norm*{\bfy - \bfb'}_{M,\bfw}^{p_M} &= \norm*{(\bfy - \bfb')\mid_{\overline{\mathcal B}}}_{M,\bfw}^{p_M} + \norm*{(\bfy - \bfb')\mid_{\mathcal B}}_{M,\bfw}^{p_M} \\
&= \norm*{(\bfy - \bfb')\mid_{\overline{\mathcal B}}}_{M,\bfw}^{p_M} + \norm*{\bfb'\mid_{\mathcal B}}_{M,\bfw}^{p_M} \pm O(\eps)\norm*{\bfb'\mid_{\mathcal B}}_{M,\bfw}^{p_M}
\end{align*}

and similarly,

\begin{align*}
\norm*{\bfy - \bfb'}_{M,\bfw'}^{p_M} &= \norm*{(\bfy - \bfb')\mid_{\overline{\mathcal B}}}_{M,\bfw'}^{p_M} + \norm*{(\bfy - \bfb')\mid_{\mathcal B}}_{M,\bfw'}^{p_M} \\
&= \norm*{(\bfy - \bfb')\mid_{\overline{\mathcal B}}}_{M,\bfw'}^{p_M} + \norm*{\bfb'\mid_{\mathcal B}}_{M,\bfw'}^{p_M} \pm O(\eps)\norm*{\bfb'\mid_{\mathcal B}}_{M,\bfw'}^{p_M}
\end{align*}

Now using that $\norm*{\bfb'}_{M,\bfw'}^{p_M} = O(1/\delta)\norm*{\bfb'}_{M,\bfw}^{p_M}$ with probability at least $1-\delta$ by Markov's inequality, we set $C = \norm*{\bfb'\mid_{\mathcal B}}_{M,\bfw}^{p_M} - \norm*{\bfb'\mid_{\mathcal B}}_{M,\bfw'}^{p_M}$ to see that

\begin{align*}
\abs*{\norm*{\bfy - \bfb'}_{M,\bfw'}^{p_M} - \norm*{\bfy - \bfb'}_{M,\bfw'}^{p_M} - C} &= \abs*{\norm*{(\bfy - \bfb')\mid_{\overline{\mathcal B}}}_{M,\bfw}^{p_M} - \norm*{(\bfy - \bfb')\mid_{\overline{\mathcal B}}}_{M,\bfw'}^{p_M}} \pm O(\eps/\delta)\norm*{\bfb'\mid_{\mathcal B}}_{M,\bfw}^{p_M} \\
&\leq O(\eps/\delta)\norm*{\bfb'\mid_{\mathcal B}}_{M,\bfw}^{p_M} = O(\eps/\delta)\OPT^{p_M}
\end{align*}

as desired.

We next show how to guarantee Equation (<ref>) using Bernstein's inequality. The proof closely follows that of Lemma <ref>, but needs a slight modification to handle the coordinates of $\bfb'$ in $\overline{\mathcal B}$.

Consider the setting of Lemma <ref>. Let $\bfA\bfx$ have $\norm*{\bfA\bfx}_{M,\bfw} = O(\OPT)$ and let $\bfy = \bfA\bfx- \bfb'$. Let $m\geq 1$ be a parameter and let $\bfw'$ be obtained from weights $\bfw$ according to Definition <ref>. Then,

\[
\Pr\braces*{\abs*{\norm*{(\bfA\bfx-\bfb')\mid_{\overline{\mathcal B}}}_{M,\bfw'}^{p_M} - \norm*{(\bfA\bfx-\bfb')\mid_{\overline{\mathcal B}}}_{M,\bfw}^{p_M}} \geq \eps\OPT^{p_M}} \leq \exp\parens*{-\Theta(1)m\eps^{2+p_M}}
\]

Let $G\subseteq \overline{\mathcal B}$ be the subset of coordinates of $\overline{\mathcal B}$ such that $\bfp_i < 1$. Define the random variable

\[
W_i \coloneqq \bfw_i' M(\bfy(i))
\]

for each $i\in G$. Then,

\[
\E\bracks*{\sum_{i\in G} W_i} = \sum_{i\in G} \frac{\bfw_i}{\bfp_i}M(\bfy(i))\cdot \bfp_i = \norm*{\bfy\mid_G}_{M,\bfw}^{p_M}.
\]

Note that

\begin{align*}
\frac{\bfw_i}{\bfp_i}M(\bfy(i)) &\leq O(1)\frac{1}{m\cdot \tilde\bfs_i^{M,\bfw}(\bfA)} (\bfw_i M([\bfA\bfx](i)) + \bfw_i M(\bfb'(i))) \\
&\leq O(1)\frac{1}{m\cdot \tilde\bfs_i^{M,\bfw}(\bfA)} \parens*{\tilde\bfs_i^{M,\bfw}(\bfA)\norm*{\bfA\bfx}_{M,\bfw}^{p_M} + \frac{\tilde\bfs_i^{M,\bfw}\OPT^{p_M}}{\eps^{p_M}}} \\
&\leq O(1)\frac{1}{m\cdot \tilde\bfs_i^{M,\bfw}(\bfA)} \parens*{\tilde\bfs_i^{M,\bfw}(\bfA)\OPT^{p_M} + \frac{\tilde\bfs_i^{M,\bfw}\OPT^{p_M}}{\eps^{p_M}}} \\
&\leq O(1)\frac{\OPT^{p_M}}{m\cdot\eps^{p_M}} \\
\end{align*}

We next bound the variance:

\begin{align*}
\Var\bracks*{\sum_{i\in G} W_i} &= \sum_{i\in G} \Var[W_i] \leq \sum_{i\in G} \frac{\bfw_i^2}{\bfp_i^2}M(\bfy(i))^2\cdot \bfp_i = \sum_{i\in G} \frac{\bfw_i}{\bfp_i}M(\bfy(i))\cdot \bfw_i M(\bfy(i)) \\
&\leq O(1)\frac{\OPT^{p_M}}{m\cdot\eps^{p_M}} \sum_{i\in G} \bfw_i M(\bfy(i)) = O(1)\frac{\OPT^{2p_M}}{m\cdot\eps^{p_M}}
\end{align*}

Then by Bernstein's inequality,

\begin{align*}
\Pr\braces*{\abs*{\sum_{i\in G} W_i - \norm*{\bfy\mid_G}_{M,\bfw}^{p_M}} > t} &\leq 2\exp\parens*{-\Theta(1)\frac{t^2}{\frac{1}{m\eps^{p_M}}\OPT^{2p_M} + \frac{1}{m\eps^{p_M}} \OPT^{p_M} t}} \\
&= 2\exp\parens*{-\Theta(1)\frac{mt^2 \eps^{p_M}}{\OPT^{p_M}(\OPT^{p_M} + t)}}.
\end{align*}

For $t = \eps\OPT^{p_M}$, this gives a bound of

\[
2\exp\parens*{-\Theta(1)\frac{m(\eps\OPT^{p_M})^2\eps^{p_M}}{\OPT^{2p_M}}} =  2\exp\parens*{-\Theta(1)m\eps^{2+p_M}}.
\]

Consider the setting of Lemma <ref>. Recall that $\bfx_c$ satisfies $\norm*{\bfA\bfx_c-\bfb}_{M,\bfw} = O(\OPT)$ and $\bfb' = \bfb-\bfA\bfx_c$. Let

\[
m = O\parens*{\frac{d}{\eps^{2+p_M}}\parens*{\log\frac1\eps}\parens*{\log\frac1\delta}}
\]

and let $\bfw'$ be as defined in Lemma <ref> with the above choice of $m$. Let $\bar\bfx$ satisfy $\norm*{\bfA\bar\bfx - \bfb'}_{M,\bfw'}^{p_M} \leq (1+\eps)\norm*{\bfA\bfx - \bfb'}_{M,\bfw'}^{p_M}$. Then, with probability at least $1-\delta$,

\[
\norm*{\bfA(\bfx_c + \bar\bfx) - \bfb}_{M,\bfw} \leq (1+O(\eps/\delta))\min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx - \bfb}_{M,\bfw}.
\]

By Lemma <ref>, we have that $\norm*{\bfA\tilde\bfx}_{M,\bfw} \leq O(\OPT)$. We may thus restrict our attention to the ball $\mathcal B_{\rho}^{M,\bfw} \subseteq\mathcal V$ of radius $\rho = O(\OPT)$, since

\[
\min_{\bfx\in\mathbb R^d} \norm*{\bfA\bfx - \bfb'}_{M,\bfw'} = \min_{\bfy\in\mathcal B_{\rho}^{M,\bfw}} \norm*{\bfy - \bfb'}_{M,\bfw'}
\]

by the above.

Next, we show Equation (<ref>) for every $\bfA\bfx\in \mathcal B_{\rho}^{M,\bfw}$. By Lemma <ref>, it suffices to show the approximation guarantee for every $\bfy\in\mathcal N$ and $\mathcal N-\bfb'$, for an $O(\eps^{p_M/q_m})\rho$-cover $\mathcal N$ of $\mathcal B_\rho^{M,\bfw}$. Let $\mathcal N$ be such a net over $\mathcal B_\rho^{M,\bfw}$, which has size at most $\abs*{\mathcal N} \leq O(d\log\frac1\eps)$ by Lemma <ref>. We may then use Lemma <ref> and a union bound over the net $\mathcal N$ to get Equation (<ref>) for every $\bfA\bfx\in \mathcal B_{\rho}^{M,\bfw}$.

Finally, note that we now have the conclusion of Lemma <ref>, since we have satisfied its assumption. Then,

\begin{align*}
\norm{\bfA \bar \bfx - \bfb'}_{M,\bfw}^{p_M} &\le \norm{\bfA \bar \bfx - \bfb'}_{M,\bfw'}^{p_M} - C + O(\epsilon/\delta) \cdot \OPT^{p_M}\nonumber\\
&\le (1+\epsilon) \cdot \min_\bfx\norm{\bfA \bfx - \bfb'}_{M,\bfw'}^{p_M} - C + O(\epsilon/\delta) \cdot \OPT^p\nonumber\\
&\le (1+\epsilon) \min_\bfx(\norm{\bfA \bfx - \bfb'}_{M,\bfw}^{p_M} + C) - C + O(\epsilon/\delta) \cdot \OPT^p \nonumber \\
&\le \min_\bfx\norm{\bfA \bfx - \bfb'}_{M,\bfw}^{p_M} + O(\epsilon/\delta) \cdot \OPT^p
\end{align*}

§.§ Nearly Input Sparsity Time Algorithm

We now combine the active regression algorithms obtained in this section with the sensitivity bounds of Section <ref>.

Input sparsity time, constant factor $M$-estimator active regression
input: Matrix $\bfA \in \R^{n \times d}$, measurement vector $\bfb\in\R^n$.

output: Constant factor approximate solution $\bfx_c\in\R^d$ to $\min_\bfx\norm*{\bfA\bfx-\bfb}_M$.
[1]
Run Algorithm <ref> to obtain approximate sensitvities $\tilde\bfs_i^M(\bfA)$ with $\tau = T$[See Theorem <ref> for the definition of this setting.]
Given weights $\mathbf{1}_n$, obtain $\bfw$ according to Definition <ref>
Run Algorithm <ref> to obtain approximate sensitvities $\tilde\bfs_i^{M,\bfw}(\bfA)$ with $\tau = d$
Given weights $\bfw$, obtain $\bfw'$ according to Definition <ref>
Let $\bfx_c$ be an approximate solution $\norm*{\bfA\bfx_c-\bfb}_{M,\bfw'} \leq O(1)\min_\bfx\norm*{\bfA\bfx-\bfb}_{M,\bfw'}$

Input sparsity time, relative error $M$-estimator active regression
input: Matrix $\bfA \in \R^{n \times d}$, measurement vector $\bfb\in\R^n$.

output: Approximate solution $\tilde\bfx\in\R^d$ to $\min_\bfx\norm*{\bfA\bfx-\bfb}_M$.
[1]
Run Algorithm <ref> to obtain a constant factor solution $\bfx_c$
Set $\bfb' \gets \bfb - \bfA\bfx_c$
Run Algorithm <ref> to obtain approximate sensitvities $\tilde\bfs_i^M(\bfA)$ with $\tau = T$
Given weights $\mathbf{1}_n$, obtain $\bfv$ according to Definition <ref>
Run Algorithm <ref> to obtain approximate sensitvities $\tilde\bfs_i^{M,\bfv}(\bfA)$ with $\tau = d$
Given weights $\bfv$, obtain $\bfv'$ according to Definition <ref>
Let $\bar\bfx$ be an approximate solution $\norm*{\bfA\bar\bfx-\bfb}_{M,\bfv'} \leq (1+\eps)\min_\bfx\norm*{\bfA\bfx-\bfb}_{M,\bfv'}$
Let $\tilde\bfx \gets \bfx_c + \bar\bfx$

There is an algorithm, Algorithm <ref>, which with probability at least $99/100$ computes $\tilde\bfx$ such that

\[
\norm*{\bfA\tilde\bfx - \bfb}_M \leq (1+\eps)\min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx - \bfb}_M
\]

in time

\[
O\parens*{\nnz(\bfA)\log^3 n + \frac{(Td^{\max\{1,p_M/2\}}+T^2)\log^4 n}{\eps^{2+p_M}}\log\frac1\eps}
\]

where constant factor $\ell_{p_M}$ Lewis weight approximation for an $n\times d$ matrix $\bfB$ takes time $O(\nnz(\bfB) + T))$ time. Furthermore, the algorithm only queries

\[
O\parens*{\frac{d^{\max\{2,p_M/2+1\}}\log^3 n}{\eps^{2+p_M}}\log\frac1\eps}
\]

entries of $\bfb$.

In order to obtain nearly input sparsity time algorithms, we will need to apply the sensitivity sampling algorithms of Theorem <ref> twice, once to reduce the number of rows to $O(\poly(d,\log n))$ rows in nearly input sparsity time, and once to reduce the number of rows all the way down to $O(d^{\max\{2,p/2+1\}}\poly\log(n))$ rows in an additional $O(\poly(d,\log n))$ time.

Constant Factor Approximation. We will first obtain a constant factor solution $\bfx_c$ such that

\[
\norm*{\bfA\bfx_c - \bfb}_M \leq O(1)\min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx - \bfb}_M.
\]

To do this, we first obtain sensitivity upper bounds $\tilde\bfs_i^M(\bfA)$ by Theorem <ref> with $\tau = T$, so that they sum to

\[
\tilde{\mathcal T}^M(\bfA) = \sum_{i=1}^n \tilde\bfs_i^M(\bfA) = O\parens*{d^{\max\{1,p_M/2\}}\log^2 n + T}
\]

and can be computed in time

\[
O\parens*{\nnz(\bfA)\log^2 n + n\log^2 n} = O\parens*{\nnz(\bfA)\log^2 n}.
\]

By Lemma <ref>, we can perform sensitivity sampling with these sensitivity estimates to obtain a set of weights $\bfw$ such that if $\bfx_c$ satisfies

\[
\norm*{\bfA\bfx_c - \bfb}_{M,\bfw} \leq O(1)\min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx - \bfb}_{M,\bfw},
\]

\[
\norm*{\bfA\bfx_c - \bfb}_{M} \leq O(1)\min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx - \bfb}_{M}.
\]

Next, we obtain weighted sensivitity upper bounds $\tilde\bfs_i^{M,\bfw}(\bfA)$ by Corollary <ref> with $\tau = d$. Note that $\norm*{\bfw}_\infty = O(n)$ since the sensitivity upper bounds are all at least $1/n$, so the weighted sensitivities sum to at most

\[
\sum_{i=1}^n \tilde\bfs_i^{M,\bfw}(\bfA) = O(d^{\max\{1,p_M/2\}}\log^3 n).
\]

Furthermore, the weights $\bfw$ only have

\[
\nnz(\bfw) = O\parens*{d\cdot \tilde{\mathcal T}^M(\bfA)} = O(d^{\max\{2,p_M/2+1\}}\log^2 n + dT)
\]

nonzero entries with constant probability by Markov's inequality, so the running time is

\begin{align*}
&O\parens*{\nnz(\bfA)\log^3(\nnz(\bfw)) + \log n\frac{\nnz(\bfw)T}{d}\log (\nnz(\bfw))} \\
\leq~&O\parens*{\nnz(\bfA)\log^3 n + (Td^{\max\{1,p_M/2\}}+T^2)\log^4 n}.
\end{align*}

By Lemma <ref>, we can perform sensitivity sampling with these sensitivity estimates to obtain a set of weights $\bfw'$. We can now define our $\bfx_c$ to be the minimizer of

\[
\min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx - \bfb}_{M,\bfw'},
\]

which means that

\[
\norm*{\bfA\bfx_c - \bfb}_{M,\bfw} \leq O(1)\min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx - \bfb}_{M,\bfw},
\]

which in turn implies that

\[
\norm*{\bfA\bfx_c - \bfb}_{M} \leq O(1)\min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx - \bfb}_{M},
\]

as desired. Computing $\bfx_c$ only required $\nnz(\bfw')$ entries of $\bfb$ to sample, which is only $O(d^{\max\{2,p_M/2+1\}}\log^3 n)$ entries.

Relative Error Approximation. With the constant factor solution $\bfx_c$ in hand, we can now apply Lemma <ref> to refine this to a relative error approximation. Following Lemma <ref>, let

\[
\bfb' = \bfb - \bfA\bfx_c
\]

Recall the sensitivity upper bounds $\tilde\bfs_i^{M}(\bfA)$ constructed earlier. This time, we obtain a different set of weights $\bfv\geq\mathbf{1}_n$ by following Lemma <ref>, so that the expected number of rows sampled is now

\begin{align*}
\E[\nnz(\bfv)] \leq O\parens*{\frac{d}{\eps^{2+p_M}}\parens*{\log\frac1\eps}\cdot \tilde{\mathcal T}^M(\bfA)} &= O\parens*{\frac{d(d^{\max\{1,p_M/2\}}\log^2 n + T)}{\eps^{2+p_M}}\parens*{\log\frac1\eps}} \\
&= O\parens*{\frac{d^{\max\{2,p_M/2+1\}}\log^2 n + dT}{\eps^{2+p_M}}\parens*{\log\frac1\eps}}.
\end{align*}

By a Markov bound, with constant probability, the number of rows sampled has the same bound, up to constant factors. We are then guaranteed by Lemma <ref> that if $\bar\bfx$ satisfies

\[
\norm*{\bfA\bar\bfx - \bfb'}_{M,\bfv} \leq (1+\eps)\min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx - \bfb'}_{M,\bfv},
\]

\[
\norm*{\bfA\bar\bfx - \bfb'}_{M} \leq (1+O(\eps))\min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx - \bfb'}_{M}.
\]

Then, as done in the previous constant factor approximation step, we perform another round of weighted sensitivity sampling. Let $\tilde\bfs_i^{M,\bfv}(\bfA)$ be weighted sensitivity upper bounds given by Corollary <ref> with $\tau = d$, which sum to at most $O(d\log^3 n)$ and be computed in time at most

\begin{align*}
&O\parens*{\nnz(\bfA)\log^3(\nnz(\bfv)) + \log n\frac{\nnz(\bfv)T}{d}\log (\nnz(\bfv))} \\
\leq~&O\parens*{\nnz(\bfA)\log^3 n + \frac{(Td^{\max\{1,p_M/2\}}+T^2)\log^4 n}{\eps^{2+p_M}}\log\frac1\eps}.
\end{align*}

We can then again do weighted sensitivity sampling according to the $\tilde\bfs_i^{M,\bfv}(\bfA)$ as in Lemma <ref> to obtain a set of weights $\bfv'$ such that if $\bar\bfx$ satisfies

\[
\norm*{\bfA\bar\bfx - \bfb'}_{M,\bfv'} \leq (1+\eps)\min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx - \bfb'}_{M,\bfv'},
\]

\[
\norm*{\bfA\bar\bfx - \bfb'}_{M,\bfv} \leq (1+O(\eps))\min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx - \bfb'}_{M,\bfv}.
\]

It then follows from the previous discussion that

\[
\norm*{\bfA\bar\bfx - \bfb'}_{M} \leq (1+O(\eps))\min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx - \bfb'}_{M}.
\]

Furthermore, computing $\bar\bfx$ only required $\nnz(\bfv')$ entries of $\bfb'$ to sample, which is only

\[
O\parens*{\frac{d\cdot\tilde{\mathcal T}^{M,\bfv}(\bfA)}{\eps^{2+p_M}}\log\frac1\eps} = O\parens*{\frac{d^{\max\{2,p_M/2+1\}}\log^3 n}{\eps^{2+p_M}}\log\frac1\eps}
\]

§ APPLICATIONS: ACTIVE REGRESSION FOR THE HUBER LOSS

Note that our sampling bound for general $M$-estimators from Section <ref> is loose by a factor of $d$ compared to the corresponding bounds for $\ell_p$. This is due to the reliance on a naïve net argument rather than using the more efficient construction of [6, 63, 46]. A natural question then is whether this can be improved or not, for loss functions other than the $\ell_p$ loss. We answer this in the affirmative for the important special case of the Huber loss.

More specificially, in this section, we obtain a sampling bound of

\[
O(d^{4-2\sqrt 2}\poly(\eps^{-1},\log n))
\]

for subspace embeddings as well as active regression for the Huber loss, where $4 - 2\sqrt 2\approx 1.17157$. This bound is substantially better compared to the $O(d^{2}\poly(\eps^{-1},\log n))$ bound for other losses, which is obtained by a naïve Bernstein bound over a net, as well as a previous subspace embedding result of [22] which achieved roughly $O(d^4)$ rows.

In what follows, we will first develop our recursive sampling algorithm for obtaining subspace embeddings for the Huber loss with the above sampling complexity, which can be used to obtain constant factor solutions in the active setting by Lemma <ref>. In Section <ref>, we then show how to modify our subspace embedding construction in order to get relative error bounds, again following Section <ref>.

§.§ Properties of the Huber Loss

We first collect properties of the Huber loss needed to apply previous results for $M$-estimators. It is known that the Huber loss is polynomially bounded above with degree $2$ and constant $1$, and polynomially bounded below with degree $1$ and constant $1$.

The Huber loss satisfies

\[
\frac{y}{x} \leq \frac{H(y)}{H(x)} \leq \frac{y^2}{x^2}
\]

for all $y > x\geq 0$.

This in turn implies that $H^{1/2}$ is subadditive, since it is known that polynomial growth with degree $1$ and constant $1$ implies subadditivity:

Let $f:\mathbb R\to\mathbb R$. Suppose that $x\mapsto x^{-1}f(x)$ is decreasing, that is,

\[
\frac{f(y)}{f(x)} \leq \frac{y}{x}
\]

for all $y > x$. Then, $f$ is subadditive.

The following is a corollary that discretizes Lemma <ref>:

Let $\bfy\in\mathbb R^n$ and let $\alpha \in [0,1/2]$ and $\gamma = n^{-\alpha}$ with $2/n \leq \gamma \leq 1$. Let

\[
T \supseteq \braces*{i\in[n] : H(\bfy_i) \leq \gamma \norm*{\bfy}_H^2},
\]

$\eps = c/\log n$ for some sufficiently small constant $c$, and let

\[
\mathcal I = \braces*{1, 1+\eps, 1+2\eps, \dots, 2 - \eps, 2} = \braces*{1 + \eps\cdot i : i\in \{0\}\cup[1/\eps]}.
\]

Then, for some constant $c>0$, either

\[
\norm*{\bfy\mid_T}_H^2 \geq c\frac1{(\gamma n)^{3-2\sqrt 2}}\min_{p\in \mathcal I}\norm*{\bfy\mid_T}_p^p
\]

\[
\norm*{\bfy}_H^2 \geq c\gamma\min_{p\in\{1,2\}}\norm*{\bfy}_p^p
\]

where $3 - 2\sqrt 2 \approx 0.17157$.

It was shown in the proof of Lemma <ref> that if $\norm*{\bfy}_1\geq 2n$, then $\norm*{\bfy}_H^2 = \Theta(\norm*{\bfy}_1)$, and if $\norm*{\bfy}_\infty \leq 1$, then $\norm*{\bfy}_H^2 = \norm*{\bfy}_2^2$. Thus, we may assume that all the entries of $\bfy$ are bounded by $O(n)$, and $\norm*{\bfy}_\infty \geq 1$. Note then that the entries that are at most $1/2n$ can make up at most an $\ell_p$ mass of $1/2$, and thus only affects the claimed inequality by a constant factor for any $p$. For all other entries, the magnitude is bounded below by $1/\poly(n)$ and above by $\poly(n)$, so an additive $O(1/\log n)$ difference in the exponent can only affect the bounds by a constant factor. We thus conclude the desired result.

§.§ Net Arguments for the Huber Loss

When constructing subspace embeddings for norms, it suffices to preserve the lengths of vectors on the unit normed sphere, by the scale invariance/homogeneity of norms. On the other hand, general $M$-estimators, including the Huber norm, do not necessarily satisfy scale invariance, which means that this strategy, as is, does not work. For the Huber norm, we show that for sufficiently small and sufficiently large scales, we have at least an $(1\pm\eps)$-approximate scale invariance, which allows us to obtain subspace embeddings by union bounding over finitely many scales. For all scales near $0$, this is because the Huber norm coincides with the $\ell_2$ norm. For larger scales, we show that the Huber norm is the $\ell_1$ norm up to $(1\pm\eps)$ factors. Then, the following result then shows that in order to prove subspace embeddings, it suffices to preserve the Huber norms of Huber balls at approximately $\log n$ radii.

The following lemma shows that for extremely large radii, the Huber norm essentially coincides with the $\ell_1$ norm.

Let $\bfw\in\mathbb R^n$ be a set of weights. Suppose $\bfy\in\mathbb R^n$ has weighted Huber norm at least $\norm*{\bfy}_{H,\bfw}^2 \geq n\norm*{\bfw}_\infty/\eps$. Then,

\[
\norm*{\bfy}_{H,\bfw}^2 = (1\pm\eps)\norm*{\bfy}_{1,\bfw}.
\]

Suppose that $\bfy\in\mathbb R^n$ has Huber norm at least $\norm*{\bfy}_H^2 \geq n\norm*{\bfw}_\infty/\eps$. Now define the set

\[
S \coloneqq \braces*{i\in[n] : \abs*{\bfy_i} \leq 1}.
\]

Then, $\norm*{\bfy\mid_S}_{1,\bfw} \leq n\norm*{\bfw}_\infty$ so

\[
\norm*{\bfy}_{1,\bfw} = \norm*{\bfy\mid_{\overline S}}_{1,\bfw} + \norm*{\bfy\mid_{S}}_{1,\bfw} \leq \norm*{\bfy\mid_{\overline S}}_{H,\bfw}^2 + n\norm*{\bfw}_\infty \leq (1+\eps)\norm*{\bfy\mid_{\overline S}}_{H,\bfw}^2.
\]

On the other hand, we have $\norm*{\bfy}_{1,\bfw} \geq \norm*{\bfy}_{H,\bfw}^2$, so we conclude.

We now show the nets we need to obtain subspace embeddings for the Huber norm.

Let $\bfA\in\mathbb R^{n\times d}$ and let $\mathcal V = \Span(\bfA)$. Let $\bfw,\bfw'\geq\mathbf{1}_n$ be two sets of weights. Let $\bfw'\geq\mathbf{1}_n$ be another set of weights. Let $\eps>0$. Let

\[
\ell \coloneqq \log_2\sqrt{\frac{n\norm*{\bfw'}_\infty}{\eps}} = O\parens*{\log\frac{n\norm*{\bfw'}_\infty}{\eps}}.
\]

For each $k\in[\ell]$, let $\rho = 2^k$ and suppose that

\[
\norm*{\bfy}_{H,\bfw'} = \norm*{\bfy}_{H,\bfw} \pm \eps \rho
\]

for every $\bfy\in \mathcal B_\rho^H$, for every $k\in[\ell]$. Then,

\[
\norm*{\bfy}_{H,\bfw'} = (1\pm O(\eps))\norm*{\bfy}_{H,\bfw}
\]

for all $\bfy\in\mathcal V$.

We handle three cases: scales near $0$, scales near infinity, and all scales in between.

Scales near $0$. Note that if the weighted Huber norm of $\bfy\in\mathbb R^n$ is at most $1$, then all entries of $\bfy$ are at most $1$ in absolute value, so $\norm*{\bfy}_{H,\bfw}^2 = \norm*{\bfy}_{2,\bfw}^2$, and thus the Huber norm is scale invariant at these scales. Thus, $(1\pm\eps)$-approximations for the unit Huber sphere implies $(1\pm\eps)$-approximations for all vectors with $\norm*{\bfy}_{H,\bfw}^2 \leq 1$.

Scales near infinity.

Suppose that we have the guarantee that

\[
\norm*{\bfx}_{H,\bfw'}^2 = (1\pm\eps)\norm*{\bfx}_{H,\bfw}^2
\]

for every $\bfx$ in the Huber sphere $\mathcal S_{n\norm*{\bfw'}_\infty/\eps}^{H,\bfw}$ of radius $n\norm*{\bfw'}_\infty/\eps$. Then, if $\bfx$ has any Huber norm at least $n\norm*{\bfw'}_\infty/\eps$, then

\begin{align*}
\norm*{\bfx}_{H,\bfw}^2 &= (1\pm\eps)\norm*{\bfx}_{1,\bfw} && \text{Lemma \ref{lem:big-huber-is-l1}} \\
&= (1\pm\eps)\frac{\eps}{n\norm*{\bfw'}_\infty}\norm*{\bfx}_{1,\bfw} \cdot \norm*{\frac{n\norm*{\bfw'}_\infty}{\eps}\frac{\bfx}{\norm*{\bfx}_{1,\bfw}}}_{1,\bfw} && \text{Scale invariance} \\
&= (1\pm\eps)^2\frac{\eps}{n\norm*{\bfw'}_\infty}\norm*{\bfx}_{1,\bfw} \cdot \norm*{\frac{n\norm*{\bfw'}_\infty}{\eps}\frac{\bfx}{\norm*{\bfx}_{1,\bfw}}}_{H,\bfw}^2 && \text{Lemma \ref{lem:big-huber-is-l1}} \\
&= (1\pm\eps)^3\frac{\eps}{n\norm*{\bfw'}_\infty}\norm*{\bfx}_{1,\bfw} \cdot \norm*{\frac{n\norm*{\bfw'}_\infty}{\eps}\frac{\bfx}{\norm*{\bfx}_{1,\bfw}}}_{H,\bfw'}^2 && \text{Approximation guarantee} \\
&= (1\pm\eps)^4\frac{\eps}{n\norm*{\bfw'}_\infty}\norm*{\bfx}_{1,\bfw} \cdot \norm*{\frac{n\norm*{\bfw'}_\infty}{\eps}\frac{\bfx}{\norm*{\bfx}_{1,\bfw}}}_{1,\bfw'} && \text{Lemma \ref{lem:big-huber-is-l1}} \\
&= (1\pm\eps)^4 \cdot \norm*{\bfx}_{1,\bfw'}^2 && \text{Scale invariance} \\
&= (1\pm\eps)^5 \cdot \norm*{\bfx}_{H,\bfw'}^2 && \text{Lemma \ref{lem:big-huber-is-l1}}
\end{align*}

Thus, it suffices to prove approximation guarantees for the Huber sphere $\mathcal S_{n\norm*{\bfw'}_\infty/\eps}^{H,\bfw}$ to handle all $\bfx\in\mathbb R^n$ with Huber norms at least $n\norm*{\bfw'}_\infty/\eps$.

Scales in between.

For the remaining scales between $1$ and $n\norm*{\bfw'}_\infty/\eps$, we work in powers of $2$ so that there are approximately

\[
\ell \coloneqq \log_{2}\sqrt{\frac{n\norm*{\bfw}_\infty}{\eps}} = O\parens*{\log\frac{n\norm*{\bfw}_\infty}{\eps}}
\]

total scales. Then for any $\rho \in [1, \sqrt{n\norm*{\bfw}_\infty/\eps}]$, there exists a $\rho \leq \rho_k \leq 2\rho$ such that we can get an additive error of $O(\eps)\rho_k = O(\eps)\rho$ for every $\bfy\in \mathcal B_{\rho_k}$, which in particular contains $\mathcal B_\rho^H$.

§.§ Compact Rounding for the Huber Loss

In Section <ref>, we have reduced our task to proving approximation guarantees for a small number of Huber spheres. We now focus on showing small sampling bounds for a single Huber sphere.

The following rounding lemma is implicit in the proof of Theorems 7.3 and 7.4 of [6], which uses these net constructions.

Let $1\leq p < \infty$ and let

\[
\mathcal B_p \coloneqq \braces*{\bfA\bfx : \norm*{\bfA\bfx}_p \leq 1}.
\]

Let $\eps \geq \gamma > 0$ and let $\mathcal N_\gamma$ be a standard $\gamma$-net on $\mathcal B_p$ with the $\ell_p$ norm, which has

\[
\log \abs*{\mathcal N_\gamma} \leq O\parens*{d\log\frac{1}{\gamma}}
\]

by Lemma <ref>. Let $\bfy\in \mathcal N_\gamma$. Let $\ell = \ceil*{\log_{1+\eps}(2^{1/p}d^{1/p\lor 1/2})}$. Then, for $k\in[\ell]$, there exists a rounding

\[
\tilde\bfy = \sum_{k=0}^\ell \tilde\bfy^{(k)}
\]

such that:

* $\abs*{\tilde\bfy(i) - \bfy(i)} \leq 3\eps \abs*{\bfy(i)}$ for all $i\in[n]$

* $\abs*{\tilde\bfy^{(k)}(i)} \leq \bfW_{i,i}^{1/p}(1+\eps)^k$

* $\tilde\bfy^{(k)}$ for $0\leq k\leq \ell$ have disjoint supports

* each $\tilde\bfy^{(k)}$ is drawn from a set of vectors $\mathcal D_k$ with size at most

\[
\log\abs*{\mathcal D_k} \leq c(p)\frac{d}{\eps^{1+\beta}(1+\eps)^{\beta k}}\parens*{\log\frac{n}{\eps} + \log\frac1\gamma}
\]

where $\beta = (2\land p)$ and $c(p) = O\parens*{\max\braces*{p, \frac1{p-1}}}$.

The lemma above is used to prove the following similar result for the Huber loss.

Let $\eps\in(0,1)$ and $\rho\geq 1$. Let $\bfA\in\mathbb R^{n\times d}$ and let $\bfy\in \mathcal B_\rho^H \subseteq \Span(\bfA)$. Let $\mathcal I\subset[1,2]$ be the set in the statement of Corollary <ref> and suppose that $\bfy$ satisfies

\[
\frac{\min_{p\in\mathcal I}\norm*{\bfy}_p^p}{\norm*{\bfy}_H^2} \leq \kappa.
\]

Then, there exists a $p\in\mathcal I$, $\ell = O((\log d)/\eps^2)$, and a rounding

\[
\tilde\bfy = \sum_{k=0}^\ell \tilde\bfy^{(k)}
\]

such that:

* $\norm*{\tilde\bfy - \bfy}_H \leq O(\eps)\norm*{\bfy}_H$

* $H\parens*{\abs*{\tilde\bfy^{(k)}(i)}} \leq O(\kappa)\parens*{\frac{\bfw_i^p(\bfA)}{d} + \frac1n} (1+\eps^2)^{pk}\norm*{\bfy}_H^2$

* $\tilde\bfy^{(k)}$ for $0\leq k\leq \ell$ have disjoint supports

* each $\tilde\bfy^{(k)}$ is drawn from a set of vectors $\mathcal D_k$ with size at most

\[
\log\abs*{\mathcal D_k} \leq O(1)\frac{d}{\eps^{2(1+p)}(1+\eps^2)^{pk}}(\log n)\parens*{\log\frac{n}{\eps}}.
\]

Let $\gamma = \poly(\eps/n)$ and let $\mathcal N_\gamma^p$ be a $\gamma$-net for the $\ell_p$ unit ball for $p\in\mathcal I$. We apply the results of Lemma <ref> to each of these nets to obtain sets of vectors $\mathcal D_k^p$, corresponding to the set of vectors $\mathcal D_k$ in the lemma statement when we apply the result with $p$.

Now let $\bfy\in\mathcal B_\rho^H$. By assumption, we may find a $p\in\mathcal I$ so that

\[
\rho^2 \leq \norm*{\bfy}_p^p \leq \kappa\rho^2.
\]

We then let $L \leq O(\log \kappa)\leq O(\log n)$ be such that

\[
\rho^2 \cdot 2^L \leq \norm*{\bfy}_p^p \leq \rho^2 \cdot 2^{L+1}.
\]

\[
\beta\coloneqq(\rho^2 \cdot 2^{L+1})^{1/p} = \Theta(\norm*{\bfy}_p).
\]

Then, $\norm*{\bfy/\beta}_p \leq 1$ so there exists a $\bfy'\in\mathcal N_\gamma^p$ such that

\[
\norm*{\bfy' - \bfy/\beta}_p \leq \gamma
\]

\[
\norm*{\beta\bfy' - \bfy}_H^2 \leq \norm*{\beta\bfy' - \bfy}_p^p \leq \gamma^p \beta^p \leq 2\gamma^p \norm*{\bfy}_p^p \leq 2\gamma^p \kappa \norm*{\bfy}_H^2 \leq \eps^2\norm*{\bfy}_H^2
\]

for an appropriate choice of $\gamma$.

We now let

\[
\tilde\bfy = \sum_{k=0}^\ell \tilde\bfy^{(k)}
\]

be the rounding given in Lemma <ref> for $\bfy'$, which is in the $\ell_p$ unit ball, with the $\eps$ in the lemma set to our $\eps^2$. Then,

\[
\abs*{\beta\bfy'(i) - \beta\tilde\bfy(i)} \leq 3\eps^2 \abs*{\beta\bfy'(i)}
\]

for all $i\in[n]$ so

\begin{align*}
\norm*{\beta\bfy' - \beta\tilde\bfy}_H^2 &= \sum_{i=1}^n H(\abs*{\beta\bfy'(i) - \beta\tilde\bfy(i)}) \\
&\leq \sum_{i=1}^n H(3\eps^2 \abs*{\beta\bfy'(i)}) \\
&\leq 3\eps^2 \sum_{i=1}^n H(\abs*{\beta\bfy'(i)}) \\
&= 3\eps^2\norm*{\beta\bfy'}_H^2 \\
&= O(\eps^2)\norm*{\bfy}_H^2.
\end{align*}

Thus by the triangle inequality for the Huber norm,

\[
\norm*{\beta\tilde\bfy - \bfy}_H \leq O(\eps)\norm*{\bfy}_H.
\]

\[
\abs*{\beta\tilde\bfy^{(k)}(i)} \leq \bfW_{i,i}^{1/p}(1+\eps^2)^k \beta
\]

\begin{align*}
H\parens*{\abs*{\beta\tilde\bfy^{(k)}(i)}} &\leq O(1)\bfW_{i,i}(1+\eps^2)^{pk}\beta^p && \text{$H(y) \leq O\parens*{\min_{p\in[1,2]} y^p}$} \\
&\leq O(1)\parens*{\frac{\bfw_i^p(\bfA)}{d} + \frac1n} (1+\eps^2)^{pk} \cdot 2\norm*{\bfy}_p^p \\
&\leq O(\kappa)\parens*{\frac{\bfw_i^p(\bfA)}{d} + \frac1n} (1+\eps^2)^{pk}\norm*{\bfy}_H^2.
\end{align*}

We conclude as desired by rescaling the $\tilde\bfy$ and $\tilde\bfy^{(k)}$ by $\beta$.

§.§ Sampling Bounds

With the above net results in hand, we can now complete the argument for one step of the sampling recursion. In order to implement a recursive sampling scheme, note that we need to handle weighted Huber norms. These weighted Huber norms will be handled by handling the weights in a small number of groups such that the weights are within constant factors of each other.

§.§.§ Sampling Guarantees for a Single Weight Class and a Single Radius

We start with a Bernstein bound using the net results, for a single group of weights and at a single radius.

Let $\bfA\in\mathbb R^{n\times d}$ and let $\mathcal V = \Span(\bfA)$. Let $\eps\in(0,1)$ be an accuracy parameter, $\delta>0$ a failure rate parameter, $\rho \geq 1$ a Huber radius, and $\kappa\geq 1$ a distortion parameter. Let

\[
m \coloneqq \frac{d}{\poly(\eps)}(\log^2 d)\parens*{\log^2\frac{n}{\eps}}\parens*{\log\frac{\log d}{\delta\eps}}.
\]

Let $\mathcal I$ be as in the statement of Corollary <ref>. Let $\bfw\geq \mathbf{1}_n$ be a set of weights such that

\begin{equation}\label{eq:bounded-weight-ratio}
\frac{\max_{i\in[n]} \bfw_i}{\min_{i\in[n]} \bfw_i}\leq 2
\end{equation}

and let $w = \min_{i\in[n]}\bfw_i$. Let $\bfw'$ be chosen randomly so that

\[
\bfw_i' \coloneqq \begin{cases}
\bfw_i/\bfp_i & \text{w.p. $\bfp_i$} \\
0 & \text{otherwise}
\end{cases}
\]

\[
\bfp_i \geq \min\braces*{1, m\cdot \kappa\bracks*{\frac1n + \sum_{p\in\mathcal I}\frac{\bfw_i^p(\bfA)}{d}}}
\]

\[
S_\kappa \coloneqq \braces*{\bfy\in\mathcal V : \frac{\min_{p\in\mathcal I}\norm*{\bfy}_p^p}{\norm*{\bfy}_H^2}\leq \kappa}
\]

Then with probability at least $1 - \delta$,

\[
\norm*{\bfy}_{H,\bfw'} = \norm*{\bfy}_{H,\bfw} \pm \eps w\rho
\]

for all $\bfy\in \mathcal N\cap S_\kappa$, where $\mathcal N$ is an $\eps^2\rho$-net over $\mathcal B_\rho^H$ given by Lemma <ref>, by setting $\eps$ in the lemma to $\eps^2$.

Let $\tilde\bfy\in\mathcal N\cap S_\kappa$ be a net vector in the result of Lemma <ref> and let $p$ be the corresponding $p\in\mathcal I$. Note that

\[
\abs*{\norm*{\tilde\bfy}_{H,\bfw'}^2 - \norm*{\tilde\bfy}_{H,\bfw}^2} \leq \sum_{k=0}^\ell \abs*{\norm*{\tilde\bfy^{(k)}}_{H,\bfw'}^2 - \norm*{\tilde\bfy^{(k)}}_{H,\bfw}^2}
\]

by the disjointness of the supports of the $\tilde\bfy^{(k)}$. Thus, it suffices to bound each term in the sum by $\eps\rho^2/(\ell+1)$.

Fix a $0\leq k\leq \ell$ and define the random variable

\[
W_i \coloneqq \bfw_i' H(\tilde\bfy^{(k)}(i))
\]

for each $i\in[n]$. Then,

\[
\E\bracks*{\sum_{i=1}^n W_i} = \sum_{i=1}^n \frac{\bfw_i}{\bfp_i}H(\tilde\bfy^{(k)}(i))\cdot \bfp_i = \norm*{\tilde\bfy^{(k)}}_{H,\bfw}^2.
\]

We next bound the variance:

\begin{align*}
\Var\bracks*{\sum_{i=1}^n W_i} &= \sum_{i=1}^n \Var[W_i] \\
&\leq \sum_{i=1}^n \frac{\bfw_i^2}{\bfp_i^2}H(\tilde\bfy^{(k)}(i))^2\cdot \bfp_i \\
&= \sum_{i=1}^n \frac1{\bfp_i}\bfw_i H(\tilde\bfy^{(k)}(i))\cdot \bfw_i H(\tilde\bfy^{(k)}(i))
\end{align*}

Note that $W_i$ is almost surely bounded by

\begin{align*}
\frac1{\bfp_i}\bfw_i H(\tilde\bfy^{(k)}(i)) &\leq \bfw_i\frac{O(\kappa)\parens*{\frac{\bfw_i^p(\bfA)}{d} + \frac1n}(1+\eps^4)^{pk}\norm*{\bfy}_H^2}{m\cdot O(\kappa)\parens*{\frac{\bfw_i^p(\bfA)}{d} + \frac1n}} && \text{Lemma \ref{lem:blm-net-for-huber}} \\
&\leq \frac2m (1+\eps^4)^{pk}\norm*{\bfy}_{H,\bfw}^2 && \text{Equation \eqref{eq:bounded-weight-ratio}}
\end{align*}

so the variance is bounded by

\[
\Var\bracks*{\sum_{i=1}^n W_i} \leq \frac2m (1+\eps^4)^{pk}\norm*{\bfy}_{H,\bfw}^2\sum_{i=1}^n \bfw_iH(\tilde\bfy^{(k)}(i)) = \frac{2}{m}(1+\eps^4)^{pk}\norm*{\bfy}_{H,\bfw}^4.
\]

Then by Bernstein's inequality,

\begin{align*}
\Pr\braces*{\abs*{\sum_{i=1}^n W_i - \norm*{\tilde\bfy^{(k)}}_{H,\bfw}^2} > t} &\leq 2\exp\parens*{-\Theta(1)\frac{t^2}{\frac{1}{m}(1+\eps^4)^{pk}\norm*{\bfy}_{H,\bfw}^4 + \frac1m (1+\eps^4)^{pk}\norm*{\bfy}_{H,\bfw}^2 t}} \\
&= 2\exp\parens*{-\Theta(1)\frac{mt^2}{(1+\eps^4)^{pk}\norm*{\bfy}_{H,\bfw}^2(\norm*{\bfy}_{H,\bfw}^2 + t)}}.
\end{align*}

For $t = \eps w^2\rho^2/(\ell+1)$, this gives a bound of

\[
2\exp\parens*{-\Theta(1)\frac{m(\eps w^2\rho^2/(\ell+1))^2}{(1+\eps^4)^{pk}w^4\rho^4}} =  2\exp\parens*{-\Theta(1)\frac{m\eps^2}{(\ell+1)^2(1+\eps^4)^{pk}}}.
\]

We then set

\begin{align*}
m &= O(1)\frac{d}{\eps^{4(1+p)+2}}(\ell+1)^2\log\frac{n}{\eps}\log\frac{\ell}\delta \\
&= O(1)\frac{d}{\eps^{4p+10}}(\log^2 d)\parens*{\log^2\frac{n}{\eps}}\parens*{\log\frac{\log d}{\delta\eps}}
\end{align*}

which is enough to union bound over the set $\mathcal D_k$ (see Lemma <ref>) of log size at most

\[
\]

with failure probability at most $\delta/(\ell+1)$. We union bound over the $\ell+1$ choices of $k$ to obtain that

\[
\abs*{\norm*{\tilde\bfy}_{H,\bfw'}^2 - \norm*{\tilde\bfy}_{H,\bfw}^2} \leq \eps w^2\rho^2.
\]

Setting $\eps$ in the above result to $\eps^2$, we get

\[
\abs*{\norm*{\tilde\bfy}_{H,\bfw'}^2 - \norm*{\tilde\bfy}_{H,\bfw}^2} \leq \eps^2 w^2\rho^2.
\]

Now note that if both $\norm*{\tilde\bfy}_{H,\bfw'}$ and $\norm*{\tilde\bfy}_{H,\bfw}$ are at most $\eps w\rho/2$, then we have that

\[
\abs*{\norm*{\tilde\bfy}_{H,\bfw'} - \norm*{\tilde\bfy}_{H,\bfw}} \leq \eps w\rho
\]

by the triangle inequality. Otherwise, we have that

\[
\abs*{\norm*{\tilde\bfy}_{H,\bfw'} - \norm*{\tilde\bfy}_{H,\bfw}} = \frac{\abs*{\norm*{\tilde\bfy}_{H,\bfw'}^2 - \norm*{\tilde\bfy}_{H,\bfw}^2}}{\abs*{\norm*{\tilde\bfy}_{H,\bfw'} + \norm*{\tilde\bfy}_{H,\bfw}}}\leq \frac{\eps^2 w^2\rho^2}{\eps w\rho} = \eps w\rho.\qedhere
\]

We can now apply the result twice, once on all of the rows and once on only the set of rows with small Huber sensitivity, in order to obtain a sampling bound for preserving a Huber sphere at a single scale.

Let $\bfA\in\mathbb R^{n\times d}$ and let $\mathcal V = \Span(\bfA)$. Let $\eps\in(0,1)$ be an accuracy parameter, $\delta>0$ a failure rate parameter, and $\rho \geq 1$ a Huber radius. Let $\beta = 3 - 2\sqrt 2$, $\gamma = n^{-\beta/(\beta+1)}$, and let $T$ be a subset

\[
T \subseteq \braces*{i\in[n] : \bfs_i^H(\bfA) \leq \gamma}
\]

Let $\kappa \coloneqq O((\gamma n)^{\beta}) = O(\gamma^{-1})$. Let $\bfw\geq \mathbf{1}_n$ be a set of weights such that

\[
\frac{\max_{i\in[n]} \bfw_i}{\min_{i\in[n]} \bfw_i}\leq 2
\]

and let $w = \min_{i\in[n]}\bfw_i$. Let $\bfp$ be the sampling probabilities given by Lemma <ref> for $\bfA$ and let $\bfq$ be the sampling probabilities given by Lemma <ref> for $\bfA\mid_T$. We then consider sampling probabilities $\bfr$ such that

\[
\bfr_i \coloneqq \begin{cases}
1 & \text{if $i\in \overline T$} \\
\min\braces*{\bfp_i + \bfq_i, 1} & \text{if $i\in T$}
\end{cases}
\]

Let $\bfw'\geq \mathbf{1}_n$ be a set of weights chosen randomly so that

\[
\bfw_i' \coloneqq \begin{cases}
\bfw_i/\bfr_i & \text{w.p. $\bfr_i$} \\
0 & \text{otherwise}
\end{cases}
\]

\[
\norm*{\bfy}_{H,\bfw'} = \norm*{\bfy}_{H,\bfw} \pm \eps w\rho.
\]

for all $\bfy\in\mathcal B_\rho^H$ with probability at least $1 - \delta$. Furthermore,

\begin{align*}
\E\nnz(\bfw') &= O(\kappa m\log n) = O\parens*{\kappa d\poly\parens*{\frac{\log n}{\eps}}\log\frac1\delta} \\
\nnz(\bfw') &= O\parens*{\E\nnz(\bfw')}
\end{align*}

with probability at least $1 - \delta$, and $\norm*{\bfw'}_\infty \leq O(n)\norm*{\bfw}_\infty$.

Note that for our choice of $\gamma$, we have that $\gamma = (\gamma n)^{-\beta}$. Then by our choice of $\kappa$ and Corollary <ref>, any vector $\bfy$ either has at least one of

\begin{equation}\label{eq:huber-distortion-choice}
\min_{p\in\mathcal I}\frac{\norm*{\bfy}_p^p}{\norm*{\bfy}_H^2} \leq \kappa \qquad \mbox{ or } \qquad \min_{p\in\mathcal I}\frac{\norm*{\bfy\vert_T}_p^p}{\norm*{\bfy\vert_T}_H^2} \leq \kappa.
\end{equation}

Now let $\mathcal N$ and $\mathcal N_T$ denote the nets obtained by applying Lemma <ref> on $\bfA$ and $\bfA\vert_T$, respectively, with radius $\rho$. Let $\bfy\in\mathcal B_\rho^H$. Then, clearly, $\bfy$ is within $\eps^2\rho$ Huber distance of some $\tilde\bfy$ belonging to either $\mathcal N$ or $\mathcal N_T$, so $\mathcal N\cup\mathcal N_T$ is an $\eps^2\rho$-cover of $\mathcal B_\rho^H$. Furthermore, by Lemma <ref>,

\[
\norm*{\bfy}_{H,\bfw'} = \norm*{\bfy}_{H,\bfw} \pm \eps w\rho
\]

for every $\bfy\in \mathcal N\cup\mathcal N_T$ with probability at least $1 - 2\delta$. Then by Lemma <ref> we have that

\[
\norm*{\bfy}_{H,\bfw'} = \norm*{\bfy}_{H,\bfw} \pm O(\eps)w\rho
\]

for all $\bfy\in\mathcal B_\rho^H$, with probability at least $1 - 2\delta$.

Finally, letting $m$ be as in the statement of Lemma <ref>,

\begin{align*}
\E\nnz(\bfw') &\leq \sum_{i=1}^n \bfp_i + \bfq_i \\
&\leq O(\kappa m)\parens*{\sum_{p\in\mathcal I}\bracks*{\frac1n + \sum_{i=1}^n \frac{\bfw_i^p(\bfA)}{d}}+\sum_{p\in\mathcal I}\bracks*{\frac1n + \sum_{i\in T} \frac{\bfw_i^p(\bfA\mid_T)}{d}}} \\
&= O(\kappa m\abs*{\mathcal I}) = O(\kappa m\log n)
\end{align*}

since Lewis weights sum to $d$. Furthermore, by Bernstein's inequality,

\[
\nnz(\bfw') = \Theta(\E\nnz(\bfw'))
\]

with probability at least $1 - \delta$. By rescaling $\epsilon$ and $\delta$ by constant factors, we conclude.

§.§.§ Sampling Guarantees for a Single Step

We now remove the assumption of bounded weights and fixed radius by union bounding over their various levels.

Let $\bfA\in\mathbb R^{n\times d}$ and let $\mathcal V = \Span(\bfA)$. Let $\eps\in(0,1)$ be an accuracy parameter, $\delta>0$ a failure rate parameter. Let $\bfw\geq\mathbf{1}_n$ be a set of weights and for each $j\in[\ceil*{\log_2(\norm*{\bfw}_\infty)}+1]$ define the sets

\[
T_j \coloneqq \braces*{i\in[n] : \bfw_i \in [2^{j-1}, 2^j)}.
\]

Let $\bfw'$ be obtained by applying Lemma <ref> on each $\bfA\vert_{T_j}$ with weights $\bfw\vert_{T_j}$, with $\delta$ set to

\[
\Theta\parens*{\frac{\delta}{\log(n\norm*{\bfw}_\infty/\eps}}.
\]

Then, for all $\bfy\in\mathcal V$,

\[
\norm*{\bfy}_{H,\bfw'} = (1\pm\eps)\norm*{\bfy}_{H,\bfw}.
\]

Furthermore, for $\beta = 3-2\sqrt 2$ and $\kappa = n^{\beta/(\beta+1)}$,

\begin{align*}
\E\nnz(\bfw') &= O\parens*{\kappa d\poly\parens*{\frac{\log n}{\eps}}\log\frac{\norm*{\bfw}_\infty}\delta} \\
\nnz(\bfw') &= O\parens*{\E\nnz(\bfw')}
\end{align*}

with probability at least $1 - \delta$.

For any $\rho\geq 1$, by a union bound over the $\ceil*{\log_2\norm*{\bfw}_\infty} + 1$ weight classes $T_j$, the approximation guarantee from Lemma <ref> holds for all classes with probability at least $\Theta(\delta/\log(n\norm*{\bfw}_\infty/\eps))$. We then union bound over the $\Theta(\log(n\norm*{\bfw'}_\infty/\eps)) = \Theta(\log(n\norm*{\bfw}_\infty/\eps))$ radius levels $\rho$ required by Lemma <ref> so that

\[
\norm*{\bfy}_{H,\bfw'} = (1\pm\eps)\norm*{\bfy}_{H,\bfw}
\]

for all $\bfy\in\mathcal V$.

We provide the algorithm for Lemma <ref> Algorithm <ref>, whose guarantee is proven in Theorem <ref>.

Huber subspace embedding
input: Matrix $\bfA \in \R^{n \times d}$, weights $\bfw$.

output: Weights $\bfw'$.
[1]
$\beta\gets 3 - 2\sqrt 2$, $\gamma\gets n^{-\beta/(1+\beta)}$, $m\gets d\poly((\log(n\norm*{\bfw}_\infty/\delta))/\eps), W\gets \ceil*{\log_2\norm*{\bfw}_\infty}+1$
$T_j \gets \braces*{i\in[n]:\bfw_i\in[2^{j-1},2^j)}$ for $j\in[W]$
$\mathcal I\gets\{1+(c/\log n)\cdot j : j\in[(\log n)/c]\cup\{0\}\}$
Compute approximate Huber sensitivities $\tilde\bfs_i^H(\bfA\vert_{T_j})$ with Algorithm <ref> with $\tau = \gamma n$
$S_j\gets \{i\in[n] : \tilde\bfs_i^H(\bfA\vert_{T_j})\leq \gamma\}$
Set $\bfw_i' = \bfw_i$ for every $i\notin S_j$
Compute approximate Lewis weights $\tilde\bfw_i^p(\bfA)$ with Theorem <ref> for every $p\in\mathcal I$
Compute approximate Lewis weights $\tilde\bfw_i^p(\bfA\vert_{S_j\cap T_j})$ with Theorem <ref> for every $p\in\mathcal I$
$i\in S_j$
Let $\bfp_i = \min\{1,\gamma^{-1}m(1/n + \sum_{p\in\mathcal I}\tilde\bfw_i^p(\bfA)/d) + \sum_{p\in\mathcal I}\tilde\bfw_i^p(\bfA\mid_T)/d)\}$
Set $\bfw_i' = \bfw_i/\bfp_i$ with probability $\bfp_i$ and $0$ otherwise.

Let $\bfA\in\mathbb R^{n\times d}$ and let $\mathcal V = \Span(\bfA)$. Let $\eps>0$ and $\delta > 0$. Let $\bfw\geq\mathbf{1}_n$ be a set of weights. Let $\bfw'$ be the weights returned by Algorithm <ref>. Then,

\[
\norm*{\bfy}_{H,\bfw}^2 = (1\pm O(\eps))\norm*{\bfy}_H^2
\]

for all $\bfy\in\mathcal V$ and

\begin{align*}
\E\nnz(\bfw') &= O\parens*{\gamma^{-1}d\poly\parens*{\frac{\log (n\norm*{\bfw}_\infty)}{\delta\eps}}} \\
\nnz(\bfw') &= O\parens*{\E\nnz(\bfw')}
\end{align*}

with probability at least $1 - \delta$. Furthermore, Algorithm <ref> runs in time

\[
O\parens*{(\nnz(\bfA) + d^{(1+1/\beta)\omega})(\log^2 n)\log\frac1\delta}.
\]

We first identify the rows with Huber sensitivity at least $\gamma$. By running Theorem <ref> $\log(1/\delta)$ times with $\tau = \gamma n = n^{1/(\beta+1)}$, we can boost the success probability to $1-\delta$ to find sensitivity upper bounds that sum to

\[
O\parens*{d(\log^2 n)\log\frac1\delta}
\]

in time

\[
O\parens*{\bracks*{\nnz(\bfA)\log^2 n + n^{1/(\beta+1)}T\log n}\log\frac1\delta}.
\]

for $T = O(d^{\omega}\log^2 d)$ by Theorem <ref>. Note that in the above running time, the term inside the square brackets is $O(\nnz(\bfA)\log^2 n)$ time for $T \leq n^{\beta/(\beta+1)}$ and $O((\nnz(\bfA)+T^{1+1/\beta})\log^2 n)$ time for $T \geq n^{\beta/(\beta+1)}$. Then, we may find a superset of rows with Huber sensitivity at least $\gamma$ of size at most $O\parens*{\gamma^{-1}d(\log^2 n)\log\frac1\delta}$. We may also approximate the $\ell_p$ Lewis weights in time $O([\nnz(\bfA) + T](\log n)\log\frac1\delta)$ by Theorem <ref>, for each $p\in\mathcal I$. The guarantees on $\nnz(\bfw)$ follow from Lemma <ref>.

§.§ Subspace Embeddings via Recursive Sampling

For large $n$, a bound of $\tilde O(n^{\beta/(\beta+1)}d)$ with $\beta = 3-2\sqrt 2$ of Theorem <ref> is much worse than the previous $O(d^2\poly\log(n))$ bound. However, by applying Theorem <ref> after reducing the number of rows to $O(d^2\poly\log(n))$, we can obtain a sampling bound of roughly $O(d^{2\beta/(\beta+1)}\poly\log(n))$, where $1 + 2\beta/(\beta+1)\approx 1.29289$. Furthermore, by recursively applying this procedure, we can improve the dependence on $d$ to $d^{1+\beta} = d^{4-2\sqrt 2}$ after only $O(\log\log(d\log n))$ many iterations, where $4-2\sqrt 2\approx 1.17157$.

Note that as long as we apply this recursion at most $\log n$ times, the number of rows decreases as $m \to O(m^{\beta/(\beta+1)}d\poly((\log n) / \eps))$ at each iteration, since $\log\norm*{\bfw}_\infty$ increases by at most $O(1)$ at each iteration. The closed form solution to this can be found via the following recurrence:

Suppose that $(a_i)_{i=0}^\infty$ satisfies the recurrence

\[
a_{i+1} = \lambda a_i + b
\]

for some $b > 0$ and $\lambda\in(0,1)$. Then,

\[
a_i = \frac{1}{1-\lambda} \parens*{b - \lambda^{i}\parens*{b - (1-\lambda) a_0}}
\]

Note that

\[
\frac{1}{1-\lambda} \parens*{b - \lambda^{0}\parens*{b - (1-\lambda) a_0}} = \frac{1}{1-\lambda} \parens*{b -b + (1-\lambda) a_0} = a_0
\]

\begin{align*}
\frac{1}{1-\lambda} \parens*{b - \lambda^{i+1}\parens*{b - (1-\lambda) a_0}} &= \lambda \cdot \frac{1}{1-\lambda} \parens*{\frac{1-\lambda}{\lambda}b + b - \lambda^{i}\parens*{b - (1-\lambda) a_0}} \\
&= \lambda \frac{1}{1-\lambda} \parens*{b - \lambda^{i}\parens*{b - (1-\lambda) a_0}} + b.
\end{align*}

By applying the above lemma on the logarithm of the number of rows after the $i$th recursive application of row sampling, where $b$ corresponds to $\log(d\poly((\log n) / \eps)))$, we obtain the following algorithm and theorem:

Recursive Huber subspace embedding
input: Matrix $\bfA \in \R^{n \times d}$, weights $\bfw$ (defaulted to $\mathbf{1}_n$).

output: Weights $\bfw''$.
[1]
$\beta\gets3-2\sqrt 2$
$\nnz(\bfw) \leq O(d^{1+\beta}\poly(\log (n/\delta)/\eps))$
Obtain weights $\bfw'$ by running Algorithm <ref> on $\bfA$ with weights $\bfw$
Obtain weights $\bfw''$ by recursively running Algorithm <ref> on $\bfA$ with weights $\bfw'$

There is an algorithm (Algorithm <ref>) which, with probability $1-\delta$, computes a set of weights $\bfw\geq\mathbf{1}_n$ with

\[
\nnz(\bfw) \leq O(d^{4-2\sqrt 2}\poly(\log \frac{n}{\delta},\eps^{-1})
\]

such that $\norm*{\bfA\bfx}_{H,\bfw}^2 = (1\pm \eps) \norm*{\bfA\bfx}_{H}^2$ for all $\bfx\in\mathbb R^d$, and runs in time

\[
O\parens*{\bracks*{\nnz(\bfA)(\log^2 n)(\log\log n) + d^{\omega\cdot(1+1/\beta)}\poly\parens*{\log n,\eps^{-1}}}\log\frac1\delta}.
\]

We apply Lemma <ref> for $O(\log\log n)$ recursive steps with $\delta = O(1/\log\log n)$ so that we have a sequence of upper bounds $m_i$ on the number of rows such that, for $\beta = 3-2\sqrt 2$,

\[
m_i = O\parens*{m_{i-1}^{\beta/(\beta+1)}d \poly\parens*{\frac{\log (n/\delta)}{\eps}}}.
\]

Note then that $a_i\coloneqq \log m_i$ satisfies

\[
a_i = \frac{\beta}{\beta+1} a_{i-1} + b
\]

\[
b = \log\bracks*{O\parens*{d\poly\parens*{\frac{\log (n/\delta)}{\eps}}}}.
\]

The closed form solution for this is given in Lemma <ref>. Note that if $i$ satisfies

\[
\parens*{\frac{\beta}{\beta+1}}^{i}\parens*{b + \frac{1}{\beta+1} a_0} \leq \frac1{\log b}
\]

\[
a_i = (1+\beta) b \pm \frac1{\log b}
\]

and thus the number of rows is at most

\[
m_i = \exp(a_i) = O(\exp(b)^{1+\beta}) = O(d^{1+\beta}\poly((\log n)/\eps)) =  O(d^{4-2\sqrt2}\poly((\log n)/\eps)).
\]

This only requires $i$ at most $O(\log(b + a_0)) = O(\log\log n)$ iterations.

§.§.§ Running Time Trade-offs

The above running time is quite a large polynomial in $d$; for the current matrix multiplication exponent of $\omega\approx 2.37286$, the exponent of $d$ is

\[
\omega\cdot \parens*{1+\frac1\beta} \approx  16.20290.
\]

However, by a further multi-step algorithm, we can further reduce the running time significantly, in a similar manner to Theorem <ref>. First note that by always only computing sensitivities which sum to at most $d^\omega\poly\log n$, we can achieve roughly $\nnz(\bfA) + d^\omega$ time to reduce the number of rows to roughly $n' = d^{\omega(1+\beta)}$ with the same proof, where we solve the recursion of reducing the dimension from $n$ to $n^{\beta/(1+\beta)}d^\omega$ at each step. We can then spend roughly $\nnz(\bfA) + n' d^\omega / d = \nnz(\bfA) + d^{\omega(2+\beta) - 1}$ in Theorem <ref> to fully reduce the dimension to $d^{1+\beta}$, where the exponent in $d$ is now only

\[
\omega\cdot\parens*{2+\beta} - 1 \approx 4.14663.
\]

This running time may still be undesirable, if one is willing to sacrifice in the dimension reduction bound; for example, this does not beat the running time $\nnz(\bfA) + d^4$ in [34], even though the dimension reduction bound is significantly better. If we wish for a running time that is strictly $\nnz(\bfA) + d^\omega$, up to polylogarithmic factors, then the previous argument shows that we obtain a dimension reduction bound of

\[
\omega\cdot(1+\beta) \approx 2.77998
\]

rows, which indeed slightly improves over their bound of $d^3$ rows in the same running time.

To find an intermediate trade-off by balancing the running time and row count, let $C \in [1, \omega]$ be a parameter. Assume we have already reduced to $n = d^{\omega(1+\beta)}$ rows in $\nnz(\bfA)+d^\omega$ time. Then, we can spend an additional

\[
\frac{nd^\omega}{d^C} = d^{\omega(2+\beta) - C}
\]

time to reduce to $d^{C(1+\beta)}$ rows. One interesting choice is to balance the running time to be equal to $d$ times the row count, which gives $C = \frac{\omega(2+\beta)-1}{2+\beta} \approx 1.91236$ which allows for reduction to $d^C$ rows in $d^{C+1}$ time.

§.§ Active Regression Algorithms

As noted before, Theorem <ref> gives constant factor solutions in the active regression setting, by Lemma <ref>. We now discuss how to build on this result to obtain a relative error solution, using the more refined techniques of Section <ref>.

Our main result of this section is the following:

Active regression with the Huber loss
input: Matrix $\bfA \in \R^{n \times d}$, $\bfb\in\mathbb R^n$, weights $\bfw\geq\mathbf{1}_n$.

output: Approximate solution $\tilde\bfx\in\R$ to $\min_\bfx\norm*{\bfA\bfx-\bfb}_H$.
[1]
$\nnz(\bfw) \leq d^{4-2\sqrt2}\poly(\log n,\eps^{-1})$
Compute an approximate solution $\tilde\bfx$ such that $\norm*{\bfA\tilde\bfx-\bfb}_{H,\bfw} \leq (1+\eps)\min_\bfx\norm*{\bfA\bfx-\bfb}_{H,\bfw}$
Let $\bfw'$ be generated by Algorithm <ref>
Let $\bfx_c$ be an approximate minimizer $\norm*{\bfA\bfx_c-\bfb}_{H,\bfw'} \leq O(1)\min_\bfx\norm*{\bfA\bfx-\bfb}_{H,\bfw'}$
Set $\bfz\gets\bfb - \bfA\bfx_c$ Implicit, for analysis
Let $\bfw''$ be generated by Algorithm <ref>, with possibly larger $\poly(\log(n/\delta)/\eps)$ terms
Recursively compute an approximate solution $\tilde\bfx$ with inputs $(\bfA,\bfz,\bfw'')$ with Algorithm <ref>
$\tilde\bfx + \bfx_c$

We first obtain an analog of Lemmas <ref> and <ref>, mirroring some of the logic from Lemma <ref>.

Let $\bfw\geq\mathbf{1}_n$ be a set of weights. Let $\bfA\in\R^{n\times d}$. Let $\bfz\in\mathbb R^n$ be such that $\norm*{\bfz}_{H,\bfw} \leq O(\OPT)$, where

\[
\OPT \coloneqq \min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx-\bfz}_{H,\bfw}.
\]

Let $\tilde\bfs_i^{H,\bfw}(\bfA)\geq\bfs_i^{H,\bfw}(\bfA)$ be an upper bound on the Huber sensitivities. Let $\eps>0$ and consider a subset

\[
\mathcal B\subseteq\braces*{i\in[n] : \bfw_i H(\bfz_i) \geq \frac{\tilde\bfs_i^{H,\bfw}(\bfA)}{\eps^2}\OPT^2}
\]

and let $\bar\bfz\in\mathbb R^n$ be equal to $\bfz$, but with all entries in $\mathcal B$ set to $0$. Then, for all $\bfx\in\mathbb R^d$ with $\norm*{\bfA\bfx}_{H,\bfw} \leq O(\OPT)$,

\[
\abs*{\norm*{\bfA\bfx-\bar\bfz}_{H,\bfw}^2-\norm*{\bfA\bfx-\bfz}_{H,\bfw}^2 - \norm*{\bfz-\bar\bfz}_{H,\bfw}^2} = O(\eps)\OPT^2.
\]

For any $\bfx\in\R^d$ such that $\norm*{\bfA\bfx}_{H,\bfw} \leq O(\OPT)$. Then, for any $i\in\mathcal B$,

\begin{equation}\label{eq:huber-small-entries}
\begin{aligned}
\bfw_i H([\bfA\bfx](i)) &\leq \tilde\bfs_i^{H,\bfw}(\bfA)\norm*{\bfA\bfx}_{H,\bfw}^2 \\
&\leq \tilde\bfs_i^{H,\bfw}(\bfA)O(\OPT^2) \\
&\leq O(\eps^2)\cdot\bfw_i H(\bfz_i).
\end{aligned}
\end{equation}

\begin{align*}
&\bfw_i H([\bfA\bfx](i)-\bfz(i)) - \bfw_i H([\bfA\bfx](i)-\bar\bfz(i)) \\
=~&\bfw_i H([\bfA\bfx](i)-\bfz(i)) - \bfw_i H([\bfA\bfx](i)) && i\in\mathcal B \\
=~&\bfw_i\parens*{H^{1/2}(\bfz(i))\pm H^{1/2}([\bfA\bfx](i))}^2 - \bfw_i H([\bfA\bfx](i)) && \text{see Lemma \ref{lem:triangle-inequality-m}} \\
=~&(1\pm O(\eps))\bfw_i H(\bfz_i) && \text{Equation \eqref{eq:huber-small-entries}}
\end{align*}

\begin{align*}
&\abs*{\norm*{\bfA\bfx-\bfz}_{H,\bfw}^2 - \norm*{\bfA\bfx-\bar\bfz}_{H,\bfw}^2 - \norm*{\bfz-\bar\bfz}_{H,\bfw}^2} \\
=~&\abs*{\sum_{i\in\mathcal B}(1\pm O(\eps))\bfw_i H(\bfz_i) - \bfw_i H(\bfz_i)} \\
=~&\abs*{\sum_{i\in\mathcal B} O(\eps)\bfw_i H(\bfz_i)} \\
=~&O(\eps)\norm*{\bfz-\bar\bfz}_{H,\bfw}^2 \leq O(\eps)\norm*{\bfz}_{H,\bfw}^2 \leq O(\eps)\OPT^2.\qedhere
\end{align*}

Next, we show that our Huber subspace embedding algorithm approximately preserves Huber norms $\norm*{\bfA\bfx-\bar\bfz}_{H,\bfw}^2$ for all $\bfx\in\mathbb R^d$ such that $\norm*{\bfA\bfx}_{H,\bfw} = O(\OPT)$.

Consider the setting of Lemma <ref>. Suppose that $\bfw'$ is obtained by performing one step of Algorithm <ref> (with the $\poly(\log(n\norm*{\bfw}_\infty/\delta)/\eps)$ term in $m$ possibly being a larger polynomial). Then, there exists a $\mathcal B$ satisfying the requirements of Lemma <ref> (see Equation (<ref>)) such that, with probability at least $1-\delta$, for all $\bfx\in\mathbb R^{d}$ with $\norm*{\bfA\bfx}_{H,\bfw} = O(\OPT)$,

\begin{equation}\label{eq:huber-additive2-guarantee}
\abs*{\norm*{\bfA\bfx-\bar\bfz}_{H,\bfw'}^2 - \norm*{\bfA\bfx-\bar\bfz}_{H,\bfw}^2} = O(\eps)\cdot\OPT^2.
\end{equation}

We will make a careful modification of our Huber rounding technique of Section <ref>. To adapt the rounding to the active setting, we will need to modify the proof for the [6] net construction rather than using them as a black box.

Note that we only need to provide Huber norm preservation guarantees for a single scale of Huber radius $\rho = O(\OPT)$ for this lemma, as opposed to $O(\log n)$ scales for the Huber subspace embedding result. Furthermore, by Lemma <ref>, it suffices to show the approximation guarantee for every $\bfy\in\mathcal N\cup(\mathcal N-\bar\bfz)$, where $\mathcal N$ is an $O(\eps^2\rho)$-cover of $\mathcal B_\rho^{H,\bfw}$. We have already shown in Theorem <ref> the preservation of $\norm*{\bfy}_{H,\bfw}^2$ for every $\bfy$ in $\mathcal N$, so we focus on showing the preservation of $\norm*{\bfy}_{H,\bfw}^2$ for every $\bfy\in\mathcal N-\bar\bfz$, which is similar but requires a slightly more involved argument to handle the $\bar\bfz$ translation.

We start off with the following claim, which is the active version of Lemma <ref>.

Let $\mathcal I\subset[1,2]$ be the set in the statement of Corollary <ref>. Let $\kappa \geq 1$ be a distortion parameter. Suppose that $\tilde\bfs_i^{H}(\bfA)\geq\bfs_i^{H}(\bfA)$ is an upper bound on the Huber sensitivities such that

\begin{equation}\label{eq:huber-lp-sensitivity-bound}
\kappa\bfw_i^p(\bfA) \leq \tilde\bfs_i^{H}(\bfA)
\end{equation}

for every $p\in\mathcal I$. Let $\Delta>0$ be a scale parameter and let $\bfv\in\mathbb R^n$ be a target vector satisfying

\[
H(\bfv(i)) \leq \frac{\tilde\bfs_i^{H}(\bfA)}{\eps^4}\Delta^2
\]

for each $i\in[n]$.
Let $\bfr = \bfy - \bfv\in (\mathcal N - \bfv)$ for for $\mathcal N$ an $O(\eps^2\rho)$-cover of $\mathcal B_\rho^H$ for $\rho = O(\Delta)$. Suppose that $\bfy$ satisfies

\[
\frac{\min_{p\in\mathcal I}\norm*{\bfy}_p^p}{\norm*{\bfy}_H^2} \leq \kappa.
\]

Then, there exists a $p\in\mathcal I$, $\ell = O((\log d)/\eps^2)$, and a rounding $\bfr' = \bfe + \sum_{k=0}^\ell \bfd_k$ such that:

* $\norm*{\bfr' - \bfr}_H \leq O(\eps)\Delta$

* $H\parens*{\abs*{\bfd_k(i)}} \leq O(1)\frac{\tilde\bfs_i^H(\bfA)\cdot(1+\eps^2)^{pk}}{\eps^4\cdot d}\Delta^2$

* $\bfe, \bfd_0, \bfd_1, \dots, \bfd_\ell$ have disjoint supports

* $\bfe$ is a single fixed vector with $H(\abs{\bfe(i)}) \leq O(1)\frac{\tilde\bfs_i^H(\bfA)}{\eps^4\cdot d}\Delta^2$

* each $\bfd_k$ is drawn from a set of vectors $\mathcal D_k$ with size at most

\[
\log\abs*{\mathcal D_k} \leq O(1)\frac{d}{\eps^{2(1+p)}(1+\eps^2)^{pk}}\log\frac{n}{\eps}.
\]

Let $p\in\mathcal I$ be the $p$ achieving the minimum in Corollary <ref> for $\bfy$. We will now define $\ell$ nets $\mathcal N_k$ for $k\in[\ell]$ as follows. Recall the following result from [6] and [63]:

Let $\bfA\in\mathbb R^{n\times d}$ and $1 \leq p < \infty$. Let $\bfw\in\mathbb R^{n\times n}$ be the diagonal matrix with

\[
\bfW_{i,i} \coloneqq \frac12\parens*{\frac{\bfw_i^p(\bfA)}{d} + \frac1n.}
\]

Let $\mathcal B_p \coloneqq \braces*{\bfA\bfx : \norm*{\bfA\bfx}_p\leq 1}$. Then for any $t > 0$, there is a net $\mathcal N_\infty\subseteq\mathbb R^n$ such that, for any $\bfy\in\mathcal B_p$, there exists a $\bfy\in\mathcal N_\infty$ with $\norm*{\bfW^{-1/p}(\bfy-\bfy')}_\infty \leq t$ and

\[
\log\abs*{\mathcal N_\infty} \leq c(p) \frac{d}{t^{2\land p}}\log n
\]

where $c(p)$ is a constant depending only on $p$ <cit.>.

We use the following definition of a sequence of nets from [6].

[Index sets [6]]
Let $\eps' > 0$. For $k\in[\ell]$, let $\mathcal N_k$ be the net from Lemma <ref> with $t = \eps'(1+\eps')^k / 3$. If $t \leq 1$, we take $\mathcal N_k = \mathcal N$. Now for each $\bfy\in\mathcal N$ and $k\in[\ell]$, let $\bff_{k,\bfy}\in\mathcal N_k$ satisfy

\[
\norm*{\bfW^{-1/p}(\bff_{k,\bfy}-\bfy)}_\infty \leq \eps'(1+\eps')^k / 3
\]

as given by Lemma <ref>. Define the index sets

\begin{align*}
C_{k,\bfy} &\coloneqq \braces*{i\in[n] : |\bfW^{-1/p} \bff_{k,\bfy}(i)| \ge (1+\epsilon')^{k-1}} \\
D_{k,\bfy} &\coloneqq C_{k,\bfy} \setminus \bigcup_{k' > k} C_{k',\bfy} \\
D_{0,\bfy} &\coloneqq [n] \setminus \bigcup_{k \ge 1} C_{k,\bfy}
\end{align*}

Then, for each $k$ we have $\log |\mathcal{N}_k| = O \left (\frac{d\log n}{(\epsilon'(1+\epsilon')^{k})^{2\land p}} \right )$, and for every $i\in D_{k,\bfy}$, we have

\begin{align}\label{yi_bound}
\frac{\bfw_i^p(\bfA)^{1/p} \cdot (1+\epsilon')^{k-2}}{d^{1/p}} \le |\bfy(i)| \le \frac{\bfw_i^p(\bfA)^{1/p} \cdot (1+\epsilon')^{k+2}}{d^{1/p}}.
\end{align}

Note that since $\norm{\bfW^{-1/p}(\bff_{k,\bfy}-\bfy)}_\infty \le \epsilon'(1+\epsilon')^k/3$, if $i \in C_{k,\bfy}$ for $1 \le k \le \ell$,

\begin{align*}
|\bfy(i)| &\ge |\bff_{k,\bfy}(i)| - \frac{\bfw^p_i(\bfA)^{1/p} \cdot \epsilon' (1+\epsilon')^k/3}{d^{1/p}}\\
&\ge \frac{\bfw^p_i(\bfA)^{1/p}}{d^{1/p}} \cdot \left [(1+\epsilon')^{k-1} - \epsilon' (1+\epsilon')^k/3 \right ]\tag{By definition, if $i \in C_{k,\bfy}$, $|\bff_{k,\bfy}(i)| \ge \frac{\bfw^p_i(\bfA)^{1/p} \cdot (1+\epsilon')^{k-1}}{d^{1/p}}$.}\\
&\ge \frac{\bfw^p_i(\bfA)^{1/p} \cdot (1+\epsilon')^{k-2}}{d^{1/p}}.
\end{align*}

If $i \notin C_{k,\bfy}$ we analogously have $|\bfy(i)| \le \frac{\bfw_i^p(\bfA)^{1/p} \cdot (1+\epsilon')^{k+2}}{d^{1/p}}$.
Thus, for $i \in D_{k,\bfy} = C_{k,\bfy} \setminus \bigcup_{k' > k} C_{k',\bfy}$ for $1 \le k < \ell$, Equation (<ref>) holds. Equation (<ref>) also holds for $k = \ell$ via Lemma <ref> – since $\norm{\bfy}_p \le 1$, $|\bfy(i)| \le d^{\max(0,1/2-1/p)} \cdot \bfw_i^p(\bfA)^{1/p}$ for all $i$. And so vacuously, since for $\ell = \log_{1+\epsilon'} d^{\max(1/p,1/2 )}$, $(1+\epsilon')^{\ell+2} \ge d^{\max(1/p,1/2)}$, we have $|\bfy(i)| \le \frac{\bfw_i^p(\bfA)^{1/p} \cdot (1+\epsilon')^{\ell + 2}}{d^{1/p}}$.

We now apply Claim <ref> with $\eps' = \eps^2$ to obtain the nets $\mathcal N_k$. We then define index sets for the target vector $\bfv$ as follows:

[Huber index sets for $\bfv$]
Define the following sets:

\begin{align*}
B_{k,\bfy} &\coloneqq \braces*{i \in D_{k,\bfy}: H(\bfv(i)) \le \frac{\tilde\bfs_i^H(\bfA) \cdot (1+\epsilon^2)^{p(k+2)}}{\epsilon^4 \cdot d}\Delta^2} \tag{$k\in[\ell]\cup\{0\}$} \\
H_k &\coloneqq \braces*{i\in[n]: \frac{\tilde\bfs_i^H(\bfA)\cdot (1+\epsilon^2)^{p(k+1)}}{\epsilon^4 \cdot d}\Delta^2 < H(\bfv(i))\le \frac{\tilde\bfs_i^H(\bfA) \cdot (1+\epsilon^2)^{p(k+2)}}{\epsilon^4 \cdot d}\Delta^2 } \tag{$k\in[\ell]$} \\
H_B &\coloneqq \braces*{ i\in[n]: H(\bfv(i)) > \frac{\tilde\bfs_i^H(\bfA) \cdot (1+\epsilon^2)^{p(\ell+2)}}{\epsilon^4 \cdot d}\Delta^2 } \\
G_{k,\bfy} &\coloneqq H_{k} \setminus \bigcup_{k' \ge k} C_{k',\bfy} \tag{$k\in[\ell]$}
\end{align*}

Then, $B_{0,\bfy},\ldots,B_{\ell,\bfy},G_{1,\bfy},\ldots,G_{\ell,\bfy}, H_B$ form a partition of $[n]$.

Define $H_{0} = \{i: H(\bfv(i)) \le \frac{\tilde\bfs_i^H(\bfA) \cdot (1+\epsilon^2)^{2p}}{\epsilon^4 \cdot d}\Delta^2 \}$, so that $H_B \cup \bigcup_{k=0}^\ell H_{k} = [n]$ is a partition the range of all possible $H(\bfv(i))$ values. Furthermore, every $i\in H_0$ lies in some $B_{k,\bfy}$.

At a high level, our argument will be as follows. $H_k$ can be furthered partitioned into the indices that lie in some $\bigcup_{k'\geq k}C_{k,\bfy}$, or those that are not in any $\bigcup_{k'\geq k} C_{k',\bfy}$. Those that lie in some $\bigcup_{k'\geq k}C_{k,\bfy}$ are then regrouped into the $B_{k,\bfy}$, which is also a valid partition by definition of the $D_{k,\bfy}$.

Any $i \in H_{k} \cap \left (\bigcup_{k' \ge k} C_{k',\bfy} \right)$ has $H(\bfv(i))\le \frac{\tilde\bfs_i^H(\bfA) \cdot (1+\epsilon^2)^{p(k+2)}}{\epsilon^4 \cdot d}\Delta^2$ and thus must lie in $B_{k',\bfy}$ for some $k' \ge k$.
Thus, $G_{k,\bfy} \cup (\bigcup_{k' \ge k} B_{k',\bfy} ) \supseteq H_{k}$ and so overall $(\bigcup_{k=1}^\ell G_{k,\bfy} \cup B_{k,\bfy}) \cup B_{0,\bfy} \cup H_B = [n]$.

Now clearly, the $H_{k}$ and in turn the $G_{k,\bfy}$ are disjoint. $H_B$ is disjoint from all $H_k$ and hence all $G_{k,\bfy}$. It is also disjoint from all $B_{k,\bfy}$ since the $H(|\bfv(i)|)$ values in $H_B$ are too large to be assigned to $B_{k,\bfy}$, even for $k = \ell$. The $D_{k,\bfy}$, and in turn the $B_{k,\bfy}$, are also disjoint by construction. Further, $\bigcup_{k' \ge k} B_{k',\bfy} \subseteq \bigcup_{k' \ge k} D_{k',\bfy} = \bigcup_{k' \ge k} C_{k',\bfy}$, and thus $G_{k,\bfy}$ is disjoint from $B_{k',\bfy}$ for all $k' \ge k$. Also, for $k \ge 1$, being in $G_{k,\bfy}$ requires $H(|\bfv(i)|) > \frac{\tilde\bfs_i^H(\bfA)\cdot (1+\epsilon^2)^{p(k+1)}}{\epsilon^4 \cdot d}\Delta^2$, while being in $B_{k',\bfy}$ for $k' < k$ requires $H(\bfv(i)) \le \frac{\tilde\bfs_i^H(\bfA) \cdot (1+\epsilon^2)^{p(k+1)}}{\epsilon^4 \cdot d}\Delta^2$. Thus $G_{k,\bfy}$ is disjoint from $B_{k',\bfy}$ for all $k' < k$. Overall, $G_{k,\bfy}$ is disjoint from all $B_{k',\bfy}$. Thus, $B_{0,\bfy},\ldots,B_{\ell,\bfy},G_{1,\bfy},\ldots,G_{\ell,\bfy}, H_B$ are all mutually disjoint and partition $[n]$.

Using the index sets of Claim <ref>, we next show how to round $\bfr = \bfy-\bfv$ to a nearby $\bfr'$ using the partition of $[n]$ defined above.

[$\ell_\infty$ Error Bound]
Define the vectors $\bfr'$, $\bfe$, and $\bfd_k$ with

\[
\bfr' = \bfe + \sum_{k=0}^\ell \bfd_k
\]

as follows:

\begin{equation}\label{eq:active-rounding-definition}
\begin{aligned}
\bfd_0(i) &\coloneqq \bfy(i)-\bfv(i) && i\in B_{0,\bfy} \\
\bfd_k(i) &\coloneqq \frac{(1+\epsilon^2)^k \cdot \bfw_i^p(\bfA)^{1/p}}{d^{1/p}} - \bfv(i) && k\in[\ell], i\in B_{{k,\bfy}} \\
\bfd_k(i) &\coloneqq \bfv(i) && i\in G_{k,\bfy} \\
\bfd_k(i) &\coloneqq 0 && \text{otherwise} \\
\bfe_i(i) &\coloneqq \bfv(i) && i\in H_B \\
\bfe_i(i) &\coloneqq 0 && \text{otherwise}
\end{aligned}
\end{equation}

$$|\bfr(i) - \bfr'(i)| \le \epsilon^2( |\bfy(i)| + |\bfv(i)|).$$

For $i \in B_{0,\bfy}$ we have $\bfr(i) = \bfr'(i)$ so the claim trivially holds. For $i \in B_{k,\bfy}$ for $k \ge 1$, by (<ref>),

\[
\left |\bfy(i) - \frac{(1+\epsilon^2)^k \cdot \bfw_i^p(\bfA)^{1/p}}{d^{1/p}} \right | \le O(\epsilon^2) \cdot |\bfy(i)|.
\]

Thus by triangle inequality, $|\bfr(i) - \bfr'(i)| \le O(\epsilon^2) \cdot |\bfy(i)|$. So the claim holds after adjusting constants on $\epsilon$.
Finally, for $i\in G_{k,\bfy}$ we have

\begin{align*}
H(\abs*{\bfy(i)}) &\leq \abs*{\bfy(i)}^p \tag{$p\in[1,2]$} \\
&\leq \frac{\bfw_i^p(\bfA)\cdot (1+\eps^2)^{p(k+2)}}{d}\norm*{\bfy}_p^p \tag{$i\notin \bigcup_{k'\geq k}C_{k',\bfy}$} \\
&\leq \kappa\frac{\bfw_i^p(\bfA)\cdot (1+\eps^2)^{p(k+2)}}{d}\norm*{\bfy}_H^2\\
&\leq \frac{\tilde\bfs_i^H(\bfA)\cdot (1+\eps^2)^{p(k+2)}}{d}\cdot O(\Delta^2) \tag{Equation \eqref{eq:huber-lp-sensitivity-bound}} \\
&\leq O(\eps^4)H(\bfv(i)) \tag{$i\in H_k$}
\end{align*}

and similarly for $i\in H_B$.

We then have that

\begin{align*}
\sum_{i=1}^n H(\abs*{\bfr'(i)-\bfr(i)}) &\leq \sum_{i=1}^n  H(\eps^2(\abs{\bfy(i)} + \abs*{\bfv(i)})) \tag{monotonicity} \\
&\leq O(\eps^2)\sum_{i=1}^n H(\abs{\bfy(i)} + \abs*{\bfv(i)}) \tag{at least linear growth} \\
&= O(\eps^2)(\norm*{\bfy}_{H} + \norm*{\bfv}_{H})^2 = O(\eps^2)\Delta^2 \tag{triangle inequality}
\end{align*}

so $\norm*{\bfr'-\bfr}_{H} \leq O(\eps)\Delta$. The disjointness of the $\bfe$ and $\bfd_k$ as well as the net size bounds carry over directly. Finally, the bounds on the coordinates of $\bfd_k$ follow from

\begin{align*}
H(\abs{\bfd_k(i)}) &\leq \abs*{\bfd_k(i)}^p \tag{$p\in[1,2]$} \\
&\leq O(1)(\abs{\bfy(i)}^p + \abs{\bfv(i)}^p) \\
&\leq O(1)\bracks*{\frac{\bfw_i^p(\bfA)}{d}(1+\eps^2)^{pk}\norm*{\bfy}_p^p + \frac{\tilde\bfs_i^H(\bfA)\cdot(1+\eps^2)^{p(k+2)}}{\eps^4\cdot d}\Delta^2} \tag{$i\in B_{k,\bfy}$} \\
&\leq O(1)\bracks*{\kappa\frac{\bfw_i^p(\bfA)}{d}(1+\eps^2)^{pk}\norm*{\bfy}_H^2 + \frac{\tilde\bfs_i^H(\bfA)\cdot(1+\eps^2)^{p(k+2)}}{\eps^4\cdot d}\Delta^2} \\
&= O(1)\frac{\tilde\bfs_i^H(\bfA)\cdot(1+\eps^2)^{pk}}{\eps^4\cdot d}\Delta^2
\end{align*}

for $i\in B_{k,\bfy}$, and similarly for other cases as well as for $\bfe$.

This completes the proof of Claim <ref>.

We return to proving Lemma <ref>. Recall that we reduced our task to proving the approximation guarantee Equation (<ref>) for all vectors $\bfr'$ belonging to an $O(\eps^2\rho)$-cover over the set $\mathcal B_\rho^H - \bar\bfz$. We will do so using the cover given by Claim <ref>, with $\eps$ set to $\eps^2$. Furthermore, we will set the Huber sensitivity upper bounds $\tilde\bfs_i^H(\bfA)$ to be at least $\kappa\sum_{p\in\mathcal I} \tilde\bfw_i^p(\bfA)$. This indeed satisfies the hypothesis (<ref>) of Claim <ref>. Then, performing a Bernstein bound gives the active version of Lemma <ref>.

Consider the setting of Claim <ref>. Let $m = d\poly((\log(n/\delta))/\eps)$. Let $\bfw\geq \mathbf{1}_n$ be a set of weights such that

\[
\frac{\max_{i\in[n]} \bfw_i}{\min_{i\in[n]} \bfw_i}\leq 2
\]

and let $w = \min_{i\in[n]}\bfw_i$. Let $\bfw'$ be chosen randomly so that

\[
\bfw_i' \coloneqq \begin{cases}
\bfw_i/\bfp_i & \text{w.p. $\bfp_i$} \\
0 & \text{otherwise}
\end{cases}
\]

\[
\bfp_i \geq \min\braces*{1, m\cdot \frac{\tilde\bfs_i^H(\bfA)}{d}}.
\]

\[
S_\kappa \coloneqq \braces*{\bfy\in \mathcal B_\rho^H-\bfv : \frac{\min_{p\in\mathcal I}\norm*{\bfy}_p^p}{\norm*{\bfy}_H^2}\leq \kappa}
\]

Then with probability at least $1 - \delta$,

\[
\norm*{\bfy}_{H,\bfw} = (1\pm\eps)\norm*{\bfy}_H
\]

for all $\bfy\in \mathcal N\cap S_\kappa$, where $\mathcal N$ is an $\eps^2\rho$-net over $\mathcal B_\rho^H-\bfv$ given by Claim <ref>, by setting $\eps$ in the lemma to $\eps^2$.

The proof follows that of Lemma <ref> almost exactly, with minimal changes. The only necessary change is that the coordinate upper bounds of Claim <ref> are now proportional to $\tilde\bfs_i^H(\bfA)$ rather than $\bfw_i^p(\bfA)$. However, this is not a problem, since we are sampling with these probabilities, and thus does not change the proof. Note also that the sensitivity upper bounds incur extra $\poly(\eps)$ factors, which changes the $m$ by $\poly(\eps)$ factors.

Next, we use Claim <ref> to obtain an active version of Lemma <ref>.

Consider the setting of Claim <ref>. Let $\bfz\in\mathbb R^n$ with $\norm*{\bfz}_H = O(\Delta)$. Let $\gamma$ be as defined in Lemma <ref> and let $\kappa = 1/\gamma$. Let $\tilde\bfs_i^H(\bfA)$ be a Huber sensitivity upper bound and let

\[
T \subseteq \braces*{i\in[n] : \tilde\bfs_i^H(\bfA) \leq \gamma}.
\]

We now define a second set of Huber sensitivity upper bounds $\overline{\bfs_i}^H(\bfA)$ via

\[
\overline{\bfs_i}^H(\bfA) \coloneqq \begin{cases}
1 & \text{if $i\in \overline T$} \\
\kappa\sum_{p\in\mathcal I} \tilde\bfw_i^p(\bfA) + \tilde\bfw_i^p(\bfA\vert_T) & \text{if $i\in T$}
\end{cases}
\]

Then consider a superset

\[
\mathcal B \supseteq \braces*{i\in[n] : H(\bfz_i) \geq \frac{\overline{\bfs_i}^H(\bfA)}{\eps^4}\Delta^2}
\]

and let $\bar\bfz\in\mathbb R^n$ be equal to $\bfz$, but with all entries in $\mathcal B$ set to $0$. Let $\bfp$ be the sampling probabilities given by Lemma <ref> for $\bfA$ and let $\bfq$ be the sampling probabilities given by Claim <ref> for $\bfA\mid_T$. We then consider sampling probabilities $\bfr$ such that

\[
\bfr_i \coloneqq \begin{cases}
1 & \text{if $i\in \overline T$} \\
\min\braces*{\bfp_i + \bfq_i, 1} & \text{if $i\in T$}
\end{cases}
\]

Let $\bfw\geq \mathbf{1}_n$ be a set of weights such that

\[
\frac{\max_{i\in[n]} \bfw_i}{\min_{i\in[n]} \bfw_i}\leq 2
\]

and let $w = \min_{i\in[n]}\bfw_i$. Let $\bfw'$ be chosen randomly so that

\[
\bfw_i' \coloneqq \begin{cases}
\bfw_i/\bfr_i & \text{w.p. $\bfr_i$} \\
0 & \text{otherwise}
\end{cases}
\]

\[
\norm*{\bfy}_{H,\bfw'} = \norm*{\bfy}_{H,\bfw} \pm \eps w\rho.
\]

for all $\bfy\in\mathcal B_\rho^H - \bar\bfz$ with probability at least $1 - \delta$. Furthermore,

\begin{align*}
\E\nnz(\bfw') &= O(\kappa m\log n) = O\parens*{\kappa d\poly\parens*{\frac{\log n}{\eps}}\log\frac1\delta} \\
\nnz(\bfw') &= O\parens*{\E\nnz(\bfw')}
\end{align*}

with probability at least $1 - \delta$, and $\norm*{\bfw'}_\infty \leq O(n)\norm*{\bfw}_\infty$.

We may check that $\bar\bfz$ satisfies the hypothesis for $\bfv$ in Claim <ref>, both for $\bfA$ and for $\bfA\vert_T$. After noticing this fact, we simply follow the proof of Lemma <ref>, combined with Claim <ref>.

Finally, we adapt Lemma <ref> to the following claim:

Let $\bfw\geq\mathbf{1}_n$ be a set of weights and for each $j\in[\ceil*{\log_2(\norm*{\bfw}_\infty)}+1]$ define the sets

\[
T_j \coloneqq \braces*{i\in[n] : \bfw_i \in [2^{j-1}, 2^j)}.
\]

Let $\bfz\in\mathbb R^n$ with $\norm*{\bfz}_{H,\bfw} = O(\OPT)$. Now define a set of weighted Huber sensitivity upper bounds $\tilde\bfs_i^{H,\bfw}(\bfA)$ via

\[
\tilde\bfs_i^{H,\bfw}(\bfA) \coloneqq 2\cdot\overline{\bfs_i}^H(\bfA\vert_{T_j})
\]

for each $i\in T_j$, where $\overline{\bfs_i}^H(\bfA\vert_{T_j})$ is as defined in Claim <ref>. Then consider the set

\begin{equation}\label{eq:huber-large-entries-active}
\mathcal B = \braces*{i\in[n] : \bfw_i H(\bfz_i) \geq \frac{\tilde\bfs_i^{H,\bfw}(\bfA)}{\eps^4}\OPT^2}
\end{equation}

and let $\bar\bfz\in\mathbb R^n$ be equal to $\bfz$, but with all entries in $\mathcal B$ set to $0$. Let $\bfw'$ be obtained by applying Claim <ref> on each $\bfA\vert_{T_j}$ with weights $\bfw\vert_{T_j}$. Then,

\[
\norm*{\bfy}_{H,\bfw'} = \norm*{\bfy}_{H,\bfw} \pm \eps\cdot\OPT
\]

for all $\bfy\in\mathcal B_\rho^H - \bar\bfz$ with probability at least $1 - \delta$. Furthermore, for $\beta = 3-2\sqrt 2$ and $\kappa = n^{\beta/(\beta+1)}$,

\begin{align*}
\E\nnz(\bfw') &= O\parens*{\kappa d\poly\parens*{\frac{\log(n\norm*{\bfw}_\infty/\delta)}{\eps}}} \\
\nnz(\bfw') &= O\parens*{\E\nnz(\bfw')}
\end{align*}

with probability at least $1 - \delta$.

First note that $\tilde\bfs_i^{H,\bfw}(\bfA)$ are valid sensitivity upper bounds, since for any $\bfx\in\mathbb R^d$,

\[
\frac{\bfw_i H([\bfA\bfx](i))}{\norm*{\bfA\bfx}_{H,\bfw}^2} \leq \frac{\bfw_i H([\bfA\bfx](i))}{2^{j-1}\norm*{\bfA\vert_{T_j}\bfx}_{H}^2} \leq 2\frac{H([\bfA\bfx](i))}{\norm*{\bfA\vert_{T_j}\bfx}_{H}^2} \leq 2\cdot\overline{\bfs_i}^H(\bfA\vert_{T_j}).
\]

Since $\bfw_i \geq1$, we have that

\[
\braces*{i\in T_j : H(\bfz_i) \geq \frac{\overline{\bfs_i}^{H}(\bfA\vert_{T_j})}{\eps^4}\OPT^2} \subseteq \mathcal B
\]

for each $j$. Thus, restricting $\bar\bfz$ to $T_j$ satisfies the hypotheses of Claim <ref>. We then set $\Delta = \OPT/2^{j-1}$ and apply Claim <ref> so that for all $\bfy\in\mathcal B_\rho^{H,\bfw}-\bar\bfz$ with $\rho = O(\OPT)$,

\[
\norm*{\bfy\vert_{T_j}}_{H,\bfw'} = \norm*{\bfy\vert_{T_j}}_{H,\bfw} \pm \eps \cdot\OPT
\]

\[
\norm*{\bfy\vert_{T_j}}_{H,\bfw'}^2 = \norm*{\bfy\vert_{T_j}}_{H,\bfw}^2 + \eps^2\cdot\OPT^2\pm 2\eps\cdot\OPT\norm*{\bfy\vert_{T_j}}_{H,\bfw}.
\]

Union bounding and summing over the $j$ gives us that

\[
\norm*{\bfy}_{H,\bfw'}^2 = \norm*{\bfy}_{H,\bfw}^2 \pm \Theta(\eps\log\norm*{\bfw}_\infty)\OPT^2
\]

with probability at least $1 - \Theta(\delta\log\norm*{\bfw}_\infty)$. Rescaling $\eps$ and $\delta$ by $\Theta(\log\norm*{\bfw}_\infty)$ gives the claim.

This completes the proof of Lemma <ref>.

Given Lemma <ref>, we now prove the active regression guarantees for a single step of the recursion.

Consider the setting of Lemma <ref>. Let $\bar\bfx$ satisfy

\[
\norm*{\bfA\bar\bfx-\bfz}_{H,\bfw'} \leq (1+\eps)\min_{\bfx}\norm*{\bfA\bfx-\bfz}_{H,\bfw'}.
\]

Then, with probability at least $1-\delta$,

\[
\norm*{\bfA\bar\bfx-\bfz}_{H,\bfw} \leq (1+O(\eps/\delta))\min_{\bfx}\norm*{\bfA\bfx-\bfz}_{H,\bfw}.
\]

We follow the proof of Theorem <ref>. By Markov's inequality, we have that

\[
\norm*{\bfz}_{H,\bfw'}^2 = O(1/\delta)\norm*{\bfz}_{H,\bfw}^2
\]

with probability at least $1-\delta$. Let $\bfA\bfx$ have $\norm*{\bfA\bfx}_{H,\bfw} = O(\OPT)$. Let $\bar\bfz$ be as defined in Lemma <ref>.
Lemma <ref>, we then have that

\begin{align*}
&\abs*{\norm*{\bfA\bfx-\bar\bfz}_{H,\bfw}^2-\norm*{\bfA\bfx-\bfz}_{H,\bfw}^2 - \norm*{\bfz-\bar\bfz}_{H,\bfw}^2} = O(\eps)\OPT^2 \\
&\abs*{\norm*{\bfA\bfx-\bar\bfz}_{H,\bfw'}^2-\norm*{\bfA\bfx-\bfz}_{H,\bfw'}^2 - \norm*{\bfz-\bar\bfz}_{H,\bfw'}^2} = O(\eps/\delta)\OPT^2.
\end{align*}

Furthermore, Lemma <ref> guarantees that

\[
\abs*{\norm*{\bfA\bfx-\bar\bfz}_{H,\bfw'}^2 - \norm*{\bfA\bfx-\bar\bfz}_{H,\bfw}^2} = O(\eps)\cdot\OPT^2.
\]

Then by the triangle inequality,

\[
\abs*{\norm*{\bfA\bfx-\bfz}_{H,\bfw'}^2-\norm*{\bfA\bfx-\bfz}_{H,\bfw}^2 + (\norm*{\bfz-\bar\bfz}_{H,\bfw'}^2 - \norm*{\bfz-\bar\bfz}_{H,\bfw}^2)} \leq O(\eps/\delta)\OPT^2.
\]

Let $C\coloneqq \norm*{\bfz-\bar\bfz}_{H,\bfw'}^2 - \norm*{\bfz-\bar\bfz}_{H,\bfw}^2$. Note that $\abs*{C} \leq O(1/\delta)\OPT^2$. Let $\bfx^*\coloneqq \arg\min_\bfx\norm*{\bfA\bfx-\bfz}_{H,\bfw}$. Then,

\begin{align*}
\norm*{\bfA\bar\bfx-\bfz}_{H,\bfw}^2 &\leq \norm*{\bfA\bar\bfx-\bfz}_{H,\bfw'}^2+C + O(\eps/\delta)\OPT^2 \\
&\leq (1+\eps)\norm*{\bfA\bfx^*-\bfz}_{H,\bfw}^2 + C + O(\eps/\delta)\OPT^2 \\
&\leq (1+\eps)(\norm*{\bfA\bfx^*-\bfz}_{H,\bfw}^2 - C) + C + O(\eps/\delta)\OPT^2 \\
&\leq \norm*{\bfA\bfx^*-\bfz}_{H,\bfw}^2+ O(\eps/\delta)\OPT^2.\qedhere
\end{align*}

We can now prove Theorem <ref>.

By an analysis identical to Theorem <ref>, the algorithm makes at most $O(\log\log n)$ recursive calls. We then scale the failure rate $\delta$ and accuracy $\eps$ in <ref> by $\Theta(1/\log\log n)$ so that we can union bound over $O(\log\log n)$ recursive calls, and so that the error incurred at each recursive call is at most $(1+O(\eps/\log\log n))$, with probability at least $99/100$. The total error over the $O(\log\log n)$ recursive calls is then at most

\[
(1+O(\eps/\log\log n))^{O(\log\log n)} \leq 1+O(\eps)
\]

so we conclude as desired.

§ APPLICATIONS: DIMENSION REDUCTION FOR GAMMA FUNCTIONS

The work of [5] introduced the $\gamma_p$ functions, which are also known as smoothed $p$-norms and mixed $\ell_2$-$\ell_p$ norms:

For $t\geq 0$ and $p\geq 1$, define

\[
\gamma_p(t, x) = \begin{cases}
\frac{p}{2} t^{p-2} x^2 & \text{if $\abs{x}\leq t$} \\
\abs*{x}^p + \parens*{\frac{p}{2}-1}t^p & \text{if $\abs{x}>t$}
\end{cases}
\]

The $\gamma_p$ functions have found many applications in fast algorithms for $\ell_p$ regression and related problems [5, 3, 1, 34]. Recall also the simpler $\ell_2$-$\ell_p$ loss defined in Lemma <ref>, which is simply $M_p(x) = \min\{\abs{x}^2,\abs{x}^p\}$. Now note that $\gamma_p(t, x) = \Theta(t^{p}) M_p(x/t)$. Thus, for the purpose of bounding the associated sensitivities up to constant factors, we may replace our discussion of $\gamma_p$ by $M_p$, for simplicity.

We now state our result for dimension reduction for $\gamma_p$ functions, which is analogous to our Huber subspace embedding result in Theorem <ref>.

For $p\geq 1$, thresholds $\bft\geq \mathbf{1}_n$, and weights $\bfw\geq\mathbf{1}_n$, define the weighted $\gamma_p$-norm of $\bfy\in\mathbb R^n$ as

\[
\norm*{\bfy}_{\gamma_p(\bft,\cdot), \bfw} \coloneqq \bracks*{\sum_{i=1}^n \bfw_i \cdot \gamma_p(\bft_i, \bfy_i)}^{1/p}.
\]

If $\bfw = \mathbf{1}_n$. We simply write $\norm*{\bfy}_{\gamma_p(\bft,\cdot)}$. Note that while we refer to this as a norm by abuse of notation, it does not satisfy the properties of a norm.

Let $p\in(1,2)$ and let $\beta$ be as defined in Lemma <ref>. Let $\bft\geq\mathbf{1}_n$. There is an algorithm which, with probability $1-\delta$, computes a set of weights $\bfw\geq\mathbf{1}_n$ with

\[
\nnz(\bfw) \leq O\parens*{d^{1+\beta}\poly\parens*{\log \frac{n}{\delta},\eps^{-1},\log\norm*{\bft}_\infty}}
\]

such that $\norm*{\bfA\bfx}_{\gamma_p(\bft,\cdot), \bfw}^2 = (1\pm \eps) \norm*{\bfA\bfx}_{\gamma_p(\bft,\cdot)}^2$ for all $\bfx\in\mathbb R^d$, and runs in time

\[
O\parens*{\bracks*{\nnz(\bfA) + d^{\omega\cdot(1+1/\beta)}}\poly\parens*{\log n,\log\frac1\delta,\eps^{-1},\log\norm*{\bft}_\infty}}.
\]

Because our proof is a simple extension of our Huber subspace embedding proof (Theorem <ref>), we only outline the proof and necessary changes from the Huber result.

We first bucket the coordinates of $\bft$ into $O(\log\norm*{\bft}_\infty)$ buckets $\braces*{i\in[n] : 2^{j-1} \leq \bft_i < 2^j}$. Note that in typical applications, $\norm*{\bft}_\infty$ is roughly $\poly(n)$ (see, e.g., the paragraph following Lemma 5.5 of [3]). A similar bucketing trick is used in [34]. Once we have done this, the $\gamma_p$ functions are within constant factors of each other, and thus we may focus on sensitivity sampling with a fixed threshold $t$, which we WLOG assume is $1$, and just refer to the corresponding loss as $\gamma_p$. Note that even with a constant factor variation in the $\gamma_p(t_i,\cdot)$ function, we can still get $(1+\eps)$ approximations around the mean by sampling (see, e.g., how we handle weights differing by constant factors in Lemma <ref>).

We first go over one step of the recursive sampling procedure, which reduces the number of rows from roughly $n$ to $n^{\beta/(1+\beta)}d$. Let $\gamma$ satisfy $1/\gamma = (\gamma n)^\beta$, or $\gamma = n^{-\beta/(1+\beta)}$. Now consider the coordinates $i\in[n]$ with $\gamma_p$-sensitivities at least $\gamma$, that is,

\[
S\coloneqq \braces*{i\in[n] : \bfs_i^{\gamma_p}(\bfA) \geq \gamma}
\]

\[
\bfs_i^{\gamma_p}(\bfA) \coloneqq \sup_{\bfx\in\mathbb R^d} \frac{\gamma_p([\bfA\bfx](i))}{\norm*{\bfA\bfx}_{\gamma_p}^p} = \sup_{\bfx\in\mathbb R^d} \frac{\gamma_p([\bfA\bfx](i))}{\sum_{i'=1}^n \gamma_p([\bfA\bfx](i'))} = \sup_{\bfx\in\mathbb R^d} \frac{\gamma_p(1, [\bfA\bfx](i))}{\sum_{i'=1}^n \gamma_p(1, [\bfA\bfx](i'))}.
\]

As done in Lemma <ref>, we may find a superset $S'\supseteq S$ of such coordinates in time at most

\[
O\parens*{\bracks*{\nnz(\bfA) + n^{1/(1+\beta)}d^\omega}\poly\parens*{\log n,\log\frac1\delta}} \leq O\parens*{\bracks*{\nnz(\bfA) + d^{\omega\cdot(1+1/\beta)}}\poly\parens*{\log n,\log\frac1\delta}}
\]

Since the sensitivity upper bounds sum to at most $O(d\poly\log(n/\delta))$, we have that
\[
\abs*{S'} \leq O(\gamma^{-1}d\poly\log(n/\delta)).
\]

Now by Lemma <ref>, any $\bfy\in\mathbb R^n$ either has distortion at most $\gamma$ in the ratio between $\norm*{\bfy}_{\gamma_p}^p$ and $\min_{q\in[p,2]}\norm*{\bfy}_{q}^q$, or the same ratio when restricted to coordinates outside of $S'$. This means that the sum of Lewis weights in the range $[p, 2]$ (appropriately discretized by increments of $1/\log n$) bounds the $\gamma_p$-sensitivities up to a factor of $\gamma$ either on all the coordinates or restricted outside of the set $S'$, which means that the net arguments for Lewis weights via [6] yields a net bound of shows that we can preserve all $\gamma_p$ norms up to $(1\pm\eps)$ factors with $O(\gamma^{-1}d\poly(\log(n/\delta), \eps^{-1})) = O(n^{\beta/(1+\beta)}d\poly(\log(n/\delta), \eps^{-1}))$ samples.

Finally, via a recursive analysis similar to that done in Theorem <ref>, we can repeat until we reach the desired number of rows in $O(\log\log n)$ iterations.

Although the running time as stated is quite large, one can further optimize this in a similar manner as discussed in Section <ref>. For example, we can reduce to $d^C$ dimensions in time $d^{C+1}$, where $C\approx 1.91236$, even for $p$ near $1$, and the trade-offs only improve as $p\to 2$.

§ APPLICATIONS: ACTIVE REGRESSION FOR THE TUKEY LOSS

An important loss function in practice is the Tukey loss:

For any $\tau \ge 0$, the Tukey loss function is given by:

\begin{align*}
T(x) = \begin{cases} \frac{\tau^2}{6} [1 - (1-x ^2/\tau^2)^3] \quad |x| \le \tau \\ \frac{\tau^2}{6} \quad \mathrm{otherwise} \end{cases}
\end{align*}

For simplicity throughout, we assume that $\tau = 1$.
For a vector $\bfy \in \R^n$ we denote $\norm{\bfy}_M = \sum_{i=1}^n M(\bfy(i))$. For a set of non-negative weights $\bfw \in \R^n$, we let $\norm{\bfy}_{M,\bfw} = \sum_{i=1}^n \bfw(i) \cdot M(\bfy(i))$.

Note that our theorem for general $M$-estimators does not capture Tukey loss, due to our assumption of at least polynomial growth, whereas the Tukey loss is flat for $\abs*{x}\geq\tau$. In this section, we show that nonetheless, our active regression techniques apply to this loss as well.

Our results will in fact hold for a more general class of Tukey-like loss functions defined in Assumption 1 of [23]. This class includes the $\ell_p$ Tukey losses, which we focus on in our proofs for simplicity.
For any $\tau \ge 0$, the $\ell_p$ Tukey loss function is given by:

\begin{align*}
M(x) = \begin{cases} |x|^p \quad |x| \le \tau \\ \tau^p \quad \mathrm{otherwise} \end{cases}
\end{align*}

For simplicity throughout, we assume that $\tau = 1$.

We can check that for $\tau = 1$ and $p = 2$, the loss in Definition <ref> is equal to the loss in Definition <ref> up to constants, so identical bounds will hold for it.

§.§ Preliminaries

We next state a few basic results regarding the $\ell_p$ Tukey loss function, shown in [23]. Mostly, these are analogs of the triangle inequality, which must be shown since $\norm{\cdot }_M$ is not a norm.

There is some fixed constant $c$ (depending on $p$) such that, for any $a, b \in \R^n$, and any non-negative weight vector $w \in \R^n$,

$$\norm{\bfa+\bfb}_{M,\bfw} \le c \cdot (\norm{\bfa}_{M,\bfw} + \norm{\bfb}_{M,\bfw}).$$

This simply uses the fact that $M(a) = |a|^p$ when $|a| \le 1$ and $M(a) = 1$ otherwise.

There is some fixed constant $c$ (depending on $p$) such that, for any $a, b \in \R$ with $|a| \le \epsilon |b|$,

$$ (1- c \epsilon) M(b) \le M(a+b) \le (1+ c \epsilon) M(b).$$

From Lemma <ref> we can also prove the following, which is used in Lemma <ref>:

There is some fixed constant $c$ such that, for any $a, b \in \R$ with $M(a) \le \epsilon^p M(b)$,

$$ (1- c \epsilon) M(b) \le M(a+b) \le (1+ c \epsilon) M(b).$$

If $M(a) \le \epsilon^p M(b)$, then since $M(a) = |a|^p$ when $|a|\le 1$ and $M(a) = 1$ otherwise, we must have $|a| = O(\epsilon) |b|$. The bound then follows from Lemma <ref>.

Relatedly we have,

There is some fixed constant $c$ (depending on $p$) such that, for any $\bfa, \bfb \in \R^n$ and non-negative weight vector $\bfw \in \R^n$ with $\norm{\bfa}_{M,\bfw} \le \epsilon^{2p+1} \norm{\bfb}_{M,\bfw}$,

$$ (1- c \epsilon) \norm{\bfb}_{M,\bfw} \le \norm{\bfa+\bfb}_{M,\bfw} \le (1+ c \epsilon) \norm{\bfb}_{M,\bfw} .$$

§.§ Constant Factor Approximation

As with $\ell_p$ regression, we start by using a subspace embedding result for Tukey regression to give a constant factor approximate regression algorithm. Unlike Theorem <ref> for the $\ell_p$ norm, our embedding will only hold for a finite net of vectors rather than the full column span of the sampled matrix. As in [23], we will later handle this assumption by assuming finite bit complexity of our inputs and outputs. See Section <ref>.

For any $\bfA \in \R^{n \times d},\bfb \in \R^n$, there is an algorithm that outputs a weight vector $\bfw \in \R^n$ such that for any set finite set of vectors $\mathcal{N}$ in the span of $[\bfA,\bfb]$, $\bfw$ only has $\tilde O(d^{\max(1,p/2)} \log(|\mathcal{N}|/\delta)/\epsilon^2)$ nonzero entries and with probability $\ge 1-\delta$, for all $\bfy \in \mathcal{N}$,

$$(1-\epsilon) \norm{\bfy}_{M} \le \norm{\bfy}_{M,\bfw} \le (1+\epsilon) \norm{\bfy}_{M}.$$

We observe that the sampling algorithm used to achieve Lemma 6.7 preserves the norm of any vector in expectation. Thus via Markov's inequality we have:

Let $\bfw$ be the weight vector produced by Theorem <ref>. Let $\bfx \in \mathbb{R}^d$ be fixed.
Then with arbitrarily large constant probability, $\norm{\bfA\bfx-\bfb}_{M,\bfw} = O(1) \norm{\bfA\bfx-\bfb}_M$.

We use the notation in Section 6.1 of [23].
Given a weight vector $\bfw$, one step of the recursive sampling algorithm of [23] is to
assign probabilities $p_1, \ldots, p_n$ to each of the $n$ rows of $\bfA$, respectively. For the purposes of this lemma, it does not
matter what these sampling probabilities are.
One then chooses a sampling and rescaling matrix $\bfS$ according to this distribution, where if row $i$ is sampled,
then its new weight $\bfw_{i}'$ is set to $\bfw_i/p_i$, otherwise its weight is set to $0$. Thus, for any fixed $\bfx \in \R^d$ which does not depend on $\bfw$:

$$\E [\|\bfA\bfx-\bfb\|_{M, \bfw'}] = \|\bfA\bfx-\bfb\|_{M, \bfw},$$

When the algorithm begins, $\bfw = \mathbf{1}_n$. It follows inductively that for the final weight vector $\bfw$, that for any fixed $\bfx$, we have:

$$\E[\|\bfA\bfx-\bfb\|_{M, \bfw}] = \|\bfA\bfx-\bfb\|_M.$$

The lemma follows by Markov's inequality.

Combining Theorem <ref> and Lemma <ref>, we have an analog to Theorem <ref> for Tukey regression.

Let $\bfw$ be a weight vector, and
suppose $\norm{\bfA (\bfx_1-\bfx _2)}_{M,\bfw} = (1 \pm \epsilon_0) \norm{\bfA (\bfx_1-\bfx _2)}_{M}$ for
all $\bfx_1,\bfx_2 \in \mathcal{N}$, where $\mathcal{N}$ is a fixed set of vectors. Let
$\tilde \bfx = \argmin_{\bfx \in \mathcal{N}} \norm{\bfA\bfx-\bfb}_{M,\bfw}$, and let
$\bfx^* = \argmin_{\bfx \in \mathcal{N}} \norm{\bfA\bfx-\bfb}_{M}$.
Suppose also that $\norm{\bfA\bfx^*-\bfb}_{M,\bfw} = O(1) \norm{\bfA\bfx^*-\bfb}_M$.
Then $\norm{\bfA \tilde \bfx-\bfb}_M = O(1) \norm{\bfA\bfx^*-\bfb}_M$.

By Theorem <ref> and Lemma <ref> there is an algorithm producing $\bfw$ satisfying the required subspace embedding and norm preservation bounds with just $\tilde O(d^{\max(1,p/2)} \log(|\mathcal{N}|/\delta)/\epsilon^2)$ non-zero entries. Thus, computing $\tilde \bfx$ requires reading just this number of entries from $\bfb$.

We have

\begin{align*}
\norm{\bfA \tilde \bfx-\bfb}_M &= O(1) \cdot (\norm{\bfA \tilde \bfx-\bfA\bfx^*}_M + \norm{\bfA\bfx^*-\bfb}_M)\tag{Approximate triangle inequality, Lem. \ref{lem:tukTri1}} \\
&\leq O(1) \cdot \left (\frac{1}{1-\epsilon_0} \cdot \norm{\bfA (\tilde \bfx-\bfx ^*)}_{M,\bfw} + \norm{\bfA\bfx^*-\bfb}_M \right )\tag{Subspace embedding assumption, since $\tilde \bfx,\bfx^* \in \mathcal{N}$}\\
&\leq O(1) \cdot \left (\frac{1}{1-\epsilon_0} \cdot \left (\norm{\bfA \tilde \bfx-\bfb}_{M,\bfw} + \norm{\bfA\bfx^*-\bfb}_{M,\bfw} \right ) + \norm{\bfA\bfx^*-\bfb}_M \right )\tag{Approximate triangle inequality, Lem. \ref{lem:tukTri1}} \\
&\leq O(1) \cdot \left (\frac{2}{1-\epsilon_0} \norm{\bfA\bfx^*-\bfb}_{M,\bfw} + \norm{\bfA\bfx^*-\bfb}_M \right )\tag{Optimality of $\tilde \bfx$ for the sampled problem}\\
& \leq O(1) \cdot \norm{\bfA\bfx^*-\bfb}_M\tag{Assumption that $\norm{\bfA\bfx^*-\bfb}_{M,\bfw} = O(1) \norm{\bfA\bfx^*-\bfb}_M$}.
\end{align*}

§.§ Net Argument

The constant factor approximation result of Theorem <ref> allows us to approximately minimize the Tukey loss regression problem over any finite net of vectors $\mathcal{N}$.
As in [23],
we need to make an assumption for the net $\mathcal{N}$ to have a small size. For simplicity, we adopt the
assumption stated after Assumption 2 in [23], that the entries of $\bfA$ and $\bfb$ are
integers of absolute values at most $n^{\poly(d)}$.
We can also assume
$\bfA$ has full rank, since removing linearly dependent columns can be done in a preprocessing stage.

Assuming that $\bfA,\bfb$ have integer entries bounded by $n^{\poly(d)}$ and $\bfA$ has full column rank, if $\bfx \in \R^d$ is such that $\norm{\bfA\bfx-\bfb}_M \le C \cdot \min_{\bfz \in \R^d} \norm{\bfA \bfz-\bfb}_M$ for some $C \leq \poly(n)$,
and we round the entries of $\bfx$ to the nearest integer multiple of $1/n^{\poly(d)}$, obtaining $\bfx'$,
then $\norm{\bfA\bfx' -\bfb}_M \le C \cdot (1+1/n^{\poly(d)}) \cdot \min_{\bfz \in \R^d} \norm{\bfA \bfz-\bfb}_M$. Further, all entries
of $\bfx$ are at most $n^{d^{\alpha}}$ for a fixed constant $\alpha > 0$.

Further, if $\norm{\bfA\bfx-\bfb}_M = 0$, then necessarily $\bfx' =\bfx$. Also, if $\norm{\bfA \bfz-\bfb}_M > 0$
for all $\bfz$, then $\min_\bfz \norm{\bfA \bfz-\bfb}_M \geq 1/n^{\poly(d)}$.

First suppose that $\min_\bfz \norm{\bfA \bfz-\bfb}_M = 0$. This implies that $\bfA\bfx = \bfb$ and thus $\bfx = (\bfA^\top \bfA)^{-1} \bfA^\top \bfb$. By Cramer's rule
and our integrality assumptions, it follows that the entries of $\bfx$ are integer multiples of $1/n^{\poly(d)}$
and bounded in absolute value by $n^{\poly(d)}$. Thus, $\bfx$ is rounded to itself, showing that $\bfx' =\bfx$
and the lemma holds.

Otherwise, suppose that $\min_\bfz \norm{\bfA\bfx-\bfb}_M > 0$. Then the claim
is that $\min_\bfz \norm{\bfA \bfz-\bfb}_M \ge 1/n^{\poly(d)}$. To see this,
note that the matrix $[\bfA, \bfb]$ has rank-$(d+1)$, since the columns of $\bfA$ are linearly independent.
Since $[\bfA, \bfb]$ is an integer matrix, it follows from
Lemma 4.1 of [20]
that its minimum singular value is at least $n^{-\textrm{poly}(d)}$.
Notice that $\min_\bfz \|\bfA\bfz-\bfb\|_2^2$ is at least the squared minimum singular value of $[\bfA,\bfb]$, since:

\begin{align*}
\min_{\bfz \in \R^d} \|\bfA\bfz-\bfb\|_2^2 = \min_{\bfz' \in \R^{d+1}: \bfz'(d+1) = -1} \norm{[\bfA,\bfb] \bfz'}_2^2 \ge \min_{\bfz' \in \R^{d+1}: \norm{\bfz'} \ge 1} \norm{[\bfA,\bfb] \bfz'}_2^2.
\end{align*}

The last equation is exactly the minimum squared singular value of $[\bfA,\bfb]$.

Consequently, the $\ell_2$-regression
cost is lower bounded by $n^{-\textrm{poly}(d)}$, which implies the $\ell_p$-regression cost
is lower bounded by $n^{-\textrm{poly}(d)}$, and the Tukey loss is lower bounded by
$\min(1, n^{-\textrm{poly}(d)}) = n^{-\textrm{poly}(d)}$ (since we assume $\tau = 1$). In this
case, if we round $\bfx$ to the nearest vector with entries that are integer multiples of $1/n^{\poly(d)}$,
then using that the entries of $\bfA$ are bounded by $n^{\poly(d)}$, we still obtain a
$C \cdot (1+1/n^{\poly(d)})$-approximate solution.

Finally, note that the entries in any $\poly(n)$-approximate solution $\bfx$ cannot be larger than
$n^{\poly(d)}$. Indeed, by Lemma 5.2 of [23] and since parts 2 and 3 of
Assumption 2 of [23] hold by our integrality assumption, it follows that part 1
of Assumption 2 of [23] also holds, and so $\|\bfx\|_2 \leq n^{d^{\alpha}}$ for a fixed
constant $\alpha > 0$.

Consider the set $\mathcal{N}$ of all vectors $\bfx \in \mathbb{R}^d$ that are integer multiples of
$1/n^{\poly(d)}$ and bounded in absolute value by $n^{\poly(d)}$, for a sufficiently large $\poly(d)$.
Then $|\mathcal{N}| \leq n^{\poly(d)}$ and we have

Suppose $\norm{\bfA\bfx'-\bfb}_{M,\bfw} = (1 \pm \epsilon)\norm{\bfA\bfx'-\bfb}_M$ for all $\bfx' \in \mathcal{N}$
for a weight vector $\bfw$ with $\|\bfw\|_{\infty} \leq n^{\poly(d)}$.
Then $\norm{\bfA\bfx-\bfb}_{M,\bfw} = (1 \pm \epsilon)\norm{\bfA\bfx-\bfb}_M$ for all $\bfx$ for which
$\|\bfx\|_2 \leq n^{d^\alpha}$ for a fixed constant $\alpha > 0$. Further,
any $C$-approximate minimizer to $\norm{\bfA\bfx-\bfb}_M$, for some $C \leq \poly(n)$,
necessarily satisfies $\|\bfx\|_2 \leq n^{d^{\alpha}}$.

Let $\bfx'$ be the rounding in $\mathcal{N}$ of such an $\bfx$, as in Lemma <ref>.
If $\min_\bfx \norm{\bfA\bfx-\bfb}_M = 0$, then necessarily $\bfx' =\bfx$ by Lemma <ref>.
Note that since the columns of $\bfA$ are linearly independent, $\bfx$ is unique, and as in
the proof of Lemma <ref>, we have $\|\bfx\|_2 \leq n^{d^{\alpha}}$ for a fixed
constant $\alpha > 0$.

Otherwise, we can observe that $\|\bfw\|_{\infty} \cdot \|\bfA(\bfx-\bfx ')\|_{\infty} = 1/n^{\poly(d)}$ by our assumption that $\norm{\bfw}_\infty \le n^{\poly(d)}$ and the fact that that the granularity of our net is $1/n^{\poly(d)}$ for a sufficiently large $\poly(d)$. This gives $\norm{\bfA (\bfx-\bfx ')}_{M,\bfw} \le 1/n^{\poly(d)}$. Further, $\norm{\bfA\bfx'-\bfb}_M \ge 1/n^{\poly(d)}$ and $\norm{\bfA\bfx-\bfb}_M \ge 1/n^{\poly(d)}$ by Lemma <ref>. In turn, $\norm{\bfA\bfx'-\bfb}_{M,\bfw} \ge 1/n^{\poly(d)}$ by our subspace embedding assumption on $\bfw$. Combining these bounds with the approximate triangle inequality of Lemma <ref> gives:

\begin{align*}
\norm{\bfA\bfx-\bfb}_{M,\bfw} & = \norm{\bfA\bfx'-\bfb + \bfA(\bfx-\bfx ')}_{M,\bfw}\\
&= (1 \pm 1/n^{\poly(d)}) \norm{\bfA\bfx'-\bfb}_{M,\bfw}\tag{Lemma \ref{lem:tukTri4}}\\
%& = & (1 \pm \epsilon) \norm{\bfA\bfx'-\bfb}_M \pm 1/n^{\poly(d)}\\
&= (1 \pm O(\epsilon)) \norm{\bfA\bfx'-\bfb}_M\tag{Subspace embedding assumption on $\bfw$ for $\bfx' \in \mathcal{N}$}\\
&= (1 \pm O(\epsilon)) \norm{\bfA\bfx-\bfb  + \bfA(\bfx'-\bfx )}_M\\
&= (1 \pm O(\epsilon)) \norm{\bfA\bfx-\bfb}_M\tag{Lemma \ref{lem:tukTri4} again}.
\end{align*}

Theorem <ref> implies that we can use Theorem <ref> with our net $\mathcal{N}$,
so that if we obtain a constant factor solution to the sampled regression problem $\min_\bfx \norm{\bfA\bfx-\bfb}_{M,\bfw}$ over $\bfx$ with $\|\bfx\|_2 \leq n^{d^{\alpha}}$, then we
obtain an initial constant factor solution to the original regression problem.

§.§ Relative Error Approximation

We now build on our constant factor approximation result and the net arguments above to give a relative error approximation result for Tukey regression, akin to Theorem <ref>. We assume that $\norm{\bfb}_M = O(\OPT)$, which we can ensure by computing an initial constant factor approximation as in Algorithm <ref>, using Theorem <ref> combined with Theorem <ref>.

Consider $\bfA \in \R^{n \times d}$, $\bfb \in \R^n$, and $p \ge 1$. Let $\OPT = \min_{\bfx \in \R^d} \norm{\bfA\bfx - \bfb}_M$ where $M$ is the $\ell_p$ Tukey loss (Def. <ref>) and assume that $\norm{\bfb}_M = O(\OPT)$. There is an algorithm that reads $\tilde O \left (\frac{d^{\max(1,p/2) + O(1)}}{\epsilon^{2+p}} \poly\log(n)\right )$ entries of $\bfb$ and outputs $\tilde \bfx$ satisfying $\norm{\bfA \tilde \bfx - \bfb}_M \le (1+\epsilon) \OPT$ with probability $99/100$.

Let $\tilde\bfs_i^{M}(\bfA)$ be as computed by the $\ell_p$ Tukey sensitivities as computed by Algorithm <ref> and Theorem <ref>. Note then that

\[
\sum_{i=1}^n \tilde\bfs_i^{M}(\bfA) \leq O(d^{\max\{1,p/2\}}\log^2 n).
\]

Note that Lemmas <ref>, <ref>, and <ref> all hold for the Tukey loss, where assumption (1) is satisfied in Lemma <ref> due to Lemma <ref>. Thus, as reasoned in Lemma <ref>, we just need to choose an $m$ that is sufficient to union bound over a net. For the Tukey loss, we set

\[
m = \frac{\poly(d) \log n}{\epsilon^{2+p}}
\]

and apply a union bound over all vectors $\bfA\bfx$ in a net $\mathcal{N}$ with $|\mathcal{N}| \le n^{\poly(d)}$ gives that this bound holds with high probability simultaneously for all vectors in the net.
Letting this net be as in Theorem <ref>, this is enough to give that the bound holds for all $\bfx$ of interest.

Note that the number of entries of $\bfb$ that are read is, in expectation,

\[
\sum_{i=1}^n \bfp_i \leq \sum_{i=1}^n m\cdot \tilde\bfs_i^M(\bfA) = O\parens*{d^{\max\{1+p/2\}+O(1)} \frac{\log^3 n}{\eps^{2+p}}}.
\]

Thus, this is, up to constant factors, the number of entries read with probability at least $99/100$ as well.

§ APPLICATIONS: KRONECKER PRODUCT REGRESSION

In the $q$-th order Kronecker product regression problem, one is
given matrices $\bfA_1, \bfA_2,\ldots, \bfA_q$, where $\bfA_i \in \mathbb{R}^{n_i \times d_i}$,
as well as a vector $\bfb \in \mathbb{R}^{n_1 n_2 \cdots n_q}$, and
the goal is to obtain a solution to the problem:

\[\min_{\bfx \in \R^{d_1 d_2 \cdots d_q}} \|( \bfA_1 \otimes \bfA_2 \cdots \otimes \bfA_q)x - \bfb\|_p, \]

where $\bfA_1 \otimes \bfA_2 \cdots \otimes \bfA_q$ is the
$\mathbb{R}^{(n_1 \times n_2 \times \cdots \times n_q)\bfx (d_1 \times d_2 \times \cdot \times d_q)}$ matrix
whose $((i_1, i_2, \ldots, i_q), (j_1, \ldots, j_q))$-th entry is equal to
$(\bfA_1)_{i_1, j_1} \cdots (\bfA_q)_{i_q, j_q}$. Typically $\prod_{i=1}^q d_i$ is much less
than $\prod_{i=1}^q n_i$, and the goal is to obtain algorithms that do not explicitly form
$\bfA_1 \otimes \bfA_2 \otimes \cdots \otimes \bfA_q$, which is too expensive.

Kronecker product regression is a special case of $\ell_p$-regression
when the design matrix is highly structured. Such matrices
arise in spline regression,
signal processing, and multivariate data fitting
[68, 69, 36].
Recent work [26, 25] designed algorithms for $p \in [1,2]$ and constant $q$
which solve this
problem with constant probability up to a $(1 \pm \epsilon)$-factor in
$\tilde{O}(\sum_{i=1}^q \nnz(\bfA_i) + \nnz(\bfb) + \poly(d_i/\epsilon))$ time;
see Theorem 1.2 of [25]. While $\sum_{i=1}^q \nnz(\bfA_i)$ is bounded by
$O(\sum_i n_i d_i)$, the main drawback of these works is that $\nnz(\bfb)$ may
be as large as $n_1 \times n_2 \times \cdots \times n_q$, which is prohibitive.
Theorem <ref> shows that for every real $p \geq 1$, only $\tilde O \left (\prod_{i=1}^q d_i^{\max(1,p/2)} \right )$
samples are needed, and our runtime is polynomial in this quantity. This should
be contrasted with prior work which required $\nnz(\bfb)$ time for every $p \neq 2$,
where $\nnz(\bfb)$ can be as large as $\prod_{i=1}^q n_i$.

Before stating our theorem, we need to understand Lewis weights of the Kronecker product
of two matrices $\bfA$ and $\bfB$. Recall from Definition <ref> that the Lewis weights of an $n \times d$ matrix $\bfA$ are
the unique numbers $\bfw^p_1(\bfA), \ldots, \bfw^p_n(\bfA)$ which satisfy
$\bfw_i^p(\bfA) = \bfs_i^2 (\bfW^{1/2-1/p} \bfA)$, where $\bfs_i^2(\bfW^{1/2-1/p} \bfA)$ is the $i^{th}$ leverage score of $\bfW^{1/2-1/p} \bfA$, and $\bfW$ is a diagonal matrix with $i^{th}$
diagonal entry equal to $\bfw_i^p(\bfA)$. We can show that the Lewis weights of $\bfA \otimes \bfB$ are just the products of those of $\bfA$ and $\bfB$.

Given matrices $\bfA,\bfB$, the $(i,j)^{th}$ Lewis weight of $\bfA \otimes \bfB$ is given by

$$\bfw_{i,j}^p(\bfA \otimes \bfB) = \bfw_i^p(\bfA) \cdot \bfw_j^p(\bfB).$$

let $\bfW_\bfA$ and $\bfW_\bfB$ denote the diagonal matrices with the Lewis weights of $\bfA$ and $\bfB$ respectively on their diagonals. To prove the lemma, we just need to verify that $\bfw_i^p(\bfA) \cdot \bfw_j^p(\bfA) = \bfs_{i,j}^2 ((\bfW_\bfA^{1/2 - 1/p} \bfA) \otimes (\bfW_\bfB^{1/2-1/p} \bfB))$.
By Proposition 3.2 of [25], the $(i,j)$-th row leverage score of the
Kronecker product of two matrices $\bfC$ and $\bfD$ is $\bfs_i^2(\bfC) \cdot \bfs_j^2(\bfD)$. For completeness: to see this,
recall the leverage scores are basis-invariant and thus a property of the column span of $\bfC \otimes \bfD$.
If $\bfR$ and $\bfS$ are square matrices for which $\bfC\bfR$ has orthonormal columns and $\bfD\bfS$ has orthonormal columns, then
one can verify that $(\bfC \otimes \bfD) \cdot (\bfR \otimes \bfS)$ has orthonormal columns, which implies the statement.

$$\bfs_{i,j}^2 ((\bfW_\bfA^{1/2 - 1/p} \bfA) \otimes (\bfW_\bfB^{1/2-1/p} \bfB)) = \bfs_i^2(\bfW_\bfA^{1/2 - 1/p} \bfA) \cdot \bfs_j^2(\bfW_\bfB^{1/2 - 1/p} \bfB) = \bfw_i^p(\bfA) \cdot \bfw_i^p(\bfB).$$

Given matrices $\bfA_1, \ldots \bfA_q$,
the $(i_1, \ldots, i_q)$-th Lewis weight of $\bfA_1 \otimes \cdots \otimes \bfA_q$ is

$$\bfw_{i_1, \ldots, i_q}^p(\bfA_1 \otimes \cdots \otimes \bfA_q) = \bfw^p_{i_1}(\bfA_1) \cdot \bfw^p_{i_2}(\bfA_2) \cdot \ldots \cdot \bfw^p_{i_q}(\bfA_q).$$

Follows from Lemma <ref> and induction.

The following theorem is notable, in that it is the first algorithm for Kronecker product regression,
for every $p \geq 1$, whose running time does not depend on $\nnz(\bfb)$.

Let $q \geq 1$ and $p \geq 1$ be constant, and $\epsilon > 0$. The Kronecker product regression
problem can be solved up to a $(1+\epsilon)$-factor with constant probability in
$\tilde{O}(\sum_{i=1}^q \nnz(\bfA_i) + \poly(\prod_{i=1}^q d_i / \epsilon))$ time.

By Theorem 6.1 of [17], in $\tilde{O}(\sum_{i=1}^q \nnz(\bfA_i))$ time
we can compute $2$-approximations $\hat{u^{\ell}_i}$ to all of the Lewis weights
of all of the input matrices $\bfA^{\ell}$, where $\{u^{\ell}_i\}_i$ is the set of Lewis weights of $\bfA^{\ell}$. By Corollary
<ref>, for each $(i_1, \ldots, i_q)$, we have that $\hat{u^{1}_{i_1}} \cdots \hat{u^q_{i_q}}$ is a $2^q = O(1)$-approximation
to the $(i_1, \ldots, i_q)$-th Lewis weight of $\bfA^1 \otimes \cdots \otimes \bfA^q$. By Theorem <ref>, it follows
that for constant $p$ if we sample $\poly(\prod_{i=1}^q d_i / \epsilon)$ rows of $\bfA^1 \otimes \cdot \otimes \bfA^q$ and
corresponding entries of $\bfb$, then we can solve for a $(1+\epsilon)$-approximate solution $\bfx$ to the $\ell_p$-Kronecker
product regression problem in $\poly(\prod_{i=1}^q d_i / \epsilon)$ time. To sample proportional to the
$\hat{u^{1}_{i_1}} \cdots \hat{u^q_{i_q}}$ values, we can first sample an $i_1$ proportional to $\hat{u^1_{i_1}}$, then independently
sample an $i_2$ proportional to $\hat{u^2_{i_2}}$, and so on. It is possible to build a data structure in $O(n)$ time, one for
each $\bfA^{\ell}$, so that a single sample can be extracted in $O(1)$ time [70, 71].
So the total time for the sampling is
$O(n) + \poly(\prod_{i=1}^q d_i / \epsilon)$, and we can assume $\tilde{O}(\sum_{i=1}^q \nnz(\bfA_i))$ is at least $n$ (otherwise we
can throw away zero rows), and so the overall time is $\tilde{O}(\sum_{i=1}^q \nnz(\bfA_i) + \poly(\prod_{i=1}^q d_i / \epsilon))$.

§ APPLICATIONS: SUBSPACE EMBEDDINGS FOR ORLICZ NORMS

In this section, we show that our sensitivity sampling techniques immediately carry over to Orlicz norms, which can be thought of as scale-invariant versions of $M$-estimators. Subspace embeddings, regression, and row sampling algorithms for Orlicz norms have recently been studied by [4] and [61], respectively.

The row sampling algorithms for Orlicz norm regression in [61] proceeds by constructing a generalization of well-conditioned bases for Orlicz norms using an oblivious subspace embedding, and then sampling the rows proportionally to an importance distribution defined using the well-conditioned basis. Due to this indirect construction of importance distributions, the algorithm of [61] incurs large polynomial factors in $d$ on the sample size (i.e., the number of rows in the subspace embedding). Specifically, their total sensitivity bound from Lemma 4 is $O(d\kappa^2)$ using a well-conditioned basis of distortion $\kappa$, where $\kappa$ given by their Corollary 12 is $O(d^{5/2}\poly\log n)$. Then, the subspace embedding they obtain in Lemma 7 has $O(d^7\poly(\log n,\eps^{-1}))$ rows with constant probability. Furthermore, small oblivious subspace embeddings do not exist when the $M$ estimator function $M$ grows faster than quadratic [8], and this constraint carries over to their construction as well (see Assumption 1 of [61]).

We overcome both of these problems by directly constructing an importance sampling distribution according to the Orlicz sensitivities (see Definition <ref>) using our Theorem <ref>, which removes the quadratic growth assumption, allowing for a polynomial growth bound of degree $p_G$, and significantly improves the dependence on $d$ in the dimension of the subspace embedding. More specificially, in Theorem <ref>, we obtain a $\tilde O(\nnz(\bfA) + \poly(d/\eps))$ time algorithm that computes a weighted subset of $O(d^2\poly(\log n, \eps^{-1}))$ rows that achieves a $(1\pm\eps)$ distortion subspace embedding for $p_G\leq 2$, as well as $O(d^{p_G/2+1}\poly(\log n, \eps^{-1}))$ rows for $p_G > 2$.

Let $G$ satisfy the hypotheses of Corollary <ref> and let $p_G$ be the polynomial upper bound exponent for $G$. There is an algorithm which, with constant probability, computes a set of weights $\bfw'$ such that

\[
\norm*{\bfA\bfx}_{G,\bfw'} = (1\pm\eps)\norm*{\bfA\bfx}_G
\]

for all $\bfx\in\mathbb R^d$ in time

\[
O\parens*{\nnz(\bfA)\log^3 n + \frac{(Td^{\max\{1,p_G/2\}} + T^2)\log^4 n}{\eps^2}\log\frac1\eps},
\]

where constant factor $\ell_{p_G}$ Lewis weight approximation for an $n\times d$ matrix $\bfB$ takes time $O(\nnz(\bfb) + T))$ time (see Theorem <ref>). Furthermore,

\[
\nnz(\bfw') =O\parens*{\frac{d^{\max\{2,p_G/2+1\}} \log^3 n}{\eps^2}\log\frac1\eps}
\]

with constant probability.

§.§ Fast Approximation of Orlicz Sensitivities

First recall the definition of Orlicz norms.