3 # Any copyright is dedicated to the Public Domain.
4 # https://creativecommons.org/publicdomain/zero/1.0/
6 # Written by Francois Fleuret <francois@fleuret.org>
9 import matplotlib.pyplot as plt
17 ######################################################################
19 def pol_value(alpha, x):
20 x_pow = x.view(-1, 1) ** torch.arange(alpha.size(0)).view(1, -1)
23 def fit_alpha(x, y, D, a = 0, b = 1, rho = 1e-12):
24 M = x.view(-1, 1) ** torch.arange(D + 1).view(1, -1)
28 q = torch.arange(2, D + 1, dtype = x.dtype).view(1, -1)
30 beta = x.new_zeros(D + 1, D + 1)
31 beta[2:, 2:] = (q-1) * q * (r-1) * r * (b**(q+r-3) - a**(q+r-3))/(q+r-3)
32 l, U = beta.eig(eigenvectors = True)
33 Q = U @ torch.diag(l[:, 0].pow(0.5))
34 B = torch.cat((B, y.new_zeros(Q.size(0))), 0)
35 M = torch.cat((M, math.sqrt(rho) * Q.t()), 0)
37 return torch.lstsq(B, M).solution[:D+1, 0]
39 ######################################################################
42 return torch.abs(torch.abs(x - 0.4) - 0.2) + x/2 - 0.1
44 ######################################################################
48 mse_train = torch.zeros(nb_runs, D_max + 1)
49 mse_test = torch.zeros(nb_runs, D_max + 1)
51 for k in range(nb_runs):
52 x_train = torch.rand(nb_train_samples, dtype = torch.float64)
53 y_train = phi(x_train)
54 if train_noise_std > 0:
55 y_train = y_train + torch.empty_like(y_train).normal_(0, train_noise_std)
56 x_test = torch.linspace(0, 1, 100, dtype = x_train.dtype)
59 for D in range(D_max + 1):
60 alpha = fit_alpha(x_train, y_train, D)
61 mse_train[k, D] = ((pol_value(alpha, x_train) - y_train)**2).mean()
62 mse_test[k, D] = ((pol_value(alpha, x_test) - y_test)**2).mean()
64 mse_train = mse_train.median(0).values
65 mse_test = mse_test.median(0).values
67 ######################################################################
68 # Plot the MSE vs. degree curves
72 ax = fig.add_subplot(1, 1, 1)
75 ax.set_xlabel('Polynomial degree', labelpad = 10)
76 ax.set_ylabel('MSE', labelpad = 10)
78 ax.axvline(x = nb_train_samples - 1, color = 'gray', linewidth = 0.5)
79 ax.plot(torch.arange(D_max + 1), mse_train, color = 'blue', label = 'Train error')
80 ax.plot(torch.arange(D_max + 1), mse_test, color = 'red', label = 'Test error')
82 ax.legend(frameon = False)
84 fig.savefig('dd-mse.pdf', bbox_inches='tight')
86 ######################################################################
87 # Plot some examples of train / test
89 torch.manual_seed(9) # I picked that for pretty
91 x_train = torch.rand(nb_train_samples, dtype = torch.float64)
92 y_train = phi(x_train)
93 if train_noise_std > 0:
94 y_train = y_train + torch.empty_like(y_train).normal_(0, train_noise_std)
95 x_test = torch.linspace(0, 1, 100, dtype = x_train.dtype)
98 for D in range(D_max + 1):
101 ax = fig.add_subplot(1, 1, 1)
102 ax.set_title(f'Degree {D}')
103 ax.set_ylim(-0.1, 1.1)
104 ax.plot(x_test, y_test, color = 'black', label = 'Test values')
105 ax.scatter(x_train, y_train, color = 'blue', label = 'Train samples')
107 alpha = fit_alpha(x_train, y_train, D)
108 ax.plot(x_test, pol_value(alpha, x_test), color = 'red', label = 'Fitted polynomial')
110 ax.legend(frameon = False)
112 fig.savefig(f'dd-example-{D:02d}.pdf', bbox_inches='tight')
114 ######################################################################