# Kernel learning example¶

Here, we imagine that we want to use linear regression with a Gaussian kernel to fit a noisy sine wave.

In [1]:
using PyPlot
using Random
using LinearAlgebra
import Statistics.mean

â”Œ Warning: PyPlot is using tkagg backend, which is known to cause crashes on MacOS (#410); use the MPLBACKEND environment variable to request a different backend.
â”” @ PyPlot /Users/cdesa/.julia/packages/PyPlot/XHEG0/src/init.jl:192

In [2]:
Random.seed!(123456);
N = 10000; # size of data set
x = 6 * rand(N);
y = sin.(x.^2) .+ 0.2 * randn(N);

In [3]:
scatter(x,y; s=1);
plot(collect(0:0.01:6), sin.(collect(0:0.01:6).^2), "--k"; linewidth=1);
xlabel("x");
ylabel("y");
title("Training Set");

In [4]:
gamma = 100;
function K(x1, x2)
return exp(-gamma*(x1 - x2)^2/2);
end

Out[4]:
K (generic function with 1 method)
In [5]:
# compute the Gram matrix
@time begin
G = zeros(N, N);
for i = 1:N
for j = 1:N
G[i,j] = K(x[i], x[j]);
end
end
end

 25.088039 seconds (889.85 M allocations: 15.497 GiB, 9.33% gc time, 0.13% compilation time)

In [6]:
# actually solve the regularized linear regression problem
@time w = (G + 0.001 * LinearAlgebra.I) \ y;

  6.881990 seconds (2.47 M allocations: 1.624 GiB, 1.11% gc time, 10.07% compilation time)

In [7]:
function predict(xt :: Float64)
rv = 0.0;
for i = 1:N
rv += w[i] * K(x[i], xt);
end
return rv;
end

Out[7]:
predict (generic function with 1 method)
In [8]:
scatter(x,y; s=1);
x_test = collect(0:0.01:6);
y_test = sin.(x_test.^2);
y_pred = predict.(x_test);
plot(x_test, y_pred, "k"; linewidth=1);
xlabel("x");
ylabel("y");
title("Learned model: Gram matrix");

println("average test error: $(sqrt(mean((y_test - y_pred).^2)))")  average test error: 0.02017090445761806  ## Now, what if we do it with random Fourier features?¶ In [9]: D = 1000; U = randn(D) * sqrt(gamma); b = 2 * pi * rand(D);  In [10]: # compute the feature map @time begin phi = zeros(N, D); for i = 1:N for j = 1:D phi[i,j] = sqrt(2)*cos(U[j] * x[i] + b[j]); end end end   2.994333 seconds (99.87 M allocations: 1.714 GiB, 13.57% gc time, 0.92% compilation time)  In [11]: # actually solve the (now not regularized) linear regression problem @time w_rff = (phi' * phi + 0.001 * LinearAlgebra.I) \ (phi' * y);   0.901588 seconds (2.66 M allocations: 166.771 MiB, 12.40% gc time, 87.94% compilation time)  In [12]: function predict_rff(xt :: Float64) phi_xt = sqrt(2)*cos.(U .* xt .+ b); return dot(phi_xt, w_rff); end  Out[12]: predict_rff (generic function with 1 method) In [13]: scatter(x,y; s=1); x_test = collect(0:0.01:6); y_test = sin.(x_test.^2); y_pred = predict_rff.(x_test); plot(x_test, y_pred, "k"; linewidth=1); xlabel("x"); ylabel("y"); title("Learned model: Random Fourier features"); println("average test error:$(sqrt(mean((y_test - y_pred).^2)))")

average test error: 0.016652825336847945

In [ ]:


In [ ]: