TODAY: Classical Variants of Linear Regression¶
- Recursive/Sequential least squares
- Total least squares
- Nonlinear least squares
Example problem¶
1D regression over $[-5, 5]$,
$$ f(x) = (1-x-2x^2)\sin(x) $$In [3]:
plt.plot(x,y,'b-')
plt.xlabel('x')
plt.ylabel('y')
Out[3]:
Text(0, 0.5, 'y')
Basic least squares¶
$$ f(x) \approx \sum_{i=1}^m w_i\phi_i(x) \equiv w^T\phi(x) $$Choose $$ \phi(x) = [1, x, x^2, \cdots, x^{m-1}] $$
In [5]:
Fp = lambda x: fphi(x, 8)
ny = 2*np.random.randn(len(y),)
w = leastSquares(x, y+ny, Fp)
yhat = Fp(x).dot(w)
plt.plot(x, y+ny, 'y.', markerfacecolor='none', label='Data')
plt.plot(x, y, 'b-', label='Truth')
plt.plot(x, yhat, 'r--', label=f'Prediction, Error={error(y, yhat)*100:3.1f}%')
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
Out[5]:
Text(0, 0.5, 'y')
Recursive/Sequential least squares¶
Still choose same polynomials.
In [7]:
Fp = lambda x: fphi(x, 8)
idx = np.random.permutation(len(x))
# ny = 0*np.random.randn(len(y),)
X, Y = x[idx], (y+ny)[idx]
ws, es = seqLeastSquares(X, Y, Fp, 20)
plt.semilogy(np.arange(len(es))+1, es)
plt.xlabel('Iteration')
plt.ylabel('Global Error')
print(es[-1])
0.11374668925285274
In [8]:
F = Fp(x)
plt.plot(x, y+ny, 'y.', markerfacecolor='none')
plt.plot(x, y, 'b-')
for _i in range(5):
yhat = F.dot(ws[_i])
plt.plot(x, yhat, 'r-', alpha=0.1*(_i+1))
yhat = F.dot(ws[len(es)-_i-1])
plt.plot(x, yhat, 'k-', alpha=0.1*(5-_i))
plt.xlabel('x')
plt.ylabel('y')
plt.ylim([-45,55])
Out[8]:
(-45.0, 55.0)
Total least squares¶
Still choose same polynomials.
In [10]:
Fp = lambda x: fphi(x, 8)
nx = 0.2*np.random.randn(len(y),)
ny = 2*np.random.randn(len(y),)
wr = leastSquares(x+nx, y+ny, Fp)
wt, Ft = totLeastSquares(x+nx, y+ny, Fp)
yr = Fp(x).dot(wr)
yt = (Fp(x+nx)+Ft).dot(wt)
# yt = (Fp(x)).dot(wt)
plt.plot(x+nx, y+ny, 'y.', markerfacecolor='none', label='Data')
plt.plot(x, y, 'b-', label='Truth')
plt.plot(x, yr, 'r--', label=f'Basic LR, Error={error(y, yr)*100:3.1f}%')
plt.plot(x, yt, 'g-.', label=f'Total LR, Error={error(y, yt)*100:3.1f}%')
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
Out[10]:
Text(0, 0.5, 'y')
Nonlinear least squares¶
- Fixed learning rate
- Newton's method
- Gauss-Newton
- Levenberg-Marquardt
Optimize model $$ f(x) = \cos(w_1 x)(1+w_2 x^2) $$
Gradient $$ g(x) = \nabla f = \begin{bmatrix} -x\sin(w_1 x)(1+w_2 x^2) \\ \cos(w_1 x)x^2 \end{bmatrix} $$
Hessian $$ h=\nabla g = \begin{bmatrix} -x^2\cos(w_1 x)(1+w_2 x^2) & -x^3\sin(w_1 x) \\ \text{sym} & 0 \end{bmatrix} $$
In [13]:
def _func(x, w):
return np.cos(w[0]*x) * (1+w[1]*x**2)
def _grad(x, w):
return np.vstack([
-x*np.sin(w[0]*x) * (1+w[1]*x**2),
np.cos(w[0]*x) * x**2
])
def _hess(x, w):
_h = np.zeros((2,2,len(x)))
_h[0,0] = -x**2*np.cos(w[0]*x) * (1+w[1]*x**2)
_h[0,1] = -x**3 * np.sin(w[0]*x)
_h[1,0] = _h[0,1]
return _h
w = [1.8, 0.5]
x = np.linspace(-1,1,21)
y = _func(x, w)
DAT = np.vstack([x,y]).T
plt.plot(x, y)
Out[13]:
[<matplotlib.lines.Line2D at 0x11919f460>]
In [28]:
# x0 = [3.0,0.5]
x0 = [1.5, -1.5]
# xs, fs, gs, hs = gradDesc(FGH, x0, fixLR(0.01))
# xs, fs, gs, hs = gradDesc(FGH, x0, LM(15.))
# xs, fs, gs, hs = gradDesc(FGH, x0, NR, ifLS=True)
xs, fs, gs, hs = gradDesc(FGHP, x0, LM(10.))
# xs, fs, gs, hs = gradDesc(FGHP, x0, GN)
f1, ax = pltFunc()
ax[0] = pltTraj(ax[0], xs)
ax[1] = pltHist(ax[1], xs)
In [ ]: