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.
