Low-rank matrix approximations

${\displaystyle {\hat {K}}}$ is a kernel matrix for some kernel method. Consider the first ${\displaystyle q<n}$ points in the training set. Then there exists a matrix ${\displaystyle {\tilde {K}}}$ of rank ${\displaystyle q}$:

${\displaystyle {\tilde {K}}={\hat {K}}_{n,q}{\hat {K}}_{q}^{-1}{\hat {K}}_{n,q}^{\text{T}}}$ , where

${\displaystyle ({\hat {K}}_{q})_{i,j}=K(x_{i},x_{j}),i,j=1,\dots ,q}$ ,

${\displaystyle {\hat {K}}_{q}}$ is invertible matrix

and

${\displaystyle ({\hat {K}}_{n,q})_{i,j}=K(x_{i},x_{j}),i=1,\dots ,n{\text{ and }}j=1,\dots ,q.}$

Applying singular-value decomposition (SVD) to matrix ${\displaystyle A}$ with dimensions ${\displaystyle p\times m}$ produces a singular system consisting of singular values ${\displaystyle \{\sigma _{j}\}_{j=1}^{k},{\text{ }}(\sigma _{j}>0{\text{ }}\forall j=1,\dots ,k),}$ vectors ${\displaystyle \{v_{j}\}_{j=1}^{m}\in \mathbb {C} ^{m}}$ and ${\displaystyle \{u_{j}\}_{j=1}^{p}\in \mathbb {C} ^{p}}$ such that they form orthonormal bases of ${\displaystyle \mathbb {C} ^{m}}$ and ${\displaystyle \mathbb {C} ^{p}}$ respectively:

$${\displaystyle {\begin{cases}A^{\text{T}}Av_{j}=\sigma _{j}v_{j},{\text{ }}j=1,\dots ,k,\\A^{\text{T}}Av_{j}=0,{\text{ }}j=k+1,\dots ,m,\\AA^{\text{T}}u_{j}=\sigma _{j}u_{j},{\text{ }}j=1,\dots ,k,\\AA^{\text{T}}u_{j}=0,{\text{ }}j=k+1,\dots ,p.\end{cases}}}$$

If ${\displaystyle U}$ and ${\displaystyle V}$ are matrices with ${\displaystyle u}$'s and ${\displaystyle v}$'s in the columns and ${\displaystyle \Sigma }$ is a diagonal ${\displaystyle p\times m}$ matrix having singular values ${\displaystyle \sigma _{i}}$ on the first ${\displaystyle k}$-entries on the diagonal (all the other elements of the matrix are zeros):

$${\displaystyle {\begin{cases}Av_{j}={\sqrt {\sigma _{j}}}u_{j},{\text{ }}j=1,\dots ,k,\\Av_{j}=0,{\text{ }}j=k+1,\dots ,m,\\A^{\text{T}}u_{j}={\sqrt {\sigma _{j}}}v_{j},{\text{ }}j=1,\dots ,k,\\A^{\text{T}}u_{j}=0,{\text{ }}j=k+1,\dots ,p,\end{cases}}}$$

then the matrix ${\displaystyle A}$ can be rewritten as:[4]

$${\displaystyle A=U\Sigma ^{1/2}V^{\text{T}}.}$$

• ${\displaystyle {\hat {X}}}$ is ${\displaystyle n\times D}$ data matrix
• ${\displaystyle {\hat {K}}={\hat {X}}{\hat {X}}^{\text{T}}}$
• ${\displaystyle {\hat {C}}={\hat {X}}^{\text{T}}{\hat {X}}}$

Applying singular-value decomposition to these matrices:

$${\displaystyle {\hat {X}}={\hat {U}}{\hat {\Sigma }}^{1/2}{\hat {V}}^{\text{T}},{\text{ }}{\hat {K}}={\hat {U}}{\hat {\Sigma }}{\hat {U}}^{T},{\text{ }}{\hat {C}}={\hat {V}}{\hat {\Sigma }}{\hat {V}}^{\text{T}}.}$$

• ${\displaystyle {\hat {X}}_{q}}$ is the ${\displaystyle q\times D}$-dimensional matrix consisting of the first ${\displaystyle q}$ rows of matrix ${\displaystyle {\hat {X}}}$
• ${\displaystyle {\hat {K}}_{q}={\hat {X}}_{q}{\hat {X}}_{q}^{\text{T}}}$
• ${\displaystyle {\hat {C}}_{q}={\hat {X}}_{q}^{\text{T}}{\hat {X}}_{q}}$

Applying singular-value decomposition to these matrices:

${\displaystyle {\hat {X}}_{q}={\hat {U}}_{q}{\hat {\Sigma }}_{q}^{1/2}{\hat {V}}_{q}^{\text{T}},{\text{ }}{\hat {K}}_{q}={\hat {U}}_{q}{\hat {\Sigma }}_{q}{\hat {U}}_{q}^{T},{\text{ }}{\hat {C}}_{q}={\hat {V}}_{q}{\hat {\Sigma }}_{q}{\hat {V}}_{q}^{\text{T}}.}$

Since ${\displaystyle {\hat {U}},{\text{ }}{\hat {V}},{\hat {U}}_{q}{\text{ and }}{\hat {V}}_{q}}$ are orthogonal matrices,

$${\displaystyle {\hat {U}}={\hat {X}}{\hat {V}}{\hat {\Sigma }}^{-1/2},{\text{ }}{\hat {V}}_{q}={\hat {X}}_{q}^{\text{T}}{\hat {U}}_{q}{\hat {\Sigma }}_{q}^{-1/2}.}$$

Replacing ${\displaystyle {\hat {V}},{\hat {\Sigma }}}$ by ${\displaystyle {\hat {V}}_{q}}$ and ${\displaystyle {\hat {\Sigma }}_{q}}$, an approximation for ${\displaystyle {\hat {U}}}$ can be obtained:

${\displaystyle {\tilde {U}}={\hat {X}}{\hat {X}}_{q}^{\text{T}}{\hat {U}}_{q}{\hat {\Sigma }}_{q}^{-1}}$ ( ${\displaystyle {\tilde {U}}}$ is not necessarily an orthogonal matrix).

However, defining ${\displaystyle {\tilde {K}}={\tilde {U}}{\hat {\Sigma }}_{q}{\tilde {U}}^{\text{T}}}$, it can be computed the next way:

$${\displaystyle {\begin{aligned}{\tilde {K}}={\tilde {U}}{\hat {\Sigma }}_{q}{\tilde {U}}^{\text{T}}={\hat {X}}{\hat {X}}_{q}^{\text{T}}{\hat {U}}_{q}{\hat {\Sigma }}_{q}^{-1}{\hat {\Sigma }}_{q}({\hat {X}}{\hat {X}}_{q}^{\text{T}}{\hat {U}}_{q}{\hat {\Sigma }}_{q}^{-1})^{\text{T}}\\={\hat {X}}{\hat {X}}_{q}^{\text{T}}{\big \{}{\hat {U}}_{q}({\hat {\Sigma }}_{q}^{-1})^{\text{T}}{\hat {U}}_{q}^{\text{T}}{\big \}}({\hat {X}}{\hat {X}}_{q}^{\text{T}})^{\text{T}}\\\end{aligned}}}$$

By the characterization for orthogonal matrix ${\displaystyle {\hat {U}}_{q}}$: equality ${\displaystyle ({\hat {U}}_{q})^{\text{T}}=({\hat {U}}_{q})^{-1}}$ holds. Then, using the formula for the inverse of matrix multiplication ${\displaystyle (AB)^{-1}=B^{-1}A^{-1}}$ for invertible matrices ${\displaystyle A}$ and ${\displaystyle B}$, the expression in braces can be rewritten as:

$${\displaystyle {\hat {U}}_{q}({\hat {\Sigma }}_{q}^{-1})^{\text{T}}{\hat {U}}_{q}^{\text{T}}=({\hat {U}}_{q}{\hat {\Sigma }}_{q}^{\text{T}}{\hat {U}}_{q}^{\text{T}})^{-1}=({\hat {K}}_{q})^{-1}.}$$

Then the expression for ${\displaystyle {\tilde {K}}}$:

$${\displaystyle {\tilde {K}}=({\hat {X}}{\hat {X}}_{q}^{\text{T}}){\hat {K}}_{q}^{-1}({\hat {X}}{\hat {X}}_{q}^{\text{T}})^{\text{T}}.}$$

Defining ${\displaystyle {\hat {K}}_{n,q}={\hat {X}}{\hat {X}}_{q}^{\text{T}}}$, the proof is finished.

For a feature map ${\displaystyle \Phi$ :{\mathcal {X}}\to {\mathcal {F}}} with associated kernel ${\displaystyle K(x,x')=\langle \Phi (x),\Phi (x')\rangle _{\mathcal {F}}}$: equality ${\displaystyle {\hat {K}}={\hat {K}}_{n,q}{\hat {K}}_{q}^{-1}{\hat {K}}_{n,q}^{\text{T}}}$ also follows by replacing ${\displaystyle {\hat {X}}}$ by the operator ${\displaystyle {\hat {\Phi }}:{\mathcal {F}}\to \mathbb {R} ^{n}}$ such that ${\displaystyle \langle {\hat {\Phi }}w\rangle _{i}=\langle \Phi (x_{i}),w\rangle _{\mathcal {F}}}$, ${\displaystyle {\text{ }}i=1,\dots ,n}$, ${\displaystyle w\in {\mathcal {F}}}$, and ${\displaystyle {\hat {X}}_{q}}$ by the operator ${\displaystyle {\hat {\Phi }}_{q}:{\mathcal {F}}\to \mathbb {R} ^{q}}$ such that ${\displaystyle \langle {\hat {\Phi }}_{q}w\rangle _{i}=\langle \Phi (x_{i}),w\rangle _{\mathcal {F}}}$, ${\displaystyle {\text{ }}i=1,\dots ,q}$, ${\displaystyle w\in {\mathcal {F}}}$ . Once again, a simple inspection shows that the feature map is only needed in the proof while the end result only depends on computing the kernel function.

In a vector and kernel notation, the problem of Regularized least squares can be rewritten as:

$${\displaystyle \min _{c\in \mathbb {R} ^{n}}{\frac {1}{n}}\|{\hat {Y}}-{\hat {K}}c\|_{\mathbb {R} ^{n}}^{2}+\lambda \langle c,{\hat {K}}c\rangle _{\mathbb {R} ^{n}}.}$$

Computing the gradient and setting in to 0, the minimum can be obtained:

$${\displaystyle {\begin{aligned}&-{\frac {1}{n}}{\hat {K}}({\hat {Y}}-{\hat {K}}c)+\lambda {\hat {K}}c=0\\\Rightarrow {}&{\hat {K}}({\hat {K}}+\lambda nI)c={\hat {K}}{\hat {Y}}\\\Rightarrow {}&c=({\hat {K}}+\lambda nI)^{-1}{\hat {Y}},{\text{ where }}c\in \mathbb {R} ^{n}\end{aligned}}}$$

The inverse matrix ${\displaystyle ({\hat {K}}+\lambda nI)^{-1}}$ can be computed using Woodbury matrix identity:

$${\displaystyle {\begin{aligned}({\hat {K}}+\lambda nI)^{-1}&={\frac {1}{\lambda n}}\left({\frac {1}{\lambda n}}{\hat {K}}+I\right)^{-1}\\&={\frac {1}{\lambda n}}\left(I+{\hat {K}}_{n,q}({\lambda n}{\hat {K}}_{q})^{-1}{\hat {K}}_{n,q}^{\text{T}}\right)^{-1}\\&={\frac {1}{\lambda n}}\left(I-{\hat {K}}_{n,q}(\lambda n{\hat {K}}_{q}+{\hat {K}}_{n,q}^{\text{T}}{\hat {K}}_{n,q})^{-1}{\hat {K}}_{n,q}^{\text{T}}\right)\end{aligned}}}$$

It has the desired storage and complexity requirements.

Random Fourier features map produces a Monte Carlo approximation to the feature map. The Monte Carlo method is considered to be randomized. These random features consists of sinusoids ${\displaystyle \cos(w^{\text{T}}\mathbf {x} +b)}$ randomly drawn from Fourier transform of the kernel to be approximated, where ${\displaystyle w\in \mathbb {R} ^{d}}$ and ${\displaystyle b\in \mathbb {R} }$ are random variables. The line is randomly chosen, then the data points are projected on it by the mappings. The resulting scalar is passed through a sinusoid. The product of the transformed points will approximate a shift-invariant kernel. Since the map is smooth, random Fourier features work well on interpolation tasks.

A random binning features map partitions the input space using randomly shifted grids at randomly chosen resolutions and assigns to an input point a binary bit string that corresponds to the bins in which it falls. The grids are constructed so that the probability that two points ${\displaystyle \mathbf {x} ,\mathbf {x'} \in \mathbb {R} ^{d}}$ are assigned to the same bin is proportional to ${\displaystyle K(\mathbf {x} ,\mathbf {x'} )}$. The inner product between a pair of transformed points is proportional to the number of times the two points are binned together, and is therefore an unbiased estimate of ${\displaystyle K(\mathbf {x} ,\mathbf {x'} )}$. Since this mapping is not smooth and uses the proximity between input points, Random Binning Features works well for approximating kernels that depend only on the ${\displaystyle L_{1}}$ distance between datapoints.