I am trying to understand either intuitively/geometricaly and/or mathematicaly why the followings are equivalent:

  • Classic Ordinay Least Squares linear regression
  • Linear-kernelized Ordinary Least Squares linear regression

A python example shows that they are numericaly equivalent: I use an OLS regression using only numpy and matrix inversion, and sklearn with classic linear regressor and KernelRidge regressor with 0 regularization. This example show that although the kernel approach uses 100 coefficients, the solution is linear.

%matplotlib qt
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.kernel_ridge import KernelRidge
from sklearn.metrics.pairwise import pairwise_kernels

n = 100
X_ = np.random.uniform(3, 10, n).reshape(-1, 1)
beta_0 = 2
beta_1 = 2
true_y = beta_1 * X_ + beta_0
noise = np.random.randn(n, 1) * 0.5 # change the scale to reduce/increase noise
y = true_y + noise

fig, axes = plt.subplots(1, 2, squeeze=False, sharex=True, sharey=True, figsize=(18, 8))
axes[0, 0].plot(X_, y, "o", label="Input data")
axes[0, 0].plot(X_, true_y, '--', label='True linear relation')
axes[0, 0].set_xlim(0, 11)
axes[0, 0].set_ylim(0, 30)
axes[0, 0].legend()

# f_0 is a column of 1s
# f_1 is the column of x1
X = np.c_[np.ones((n, 1)), X_]

beta_OLS_scratch = np.linalg.inv(X.T @ X) @ X.T @ y
lr = LinearRegression(
    fit_intercept=False, # do not fit intercept independantly, since we added the 1 column for this purpose
).fit(X, y)
krr = KernelRidge(
    alpha=0, # no regularization
).fit(X, y)

new_X = np.linspace(0, 15, 50).reshape(-1, 1)
new_X = np.c_[np.ones((50, 1)), new_X]
new_y_OLS_scratch = new_X @ beta_OLS_scratch 
new_y_lr = lr.predict(new_X)
new_y_krr = krr.predict(new_X)
axes[0, 1].plot(X_, y, 'o', label='Input data')
axes[0, 1].plot(new_X[:, 1], new_y_OLS_scratch, '-o', alpha=0.5,  label=r"OLS scratch solution")# $\beta=[" + f'{beta_OLS_scratch.flatten()}'+r"]$")
axes[0, 1].plot(new_X[:, 1], new_y_lr, '-*', alpha=0.5, label=r"sklearn.lr OLS solution")# $\beta=["+f"{lr.coef_.flatten()}" + r"]$")
axes[0, 1].plot(new_X[:, 1], new_y_krr, '-*', alpha=0.5, label=r"sklearn.krr linear kernel")# $\beta=["+f"{lr.coef_.flatten()}" + r"]$")
axes[0, 1].legend()

All regressions lead to the same results

Why it make sense that they are equivalent: Mathematicaly, I feel like since the kernel uses only linear operations, it kinda make sense that we get a linear solution, but this feels not enough justification compared to the fact that it creates a huge dimension space. Also, it kinda make sence since the associated transformation to linear kernel is the identity: $$K(x,y) = \left< \phi(x),\phi(y)\right> = \phi(x)^T \phi(y) = x^T y$$ since for linear kernel $\phi(x)=x$.

But also why it doesn't make sense: Intuitively, I would've expected that since using a kernel expands the feature space to lots of dimensions (given by the number of inputs in the traning set, here $m'=n=100$ dimensions), the solution (projected back onto the original space) would've been non-linear (and probably expose overfitting). But no, the final solution is perfectly linear AND identical to classic OLS regression.

My take at prooving this mathematicaly


  • $X_{ij}$ row-i, column-j element of input matrix $X$ of shape n by m.
  • $i=1...n$, with $n$ the number of rows in the input data matrix
  • $j=1...m$, with $m$ the number of columns in the input data matrix
  • $j^{'}=1...m'=n$, with $m'$ the number of columns of the linear kernel approach, that is $n$
  • $\beta_j$ the regression coefficient of the standard OLS regression, $j=1...m$
  • $\beta^{'}_{j^{'}}$ the regression coefficient of the kernel approach, $j^{'}=1...m^{'}=n$

The initial linear regression is: $$Y_n = \begin{bmatrix} y_1\\ ...\\ y_n \end{bmatrix}_{[n]} = X_{[n\times m]} \begin{bmatrix} \beta_0\\ ...\\ \beta_m \end{bmatrix}_{[m]} + \epsilon_{[n]}$$ $$y_i = \sum_j^m X_{ij}\beta_j + \epsilon_i$$ with solution: $$\beta = (X^TX)^{-1}X^T y$$

The linear-kernel transformed problem becomes with kernel: $$K = k_{ij} = X_i \cdot X_j = \sum_k^m X_{ik} X_{jk}$$ we get $$Y_n = \begin{bmatrix} y_1\\ ...\\ y_n\\ \end{bmatrix}= K_{[n\times m^{'}=n]} \begin{bmatrix} \beta_0\\ ...\\ \beta_{m^{'}}\\ \end{bmatrix}_{[m^{'}=n]} + \epsilon^{'}_{[n]} $$ $$y_i = \sum_{j^{'}}^{m^{'}=n} K_{ij^{‘}} \beta^{'}_{j^{'}} + \epsilon^{'}_{i}$$ $$y_i = \sum_{j^{'}}^{m^{'}=n} \left( \sum_k^m X_{ik}X_{j'k} \right) \beta^{'}_{j^{'}} + \epsilon^{'}_{i}$$ $$y_i = \sum_{j^{'}}^{m^{'}=n} \overrightarrow{X_i} \cdot \overrightarrow{X_{j'}} \beta^{'}_{j^{'}} + \epsilon^{'}_{i}$$ with solution: $$\beta^{'} = (K^TK)^{—1}K^Ty$$ Since $$K = (\overrightarrow{X_i}\cdot \overrightarrow{X_j})_{ij} = XX^T$$ $$\beta^{'} = ((XX^T)^TXX^T)^{—1}(XX^T)^Ty \text{(1)}$$

So I guess to show that both approaches are equivalent, I would have to proove that any new dataset $X_{\text{new}}$ would give the new y values, that is: $$Y_{\text{new}}= X_{\text{new}} \beta= X_{\text{new}} X^T \beta^{'} $$ Am I right ? is this possible ?

Conclusion Although some arguments can be accounted for the fact that linear regression and linear-kernel regression lead to the same solution, I cannot understand why using a high-dimension space still gives back a linear solution once projected back onto the low-dimension space. Also, I cannot proove that both approaches are equivalent mathematicaly, although it seems possible.

[EDIT] With this post, I ended up showing exactly the mathematical proof I was looking for:
Continuing from $\text{(1)}$ and using $(XX^T)^T=X^{TT}X^T=XX^T$: $$\beta^{'} = ((XX^T)^TXX^T)^{—1}(XX^T)^Ty \text{(1)}$$ $$\beta^{'} = (XX^TXX^T)^{—1}XX^Ty$$ $$\beta^{'} = X^{-T}(X^TX)^{-1}X^{—1}XX^Ty$$ $$\beta^{'} = X^{-T}(X^TX)^{-1}X^Ty$$ so $$X_{\text{new}} X^T \beta^{'} = X_{\text{new}} X^T X^{-T}(X^TX)^{-1}X^Ty = X_{\text{new}} (X^TX)^{-1}X^Ty = X_{\text{new}} \beta $$ I am still looking for a geometric/intuitive proof.


1 Answer 1


Well, this might or might not be an intuition, but the simplest explanation that I can come up with is the following.

take a linear model with $w$ as weights, we can prove that $$ w = \sum_i \alpha_i x_i $$

that is, the weights must be a linear combination of the training set.

Then, we can simply rewrite the linear model loss function with this new formalization:
$$ L(\theta) = \sum_j(w^Tx_j- y_j)^2 = \sum_j \left[(\sum_i \alpha_i x_i)x_j - y_i\right]^2 $$

And now you can clearly see that this is literally kernel regression with linear kernel


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.