### Regularisation

Polynomial examples



In [58]:

import numpy as np
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('png', 'pdf')




In [59]:

def fmap_poly(tobs,n):
A = np.zeros((len(tobs),n))
for i, ti in enumerate(tobs):
#A[i,:] = np.array([1,10*ti,-0.5*ti**2])
A[i,:] = np.array([ti**j for j in range(0,n)])
return A




In [60]:

#fine time grid
t = np.linspace(0,1,1000)

#observation times
#tobs = np.array([1,13]) #under-determined
tobs = np.array([0.0,0.1,0.4,0.7,1.0]) #over-determined

order = 10
Ap_obs = fmap_poly(tobs,n=order)
Ap_fine = fmap_poly(t,n=order)




In [61]:

#parameters
thetap_true = np.zeros(order)
thetap_true[1] = 1
thetap_true[2] = 10
thetap_true




Out[61]:

array([  0.,   1.,  10.,   0.,   0.,   0.,   0.,   0.,   0.,   0.])




In [62]:

#fine and obs
plt.plot(tobs,np.dot(Ap_obs,thetap_true),'ro',markersize=10)
#plt.plot(tobs,yp_obs,'ro',markersize=10)

plt.plot(t,np.dot(Ap_fine,thetap_true),'r')
#plt.plot(t,np.dot(Atrue,theta_min),'b')
plt.show()







In [87]:

#observed data
#yp_obs = np.dot(Ap_obs,thetap_true)
yp_obs = np.dot(Ap_obs,thetap_true) + np.random.normal(0,1.0,size=len(tobs)) #with noise

#invert
Ap_obs_pinv = np.linalg.pinv(Ap_obs)
thetap_min = np.dot(Ap_obs_pinv,yp_obs)




In [88]:

print(thetap_min)
print(thetap_true)




[  0.43458213 -24.81296717  80.03423016  -8.92353134 -34.86760772
-28.37125087 -12.94825018   2.29933733  14.76337007  24.20549425]
[  0.   1.  10.   0.   0.   0.   0.   0.   0.   0.]




In [89]:

print(np.linalg.norm(thetap_min,2))
print(np.linalg.norm(thetap_true,2))




100.491117703
10.0498756211




In [90]:

#key figure: simplest, within perfect fit class...
plt.plot(tobs,np.dot(Ap_obs,thetap_true),'ro',markersize=10)
plt.plot(tobs,yp_obs,'bo',markersize=10)

plt.plot(t,np.dot(Ap_fine,thetap_true),'r')
plt.plot(t,np.dot(Ap_fine,thetap_min),'b')
plt.show()






#### Tikhonov



In [91]:

yp_obs = np.dot(Ap_obs,thetap_true) + np.random.normal(0,0.8,size=len(tobs)) #with noise




In [129]:

alpha = 0.4
#alpha = 0.0




In [130]:

b = np.zeros(len(thetap_true))




In [131]:

yp_aug = np.hstack([yp_obs,b])




In [132]:

I = np.eye(len(thetap_true))
#I[0,0] = 0




In [133]:

#augment
Ap_aug = np.vstack([Ap_obs,alpha*I])

#invert
Ap_aug_pinv = np.linalg.pinv(Ap_aug)
thetap_aug_min = np.dot(Ap_aug_pinv,yp_aug)




In [134]:

print(np.linalg.norm(thetap_aug_min[1:],2))
print(np.linalg.norm(thetap_true[1:],2))




3.59510693991
10.0498756211




In [135]:

#key figure: simplest, within perfect fit class...
plt.plot(tobs,np.dot(Ap_obs,thetap_true),'ro',markersize=10)
#plt.plot(tobs,np.dot(Ap_obs,thetap_aug_min),'bo',markersize=10)
plt.plot(tobs,yp_obs,'bo',markersize=10)

plt.plot(t,np.dot(Ap_fine,thetap_true),'r')
#plt.plot(t,np.dot(Ap_fine,thetap_min),'b--')
plt.plot(t,np.dot(Ap_fine,thetap_aug_min),'b',linewidth=2)
plt.show()






#### Loop



In [140]:

alphav = np.logspace(-2.0, 0.5, num=1000)




In [141]:

thetas = np.zeros((len(alphav),len(thetap_true)))
theta_norms = np.zeros(len(alphav))
data_norms = np.zeros(len(alphav))
I = np.eye(len(thetap_true))
#I[0,0] = 0

for i, alphai in enumerate(alphav):

#augment
Ap_aug = np.vstack([Ap_obs,alphai*I])

#invert
Ap_aug_pinv = np.linalg.pinv(Ap_aug)
thetas[i,:] = np.dot(Ap_aug_pinv,yp_aug)

#calc norms
theta_norms[i] = np.linalg.norm(thetas[i,:],2)
data_norms[i] = np.linalg.norm(yp_obs - np.dot(Ap_obs,thetas[i,:]))




In [142]:

plt.plot(theta_norms**2,data_norms**2,linewidth=5)
plt.ylabel(r'$||A\theta-y||^2$',fontsize=14)
plt.xlabel(r'$||\theta||^2$',fontsize=14)
plt.show()







In [143]:

plt.loglog(theta_norms**2,data_norms**2,linewidth=5)
plt.ylabel(r'$||A\theta-y||^2$',fontsize=14)
plt.xlabel(r'$||\theta||^2$',fontsize=14)
plt.show()







In [144]:

plt.plot(np.log10(theta_norms**2),np.log10(data_norms**2),linewidth=5)
plt.ylabel(r'$\log_{10}||A\theta-y||^2$',fontsize=14)
plt.xlabel(r'$\log_{10}||\theta||^2$',fontsize=14)
plt.show()







In [145]:

#plt.plot(alphav,data_norms**2,linewidth=5)
plt.plot(-np.log10(alphav),np.log10(data_norms**2),linewidth=5)
plt.ylabel(r'$||A\theta-y||^2$',fontsize=14)
plt.xlabel(r'$-\log_{10}(\alpha)$',fontsize=14)
#plt.xlabel(r'$\alpha$',fontsize=14)
plt.show()