From: Francois Fleuret Date: Fri, 5 Jun 2020 22:58:01 +0000 (+0200) Subject: Initial commit. X-Git-Url: https://fleuret.org/cgi-bin/gitweb/gitweb.cgi?a=commitdiff_plain;h=606f9f7294872c8291ea6449d3b88f474b21adc9;p=pytorch.git Initial commit. --- diff --git a/ddpol.py b/ddpol.py new file mode 100755 index 0000000..97d2ff5 --- /dev/null +++ b/ddpol.py @@ -0,0 +1,109 @@ +#!/usr/bin/env python + +# Any copyright is dedicated to the Public Domain. +# https://creativecommons.org/publicdomain/zero/1.0/ + +# Written by Francois Fleuret + +import math +import matplotlib.pyplot as plt +import torch + +###################################################################### + +def compute_alpha(x, y, D, a = 0, b = 1, rho = 1e-11): + M = x.view(-1, 1) ** torch.arange(D + 1).view(1, -1) + B = y + + if D+1 > 2: + q = torch.arange(2, D + 1).view( 1, -1).to(x.dtype) + r = q.view(-1, 1) + beta = x.new_zeros(D + 1, D + 1) + beta[2:, 2:] = (q-1) * q * (r-1) * r * (b**(q+r-3) - a**(q+r-3))/(q+r-3) + l, U = beta.eig(eigenvectors = True) + Q = U @ torch.diag(l[:, 0].pow(0.5)) + B = torch.cat((B, y.new_zeros(Q.size(0))), 0) + M = torch.cat((M, math.sqrt(rho) * Q.t()), 0) + + alpha = torch.lstsq(B, M).solution.view(-1)[:D+1] + + return alpha + +###################################################################### + +def phi(x): + return 4 * (x - 0.5) ** 2 * (x >= 0.5) + +###################################################################### + +torch.manual_seed(0) + +nb_train_samples = 7 +D_max = 16 +nb_runs = 250 + +mse_train = torch.zeros(nb_runs, D_max + 1) +mse_test = torch.zeros(nb_runs, D_max + 1) + +for k in range(nb_runs): + x_train = torch.rand(nb_train_samples, dtype = torch.float64) + y_train = phi(x_train) + y_train = y_train + torch.empty(y_train.size(), dtype = y_train.dtype).normal_(0, 0.1) + x_test = torch.linspace(0, 1, 100, dtype = x_train.dtype) + y_test = phi(x_test) + + for D in range(D_max + 1): + alpha = compute_alpha(x_train, y_train, D) + X_train = x_train.view(-1, 1) ** torch.arange(D + 1).view(1, -1) + X_test = x_test.view(-1, 1) ** torch.arange(D + 1).view(1, -1) + mse_train[k, D] = ((X_train @ alpha - y_train)**2).mean() + mse_test[k, D] = ((X_test @ alpha - y_test)**2).mean() + +mse_train = mse_train.median(0).values +mse_test = mse_test.median(0).values + +###################################################################### + +torch.manual_seed(4) # I picked that for pretty + +x_train = torch.rand(nb_train_samples, dtype = torch.float64) +y_train = phi(x_train) +y_train = y_train + torch.empty(y_train.size(), dtype = y_train.dtype).normal_(0, 0.1) +x_test = torch.linspace(0, 1, 100, dtype = x_train.dtype) +y_test = phi(x_test) + +for D in range(D_max + 1): + fig = plt.figure() + + ax = fig.add_subplot(1, 1, 1) + ax.set_title(f'Degree {D}') + ax.set_ylim(-0.1, 1.1) + ax.plot(x_test, y_test, color = 'blue', label = 'Test values') + ax.scatter(x_train, y_train, color = 'blue', label = 'Training examples') + + alpha = compute_alpha(x_train, y_train, D) + X_test = x_test.view(-1, 1) ** torch.arange(D + 1).view(1, -1) + ax.plot(x_test, X_test @ alpha, color = 'red', label = 'Fitted polynomial') + + ax.legend(frameon = False) + + fig.savefig(f'dd-example-{D:02d}.pdf', bbox_inches='tight') + +###################################################################### + +fig = plt.figure() + +ax = fig.add_subplot(1, 1, 1) +ax.set_yscale('log') +ax.set_xlabel('Polynomial degree', labelpad = 10) +ax.set_ylabel('MSE', labelpad = 10) + +ax.axvline(x = nb_train_samples - 1, color = 'gray', linewidth = 0.5) +ax.plot(torch.arange(D_max + 1), mse_train, color = 'blue', label = 'Train error') +ax.plot(torch.arange(D_max + 1), mse_test, color = 'red', label = 'Test error') + +ax.legend(frameon = False) + +fig.savefig('dd-mse.pdf', bbox_inches='tight') + +######################################################################