Tensor sketch

The tensor sketch of Pham and Pagh[4] computes ${\displaystyle C^{(1)}x\ast C^{(2)}y}$, where ${\displaystyle C^{(1)}}$ and ${\displaystyle C^{(2)}}$ are independent count sketch matrices and ${\displaystyle \ast }$ is vector convolution. They show that, amazingly, this equals ${\displaystyle C(x\otimes y)}$ – a count sketch of the tensor product!

It turns out that this relation can be seen in terms of the face-splitting product as

${\displaystyle C^{(1)}x\ast C^{(2)}y={\mathcal {F}}^{-1}({\mathcal {F}}C^{(1)}x\circ {\mathcal {F}}C^{(2)}y)}$, where ${\displaystyle {\mathcal {F}}}$ is the Fourier transform matrix.

Since ${\displaystyle {\mathcal {F}}}$ is an orthonormal matrix, ${\displaystyle {\mathcal {F}}^{-1}}$ doesn't impact the norm of ${\displaystyle Cx}$ and may be ignored. What's left is that ${\displaystyle C\sim {\mathcal {C}}^{(1)}\bullet {\mathcal {C}}^{(2)}}$.

On the other hand,

${\displaystyle {\mathcal {F}}(C^{(1)}x\ast C^{(2)}y)={\mathcal {F}}C^{(1)}x\circ {\mathcal {F}}C^{(2)}y=({\mathcal {F}}C^{(1)}\bullet {\mathcal {F}}C^{(2)})(x\otimes y)}$.

The problem with the original tensor sketch algorithm was that it used count sketch matrices, which aren't always very good dimensionality reductions.

In 2020[15] it was shown that any matrices with random enough independent rows suffice to create a tensor sketch. This allows using matrices with stronger guarantees, such as real Gaussian Johnson Lindenstrauss matrices.

In particular, we get the following theorem

Consider a matrix ${\displaystyle T}$ with i.i.d. rows ${\displaystyle T_{1},\dots ,T_{m}\in \mathbb {R} ^{d}}$, such that ${\displaystyle E[(T_{1}x)^{2}]=\|x\|_{2}^{2}}$ and ${\displaystyle E[(T_{1}x)^{p}]^{1/p}\leq {\sqrt {ap}}\|x\|_{2}}$. Let ${\displaystyle T^{(1)},\dots ,T^{(c)}}$ be independent consisting of ${\displaystyle T}$ and ${\displaystyle M=T^{(1)}\bullet \dots \bullet T^{(c)}}$.

Then ${\displaystyle |\|Mx\|_{2}-\|x\|_{2}|<\varepsilon \|x\|_{2}}$ with probability ${\displaystyle 1-\delta }$ for any vector ${\displaystyle x}$ if

${\displaystyle m=(4a)^{2c}\varepsilon ^{-2}\log 1/\delta +(2ae)\varepsilon ^{-1}(\log 1/\delta )^{c}}$.

In particular, if the entries of ${\displaystyle T}$ are ${\displaystyle \pm 1}$ we get ${\displaystyle m=O(\varepsilon ^{-2}\log 1/\delta +\varepsilon ^{-1}({\tfrac {1}{c}}\log 1/\delta )^{c})}$ which matches the normal Johnson Lindenstrauss theorem of ${\displaystyle m=O(\varepsilon ^{-2}\log 1/\delta )}$ when ${\displaystyle \varepsilon }$ is small.

The paper[15] also shows that the dependency on ${\displaystyle \varepsilon ^{-1}({\tfrac {1}{c}}\log 1/\delta )^{c}}$ is necessary for constructions using tensor randomized projections with Gaussian entries.

Because of the exponential dependency on ${\displaystyle c}$ in tensor sketches based on the face-splitting product, a different approach was developed in 2020[15] which applies

${\displaystyle M(x\otimes y\otimes \cdots )=M^{(1)}(x\otimes (M^{(2)}y\otimes \cdots ))}$

We can achieve such an ${\displaystyle M}$ by letting

${\displaystyle M=M^{(c)}(M^{(c-1)}\otimes I_{d})(M^{(c-2)}\otimes I_{d^{2}})\cdots (M^{(1)}\otimes I_{d^{c-1}})}$.

With this method, we only apply the general tensor sketch method to order 2 tensors, which avoids the exponential dependency in the number of rows.

It can be proved[15] that combining ${\displaystyle c}$ dimensionality reductions like this only increases ${\displaystyle \varepsilon }$ by a factor ${\displaystyle {\sqrt {c}}}$.

The fast Johnson–Lindenstrauss transform is a dimensionality reduction matrix

Given a matrix ${\displaystyle M\in \mathbb {R} ^{k\times d}}$, computing the matrix vector product ${\displaystyle Mx}$ takes ${\displaystyle kd}$ time. The Fast Johnson Lindenstrauss Transform (FJLT),[16] was introduced by Ailon and Chazelle in 2006.

A version of this method takes ${\displaystyle M=\operatorname {SHD} }$ where

1. ${\displaystyle D}$ is a diagonal matrix where each diagonal entry ${\displaystyle D_{i,i}}$ is ${\displaystyle \pm 1}$ independently.

The matrix-vector multiplication ${\displaystyle Dx}$ can be computed in ${\displaystyle O(d)}$ time.

1. ${\displaystyle H}$ is a Hadamard matrix, which allows matrix-vector multiplication in time ${\displaystyle O(d\log d)}$
2. ${\displaystyle S}$ is a ${\displaystyle k\times d}$ sampling matrix which is all zeros, except a single 1 in each row.

If the diagonal matrix is replaced by one which has a tensor product of ${\displaystyle \pm 1}$ values on the diagonal, instead of being fully independent, it is possible to compute ${\displaystyle \operatorname {SHD} (x\otimes y)}$ fast.

For an example of this, let ${\displaystyle \rho ,\sigma \in \{-1,1\}^{2}}$ be two independent ${\displaystyle \pm 1}$ vectors and let ${\displaystyle D}$ be a diagonal matrix with ${\displaystyle \rho \otimes \sigma }$ on the diagonal. We can then split up ${\displaystyle \operatorname {SHD} (x\otimes y)}$ as follows:

${\displaystyle {\begin{aligned}&\operatorname {SHD} (x\otimes y)\\&\quad ={\begin{bmatrix}1&0&0&0\\0&0&1&0\\0&1&0&0\end{bmatrix}}{\begin{bmatrix}1&1&1&1\\1&-1&1&-1\\1&1&-1&-1\\1&-1&-1&1\end{bmatrix}}{\begin{bmatrix}\sigma _{1}\rho _{1}&0&0&0\\0&\sigma _{1}\rho _{2}&0&0\\0&0&\sigma _{2}\rho _{1}&0\\0&0&0&\sigma _{2}\rho _{2}\\\end{bmatrix}}{\begin{bmatrix}x_{1}y_{1}\\x_{2}y_{1}\\x_{1}y_{2}\\x_{2}y_{2}\end{bmatrix}}\\[5pt]&\quad =\left({\begin{bmatrix}1&0\\0&1\\1&0\end{bmatrix}}\bullet {\begin{bmatrix}1&0\\1&0\\0&1\end{bmatrix}}\right)\left({\begin{bmatrix}1&1\\1&-1\end{bmatrix}}\otimes {\begin{bmatrix}1&1\\1&-1\end{bmatrix}}\right)\left({\begin{bmatrix}\sigma _{1}&0\\0&\sigma _{2}\\\end{bmatrix}}\otimes {\begin{bmatrix}\rho _{1}&0\\0&\rho _{2}\\\end{bmatrix}}\right)\left({\begin{bmatrix}x_{1}\\x_{2}\end{bmatrix}}\otimes {\begin{bmatrix}y_{1}\\y_{2}\end{bmatrix}}\right)\\[5pt]&\quad =\left({\begin{bmatrix}1&0\\0&1\\1&0\end{bmatrix}}\bullet {\begin{bmatrix}1&0\\1&0\\0&1\end{bmatrix}}\right)\left({\begin{bmatrix}1&1\\1&-1\end{bmatrix}}{\begin{bmatrix}\sigma _{1}&0\\0&\sigma _{2}\\\end{bmatrix}}{\begin{bmatrix}x_{1}\\x_{2}\end{bmatrix}}\,\otimes \,{\begin{bmatrix}1&1\\1&-1\end{bmatrix}}{\begin{bmatrix}\rho _{1}&0\\0&\rho _{2}\\\end{bmatrix}}{\begin{bmatrix}y_{1}\\y_{2}\end{bmatrix}}\right)\\[5pt]&\quad ={\begin{bmatrix}1&0\\0&1\\1&0\end{bmatrix}}{\begin{bmatrix}1&1\\1&-1\end{bmatrix}}{\begin{bmatrix}\sigma _{1}&0\\0&\sigma _{2}\\\end{bmatrix}}{\begin{bmatrix}x_{1}\\x_{2}\end{bmatrix}}\,\circ \,{\begin{bmatrix}1&0\\1&0\\0&1\end{bmatrix}}{\begin{bmatrix}1&1\\1&-1\end{bmatrix}}{\begin{bmatrix}\rho _{1}&0\\0&\rho _{2}\\\end{bmatrix}}{\begin{bmatrix}y_{1}\\y_{2}\end{bmatrix}}.\end{aligned}}}$

In other words, ${\displaystyle \operatorname {SHD} =S^{(1)}HD^{(1)}\bullet S^{(2)}HD^{(2)}}$, splits up into two Fast Johnson–Lindenstrauss transformations, and the total reduction takes time ${\displaystyle O(d_{1}\log d_{1}+d_{2}\log d_{2})}$ rather than ${\displaystyle d_{1}d_{2}\log(d_{1}d_{2})}$ as with the direct approach.

The same approach can be extended to compute higher degree products, such as ${\displaystyle \operatorname {SHD} (x\otimes y\otimes z)}$

Ahle et al.[15] shows that if ${\displaystyle \operatorname {SHD} }$ has ${\displaystyle \varepsilon ^{-2}(\log 1/\delta )^{c+1}}$ rows, then ${\displaystyle |\|\operatorname {SHD} x\|_{2}-\|x\||\leq \varepsilon \|x\|_{2}}$ for any vector ${\displaystyle x\in \mathbb {R} ^{d^{c}}}$ with probability ${\displaystyle 1-\delta }$, while allowing fast multiplication with degree ${\displaystyle c}$ tensors.

Jin et al.,[17] the same year, showed a similar result for the more general class of matrices call RIP, which includes the subsampled Hadamard matrices. They showed that these matrices allow splitting into tensors provided the number of rows is ${\displaystyle \varepsilon ^{-2}(\log 1/\delta )^{2c-1}\log d}$. In the case ${\displaystyle c=2}$ this matches the previous result.

These fast constructions can again be combined with the recursion approach mentioned above, giving the fastest overall tensor sketch.

Kernel methods are popular in machine learning as they give the algorithm designed the freedom to design a "feature space" in which to measure the similarity of their data points. A simple kernel-based binary classifier is based on the following computation:

${\displaystyle {\hat {y}}(\mathbf {x'} )=\operatorname {sgn} \sum _{i=1}^{n}y_{i}k(\mathbf {x} _{i},\mathbf {x'} ),}$

where ${\displaystyle \mathbf {x} _{i}\in \mathbb {R} ^{d}}$ are the data points, ${\displaystyle y_{i}}$ is the label of the ${\displaystyle i}$th point (either −1 or +1), and ${\displaystyle {\hat {y}}(\mathbf {x'} )}$ is the prediction of the class of ${\displaystyle \mathbf {x'} }$. The function ${\displaystyle k:\mathbb {R} ^{d}\times \mathbb {R} ^{d}\to \mathbb {R} }$ is the kernel. Typical examples are the radial basis function kernel, ${\displaystyle k(x,x')=\exp(-\|x-x'\|_{2}^{2})}$, and polynomial kernels such as ${\displaystyle k(x,x')=(1+\langle x,x'\rangle )^{2}}$.

When used this way, the kernel method is called "implicit". Sometimes it is faster to do an "explicit" kernel method, in which a pair of functions ${\displaystyle f,g:\mathbb {R} ^{d}\to \mathbb {R} ^{D}}$ are found, such that ${\displaystyle k(x,x')=\langle f(x),g(x')\rangle }$. This allows the above computation to be expressed as

${\displaystyle {\hat {y}}(\mathbf {x'} )=\operatorname {sgn} \sum _{i=1}^{n}y_{i}\langle f(\mathbf {x} _{i}),g(\mathbf {x'} )\rangle =\operatorname {sgn} \left\langle \left(\sum _{i=1}^{n}y_{i}f(\mathbf {x} _{i})\right),g(\mathbf {x'} )\right\rangle ,}$

where the value ${\displaystyle \sum _{i=1}^{n}y_{i}f(\mathbf {x} _{i})}$ can be computed in advance.

The problem with this method is that the feature space can be very large. That is ${\displaystyle D>>d}$. For example, for the polynomial kernel ${\displaystyle k(x,x')=\langle x,x'\rangle ^{3}}$ we get ${\displaystyle f(x)=x\otimes x\otimes x}$ and ${\displaystyle g(x')=x'\otimes x'\otimes x'}$, where ${\displaystyle \otimes }$ is the tensor product and ${\displaystyle f(x),g(x')\in \mathbb {R} ^{D}}$ where ${\displaystyle D=d^{3}}$. If ${\displaystyle d}$ is already large, ${\displaystyle D}$ can be much larger than the number of data points (${\displaystyle n}$) and so the explicit method is inefficient.

The idea of tensor sketch is that we can compute approximate functions ${\displaystyle f',g':\mathbb {R} ^{d}\to \mathbb {R} ^{t}}$ where ${\displaystyle t}$ can even be smaller than ${\displaystyle d}$, and which still have the property that ${\displaystyle \langle f'(x),g'(x')\rangle \approx k(x,x')}$.

This method was shown in 2020[15] to work even for high degree polynomials and radial basis function kernels.

Assume we have two large datasets, represented as matrices ${\displaystyle X,Y\in \mathbb {R} ^{n\times d}}$, and we want to find the rows ${\displaystyle i,j}$ with the largest inner products ${\displaystyle \langle X_{i},Y_{j}\rangle }$. We could compute ${\displaystyle Z=XY^{T}\in \mathbb {R} ^{n\times n}}$ and simply look at all ${\displaystyle n^{2}}$ possibilities. However, this would take at least ${\displaystyle n^{2}}$ time, and probably closer to ${\displaystyle n^{2}d}$ using standard matrix multiplication techniques.

The idea of Compressed Matrix Multiplication is the general identity

${\displaystyle XY^{T}=\sum _{i=1}^{d}X_{i}\otimes Y_{i}}$

where ${\displaystyle \otimes }$ is the tensor product. Since we can compute a (linear) approximation to ${\displaystyle X_{i}\otimes Y_{i}}$ efficiently, we can sum those up to get an approximation for the complete product.

Tensor sketches can be used to decrease the number of variables needed when implementing Bilinear Pooling in a neural network.

Bilinear pooling is the technique of taking two input vectors, ${\displaystyle x,y}$ from different sources, and using the tensor product ${\displaystyle x\otimes y}$ as the input layer to a neural network.

In[19] the authors considered using tensor sketch to reduce the number of variables needed.

In 2017 another paper[20] takes the FFT of the input features, before they are combined using the element-wise product. This again corresponds to the original tensor sketch.