Improved Subsampled Randomized Hadamard Transform for Linear SVM

Given a data matrix $\mathbf{X}\in\mathbb{R}^{n\times d}$, random projection methods reduce the dimensionality of $\mathbf{X}$ by multiplying it by a random orthonormal matrix $\mathbf{R}\in\mathbb{R}^{d\times r}$ with parameter $r\ll d$. The projected data in the low-dimensional space is

Note that the matrix $\mathbf{R}$ is randomly generated and is independent of the input data $\mathbf{X}$. The theoretical foundation of random projection is the Johnson-Lindenstrauss Lemma (JL lemma) as shown in following,

JL lemma proves that in Euclidean space, high-dimensional data can be randomly projected into lower dimensional space while the pairwise distances between all data points are well preserved with high probability. There are several different ways to construct the random projection matrix $\mathbf{R}$. In this paper, we focus on SRHT which uses a structured orthonormal matrix for random projection. The details about SRHT are as follows.

For $d=2^{q}$where $q$ is any positive integer¹¹1We can ensure this by padding zeros to original data , SRHT defines a $d\times r$ matrix as:

where

• •

  $\mathbf{D}\in\mathbb{R}^{d\times d}$ is a diagonal matrix whose elements are independent random signs {1,$-1$};

• •

  $\mathbf{H}\in\mathbb{R}^{d\times d}$ is a normalized Walsh-Hadamard matrix. The Walsh-Hadamard matrix is defined recursively as:
  $\mathbf{H}_{d}=\left[\begin{array}[]{cc}\mathbf{H}_{d/2}&\mathbf{H}_{d/2}\\ \mathbf{H}_{d/2}&-\mathbf{H}_{d/2}\end{array}\right]$ with $\mathbf{H}_{2}=\left[\begin{array}[]{cc}1&1\\ 1&-1\end{array}\right],$ and $\mathbf{H}=\frac{1}{\sqrt{d}}\mathbf{H}_{d}\in\mathbb{R}^{d\times d}$;

• •

  $\mathbf{S}\in\mathbb{R}^{d\times r}$ is a subset of randomly sampling $r$ columns from the $d\times d$ identity matrix. The purpose of multiplying $\mathbf{S}$ is to uniformly sample $r$ columns from the rotated data matrix $\mathbf{X}_{r}=\mathbf{XDH}$.

There are several advantages of using SRHT for random projection. Firstly, we do not need to explicitly represent $\mathbf{H}$ in memory and only need to store a diagonal random sign matrix $\mathbf{D}$ and a sampling matrix $\mathbf{S}$ with $r$ non-zero values. Therefore, the memory cost of storing $\mathbf{R}$ is only $O(d+r)$. Secondly, due to the recursive structure of Hadamard matrix $\mathbf{H}$, matrix multiplication $\mathbf{XR}$ only takes $O(nd\log(d))$ time by using Fast Fourier Transform (FFT) for matrix multiplication (?). ? (?) further improve the time complexity of SRHT to $O(nd\log(r))$ if only $r$ columns in $\mathbf{X}_{r}$ are needed.

Despite the advantages of SRHT as mentioned in the previous section, one limitation of SRHT is that the random matrix $\mathbf{R}$ is constructed without considering any specific properties of input data $\mathbf{X}$. Based on the definition in (2), SRHT works by first rotating the input data matrix $\mathbf{X}$ by multiplying it with $\mathbf{DH}$ and then uniformly random sampling $r$ columns from the rotated matrix $\mathbf{X}_{r}=\mathbf{XDH}$. Since this uniform column sampling method does not consider the underlying data distribution of $\mathbf{X}_{r}$, it would result in low and unstable accuracy when used as a dimensionality reduction method in a specific classification problem. To illustrate the limitation of uniformly random sampling step in SRHT, we generate a simple two-dimensional synthetic data as shown in Figure 1(a). The first class (i.e., class label ‘+’) was generated from $N$($\mu_{1}=[3,3]$, $\Sigma=\left[\begin{array}[]{cc}1&1.75\\ 1.75&1\end{array}\right]$). The second class (i.e., class label ‘$\bullet$’) was generated from $N$ ($\mu_{2}=[-3,-3]$, $\Sigma=\left[\begin{array}[]{cc}1&1.75\\ 1.75&1\end{array}\right]$). After multiplying with matrix $\mathbf{DH}$, the original data was rotated into new feature space as shown in Figure 1(b). In this simple illustration example, SRHT might choose the feature that can not distinguish these two classes (i.e, the vertical feature in Figure 1(b)) and then leads to low classification accuracy. In Figure 1(c), we shown the classification accuracy of using SRHT on a real-life dataset for 15 different runs. It clearly shows that this data-independent random projection method could result in low and unstable performance on this synthetic dataset. In the following sections, we will describe our proposed methods to construct the data-dependent sampling matrix $\mathbf{S}$ in (2) by exploiting the underlying data properties from both unsupervised and supervised perspectives.

First, we present our proposed non-uniform sampling methods from unsupervised perspective. In this paper, we study SRHT in the context of using linear SVM for classification. Assume we are given a training data set $D=\{\mathbf{x}_{i},y_{i}\}_{i=1}^{n}$, where $\mathbf{x}_{i}$ is a $d$-dimensional input feature vector, $y_{i}$ is its corresponding class label. The dual problem of SVM can be written in matrix format as shown in following,

where $\mathbf{X}\in\mathbb{R}^{n\times d}$ is the input feature matrix, and $\mathbf{Y}$ is a diagonal matrix where the elements in the diagonal are class labels, $\mathbf{Y}_{ii}=y_{i}$, $C$ is the regularization parameter.

By applying SRHT to project original $\mathbf{X}$ from $d$-dimensional space into $r$-dimensional space, i.e., $\widetilde{\mathbf{X}}=\sqrt{\frac{d}{r}}\mathbf{XDHS}$, the SVM optimization problem on projected data $\widetilde{\mathbf{X}}$ will be

Comparing the dual problem (3) in the original space and the dual problem (4) in the projected low-dimensional space, the only difference is that $\mathbf{X}\mathbf{X}^{T}$ in (3) is replaced by $\widetilde{\mathbf{X}}\widetilde{\mathbf{X}}^{T}$ in (4). The constraints on optimal solutions (i.e., $\boldsymbol{\alpha}^{*}$ of (3) and $\boldsymbol{\widetilde{\alpha}}^{*}$ of (4) ) do not depend on input data matrix. Therefore, $\boldsymbol{\widetilde{\alpha}}^{*}$ is expected to be close to $\boldsymbol{\alpha}^{*}$ when $\widetilde{\mathbf{X}}\widetilde{\mathbf{X}}^{T}$ is close to $\mathbf{X}\mathbf{X}^{T}$.

Let us rewrite $\widetilde{\mathbf{X}}=\mathbf{X}_{r}\mathbf{S}\sqrt{\frac{d}{r}}$ where $\mathbf{X}_{r}=\mathbf{XDH}$, we can view the effect of $\mathbf{X}_{r}\mathbf{S}\sqrt{\frac{d}{r}}$ is to uniformly sample $r$ columns from $\mathbf{X}_{r}$ and then re-scale each selected column by a scaling factor $\sqrt{\frac{d}{r}}$. It corresponds to use uniform column sampling for approximating matrix multiplication (?). It has been proved in (?) that ($i,j$)-th element in $\widetilde{\mathbf{X}}\widetilde{\mathbf{X}}^{T}$ equals to the ($i,j$)-th element of exact product $\mathbf{X}_{r}\mathbf{X}_{r}^{T}$ in expectation regardless of the sampling probabilities. However, the variance of the ($i,j$)-th element in $\widetilde{\mathbf{X}}\widetilde{\mathbf{X}}^{T}$ depends on the choice of sampling probabilities. Therefore, simply using uniform sampling would result in large variance of each element in $\widetilde{\mathbf{X}}\widetilde{\mathbf{X}}^{T}$ and then cause a large approximation error between $\widetilde{\mathbf{X}}\widetilde{\mathbf{X}}^{T}$ and $\mathbf{X}_{r}\mathbf{X}_{r}^{T}$.

Motivated by this observation, we first propose to improve SRHT by employing importance sampling for randomized matrix multiplication. Specifically, we would like to design a method to construct $\widetilde{\mathbf{S}}$ based on specific data properties of $\mathbf{X}_{r}$ instead of using random sampling matrix $\mathbf{S}$ in original SRHT. The idea is based on importance sampling for approximating matrix multiplication: (1) we assign each column $i$ in $\mathbf{X}_{r}$ a probability $p_{i}$ such that $p_{i}\geq 0$ and $\sum_{i=1}^{d}p_{i}=1$; (2) we then select $r$ columns from $\mathbf{X}_{r}$ based on the probability distribution $\{p_{i}\}_{i=1}^{d}$ to form $\widetilde{\mathbf{X}}$; (3) a column in $\widetilde{\mathbf{X}}$ is formed as $\sqrt{\frac{1}{rp_{i}}}\mathbf{X}_{r({:,i}})$ if this column is chosen from the $i$-th column in $\mathbf{X}_{r}$ where $\sqrt{\frac{1}{rp_{i}}}$ is the re-scaling factor. In other word, the matrix $\widetilde{\mathbf{S}}\in\mathbb{R}^{d\times r}$ is constructed as: for every column $j$ in $\widetilde{\mathbf{S}}$ only one entry $i$ ($i$ is a random number from $\{1,\dots,d\}$) is selected by following the probability distribution $\{p_{i}\}_{i=1}^{d}$ and a non-zero value is assigned to this entry. The non-zero value is set as

In the uniform sampling setting, $p_{i}=\frac{1}{d}$, therefore $\widetilde{\mathbf{S}}=\sqrt{\frac{d}{r}}\mathbf{S}$ which is consistent with the setting in original SRHT as shown in (2).

Now, the question is what choice of probability $\{p_{i}\}_{i=1}^{d}$ can optimally reduce the expected approximation error between $\widetilde{\mathbf{X}}\widetilde{\mathbf{X}}^{T}$ and $\mathbf{X}\mathbf{X}^{T}$. We use Frobenius norm to measure the approximation error between $\widetilde{\mathbf{X}}\widetilde{\mathbf{X}}^{T}$ and $\mathbf{X}\mathbf{X}^{T}$.

Let us use $\mathbf{X}_{r(:,j)}$ to denote the $j$-th column of $\mathbf{X}_{r}$, as shown in the Proposition 1, the optimal sampling probability to minimize $\mathbf{E}[\|\widetilde{\mathbf{X}}\widetilde{\mathbf{X}}^{T}-\mathbf{X}\mathbf{X}^{T}\|_{F}]$ is

To prove Proposition 1, we first use the fact $\mathbf{X}\mathbf{X}^{T}=\mathbf{X}_{r}\mathbf{X}_{r}^{T}$ since $\mathbf{D}\mathbf{H}\mathbf{H}^{T}\mathbf{D}^{T}=\mathbf{I}$. Then, we apply the Lemma 4 from (?) to obtain the optimal sampling probability.

In addition, we also propose a deterministic top-$r$ sampling, which select the $r$ columns with largest Euclidean norms from $\mathbf{X}_{r}$. As shown in the Proposition 2, this deterministic top-$r$ sampling method directly optimizes an upper bound of $\|\widetilde{\mathbf{X}}\widetilde{\mathbf{X}}^{T}-\mathbf{X}\mathbf{X}^{T}\|_{F}$ without using re-scaling factor.

In this section, we describe our proposed non-uniform sampling method from supervised perspective. We propose to construct the sampling matrix $\widetilde{\mathbf{S}}$ by incorporating label information. Our proposed method is based on the idea of metric learning (?). The idea is to select $r$ columns from $\mathbf{X}_{r}$ for improving the inter-class separability along with the intra-class compactness as used in (?).

Let $\mathbf{z}_{i}$ be new feature representation of the $i$-th data sample in $\mathbf{X}_{r}$. Then, we measure the interclass separability by the sum of squared pair-wise distances between $\mathbf{z}_{i}$’s from different classes as $\sum_{y_{i}\neq y_{j}}\|\mathbf{z}_{i}-\mathbf{z}_{j}\|^{2}$. Similarly, we measure the intra-class tightness by the sum of squared pair-wise distances between $\mathbf{z}_{i}$’s from the same class as $\sum_{y_{i}=y_{j}}\|\mathbf{z}_{i}-\mathbf{z}_{j}\|^{2}$. We formulate our objective as a combination of these two terms,

where $\mathbf{A}_{ij}=1$ if $y_{i}=y_{j}$, and $\mathbf{A}_{ij}=-a$ if $y_{i}\neq y_{j}$, parameter $a$ is used to balance the tradeoff between inter-class separability and intra-class tightness, $\mathbf{L}$ is the Laplacian matrix of $\mathbf{A}$, defined as $\mathbf{L}=\mathbf{D}-\mathbf{A}$, and $\mathbf{D}$ is a diagonal degree matrix of $\mathbf{A}$ such that $\mathbf{D}_{ii}=\sum_{j}\mathbf{A}_{ij}$. $\widetilde{\mathbf{S}}$ is the sampling matrix whose elements are 0 or 1. We want to learn $\widetilde{\mathbf{S}}$ by minimizing (8).

Without losing any information, let us rewrite our variable $\widetilde{\mathbf{S}}$ as a $d\times d$ diagonal matrix $\mathbf{P}$ where the diagonal element $\mathbf{P}_{ii}\in\{0,1\}$. $\mathbf{P}_{ii}=1$ means the $i$-th column in $\mathbf{X}_{r}$ is selected and $\mathbf{P}_{ii}=0$ means the $i$-th column is not selected. Then $\widetilde{\mathbf{S}}$ can be viewed as a submatrix of $\mathbf{P}$ that selects the columns with $\mathbf{P}_{ii}=1$ from $\mathbf{P}$.

It is straight-forward to verify that $\widetilde{\mathbf{S}}\widetilde{\mathbf{S}}^{\top}$ in (8) equals to $\mathbf{P}$. Therefore, our optimization is formulated as

Let vector $\mathbf{b}$ denote the diagonal of the matrix ${\mathbf{X}_{r}}^{\top}\mathbf{L}\mathbf{X}_{r}$, i.e $b_{i}=({\mathbf{X}_{r}}^{\top}\mathbf{L}\mathbf{X}_{r})_{ii}$, vector $\boldsymbol{p}$ denote the diagonal of $\mathbf{P}$. Then, the objective (9) is simplified to

The optimal solution of (10) is just to select the largest $r$ elements in $\mathbf{b}_{i}$ and set the corresponding entries in $\mathbf{p}_{i}$ to 1. Then $\widetilde{\mathbf{S}}$ can easily obtained based on non-zero values in $\mathbf{p}$.

We summarize our proposed algorithm in algorithm 1 and name it as Improved Subsampled Randomized Hadamard Transform (ISRHT). With respect to training time complexity, step 2 needs $O(ndlog(d))$ time by using FFT. In step 3, for unsupervised column selection, we need $O(nd)$ time to compute the Euclidean norm for each column. For supervised column selection, $b_{i}$ in (10) are computed as $b_{i}=({\mathbf{X}_{r}}^{\top}\mathbf{L}\mathbf{X}_{r})_{ii}$ and the Laplacian matrix L is equal to $\mathbf{D-yy^{T}}$, therefore, the required time is still $O(nd)$. Step 4 needs $O(r)$ time and step 5 needs $O(tnr)$ where $t$ is the number of iterations needed for training linear SVM. Therefore, the overall time complexity of our proposed algorithm is $O(ndlog(d)+nd+tnr)$. With respect to prediction time complexity, our proposed method needs $O(d\log(d)+r)$ time to make a prediction for a new test sample. With respect to space complexity, $\mathbf{D}$ is a diagonal matrix, matrix $\widetilde{\mathbf{S}}$ is a sampling matrix with $r$ non-zero entries. Therefore, the required space for storing the projection matrix $\mathbf{R}$ of our proposed methods is $O(d+r)$. Our proposed methods achieve very efficient time and space complexities for model deployment. We compare the time and space complexity of different random projection methods in Table 1.

In this section, we describe an extension of our algorithm on high-dimensional sparse data. One disadvantage of applying Hadamard transform on sparse data is that it produces a dense matrix. As shown in the step 2 in Algorithm 1, the memory cost is increased from $O(nnz(\mathbf{X}))$ to $O(nd)$. To tackle this challenge, we use a combination of sparse embedding and SRHT as proposed in (?) for memory-efficient projection. We improve their work by replacing the original SRHT by our proposed ISRHT methods. We expect that our proposed data-dependent projection methods can obtain higher accuracies than data-independent projection methods on high-dimensional sparse data.

Specifically, we first use sparse embedding matrix $\mathbf{R}_{sparse}\in\mathbb{R}^{d\times r^{\prime}}$ to project $\mathbf{X}$ into $r^{\prime}$-dimensional space and then use our improved SRHT methods to project into $r$ dimensional space. The sparse embedding matrix $\mathbf{R}_{sparse}$ is generated as follows, for each row in random matrix $\mathbf{R}_{sparse}$, randomly selected one column and assign either 1 or $-1$ with probability 0.5 to this entry. All other entries are 0s. Due to the sparse structure of $\mathbf{R}_{sparse}$, matrix multiplication $\mathbf{X}\mathbf{R}_{sparse}$ takes $O(nnz(\mathbf{X}))$ time. By doing this, the time complexity of the proposed method on high-dimensional sparse data was reduced to $O(nnz(\mathbf{X})+nr^{\prime}log(r^{\prime})+nr^{\prime})$ and the memory cost was reduced to $O(nnz(\mathbf{X})+nr^{\prime})$. In our experimental setting, we set $r^{\prime}=2r$ as suggested in (?).