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] http://scikit-learn.org/stable/modules/kernel_approximation.html#nystroem-kernel-approx

[2] http://www.jmlr.org/papers/volume13/kumar12a/kumar12a.pdf


1 Answer 1


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.

  • 1
    $\begingroup$ 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? $\endgroup$ Commented Feb 12, 2017 at 11:44
  • 1
    $\begingroup$ 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. $\endgroup$
    – Danica
    Commented Feb 12, 2017 at 12:07
  • 1
    $\begingroup$ 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$. $\endgroup$
    – Danica
    Commented Feb 12, 2017 at 12:43
  • 1
    $\begingroup$ 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$. $\endgroup$
    – Danica
    Commented Feb 12, 2017 at 12:54
  • 1
    $\begingroup$ "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$. $\endgroup$
    – Newton
    Commented Nov 9, 2023 at 18:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.