From 2545aec1eb1465646a72b264176c043a72feb68b Mon Sep 17 00:00:00 2001 From: Morten Hjorth-Jensen Date: Thu, 17 Sep 2020 11:14:54 +0200 Subject: [PATCH] adding material to week 38 --- doc/src/week38/DataFiles/chddata.csv | 1 + doc/src/week38/DataFiles/chddata.csv~ | 101 --- doc/src/week38/week38.do.txt | 908 +++++++++++++++++++++++++- 3 files changed, 906 insertions(+), 104 deletions(-) delete mode 100644 doc/src/week38/DataFiles/chddata.csv~ diff --git a/doc/src/week38/DataFiles/chddata.csv b/doc/src/week38/DataFiles/chddata.csv index 591f7ad32..9c52675fd 100644 --- a/doc/src/week38/DataFiles/chddata.csv +++ b/doc/src/week38/DataFiles/chddata.csv @@ -1,3 +1,4 @@ +ID Age Agegroup CHD 1 21 1 0 2 23 1 0 3 25 1 1 diff --git a/doc/src/week38/DataFiles/chddata.csv~ b/doc/src/week38/DataFiles/chddata.csv~ deleted file mode 100644 index 9c52675fd..000000000 --- a/doc/src/week38/DataFiles/chddata.csv~ +++ /dev/null @@ -1,101 +0,0 @@ -ID Age Agegroup CHD -1 21 1 0 -2 23 1 0 -3 25 1 1 -4 29 1 0 -5 21 1 0 -6 24 1 0 -7 27 1 0 -8 29 1 0 -9 28 1 0 -10 26 1 0 -11 30 2 0 -12 31 2 0 -13 31 2 0 -14 31 2 1 -15 32 2 0 -16 34 2 0 -17 34 2 0 -18 31 2 0 -19 32 2 0 -20 32 2 0 -21 33 2 0 -22 34 2 0 -23 31 2 1 -24 30 2 0 -25 33 2 0 -26 36 3 1 -27 35 3 0 -28 35 3 0 -29 38 3 0 -30 37 3 1 -31 36 3 0 -32 35 3 0 -33 39 3 0 -34 39 3 0 -35 38 3 1 -36 37 3 0 -37 37 3 0 -38 40 4 0 -39 41 4 1 -40 44 4 0 -41 44 4 0 -42 43 4 1 -43 42 4 0 -44 41 4 0 -45 40 4 1 -46 42 4 0 -47 42 4 0 -48 43 4 0 -49 44 4 1 -50 44 4 0 -51 42 4 0 -52 41 4 1 -53 45 5 0 -54 45 5 1 -55 49 5 0 -56 48 5 1 -57 47 5 0 -58 49 5 1 -59 46 5 1 -60 45 5 0 -61 49 5 1 -62 48 5 0 -63 47 5 1 -64 46 5 0 -65 47 5 0 -66 50 6 1 -67 51 6 1 -68 51 6 0 -69 54 6 1 -70 53 6 1 -71 51 6 0 -72 52 6 1 -73 54 6 0 -74 55 7 1 -75 56 7 1 -76 58 7 0 -77 59 7 1 -78 59 7 1 -79 58 7 0 -80 55 7 1 -81 56 7 1 -82 57 7 1 -83 58 7 1 -84 59 7 0 -85 55 7 1 -86 56 7 1 -87 57 7 1 -88 58 7 0 -89 59 7 1 -90 56 7 1 -91 60 8 1 -92 65 8 1 -93 67 8 1 -94 66 8 0 -95 63 8 1 -96 61 8 1 -97 69 8 1 -98 65 8 1 -99 64 8 1 -100 63 8 0 \ No newline at end of file diff --git a/doc/src/week38/week38.do.txt b/doc/src/week38/week38.do.txt index 3c8d4b781..fa4679c96 100644 --- a/doc/src/week38/week38.do.txt +++ b/doc/src/week38/week38.do.txt @@ -11,12 +11,913 @@ DATE: today !split -===== Thursday: ===== +===== Thursday September 17 ===== -Material will added during Wednesday 16 !split -===== Friday: Intro to Logistic Regression ===== +===== Ridge and LASSO Regression, reminder ===== + +The expression for the standard Mean Squared Error (MSE) which we used to define our cost function and the equations for the ordinary least squares (OLS) method, that is +our optimization problem is +!bt +\[ +{\displaystyle \min_{\bm{\beta}\in {\mathbb{R}}^{p}}}\frac{1}{n}\left\{\left(\bm{y}-\bm{X}\bm{\beta}\right)^T\left(\bm{y}-\bm{X}\bm{\beta}\right)\right\}. +\] +!et +or we can state it as +!bt +\[ +{\displaystyle \min_{\bm{\beta}\in +{\mathbb{R}}^{p}}}\frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\frac{1}{n}\vert\vert \bm{y}-\bm{X}\bm{\beta}\vert\vert_2^2, +\] +!et +where we have used the definition of a norm-2 vector, that is +!bt +\[ +\vert\vert \bm{x}\vert\vert_2 = \sqrt{\sum_i x_i^2}. +\] +!et + +By minimizing the above equation with respect to the parameters +$\bm{\beta}$ we could then obtain an analytical expression for the +parameters $\bm{\beta}$. We can add a regularization parameter $\lambda$ by +defining a new cost function to be optimized, that is + +!bt +\[ +{\displaystyle \min_{\bm{\beta}\in +{\mathbb{R}}^{p}}}\frac{1}{n}\vert\vert \bm{y}-\bm{X}\bm{\beta}\vert\vert_2^2+\lambda\vert\vert \bm{\beta}\vert\vert_2^2 +\] +!et + +which leads to the Ridge regression minimization problem where we +require that $\vert\vert \bm{\beta}\vert\vert_2^2\le t$, where $t$ is +a finite number larger than zero. By defining + +!bt +\[ +C(\bm{X},\bm{\beta})=\frac{1}{n}\vert\vert \bm{y}-\bm{X}\bm{\beta}\vert\vert_2^2+\lambda\vert\vert \bm{\beta}\vert\vert_1, +\] +!et + +we have a new optimization equation +!bt +\[ +{\displaystyle \min_{\bm{\beta}\in +{\mathbb{R}}^{p}}}\frac{1}{n}\vert\vert \bm{y}-\bm{X}\bm{\beta}\vert\vert_2^2+\lambda\vert\vert \bm{\beta}\vert\vert_1 +\] +!et +which leads to Lasso regression. Lasso stands for least absolute shrinkage and selection operator. + +Here we have defined the norm-1 as +!bt +\[ +\vert\vert \bm{x}\vert\vert_1 = \sum_i \vert x_i\vert. +\] +!et + + + +!split +===== Various steps in cross-validation ===== + +When the repetitive splitting of the data set is done randomly, +samples may accidently end up in a fast majority of the splits in +either training or test set. Such samples may have an unbalanced +influence on either model building or prediction evaluation. To avoid +this $k$-fold cross-validation structures the data splitting. The +samples are divided into $k$ more or less equally sized exhaustive and +mutually exclusive subsets. In turn (at each split) one of these +subsets plays the role of the test set while the union of the +remaining subsets constitutes the training set. Such a splitting +warrants a balanced representation of each sample in both training and +test set over the splits. Still the division into the $k$ subsets +involves a degree of randomness. This may be fully excluded when +choosing $k=n$. This particular case is referred to as leave-one-out +cross-validation (LOOCV). + +!split +===== How to set up the cross-validation for Ridge and/or Lasso ===== + +* Define a range of interest for the penalty parameter. + +* Divide the data set into training and test set comprising samples $\{1, \ldots, n\} \setminus i$ and $\{ i \}$, respectively. + +* Fit the linear regression model by means of ridge estimation for each $\lambda$ in the grid using the training set, and the corresponding estimate of the error variance $\bm{\sigma}_{-i}^2(\lambda)$, as +!bt +\begin{align*} +\bm{\beta}_{-i}(\lambda) & = ( \bm{X}_{-i, \ast}^{T} +\bm{X}_{-i, \ast} + \lambda \bm{I}_{pp})^{-1} +\bm{X}_{-i, \ast}^{T} \bm{y}_{-i} +\end{align*} +!et + +* Evaluate the prediction performance of these models on the test set by $\log\{L[y_i, \bm{X}_{i, \ast}; \bm{\beta}_{-i}(\lambda), \bm{\sigma}_{-i}^2(\lambda)]\}$. Or, by the prediction error $|y_i - \bm{X}_{i, \ast} \bm{\beta}_{-i}(\lambda)|$, the relative error, the error squared or the R2 score function. + +* Repeat the first three steps such that each sample plays the role of the test set once. + +* Average the prediction performances of the test sets at each grid point of the penalty bias/parameter. It is an estimate of the prediction performance of the model corresponding to this value of the penalty parameter on novel data. It is defined as +!bt +\begin{align*} +\frac{1}{n} \sum_{i = 1}^n \log\{L[y_i, \mathbf{X}_{i, \ast}; \bm{\beta}_{-i}(\lambda), \bm{\sigma}_{-i}^2(\lambda)]\}. +\end{align*} +!et + +!split +===== Cross-validation in brief ===== + +For the various values of $k$ + +o shuffle the dataset randomly. +o Split the dataset into $k$ groups. +o For each unique group: + o Decide which group to use as set for test data + o Take the remaining groups as a training data set + o Fit a model on the training set and evaluate it on the test set + o Retain the evaluation score and discard the model +o Summarize the model using the sample of model evaluation scores + + + +!split +===== Code Example for Cross-validation and $k$-fold Cross-validation ===== + +The code here uses Ridge regression with cross-validation (CV) resampling and $k$-fold CV in order to fit a specific polynomial. +!bc pycod +import numpy as np +import matplotlib.pyplot as plt +from sklearn.model_selection import KFold +from sklearn.linear_model import Ridge +from sklearn.model_selection import cross_val_score +from sklearn.preprocessing import PolynomialFeatures + +# A seed just to ensure that the random numbers are the same for every run. +# Useful for eventual debugging. +np.random.seed(3155) + +# Generate the data. +nsamples = 100 +x = np.random.randn(nsamples) +y = 3*x**2 + np.random.randn(nsamples) + +## Cross-validation on Ridge regression using KFold only + +# Decide degree on polynomial to fit +poly = PolynomialFeatures(degree = 6) + +# Decide which values of lambda to use +nlambdas = 500 +lambdas = np.logspace(-3, 5, nlambdas) + +# Initialize a KFold instance +k = 5 +kfold = KFold(n_splits = k) + +# Perform the cross-validation to estimate MSE +scores_KFold = np.zeros((nlambdas, k)) + +i = 0 +for lmb in lambdas: + ridge = Ridge(alpha = lmb) + j = 0 + for train_inds, test_inds in kfold.split(x): + xtrain = x[train_inds] + ytrain = y[train_inds] + + xtest = x[test_inds] + ytest = y[test_inds] + + Xtrain = poly.fit_transform(xtrain[:, np.newaxis]) + ridge.fit(Xtrain, ytrain[:, np.newaxis]) + + Xtest = poly.fit_transform(xtest[:, np.newaxis]) + ypred = ridge.predict(Xtest) + + scores_KFold[i,j] = np.sum((ypred - ytest[:, np.newaxis])**2)/np.size(ypred) + + j += 1 + i += 1 + + +estimated_mse_KFold = np.mean(scores_KFold, axis = 1) + +## Cross-validation using cross_val_score from sklearn along with KFold + +# kfold is an instance initialized above as: +# kfold = KFold(n_splits = k) + +estimated_mse_sklearn = np.zeros(nlambdas) +i = 0 +for lmb in lambdas: + ridge = Ridge(alpha = lmb) + + X = poly.fit_transform(x[:, np.newaxis]) + estimated_mse_folds = cross_val_score(ridge, X, y[:, np.newaxis], scoring='neg_mean_squared_error', cv=kfold) + + # cross_val_score return an array containing the estimated negative mse for every fold. + # we have to the the mean of every array in order to get an estimate of the mse of the model + estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds) + + i += 1 + +## Plot and compare the slightly different ways to perform cross-validation + +plt.figure() + +plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score') +plt.plot(np.log10(lambdas), estimated_mse_KFold, 'r--', label = 'KFold') + +plt.xlabel('log10(lambda)') +plt.ylabel('mse') + +plt.legend() + +plt.show() + +!ec + + + +!split +===== Bias-Variance tradeoff with Bootstrap ===== +!bc pycod +import matplotlib.pyplot as plt +import numpy as np +from sklearn.linear_model import LinearRegression, Ridge, Lasso +from sklearn.preprocessing import PolynomialFeatures +from sklearn.model_selection import train_test_split +from sklearn.pipeline import make_pipeline +from sklearn.utils import resample + +np.random.seed(2018) + +n = 40 +n_boostraps = 100 +maxdegree = 14 + + +# Make data set. +x = np.linspace(-3, 3, n).reshape(-1, 1) +y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape) +error = np.zeros(maxdegree) +bias = np.zeros(maxdegree) +variance = np.zeros(maxdegree) +polydegree = np.zeros(maxdegree) +x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2) + +for degree in range(maxdegree): + model = make_pipeline(PolynomialFeatures(degree=degree), LinearRegression(fit_intercept=False)) + y_pred = np.empty((y_test.shape[0], n_boostraps)) + for i in range(n_boostraps): + x_, y_ = resample(x_train, y_train) + y_pred[:, i] = model.fit(x_, y_).predict(x_test).ravel() + + polydegree[degree] = degree + error[degree] = np.mean( np.mean((y_test - y_pred)**2, axis=1, keepdims=True) ) + bias[degree] = np.mean( (y_test - np.mean(y_pred, axis=1, keepdims=True))**2 ) + variance[degree] = np.mean( np.var(y_pred, axis=1, keepdims=True) ) + print('Polynomial degree:', degree) + print('Error:', error[degree]) + print('Bias^2:', bias[degree]) + print('Var:', variance[degree]) + print('{} >= {} + {} = {}'.format(error[degree], bias[degree], variance[degree], bias[degree]+variance[degree])) + +plt.plot(polydegree, error, label='Error') +plt.plot(polydegree, bias, label='bias') +plt.plot(polydegree, variance, label='Variance') +plt.legend() +plt.show() + + + + +!ec + + +!split +===== Another Example from Scikit-Learn's Repository ===== +!bc pycod +""" +============================ +Underfitting vs. Overfitting +============================ + +This example demonstrates the problems of underfitting and overfitting and +how we can use linear regression with polynomial features to approximate +nonlinear functions. The plot shows the function that we want to approximate, +which is a part of the cosine function. In addition, the samples from the +real function and the approximations of different models are displayed. The +models have polynomial features of different degrees. We can see that a +linear function (polynomial with degree 1) is not sufficient to fit the +training samples. This is called **underfitting**. A polynomial of degree 4 +approximates the true function almost perfectly. However, for higher degrees +the model will **overfit** the training data, i.e. it learns the noise of the +training data. +We evaluate quantitatively **overfitting** / **underfitting** by using +cross-validation. We calculate the mean squared error (MSE) on the validation +set, the higher, the less likely the model generalizes correctly from the +training data. +""" + +print(__doc__) + +import numpy as np +import matplotlib.pyplot as plt +from sklearn.pipeline import Pipeline +from sklearn.preprocessing import PolynomialFeatures +from sklearn.linear_model import LinearRegression +from sklearn.model_selection import cross_val_score + + +def true_fun(X): + return np.cos(1.5 * np.pi * X) + +np.random.seed(0) + +n_samples = 30 +degrees = [1, 4, 15] + +X = np.sort(np.random.rand(n_samples)) +y = true_fun(X) + np.random.randn(n_samples) * 0.1 + +plt.figure(figsize=(14, 5)) +for i in range(len(degrees)): + ax = plt.subplot(1, len(degrees), i + 1) + plt.setp(ax, xticks=(), yticks=()) + + polynomial_features = PolynomialFeatures(degree=degrees[i], + include_bias=False) + linear_regression = LinearRegression() + pipeline = Pipeline([("polynomial_features", polynomial_features), + ("linear_regression", linear_regression)]) + pipeline.fit(X[:, np.newaxis], y) + + # Evaluate the models using crossvalidation + scores = cross_val_score(pipeline, X[:, np.newaxis], y, + scoring="neg_mean_squared_error", cv=10) + + X_test = np.linspace(0, 1, 100) + plt.plot(X_test, pipeline.predict(X_test[:, np.newaxis]), label="Model") + plt.plot(X_test, true_fun(X_test), label="True function") + plt.scatter(X, y, edgecolor='b', s=20, label="Samples") + plt.xlabel("x") + plt.ylabel("y") + plt.xlim((0, 1)) + plt.ylim((-2, 2)) + plt.legend(loc="best") + plt.title("Degree {}\nMSE = {:.2e}(+/- {:.2e})".format( + degrees[i], -scores.mean(), scores.std())) +plt.show() +!ec + + + +!split +===== Cross-validation with Ridge ===== +!bc pycod +import numpy as np +import matplotlib.pyplot as plt +from sklearn.model_selection import KFold +from sklearn.linear_model import Ridge +from sklearn.model_selection import cross_val_score +from sklearn.preprocessing import PolynomialFeatures + +# A seed just to ensure that the random numbers are the same for every run. +np.random.seed(3155) +# Generate the data. +n = 100 +x = np.linspace(-3, 3, n).reshape(-1, 1) +y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape) +# Decide degree on polynomial to fit +poly = PolynomialFeatures(degree = 10) + +# Decide which values of lambda to use +nlambdas = 500 +lambdas = np.logspace(-3, 5, nlambdas) +# Initialize a KFold instance +k = 5 +kfold = KFold(n_splits = k) +estimated_mse_sklearn = np.zeros(nlambdas) +i = 0 +for lmb in lambdas: + ridge = Ridge(alpha = lmb) + estimated_mse_folds = cross_val_score(ridge, x, y, scoring='neg_mean_squared_error', cv=kfold) + estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds) + i += 1 +plt.figure() +plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score') +plt.xlabel('log10(lambda)') +plt.ylabel('MSE') +plt.legend() +plt.show() + + +!ec + + + + + + + + + + +!split +===== The Ising model ===== + +The one-dimensional Ising model with nearest neighbor interaction, no +external field and a constant coupling constant $J$ is given by + +!bt +\begin{align} + H = -J \sum_{k}^L s_k s_{k + 1}, +\end{align} +!et + +where $s_i \in \{-1, 1\}$ and $s_{N + 1} = s_1$. The number of spins +in the system is determined by $L$. For the one-dimensional system +there is no phase transition. + +We will look at a system of $L = 40$ spins with a coupling constant of +$J = 1$. To get enough training data we will generate 10000 states +with their respective energies. + + +!bc pycod +import numpy as np +import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import make_axes_locatable +import seaborn as sns +import scipy.linalg as scl +from sklearn.model_selection import train_test_split +import tqdm +sns.set(color_codes=True) +cmap_args=dict(vmin=-1., vmax=1., cmap='seismic') + +L = 40 +n = int(1e4) + +spins = np.random.choice([-1, 1], size=(n, L)) +J = 1.0 + +energies = np.zeros(n) + +for i in range(n): + energies[i] = - J * np.dot(spins[i], np.roll(spins[i], 1)) +!ec + +Here we use ordinary least squares +regression to predict the energy for the nearest neighbor +one-dimensional Ising model on a ring, i.e., the endpoints wrap +around. We will use linear regression to fit a value for +the coupling constant to achieve this. + +!split +===== Reformulating the problem to suit regression ===== + +A more general form for the one-dimensional Ising model is + +!bt +\begin{align} + H = - \sum_j^L \sum_k^L s_j s_k J_{jk}. +\end{align} +!et + +Here we allow for interactions beyond the nearest neighbors and a state dependent +coupling constant. This latter expression can be formulated as +a matrix-product +!bt +\begin{align} + \bm{H} = \bm{X} J, +\end{align} +!et + +where $X_{jk} = s_j s_k$ and $J$ is a matrix which consists of the +elements $-J_{jk}$. This form of writing the energy fits perfectly +with the form utilized in linear regression, that is + +!bt +\begin{align} + \bm{y} = \bm{X}\bm{\beta} + \bm{\epsilon}, +\end{align} +!et + +We split the data in training and test data as discussed in the previous example + +!bc pycod +X = np.zeros((n, L ** 2)) +for i in range(n): + X[i] = np.outer(spins[i], spins[i]).ravel() +y = energies +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) +!ec + +!split +===== Linear regression ===== + +In the ordinary least squares method we choose the cost function + +!bt +\begin{align} + C(\bm{X}, \bm{\beta})= \frac{1}{n}\left\{(\bm{X}\bm{\beta} - \bm{y})^T(\bm{X}\bm{\beta} - \bm{y})\right\}. +\end{align} +!et + +We then find the extremal point of $C$ by taking the derivative with respect to $\bm{\beta}$ as discussed above. +This yields the expression for $\bm{\beta}$ to be + +!bt +\[ + \bm{\beta} = \frac{\bm{X}^T \bm{y}}{\bm{X}^T \bm{X}}, +\] +!et + +which immediately imposes some requirements on $\bm{X}$ as there must exist +an inverse of $\bm{X}^T \bm{X}$. If the expression we are modeling contains an +intercept, i.e., a constant term, we must make sure that the +first column of $\bm{X}$ consists of $1$. We do this here + +!bc pycod +X_train_own = np.concatenate( + (np.ones(len(X_train))[:, np.newaxis], X_train), + axis=1 +) +X_test_own = np.concatenate( + (np.ones(len(X_test))[:, np.newaxis], X_test), + axis=1 +) +!ec + +!bc pycod +def ols_inv(x: np.ndarray, y: np.ndarray) -> np.ndarray: + return scl.inv(x.T @ x) @ (x.T @ y) +beta = ols_inv(X_train_own, y_train) +!ec + + +!split +===== Singular Value decomposition ===== + +Doing the inversion directly turns out to be a bad idea since the matrix +$\bm{X}^T\bm{X}$ is singular. An alternative approach is to use the _singular +value decomposition_. Using the definition of the Moore-Penrose +pseudoinverse we can write the equation for $\bm{\beta}$ as + +!bt +\[ + \bm{\beta} = \bm{X}^{+}\bm{y}, +\] +!et + +where the pseudoinverse of $\bm{X}$ is given by + +!bt +\[ + \bm{X}^{+} = \frac{\bm{X}^T}{\bm{X}^T\bm{X}}. +\] +!et + +Using singular value decomposition we can decompose the matrix $\bm{X} = \bm{U}\bm{\Sigma} \bm{V}^T$, +where $\bm{U}$ and $\bm{V}$ are orthogonal(unitary) matrices and $\bm{\Sigma}$ contains the singular values (more details below). +where $X^{+} = V\Sigma^{+} U^T$. This reduces the equation for +$\omega$ to +!bt +\begin{align} + \bm{\beta} = \bm{V}\bm{\Sigma}^{+} \bm{U}^T \bm{y}. +\end{align} +!et + +Note that solving this equation by actually doing the pseudoinverse +(which is what we will do) is not a good idea as this operation scales +as $\mathcal{O}(n^3)$, where $n$ is the number of elements in a +general matrix. Instead, doing $QR$-factorization and solving the +linear system as an equation would reduce this down to +$\mathcal{O}(n^2)$ operations. + + +!bc pycod +def ols_svd(x: np.ndarray, y: np.ndarray) -> np.ndarray: + u, s, v = scl.svd(x) + return v.T @ scl.pinv(scl.diagsvd(s, u.shape[0], v.shape[0])) @ u.T @ y +!ec + +!bc pycod +beta = ols_svd(X_train_own,y_train) +!ec + +When extracting the $J$-matrix we need to make sure that we remove the intercept, as is done here + +!bc pycod +J = beta[1:].reshape(L, L) +!ec + +A way of looking at the coefficients in $J$ is to plot the matrices as images. + + +!bc pycod +fig = plt.figure(figsize=(20, 14)) +im = plt.imshow(J, **cmap_args) +plt.title("OLS", fontsize=18) +plt.xticks(fontsize=18) +plt.yticks(fontsize=18) +cb = fig.colorbar(im) +cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18) +plt.show() +!ec +It is interesting to note that OLS +considers both $J_{j, j + 1} = -0.5$ and $J_{j, j - 1} = -0.5$ as +valid matrix elements for $J$. +In our discussion below on hyperparameters and Ridge and Lasso regression we will see that +this problem can be removed, partly and only with Lasso regression. + +In this case our matrix inversion was actually possible. The obvious question now is what is the mathematics behind the SVD? + + + + + +!split +===== The one-dimensional Ising model ===== + +Let us bring back the Ising model again, but now with an additional +focus on Ridge and Lasso regression as well. We repeat some of the +basic parts of the Ising model and the setup of the training and test +data. The one-dimensional Ising model with nearest neighbor +interaction, no external field and a constant coupling constant $J$ is +given by + +!bt +\begin{align} + H = -J \sum_{k}^L s_k s_{k + 1}, +\end{align} +!et +where $s_i \in \{-1, 1\}$ and $s_{N + 1} = s_1$. The number of spins in the system is determined by $L$. For the one-dimensional system there is no phase transition. + +We will look at a system of $L = 40$ spins with a coupling constant of $J = 1$. To get enough training data we will generate 10000 states with their respective energies. + + +!bc pycod +import numpy as np +import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import make_axes_locatable +import seaborn as sns +import scipy.linalg as scl +from sklearn.model_selection import train_test_split +import sklearn.linear_model as skl +import tqdm +sns.set(color_codes=True) +cmap_args=dict(vmin=-1., vmax=1., cmap='seismic') + +L = 40 +n = int(1e4) + +spins = np.random.choice([-1, 1], size=(n, L)) +J = 1.0 + +energies = np.zeros(n) + +for i in range(n): + energies[i] = - J * np.dot(spins[i], np.roll(spins[i], 1)) +!ec + +A more general form for the one-dimensional Ising model is + +!bt +\begin{align} + H = - \sum_j^L \sum_k^L s_j s_k J_{jk}. +\end{align} +!et + +Here we allow for interactions beyond the nearest neighbors and a more +adaptive coupling matrix. This latter expression can be formulated as +a matrix-product on the form +!bt +\begin{align} + H = X J, +\end{align} +!et + +where $X_{jk} = s_j s_k$ and $J$ is the matrix consisting of the +elements $-J_{jk}$. This form of writing the energy fits perfectly +with the form utilized in linear regression, viz. +!bt +\begin{align} + \bm{y} = \bm{X}\bm{\beta} + \bm{\epsilon}. +\end{align} +!et +We organize the data as we did above +!bc pycod +X = np.zeros((n, L ** 2)) +for i in range(n): + X[i] = np.outer(spins[i], spins[i]).ravel() +y = energies +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.96) + +X_train_own = np.concatenate( + (np.ones(len(X_train))[:, np.newaxis], X_train), + axis=1 +) + +X_test_own = np.concatenate( + (np.ones(len(X_test))[:, np.newaxis], X_test), + axis=1 +) +!ec + +We will do all fitting with _Scikit-Learn_, + +!bc pycod +clf = skl.LinearRegression().fit(X_train, y_train) +!ec +When extracting the $J$-matrix we make sure to remove the intercept +!bc pycod +J_sk = clf.coef_.reshape(L, L) +!ec +And then we plot the results +!bc pycod +fig = plt.figure(figsize=(20, 14)) +im = plt.imshow(J_sk, **cmap_args) +plt.title("LinearRegression from Scikit-learn", fontsize=18) +plt.xticks(fontsize=18) +plt.yticks(fontsize=18) +cb = fig.colorbar(im) +cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18) +plt.show() +!ec +The results perfectly with our previous discussion where we used our own code. + +!split +===== Ridge regression ===== + +Having explored the ordinary least squares we move on to ridge +regression. In ridge regression we include a _regularizer_. This +involves a new cost function which leads to a new estimate for the +weights $\bm{\beta}$. This results in a penalized regression problem. The +cost function is given by + +!bt +\begin{align} + C(\bm{X}, \bm{\beta}; \lambda) = (\bm{X}\bm{\beta} - \bm{y})^T(\bm{X}\bm{\beta} - \bm{y}) + \lambda \bm{\beta}^T\bm{\beta}. +\end{align} +!et +!bc pycod +_lambda = 0.1 +clf_ridge = skl.Ridge(alpha=_lambda).fit(X_train, y_train) +J_ridge_sk = clf_ridge.coef_.reshape(L, L) +fig = plt.figure(figsize=(20, 14)) +im = plt.imshow(J_ridge_sk, **cmap_args) +plt.title("Ridge from Scikit-learn", fontsize=18) +plt.xticks(fontsize=18) +plt.yticks(fontsize=18) +cb = fig.colorbar(im) +cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18) + +plt.show() +!ec + +!split +===== LASSO regression ===== + +In the _Least Absolute Shrinkage and Selection Operator_ (LASSO)-method we get a third cost function. + +!bt +\begin{align} + C(\bm{X}, \bm{\beta}; \lambda) = (\bm{X}\bm{\beta} - \bm{y})^T(\bm{X}\bm{\beta} - \bm{y}) + \lambda \sqrt{\bm{\beta}^T\bm{\beta}}. +\end{align} +!et + +Finding the extremal point of this cost function is not so straight-forward as in least squares and ridge. We will therefore rely solely on the function ``Lasso`` from _Scikit-Learn_. + +!bc pycod +clf_lasso = skl.Lasso(alpha=_lambda).fit(X_train, y_train) +J_lasso_sk = clf_lasso.coef_.reshape(L, L) +fig = plt.figure(figsize=(20, 14)) +im = plt.imshow(J_lasso_sk, **cmap_args) +plt.title("Lasso from Scikit-learn", fontsize=18) +plt.xticks(fontsize=18) +plt.yticks(fontsize=18) +cb = fig.colorbar(im) +cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18) + +plt.show() +!ec + +It is quite striking how LASSO breaks the symmetry of the coupling +constant as opposed to ridge and OLS. We get a sparse solution with +$J_{j, j + 1} = -1$. + + + +!split +===== Performance as function of the regularization parameter ===== + +We see how the different models perform for a different set of values for $\lambda$. + + +!bc pycod +lambdas = np.logspace(-4, 5, 10) + +train_errors = { + "ols_sk": np.zeros(lambdas.size), + "ridge_sk": np.zeros(lambdas.size), + "lasso_sk": np.zeros(lambdas.size) +} + +test_errors = { + "ols_sk": np.zeros(lambdas.size), + "ridge_sk": np.zeros(lambdas.size), + "lasso_sk": np.zeros(lambdas.size) +} + +plot_counter = 1 + +fig = plt.figure(figsize=(32, 54)) + +for i, _lambda in enumerate(tqdm.tqdm(lambdas)): + for key, method in zip( + ["ols_sk", "ridge_sk", "lasso_sk"], + [skl.LinearRegression(), skl.Ridge(alpha=_lambda), skl.Lasso(alpha=_lambda)] + ): + method = method.fit(X_train, y_train) + + train_errors[key][i] = method.score(X_train, y_train) + test_errors[key][i] = method.score(X_test, y_test) + + omega = method.coef_.reshape(L, L) + + plt.subplot(10, 5, plot_counter) + plt.imshow(omega, **cmap_args) + plt.title(r"%s, $\lambda = %.4f$" % (key, _lambda)) + plot_counter += 1 + +plt.show() +!ec + +We see that LASSO reaches a good solution for low +values of $\lambda$, but will "wither" when we increase $\lambda$ too +much. Ridge is more stable over a larger range of values for +$\lambda$, but eventually also fades away. + +!split +===== Finding the optimal value of $\lambda$ ===== + +To determine which value of $\lambda$ is best we plot the accuracy of +the models when predicting the training and the testing set. We expect +the accuracy of the training set to be quite good, but if the accuracy +of the testing set is much lower this tells us that we might be +subject to an overfit model. The ideal scenario is an accuracy on the +testing set that is close to the accuracy of the training set. + + +!bc pycod +fig = plt.figure(figsize=(20, 14)) + +colors = { + "ols_sk": "r", + "ridge_sk": "y", + "lasso_sk": "c" +} + +for key in train_errors: + plt.semilogx( + lambdas, + train_errors[key], + colors[key], + label="Train {0}".format(key), + linewidth=4.0 + ) + +for key in test_errors: + plt.semilogx( + lambdas, + test_errors[key], + colors[key] + "--", + label="Test {0}".format(key), + linewidth=4.0 + ) +plt.legend(loc="best", fontsize=18) +plt.xlabel(r"$\lambda$", fontsize=18) +plt.ylabel(r"$R^2$", fontsize=18) +plt.tick_params(labelsize=18) +plt.show() +!ec + +From the above figure we can see that LASSO with $\lambda = 10^{-2}$ +achieves a very good accuracy on the test set. This by far surpasses the +other models for all values of $\lambda$. + + + + + + + + + + +!split +===== Friday September 18: Intro to Logistic Regression ===== !split @@ -550,3 +1451,4 @@ plt.show() +