From 1925c4cce4a5c6d9d55714ef75bd2e691247f911 Mon Sep 17 00:00:00 2001 From: mhjensen Date: Thu, 10 Sep 2020 11:30:21 +0200 Subject: [PATCH] added dot files --- .../Regression/html/._Regression-bs115.html | 540 +++++++++++++++ .../Regression/html/._Regression-bs116.html | 578 ++++++++++++++++ .../Regression/html/._Regression-bs117.html | 624 ++++++++++++++++++ .../Regression/html/._Regression-bs118.html | 524 +++++++++++++++ .../Regression/html/._Regression-bs119.html | 526 +++++++++++++++ .../Regression/html/._Regression-bs120.html | 541 +++++++++++++++ .../Regression/html/._Regression-bs121.html | 538 +++++++++++++++ 7 files changed, 3871 insertions(+) create mode 100644 doc/pub/Regression/html/._Regression-bs115.html create mode 100644 doc/pub/Regression/html/._Regression-bs116.html create mode 100644 doc/pub/Regression/html/._Regression-bs117.html create mode 100644 doc/pub/Regression/html/._Regression-bs118.html create mode 100644 doc/pub/Regression/html/._Regression-bs119.html create mode 100644 doc/pub/Regression/html/._Regression-bs120.html create mode 100644 doc/pub/Regression/html/._Regression-bs121.html diff --git a/doc/pub/Regression/html/._Regression-bs115.html b/doc/pub/Regression/html/._Regression-bs115.html new file mode 100644 index 000000000..e53c538b1 --- /dev/null +++ b/doc/pub/Regression/html/._Regression-bs115.html @@ -0,0 +1,540 @@ + + + + + + + + +Data Analysis and Machine Learning: Linear Regression and more Advanced Regression Analysis + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ +

 

 

 

+ + + + +

Linear regression

+ +

+In the ordinary least squares method we choose the cost function + +$$ +\begin{align} + C(\boldsymbol{X}, \boldsymbol{\beta})= \frac{1}{n}\left\{(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})^T(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})\right\}. +\tag{25} +\end{align} +$$ + +

+We then find the extremal point of \( C \) by taking the derivative with respect to \( \boldsymbol{\beta} \) as discussed above. +This yields the expression for \( \boldsymbol{\beta} \) to be + +$$ + \boldsymbol{\beta} = \frac{\boldsymbol{X}^T \boldsymbol{y}}{\boldsymbol{X}^T \boldsymbol{X}}, +$$ + +

+which immediately imposes some requirements on \( \boldsymbol{X} \) as there must exist +an inverse of \( \boldsymbol{X}^T \boldsymbol{X} \). If the expression we are modeling contains an +intercept, i.e., a constant term, we must make sure that the +first column of \( \boldsymbol{X} \) consists of \( 1 \). We do this here + +

+ + +

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
+)
+
+

+ + +

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)
+
+

+

+ +

+ + +
+ + + + + + + +
+ +
+ + + + + + diff --git a/doc/pub/Regression/html/._Regression-bs116.html b/doc/pub/Regression/html/._Regression-bs116.html new file mode 100644 index 000000000..862e578f9 --- /dev/null +++ b/doc/pub/Regression/html/._Regression-bs116.html @@ -0,0 +1,578 @@ + + + + + + + + +Data Analysis and Machine Learning: Linear Regression and more Advanced Regression Analysis + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ +

 

 

 

+ + + + +

Singular Value decomposition

+ +

+Doing the inversion directly turns out to be a bad idea since the matrix +\( \boldsymbol{X}^T\boldsymbol{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 \( \boldsymbol{\beta} \) as + +$$ + \boldsymbol{\beta} = \boldsymbol{X}^{+}\boldsymbol{y}, +$$ + +

+where the pseudoinverse of \( \boldsymbol{X} \) is given by + +$$ + \boldsymbol{X}^{+} = \frac{\boldsymbol{X}^T}{\boldsymbol{X}^T\boldsymbol{X}}. +$$ + +

+Using singular value decomposition we can decompose the matrix \( \boldsymbol{X} = \boldsymbol{U}\boldsymbol{\Sigma} \boldsymbol{V}^T \), +where \( \boldsymbol{U} \) and \( \boldsymbol{V} \) are orthogonal(unitary) matrices and \( \boldsymbol{\Sigma} \) contains the singular values (more details below). +where \( X^{+} = V\Sigma^{+} U^T \). This reduces the equation for +\( \omega \) to +$$ +\begin{align} + \boldsymbol{\beta} = \boldsymbol{V}\boldsymbol{\Sigma}^{+} \boldsymbol{U}^T \boldsymbol{y}. +\tag{26} +\end{align} +$$ + +

+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. + +

+ + +

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
+
+

+ + +

beta = ols_svd(X_train_own,y_train)
+
+

+When extracting the \( J \)-matrix we need to make sure that we remove the intercept, as is done here + +

+ + +

J = beta[1:].reshape(L, L)
+
+

+A way of looking at the coefficients in \( J \) is to plot the matrices as images. + +

+ + +

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()
+
+

+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? + +

+

+ +

+ + +
+ + + + + + + +
+ +
+ + + + + + diff --git a/doc/pub/Regression/html/._Regression-bs117.html b/doc/pub/Regression/html/._Regression-bs117.html new file mode 100644 index 000000000..01a436609 --- /dev/null +++ b/doc/pub/Regression/html/._Regression-bs117.html @@ -0,0 +1,624 @@ + + + + + + + + +Data Analysis and Machine Learning: Linear Regression and more Advanced Regression Analysis + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ +

 

 

 

+ + + + +

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 + +$$ +\begin{align} + H = -J \sum_{k}^L s_k s_{k + 1}, +\tag{27} +\end{align} +$$ + +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. + +

+ + +

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))
+
+

+A more general form for the one-dimensional Ising model is + +$$ +\begin{align} + H = - \sum_j^L \sum_k^L s_j s_k J_{jk}. +\tag{28} +\end{align} +$$ + +

+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 +$$ +\begin{align} + H = X J, +\tag{29} +\end{align} +$$ + +

+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. +$$ +\begin{align} + \boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}. +\tag{30} +\end{align} +$$ + +We organize the data as we did above +

+ + +

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
+)
+
+

+We will do all fitting with Scikit-Learn, + +

+ + +

clf = skl.LinearRegression().fit(X_train, y_train)
+
+

+When extracting the \( J \)-matrix we make sure to remove the intercept +

+ + +

J_sk = clf.coef_.reshape(L, L)
+
+

+And then we plot the results +

+ + +

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()
+
+

+The results perfectly with our previous discussion where we used our own code. + +

+

+ +

+ + +
+ + + + + + + +
+ +
+ + + + + + diff --git a/doc/pub/Regression/html/._Regression-bs118.html b/doc/pub/Regression/html/._Regression-bs118.html new file mode 100644 index 000000000..9c35368ee --- /dev/null +++ b/doc/pub/Regression/html/._Regression-bs118.html @@ -0,0 +1,524 @@ + + + + + + + + +Data Analysis and Machine Learning: Linear Regression and more Advanced Regression Analysis + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ +

 

 

 

+ + + + +

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 \( \boldsymbol{\beta} \). This results in a penalized regression problem. The +cost function is given by + +$$ +\begin{align} + C(\boldsymbol{X}, \boldsymbol{\beta}; \lambda) = (\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})^T(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y}) + \lambda \boldsymbol{\beta}^T\boldsymbol{\beta}. +\tag{31} +\end{align} +$$ + +

+ + +

_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()
+
+

+

+ +

+ + +
+ + + + + + + +
+ +
+ + + + + + diff --git a/doc/pub/Regression/html/._Regression-bs119.html b/doc/pub/Regression/html/._Regression-bs119.html new file mode 100644 index 000000000..a8085cc35 --- /dev/null +++ b/doc/pub/Regression/html/._Regression-bs119.html @@ -0,0 +1,526 @@ + + + + + + + + +Data Analysis and Machine Learning: Linear Regression and more Advanced Regression Analysis + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ +

 

 

 

+ + + + +

LASSO regression

+ +

+In the Least Absolute Shrinkage and Selection Operator (LASSO)-method we get a third cost function. + +$$ +\begin{align} + C(\boldsymbol{X}, \boldsymbol{\beta}; \lambda) = (\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})^T(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y}) + \lambda \sqrt{\boldsymbol{\beta}^T\boldsymbol{\beta}}. +\tag{32} +\end{align} +$$ + +

+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. + +

+ + +

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()
+
+

+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 \). + +

+

+ +

+ + +
+ + + + + + + +
+ +
+ + + + + + diff --git a/doc/pub/Regression/html/._Regression-bs120.html b/doc/pub/Regression/html/._Regression-bs120.html new file mode 100644 index 000000000..2756d5d50 --- /dev/null +++ b/doc/pub/Regression/html/._Regression-bs120.html @@ -0,0 +1,541 @@ + + + + + + + + +Data Analysis and Machine Learning: Linear Regression and more Advanced Regression Analysis + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ +

 

 

 

+ + + + +

Performance as function of the regularization parameter

+ +

+We see how the different models perform for a different set of values for \( \lambda \). + +

+ + +

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()
+
+

+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. + +

+

+ +

+ + +
+ + + + + + + +
+ +
+ + + + + + diff --git a/doc/pub/Regression/html/._Regression-bs121.html b/doc/pub/Regression/html/._Regression-bs121.html new file mode 100644 index 000000000..ff734262c --- /dev/null +++ b/doc/pub/Regression/html/._Regression-bs121.html @@ -0,0 +1,538 @@ + + + + + + + + +Data Analysis and Machine Learning: Linear Regression and more Advanced Regression Analysis + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ +

 

 

 

+ + + + +

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. + +

+ + +

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()
+
+

+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 \). + +

+ +

+ +

+ + +
+ + + + + + + +
+ +
+ + + + + +