ddpol.py

   1 #!/usr/bin/env python
   2
   3 # Any copyright is dedicated to the Public Domain.
   4 # https://creativecommons.org/publicdomain/zero/1.0/
   5
   6 # Written by Francois Fleuret <francois@fleuret.org>
   7
   8 import math
   9 import matplotlib.pyplot as plt
  10 import torch
  11
  12 ######################################################################
  13
  14 def compute_alpha(x, y, D, a = 0, b = 1, rho = 1e-11):
  15     M = x.view(-1, 1) ** torch.arange(D + 1).view(1, -1)
  16     B = y
  17
  18     if D+1 > 2:
  19         q = torch.arange(2, D + 1).view( 1, -1).to(x.dtype)
  20         r = q.view(-1,  1)
  21         beta = x.new_zeros(D + 1, D + 1)
  22         beta[2:, 2:] = (q-1) * q * (r-1) * r * (b**(q+r-3) - a**(q+r-3))/(q+r-3)
  23         l, U = beta.eig(eigenvectors = True)
  24         Q = U @ torch.diag(l[:, 0].pow(0.5))
  25         B = torch.cat((B, y.new_zeros(Q.size(0))), 0)
  26         M = torch.cat((M, math.sqrt(rho) * Q.t()), 0)
  27
  28     alpha = torch.lstsq(B, M).solution.view(-1)[:D+1]
  29
  30     return alpha
  31
  32 ######################################################################
  33
  34 def phi(x):
  35     return 4 * (x - 0.5) ** 2 * (x >= 0.5)
  36
  37 ######################################################################
  38
  39 torch.manual_seed(0)
  40
  41 nb_train_samples = 7
  42 D_max = 16
  43 nb_runs = 250
  44
  45 mse_train = torch.zeros(nb_runs, D_max + 1)
  46 mse_test = torch.zeros(nb_runs, D_max + 1)
  47
  48 for k in range(nb_runs):
  49     x_train = torch.rand(nb_train_samples, dtype = torch.float64)
  50     y_train = phi(x_train)
  51     y_train = y_train + torch.empty(y_train.size(), dtype = y_train.dtype).normal_(0, 0.1)
  52     x_test = torch.linspace(0, 1, 100, dtype = x_train.dtype)
  53     y_test = phi(x_test)
  54
  55     for D in range(D_max + 1):
  56         alpha = compute_alpha(x_train, y_train, D)
  57         X_train = x_train.view(-1, 1) ** torch.arange(D + 1).view(1, -1)
  58         X_test = x_test.view(-1, 1) ** torch.arange(D + 1).view(1, -1)
  59         mse_train[k, D] = ((X_train @ alpha - y_train)**2).mean()
  60         mse_test[k, D] = ((X_test @ alpha - y_test)**2).mean()
  61
  62 mse_train = mse_train.median(0).values
  63 mse_test = mse_test.median(0).values
  64
  65 ######################################################################
  66
  67 torch.manual_seed(4) # I picked that for pretty
  68
  69 x_train = torch.rand(nb_train_samples, dtype = torch.float64)
  70 y_train = phi(x_train)
  71 y_train = y_train + torch.empty(y_train.size(), dtype = y_train.dtype).normal_(0, 0.1)
  72 x_test = torch.linspace(0, 1, 100, dtype = x_train.dtype)
  73 y_test = phi(x_test)
  74
  75 for D in range(D_max + 1):
  76     fig = plt.figure()
  77
  78     ax = fig.add_subplot(1, 1, 1)
  79     ax.set_title(f'Degree {D}')
  80     ax.set_ylim(-0.1, 1.1)
  81     ax.plot(x_test, y_test, color = 'blue', label = 'Test values')
  82     ax.scatter(x_train, y_train, color = 'blue', label = 'Training examples')
  83
  84     alpha = compute_alpha(x_train, y_train, D)
  85     X_test = x_test.view(-1, 1) ** torch.arange(D + 1).view(1, -1)
  86     ax.plot(x_test, X_test @ alpha, color = 'red', label = 'Fitted polynomial')
  87
  88     ax.legend(frameon = False)
  89
  90     fig.savefig(f'dd-example-{D:02d}.pdf', bbox_inches='tight')
  91
  92 ######################################################################
  93
  94 fig = plt.figure()
  95
  96 ax = fig.add_subplot(1, 1, 1)
  97 ax.set_yscale('log')
  98 ax.set_xlabel('Polynomial degree', labelpad = 10)
  99 ax.set_ylabel('MSE', labelpad = 10)
 100
 101 ax.axvline(x = nb_train_samples - 1, color = 'gray', linewidth = 0.5)
 102 ax.plot(torch.arange(D_max + 1), mse_train, color = 'blue', label = 'Train error')
 103 ax.plot(torch.arange(D_max + 1), mse_test, color = 'red', label = 'Test error')
 104
 105 ax.legend(frameon = False)
 106
 107 fig.savefig('dd-mse.pdf', bbox_inches='tight')
 108
 109 ######################################################################