{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"toc": true
},
"source": [
"Table of Contents
\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Why is gradient descent so bad at optimizing polynomial regression?# \n",
"\n",
"Question from Stackexchange:\n",
"https://stats.stackexchange.com/questions/350130/why-is-gradient-descent-so-bad-at-optimizing-polynomial-regression\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Linear regression\n",
"\n",
"#### Cost function\n",
"$J(\\theta) = \\frac{1}{2m}\\sum_{i=1}^m(h_\\theta(x^{(i)}) - y^{(i)})^2 $\n",
"\n",
"$J(\\theta) = \\frac{1}{2m}(X\\theta - y)^T(X\\theta - y) $ (vectorized version)\n",
"\n",
"#### Gradient \n",
"$\\frac{\\partial J(\\theta)}{\\partial \\theta} = \\frac{1}{m}X^T(X\\theta - y) $ \n",
"\n",
"##### Hessian\n",
"$\\frac{\\partial^2 J(\\theta)}{\\partial \\theta^2} = \\frac{1}{m}X^T X $ \n",
"\n",
"## Polynomial regression\n",
"The design matrix is of the form: \n",
"\n",
"$ \\mathbf{X = [1 , x , x^2 , x^3 , ... , x}^n]$\n",
"\n",
"### Libraries"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/anaconda3/lib/python3.6/site-packages/statsmodels/compat/pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.\n",
" from pandas.core import datetools\n"
]
}
],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import scipy.optimize as opt\n",
"from sklearn import linear_model\n",
"import statsmodels.api as sm\n",
"\n",
"from matplotlib import pyplot as plt\n",
"from mpl_toolkits.mplot3d import Axes3D\n",
"%matplotlib inline\n",
"plt.style.use('seaborn-white')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Helper Functions "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def costfunction(theta,X,y):\n",
" m = np.size(y)\n",
" theta = theta.reshape(-1,1)\n",
" \n",
" #Cost function in vectorized form\n",
" h = X @ theta\n",
" J = float((1./(2*m)) * (h - y).T @ (h - y)); \n",
" return J;\n",
"\n",
"def gradient_descent(theta,X,y,alpha = 0.0005,num_iters=1000):\n",
" m = np.size(y)\n",
" J_history = np.empty(num_iters)\n",
" count_history = np.empty(num_iters)\n",
" theta_1_hist, theta_2_hist = [], []\n",
" \n",
" for i in range(num_iters):\n",
" #Grad function in vectorized form\n",
" h = X @ theta\n",
" theta = theta - alpha * (1/m)* (X.T @ (h-y))\n",
" \n",
" #Tracker values for plotting\n",
" J_history[i] = costfunction(theta,X,y)\n",
" count_history[i] = i\n",
" theta_1_hist.append(theta[0,0])\n",
" theta_2_hist.append(theta[1,0])\n",
" \n",
" return theta, J_history,count_history, theta_1_hist, theta_2_hist\n",
"\n",
"def grad(theta,X,y):\n",
" #Initializations\n",
" theta = theta[:,np.newaxis]\n",
" m = len(y)\n",
" grad = np.zeros(theta.shape)\n",
" \n",
" #Computations\n",
" h = X @ theta\n",
" \n",
" grad = (1./m)*(X.T @ ( h - y))\n",
" \n",
" return (grad.flatten())\n",
"\n",
"\n",
"def polynomial_features(data, deg):\n",
" data_copy=data.copy()\n",
" \n",
" for i in range(1,deg):\n",
" data_copy['X'+str(i+1)]=data_copy['X'+str(i)]*data_copy['X1']\n",
"\n",
" return data_copy\n",
"\n",
"def hessian(theta,X,y):\n",
" m,n = X.shape\n",
" X = X.values\n",
" return ((1./m)*(X.T @ X))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# : Comparing results for high order polynomial regression "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Initializing the data"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"#Create data from sin function with uniform noise\n",
"x = np.linspace(0.1,1,40)\n",
"noise = np.random.uniform( size = 40)\n",
"y = np.sin(x * 1.5 * np.pi ) \n",
"y_noise = (y + noise).reshape(-1,1)\n",
"y_noise = y_noise - y_noise.mean() #Centering the data\n",
"\n",
"degree = 7\n",
"X_d = polynomial_features(pd.DataFrame({'X0':1,'X1': x}),degree)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Closed form solution "
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[-2.95437720e+00],\n",
" [ 6.50561451e+01],\n",
" [-5.29277649e+02],\n",
" [ 2.29876198e+03],\n",
" [-5.44238666e+03],\n",
" [ 6.98373410e+03],\n",
" [-4.56867395e+03],\n",
" [ 1.19415965e+03]])"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def closed_form_solution(X,y):\n",
" return np.linalg.inv(X.T @ X) @ X.T @ y\n",
"\n",
"coefs = closed_form_solution(X_d.values,y_noise)\n",
"coefs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Numpy only fit"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 0.46548572],\n",
" [ 1.24606966],\n",
" [-1.36884742],\n",
" [-1.79417815],\n",
" [-1.19320446],\n",
" [-0.36286233],\n",
" [ 0.39300774],\n",
" [ 0.99476541]])"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"stepsize = .1\n",
"theta_result_1,J_history_1, count_history_1, theta_1_hist, theta_2_hist = gradient_descent(np.zeros((len(X_d.T),1)).reshape(-1,1), X_d,y_noise,alpha = stepsize,num_iters=5000)\n",
"display(theta_result_1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Sciy optimize fit using first order derivative only\n",
"#### Comment: BFGS does very well but requires adjustment of options\n",
"In particular tolerance must be made smaller as the cost function is very flat near the global minimum"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.036158\n",
" Iterations: 152\n",
" Function evaluations: 156\n",
" Gradient evaluations: 156\n"
]
},
{
"data": {
"text/plain": [
"array([-2.95437890e+00, 6.50561760e+01, -5.29277889e+02, 2.29876294e+03,\n",
" -5.44238875e+03, 6.98373664e+03, -4.56867556e+03, 1.19416006e+03])"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import scipy.optimize as opt\n",
"\n",
"theta_init = np.ones((len(X_d.T),1)).reshape(-1,1)\n",
"model_t = opt.minimize(fun = costfunction, x0 = theta_init , args = (X_d, y_noise), \n",
" method = 'BFGS', jac = grad, options={'maxiter':1000, 'gtol': 1e-10, 'disp' : True})\n",
"model_t.x"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Sciy optimize fit using hessian matrix\n",
"#### As expected, 2nd order information allows to converge much faster"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.036158\n",
" Iterations: 20\n",
" Function evaluations: 21\n",
" Gradient evaluations: 21\n",
" Hessian evaluations: 20\n"
]
},
{
"data": {
"text/plain": [
"array([-2.95437904e+00, 6.50561791e+01, -5.29277914e+02, 2.29876304e+03,\n",
" -5.44238897e+03, 6.98373692e+03, -4.56867574e+03, 1.19416011e+03])"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import scipy.optimize as opt\n",
"\n",
"theta_init = np.ones((len(X_d.T),1)).reshape(-1,1)\n",
"model_t = opt.minimize(fun = costfunction, x0 = theta_init , args = (X_d, y_noise), \n",
" method = 'dogleg', jac = grad, hess= hessian, options={'maxiter':1000, 'disp':True})\n",
"model_t.x\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Sklearn fit "
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[-2.95437904e+00, 6.50561790e+01, -5.29277914e+02,\n",
" 2.29876304e+03, -5.44238897e+03, 6.98373691e+03,\n",
" -4.56867573e+03, 1.19416011e+03]])"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from sklearn import linear_model\n",
"model_d = linear_model.LinearRegression(fit_intercept=False)\n",
"model_d.fit(X_d,y_noise)\n",
"model_d.coef_"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Statsmodel fit "
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: y R-squared: 0.891\n",
"Model: OLS Adj. R-squared: 0.868\n",
"Method: Least Squares F-statistic: 37.49\n",
"Date: Thu, 07 Jun 2018 Prob (F-statistic): 1.16e-13\n",
"Time: 10:21:12 Log-Likelihood: -4.2236\n",
"No. Observations: 40 AIC: 24.45\n",
"Df Residuals: 32 BIC: 37.96\n",
"Df Model: 7 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"X0 -2.9544 4.220 -0.700 0.489 -11.551 5.642\n",
"X1 65.0562 86.252 0.754 0.456 -110.633 240.745\n",
"X2 -529.2779 673.778 -0.786 0.438 -1901.719 843.163\n",
"X3 2298.7630 2651.280 0.867 0.392 -3101.718 7699.244\n",
"X4 -5442.3890 5762.460 -0.944 0.352 -1.72e+04 6295.358\n",
"X5 6983.7369 7007.140 0.997 0.326 -7289.340 2.13e+04\n",
"X6 -4568.6757 4460.607 -1.024 0.313 -1.37e+04 4517.283\n",
"X7 1194.1601 1156.577 1.032 0.310 -1161.709 3550.029\n",
"==============================================================================\n",
"Omnibus: 1.721 Durbin-Watson: 2.129\n",
"Prob(Omnibus): 0.423 Jarque-Bera (JB): 1.447\n",
"Skew: 0.306 Prob(JB): 0.485\n",
"Kurtosis: 2.297 Cond. No. 3.02e+05\n",
"==============================================================================\n",
"\n",
"Warnings:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"[2] The condition number is large, 3.02e+05. This might indicate that there are\n",
"strong multicollinearity or other numerical problems.\n"
]
}
],
"source": [
"import statsmodels.api as sm\n",
"model_sm = sm.OLS(y_noise, X_d)\n",
"res = model_sm.fit()\n",
"print(res.summary())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# : Repeating the experiment with 2 polynomial variables and visualizing the results\n",
"Here we will focus on a 2-D design matrix with $x$ and $x^2$ values. The y values have been centered so we will ignore the constant term and y-intercept"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Initializing the data"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"#Create data from sin function with uniform noise\n",
"x = np.linspace(0.1,1,40) #Adjusting the starting point to reduce numerical instability\n",
"noise = np.random.uniform( size = 40)\n",
"y = np.sin(x * 1.5 * np.pi ) \n",
"y_noise = (y + noise).reshape(-1,1)\n",
"y_noise = y_noise - y_noise.mean() #Centering the data\n",
"\n",
"#2nd order polynomial only\n",
"degree = 2\n",
"X_d = polynomial_features(pd.DataFrame({'X1': x}),degree)\n",
"\n",
"#Starting point for gradient descent - see later diagrams\n",
"initial_theta = np.array([0,-2]).reshape(-1,1)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" X1 | \n",
" X2 | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" 0.100000 | \n",
" 0.010000 | \n",
"
\n",
" \n",
" | 1 | \n",
" 0.123077 | \n",
" 0.015148 | \n",
"
\n",
" \n",
" | 2 | \n",
" 0.146154 | \n",
" 0.021361 | \n",
"
\n",
" \n",
" | 3 | \n",
" 0.169231 | \n",
" 0.028639 | \n",
"
\n",
" \n",
" | 4 | \n",
" 0.192308 | \n",
" 0.036982 | \n",
"
\n",
" \n",
"
\n",
"