In [1]:
%matplotlib inline
In [2]:
import numpy as np
import matplotlib.pyplot as plt

Generate $x$ and $y$

In [70]:
n=100
x = np.arange(1,n+1,1)
x = x/n  
y = 2.3 + 20*x + 100*x**2 + 100*x**3 +np.random.randn(n)*10
In [71]:
plt.plot(x,y)
Out[71]:
[<matplotlib.lines.Line2D at 0x128d6c2b0>]

Truncated Vandermonde matrix

In [86]:
A = np.vander(x,increasing=True)[:,:4]
In [87]:
plt.imshow(A,aspect="auto")
Out[87]:
<matplotlib.image.AxesImage at 0x1295d4f60>

Compute the SVD of $A$

In [88]:
U, s, Vs = np.linalg.svd(A, full_matrices=False)
V = Vs.T
U.shape,s.shape,Vs.shape
Out[88]:
((100, 4), (4,), (4, 4))
In [89]:
plt.imshow(U,aspect="auto")
Out[89]:
<matplotlib.image.AxesImage at 0x1296bee80>
In [90]:
plt.imshow(V,aspect="auto")
Out[90]:
<matplotlib.image.AxesImage at 0x1297aacc0>
In [91]:
plt.imshow(U.T.dot(U))
Out[91]:
<matplotlib.image.AxesImage at 0x12989f748>
In [92]:
plt.imshow(U.dot(U.T))
Out[92]:
<matplotlib.image.AxesImage at 0x129986fd0>

Compute $c$

In [93]:
U_y = U.T.dot(y)
U_y_div_sigma = U_y*1./s
c = V.dot(U_y_div_sigma)
In [94]:
plt.plot(c)
Out[94]:
[<matplotlib.lines.Line2D at 0x129a602e8>]
In [95]:
plt.plot(A.dot(c))
Out[95]:
[<matplotlib.lines.Line2D at 0x129abb668>]

Generate new values, compute $A$ and apply $c$

In [96]:
n=200
x_ = np.arange(1,n+1,1)
x_ = x_/n
x_ = x_[np.logical_and(x_>=x.min(), x_<=x.max())]
In [97]:
A_ = np.vander(x_,increasing=True)
A_ = A_[:,:len(c)]
In [98]:
plt.imshow(A_,aspect="auto")
Out[98]:
<matplotlib.image.AxesImage at 0x129b10e80>
In [99]:
plt.plot(x_,A_.dot(c))
plt.plot(x,y,"o")
Out[99]:
[<matplotlib.lines.Line2D at 0x129b4fd30>]
In [ ]:
In [ ]: