Sketching for Kronecker Product Regression and P-splines

Huaian Diao<EMAIL_ADDRESS>Zhao Song<EMAIL_ADDRESS>Wen Sun
<EMAIL_ADDRESS>David P. Woodruff<EMAIL_ADDRESS>
Northeast Normal University Harvard U & UT-Austin Carnegie Mellon University
Carnegie Mellon University

###### Abstract

TensorSketch is an oblivious linear sketch introduced in (Pagh, 2013) and
later used in (Pham and Pagh, 2013) in the context of SVMs for polynomial
kernels. It was shown in (Avron et al., 2014) that TensorSketch provides a
subspace embedding, and therefore can be used for canonical correlation
analysis, low rank approximation, and principal component regression for the
polynomial kernel. We take TensorSketch outside of the context of polynomials
kernels, and show its utility in applications in which the underlying design
matrix is a Kronecker product of smaller matrices. This allows us to solve
Kronecker product regression and non-negative Kronecker product regression, as
well as regularized spline regression. Our main technical result is then in
extending TensorSketch to other norms. That is, TensorSketch only provides
input sparsity time for Kronecker product regression with respect to the
$2$-norm. We show how to solve Kronecker product regression with respect to
the $1$-norm in time sublinear in the time required for computing the
Kronecker product, as well as for more general $p$-norms.

##  1 INTRODUCTION

In the overconstrained least squares regression problem, we are given an
$n\times d$ matrix $A$ called the “design matrix”, $n\gg d$, and an $n\times
1$ vector $b$, and the goal is to find an $x$ which minimizes $\|Ax-b\|_{2}$.
There are many variants to this problem, such as regularized versions of the
problem in which one seeks an $x$ so as to minimize
$\|Ax-b\|_{2}^{2}+\|Lx\|_{2}^{2}$ for a matrix $L$, or regression problems
which seek to minimize more robust loss functions, such as
$\ell_{1}$-regression $\|Ax-b\|_{1}$.

In the era of big data, large scale matrix computations have attracted
considerable interest. To obtain reasonable computational and time
complexities for large scale matrix computations, a number of randomized
approximation algorithms have been developed. For example, in (Woodruff,
2014), it was shown how to output a vector $x^{\prime}\in\mathbb{R}^{d}$ for
which
$\|Ax^{\prime}-b\|_{2}\leq(1+\epsilon)\min_{x\in\mathbb{R}^{d}}\|Ax-b\|_{2}$
in $\operatorname{\mathtt{nnz}}(A)+\operatorname{poly}(d/\epsilon)$ time,
where $\operatorname{\mathtt{nnz}}(A)$ denotes the number of non-zero entries
of matrix $A$. We refer the reader to the recent surveys (Kannan and Vempala,
2009; Mahoney, 2011; Woodruff, 2014) for a detailed treatment of this topic.

In this work we focus on regression problems for which the design matrix is a
Kronecker product matrix, that is, it has the form $A\otimes B$ for
$A\in\mathbb{R}^{n_{1}\times d_{1}}$ and $B\in\mathbb{R}^{n_{2}\times d_{2}}$.
Also, $b\in\mathbb{R}^{n_{1}n_{2}}$, and one seeks to solve the problem
$\min_{x\in\mathbb{R}^{d_{1}d_{2}}}\|(A\otimes B)x-b\|_{2}$ in the standard
setting, which can also be generalized to regularized and robust variants. One
can also ask the question when the design matrix is a Kronecker product of
more than two matrices.

Kronecker product matrices have many applications in statistics, linear system
theory, signal processing, photogrammetry, multivariate data fitting, etc.;
see (Golub and Van Loan, 2013; Van Loan and Pitsianis, 1993; Van Loan, 1992).
The linear least squares problem involving the Kronecker product arises in
many applications, such as structured linear total least norm on blind
deconvolution problems (Oh and Yun, 2005), constrained linear least squares
problems with a Kronecker product structure, the bivariate problem of surface
fitting and multidimensional density smoothing (Eilers and Marx, 2006).

One way to solve Kronecker product regression is to form the matrix
$C=A\otimes B$ explicitly, where $A\in\mathbb{R}^{n_{1}\times
d_{1}},B\in\mathbb{R}^{n_{2}\times d_{2}}$, and then apply a randomized
algorithm to $C$. However, this takes at least
$\operatorname{\mathtt{nnz}}(A)\cdot\operatorname{\mathtt{nnz}}(B)$ time and
space. It is natural to ask if it is possible to solve the Kronecker product
regression problem in time faster than computing $A\otimes B$. This is in fact
the case, as Fausett and Fulton (Fausett and Fulton, 1994) show that one can
solve Kronecker product regression in $O(n_{1}d_{1}^{2}+n_{2}d_{2}^{2})$ time;
indeed, the solution vector
$x=\textrm{vec}((B^{\perp})^{\top}D^{-1}A^{\perp})$, where $D=\textrm{vec}(b)$
and vec$(E)$ for a matrix $E$ denotes the operation of stacking the columns of
$E$ into a single long vector. While such a computation does not involve
computing $A\otimes B$, it is more expensive than what one would like.

A natural question is if one can approximately solve Kronecker product
regression in
$\operatorname{\mathtt{nnz}}(A)+\operatorname{\mathtt{nnz}}(B)+\operatorname{poly}(d/\epsilon)$
time, which, up to the $\operatorname{poly}(d/\epsilon)$ term, would match the
description size of the input. Another natural question is Kronecker product
regression with regularization. Such regression problems arise frequently in
the context of splines (Eilers and Marx, 2006). Finally, what about Kronecker
product regression with other, more robust, norms such as the $\ell_{1}$-norm?

###  1.1 Our Contributions

We first observe that the random linear map TensorSketch, introduced in the
context of problems for the polynomial kernel by Pagh (2013), Pham and Pagh
(2013), is exactly suited for this task. Namely, in (Avron et al., 2014) it
was shown that for a $d$-dimensional subspace $C$ of $\mathbb{R}^{n}$,
represented as an $n\times d$ matrix, there is a distribution on linear maps
$S$ with $O(d^{2}/\epsilon^{2})$ rows such that with constant probability,
simultaneously for all $x\in\mathbb{R}^{d}$,
$\|SCx\|_{2}=(1\pm\epsilon)\|Cx\|_{2}$. That is, $S$ provides an Oblivious
Subspace Embedding (OSE) for the column span of $C$. Further, it is known that
if $b$ is an $n$-dimensional vector, then one has that
$\|S[C,b]x\|_{2}=(1\pm\epsilon)\|[C,b]x\|_{2}$ for all $x\in\mathbb{R}^{d+1}$.
Consequently, to solve the regression problem
$\min_{x\in\mathbb{R}^{d}}\|Cx-b\|_{2}$, one can instead solve the much
smaller problem $\min_{x\in\mathbb{R}^{d}}\|SCx-Sb\|_{2}$, and the solution
$x^{\prime}$ to the latter problem will satisfy
$\|Cx^{\prime}-b\|_{2}\leq(1+\epsilon)\min_{x\in\mathbb{R}^{d}}\|Cx-b\|_{2}$.

Importantly, if $n=n_{1}\cdot n_{1}$ and there is a basis for the column span
of $C$ of the form $A_{1}\otimes A_{1},A_{2}\otimes A_{2},\ldots,A_{d}\otimes
A_{d}$, then given $A_{1},\ldots,A_{d}$, it holds that $SC$ can be computed in
$\operatorname{\mathtt{nnz}}(A)$ time, where $A$ is the $n_{1}\times d$ matrix
whose $i$-th column is $A_{i}$. Further, given a vector $b$ with
$\operatorname{\mathtt{nnz}}(b)$ non-zero entries, one can compute $Sb$ in
$\operatorname{\mathtt{nnz}}(b)$ time. Thus, one obtains a
$(1+\epsilon)$-approximate solution to the regression
$\min_{x\in\mathbb{R}^{d}}\|Cx-b\|_{2}$ in
$\operatorname{\mathtt{nnz}}(A)+\operatorname{\mathtt{nnz}}(b)+\operatorname{poly}(d/\epsilon)$
time.

While not immediately useful for our problem, we show that via simple
modifications, the claim about TensorSketch above can be generalized to the
case when there is a basis for the column span of $C$ of the form
$A_{i}\otimes B_{j}$ for arbitrary vectors
$A_{1},\ldots,A_{d_{1}}\in\mathbb{R}^{n_{1}}$ and vectors
$B_{1},\ldots,B_{d_{2}}\in\mathbb{R}^{n_{2}}$. That is, we observe that in
this case $SC$ can be computed in
$\operatorname{\mathtt{nnz}}(A)+\operatorname{\mathtt{nnz}}(B)$ time, where
$A$ is the $n_{1}\times d_{1}$ matrix whose $i$-th column is $A_{i}$, and $B$
is the $n_{2}\times d_{2}$ matrix whose $i$-th column is $B_{i}$. In this case
$C=A\otimes B$, which is exactly the case of Kronecker product regression.
Using the above connection to regression, we obtain an algorithm for Kronecker
product regression in
$\operatorname{\mathtt{nnz}}(A)+\operatorname{\mathtt{nnz}}(B)+\operatorname{\mathtt{nnz}}(b)+\operatorname{poly}(d_{1}d_{2}/\epsilon)$
time. Using the fact that $\|SCx-Sb\|_{2}=(1\pm\epsilon)\|Cx-b\|_{2}$ for all
$x$, we have in particular that this holds for all non-negative $x$, and so
also obtain the same reduction in problem size for non-negative Kronecker
product regression, which occurs often in image and text data; see e.g., (Chen
and Plemmons, 2010).

The above observation allows us to extend many existing variants of least
squares regression to the case when $C$ is a Kronecker product of two
matrices. For example, the results in (Avron et al., 2016) for the ridge
regression problem
$\min_{x\in\mathbb{R}^{d}}\|Ax-b\|_{2}^{2}+\lambda\|x\|_{2}^{2}$ immediately
hold when $C$ is a Kronecker product matrix, since the conditions needed for
the oblivious embedding in (Avron et al., 2016) come down to a subspace
embedding and an approximate matrix product condition, both of which are known
to hold for TensorSketch (Avron et al., 2014). More interestingly we are able
to extend the results in (Avron et al., 2016) for Kronecker ridge regression
to the case when the regularizer is a general matrix, namely, to the problem
$\min_{x\in\mathbb{R}^{d}}\|Ax-b\|_{2}^{2}+\lambda\|Lx\|_{2}^{2}$, where $L$
is an arbitrary matrix. Such problems occur in the context of spline
regression (Eilers and Marx, 1996, 2006; Eilers et al., 2015). The number of
rows of TensorSketch depends on a generalized notion of statistical dimension
depending on the generalized singular values $\gamma_{i}$ of $[A;L]$, and may
be much smaller than $d_{1}$ or $d_{2}$. In (Avron et al., 2016), only results
for $L$ equal to the identity were obtained.

Finally, our main technical result is to extend our results to least absolute
deviation Kronecker product regression
$\min_{x\in\mathbb{R}^{d}}\|Cx-b\|_{1}$, which, since it involves the
$1$-norm, is often considered more robust than least squares regression
(Rousseeuw and Leroy, 2005) and has been widely used in applications such as
computer vision (Zheng et al., 2012). We in fact extend this to general
$p$-norms but focus the discussion on $p=1$. We give the first algorithms for
this problem that are faster than computing $C=A\otimes B$ explicitly, which
would take at least
$\operatorname{\mathtt{nnz}}(A)\cdot\operatorname{\mathtt{nnz}}(B)\geq
n_{1}n_{2}$ time. Namely, and for simplicitly focusing on the case when
$n_{1}=n_{2}$, for which the goal is to do better than $n_{1}^{2}$ time, we
show how to output an $x\in\mathbb{R}^{d_{1}\times d_{2}}$ for which
$\|Cx^{\prime}-b\|_{1}\leq(1+\epsilon)\min_{x\in\mathbb{R}^{d_{1}d_{2}}}\|Cx-b\|_{1}$
in time $n^{3/2}\operatorname{poly}(d_{1}d_{2}/\epsilon)$. While this is more
expensive than solving Kronecker product least squares, the $1$-norm may lead
to significantly more robust solutions. From a technical perspective,
TensorSketch when applied to a vector actually destroys the $1$-norm of a
vector, preserving only its $2$-norm, so new ideas are needed. We show how to
use multiple TensorSketch matrices to implicitly obtain very crude estimates
to the so-called $\ell_{1}$-leverage scores of $C$, which can be interpreted
as probabilities of sampling rows of $C$ in order to obtain an
$\ell_{1}$-subspace embedding of the column span of $C$ (see, e.g., (Woodruff,
2014)). However, since $C$ has $n_{1}n_{2}$ rows, we cannot even afford to
write down the $\ell_{1}$-leverage scores of $C$. We show how to implicitly
represent such leverage scores and to sample from them without ever explicitly
writing them down. Balancing the phases of our algorithm leads to our overall
time.

###  1.2 Notation

We consider Kronecker Product of $q$ 2-d matrices $A_{1}\otimes
A_{2}\otimes...\otimes A_{q}$, where each $A_{i}\in\mathbb{R}^{n_{i}\times
d_{i}}$. We denote $\mathcal{A}=A_{1}\otimes A_{2}\otimes...\otimes A_{q}$,
$n=\prod_{i=1}^{q}n_{i}$, and $d=\prod_{i=1}^{q}d_{i}$. We denote $[n]$ as the
set $\\{1,2,3,\cdots,n\\}$. The $\ell_{p}$ norm for a vector
$x\in\mathbb{R}^{d}$ is defined as
$\|x\|_{p}=(\sum_{i=1}^{d}|x_{i}|^{p})^{1/p}$, where $x_{i}$ stands for the
$i$’th entry of the vector $x$. For any matrix $M$, we use $M_{i,*}$ to
represent the $i$’th row of $M$ and $M_{*,j}$ as the $j$’th column of $M$.

We define a _Well-Conditioned Basis_ and _Statistical Dimension_ as follow
(similar definitions can be found in Clarkson (2005); Dasgupta et al. (2009);
Sohler and Woodruff (2011); Meng and Mahoney (2013); Song et al. (2017a, b)
and Avron et al. (2016)):

######  Definition 1.1 (Well-Conditioned Basis).

Let $A$ be an $n\times m$ matrix with rank $d$, let $p\in[1,\infty)$, and let
$\|\cdot\|_{q}$ be the dual norm of $\|\cdot\|_{p}$, i.e., $1/p+1/q=1$. Then
an $n\times d$ matrix $U$ is an $(\alpha,\beta,p)$-well-conditioned basis for
the column space of $A$, if the columns of $U$ span the column space of $A$
and (1) $\|U\|_{p}\leq\alpha$, and (2) for all $z\in\mathbb{R}^{d}$,
$\|z\|_{q}\leq\beta\|Uz\|_{p}$.

Consider the classic Ridge Regression:
$\min_{x}\|Ax-b\|^{2}_{2}+\lambda\|x\|_{2}^{2}$. The _Statistical Dimension_
is defined as:

######  Definition 1.2 (Statistical Dimension).

Let $A$ be an $n\times m$ matrix with rank $d$ and singular values
$\sigma_{i},i\in[d]$. For $\lambda\geq 0$, the statistical dimension
$\operatorname{\mathtt{sd}}_{\lambda}(A)$ is defined as the quantity
$\operatorname{\mathtt{sd}}_{\lambda}(A)=\sum_{i\in[d]}1/(1+\lambda/\sigma_{i}^{2})$.

##  2 BACKGROUND:TensorSketch

We briefly introduce TensorSketch (Pagh, 2013; Avron et al., 2014) and how to
apply TensorSketch to the Kronecker product of multiple matrices efficiently
without explicitly computing the tensor product.111We refer readers to (Avron
et al., 2014) for more details about TensorSketch.

We want to find a oblivious subspace embedding $S$ such that for any
$x\in\mathbb{R}^{n}$, we have
$\|S\mathcal{A}x\|_{2}=(1\pm\epsilon)\|\mathcal{A}x\|_{2}$, where the notation
$a=(1\pm\epsilon b)$ stands for $(1-\epsilon)b\leq a\leq(1+\epsilon)b$, for
any $a,b\in\mathbb{R}$. Consider the $(i_{1},i_{2},...,i_{q})$’th column of
$\mathcal{A}$ ($i_{j}\in[d_{j}]$): $A_{1{{}_{*,i_{1}}}}\otimes
A_{2_{*,i_{2}}}\otimes\cdots\otimes A_{{q}_{*,i_{q}}}$. Assume the sketching
target dimension is $m$. TensorSketch is defined using $q$ 3-wise indepedent
hash functions $h_{i}:[n_{i}]\to[m]$, and $q$ 4-wise independent sign
functions $s_{i}:[n_{i}]\to\\{-1,1\\}$, $\forall i\in[q]$. Define hash
function $H:[n_{1}]\times[n_{2}]\times\cdots\times[n_{q}]\to[m]$ as
$H(i_{1},i_{2},...,i_{q})=((\sum_{k=1}^{q}h_{k}(i_{k}))\mod m)$, and sign
function $S:[n_{1}]\times[n_{2}]\times\cdots\times[n_{q}]\to\\{-1,1\\}$ as
$S(i_{1},i_{2},...,i_{q})=\prod_{k=1}^{q}s_{k}(i_{k})$. Applying Tensorsketch
to the Kronecker product of vectors $A_{1_{*,i_{1}}},...,A_{q_{*,i_{q}}}$ is
equivalent to applying CountSketch (Charikar et al., 2004) defined with $H$
and $S$ to the vector $(A_{1_{*,i_{1}}}\otimes...\otimes A_{q_{*,i_{q}}})$. To
apply Tensorsketch to $\mathcal{A}$, we just need to apply CountSketch defined
with $H$ and $S$ to the columns of $\mathcal{A}$ one by one.

Applying tensorsketch to the Kronecker product of
$(A_{1_{*,i_{1}}},...,A_{q_{*,i_{q}}})$ naively would require at least $O(n)$
time. Pagh (2013) shows that one can apply tensorsketch to the Kronecker
product of $(A_{1_{*,i_{1}}},...,A_{q_{*,i_{q}}})$ without explicitly
computing the Kronecker product of these vectors using the Fast Fourier
Transformation. Particularly, Pham and Pagh (2013) show that one only needs
$O(\sum_{j=1}^{q}\operatorname{\mathtt{nnz}}(A_{j_{*,i_{j}}})+qm\log(m))$ time
to compute $S(A_{1_{*,i_{1}}}\otimes...\otimes A_{q_{*,i_{q}}})$ where
$A_{i_{*,j}}$ stands for the $j$’th column of $A_{i}$. As $\mathcal{A}$ has
$d$ columns, computing $S\mathcal{A}$ takes
$O(d(\sum_{i=1}^{q}\operatorname{\mathtt{nnz}}(A_{i}))+dqm\log(m))$ time,
which is much smaller than $O(\prod_{i=1}^{q}n_{i})$, which is the least
amount of time one needs for explicitly computing $\mathcal{A}$. In the rest
of the paper, we assume that we compute $S\mathcal{A}$ using the above
efficient procedure without explicitly computing $A_{1}\otimes\ldots\otimes
A_{q}$.

##  3 TENSOR PRODUCT LEAST SQUARES REGRESSION

Consider the tensor product least squares regression problem
$\min_{x}\|\mathcal{A}x-b\|_{2}$ where $\mathcal{A}\in\mathbb{R}^{n\times d}$
and $b\in{\mathbb{R}}^{n}$. Let $S\in{\mathbb{R}}^{m\times n}$ be the matrix
form of TensorSketch of Section 2. We propose a TensorSketch -type Algorithm 1
for the tensor product regression problem. The following theorem shows that
the solution obtained from Alg. 1 is a good approximation of the optimal
solution of the original tensor product regression.

Let us define $\operatorname{OPT}$ to be the optimal cost of the optimization
problem, e.g., $\operatorname{OPT}=\min_{x}\|{\cal A}x-b\|_{2}$. The following
theorem shows that Alg. 1 computes a good approximate solution.

Algorithm 1 Tensor product regression

1:procedure TRegression($A,b,\epsilon,\delta$)

2: $m\leftarrow(d_{1}d_{2}\cdots d_{q}+1)^{2}(2+3^{q})/(\epsilon^{2}\delta)$

3: Choose $S$ to be an $m\times(n_{1}n_{2}\cdots n_{q})$ TensorSketch matrix

4: Compute $S(A_{1}\otimes A_{2}\otimes\cdots\otimes A_{q})$ and $Sb$

5: $\widetilde{x}\leftarrow\min_{x}\|S(A_{1}\otimes A_{2}\otimes\cdots\otimes
A_{q})x-Sb\|_{2}$

6: return $\widetilde{x}$

7:end procedure

######  Theorem 3.1.

(Tensor regression) Suppose $\widetilde{x}$ is the output of Algorithm 1 with
tensorsketch $S\in\mathbb{R}^{m\times n}$, where $m=8(d_{1}d_{2}\cdots
d_{q}+1)^{2}(2+3^{q})/(\epsilon^{2}\delta)$. Then the following approximation
$\|(A_{1}\otimes A_{2}\otimes\cdots\otimes
A_{q})\widetilde{x}-b\|_{2}\leq(1+\epsilon)\operatorname{OPT},$ holds with
probability at least $1-\delta$.

The proof of Theorem 3.1 can be found in Appendix C.1. Theorem 3.1 shows that
we can achieve an $\epsilon$-close solution by solving a much smaller
regression problem with a number of samples of order
$O(\operatorname{poly}(d/\epsilon))$, which is independent of the large
dimension $n$. Using the technique we introduced in Sec. 2, we can also
compute $S{\cal A}$ without explicitly computing the tensor product.

We can extend Theorem 3.1 to the nonnegative tensor product regression problem
$\min_{x\geq 0}\|(A_{1}\otimes A_{2}\otimes\cdots\otimes A_{q})x-b\|_{2},$
where $A_{i}\in{\mathbb{R}}^{n_{i}\times d_{i}},\,i=1,\ldots,q$ and
$b\in{\mathbb{R}}^{n_{1}n_{2}\cdots n_{q}}$. Suppose $x$ is the optimal
solution. Similarly, let $S\in{\mathbb{R}}^{m\times(n_{1}n_{2}\cdots n_{q})}$
be the matrix form of TensorSketch of Section 2. If $\widetilde{x}$ is the
optimal solution to $\min_{x\geq 0}\|S(A_{1}\otimes A_{2}\otimes\cdots\otimes
A_{q})x-Sb\|_{2}.$, we have the following:

######  Corollary 3.2.

(Sketch for tensor nonnegative regression) Suppose $\tilde{x}=\min_{x\geq
0}\|S\mathcal{A}x-Sb\|_{2}$ with tensorsketch $S\in\mathbb{R}^{m\times n}$,
where $m=8(d_{1}d_{2}\cdots d_{q}+1)^{2}(2+3^{q})/(\epsilon^{2}\delta)$. Then
the following approximation $\|(A_{1}\otimes A_{2}\otimes\cdots\otimes
A_{q})\widetilde{x}-b\|_{2}\leq(1+\epsilon)\operatorname{OPT}$ holds with
probability at least $1-\delta$, where $\operatorname{OPT}=\min_{x\geq
0}\left\|(A_{1}\otimes A_{2}\otimes\cdots\otimes A_{q}){x}-b\right\|_{2}$.

The proof of Corollary 3.2 can be found in Appendix C.1.

##  4  P-SPLINES

B-splines are local basis functions, consisting of low degree (e.g.,
quadratic, cubic) polynomial segments. The positions where the segments join
are called the knots. B-splines have local support and are thus suitable for
smoothing and interpolating data with complex patterns. Unfortunately, control
over smoothness is limited: one can only change the number and positions of
the knots. If there are no reasons to assume that smoothness is non-uniform,
the knots will be equally spaced and the only tuning parameter is their
(discrete) number. In contrast P-spline (Eilers and Marx, 1996) equally spaces
B-splines, discards the derivative completely, and controls smoothness by
regularizing the sum of squares of differences of coefficients. Specifically
Eilers and Marx Eilers and Marx (1996) proposed the P-spline recipe: (1) use a
(quadratic or cubic) B-spline basis with a large number of knots, say 10-50;
(2) introduce a penalty on (second or third order) differences of the B-spline
coefficients; (3) minimize the resulting penalized likelihood function; (4)
tune smoothness with the weight of the penalty, using cross-validation or AIC
to determine the optimal weight.

We give a brief overview of B-Splines and P-Splines below. Let $b$ and $u$,
each vectors of length $n$, represent the observed and explanatory variables,
respectively. Once a set of knots is chosen, the B-spline basis $A$ follows
from $u$. If there are $d$ basis functions then $A$ is $n\times d$. In the
case of normally distributed observations the model is $b=Ax+e$, with
independent errors $e$. In the case of B-spline regression the sum of squares
of residuals $\left\|b-Ax\right\|_{2}$ is minimized and the normal equations
$A^{\top}A\hat{x}=A^{\top}b$ are obtained; the explicit solution
$\hat{x}=(A^{\top}A)^{-1}A^{\top}b$ results. The P-spline approach minimizes
the penalized least-squares function

$\left\|b-Ax\right\|_{2}^{2}+\lambda\|Lx\|_{2}^{2},$ (1)

where $L\in{\mathbb{R}}^{p\times n}$ is a matrix that forms differences of
order $\ell$, i.e., $L_{\ell}x=\Delta^{\ell}x$. Examples of this matrix, for
$\ell=1$ and $\ell=2$ are :

$\displaystyle L_{1}=\begin{bmatrix}-1&1&0&0\\\ 0&-1&1&0\\\
0&0&-1&1\end{bmatrix},L_{2}=\begin{bmatrix}1&-2&1&0&0\\\ 0&1&-2&1&0\\\
0&0&1&-2&1\end{bmatrix}.$

The parameter $\lambda$ determines the influence of the penalty. If $\lambda$
is zero, we are back to B-spline regression; increasing $\lambda$ makes
$\hat{x}$, and hence $\hat{b}=A\hat{x}$, smoother.

Let $x^{*}$ denote $\operatorname*{argmin}_{x\in{\mathbb{R}}^{d}}\lVert
Ax-b\rVert_{2}^{2}+\lambda\lVert Lx\rVert_{2}^{2}$, and $\operatorname{OPT}$
denote $\lVert Ax^{*}-b\rVert_{2}^{2}+\lambda\lVert Lx^{*}\rVert_{2}^{2}$. In
general $x^{*}=(A^{\top}A+\lambda
L^{\top}L)^{-1}A^{\top}b=A^{\top}(AA^{\top}+\lambda LL^{\top})^{-1}b$, so
$x^{*}$ can be found in $O(\operatorname{\mathtt{nnz}}(A)\min(n,d))$ time
using an iterative method (e.g., LSQR). Our first goal in this section is to
design faster algorithms that find an approximate $\widetilde{x}$ in the
following sense:

$\lVert A\widetilde{x}-b\rVert_{2}^{2}+\lambda\lVert
L\widetilde{x}\rVert_{2}^{2}\leq(1+\epsilon)\operatorname{OPT}.\vspace{-4mm}$
(2)

###  4.1 Sketching for P-Spline

We first introduce a new definition of _Statistical Dimension_ that extends
the statistical dimension defined for Ridge Regression (i.e., $L$ is an
identity matrix in Eq. 1) (Avron et al., 2016) to P-Spline regression.

The problem (1) can also be analyzed by generalized singular value
decomposition (GSVD) Golub and Van Loan (2013). For matrices
$A\in{\mathbb{R}}^{n\times d}$ and $L\in{\mathbb{R}}^{p\times d}$ with
$\mathop{\rm rank}\nolimits(L)=p$ and $\mathop{\rm
rank}\nolimits\left(\left[\begin{array}[]{c}A\\\
L\end{array}\right]\right)=d$, the GSVD of $(A,L)$ is given by the pair of
factorizations $A=U\begin{bmatrix}\Sigma&{0}_{p\times(n-p)}\cr{0}_{(n-p)\times
p}&I_{d-p}\end{bmatrix}RQ^{\top}\quad\hbox{and}\quad
L=V\begin{bmatrix}\Omega&{0}_{p\times(n-p)}\end{bmatrix}RQ^{\top},$ where
$U\in{\mathbb{R}}^{m\times n}$ has orthonormal columns,
$V\in{\mathbb{R}}^{p\times p}$, $Q\in{\mathbb{R}}^{d\times d}$ are orthogonal,
$R\in\mathbb{R}^{d\times d}$ is upper triangular and nonsingular, and $\Sigma$
and $\Omega$ are $p\times p$ diagonal matrices: $\Sigma=\mathop{\rm
diag}\nolimits(\sigma_{1},\sigma_{2},\ldots,\sigma_{p})$ and
$\Omega=\mathop{\rm diag}\nolimits(\mu_{1},\mu_{2},\ldots,\mu_{p})$ with
$0\leq\sigma_{1}\leq\sigma_{2}\leq\ldots\leq\sigma_{p}<1\quad\hbox{and}\quad
1\geq\mu_{1}\geq\mu_{2}\geq\ldots\geq\mu_{p}>0,$ satisfying
$\Sigma^{\top}\Sigma+\Omega^{\top}\Omega=I_{p}$. The generalized singular
values $\gamma_{i}$ of $(A,L)$ are defined by the ratios
$\gamma_{i}=\sigma_{i}/\mu_{i}$ ($i=[p]$).

In this section we design an algorithm that is aimed at the case when $n\gg
d$. The general strategy is to design a distribution on matrices of size
$m$-by-$n$ ($m$ is a parameter), sample an $S$ from that distribution, and
solve $\widetilde{x}\equiv\operatorname*{argmin}_{x\in{\mathbb{R}}^{d}}\lVert
S(Ax-b)\rVert_{2}^{2}+\lambda\lVert Lx\rVert_{2}^{2}\,$.

The following lemma defines conditions on the distribution of $S$ that
guarantees Eq. (2) holds with constant probability (which can be boosted to
high probability by repetition and taking the minimum objective).

######  Lemma 4.1.

Let $x^{*}\in{\mathbb{R}}^{d}$, $A\in{\mathbb{R}}^{n\times d}$ and
$b\in{\mathbb{R}}^{n}$ as above. Let $U_{1}\in{\mathbb{R}}^{n\times d}$ denote
the first $n$ rows of an orthogonal basis for $\left[\begin{smallmatrix}A\\\
\sqrt{\lambda}L\end{smallmatrix}\right]\in{\mathbb{R}}^{(n+p)\times d}$. Let
sketching matrix $S\in{\mathbb{R}}^{m\times n}$ have a distribution such that
with constant probability

$\displaystyle\mathrm{(\@slowromancap i@)}~{}\lVert
U_{1}^{\top}S^{\top}SU_{1}-U_{1}^{\top}U_{1}\rVert_{2}\leq 1/4,$

and

$\displaystyle\mathrm{(\@slowromancap
ii@)}~{}\|U_{1}^{\top}(S^{\top}S-I)(b-Ax^{*})\|_{2}\leq\sqrt{\epsilon\operatorname{OPT}/2}.$

Let $\widetilde{x}$ denote
$\operatorname*{argmin}_{x\in{\mathbb{R}}^{d}}\lVert
S(Ax-b)\rVert_{2}^{2}+\lambda\lVert Lx\rVert_{2}^{2}$. Then with probability
at least $9/10$,

$\displaystyle\lVert A\widetilde{x}-b\rVert_{2}^{2}+\lambda\lVert
L\widetilde{x}\rVert_{2}^{2}\leq(1+\epsilon)\operatorname{OPT}.$

Define the _statistical dimension_ for P-Splines as follows:

######  Definition 4.2 (Statistical Dimension for P-Splines).

For S-Spline in Eq. (1), the statistical dimension is defined as
$\operatorname{\mathtt{sd}}_{\lambda}(A,L)=\sum_{i}1/(1+\lambda/\gamma_{i}^{2})+d-p$.

The following theorem shows that there is a sparse subspace embedding matrix
$S\in\mathbb{R}^{m\times n}$ (e.g., CountSketch), with $m\geq
K(\operatorname{\mathtt{sd}}_{\lambda}(A,L)/\epsilon+\operatorname{\mathtt{sd}}_{\lambda}(A,L)^{2})$,
that satisfies Property (i@) and (ii@) of Lemma 4.1, and hence achieves an
$\epsilon-$approximation solution to problem 1:

######  Theorem 4.3.

(P-Spline regression) There is a constant $K>0$ such that for $m\geq
K(\epsilon^{-1}\operatorname{\mathtt{sd}}_{\lambda}(A,L)+\operatorname{\mathtt{sd}}_{\lambda}(A,L)^{2})$
and $S\in{\mathbb{R}}^{m\times n}$ a sparse embedding matrix (e.g.,
Countsketch) with $SA$ computable in $O(\operatorname{\mathtt{nnz}}(A))$ time,
Property (i@) and (ii@) of Lemma 4.1 apply, and with constant probability the
corresponding
$\widetilde{x}=\operatorname*{argmin}_{x\in{\mathbb{R}}^{d}}\lVert
S(Ax-b)\rVert_{2}+\lambda\lVert Lx\rVert_{2}^{2}$ is an $\epsilon$-approximate
solution to $\min_{x\in{\mathbb{R}}^{d}}\lVert
b-Ax\rVert_{2}^{2}+\lambda\lVert Lx\rVert_{2}^{2}$.

Note $\operatorname{\mathtt{sd}}_{\lambda}(A,L)$ is upper bounded by $d$. The
above theorem shows that the statistical dimension allows us to design smaller
sketch matrices whose size only depends on
$O(\operatorname{poly}(\operatorname{\mathtt{sd}}_{\lambda}(A,L)/\epsilon))$
instead of $O(\operatorname{poly}(d/\epsilon))$, without sacrificing the
approximation accuracy.

###  4.2 Tensor Sketching for Multi-Dimensional P-Spline

Algorithm 2 P-Spline Tensor product regression

1:procedure PTRegression($A,b,K,L,\epsilon,\delta$)

2: $m\leftarrow
K(\epsilon^{-1}\operatorname{\mathtt{sd}}_{\lambda}(\mathcal{A},L)+\operatorname{\mathtt{sd}}_{\lambda}(\mathcal{A},L)^{2})$

3: Choose $S$ to be a $m\times n$ TensorSketch matrix

4: Compute $S(A_{1}\otimes A_{2}\otimes\cdots\otimes A_{q})$ and $Sb$

5: $\widetilde{x}=\operatorname*{argmin}_{x\in{\mathbb{R}}^{d}}\lVert
S(Ax-b)\rVert_{2}^{2}+\lambda\lVert Lx\rVert_{2}^{2}$.

6: return $\widetilde{x}$

7:end procedure

Tensor products allow a natural extension of one-dimensional P-spline
smoothing to multi-dimensional P-Spline. We focus on 2-dimensional P-Spline
but our results can be generalized to the multi-dimensional setting. Assume
that in addition to $u$ we have a second explanatory variable $v$. We have
data triples $(u_{i},v_{j},b_{{(i-1)\cdot n_{2}+j}})$ for $i=1,\ldots,n_{1}$
and $j=1,\ldots,n_{2}$. We seek a smooth surface $f(u,v)$ which gives a good
approximation to the response $b$. Let $A_{1}$, $n\times d_{1}$, be a B-spline
basis along $u$, and $A_{2}$, $n\times d_{2}$, be a B-spline basis along $v$.
We form the tensor product basis as $A_{1}\otimes A_{2}$. When $n_{1}$ and
$n_{2}$ are large, we do not compute $A_{1}\otimes A_{2}$. We apply
tensorsketch here to avoid explicitly forming $A_{1}\otimes A_{2}$ to speed up
computation.

Let $X=[x_{kl}]$ be a $d_{1}\times d_{2}$ matrix of coefficients. Then, for
given $X$, the value fit at $(u,v)$ is
$f(u,v)=\sum_{k}\sum_{l}A_{2,k}(v)A_{1,l}(u)x_{kl}$ and so $X$ may be chosen
using least squares by minimizing $\sum_{i,j}[b_{{(i-1)\cdot
n_{2}+j}}-f(u_{i},v_{i})]^{2}=\sum_{i}\left[b_{{(i-1)\cdot
n_{2}+j}}-\sum_{k}\sum_{l}A_{2,k}(v_{i})A_{1,l}(u_{i})x_{kl}\right]^{2}.$
Using Kronecker product, the above minimization can be written in the form
$\min\left\|b-{\cal A}x\right\|_{2},$ where ${\cal A}=A_{1}\otimes
A_{2}\in{\mathbb{R}}^{n_{1}n_{2}\times d_{1}d_{2}}$ and $x=\textrm{ vec}(X)$.
Again the P-spline approach minimizes the penalized least-squares function
$\left\|b-(A_{1}\otimes A_{2})x\right\|_{2}^{2}+\lambda\|Lx\|_{2}^{2}.$
Consider the tensor p-spline regression problem $\min_{x}\|(A_{1}\otimes
A_{2}\otimes\cdots\otimes A_{q})x-b\|_{2}^{2}+\lambda\lVert Lx\rVert_{2}^{2},$
where $L\in{\mathbb{R}}^{p\times n}$, $A_{i}\in{\mathbb{R}}^{n_{i}\times
d_{i}},\,i=1,\ldots,q$ and $b\in{\mathbb{R}}^{n}$. Let
$S\in{\mathbb{R}}^{m\times n}$ be the matrix form of TensorSketch of Section
2. Algorithm 2 summarizes the procedure for efficiently solving multi-
dimensional P-Spline.

Let ${\cal A}=A_{1}\otimes A_{2}\otimes\cdots\otimes A_{q}$. Replacing the
matrix $A$ in Theorem 4.3 with $\mathcal{A}$, we have the following corollary
for multi-dimensional P-Spline:

######  Corollary 4.4 (P-Spline tensor regression).

Suppose $\lambda\leq\sigma_{1}^{2}/\epsilon$. There is a constant $K>0$ such
that for $m\geq K(\epsilon^{-1}\operatorname{\mathtt{sd}}_{\lambda}({\cal
A},L)+\operatorname{\mathtt{sd}}_{\lambda}({\cal A},L)^{2})$ and
$S\in{\mathbb{R}}^{m\times n}$ a TensorSketch matrix with $S{\cal A}$
computable in $O(\operatorname{\mathtt{nnz}}({\cal A}))$ time, Property (i@)
and Property (ii@) of Lemma 4.1 apply, and with constant probability the
corresponding
$\widetilde{x}=\operatorname*{argmin}_{x\in{\mathbb{R}}^{d}}\lVert S({\cal
A}x-b)\rVert_{2}+\lambda\lVert Lx\rVert_{2}^{2}$ is an $\epsilon$-approximate
solution to $\min_{x\in{\mathbb{R}}^{d}}\lVert b-{\cal
A}x\rVert_{2}^{2}+\lambda\lVert Lx\rVert_{2}^{2}$.

##  5 TENSOR PRODUCT ABSOLUTE DEVIATION REGRESSION

We extend our previous results for the $\ell_{2}$ norm (i.e., least squares
regression) to general $\ell_{p}$ norms, with a focus on $p=1$ (i.e., absolute
deviation regression). Specifically we consider $\min_{x}\|{\cal A}x-b\|_{p}$,
where ${\cal A}=A_{1}\otimes A_{2}\otimes\cdots\otimes A_{q}$. We will show in
this section that with probability at least $2/3$, we can quickly find an
$\tilde{x}$ for which $\|{\cal
A}\tilde{x}-b\|_{1}\leq(1+\epsilon)\min_{x}\|{\cal A}x-b\|_{1}.$

As in Clarkson and Woodruff (2013), for each $i$, in
$O(\operatorname{\mathtt{nnz}}(A_{i})\log(d_{i}))+O(r_{i}^{3})$ time we can
replace the input matrix $A_{i}\in{\mathbb{R}}^{n_{i}\times d_{i}}$ with a new
matrix with the same column space of $A_{i}$ and full column rank. We
therefore assume ${\cal A}$ has full rank in what follows.

Suppose $S$ is the TensorSketch matrix defined in Section 2. Let
$w_{i}\in\mathbb{N}$ and assume $w_{i}\mid n_{i}$. Split $A_{i}$ into
$n_{i}/w_{i}$ matrices $A_{1}^{(i)},\ldots,A_{n_{i}/w_{i}}^{(i)}$, each
$w_{i}\times d_{i}$, so that $A_{j}^{(i)}$ is the submatrix of $A_{i}$ indexed
by the $j$-th block of $w_{i}$ rows. Note ${\cal A}$ can be written as:

$\displaystyle\begin{bmatrix}A_{1}^{(1)}\otimes
A_{1}^{(2)}\otimes\cdots\otimes A_{1}^{(q-1)}\otimes A_{1}^{(q)}\\\
A_{1}^{(1)}\otimes A_{1}^{(2)}\otimes\cdots\otimes A_{1}^{(q-1)}\otimes
A_{2}^{(q)}\\\ \vdots\\\ A_{n_{1}/w_{1}}^{(1)}\otimes
A_{n_{2}/w_{2}}^{(2)}\otimes\cdots\otimes A_{n_{q-1}/w_{q-1}}^{(q-1)}\otimes
A_{n_{q}/w_{q}}^{(q)}.\end{bmatrix}.$

For each $A_{i_{1}}^{(1)}\otimes A_{i_{2}}^{(2)}\otimes\cdots\otimes
A_{i_{q}}^{(q)}$, we can use the TensorSketch matrix $S_{i_{1}i_{2}\ldots
i_{q}}\in{\mathbb{R}}^{m\times\prod_{i=1}^{q}w_{i}}$, where we set $m\geq
100\prod_{i=1}^{q}d_{i}^{2}(2+3^{q})/\epsilon^{2}$, such that with probability
at least $.99$, $\|S_{i_{1}i_{2}\ldots i_{q}}A_{i_{1}}^{(1)}\otimes
A_{i_{2}}^{(2)}\otimes\cdots\otimes
A_{i_{q}}^{(q)}x\|_{2}=(1\pm\epsilon)\|A_{i_{1}}^{(1)}\otimes
A_{i_{2}}^{(2)}\otimes\cdots\otimes A_{i_{q}}^{(q)}x\|_{2}$ simultaneously for
all $x\in{\mathbb{R}}^{d}$ as $S_{i_{1}i_{2}\ldots i_{q}}$ is an oblivious
subspace embedding (Lemma. B.3) . Now we use Algorithm 4 from (Liang et al.,
2014) to boost the success probability by computing $t=O(\log(1/\delta))$
independent TensorSketch products $S^{(j)}_{i_{1}i_{2}\ldots
i_{q}}A_{i_{1}}^{(1)}\otimes A_{i_{2}}^{(2)}\otimes\cdots\otimes
A_{i_{q}}^{(q)}$, $j=[t]$, each with only constant success probability, and
then running a cross validation procedure similar to that in Algorithm 4 of
Liang et al. (2014), to find one which succeeds with probability at least
$1-\delta$:

######  Lemma 5.1 ((Liang et al., 2014)).

For $\delta,\epsilon\in(0,1)$, let $t=O(\log 1/\delta)$ and $m\geq
100\prod_{i=1}^{q}d_{i}^{2}(2+3^{q})/\epsilon^{2}$. Running algorithm
Algorithm 4 in (Liang et al., 2014) with parameters $t,m$, we can obtain a
tensorsketch $S\in\mathbb{R}^{m\times\prod_{i=1}^{q}w_{i}}$ such that with
probability at least $1-\delta$, $\|SA_{i_{1}}^{(1)}\otimes
A_{i_{2}}^{(2)}\otimes\cdots\otimes
A_{i_{q}}^{(q)}x\|_{2}=(1\pm\epsilon)\|A_{i_{1}}^{(1)}\otimes
A_{i_{2}}^{(2)}\otimes\cdots\otimes A_{i_{q}}^{(q)}x\|_{2}$ for all
$x\in\mathbb{R}^{d}$.

After computing a tensorsketch $S_{i_{1}i_{2}\ldots
i_{q}}\in\mathbb{R}^{m\times\prod_{i=1}^{q}w_{i}}$ using lemma 5.1 for each
row-block $A_{i_{1}}^{(1)}\otimes...\otimes A_{i_{q}}^{(q)}$, we compose a
single tensorsketch for $\mathcal{A}$ as
$S={\rm{diag}}(S_{1,\cdots,1},\cdots,S_{i_{1},\cdots,i_{q}},\cdots,S_{(n_{1}/w_{1}),\cdots,(n_{q}/w_{q})})\in\mathbb{R}^{m\prod_{i=1}^{q}n_{i}/w_{i}\times\prod_{i=1}^{q}{n_{i}}}$,
which is defined as:

$\displaystyle\begin{bmatrix}S_{1,\cdots,1}&&&\\\ &\ddots&&&\\\
&&S_{i_{1},\cdots,i_{q}}&&\\\ &&&\ddots&\\\
&&&&S_{(n_{1}/w_{1}),\cdots,(n_{q}/w_{q})}\end{bmatrix},$

where each block on the diagonal is from Lemma 5.1. Note that ${\cal A}$ has
in total $\prod_{i=1}^{q}(n_{i}/w_{i})$ many blocks. Using Lemma 5.1 with a
union bound over all blocks of ${\cal A}$, we have the following theorem which
shows $S$ is an oblivious subspace embedding for ${\cal A}$ in $\ell_{2}$
norm:

######  Theorem 5.2 ($\ell_{2}$ OSE for tensor matrices).

Given $\delta,\epsilon\in(0,1)$, let
$S\in\mathbb{R}^{m\prod_{i=1}^{q}n_{i}/w_{i}\times\prod_{i=1}^{q}n_{i}}$
denote the matrix that has $\prod_{i=1}^{q}n_{i}/w_{i}$ diagonal block
matrices where each diagonal block $S_{i_{1}i_{2}\ldots
i_{q}}\in\mathbb{R}^{m\times\prod_{i=1}^{q}w_{i}}$ is from Lemma 5.1. With
probability at least $1-\prod_{i=1}^{q}(n_{i}/w_{i})\delta$, $\|S{\cal
A}x\|_{2}=(1\pm\epsilon)\|{\cal A}x\|_{2},\forall x\in\mathbb{R}^{d}$.

It is known that for any matrix $A\in\mathbb{R}^{n\times r}$, we can compute a
change of basis $U\in\mathbb{R}^{r\times r}$ such that $AU$ is an
$(\alpha,\beta,p)$ well-conditioned basis of $A$ (see Definition 1.1), in time
polynomial with respect to $n,r$ (Dasgupta et al., 2009). Specifically,
Theorem 5 in (Dasgupta et al., 2009) shows that we can compute a change of
basis $U$ for which $AU$ is a well-conditioned basis of $A$ in time
$O(nr^{5}\log(n))$ time. However we cannot afford to directly use the results
from (Dasgupta et al., 2009) to compute a well-conditioned basis for
$\mathcal{A}$, which requires time at least $\Omega(n)$. Instead we compute a
well-conditioned basis for $\mathcal{A}$ through a sketch $S{\cal A}$ where
$S$ is the tensorsketch from Theorem 5.2.

Specifically we define the following procedure Condition(${\cal A}$) for
computing a well-conditioned basis for ${\cal A}$. Given ${\cal
A}=A_{1}\otimes\cdots\otimes A_{q}\in\mathbb{R}^{n\times r}$, 1) Compute
$S{\cal A}$; 2) Compute a $d\times d$ change of basis matrix $U$ so that
$S{\cal A}U$ is an $(\alpha,\beta,p)$-well-conditioned basis of the column
space of $S{\cal A}$; 3) Output ${\cal A}U/(d\gamma_{p})$, where
$\gamma_{p}\equiv\sqrt{2}t^{1/p-1/2}$ for $p\leq 2$, and
$\gamma_{p}\equiv\sqrt{2}w^{1/2-1/p}$ for $p\geq 2$, where
$t=\prod_{i=1}^{q}m_{i}$ and $w=\prod_{i=1}^{q}w_{i}$. The following Lemma 5.3
(the proof can be found in the Appendix) is the analogue of that in Clarkson
et al. (2013) proved for the Fast Johnson Lindenstauss Transform.

######  Lemma 5.3.

For any $p\geq 1$. Condition$({\cal A})$ computes ${\cal A}U/(d\gamma_{p})$
which is an $(\alpha,\beta\sqrt{3}d(tw)^{|1/p-1/2|},p)$-well-conditioned basis
of ${\cal A}$, with probability at least
$1-\prod_{i=1}^{q}(n_{i}/w_{i})\delta$.

Lemma 5.3 indicates that Condition$({\cal A})$ computes a
$(\alpha,\beta\cdot\operatorname{poly}(\max(d,\log n)),p)$-well-conditioned
basis. A well-conditioned basis can be used to solve $\ell_{p}$ regression
problems, via sampling a subset of rows of the well-conditioned basis
$\mathcal{A}U$ with probabilities proportional to the $p$-th power of the
$\ell_{p}$ norm of the rows (Woodruff, 2014). However the first issue for
sampling is that we cannot afford to compute $\mathcal{A}U$ as this requires
$O(\operatorname{\mathtt{nnz}}(\mathcal{A})d)$ time. To fix this, we apply a
Gaussian sketch matrix $G\in\mathbb{R}^{d\times\log(n)}$ with i.i.d normal
random variables (Drineas et al., 2011) on the right hand side of
$\mathcal{A}U$. Note that $\mathcal{A}UG$ can be computed efficiently by first
compuing $UG$ and then $\mathcal{A}(UG)$ in time
$O(d^{2}\log(n)+nnz(\mathcal{A})\log(n))$. The second issue is that even with
$\mathcal{A}UG$, computing the $\ell_{p}$ norm of each row of $\mathcal{A}UG$
takes $O(n)$ time, but we want sublinear time. This leads us to the following
sampling technique.

The high level idea is that since $\mathcal{A}UG$ only has $O(\log(n))$
columns, we can afford to sample columns of $\mathcal{A}UG$ with probability
proportional to the $p$-th power of the $\ell_{p}$ norms of the columns, if we
can efficiently estimate the $\ell_{p}$ norms of the columns (note that
naïvely computing the $\ell_{p}$ norm of a column also takes $O(n)$ time). Let
us denote the first of column of $UG$ as $e\in\mathbb{R}^{\log(n)}$. We focus
on how to efficiently estimate the $\ell_{p}$ norm of the first column of
$\mathcal{A}UG$, which is $\mathcal{A}e$, and all the left columns can be
estimated in the same way. We first reshape the vector $\mathcal{A}e$ into a
2-d matrix

$\displaystyle M=A_{1}\otimes...\otimes
A_{q_{1}}E(A_{q_{1}+1}\otimes...\otimes A_{q})^{\top},\vspace{-2mm}$ (3)

where $M\in\mathbb{R}^{(n_{1}n_{2}...n_{q_{1}})\times(n_{q_{1}+1}...n_{q}})$,
$E$ is obtained from reshaping $e$ into a
$(d_{1}d_{2}...d_{q_{1}})\times(d_{q_{1}+1}..d_{q})$ matrix and
$q_{1}\in[1,q]$ is chosen such that
$(n_{1}n_{2}...n_{q_{1}})\approx(n_{q_{1}+1}...n_{q})$. Namely we reshape the
column $\mathcal{A}e$ into a (nearly) square matrix. Focusing now on $p=1$,
note that
$\|\mathcal{A}e\|_{1}=\sum_{i=1}^{n_{q_{1}+1}...n_{q}}\|M_{*,i}\|_{1}$. Hence
to estimate $\|\mathcal{A}e\|_{1}$, we only need to estimate the $\ell_{1}$
norm of the columns of $M$.

Let us apply a sketch matrix
${R}\in\mathbb{R}^{O(\log(n))\times\prod_{i=1}^{q_{1}}n_{i}}$, whose entries
are sampled i.i.d. from the Cauchy distribution, to the left hand side of $M$.
For the $i$’th column of $M$, let us define random variables
$z_{l}^{i}=R_{l,*}^{\top}M_{i}$, for $l\in[O(\log(n))]$, where $R_{l,*}$ is
the $l$’th row of $R$. Due to the 1-stability property of the Cauchy
distribution, we have that $\\{z_{l}^{i}\\}_{l=1}^{O(\log n)}$ are $O(\log n)$
independent Cauchys scaled by $\|M_{i}\|_{1}$. Applying a Chernoff bound to
independent half-Cauchys (see Claims 1, 2 and Lemmas 1, 2 in Indyk (2006)), we
have $0.5\|M_{i}\|_{1}\leq\texttt{median}_{l\in[O(\log
n)]}\\{|z_{j}^{i}|\\}\leq 1.5\|M_{i}\|_{1}$ with probability at least
$1-2e^{-cO(\log(n))}$, with constant $c\geq 0.07$. Denote the median of
$\\{|z_{l}^{i}|\\}_{l=1}^{O(\log(n))}$ as $\lambda_{i}$. By a union bound over
all columns of $M$, we have that with probability at least
$1-n_{[q]\backslash[q_{1}]}2e^{-cO(\log(n))}$:

$\displaystyle\lambda_{i}=(1\pm 0.5)\|M_{i}\|_{1},\forall
i\in[n_{[q]\backslash[q_{1}]}],$ (4)
$\displaystyle\lambda_{e}=\sum_{i=1}^{n_{[q]\backslash[q_{1}]}}\lambda_{i}=(1\pm
0.5)\|\mathcal{A}e\|_{1}$ (5)

where $n_{[q]\backslash[q_{1}]}=\prod_{i=q_{1}+1}^{q}n_{i}$.

Algorithm 3 $\ell_{1}$ tensor product regression

1:procedure L1TRregression($A,b,\epsilon,\delta$)

2: Construct a tensorsketch
$S\in\mathbb{R}^{(m\prod_{i=1}^{q}\frac{n_{i}}{w_{i}})\times n}$.

3: Run Condition(${\cal A}$) using $S$ to compute $U/(d\gamma_{p})$.

4: Generate a Gaussian matrix $G\in\mathbb{R}^{d\times O(\log(n))}$ and a
Cauchy sketch matrix $R\in\mathbb{R}^{\log(n)\times n}$.

5: for for each column $e$ in $UG$ do

6: Reshape ${\cal A}e$ to $M$ (Eq. 3).

7: Compute $\lambda_{i}$ and $\lambda_{e}=\sum_{i}\lambda_{i}$ (Eq 4 and 5).

8: end for

9: for $i\in[\sqrt{\prod_{i=1}^{q}w_{i}}\operatorname{poly}(d)$] do

10: Sample a column $(AUG)_{*,e}$ with probability proportional to
$\lambda_{e}$. . $\triangleright$ $e\in{[O(\log(n))]}$

11: Reshape $(AUG)_{*,e}$ to $M$, sample a column $M_{*,j}$ with probability
proportional to $\lambda_{j}$.

12: Sample an entry $M_{k,j}$ with probability proportional to $|M_{k,j}|$.

13: Convert $(k,j)$ back to the corresponding row index in ${\cal A}$, denoted
as $r_{i}$.

14: end for

15: $\widetilde{x}\leftarrow\min_{x}\|D(A_{1}\otimes A_{2}\otimes\cdots\otimes
A_{q})x-Db\|_{1},$ where $D$ is a diagonal matrix to select the
$r_{1},r_{2},...,r_{N}$-th rows of ${\cal A}$ and $b$ with
$N=\sqrt{\prod_{i=1}^{q}w_{i}\operatorname{poly}(d)}$.

16: return $\widetilde{x}$ $\triangleright$
$\widetilde{x}\in{\mathbb{R}}^{d_{1}d_{2}\cdots d_{q}}$

17:end procedure

Note that since $\mathcal{A}UG$ only has $O(\log n)$ columns, we can afford to
compute the $\ell_{1}$ norm of all the columns using the above procedure. Let
us denote the $\ell_{1}$ norms of the columns of $AUG$ by
$\lambda_{1},\lambda_{2},...,\lambda_{O(\log n)}$. We can sample a column
$(\mathcal{A}UG)_{*,i}$ with probability proportional to $\lambda_{i}$ (Line
10 in Alg. 3). Once we sample a column $\mathcal{A}UG_{*,i}$, we need to
sample an entry $j\in[n]$ with probability proportional to the absolute value
of the entry $|(\mathcal{A}UG)_{j,i}|$. As we cannot afford to compute
$|(\mathcal{A}UG)_{j,i}|$ for all $j\in[n]$, we use the reshaped 2-d matrix
$M$ of $({\cal A}UG)_{*,i}$. Note that sampling an entry in $M$ with
probability proportional to the absolute values of entries of $M$ is
equivalent to sampling an entry $j\in[n]$ from $({\cal A}UG)_{*,i}$ with
probability proportional to the absolute value of the entries of $({\cal
A}UG)_{*,i}$. We first sample a column $M_{*,j}$ from all the columns of $M$
with probability proportional to $\lambda_{j}$ for
$j\in[\prod_{i=q_{1}+1}^{q}n_{i}]$. We then sample an entry $M_{k,j}$ with
probability proportional to $|M_{k,j}|$ for $k\in[\prod_{i=1}^{q_{1}}n_{i}]$.
Noting that $k\in[\prod_{i=1}^{q_{1}}n_{i}]$ and
$j\in[\prod_{i=q_{1}+1}^{q}n_{i}]$, the pair $(k,j)$ uniquely determines a
corresponding row index $r$ in ${\cal A}$, for some
$r\in[\prod_{i=1}^{q}n_{i}]$. Hence we successfully sample a row from ${\cal
A}UG$ without ever computing the $\ell_{p}$ norm of the rows. The above
sampling procedure is summarized in Line 10 to Line 13 in Alg. 3. We use the
above procedure to sample $\sqrt{\prod_{i=1}^{q}w_{i}}\operatorname{poly}(d)$
rows of $\mathcal{A}$. Let $D$ be a diagonal matrix that selects the
corresponding sampled rows from $\mathcal{A}$. We can now solve a smaller ADL
problem as $\min_{x}\|D\mathcal{A}-Db\|_{1}$. Note that our analysis focuses
on the $\ell_{1}$ norm. We can extend the analysis to general $\ell_{p}$ norms
by using a sketching matrix $R\in\mathbb{R}^{O(\log
n)\times\prod_{i=1}^{q_{1}}n_{i}}$ with entries sampled i.i.d. from a
$p$-stable distribution for $p\in[1,2]$.

We now present our main theorem and defer the proofs to the appendix:

######  Theorem 5.4.

(Main result) Given $\epsilon\in(0,1)$, ${\cal A}\in{\mathbb{R}}^{n\times d}$
and $b\in{\mathbb{R}}^{n}$, Alg. 3 computes $\widehat{x}$ such that with
probability at least $1/2$, $\|{\cal
A}\hat{x}-b\|_{1}\leq(1+\epsilon)\min_{x\in{\mathbb{R}}^{d}}\|{\cal
A}x-b\|_{1}$. For the special case when $q=2$, $n_{1}=n_{2}$, the algorithm’s
running time is
$O({n_{1}}^{3/2}\operatorname{poly}(\prod_{i=1}^{2}d_{i}/\epsilon))$.

For the special case where $q=2$ and $n_{1}=n_{2}$, we can see from theorem
5.4 our algorithm computes $\widehat{x}$ in time
$O({n_{1}}^{3/2}\operatorname{poly}(\hat{d}))$, which is faster than
$O({n_{1}}^{2})$—the time needed for forming $A_{1}\otimes A_{2}$. Note that
we can run Alg. 3 $O(\log(1/\delta))$ times independently and pick the best
solution among these independent runs to boost the success probability to be
$1-\delta$, for $\delta\in[0,1)$.

##  6 NUMERICAL EXPERIMENTS

We generate matrices $A_{1}$, $A_{2}$ and $b$ with all entries sampled i.i.d
from a normal distribution. The baseline we compared to is directly solving
regression without sketching. We let ${\tt T_{1}}$ be the time for directly
solving the regression problem, and ${\tt T_{2}}$ be the time of our
algorithm. The time ratio is $r_{t}={\tt T_{2}}/{\tt T_{1}}$. The relative
residual percentage is defined by $r_{e}=\frac{100\left|\,\left\|(A_{1}\otimes
A_{2})\widetilde{x}-b\right\|_{2}-\left\|(A_{1}\otimes
A_{2})x^{*}-b\right\|_{p}\,\right|}{\left\|(A_{1}\otimes
A_{2})x^{*}-b\right\|_{p}},$ where $\widetilde{x}$ is the output of our
algorithms and $x^{*}$ is the optimal solution. Throughout the simulations, we
use a moderate input matrix size in order to accommodate the brute force
algorithm and to compare to the exact solution.

######  Example 6.1 ($\ell_{2}$ Regression).

We create a design matrix with moderate size by fixing $n_{1}=n_{2}=300$ and
$d_{1}=d_{2}=15$. Thus ${\cal A}\in{\mathbb{R}}^{90000\times 225}$. We do 10
rounds to compute the mean values of $r_{e}$ and $r_{t}$, which is reported in
Table 1 (Left).

From Table 1 (left), we can see that when the number $m$ of sampled rows in
Algorithm 1 increases from 8000 to 12000, the mean values of $r_{e}$ decrease
while the mean values of $r_{t}$ increase. In general we can see that we can
achieve around $1\%$ relative error while being 5 times faster than the direct
method.

Table 1: Examples 6.1 and 6.2: the values of $r_{e}$ and $r_{t}$ with respect
to different sampling parameters $m$.

| $m$ | $r_{e}$ | $r_{t}$ |  | $m$ | $r_{e}$ | $r_{t}$
---|---|---|---|---|---|---|---
$\ell_{2}$ | 8000 | 1.79% | 0.11 | $\ell_{1}$ | 8000 | 1.89% | 0.06
12000 | 1.24% | 0.18 | 12000 | 1.33% | 0.11
16000 | 1.01% | 0.25 | 16000 | 0.992% | 0.18

######  Example 6.2 ($\ell_{1}$ Regression).

We set $n_{1}=n_{2}=300$, $d_{1}=d_{2}=15$. For $\min_{x}\|Ax-b\|_{1}$, we
solve it by a Linear Programming solver in Gurobi (Gurobi Optimization, 2016).
We tested different numbers of sampled rows $m$ (the number of rows in $D$ in
Alg. 3). The results are summarized in Table 1 (Right).

As we can see from Table 1 (right), our method is around 10 times faster than
directly solving the problem, with relative error only around $1\%$.

######  Example 6.3 (P-Spline Regression).

For P-spline regression, $L_{3}$ is fixed. We use 30 knots and cubic
B-splines. The data
$u=(u_{i})\in{\mathbb{R}}^{n_{1}},v=(v_{j})\in{\mathbb{R}}^{n_{2}}$ and
$b\in{\mathbb{R}}^{n_{1}n_{2}}$ are generated i.i.d from the normal
distribution. We compute the B-spline basis matrices
$A_{1}\in{\mathbb{R}}^{n_{1}\times d_{1}}$ and
$A_{2}\in{\mathbb{R}}^{n_{2}\times d_{2}}$ separately, if there are $d_{1}$
and $d_{2}$ basis functions for the B-spline. We set $n_{1}=n_{2}$,
$d_{1}=d_{2}$, with $n_{1}n_{2}=10^{4}$ and $d_{1}d_{2}=529$. The original and
sketched P-spline regression problem are both solved by Regularization Tools
(Hansen, 1994) via computing their GSVDs. We test different choices of
$\lambda$. The results are shown in Table 2.

From Table 2, sampling only 20% of the rows can give around $0.05\%$ relative
error, while the computation time is half of the time of directly solving the
p-spline.

Table 2: Example 6.3: the mean values of $r_{e}$ and $r_{t}$ with respect to
different sampling parameters $m$ and regularization parameters $\lambda$.

$\lambda$ | $m$ | $r_{e}$ | $r_{t}$
---|---|---|---
1 | 2000 | 4.43e-2% | 0.52
| 4000 | 2.99e-2% | 0.70
| 6000 | 7.92e-2% | 0.91
0.1 | 2000 | 6.98e-2% | 0.47
| 4000 | 4.07e-2% | 0.72
| 6000 | 5.41e-2% | 0.94
0.01 | 2000 | 3.78e-2% | 0.46
| 4000 | 1.07e-1% | 0.71
| 6000 | 2.97e-2% | 0.95

##  7 CONCLUSION

We propose algorithms for efficiently solving tensor product least squares
regression, regularized P-splines, as well as tensor product least absolute
deviation regression, using sketching techniques. Our main contributions are:
(1) we apply tensorsketch to least square regression problems, (2) we propose
new statistical dimension measures for P-splines, extending the previous
statistical dimension defined only for classic Ridge regression, and (3) we
extend tensorsketch to $\ell_{p}$ norms and propose an algorithm that can
solve tensor product $\ell_{1}$ regression in time sublinear in the time for
explicitly computing the tensor product. Simulation results support our
theorems and demonstrate that our algorithms are much faster than brute-force
algorithms and can achieve approximate solutions that are close to optimal.

## ACKNOWLEDGEMENT

Huaian Diao is supported in part by the Fundamental Research Funds for the
Central Universities under the grant 2412017FZ007. Wen Sun is supported in
part by Office of Naval Research contract N000141512365.

## References

  * Avron et al. [2014] H. Avron, H. Nguyen, and D. Woodruff.  Subspace embeddings for the polynomial kernel.  In _Advances in Neural Information Processing Systems(NIPS)_ , pages 2258–2266, 2014.
  * Avron et al. [2016] H. Avron, K. L. Clarkson, and D. P. Woodruff.  Sharper bounds for regression and low-rank approximation with regularization.  _CoRR_ , abs/1611.03225, 2016.
  * Carter and Wegman [1979] L. Carter and M. N. Wegman.  Universal classes of hash functions.  _J. Comput. Syst. Sci._ , 18(2):143–154, 1979\.
  * Charikar et al. [2004] M. Charikar, K. Chen, and M. Farach-Colton.  Finding frequent items in data streams.  _Theor. Comput. Sci._ , 312(1):3–15, 2004.
  * Chen and Plemmons [2010] D. Chen and R. J. Plemmons.  Nonnegativity constraints in numerical analysis.  In _The birth of numerical analysis_ , pages 109–139. World Sci. Publ., Hackensack, NJ, 2010.
  * Clarkson et al. [2013] K. Clarkson, P. Drineas, M. Magdon-Ismail, M. Mahoney, X. Meng, and D. P. Woodruff.  The fast Cauchy transform and faster robust linear regression.  In _SODA_ , 2013.
  * Clarkson [2005] K. L. Clarkson.  Subgradient and sampling algorithms for $\ell_{1}$ regression.  In _Proceedings of the sixteenth annual ACM-SIAM symposium on Discrete algorithms (SODA)_ , pages 257–266, 2005.
  * Clarkson and Woodruff [2013] K. L. Clarkson and D. P. Woodruff.  Low rank approximation and regression in input sparsity time.  In _Symposium on Theory of Computing Conference, STOC’13, Palo Alto, CA, USA, June 1-4, 2013_ , pages 81–90. https://arxiv.org/pdf/1207.6365, 2013.
  * Dasgupta et al. [2009] A. Dasgupta, P. Drineas, B. Harb, R. Kumar, and M. W. Mahoney.  Sampling algorithms and coresets for $\ell_{p}$ regression.  _SIAM J. Comput._ , 38(5):2060–2078, 2009.
  * Drineas et al. [2011] P. Drineas, M. Magdon-Ismail, M. W. Mahoney, and D. P. Woodruff.  Fast approximation of matrix coherence and statistical leverage.  _CoRR_ , abs/1109.3843, 2011.
  * Eilers and Marx [2006] P. H. Eilers and B. D. Marx.  Multidimensional density smoothing with p-splines.  In _Proceedings of the 21st international workshop on statistical modelling_ , 2006.
  * Eilers and Marx [1996] P. H. C. Eilers and B. D. Marx.  Flexible smoothing with $B$-splines and penalties.  _Statist. Sci._ , 11(2):89–121, 1996.
  * Eilers et al. [2015] P. H. C. Eilers, B. D. Marx, and M. Durbán.  Twenty years of P-splines.  _SORT_ , 39(2):149–186, 2015.  ISSN 1696-2281.
  * Fausett and Fulton [1994] D. W. Fausett and C. T. Fulton.  Large least squares problems involving Kronecker products.  _SIAM J. Matrix Anal. Appl._ , 15(1):219–227, 1994.
  * Golub and Van Loan [2013] G. H. Golub and C. F. Van Loan.  _Matrix computations_.  Johns Hopkins Studies in the Mathematical Sciences. Johns Hopkins University Press, Baltimore, MD, 2013.
  * Gurobi Optimization [2016] I. Gurobi Optimization.  Gurobi optimizer reference manual, 2016.  URL http://www.gurobi.com.
  * Hansen [1994] P. C. Hansen.  Regularization tools: A matlab package for analysis and solution of discrete ill-posed problems.  _Numerical Algorithms_ , 6(1):1–35, 1994.
  * Indyk [2006] P. Indyk.  Stable distributions, pseudorandom generators, embeddings, and data stream computation.  _J. ACM_ , 53(3):307–323, 2006.
  * Kannan and Vempala [2009] R. Kannan and S. Vempala.  Spectral algorithms.  _Foundations and Trends in Theoretical Computer Science_ , 4(3-4):157–288, 2009.
  * Liang et al. [2014] Y. Liang, M.-F. F. Balcan, V. Kanchanapally, and D. Woodruff.  Improved distributed principal component analysis.  In _Advances in Neural Information Processing Systems_ , pages 3113–3121, 2014.
  * Mahoney [2011] M. W. Mahoney.  Randomized algorithms for matrices and data.  _Foundations and Trends in Machine Learning_ , 3(2):123–224, 2011.
  * Meng and Mahoney [2013] X. Meng and M. W. Mahoney.  Low-distortion subspace embeddings in input-sparsity time and applications to robust linear regression.  In _Proceedings of the forty-fifth annual ACM symposium on Theory of computing_ , pages 91–100. ACM, https://arxiv.org/pdf/1210.3135, 2013.
  * Nelson and Nguyên [2013] J. Nelson and H. L. Nguyên.  Osnap: Faster numerical linear algebra algorithms via sparser subspace embeddings.  In _2013 IEEE 54th Annual Symposium on Foundations of Computer Science (FOCS)_ , pages 117–126. IEEE, https://arxiv.org/pdf/1211.1002, 2013\.
  * Oh and Yun [2005] S. Oh, S. Kwon and J. Yun.  A method for structured linear total least norm on blind deconvolution problem.  _Applied Mathematics and Computing_ , 19:151–164, 2005\.
  * Pagh [2013] R. Pagh.  Compressed matrix multiplication.  _ACM Trans. Comput. Theory_ , 5(3):9:1–9:17, 2013\.
  * Patrascu and Thorup [2012] M. Patrascu and M. Thorup.  The power of simple tabulation hashing.  _J. ACM_ , 59(3):14, 2012.
  * Pham and Pagh [2013] N. Pham and R. Pagh.  Fast and scalable polynomial kernels via explicit feature maps.  In _Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining(KDD)_ , pages 239–247. ACM, 2013.
  * Rousseeuw and Leroy [2005] P. J. Rousseeuw and A. M. Leroy.  _Robust regression and outlier detection_ , volume 589.  John wiley & sons, 2005.
  * Sohler and Woodruff [2011] C. Sohler and D. P. Woodruff.  Subspace embeddings for the $\ell_{1}$-norm with applications.  In _Proceedings of the forty-third annual ACM symposium on Theory of computing (STOC)_ , pages 755–764. ACM, 2011.
  * Song et al. [2017a] Z. Song, D. P. Woodruff, and P. Zhong.  Low rank approximation with entrywise $\ell_{1}$-norm error.  In _Proceedings of the 49th Annual Symposium on the Theory of Computing (STOC)_. ACM, https://arxiv.org/pdf/1611.00898, 2017a.
  * Song et al. [2017b] Z. Song, D. P. Woodruff, and P. Zhong.  Relative error tensor low rank approximation.  _arXiv preprint arXiv:1704.08246_ , 2017b.
  * Van Loan [1992] C. Van Loan.  _Computational frameworks for the fast Fourier transform_ , volume 10 of _Frontiers in Applied Mathematics_.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1992.
  * Van Loan and Pitsianis [1993] C. F. Van Loan and N. Pitsianis.  Approximation with Kronecker products.  In _Linear algebra for large scale and real-time applications (Leuven, 1992)_ , volume 232 of _NATO Adv. Sci. Inst. Ser. E Appl. Sci._ , pages 293–314. Kluwer Acad. Publ., Dordrecht, 1993.
  * Woodruff [2014] D. P. Woodruff.  Sketching as a tool for numerical linear algebra.  _Foundations and Trends in Theoretical Computer Science_ , 10(1-2):1–157, 2014.
  * Zheng et al. [2012] Y. Zheng, G. Liu, S. Sugimoto, S. Yan, and M. Okutomi.  Practical low-rank matrix approximation under robust l 1-norm.  In _Computer Vision and Pattern Recognition (CVPR), 2012 IEEE Conference on_ , pages 1410–1417. IEEE, 2012.

##  Appendix A Background: CountSketch and TensorSketch

We start by describing the CountSketch transform Charikar et al. [2004]. Let
$m$ be the target dimension. When applied to $n$-dimensional vectors, the
transform is specified by a $2$-wise independent hash function
$h:[n]\rightarrow[m]$ and a $2$-wise independent sign function
$s:[n]\rightarrow\\{-1,+1\\}$. When applied to $v$, the value at coordinate
$i$ of the output, $i=1,2,\ldots,m$ is $\sum_{j\mid h(j)=i}s(j)v_{j}$. Note
that CountSketch can be represented as an $m\times n$ matrix in which the
$j$-th column contains a single non-zero entry $s(j)$ in the $h(j)$-th row.

We now describe the TensorSketch transform Pagh [2013]. Suppose we are given
points $v_{i}\in\mathbb{R}^{n_{i}}$, where $i=1,\ldots,q$ and so
$\phi(v_{1},\ldots,v_{q})=v_{1}\otimes v_{2}\otimes\cdots\otimes
v_{q}\in\mathbb{R}^{n_{1}n_{2}\cdots n_{q}}$, and the target dimension is
again $m$. The transform is specified using $q$ $3$-wise independent hash
functions $h_{i}:[n_{i}]\rightarrow[m]$, and $q$ $4$-wise independent sign
functions $s_{i}:[n_{i}]\rightarrow\\{+1,-1\\}$, where $i=1,\ldots,q$.
TensorSketch applied to $v_{1},\ldots\otimes v_{q}$ is then CountSketch
applied to $\phi(v_{1},\ldots,v_{q})$ with hash function $H:[n_{1}n_{2}\cdots
n_{q}]\rightarrow[m]$ and sign function $S:[n_{1}n_{2}\cdots
n_{q}]\rightarrow\\{+1,-1\\}$ defined as follows:

$H(i_{1},\ldots,i_{q})=h_{1}(i_{1})+h_{2}(i_{2})+\cdots+h_{q}(i_{q})\bmod m,$

and

$S(i_{1},\ldots,i_{q})=s_{1}(i_{1})\cdot s_{2}(i_{2})\cdots s_{q}(i_{q}),$

where $i_{j}\in[n_{j}]$. It is well-known that if $H$ is constructed this way,
then it is $3$-wise independent Carter and Wegman [1979], Patrascu and Thorup
[2012]. Unlike the work of Pham and Pagh Pham and Pagh [2013], which only used
that $H$ was $2$-wise independent, our analysis needs this stronger property
of $H$.

The TensorSketch transform can be applied to $v_{1},\ldots,v_{q}$ without
computing $\phi(v_{1},\ldots,v_{q})$ as follows. Let
$v_{j}=(v_{j_{\ell}})\in\mathbb{R}^{n_{j}}$. First, compute the polynomials

$p_{\ell}(x)=\sum_{i=0}^{B-1}x^{i}\sum_{j_{\ell}\mid
h_{\ell}(j_{\ell})=i}v_{j_{\ell}}\cdot s_{\ell}(j_{\ell}),$

for $\ell=1,2,\ldots,q$. A calculation Pagh [2013] shows

$\prod_{\ell=1}^{q}p_{\ell}(x)\bmod(x^{B}-1)=\sum_{i=0}^{B-1}x^{i}\sum_{(j_{1},\ldots,j_{q})\mid
H(j_{1},\ldots,j_{q})=i}v_{j_{1}}\cdots v_{j_{q}}S(j_{1},\ldots,j_{q}),$

that is, the coefficients of the product of the $q$ polynomials $\bmod\
(x^{m}-1)$ form the value of TensorSketch($v_{1},\ldots,v_{q}$). Pagh observed
that this product of polynomials can be computed in $O(qm\log m)$ time using
the Fast Fourier Transform. As it takes
$O(q\max(\operatorname{\mathtt{nnz}}(v_{i})))$ time to form the $q$
polynomials, the overall time to compute TensorSketch$(v)$ is
$O(q(\max(\operatorname{\mathtt{nnz}}(v_{i}))+m\log m))$.

##  Appendix B TensorSketch is an Oblivious Subspace Embedding (OSE)

Let $S$ be the $m\times(n_{1}n_{2}\cdots n_{q})$ matrix such that TensorSketch
$(v_{1},\ldots,v_{q})$ is $S\cdot\phi(v_{1},\ldots,v_{q})$ for a randomly
selected TensorSketch. Notice that $S$ is a random matrix. In the rest of the
paper, we refer to such a matrix as a TensorSketch matrix with an appropriate
number of rows, i.e., the number of hash buckets. We will show that $S$ is an
oblivious subspace embedding for subspaces in ${\mathbb{R}}^{n_{1}n_{2}\cdots
n_{q}}$ for appropriate values of $m$. Notice that $S$ has exactly one non-
zero entry per column. The index of the non-zero in the column
$(i_{1},\ldots,i_{q})$ is
$H(i_{1},\ldots,i_{q})=\sum_{j=1}^{q}h_{j}(i_{j})\bmod{m}$. Let $\delta_{a,b}$
be the indicator random variable of whether $S_{a,b}$ is non-zero. The sign of
the non-zero entry in column $(i_{1},\ldots,i_{q})$ is
$S(i_{1},\ldots,i_{q})=\prod_{j=1}^{q}s_{j}(i_{j})$. We show that the
embedding matrix $S$ of TensorSketch can be used to approximate matrix product
and is an oblivious subspace embedding (OSE).

######  Theorem B.1.

Let $S$ be the $m\times(n_{1}n_{2}\cdots n_{q})$ matrix such that

$\mbox{ {\sc TensorSketch}}(v_{1},\ldots,v_{q})$

is $S\cdot\phi(v_{1},\ldots,v_{q})$ for a randomly selected TensorSketch. The
matrix $S$ satisfies the following two properties.

  1. 1.

($\mathrm{Approximate~{}Matrix~{}Product:}$) Let $A$ and $B$ be matrices with
$n_{1}n_{2}\cdots n_{q}$ rows. For $m\geq(2+3^{q})/(\epsilon^{2}\delta)$, we
have

$\Pr_{S}\left[\|A^{\top}S^{\top}SB-A^{\top}B\|_{F}^{2}\leq\epsilon^{2}\|A\|_{F}^{2}\|B\|_{F}^{2}\right]\geq
1-\delta.$

  2. 2.

($\mathrm{Subspace~{}Embedding:}$) Consider a fixed $k$-dimensional subspace
$V$. If $m\geq k^{2}(2+3^{q})/(\epsilon^{2}\delta)$, then with probability at
least $1-\delta$, $\|Sx\|=(1\pm\epsilon)\|x\|$ simultaneously for all $x\in
V$.

We establish the theorem via two lemmas as in Avron et al. [2016]. The first
lemma proves the approximate matrix product property via a careful second
moment analysis.

######  Lemma B.2 (Approximate matrix product).

Let $A$ and $B$ be matrices with $n_{1}n_{2}\cdots n_{q}$ rows. For
$m\geq(2+3^{q})/(\epsilon^{2}\delta)$, we have

$\displaystyle\Pr_{S}\left[\|A^{\top}S^{\top}SB-A^{\top}B\|_{F}^{2}\leq\epsilon^{2}\|A\|_{F}^{2}\|B\|_{F}^{2}\right]\geq
1-\delta.$

###### Proof.

The proof follows that in Avron et al. [2016]. Let $C=A^{\top}S^{\top}SB$. We
have

$C_{u,u^{\prime}}=\sum_{t=1}^{m}\sum_{i,j\in{[n_{1}n_{2}\cdots
n_{q}]}}S(i)S(j)\delta_{t,i}\delta_{t,j}A_{i,u}B_{j,u^{\prime}}=\sum_{t=1}^{m}\sum_{i\neq
j\in{[n_{1}n_{2}\cdots
n_{q}]}}S(i)S(j)\delta_{t,i}\delta_{t,i}A_{i,u}B_{j,u^{\prime}}+(A^{\top}B)_{u,u^{\prime}}$

Thus, ${\bf{E}}[C_{u,u^{\prime}}]=(A^{\top}B)_{u,u^{\prime}}$.

Next, we analyze ${\bf{E}}[((C-A^{\top}B)_{u,u^{\prime}})^{2}]$. We have

$\displaystyle((C-A^{\top}B)_{u,u^{\prime}})^{2}=\sum_{t_{1},t_{2}=1}^{m}\sum_{i_{1}\neq
j_{1},i_{2}\neq j_{2}\in{[n_{1}n_{2}\cdots n_{q}]}}$ $\displaystyle
S(i_{1})S(i_{2})S(j_{1})S(j_{2})\cdot\delta_{t_{1},i_{1}}\delta_{t_{1},j_{1}}\delta_{t_{2},i_{2}}\delta_{t_{2},j_{2}}\cdot
A_{i_{1},u}A_{i_{2},u}B_{j_{1},u^{\prime}}B_{j_{2},u^{\prime}}$

For a term in the summation on the right hand side to have a non-zero
expectation, it must be the case that
${\bf{E}}[S(i_{1})S(i_{2})S(j_{1})S(j_{2})]\neq 0$. Note that
$S(i_{1})S(i_{2})S(j_{1})S(j_{2})$ is a product of random signs (possibly with
multiplicities) where the random signs in different coordinates in
$\\{1,\ldots,q\\}$ are independent and they are 4-wise independent within each
coordinate. Thus, ${\bf{E}}[S(i_{1})S(i_{2})S(j_{1})S(j_{2})]$ is either $1$
or $0$. For the expectation to be $1$, all random signs must appear with even
multiplicities. In other words, in each of the $q$ coordinates, the 4
coordinates of $i_{1},i_{2},j_{1},j_{2}$ must be the same number appearing 4
times or 2 distinct numbers, each appearing twice. All the subsequent claims
in the proof regarding $i_{1},i_{2},j_{1},j_{2}$ agreeing on some coordinates
follow from this property.

Let $S_{1}$ be the set of coordinates where $i_{1}$ and $i_{2}$ agree. Note
that $j_{1}$ and $j_{2}$ must also agree in all coordinates in $S_{1}$ by the
above argument. Let $S_{2}\subset[q]\setminus S_{1}$ be the coordinates among
the remaining where $i_{1}$ and $j_{1}$ agree. Finally, let
$S_{3}=[q]\setminus(S_{1}\cup S_{2})$. All coordinates in $S_{3}$ of $i_{1}$
and $j_{2}$ must agree. Similarly as before, note that $i_{2}$ and $j_{2}$
agree on all coordinates in $S_{2}$ and $i_{2}$ and $j_{1}$ agree on all
coordinates in $S_{3}$. We can rewrite
$i_{1}=(a,b,c),i_{2}=(a,e,f),j_{1}=(g,b,f),j_{2}=(g,e,c)$ where
$a=(a_{\ell}),g=(g_{\ell})$ with $\ell\in S_{1}$, $b=(b_{\ell}),e=(e_{\ell})$
with $\ell\in S_{2}$ and $c=(c_{\ell}),f=(f_{\ell})$ with $\ell\in S_{3}$.

First we show that the contribution of the terms where $i_{1}=i_{2}$ or
$i_{1}=j_{2}$ is bounded by
$\frac{2\|A_{u}\|_{2}^{2}\|B_{u^{\prime}}\|_{2}^{2}}{m}$, where $A_{u}$ is the
$u$th column of $A$ and $B_{u^{\prime}}$ is the $u^{\prime}$th column of $B$.
Indeed, consider the case $i_{1}=i_{2}$. As observed before, we must have
$j_{1}=j_{2}$ to get a non-zero contribution. Note that if $t_{1}\neq t_{2}$,
we always have $\delta_{t_{1},i_{1}}\delta_{t_{2},i_{2}}=0$ as $H(i_{1})$
cannot be equal to both $t_{1}$ and $t_{2}$. Thus, for fixed
$i_{1}=i_{2},j_{1}=j_{2}$,

$\displaystyle~{}{\bf{E}}\left[\sum_{t_{1},t_{2}=1}^{m}S(i_{1})S(i_{2})S(j_{1})S(j_{2})\cdot\delta_{t_{1},i_{1}}\delta_{t_{1},j_{1}}\delta_{t_{2},i_{2}}\delta_{t_{2},j_{2}}\cdot
A_{i_{1},u}A_{i_{2},u}B_{j_{1},u^{\prime}}B_{j_{2},u^{\prime}}\right]$
$\displaystyle=$
$\displaystyle~{}{\bf{E}}\left[\sum_{t_{1}=1}^{m}\delta_{i_{1},t_{1}}^{2}\delta_{j_{1},t_{1}}^{2}A_{i_{1},u}^{2}B_{j_{1},u^{\prime}}^{2}\right]$
$\displaystyle=$
$\displaystyle~{}\frac{A_{i_{1},u}^{2}B_{j_{1},u^{\prime}}^{2}}{m}$

Summing over all possible values of $i_{1},j_{1}$, we get the desired bound of
$\frac{\|A_{u}\|_{2}^{2}\|B_{u^{\prime}}\|_{2}^{2}}{m}$. The case
$i_{1}=j_{2}$ is analogous.

Next we compute the contribution of the terms where $i_{1}\neq
i_{2},j_{1},j_{2}$ i.e., there are at least 3 distinct numbers among
$i_{1},i_{2},j_{1},j_{2}$. Notice that
${\bf{E}}[\delta_{t_{1},i_{1}}\delta_{t_{1},j_{1}}\delta_{t_{2},i_{2}}\delta_{t_{2},j_{2}}]\leq\frac{1}{m^{3}}$
because the $\delta_{t,i}$’s are 3-wise independent. For fixed
$i_{1},j_{1},i_{2},j_{2}$, there are $m^{2}$ choices of $t_{1},t_{2}$ so the
total contribution to the expectation from terms with the same
$i_{1},j_{1},i_{2},j_{2}$ is bounded by
$m^{2}\cdot\frac{1}{m^{3}}\cdot|A_{i_{1},u}A_{i_{2},u}B_{j_{1},u^{\prime}}B_{j_{2},u^{\prime}}|=\frac{1}{m}|A_{i_{1},u}A_{i_{2},u}B_{j_{1},u^{\prime}}B_{j_{2},u^{\prime}}|$.

Therefore,

$\displaystyle~{}{\bf{E}}[((C-A^{\top}B)_{u,u^{\prime}})^{2}]$
$\displaystyle\leq$
$\displaystyle~{}\frac{2\|A_{u}\|_{2}^{2}\|B_{u^{\prime}}\|_{2}^{2}}{m}+\frac{1}{m}\sum_{\textnormal{partition
}S_{1},S_{2},S_{3}}\sum_{a,g,b,e,c,f}|A_{(a,b,c),u}B_{(g,b,f),u^{\prime}}A_{(a,e,f),u}B_{(g,e,c),u^{\prime}}|$
$\displaystyle\leq$
$\displaystyle~{}\frac{2\|A_{u}\|_{2}^{2}\|B_{u^{\prime}}\|_{2}^{2}}{m}+\frac{3^{q}}{m}\sum_{a,b,c,g,e,f}|A_{(a,b,c),u}B_{(g,b,f),u^{\prime}}A_{(a,e,f),u}B_{(g,e,c),u^{\prime}}|$
$\displaystyle\leq$
$\displaystyle~{}\frac{2\|A_{u}\|_{2}^{2}\|B_{u^{\prime}}\|_{2}^{2}}{m}+\frac{3^{q}}{m}\sum_{g,e,f}\biggl{(}\sum_{a,b,c}A_{(a,b,c),u}^{2}\biggr{)}^{1/2}\biggl{(}\sum_{a,b,c}B_{(g,b,f),u^{\prime}}^{2}A_{(a,e,f),u}^{2}B_{(g,e,c),u^{\prime}}^{2}\biggr{)}^{1/2}$
$\displaystyle=$
$\displaystyle~{}\frac{2\|A_{u}\|_{2}^{2}\|B_{u^{\prime}}\|_{2}^{2}}{m}+\frac{3^{q}\|A_{u}\|}{m}\sum_{g,e,f}\biggl{(}\sum_{b}B_{(g,b,f),u^{\prime}}^{2}\biggr{)}^{1/2}\biggl{(}\sum_{a,c}A_{(a,e,f),u}^{2}B_{(g,e,c),u^{\prime}}^{2}\biggr{)}^{1/2}$
$\displaystyle\leq$
$\displaystyle~{}\frac{2\|A_{u}\|_{2}^{2}\|B_{u^{\prime}}\|_{2}^{2}}{m}+\frac{3^{q}\|A_{u}\|}{m}\sum_{e}\biggl{(}\sum_{b,g,f}B_{(g,b,f),u^{\prime}}^{2}\biggr{)}^{1/2}\biggl{(}\sum_{a,c,g,f}A_{(a,e,f),u}^{2}B_{(g,e,c),u^{\prime}}^{2}\biggr{)}^{1/2}$
$\displaystyle=$
$\displaystyle~{}\frac{2\|A_{u}\|_{2}^{2}\|B_{u^{\prime}}\|_{2}^{2}}{m}+\frac{3^{q}\|A_{u}\|\cdot\|B_{u^{\prime}}\|}{m}\sum_{e}\biggl{(}\sum_{a,f}A_{(a,e,f),u}^{2}\biggr{)}^{1/2}\biggl{(}\sum_{g,c}B_{(g,e,c),u^{\prime}}^{2}\biggr{)}^{1/2}$
$\displaystyle\leq$
$\displaystyle~{}\frac{2\|A_{u}\|_{2}^{2}\|B_{u^{\prime}}\|_{2}^{2}}{m}+\frac{3^{q}\|A_{u}\|\cdot\|B_{u^{\prime}}\|}{m}\biggl{(}\sum_{a,e,f}A_{(a,e,f),u}^{2}\biggr{)}^{1/2}\biggl{(}\sum_{g,e,c}B_{(g,e,c),u^{\prime}}^{2}\biggr{)}^{1/2}$
$\displaystyle=$
$\displaystyle~{}\frac{(2+3^{q})\|A_{u}\|_{2}^{2}\|B_{u^{\prime}}\|_{2}^{2}}{m},$

where the second inequality follows from the fact that there are at most
$3^{q}$ partitions of $[q]$ into 3 sets. The other inequalities are from
Cauchy-Schwarz.

Combining the above bounds, we have
${\bf{E}}[((C-A^{\top}B)_{u,u^{\prime}})^{2}]\leq\frac{(2+3^{q})\|A_{u}\|_{2}^{2}\|B_{u^{\prime}}\|_{2}^{2}}{m}.$
For $m\geq(2+3^{q})/(\epsilon^{2}\delta)$, by the Markov inequality,
$\|A^{\top}S^{\top}SB-A^{\top}B\|_{F}^{2}\leq\epsilon^{2}\|A\|_{F}^{2}\|B\|_{F}^{2}$
with probability $1-\delta$.∎∎

The second lemma proves that the subspace embedding property follows from the
approximate matrix product property.

######  Lemma B.3 (Oblivious subspace embeddings).

Consider a fixed $k$-dimensional subspace
$V\subset{\mathbb{R}}^{n_{1}n_{2}\cdots n_{q}}$. If $m\geq
k^{2}(2+3^{q})/(\epsilon^{2}\delta)$, then with probability at least
$1-\delta$, $\|Sx\|_{2}=(1\pm\epsilon)\|x\|_{2}$ simultaneously for all $x\in
V$.

###### Proof.

Let $B$ be a $(n_{1}n_{2}\cdots n_{q})\times k$ matrix whose columns form an
orthonormal basis of $V$. Thus, we have $B^{\top}B=I_{k}$ and
$\|B\|_{F}^{2}=k$. The condition that $\|Sx\|_{2}=(1\pm\epsilon)\|x\|_{2}$
simultaneously for all $x\in V$ is equivalent to the condition that the
singular values of $SB$ are bounded by $1\pm\epsilon$. By Lemma B.2, for
$m\geq(2+3^{q})/((\epsilon/k)^{2}\delta)$, with probability at least
$1-\delta$, we have

$\|B^{\top}S^{\top}SB-B^{\top}B\|_{F}^{2}\leq(\epsilon/k)^{2}\|B\|_{F}^{4}=\epsilon^{2}$

Thus, we have $\|B^{\top}S^{\top}SB-I_{k}\|_{2}\leq\|B^{\top}S^{\top}SB-
I_{k}\|_{F}\leq\epsilon$. In other words, the squared singular values of $SB$
are bounded by $1\pm\epsilon$, implying that the singular values of $SB$ are
also bounded by $1\pm\epsilon$. Note that $\|A\|_{2}$ for a matrix $A$ denotes
its operator norm.∎∎

##  Appendix C Missing Proofs

###  C.1 Proofs for Tensor Product Least Square Regression

######  Theorem 3.1.

(Tensor regression) Suppose $\widetilde{x}$ is the output of Algorithm 1 with
tensorsketch $S\in\mathbb{R}^{m\times n}$, where $m=8(d_{1}d_{2}\cdots
d_{q}+1)^{2}(2+3^{q})/(\epsilon^{2}\delta)$. Then the following approximation
$\|(A_{1}\otimes A_{2}\otimes\cdots\otimes
A_{q})\widetilde{x}-b\|_{2}\leq(1+\epsilon)\operatorname{OPT},$ holds with
probability at least $1-\delta$.

###### Proof.

It is easy to see that

$\|(A_{1}\otimes A_{2}\otimes\cdots\otimes
A_{q})x-b\|_{2}=\left\|\begin{bmatrix}(A_{1}\otimes A_{2}\otimes\cdots\otimes
A_{q})&b\end{bmatrix}\begin{bmatrix}x\cr-1\end{bmatrix}\right\|_{2},$

and identifying

$y=\begin{bmatrix}(A_{1}\otimes A_{2}\otimes\cdots\otimes
A_{q})&b\end{bmatrix}\begin{bmatrix}x\cr-1\end{bmatrix}\in{\mathbb{R}}^{n_{1}n_{2}\cdots
n_{q}}$

and $y$ is a vector of a subspace $V\subset{\mathbb{R}}^{n_{1}n_{2}\cdots
n_{q}}$ with dimension at most $d_{1}d_{2}\cdots d_{q}+1$, we can use Lemma
B.3 to conclude that

$\displaystyle\Pr\left[\left|\left\|Sy\right\|_{2}-\|y\|_{2}\right|\leq\epsilon\|y\|_{2}\right]\geq
1-\delta$

when $m=(d_{1}d_{2}\cdots d_{q}+1)^{2}(2+3^{q})/(\epsilon^{2}\delta)$.

Thus we have

$\displaystyle\|(A_{1}\otimes A_{2}\otimes\cdots\otimes
A_{q})\widetilde{x}-b\|_{2}\leq\frac{1}{1-\epsilon}\|S(A_{1}\otimes
A_{2}\otimes\cdots\otimes A_{q})\widetilde{x}-Sb\|_{2}$

and

$\displaystyle\|S(A_{1}\otimes A_{2}\otimes\cdots\otimes
A_{q}){x}-Sb\|_{2}\leq(1+\epsilon)\|(A_{1}\otimes A_{2}\otimes\cdots\otimes
A_{q}){x}-b\|_{2}$

hold with probability at least $1-\delta$. Then using a union bound, we have

$\displaystyle~{}\|(A_{1}\otimes A_{2}\otimes\cdots\otimes
A_{q})\widetilde{x}-b\|_{2}$ $\displaystyle\leq$
$\displaystyle~{}\frac{1}{1-\epsilon}\|S(A_{1}\otimes
A_{2}\otimes\cdots\otimes A_{q})\widetilde{x}-Sb\|_{2}$ $\displaystyle\leq$
$\displaystyle~{}\frac{1}{1-\epsilon}\|S(A_{1}\otimes
A_{2}\otimes\cdots\otimes A_{q}){x}-Sb\|_{2}$ $\displaystyle\leq$
$\displaystyle~{}\frac{1+\epsilon}{1-\epsilon}\|(A_{1}\otimes
A_{2}\otimes\cdots\otimes A_{q}){x}-b\|_{2}$

holds with probability at least $1-2\delta$. ∎

######  Corollary 3.2.

(Sketch for tensor nonnegative regression) Suppose $\tilde{x}=\min_{x\geq
0}\|S\mathcal{A}x-Sb\|_{2}$ with tensorsketch $S\in\mathbb{R}^{m\times n}$,
where $m=8(d_{1}d_{2}\cdots d_{q}+1)^{2}(2+3^{q})/(\epsilon^{2}\delta)$. Then
the following approximation $\|(A_{1}\otimes A_{2}\otimes\cdots\otimes
A_{q})\widetilde{x}-b\|_{2}\leq(1+\epsilon)\operatorname{OPT}$ holds with
probability at least $1-\delta$, where $\operatorname{OPT}=\min_{x\geq
0}\left\|(A_{1}\otimes A_{2}\otimes\cdots\otimes A_{q}){x}-b\right\|_{2}$.

###### Proof.

The proof of Theorem. 3.2 is similar to the proof of theorem 3.1. Denote
$\tilde{x}=\min_{x\geq 0}\|S\mathcal{A}x-Sb\|_{2}$ and $x^{*}=\min_{x\geq
0}\|\mathcal{A}x-b\|_{2}$. Using Lemma. B.3, we have:

$\displaystyle\|\mathcal{A}\tilde{x}-b\|_{2}\leq\frac{1}{1-\epsilon}\|S\mathcal{A}\tilde{x}-Sb\|_{2},$
(6)

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

$\displaystyle\|S\mathcal{A}{x}^{*}-Sb\|_{2}\leq(1+\epsilon)\|\mathcal{A}{x}^{*}-b\|_{2},$
(7)

with probability at least $1-\delta$. Hence applying a union bound we have:

$\displaystyle~{}\|\mathcal{A}\tilde{x}-b\|_{2}$ (8) $\displaystyle\leq$
$\displaystyle~{}\frac{1}{1-\epsilon}\|S\mathcal{A}\tilde{x}-Sb\|_{2}$
$\displaystyle\leq$
$\displaystyle~{}\frac{1}{1-\epsilon}\|S\mathcal{A}{x}^{*}-Sb\|_{2}$
$\displaystyle\leq$
$\displaystyle~{}\frac{1+\epsilon}{1-\epsilon}\|\mathcal{A}{x}^{*}-b\|_{2},$
(9)

with probability at least $1-2\delta$. ∎

###  C.2 Proofs for P-Splines

######  Lemma 4.1.

Let $x^{*}\in{\mathbb{R}}^{d}$, $A\in{\mathbb{R}}^{n\times d}$ and
$b\in{\mathbb{R}}^{n}$ as above. Let $U_{1}\in{\mathbb{R}}^{n\times d}$ denote
the first $n$ rows of an orthogonal basis for $\left[\begin{smallmatrix}A\\\
\sqrt{\lambda}L\end{smallmatrix}\right]\in{\mathbb{R}}^{(n+p)\times d}$. Let
sketching matrix $S\in{\mathbb{R}}^{m\times n}$ have a distribution such that
with constant probability

$\displaystyle\mathrm{(\@slowromancap i@)}~{}\lVert
U_{1}^{\top}S^{\top}SU_{1}-U_{1}^{\top}U_{1}\rVert_{2}\leq 1/4,$

and

$\displaystyle\mathrm{(\@slowromancap
ii@)}~{}\|U_{1}^{\top}(S^{\top}S-I)(b-Ax^{*})\|_{2}\leq\sqrt{\epsilon\operatorname{OPT}/2}.$

Let $\widetilde{x}$ denote
$\operatorname*{argmin}_{x\in{\mathbb{R}}^{d}}\lVert
S(Ax-b)\rVert_{2}^{2}+\lambda\lVert Lx\rVert_{2}^{2}$. Then with probability
at least $9/10$,

$\displaystyle\lVert A\widetilde{x}-b\rVert_{2}^{2}+\lambda\lVert
L\widetilde{x}\rVert_{2}^{2}\leq(1+\epsilon)\operatorname{OPT}.$

###### Proof.

Let $\hat{A}\in{\mathbb{R}}^{(n+d)\times d}$ have orthonormal columns with
$\mathop{\rm range}\nolimits(\hat{A})=\mathop{\rm
range}\nolimits(\left[\begin{smallmatrix}A\\\
\sqrt{\lambda}L\end{smallmatrix}\right])$. (An explicit expression for one
such $\hat{A}$ is given below.) Let
$\hat{b}\equiv\left[\begin{smallmatrix}b\\\ 0_{d}\end{smallmatrix}\right]$. We
have

$\min_{y\in{\mathbb{R}}^{d}}\lVert\hat{A}y-\hat{b}\rVert_{2}$ (10)

equivalent to $\left\|b-Ax\right\|_{2}^{2}+\lambda\|Lx\|_{2}^{2}$, , in the
sense that for any $\hat{A}y\in\mathop{\rm range}\nolimits(\hat{A})$, there is
$x\in{\mathbb{R}}^{d}$ with $\hat{A}y=\left[\begin{smallmatrix}A\\\
\sqrt{\lambda}L\end{smallmatrix}\right]x$, so that
$\lVert\hat{A}y-\hat{b}\rVert_{2}^{2}=\lVert\left[\begin{smallmatrix}A\\\
\sqrt{\lambda}L\end{smallmatrix}\right]x-\hat{b}\rVert_{2}^{2}=\lVert
b-Ax\rVert_{2}^{2}+\lambda\lVert Lx\rVert_{2}^{2}$. Let
$y^{*}=\operatorname*{argmin}_{y\in{\mathbb{R}}^{d}}\lVert\hat{A}y-\hat{b}\rVert_{2}$,
so that $\hat{A}y^{*}=\left[\begin{smallmatrix}Ax^{*}\\\
\sqrt{\lambda}Lx^{*}\end{smallmatrix}\right]$.

Let $\hat{A}=\left[\begin{smallmatrix}U_{1}\\\ U_{2}\end{smallmatrix}\right]$,
where $U_{1}\in{\mathbb{R}}^{n\times d}$ and $U_{2}\in{\mathbb{R}}^{d\times
d}$, so that $U_{1}$ is as in the lemma statement.

We define $\hat{S}$ to be $\left[\begin{smallmatrix}S&0_{m\times d}\\\
0_{d\times n}&{I}_{d}\end{smallmatrix}\right]$ and $\hat{S}$ satisfies
Property (i@) and (ii@) of Lemma 4.1.

Using $\lVert U_{1}^{\top}S^{\top}SU_{1}-U_{1}^{\top}U_{1}\rVert_{2}\leq 1/4$,
with constant probability

$\lVert\hat{A}^{\top}\hat{S}^{\top}\hat{S}\hat{A}-{I}_{d}\rVert_{2}=\lVert
U_{1}^{\top}S^{\top}SU_{1}+U_{2}^{\top}U_{2}-{I}_{d}\rVert_{2}=\lVert
U_{1}^{\top}S^{\top}SU_{1}-U_{1}^{\top}U_{1}\rVert_{2}\leq 1/4.$ (11)

Using the normal equations for Eq. (10), we have

$0=\hat{A}^{\top}(\hat{b}-\hat{A}y^{*})=U_{1}^{\top}(b-Ax^{*})-\sqrt{\lambda}U_{2}^{\top}x^{*},$

and so

$\hat{A}^{\top}\hat{S}^{\top}\hat{S}(\hat{b}-\hat{A}y^{*})=U_{1}^{\top}S^{\top}S(b-Ax^{*})-\sqrt{\lambda}U_{2}^{\top}x^{*}=U_{1}^{\top}S^{\top}S(b-Ax^{*})-U_{1}^{\top}(b-Ax^{*}).$

Using Property (ii@) of Lemma 4.1, with constant probability

$\displaystyle~{}\lVert\hat{A}^{\top}\hat{S}^{\top}\hat{S}(\hat{b}-\hat{A}y^{*})\rVert_{2}$
$\displaystyle=$ $\displaystyle~{}\lVert
U_{1}^{\top}S^{\top}S(b-Ax^{*})-U_{1}^{\top}(b-Ax^{*})\rVert_{2}$
$\displaystyle\leq$ $\displaystyle~{}\sqrt{\epsilon\operatorname{OPT}/2}$
$\displaystyle=$
$\displaystyle~{}\sqrt{\epsilon/2}\lVert\hat{b}-\hat{A}y^{*}\rVert_{2}.$ (12)

It follows by a standard result from (11) and (12) that the solution
$\tilde{y}\equiv\operatorname*{argmin}_{y\in{\mathbb{R}}^{d}}\lVert\hat{S}(\hat{A}y-\hat{b})\rVert_{2}$
has
$\lVert\hat{A}\tilde{y}-\hat{b}\rVert_{2}\leq(1+\epsilon)\min_{y\in{\mathbb{R}}^{d}}\lVert\hat{A}y-\hat{b}\rVert_{2}$,
and therefore that $\tilde{x}$ satisfies the claim of the theorem.

For convenience we give the proof of the standard result: (11) implies that
$\hat{A}^{\top}\hat{S}^{\top}\hat{S}\hat{A}$ has smallest singular value at
least $3/4$. The normal equations for the unsketched and sketched problems are

$\hat{A}^{\top}(\hat{b}-\hat{A}y^{*})=0=\hat{A}^{\top}\hat{S}^{\top}\hat{S}(\hat{b}-\hat{A}\tilde{y}).$

The normal equations for the unsketched case imply
$\lVert\hat{A}\tilde{y}-\hat{b}\rVert_{2}^{2}=\lVert\hat{A}(\tilde{y}-y^{*})\rVert_{2}^{2}+\lVert\hat{b}-\hat{A}y^{*}\rVert_{2}^{2}$,
so it is enough to show that
$\lVert\hat{A}(\tilde{y}-y^{*})\rVert_{2}^{2}=\lVert\tilde{y}-y^{*}\rVert_{2}^{2}\leq\epsilon\operatorname{OPT}$.
We have

$\displaystyle(3/4)\lVert\tilde{y}-y^{*}\rVert_{2}$
$\displaystyle\leq\lVert\hat{A}^{\top}\hat{S}^{\top}\hat{S}\hat{A}(\tilde{y}-y^{*})\rVert_{2}$
by Eq. (11)
$\displaystyle=\lVert\hat{A}^{\top}\hat{S}^{\top}\hat{S}\hat{A}(\tilde{y}-y^{*})-\hat{A}^{\top}\hat{S}^{\top}\hat{S}(\hat{b}-\hat{A}\tilde{y})\rVert_{2}$
by Normal Equation
$\displaystyle=\lVert\hat{A}^{\top}\hat{S}^{\top}\hat{S}(\hat{b}-\hat{A}y^{*})\rVert_{2}$
$\displaystyle\leq\sqrt{\epsilon\operatorname{OPT}/2}$
$\displaystyle\text{by~{}Eq.~{}\eqref{eq prod2}},$

so that
$\lVert\tilde{y}-y^{*}\rVert_{2}^{2}\leq(4/3)^{2}\epsilon\operatorname{OPT}/2\leq\epsilon\operatorname{OPT}$.
The lemma follows.

∎

The following lemma computes the statistical dimension
$\operatorname{\mathtt{sd}}_{\lambda}(A,L)$ that will be used for computing
the number of rows of sketching matrix $S$.

######  Lemma C.1.

For $U_{1}$ as in Lemma 4.1,
${\|U_{1}\|}_{F}^{2}=\operatorname{\mathtt{sd}}_{\lambda}(A,L)=\sum_{i}1/(1+\lambda/\gamma_{i}^{2})+d-p$,
where $A$ has singular values $\sigma_{i}$. Also $\lVert
U_{1}\rVert_{2}=\max\\{1/\sqrt{1+\lambda/\gamma_{1}^{2}},\,1\\}$.

###### Proof.

Suppose we have the GSVD of $(A,\,L)$. Let

$D\equiv\begin{bmatrix}\Sigma^{\top}\Sigma+\lambda\Omega^{\top}\Omega&{0}_{p\times(n-p)}\cr{0}_{(n-p)\times
p}&I_{d-p}\end{bmatrix}^{-1/2}.$

Then

$\hat{A}=\left[\begin{smallmatrix}U\begin{bmatrix}\Sigma&{0}_{p\times(n-p)}\cr{0}_{(n-p)\times
p}&I_{d-p}\end{bmatrix}D\\\
\sqrt{\lambda}V\begin{bmatrix}\Omega&{0}_{p\times(n-p)}\end{bmatrix}D\end{smallmatrix}\right]$

has $\hat{A}^{\top}\hat{A}={I}_{d}$, and for given $x$, there is
$y=D^{-1}RQ^{\top}x$ with $\hat{A}y=\left[\begin{smallmatrix}A\\\
\sqrt{\lambda}L\end{smallmatrix}\right]x$. We have
${\|U_{1}\|}_{F}^{2}=\left\|{U\begin{bmatrix}\Sigma&{0}_{p\times(n-p)}\cr{0}_{(n-p)\times
p}&I_{d-p}\end{bmatrix}D}\right\|_{F}^{2}=\left\|{\begin{bmatrix}\Sigma&{0}_{p\times(n-p)}\cr{0}_{(n-p)\times
p}&I_{d-p}\end{bmatrix}D}\right\|_{F}^{2}=\sum_{i=1}^{p}1/(1+\lambda/\gamma_{i}^{2})+d-p$
as claimed. ∎

######  Theorem 4.3.

(P-Spline regression) There is a constant $K>0$ such that for $m\geq
K(\epsilon^{-1}\operatorname{\mathtt{sd}}_{\lambda}(A,L)+\operatorname{\mathtt{sd}}_{\lambda}(A,L)^{2})$
and $S\in{\mathbb{R}}^{m\times n}$ a sparse embedding matrix (e.g.,
Countsketch) with $SA$ computable in $O(\operatorname{\mathtt{nnz}}(A))$ time,
Property (i@) and (ii@) of Lemma 4.1 apply, and with constant probability the
corresponding
$\widetilde{x}=\operatorname*{argmin}_{x\in{\mathbb{R}}^{d}}\lVert
S(Ax-b)\rVert_{2}+\lambda\lVert Lx\rVert_{2}^{2}$ is an $\epsilon$-approximate
solution to $\min_{x\in{\mathbb{R}}^{d}}\lVert
b-Ax\rVert_{2}^{2}+\lambda\lVert Lx\rVert_{2}^{2}$.

###### Proof.

Recall that $\operatorname{\mathtt{sd}}_{\lambda}(A,L)={\|U_{1}\|}_{F}^{2}$.
Sparse embedding distributions satisfy the bound for approximate matrix
multiplication

${\|W^{\top}S^{\top}SH-W^{\top}H\|}_{F}\leq C{\|W\|}_{F}{\|H\|}_{F}/\sqrt{m},$

for a constant $C$ [Clarkson and Woodruff, 2013, Meng and Mahoney, 2013,
Nelson and Nguyên, 2013]; this is also true of OSE matrices. We set
$W=H=U_{1}$ and use $\lVert X\rVert_{2}\leq{\|X\|}_{F}$ for all $X$ and $m\geq
K{\|U_{1}\|}_{F}^{4}$ to obtain Property (i@) of Lemma 4.1, and set $W=U_{1}$,
$H=b-Ax^{*}$ and use $m\geq K{\|U_{1}\|}_{F}^{2}/\epsilon$ to obtain Property
(ii@) of Lemma 4.1. (Here the bound is slightly stronger than Property (ii@),
holding for $\lambda=0$.) With Property (i@) and Property (ii@), the claim for
$\tilde{x}$ from a sparse embedding follows using Lemma 4.1. ∎

###  C.3 Proofs for Tensor Product $\ell_{1}$ Regression

######  Lemma 5.3.

For any $p\geq 1$. Condition$({\cal A})$ computes ${\cal A}U/(d\gamma_{p})$
which is an $(\alpha,\beta\sqrt{3}d(tw)^{|1/p-1/2|},p)$-well-conditioned basis
of ${\cal A}$, with probability at least
$1-\prod_{i=1}^{q}(n_{i}/w_{i})\delta$.

###### Proof.

This lemma is similar to arguments in Clarkson et al. [2013], we simply adjust
notation and parameters for completeness. Applying Theorem 5.2, we have that
with probability at least $1-\prod_{i=1}^{q}(n_{i}/w_{i})\delta$, for all
$x\in\mathbb{R}^{r}$, if we consider $y={\cal A}x$ and write
$y^{\top}=[z_{1}^{\top},z_{2}^{\top},\ldots,z_{\prod_{i=1}^{q}n_{i}/w_{i}}]^{\top}$,
then for all $i\in[\prod_{i=1}^{q}n_{i}/w_{i}]$,

$\displaystyle\sqrt{\textstyle\frac{1}{2}}\|{z_{i}}\|_{2}\leq\|S_{i}z_{i}\|_{2}\leq\sqrt{\textstyle\frac{3}{2}}\|{z_{i}}\|_{2},$

where $S_{i}\in{\mathbb{R}}^{m_{i}\times\prod_{i=1}^{q}w_{i}}$. In the
following, suppose $m_{i}=t$. By relating the $2$-norm and the $p$-norm, for
$1\leq p\leq 2$, we have

$\|{S_{i}z_{i}}\|_{p}\leq t^{1/p-1/2}\|{Sz_{i}}\|_{2}\leq
t^{1/p-1/2}\sqrt{\textstyle\frac{3}{2}}\|{z_{i}}\|_{2}\leq
t^{1/p-1/2}\sqrt{\textstyle\frac{3}{2}}\|{z_{i}}\|_{p},$

and similarly,

$\|{S_{i}z_{i}}\|_{p}\geq\|{S_{i}z_{i}}\|_{2}\geq\sqrt{\textstyle\frac{1}{2}}\|{z_{i}}\|_{2}\geq\sqrt{\textstyle\frac{1}{2}}w^{1/2-1/p}\|{z_{i}}\|_{p},\,w=\prod_{j=1}^{q}w_{j}.$

If $p>2$, then

$\|{S_{i}z_{i}}\|_{p}\leq\|{S_{i}z_{i}}\|_{2}\leq\sqrt{\textstyle\frac{3}{2}}\|{z_{i}}\|_{2}\leq\sqrt{\textstyle\frac{3}{2}}w^{1/2-1/p}\|{z_{i}}\|_{p},$

and similarly,

$\|{S_{i}z_{i}}\|_{p}\geq t^{1/p-1/2}\|{S_{i}z_{i}}\|_{2}\geq
t^{1/p-1/2}\sqrt{\textstyle\frac{1}{2}}\|{z_{i}}\|_{2}\geq
t^{1/p-1/2}\sqrt{\textstyle\frac{1}{2}}\|{z_{i}}\|_{p}.$

Since $\|{{\cal A}x}\|_{p}^{p}=\|{y}\|_{p}^{p}=\sum_{i}\|{z_{i}}\|_{p}^{p}$
and $\|{S{\cal A}x}\|_{p}^{p}=\sum_{i}\|{S_{i}z_{i}}\|_{p}^{p}$, for
$p\in[1,2]$ we have with probability $1-\prod_{i=1}^{q}(n_{i}/w_{i})\delta$

$\sqrt{\textstyle\frac{1}{2}}w^{1/2-1/p}\|{{\cal A}x}\|_{p}\leq\|{S{\cal
A}x}\|_{p}\leq\sqrt{\textstyle\frac{3}{2}}t^{1/p-1/2}\|{{\cal A}x}\|_{p},$

and for $p\in[2,\infty)$ with probability
$1-\prod_{i=1}^{q}(n_{i}/w_{i})\delta$

$\sqrt{\textstyle\frac{1}{2}}t^{1/p-1/2}\|{{\cal A}x}\|_{p}\leq\|{S{\cal
A}x}\|_{p}\leq\sqrt{\textstyle\frac{3}{2}}w^{1/2-1/p}\|{{\cal A}x}\|_{p}.$

In either case,

$\|{{\cal A}x}\|_{p}\leq\gamma_{p}\|{S{\cal
A}x}\|_{p}\leq\sqrt{3}(tw)^{|1/p-1/2|}\|{{\cal A}x}\|_{p}.$ (13)

We have, from the definition of an $(\alpha,\beta,p)$-well-conditioned basis,
that

$\displaystyle\|S{\cal A}U\|_{p}\leq\alpha$ (14)

and for all $x\in\mathbb{R}^{d}$,

$\displaystyle\|x\|_{q}\leq\beta\|S{\cal A}Ux\|_{p}.$ (15)

Combining (13) and (14), we have that with probability at least
$1-\prod_{i=1}^{q}(n_{i}/w_{i})\delta$,

$\displaystyle\|{\cal A}U/(r\gamma_{p})\|_{p}\leq\sum_{i}\|{\cal
A}U_{i}/r\gamma_{p}\|_{p}\leq\sum_{i}\|S{\cal A}U_{i}/r\|_{p}\leq\alpha.$

Combining (13) and (15), we have that with probability at least
$1-\prod_{i=1}^{q}(n_{i}/w_{i})\delta$, for all $x\in\mathbb{R}^{r}$,

$\displaystyle\|x\|_{q}\leq\beta\|S{\cal
A}Ux\|_{p}\leq\beta\sqrt{3}r(tw)^{|1/p-1/2|}\|{\cal
A}U\frac{1}{r\gamma_{p}}x\|_{p}.$

Hence ${\cal A}U/(r\gamma_{p})$ is an
$(\alpha,\beta\sqrt{3}r(tw)^{|1/p-1/2|},p)$-well-conditioned basis. ∎

######  Theorem 5.4.

(Main result) Given $\epsilon\in(0,1)$, ${\cal A}\in{\mathbb{R}}^{n\times d}$
and $b\in{\mathbb{R}}^{n}$, Alg. 3 computes $\widehat{x}$ such that with
probability at least $1/2$, $\|{\cal
A}\hat{x}-b\|_{1}\leq(1+\epsilon)\min_{x\in{\mathbb{R}}^{d}}\|{\cal
A}x-b\|_{1}$. For the special case when $q=2$, $n_{1}=n_{2}$, the algorithm’s
running time is
$O({n_{1}}^{3/2}\operatorname{poly}(\prod_{i=1}^{2}d_{i}/\epsilon))$.

###### Proof.

For notational simplicity, let us denote
$n_{[q_{1}]}=\prod_{i=1}^{q_{1}}n_{i}$,
$n_{[q]\backslash[q_{1}]}=\prod_{i=q_{1}+1}^{q}n_{1}$,
$d_{[q_{1}]}=\prod_{i=1}^{q_{1}}d_{i}$, and
$d_{[q]\backslash[q_{1}]}=\prod_{i=q_{1}+1}^{q}d_{i}$. For any row-block
$A_{i_{1}}^{1}\otimes\ldots\otimes A_{i_{q}}^{(q)}$, computing
$S_{i1i2...i_{q}}(A_{i_{1}}^{1}\otimes\ldots\otimes A_{i_{q}}^{(q-1)})$ takes
$O(d(\sum_{k=1}^{q}nnz(A_{i_{k}}^{(k)}))+dqm\log(m))$ (see Sec 2). Hence for
$S{\cal A}$, it takes:

$\displaystyle\left(d\sum_{k=1}^{q}\operatorname{\mathtt{nnz}}(A_{k})\prod_{i\in[q]\backslash\\{k\\}}^{q}n_{i}/w_{i}\right)+\left(dqm\log(m)\prod_{i=1}^{q}n_{i}/w_{i}\right)$

where
$S\in{\mathbb{R}}^{(m\prod_{i=1}^{q}(n_{i}/w_{i}))\times\prod_{i=1}^{q}w_{i}}$
and $m\geq
100\prod_{i=1}^{q}d_{i}^{2}(2+3^{q})/\epsilon^{2}=O(\operatorname{poly}(d/\epsilon))$.
We need to compute an orthogonal factorization $S{\cal A}=QR_{\cal A}$ in
$O(qmd^{2})$ and then compute $U=R_{\cal A}^{-1}$ in $O(d^{3})$ time. Hence
the total running time of Algorithm Condition(${\cal A}$) is
$O(qmd^{2}+d^{3})$. Thus the total running time of computing $S{\cal A}$ and
Condition(A) is

$\displaystyle
O\left(\left(\sum_{k=1}^{q}\operatorname{\mathtt{nnz}}(A_{k})\prod_{i\in[q]\backslash\\{k\\}}^{q}n_{i}/w_{i}\right)+\left(\prod_{i=1}^{q}n_{i}/w_{i}\right)\operatorname{poly}(d/\epsilon)+qmd^{2}+d^{3}\right),$

We will compute $UG$ in $O(d^{2}\log n)$ time. We compute
$\widetilde{E}=E(A_{q_{1}+1}\otimes\ldots\otimes A_{q})^{T}$ in
$O(dn_{[q]\backslash[q_{1}]})$ time.

Then we can compute $R(A_{1}\otimes\cdots\otimes A_{q_{1}})\widetilde{E}_{j}$
in $O(n_{[q_{1}]}d_{[q_{1}]}\log n+d_{[q_{1}]}n_{[q]\backslash[q_{1}]}\log n)$
time.

Since computation of the median $\lambda_{i}$ takes $O(\log n)$ time,
computing all $\lambda_{i}$ and then $\lambda_{e}$ takes
$O(n_{[q]\backslash[q_{1}]}\log n)$ time.

As ${\cal A}UG$ has $O(\log n)$ columns, we need to compute $\lambda_{e}$ for
each ${\cal A}UG$ using the above procedure and hence it takes in total
$O(d(n_{[q_{1}]}+n_{[q]\backslash[q_{1}]})\log^{2}n)$ time.

Sampling a column of ${\cal A}UG$ using $\lambda_{e}$ takes $O(\log n)$ time,
sampling an entry in $M$ takes in total
$O(n_{[q_{1}]}+n_{[q]\backslash[q_{1}]})$ time.

Since we need $\sqrt{\prod_{k=1}^{q}w_{k}}\operatorname{poly}(r)$ samples to
select rows, the running time is
$d(n_{[q_{1}]}+n_{[q]\backslash[q_{1}]})\log^{2}n\cdot\sqrt{\prod_{k=1}^{q}w_{k}}\operatorname{poly}(r)$.

Now for simplicity, we set $q=2$, $n_{i}=n_{0}$ for $i\in[2]$. Note that it is
optimal to choose $w_{i}=w$ for $i\in[2]$. Substituting $q=2,n_{i}=n_{0}$ and
$w_{i}=w$, we that the total running time of Alg. 3:

$\displaystyle
O\left(dw^{-1}n_{0}(\operatorname{\mathtt{nnz}}(A_{1})+\operatorname{\mathtt{nnz}}(A_{2}))+w^{-2}n_{0}^{2}\operatorname{poly}(d/\epsilon)+wn_{0}\operatorname{poly}(d)\log(n)\right).$

For dense $A_{1}$ and $A_{2}$,
$\operatorname{\mathtt{nnz}}(A_{1})+\operatorname{\mathtt{nnz}}(A_{2})=O(n_{0})$
time, and so ignoring $\operatorname{poly}$ and $\log$ terms that do not
depend on $n_{0}$, the total running time can be simplified to:

$\displaystyle O(w^{-1}n_{0}^{2}+wn_{0}).$

Setting $w=\sqrt{n_{0}}$, we can minimize the above running time to
$O(n_{0}^{3/2})$, which is faster than the $n_{0}^{2}$ time for solving the
problem by forming $A_{1}\otimes A_{2}$. ∎