# Nystroem Method for Kernel Approximation

I have been reading about the Nyström method for low-rank kernel aproximation. This method is implemented in scikit-learn [1] as a method to project data samples to a low-rank approximation of the kernel feature mapping.

To the best of my knowledge, given a training set $\{x_i\}_{i=1}^n$ and a kernel function, it generates a low-rank approximation of the $n \times n$ kernel matrix $K$ by applying SVD to $W$ and $C$.

$K = \left [ \begin{array}{cc} W & K_{21}^T \\ K_{21} & K_{22} \end{array} \right ]$ $C = \left [\begin{array}{cc} W \\ K_{21} \end{array}\right ]$, $W \in \mathbb{R}^{l\times l}$

However, I do not understand how the low-rank approximation of the Kernel matrix can be used to project new samples to the approximated kernel feature space. The papers I have found (e.g. [2]) are not of great help, for they are little didactic.

Also, I am curious about the computational complexity of this method, both in training and test phases.

## 1 Answer

Let's derive the Nyström approximation in a way that should make the answers to your questions clearer.

The key assumption in Nyström is that the kernel function is of rank $$m$$. (Really we assume that it's approximately of rank $$m$$, but for simplicity let's just pretend it's exactly rank $$m$$ for now.) That means that any kernel matrix is going to have rank at most $$m$$, and in particular $$K = \begin{bmatrix} k(x_1, x_1) & \dots & k(x_1, x_n) \\ \vdots & \ddots & \vdots \\ k(x_n, x_1) & \dots & k(x_n, x_n) \end{bmatrix} ,$$ is rank $$m$$. Therefore there are $$m$$ nonzero eigenvalues, and we can write the eigendecomposition of $$K$$ (ignoring the zero eigenvalues) as $$K = U \Lambda U^T$$ with eigenvectors stored in $$U$$, of shape $$n \times m$$, and eigenvalues arranged in $$\Lambda$$, an $$m \times m$$ diagonal matrix whose diagonal entries are all positive.

So, let's pick $$m$$ elements, usually uniformly at random but possibly according to other schemes – all that matters in this simplified version is that $$K_{11}$$ be of rank $$m$$. Once we do, just relabel the points so that we end up with the kernel matrix in blocks: $$K = \begin{bmatrix} K_{11} & K_{21}^T \\ K_{21} & K_{22} \end{bmatrix} ,$$ where we evaluate each entry in $$K_{11}$$ (which is $$m \times m$$) and $$K_{21}$$ ($$(n-m) \times m$$), but don't want to evaluate any entries in $$K_{22}$$.

Now, we can split up the eigendecomposition according to this block structure too: \begin{align} K &= U \Lambda U^T \\&= \begin{bmatrix}U_1 \\ U_2\end{bmatrix} \Lambda \begin{bmatrix}U_1 \\ U_2\end{bmatrix}^T \\&= \begin{bmatrix} U_1 \Lambda U_1^T & U_1 \Lambda U_2^T \\ U_2 \Lambda U_1^T & U_2 \Lambda U_2^T \end{bmatrix} ,\end{align} where $$U_1$$ is $$m \times m$$ and $$U_2$$ is $$(n-m) \times m$$. So, we have $$K_{11} = U_1 \Lambda U_1^T$$ (although notice that this is not an eigendecomposition of $$K_{11}$$), and $$K_{21} = U_2 \Lambda U_1^T$$. If $$U_1$$ is invertible, then $$K_{11}^{-1} = U_1^{-T} \Lambda^{-1} U_1^{-1}$$, and then we have \begin{align} K_{21} K_{11}^{-1} K_{21}^T &= U_2 \Lambda U_1^T U_1^{-T} \Lambda^{-1} U_1^{-1} U_1 \Lambda U_2^T \\&= U_2 \Lambda U_2^T \\&= K_{22} .\end{align} This tells us what $$K_{22}$$ must be under this rank assumption, as long as $$U_1$$ is invertible. $$U_1$$ should be almost surely invertible if $$X$$ is sampled from a continuous distribution; it is possible for it to not be invertible (consider the full $$n \times n$$ eigenvectors being a permuted identity such that $$U_1$$ has an all-zero row), but in that case you can follow another path that I find less intuitive to the same result.

But we don't want to actually construct $$K_{22}$$, because that's expensive. Instead, notice that $$K_{21} K_{11}^{-1} K_{21}^T = \left( K_{21} K_{11}^{-\frac12} \right) \left( K_{21} K_{11}^{-\frac12} \right)^T .$$ We can see then that using the feature matrix $$K_{21} K_{11}^{-\frac12}$$, of shape $$(n-m) \times m$$, corresponds to these imputed kernel values. If we use $$K_{11} K_{11}^{-\frac12} = K_{11}^{\frac12}$$ for the points in the first partition, we have a set of $$m$$-dimensional features $$\Phi = \begin{bmatrix} K_{11}^{\frac12} \\ K_{21} K_{11}^{-\frac12} \end{bmatrix} .$$ We can just quickly verify that $$\Phi$$ corresponds to the correct kernel matrix: \begin{align} \Phi \Phi^T &= \begin{bmatrix} K_{11}^{\frac12} \\ K_{21} K_{11}^{-\frac12} \end{bmatrix} \begin{bmatrix} K_{11}^{\frac12} \\ K_{21} K_{11}^{-\frac12} \end{bmatrix}^T \\&=\begin{bmatrix} K_{11}^{\frac12} K_{11}^{\frac12} & K_{11}^{\frac12} K_{11}^{-\frac12} K_{21}^T \\ K_{21} K_{11}^{-\frac12} K_{11}^{\frac12} & K_{21} K_{11}^{-\frac12} K_{11}^{-\frac12} K_{21}^T \end{bmatrix} \\&=\begin{bmatrix} K_{11} & K_{21}^T \\ K_{21} & K_{21} K_{11}^{-1} K_{21}^T \end{bmatrix} \\&= K .\end{align}

So, all we need to do is train our regular learning model with the $$m$$-dimensional features $$\Phi$$. This will be exactly the same (under the assumptions we've made) as the kernelized version of the learning problem with $$K$$.

Now, for an individual data point $$x$$, the features in $$\Phi$$ correspond to $$\phi(x) = \begin{bmatrix} k(x, x_1) & \dots & k(x, x_m) \end{bmatrix} K_{11}^{-\frac12} .$$ This vector is just the relevant row of either $$K_{11}$$ or $$K_{21}$$, so $$\phi(x)$$ that multiplies that row by $$K_{11}^{-\frac12}$$ is indeed the corresponding row of $$\Phi$$.

So...this is still true for a novel test point. You just do the same thing: $$\Phi_\text{test} = K_{\text{test},1} K_{11}^{-\frac12} .$$ Because we assumed the kernel is rank $$m$$, the matrix $$\begin{bmatrix}K_{\text{train}} & K_{\text{train,test}} \\ K_{\text{test,train}} & K_{\text{test}} \end{bmatrix}$$ is also of rank $$m$$, and the reconstruction of $$K_\text{test}$$ is still exact by exactly the same logic as for $$K_{22}$$.

Above, we assumed that the kernel matrix $$K$$ was *exactly* rank $$m$$. This is not usually going to be the case; for a Gaussian kernel, for example, $$K$$ is always rank $$n$$, but the smaller eigenvalues typically drop off pretty quickly, so it's going to be *close to* a matrix of rank $$m$$, and our reconstructions of $$K_{21}$$ or $$K_{\text{test},1}$$ are going to be *close to* the true values (but not exactly the same). They'll be better reconstructions the closer the eigenspace of $$K_{11}$$ gets to that of $$K$$ overall, which is why choosing the right $$m$$ points is important in practice.

If $$K_{11}$$ has any zero eigenvalues, you can replace inverses with pseudoinverses and everything still works; you just replace $$K_{21}$$ in the reconstruction with $$K_{21} K_{11}^\dagger K_{11}$$.

You can use the SVD instead of the eigendecomposition if you'd like; since $$K$$ is psd, they're the same thing, but the SVD might be a little more robust to slight numerical error in the kernel matrix, and that's what scikit-learn does. scikit-learn's actual implementation does this, though it uses $$\max(\lambda_i, 10^{-12})$$ in the inverse instead of the pseudoinverse. This is related to avoiding numerical error, but I'm not sure why they make small eigenvalues $$10^{12}$$ instead of zero like np.linalg.pinv does.

• Thank you for your great answer. It was of great help. However, I still do not fully understand how the computation of $K_{11}^{-1/2}$ is carried out in scikit-learn. Can you provide some reference of how this is achieved through SVD? Commented Feb 12, 2017 at 11:44
• When $A$ is positive semidefinite, the eigendecomposition $U \Lambda U^T$ coincides with the SVD. scikit-learn, because due to numerical error $A$ might be slightly non-psd, instead computes $U \Sigma V^T$, and uses $A^{-\frac12} = V \Sigma^{-\frac12} V^T$, so that $A$'s features become $A V \Sigma^{-\frac12} V^T = U \Sigma V^T V \Sigma^{-\frac12} V^T = U \Sigma^{\frac12} V^T = A^{\frac12}$. It's the same thing, basically. Commented Feb 12, 2017 at 12:07
• Whoops, sorry, yeah they use $U \Sigma^{-\frac12} V^T = K^{-\frac12}$. It all doesn't really matter since $U \approx V$, but since they do the transpose the features for $K_{11}$ end up as $U\Sigma V^T V \Sigma^{-\frac12} U^T = U \Sigma^{\frac12} U^T$. Commented Feb 12, 2017 at 12:43
• Raising a diagonal matrix to a power is the same as raising each element to a power, and $x^{-\frac12} = 1 / \sqrt x$. In numpy broadcasting notation, elementwise multiplication by a vector is the same as right-multiplying by a diagonal matrix. Also, that code uses $V$ to mean what I was calling $V^T$. Commented Feb 12, 2017 at 12:54
• "right-multiply both sides by ..." - Why would $U_1$ be an orthogonal matrix? I get that $U$ is orthogonal, but $U_1$ is a sub-matrix of $U$. Commented Nov 9, 2023 at 18:31