Recycling Randomness with Structure for Sublinear time Kernel Expansions

Introduction

Consider a $k$-dimensional feature map of the form,

In recent years, such random feature maps have been used to dramatically accelerate the training time and inference speed of kernel methods (Schölkopf & Smola, 2002) across a variety of statistical modeling problems (Rahimi & Recht, 2007; Xie et al., 2015) and applications (Huang et al., 2014; Vedaldi & Zisserman, 2012). Standard linear techniques applied to random nonlinear embeddings of data are equivalent to learning with approximate kernels. To quantify the benefits, consider solving a kernel ridge regression task given $l$ training examples. With traditional kernel methods, dense linear algebra operations on the Gram matrix associated with the exact kernel function imply that the training complexity grows as $O(l^{3}+l^{2}n)$ and the time to make a prediction on a test sample grows as $O(ln)$. By contrast, random feature approximations reduce training complexity to $O(lk^{2}+lkn)$ and test speed to $O(kn)$. This is a major win on big datasets where $l$ is very large, provided that a small value of $k$ can provide a good approximation to the kernel function.

In practice, though, the optimal value of $k$ is often large, albeit still much smaller than $l$. For example, in a speech recognition application (Huang et al., 2014) involving around two million training examples, about hundred thousand random features are required to achieve state of the art results. In such settings, the time to construct the random feature map is dominated by matrix multiplication against the dense Gaussian random matrix, which becomes the new computational bottleneck. To alleviate this bottleneck, (Le et al., 2013) introduce the “Fastfood” approach where Gaussian random matrices are replaced by Hadamard matrices combined with diagonal matrices with Gaussian distributed diagonal entries. It was shown in (Le et al., 2013) that for the specific case of the complex exponential nonlinearity, the Fastfood feature maps provide unbiased estimates for the Gaussian kernel function, at the expense of additional statistical variance, but with the computational benefit of reducing the feature map construction time from $O(kn)$ to $O(k\log~{}n)$ by using the Fast Walsh-Hadamard transform for matrix multiplication. The Fastfood construction for kernel approximations is akin to the use of structured matrices - in lieu of Gaussian random matrices - in Fast Johnson-Lindenstrauss transform (FJLT) (Alon & Chazelle, 2009) for dimensionality reduction, fast compressed sensing (Bajwa et al., 2007; Rauhut et al., 2012), and randomized numerical linear algebra techniques (Alon & Chazelle, 2011; Mahoney, 2011) Specific structured matrices were recently applied for approximating angular kernels (Choromanska et al., 2016). Some heuristic results for approximating kernels with circulant matrices were given in (Yu et al., 2015).

Our contributions in this paper are as follows:

We study a general family of structured random matrices that can be constructed by recycling a Gaussian random vector using a sequence of elementary generator matrices (introduced in Section 3). This family includes Circulant, Toeplitz and Hankel matrices. It also includes the Fastfood construction of (Le et al., 2013) as a special case. We show that fast sublinear time random feature maps obtained from these matrices provide unbiased estimates of the exact kernel, with variance comparable to the fully unstructured Gaussian case (Section 4). We introduce various structural coherence and graph-theoretic constants that control the quality of randomness we get from our model. Our approach generalizes across various choices of nonlinearities and kernel functions.

Of particular interest for us is the class of generalized structured matrices that have low-displacement rank (Pan, 2001; Sindhwani et al., 2015). Such matrices span an increasingly rich class of structures as the displacement rank is increased: from Circulant and Toeplitz matrices, to inverses and products of Toeplitz matrices, and more. The displacement rank provides a knob with which the degree of structure and randomness can be controlled to tradeoff computational and storage requirements against statistical variance.

We provide empirical support for our theoretical results (Section 5). In particular, we show that Circulant, Fastfood and low-displacement Toeplitz-like matrices provide high quality sublinear-time feature maps for approximating various kernels. With increasing displacement rank, the quality of the approximation approaches that of the fully Gaussian random matrix.

Background and Preliminaries

We start by giving a brisk background on random feature maps and structured matrices.

Random feature maps may be viewed as arising from Monte-Carlo approximations to integral representations of kernel functions. The original construction by Rahimi & Recht (2007) was motivated by a classical result that characterizes the class of shift-invariant positive definite functions.

While studying synergies between kernel methods and deep learning, (Cho & Saul, 2009) introduce $b^{th}$-order arc-cosine kernels via the following integral representation:

where $i(\cdot)$ is the step function, i.e. $i(x)=1$ if $x>0$ and otherwise; and the density $p$ is chosen to be standard Gaussian. These kernels evaluate inner products in the representation induced by an infinitely wide single hidden layer neural network with random Gaussian weights, and admit closed form expressions in terms of the angle $\theta=cos^{-1}(\frac{\mathbf{x}^{T}\mathbf{z}}{\|\mathbf{x}\|_{2}\|\mathbf{z}\|_{2}})$ between $\mathbf{x}$ and $\mathbf{z}$:

where $\|\cdot\|_{2}$ denotes $l_{2}$ norm.

Monte Carlo approximations to the integral representations above lead to the following,

where the feature map $\Psi(\mathbf{x})$ has the form given in Eqn. 1, with rows of $\mathbf{M}$, i.e. $\mathbf{w}_{j}$ vectors, drawn from the Gaussian density, and the nonlinearity $s$ set to the following: complex exponential, $s(x)=e^{i\frac{x}{\sigma}}$, for the Gaussian kernel with bandwidth $\sigma$; hard-thresholding, $s(x)=i(x)$, for the angular similarity kernel in Eqn. 2; and ReLU activation, $s(x)=\max(x,0)$, for the first order arc-cosine kernel in Eqn. 3.

2 Structured Matrices

A $m\times n$ matrix is called a structured matrix if it satisfies the following two properties: (1) it has much fewer degrees of freedom than $mn$ independent entries, and hence can be implicitly stored more efficiently than general matrices, and (2) the structure in the matrix can be exploited for fast linear algebra operations such as fast matrix-vector multiplication. Examples include the Discrete Fourier Transform (DFT), the Discrete Cosine Transform (DCT) and the Walsh-Hadamard Transform (WHT) matrices. Here, we give other examples particularly relevant to this paper. The matrices described below are square. Rectangular matrices can be obtained by appropriately selecting rows or columns.

Circulant Matrices: These matrices are intimately associated with circular convolutions and have been used for fast compressed sensing in (Rauhut et al., 2012). A $n\times n$ Circulant matrix is completely determined by its first column/row, i.e., $n$ parameters. Each column/row of a Circulant matrix is generated by cyclically down/right-shifting the previous column/row. A skew-Circulant matrix has identical structure to Circulant, except that the upper triangular part of the matrix is negated. This general structure looks like,

with $f=1$ for Circulant and $f=-1$ for skew-Circulant matrix. Both these matrices admit $O(n~{}\log~{}n)$ matrix-vector multiplication as they are diagonalized by the DFT matrix (Pan, 2001). We will use the notation $\texttt{circ}[\mathbf{g}]$ and $\texttt{scirc}[\mathbf{g}]$ for Circulant and skew-Circulant matrices respectively.

Toeplitz and Hankel Matrices: These matrices implement discrete linear convolution and arise naturally in dynamical systems and time series analysis. Toeplitz matrices are characterized by constant diagonals as follows,

Closely related Hankel matrices have constant anti-diagonals. Toeplitz-vector multiplication can be reduced to $O(n~{}\log~{}n)$ Circulant-vector multiplication. For detailed properties of Circulant and Toeplitz matrices, we point the reader to (Gray, 2006)

For a given displacement rank parameter $r$, the class of matrices for which the rank of $L[\mathbf{T}]$ is at most $r$ is called Toeplitz-like. Remarkably, this class of matrices admits a closed-form parameterization in terms of the low-rank factorization of $L[\mathbf{T}]$:

: If an $n\times n$ matrix $\mathbf{T}$ satisfies $rank(\mathbf{Z}_{1}\mathbf{T}-\mathbf{T}\mathbf{Z}_{-1})\leq r$, then it can be written as,

The family of matrices expressible by Eqn. 15 is very rich (Pan, 2001), i.e., it covers (i) all Circulant and Skew-circulant matrices for $r=1$, (ii) all Toeplitz matrices and their inverses for $r=2$, (iii) Products, inverses, linear combinations of distinct Toeplitz matrices with increasing $r$, and (iv) all $n\times n$ matrices for $r=n$. Since Toeplitz-like matrices under the parameterization of Eqn. 15 are a sum of products between Circulant and Skew-circulant matrices, they inherit fast FFT based matrix-vector multiplication with cost $O(nrlog~{}n)$, where $r$ is the displacement rank. Hence, $r$ provides a knob on the degree of structure imposed on the matrix with which storage requirements, computational constraints and statistical capacity can be explicitly controlled. Recently such matrices were used in the context of learning mobile-friendly neural networks in (Sindhwani et al., 2015). We note in passing that the displacement rank framework generalizes to other types of base structures (e.g. Vandermonde); see (Pan, 2001).

3 FastFood

In the context of fast kernel approximations, (Le et al., 2013) introduce the Fastfood technique where the matrix $\mathbf{M}$ in Eqn. 1 is parameterized by a product of diagonal and simple matrices as follows:

Here, $\mathbf{S},\mathbf{G},\mathbf{B}$ are diagonal random matrices, $\mathbf{P}$ is a permutation matrix and $\mathbf{H}$ is the Walsh-Hadamard matrix. The $k\times n$ matrix $\mathbf{M}$ is obtained by vertically stacking $k/n$ independent copies of the $n\times n$ matrix $\mathbf{F}$. Multiplication against such a matrix can be performed in time $O(k\log~{}n)$. The authors prove that (1) the Fastfood approximation is unbiased, (2) its variance is at most the variance of standard Gaussian random features with an additional $O(\frac{1}{k})$ term, and (3) for a given error probability $\delta$, the pointwise approximation error of a $n\times n$ block of Fastfood is at most $O(\sqrt{\log(n/\delta)})$ larger than that of standard Gaussian random features. However, note that the Fastfood analysis is limited to the Gaussian kernel and their variance bound uses properties of the complex exponential. The authors also conjecture that the Hadamard matrix $\mathbf{H}$ above, can be replaced by any matrix $\mathbf{T}$ such that $\mathbf{T}/\sqrt{n}$ is orthonormal, the maximum entry in $\mathbf{T}$ is small, and matrix-vector product against $\mathbf{T}$ can be computed in $O(n\log~{}n)$ time.

Structured Matrices from Gaussian Vectors

In this section, we present a general structured matrix model that allows a small Gaussian vector to be recycled in order to mimic the properties of a Gaussian random matrix suitable for generating random features. We first introduce some basic concepts in our construction. Note that we emphasize intuitions in our exposition - formal proofs are provided in our supplementary material.

Budget of Randomness: Let $t$ be some given parameter. Consider the column vector $\mathbf{g}=(g_{1},...,g_{t})^{T}$, where each entry is an independent Gaussian taken from $\mathcal{N}(0,1)$. This vector stands for the “budget of randomness” used in our structured matrix construction scheme.

Our goal is to recycle the Gaussian vector $\mathbf{g}$ to construct random matrices with desirable properties. This is accomplished using a sequence of matrices which we call the $\mathcal{P}$-model.

where $\mathbf{g}$ is a Gaussian random vector of length $t$.

In the constructions of interest to us, the sequence $\mathcal{P}$ is designed to separate structure from Gaussian randomness; though elements of $\mathcal{P}$ can be deterministic or itself random, Gaussianity is restricted to the vector $\mathbf{g}$. The ability of $\mathcal{P}$ to recycle a Gaussian vector effectively depends on certain structural constants that we now define.

For $\mathcal{P}=\{\mathbf{P}_{i}\}_{i=1}^{m}$, let $\mathbf{P}_{ij}$ denote the $j^{th}$ column of the $i^{th}$ matrix. The coherence of a $\mathcal{P}$-model is defined as,

Note that $\mu[\mathcal{P}]$ is a maximum over all pairs of rows $1\leq i\leq j\leq m$ of the rescaled sums of cross-correlations $\mathbf{P}^{T}_{i,n_{1}}\mathbf{P}_{j,n_{2}}$ for all pairs of different column indices $n_{1},n_{2}$. Lower values of $\mu[\mathcal{P}]$ will lead to better quality models. In practice, as we will see in subsequent analysis, it suffices if $\mu[\mathcal{P}]=O(poly(\log(n)))$ which is the case for instance for Toeplitz and Circulant matrices.

The coherence of the $\mathcal{P}$-model is an extremal statistic of pairwise correlations. We couple it with another set of objects describing global structural properties of the model, namely the coherence graphs.

Let $1\leq i,j\leq m$. We define by $\mathcal{G}_{i,j}$ an undirected graph with the set of vertices $V(\mathcal{G}_{i,j})=\{\{n_{1},n_{2}\}:1\leq n_{1}\neq n_{2}\leq n$ and $\mathbf{P}^{T}_{i,n_{1}}\mathbf{P}_{j,n_{2}}\neq 0\}$ and the set of edges $E(\mathcal{G}_{i,j})=\{\{\{n_{1},n_{2}\},\{n_{2},n_{3}\}\}:\{n_{1},n_{2}\},\{n_{2},n_{3}\}\in V(\mathcal{G}_{i,j})\}$. In other words, edges are between these vertices such that their corresponding $2$-element subsets intersect. The chromatic number $\chi(i,j)$ of a graph $\mathcal{G}_{i,j}$ is the smallest number of colors that can be used to color all vertices of $\mathcal{G}_{i,j}$ in such a way that no two adjacent vertices share the same color.

The chromatic number of a $\mathcal{P}$-model is defined as follows:

The chromatic number $\chi[\mathcal{P}]$ of a $\mathcal{P}$-model is given as:

where $\mathcal{G}_{i,j}$ are associated coherence graphs.

As it was the case for the coherence $\mu[\mathcal{P}]$, smaller values of the chromatic number $\chi[\mathcal{P}]$ lead to better theoretical results regarding the quality of the model. Intuitively speaking, coherence graphs encode in a compact combinatorial way correlations between different rows of the structured matrix produced by the $\mathcal{P}$-model. The chromatic number $\chi[\mathcal{P}]$ is a single combinatorial parameter measuring quantitatively these dependencies. It can be easily computed or at least upper-bounded (which is enough for us) for $\mathcal{P}$-models related to all structured matrices considered in this paper. The following is a well-known fact from graph theory:

The chromatic number $\chi(G)$ of an undirected graph $G$ with maximum degree $d_{max}$ satisfies: $\chi(G)\leq d_{max}+1$.

For all instantiations of $\mathcal{P}$-models considered in this paper leading to various structured matrices, the vertices of associated coherence graphs will turn out to have small degrees and hence, by Lemma 3.1, small chromatic numbers.

We will introduce one more structural parameter of the $\mathcal{P}$-model, depending on whether it is specified deterministically or randomly.

2 Examples of 𝒫𝒫\mathcal{P}-model structured matrices

Below we observe that various structured random matrices can be constructed according to the $\mathcal{P}$-model, i.e. by specifying a sequence of matrices $\mathbf{P}_{i}$ in Eqn. 17. We note that chromatic numbers and coherence values of these $\mathcal{P}$-models are low. In the next section, we show that this implies that we can get unbiased, low-variance kernel approximations from these matrices, for various choices of nonlinearities. Here we consider square structured matrices for which $m=n$, or rectangular matrices with $m<n$ obtained by selecting first $m$ rows of a structured matrix.

2.2 Toeplitz and Hankel matrices

2.3 Fastfood matrices

2.4 Toeplitz-like semi-Gaussian matrices

Random discretized vectors $\mathbf{h}^{i}$: Each dimension of each $\mathbf{h}^{i}$ is chosen independently at random from the binary set $\{-\frac{1}{\sqrt{nr}},\frac{1}{\sqrt{nr}}\}$.

Sparse setting: Each $\mathbf{h}^{i}$ is sparse (but nonzero), i.e. has only few nonzero entries. Furthermore, the sign of each $\mathbf{h}^{i}_{j}$ is chosen independently at random and the following holds: $\|\mathbf{h}^{1}\|^{2}+...+\|\mathbf{h}^{r}\|^{2}=1$. This setting is characterized by a parameter $\kappa$ defining the size of the set of dimensions that are nonzero for at least one $\mathbf{h}^{i}$.

We refer to such matrices as Toeplitz-like semi-Gaussian matrices. We now sketch how they can be obtained from the $\mathcal{P}$-model. We take $t=nr$ and $\mathbf{g}=(g^{1}_{1},...,g^{1}_{n},...,g^{r}_{1},...,g^{r}_{n})^{T}$. The matrix $\mathbf{P}_{1}$ is constructed by vertically stacking $r$ matrices $\mathbf{S}_{j}$ for $j=1,...,r$, where each $\mathbf{S}_{j}$ is constructed as follows. The first column of $\mathbf{S}_{j}$ is $\mathbf{h}^{j}$ and the subsequent columns are obtained from previous by skew-Circulant downward shifts. Matrix $\mathbf{P}_{i}$ for $i>1$ is obtained from $\mathbf{P}_{i-1}$ by upward Circulant shifts, independently for each column at each block $\mathbf{S}_{j}$.

Matrices constructed according to this procedure satisfy conditions regarding certain structural parameters of the $\mathcal{P}$-model (see: Theorem 4.4). In particular, in the sparse semi-Gaussian setting the corresponding coherence graphs have vertices of degrees bounded by a constant; thus, by Lemma 3.1 the $\mathcal{P}$-models associated with them have low chromatic numbers.

3 Construction of Random Feature Maps

Given $S[\mathcal{P}]$, the $m\times n$ structured random matrix defined by a $\mathcal{P}$-model, in lieu of using the $k\times n$ Gaussian random matrix $\mathbf{M}$ in Eqn. 1, the feature map for a data vector $\mathbf{x}$ is constructed as follows.

Return $\Psi(\mathbf{x})=\frac{1}{\sqrt{k}}s(\mathbf{\bar{x}})$

Note that the displacement rank $r$ for low displacement rank matrices and the number of rows $m$ of a single structured block can be used to control the “budget of randomness”; $m=1$ reduces to a completely unstructured matrix.

Theoretical results

In this section we provide concentration results regarding $\mathcal{P}$-model for Gaussian and arc-cosine kernels, showing in particular that the variance of the computed structured approximation of the kernel is close to the unstructured one. We also present results targeting specifically low displacement rank structured matrices, and show how the displacement rank knob can be used to increase the budget of randomness and reduce the variance.

The orthogonality condition $\mathbf{P}^{T}_{i,j}\mathbf{P}_{i,k}=0$ is trivially satisfied by Hankel, circulant or Toeplitz structured matrices produced by the $\mathcal{P}$-model as well as Toeplitz-like semi-Gaussian matrices, where each $\mathbf{h}^{i}$ has one nonzero entry. It is also satisfied in expectation (which in practice suffices) for all presented Toeplitz-like semi-Gaussian matrices.

For a $\mathcal{P}$-model, where matrices $\mathbf{P}_{i}$ were chosen randomly we denote as $\eta[\mathcal{P}]$ the maximum possible value that a random variable $(\mathbf{P}_{i,n_{1}}^{T}\mathbf{P}_{j,n_{1}})^{2}$ can take for $1\leq i<j\leq m,1\leq n_{1}\leq n$. Without loss of generality we will assume that data vectors are drawn from the ball $\mathcal{B}(0,1)$ centered at of unit $l_{2}$ norm. Below we state results regarding $d^{th}$ moments of the obtained kernel’s approximation via the $\mathcal{P}$-model that lead to the concentration results.

and expectations are taken in respect to random choice for a Gaussian vector $\mathbf{g}$. If $\mathbf{P}_{i}s$ are chosen from the probabilistic model then the above holds with probability at least $1-p_{wrong}$ in respect to random choices of $\mathbf{P}_{i}s$, where

Let us comment on the result above. The upper bound is built from two main components: $p_{gen}$ and $p_{struct}$. The first one depends on the general parameters of the setting: dimensionality of the data $n$ and order of the computed moment $d$. The second one is crucial to understand how the structure of the matrix influences the quality of the model. We can immediately see that low chromatic numbers $\chi(i,j)$ (see: Section 3.1) improve quality since they decrease computed upper bound. Furthermore, low values of the coherence $\mu[\mathcal{P}]$ and chromatic number $\chi[\mathcal{P}]$ also lead to stronger concentration results. Both observations were noticed by us before, but now we see how they are implied by general theoretical results. Finally, for all considered settings, where matrices $\mathbf{P}_{i}$ are constructed randomly parameter $\eta[\mathcal{P}]$ is of order $O(1)$ thus $p_{wrong}$ in negligibly small.

In particular, if both the chromatic number $\chi[\mathcal{P}]$ and the coherence $\mu[\mathcal{P}]$ are of the order $O(poly(\log(n)))$ then $p_{struct}$ if inversely proportional to the superpolynomial function of $n$ thus is negligible in practice. That, as we will see soon, will be the case for proposed Toeplitz-like semi-Gaussian matrices with sparse vectors $h^{i}$.

Theorem 4.1 implies also that variances of the kernel approximation for the structured $\mathcal{P}$-model case and unstructured setting are very similar (we borrow denotation from Theorem 4.1).

Consider the setting as in Theorem 4.1. If matrices $\mathbf{P}_{i}$ are chosen deterministically then for any $T,\epsilon>0$ the following is true for $n$ large enough:

where $Var$ stands for the variance and $\Delta=p_{gen}(T)+p_{struct}+\epsilon$. If $\mathbf{P}_{i}s$ are chosen from the probabilistic model then the above holds with probability at least $1-p_{wrong}$, where $p_{wrong}$ is as in Theorem 4.1.

Note that in practice it means that the variance in the structured and unstructured setting is similar. In particular, choosing $\epsilon=O(\frac{1}{m^{2}})$, $T>7\log(m)$, one can deduce that the variance in the structured setting is of the order $O(\frac{1}{m})$ for $n$ large enough (the well known fact is that the unstructured variance is of the order $O(\frac{1}{m})$). Note also that as expected, for $m=1$ the structured setting becomes an unstructured one, since each structured block consists of just one row and different blocks are constructed independently.

Toeplitz-like semi-Gaussian Low-displacement rank matrices: Note that the structure of a matrix affects only the $p_{struct}$ factor in the statements above. Thus, we will focus on the structured parameters of the $\mathcal{P}$-model. We will show that Toeplitz-like semi-Gaussian matrices can be set up so that the above parameters are of required order.

The richness of the low displacement rank mechanism comes from the fact that the budget of randomness can be controlled by the rank parameter $r$ and increasing $r$ leads to better quality approximations. In particular, we have:

Note that increasing rank $r$ leads to sharper upper bounds on the coherence $\mu[\mathcal{P}]$ (in practice $r$ polynomial in $log(n)$ suffices) and thus, from what we have said so far, to better concentration results for the entire structured scheme. Analogous variance bounds can also be derived for Toeplitz-like semi-Gaussian matrices where the $\mathbf{h}^{i}$ vectors are chosen to be dense. But due to lack of space, these results are included in our supplementary material.

Empirical Support

In this section, we compare feature maps obtained with fully Gaussian, Fastfood, Circulant, and Toeplitz-like matrices with increasing displacement rank. Our goal is to lend support to the theoretical contributions of this paper by showing that high-quality feature maps can be constructed from a broad class of structured matrices as instantiations of the proposed $\mathcal{P}$-model.

Results on publicly available real-world classification datasets, averaged over $100$ runs, are reported in Table 1 for complex exponential nonlinearity (Gaussian kernel). Results with ReLU (arc-cosine) are similar but not shown for lack of space. As observed in previous papers, better Gram matrix approximation is not often correlated with higher classification accuracy. Nonetheless, it is clear that the design of space of valid feature map constructions based on structured matrices is much larger than what has so far been explored in the literature: Circulant and Toeplitz-like matrices are very competitive with Fastfood, and sometimes give better results particularly with increasing displacement rank. The effectiveness of such feature maps for nonlinearities other than the complex exponential also validates our theoretical contributions. Among the unstructured baselines, we also include Quasi-Monte Carlo (QMC) feature maps of (Yang et al., 2014) using Halton low-discrepancy sequences. The use of structured matrices to accelerate QMC techniques building on (Dick et al., 2015) is of interest for future work.

Speedups: Figure 3 shows the speedup obtained in featuremap construction time using structured matrices relative to using unstructured Gaussian random matrices (on a 6-core 32-GB Intel(R) Xeon(R) machine running Matlab R2014a). The benefits of sub-quadratic matrix-vector multiplication with FFT-variations tend to show up beyond $1024$ dimensions. Circulant-based feature maps are the fastest to compute. Fastfood (with DCT instead of Hadamard matrices) is about as fast as Toeplitz-like matrices with displacement rank 1 or 2. Higher displacement rank matrices show speedups at higher dimensions as expected. Fastfood with inbuilt fwht routine in Matlab performed poorly in our experiments.

Conclusions

We have theoretically justified and empirically validated the use of a broad family of structured matrices for accelerating the construction of random embeddings for approximating various kernel functions. In particular, the class of Toeplitz-like semi-Gaussian matrices allows our construction to span highly compact to fully random matrices.

References

Appendix

We now prove all theoretical results of the paper. We need to introduce some technical denotation.

where $g^{1},...,g^{m}$ is the set of $m$ gaussian vectors forming gaussian matrix $G$, each obtained by sampling independently $n$ values from the distribution $\mathcal{N}(0,1)$ and $\phi_{i}$s differ by the choice of nonlinear mapping $f_{i}\in\mathcal{F}$, can be accurately approximated by its structured version $T^{A,d}_{v^{1},v^{2}}((\mathcal{R},\alpha_{1},...,\alpha_{r})$ which is of the form:

where $a^{1},...,a^{m}$ are rows of the structured matrix $A=G_{struct}^{i}$. The importance of $T^{G,d}_{v^{1},v^{2}}(\mathcal{R},\alpha_{1},...,\alpha_{r})$ and $T^{A,d}_{v^{1},v^{2}}(\mathcal{R},\alpha_{1},...,\alpha_{r})$ lies in the fact that $d^{th}$ moments of the random variables approximating considered kernels in the unstructured and structured mechanism can be expressed as weighted sums of the expressions of the form $T^{G,d}_{v^{1},v^{2}}(\alpha_{1},...,\alpha_{r})$ and $T^{A,d}_{v^{1},v^{2}}(\alpha_{1},...,\alpha_{r})$ respectively if $\Psi(x_{1},...,x_{r})=x_{1}\cdot...\cdot x_{r}$. Thus if $T^{A,d}_{v^{1},v^{2}}(\alpha_{1},...,\alpha_{r})$ closely approximates $T^{G,d}_{v^{1},v^{2}}(\alpha_{1},...,\alpha_{r})$ then the corresponding moments are similar. That, as we will see soon, implies several theoretical guarantees for the structured method. In particular, this means that the variances are similar. Since in the unstructured setting the variance is of the order $O(\frac{1}{m})$, that will be also the case for the structured setting. This in turn will imply concentration results providing theoretical explanation for the observations from the experimental section that show the quality of the proposed structured setting.

Note that the value of the function $\phi_{i}(v^{1},v^{2},g^{i})^{\alpha_{i}}$ depends only on the projection $g^{i}_{proj}$ of $g^{i}$ on the $2$-dimensional space spanned by $v^{1}$ and $v^{2}$. Thus for a given pair $v^{1},v^{2}$ function $\phi$ is in fact a function $B^{v^{1},v^{2}}_{i}$ of this projection.

where the supremum is taken over all indices $i=1,...,m$, all pairs of linearly independent vectors from the domain, all coordinate systems in $span(v^{1},v^{2})$ and vectors $\zeta$ of $L_{1}$-norm at most $\epsilon$ in some of these coordinate systems.

We will use the following notation: $\sigma_{i,j}(n_{1},n_{2})=\mathbf{P}_{i,n_{1}}^{T}\mathbf{P}_{j,n_{2}}$. To compress the statements of our theoretical results, we will use also the following notation:

Note first that the preprocessing step preserves kernels’ values since transformation $HD_{0}$ is an isometry and considered kernels are spherically-invariant. We start with Lemma 4.1.

where $G$ is the unstructured gaussian matrix. Let $g^{i,j}_{struct}$ be the $j^{th}$ row of $G^{i}_{struct}$ and let $g^{j}$ be the $j^{th}$ row of $G$. Note that we have:

The latter follows from the fact that $g^{i,j}_{struct}$ has the same distribution as $g$. To see this note that $g^{i,j}_{struct}=g\cdot P_{i}$. Thus dimensions of $g^{i,j}_{struct}$ are projections of $g$ onto columns of $P_{i}$. Each projection is trivially gaussian from $\mathcal{N}(0,1)$ (that is implied by the fact that each column is normalized). The independence of different dimensions of $g^{i,j}_{struct}$ comes from the observation that different columns are orthogonal. Thus we can use a simple property of gaussian vectors stating that the projections of a gaussian vector on mutually orthogonal directions are independent. The equation 25 implies equation 26 by the linearity of expectation and that completes the proof. ∎

Now we prove Theorem 4.1. This one is easily implied by a more general result that we state below. We will assume that function $\Psi$ from equations: 22, 23 is $M$-bounded for some given $M>0$. We will assume that expected values defining $T^{A,d}$ are not with respect to the random choices determining $P_{i}s$.

If $P_{i}s$ are chosen from the probabilistic model then the above holds with probability at least $1-p_{wrong}$, where $p_{wrong}=2\sum_{i\leq i_{1}<i_{2}\leq m}e^{-\frac{n^{2}}{8\log^{6}(n)\sum_{j=1}^{n}(\sigma_{i_{1},i_{2}}(j,j))^{2}}}$.

We will show that $s^{i,j}$s, even though not necessarily pairwise orthogonal, are close to be pairwise orthogonal with high probability. Let us assume now that vectors $s^{i,j}$ can be chosen in such a way that each $s^{i,j}$ satisfies: $s^{i,j}=w^{i,j}+\rho(i,j)$, where vectors $w^{i,j}$ are mutually orthogonal, we have $\|s^{i,j}\|_{2}=\|w^{i,j}\|_{2}$ and furthermore $\|\rho(i,j)\|_{2}\leq\rho$ for some given $\rho>0$. We call this property the $\rho$-orthogonality property. We will later show that the $\rho$-orthogonality property depends on the random diagonal matrix $D_{1}$.

Assume now that the $\rho$-orthogonality property is satisfied. Denote by $g^{\mathcal{H}}$ the projection of the “budget-of-randomness” vector $g$ onto $2r$-dimensional linear space $\mathcal{H}$ spanned by vectors from $\{s^{i,j}\}$. Note that then the coordinates of $a^{i}_{proj}s$ in $\mathcal{B}$ can be rewritten as $g\cdot w^{i,j}+\epsilon(i,j)$, where $|\epsilon(i,j)|\leq\epsilon$ and $\epsilon=\|g^{\mathcal{H}}\|_{2}\rho$. Thus each $\psi_{i}$ in the formula from equation 23 can be then expressed as $B_{i}^{v^{1},v^{2}}(g^{i}_{proj}+\epsilon(i))$, where $g^{i}_{proj}s$ stand for the projections onto $2$-dimensional linear space spanned by $v^{1}$ and $v^{2}$ of independent copies of gaussian vectors $g^{i}$. Each $g^{i}$ is of the same distribution as the corresponding structured vector $a^{i}$ and $\epsilon(i)s$ are vectors with the $L_{1}$-norm satisfying $\|\epsilon(i)\|\leq\epsilon$. The independence comes from the fact that variables of the form $g\cdot w^{i,j}$ are independent. That, as in the proof of Lemma 4.1 is implied by the well known fact that dot products of a given gaussian vector with orthogonal vectors are independent. Note that if not the term $\epsilon(i)$ then the formula for $T^{A,d}$ would collapse to its unstructured counterpart $T^{G,d}$. We will argue that both expressions are still close to each other if $\epsilon(i)$ have small $L_{1}$-norm.

Let us fix $\lambda>0$. Our goal is to count these indices $i$ that satisfy the following: $|\psi_{i}(v^{1},v^{2},g^{i})^{\alpha^{i}}-\psi_{i}(v^{1},v^{2},g^{i})^{\alpha^{i}}|>\lambda$, where $g^{i}s$ corresponds to the aforementioned independent counterparts of $a^{i}s$. We call them bad indices. Based on what we have said so far, we can conclude that the latter inequality can be expressed as $|B_{i}^{v^{1},v^{2}}(g^{i}_{proj}+\epsilon(i))-B_{i}^{v^{1},v^{2}}(g^{i}_{proj})|>\lambda$. Let us first find the upper bound on the probability of the event that the number of bad indices is $j$ for some fixed $1\leq j\leq d$. Note that since $g^{i}s$ are independent, we can use Bernoulli scheme to find that upped bound. Using the definition of $p_{\lambda,\epsilon}$ we obtain an upper bound of the form $p_{upper}\leq{d\choose j}(p_{\lambda,\epsilon})^{j}$. If the number of bad indices is $j$ then by the definition of $M$ and $\Delta^{\Psi}_{\lambda}$ we see that $T^{A,d}$ differs from $T^{G,d}$ by at most $iM+(d-i)\Delta^{\Psi}_{\lambda}$. Summing up over all indices $j$ we get the second term of the upper bound on $q^{d}_{v^{1},v^{2}}$ from the statement of the theorem.

However the $\rho$-orthogonality does not have to hold. Note that (by the definition of $\Psi_{max}$) to finish the proof of the theorem it suffices to show that the probability of $\rho$-orthogonality not to hold is at most $p_{gen}+p_{struct}$.

The $\rho$-orthogonality property holds with probability at least $1-(p_{gen}+p_{struct})$.

Let $x=(x_{1},...,x_{n})$ be a vector with $\|x\|_{2}=1$. We say that $x$ is $\theta$-balanced if $|x_{i}|\leq\frac{\theta}{\sqrt{n}}$ for $i=1,...,n$.

We will use the following concentration inequality, calles Azuma’s inequality

Let $X_{1},...,X_{n}$ be a martingale and assume that $-\alpha_{i}\leq X_{i}\leq\beta_{i}$ for some positive constants $\alpha_{1},...,\alpha_{n},\beta_{1},...,\beta_{n}$. Denote $X=\sum_{i=1}^{n}X_{i}$. Then the following is true:

for $l=1,...,t$, where $s^{i,j}_{l}$ stands for the $l^{th}$ dimension of $s^{i,j}$, $p^{i}_{l,k}$ is the entry in the $l^{th}$ row and $k^{th}$ column of $P_{i}$ and $d_{r}s$ are the values on the diagonal of the matrix $D_{0}$. As we noted earlier, we want to show that $s^{i,j}s$ are close to be mutually orthogonal. To do it, we will compute dot products $s^{i_{1},j_{1}}\cdot s^{i_{2},j_{2}}$. We will first do it for $i_{1}=i_{2}$. We have:

Now we take advantage of the normalization property of the matrices $P_{i}$ and the fact that $x^{1}$ is orthogonal to $x^{2}$ and conclude that the first term on the RHS of the equation above is equal to . Thus we have:

Note that if for any fixed $P_{i}$ any two different columns of $P_{i}$ are orthogonal then $\sigma_{i_{1},i_{1}}(n_{1},n_{2})=0$ and thus $s^{i_{1},j_{1}}\cdot s^{i_{1},j_{2}}=0$. This is the case for many structured matrices constructed according to the $\mathcal{P}$-model, for instance circulant, Toeplitz or Hankel matrices.

Let us consider now $s^{i_{1},j_{1}}\cdot s^{i_{2},j_{2}}$ for $i_{1}\neq i_{2}$. By the previous analysis, we obtain:

This time in general we cannot get rid of the first term in the RHS expression. This can be done if columns of the same indices in different $P_{i}s$ are orthogonal. This is in fact again the case for circulant, Toeplitz or Hankel matrices.

For $\{n_{1},n_{2}\}$ such that $n_{1}\neq n_{2}$ and $\sigma_{i_{1},i_{1}}(n_{1},n_{2})\neq 0$ let us now consider random variables $Y_{n_{1},n_{2}}$ that are defined as follows

From the definition of the chromatic number $\chi(i_{1},i_{1})$ we can deduce that the set of all this random variables can be partitioned into at most $\chi(i_{1},i_{1})$ subsets such that random variables in each subset are independent. Let us denote these subsets as: $\mathcal{L}_{1},...,\mathcal{L}_{r}$, where $r\leq\chi(i_{1},i_{1})$. Note that an event $\{|\sum_{1\leq n_{1}<n_{2}\leq n}d_{n_{1}}d_{n_{2}}x^{j_{1}}_{n_{1}}x^{j_{2}}_{n_{2}}2\sigma_{i_{1},i_{1}}(n_{1},n_{2})|>\kappa\}$ is contained in the sum of the events: $\mathcal{E}=\mathcal{E}_{1}\cup...\cup\mathcal{E}_{r}$, where each $\mathcal{E}_{j}$ is defined as follows:

Now we can use Azuma’s inequality to find an upper bound on $\mathcal{P}[\mathcal{E}_{i}]$ and we obtain:

We can conclude, using the union bound again, that for a $\log(n)$-balanced basis $\mathcal{B}$ the probability that there exist $i_{1},j_{1},j_{2}$ such that: $|s^{i_{1},j_{1}}\cdot s^{i_{1},j_{2}}|>\kappa$ is at most

Assume first that $P_{i}s$ are chosen deterministically. Note that by $\log(n)$-balanceness, we have:

Thus, by the triangle inequality, we have:

Using the same analysis as before, we then obtain the following bound on $p_{bad}(\kappa,\theta)$:

We can conclude that in the setting where $P_{i}s$ are chosen deterministically, under our assumptions on $\lambda(i_{1},i_{2})$, for $\kappa>0$ that does not depend on $n$ and $n$ large enough the following is true. The probability that there exist two different vector $s^{i_{1},j_{1}}$, $s^{i_{2},j_{2}}$ such that $|s^{i_{1},j_{1}}\cdot s^{i_{2},j_{2}}|>\kappa$ satisfies:

Now let us assume that $P_{i}s$ are chosen probabilistically. In that setting we also assume that columns of different indices are chosen independently (this is the case for instance for the FastFood Transform). Let us now denote:

If we now take $a=\frac{1}{\log(n)}$ and under $\log(n)$-balanceness assumption, we obtain:

Thus, using our bound on $Y$ for a fixed $\kappa$ and $n$ large enough we can repeat previous analysis and conclude that in the probabilistic setting of $P_{i}s$ the following is true:

Thus we can conclude that in both the deterministic and probabilistic setting for $P_{i}s$ we get:

Now we will show that the squared lengths of vectors $s^{i,j}$ are well concentrated around their means and that these means are equal to $1$. Let us remind that we have:

where the last inequality comes from the fact that each column of each $P_{i}$ has $l_{2}$-norm equal to $1$.

We can again apply Azuma’s inequality and the union bound as we did before and obtain:

where $p_{s}=4\sum_{i=1}^{m}\chi(i,i)e^{-\frac{1}{2\xi^{2}(i,i)\log^{2}(n)}\frac{n^{2}}{\log^{4}(n)}}$.

We will assume now that all $s^{i,j}$ satisfy: $|\|s^{i,j}\|^{2}_{2}-1|\leq\frac{1}{\log(n)}$, in particular:

Let us assume right now that the above inequality holds. Let $\{w^{i,j}\}$ be a set of vectors obtained from $\{s^{i,j}\}$ by the Gram-Schmidt process. Without loss of generality we can assume that $\|w^{i,j}\|_{2}=\|s^{i,j}\|_{2}$. Note that the size of the set $\{s^{i,j}\}$ is in fact not $2m$, but $2r$ and in all practical application $r\ll m$. Assume now that $|s^{i_{1},j_{1}}\cdot s^{i_{2},j_{2}}|\leq\kappa$ for any two different vectors $s^{i_{1},_{j_{1}}},s^{i_{2},j_{2}}$ and some fixed $\kappa>0$. Now, one can easily note that directly from the description of the Gram-Schmidt process that it leads to the set of vectors $\{w^{i,j}\}$ such that $\|s^{i,j}-w^{i,j}\|_{2}\leq\kappa\Gamma(2r)$, where $\Gamma$ is some constant that depends just on the size of the set $\{s^{i,j}\}$. Thus if we want $\rho$-orthogonality with $\rho=\frac{\epsilon}{\|g^{\mathcal{H}}\|_{2}}$, where $g^{\mathcal{H}}$ stands for the random projection of a vector $g$ onto $2r$-dimensional linear space spanned by vectors from $\{s^{i,j}\}$, then we want to have:

Thus we can conclude that the probability that $g^{\mathcal{H}}$ has $l_{2}$ norm larger than $\sqrt{2r}\cdot\sqrt{T}$ is at most $p_{gauss}(T)\leq\frac{4r}{\sqrt{2\pi T}}$. In such a case we need to take $\kappa$ of the form:

We are ready to finish the proof of Lemma 7.1. Take $\kappa=\frac{\epsilon}{\Gamma(2r)\sqrt{2r}\sqrt{T}}$. Let us first take the setting where $P_{i}s$ are chosen deterministically. Take an event $\mathcal{E}_{bad}$ which is the sum of the events which probabilisites are upper-bounded by $p_{gauss}(T)$, $1-p_{balanced}$, $p_{bad}(\kappa)$ and $p_{s}$. By the union bound, the probability of that event is at most $p_{gauss}+(1-p_{balanced})+p_{bad}(\kappa)+p_{s}$ which is upper-bounded by $p_{gen}+p_{struct}$ for $n$ large enough. Note that if $\mathcal{E}_{bad}$ does not hold then $\rho$-orthogonality is satisfied. Now let us take the probabilistic setting for choosing $P_{i}s$. We proceed similarly. The only difference is that right now we need to assume that the event upper-bounded by $p_{wrong}$ does not hold (this one depends only on the random choices for setting up $P_{i}s$). Thus again we get the statement of the lemma. That completes the proof of Lemma 7.1. ∎

As mentioned above, the proof of Lemma 7.1 completes the proof of the theorem. ∎

The last inequality in Eqn.62 is implied by the fact that different blocks of the structured matrix are computed independently and thus covariance related to rows from different blocks is .

where $X^{\mathcal{P}}_{i,j}$ stands for the version of $X_{i,j}$ if $\mathbf{A}$ was costructed via the $\mathcal{P}$-model and $X^{\mathbf{G}}_{i,j}$ stands for the fully unstructured one.

where the latter inequality is implied by the fact that different blocks are constructed independently.

where $\beta$ is an upper bound as in Theorem 7.1 for $d=2$. Now we can proceed in the same way as in the proof of Theorem 4.1 and the proof is completed. ∎

The fact that $\mu[\mathcal{P}]\leq\kappa$ comes directly from the definition of the coherence number and the sparse setting of semi-gaussian matrices. To see that, note that any given column $col$ of any matrix $\mathbf{P}_{i}$ in the related $\mathcal{P}$-model has a nonzero dot-product with at most $\kappa^{2}$ other columns of any matrix $\mathbf{P}_{j}$. This in turn is implied by the fact that different columns are obtained by applying skew-circulant shifts blockwise, thus the number of columns from $\mathbf{P}_{j}$ that have nonzero dot product with $col$ is at most the product of the number of nonzero dimensions of $col$ and $\mathbf{P}_{j}$. This is clearly upper bounded by $\kappa^{2}$. This leads to the upper bound on the coherence $\mu[\mathcal{P}]$.

The bound regarding the chromatic number is implied by the observation that each coherence graph in the corresponding $\mathcal{P}$-model has degree at most $\kappa^{2}$. That follows directly from the observation we used to prove the upper bound on $\mu[\mathcal{P}]$. But now we can use Lemma 3.1 and that completes the proof of Theorem 4.3. ∎

Below we present the proof of Theorem 4.4.

Fix two columns $\mathbf{P}_{i,n_{1}}$ and $\mathbf{P}_{j,n_{2}}$ and consider the expression $\mathbf{P}_{i,n_{1}}^{T}\mathbf{P}_{j,n_{2}}$. We have already mentioned in the previous proof the right approach to finding strong upper bound on $|\mathbf{P}_{i,n_{1}}^{T}\mathbf{P}_{j,n_{2}}|$. We first note that $\mathbf{P}_{i,n_{1}}^{T}\mathbf{P}_{j,n_{2}}$ can be written as a sum $w_{1}+...+w_{nr}$, where $w_{i}s$ are not necessarily independent but can be partitioned into at most three sets such that wariables in each of these sets are independent. This is true since $G^{i}_{struct}$ is produced by skew-circulant shifts and the corresponding coherence graphs has verrtices of degree at most $2$. Note also that each $w_{k}$ satisfies: $|w_{k}|\leq\frac{1}{\alpha r}$. In each of the sum we get rid of these $w_{i}s$ that are equal to . Then, by applying Azuma’s inequality independently on each of these subsets and taking union bound over these subsets, we conclude that for any $a>0$:

Now we can take the union bound over all pairs of columns and notice that for every columcn $col$ in $\mathbf{P}_{i}$ and any $\mathbf{P}_{j}$ there exists at most $\kappa$ columns in $\mathbf{P}_{j}$ that have nonzero dot product with $col$. We can then take $a=\frac{\tau}{\kappa}$ and the proof is completed. ∎

Let us now switch to dense semi-gaussian matrices. The following is true.

Consider the setting as in Theorem 4.1. Assume that entries of any fixed column of $P_{i}$ are chosen independently at random. Assume also that for any $1\leq i\leq j\leq m$ and any fixed column $col$ of $P_{i}$ each column of $P_{j}$ is a downward shift of $col$ by $b$ entries (possibly with signs of dimensions swapped) and that $b=0$ for $O(1)$ columns in $P_{j}$. Then for and $T>0$ and $n$ large enough the following holds:

where $\Delta=p_{gen}(T)+p_{struct}(T)+d\epsilon+e^{-n^{\frac{1}{3}}}$ and

The proof of this result follows along the lines of the proof of Theorem 4.1 and Theorem 4.2. Take the formulas for $s^{i_{1},j_{1}}\cdot s^{i_{2},j_{2}}$ derived in the proof of Theorem 7.1. Note that we want to have: $|s^{i_{1},j_{1}}\cdot s^{i_{2},j_{2}}|\leq\frac{\epsilon}{\Gamma(2d)\|g^{\mathcal{H}}\|_{2}}$, where $\Gamma$ is a constant that depends only on the degree $d$. Each $s^{i_{1},j_{1}}\cdot s^{i_{2},j_{2}}$ is a sum of random variables that can be decoupled into $O(1)$ subsums such that variables in each subsum are independent (here we use exactly the same trick as in the proof of Theorem 4.1). In each subsum we apply Azuma’s inequality. Straightforward computations lead to the conclusion that if one sets up $\epsilon$ as in the statement of Theorem 7.2 then the probability that there exist different $s^{i_{1},j_{1}}$, $s^{i_{2},j_{2}}$ such that $|s^{i_{1},j_{1}}\cdot s^{i_{2},j_{2}}|>\frac{\epsilon}{\Gamma(2d)\|g^{\mathcal{H}}\|_{2}}$ is of the order $e^{-n^{\frac{1}{3}}}$ for $n$ large enough. That is the extra term in the formula for $\Delta$ that was not present in the staement of Theorem 4.1. The variance results follows immediately by exactly the same analysis as in the proof of Theorem 4.2. ∎

Note that introduced dense semi-gaussian matrices trivially satisfy conditions of Theorem 7.2 (look for the description of matrices $\mathbf{P}_{i}$ from Subsection: 3.2.4). The role of rank is similar as in the sparse setting, i.e. larger values of $r$ lead to sharper concentration results. Theorem 7.2 can be applied to classes of matrices for which $|\sum_{1\leq n_{1}<n_{2}\leq n}\mathbf{P}_{i,n_{1}}^{T}\mathbf{P}_{j,n_{2}}|$ is small and random dense semi-gaussian matrices satisfy this condition with high probability.