updating many files

doc/LectureNotes/_build/html/_sources/schedule.md

@@ -78,6 +78,7 @@ For the reading assignments we use the following abbreviations:
 ###  Week 38 September 20-24
 - Lab Wednesday: Work on  Project 1
 - Lecture Thursday:  Classification problems and Logistic Regression, from binary cases to several categories
+  - Video of Lecture at https://www.uio.no/studier/emner/matnat/fys/FYS-STK4155/h21/forelesningsvideoer/LectureSeptember23.mp4?vrtx=view-as-webpage
 - Lecture Friday: Logistic Regression and gradient optimization
 
 - Reading recommendations:

doc/LectureNotes/_build/html/schedule.html

@@ -489,7 +489,11 @@
 <h3>Week 38 September 20-24<a class="headerlink" href="#week-38-september-20-24" title="Permalink to this headline">¶</a></h3>
 <ul class="simple">
 <li><p>Lab Wednesday: Work on  Project 1</p></li>
-<li><p>Lecture Thursday:  Classification problems and Logistic Regression, from binary cases to several categories</p></li>
+<li><p>Lecture Thursday:  Classification problems and Logistic Regression, from binary cases to several categories</p>
+<ul>
+<li><p>Video of Lecture at <a class="reference external" href="https://www.uio.no/studier/emner/matnat/fys/FYS-STK4155/h21/forelesningsvideoer/LectureSeptember23.mp4?vrtx=view-as-webpage">https://www.uio.no/studier/emner/matnat/fys/FYS-STK4155/h21/forelesningsvideoer/LectureSeptember23.mp4?vrtx=view-as-webpage</a></p></li>
+</ul>
+</li>
 <li><p>Lecture Friday: Logistic Regression and gradient optimization</p></li>
 <li><p>Reading recommendations:</p>
 <ul>

doc/LectureNotes/_build/jupyter_execute/chapter3.py

@@ -1167,6 +1167,969 @@ plt.ylabel('MSE')
 plt.legend()
 plt.show()
 
+## More on Rescaling data
+
+We end this chapter by adding some words on scaling and how to deal with the intercept for regression cases.
+
+When you are comparing your own code with for example **Scikit-Learn**'s
+library, there are some technicalities to keep in mind.  The examples
+here demonstrate some of these aspects with potential pitfalls.
+
+The discussion here focuses on the role of the intercept, how we can
+set up the design matrix, what scaling we should use and other topics
+which tend  confuse us.
+
+The intercept can be interpreted as the expected value of our
+target/output variables when all other predictors are set to zero.
+Thus, if we cannot assume that the expected outputs/targets are zero
+when all predictors are zero (the columns in the design matrix), it
+may be a bad idea to implement a model which penalizes the intercept.
+Furthermore, in for example Ridge and Lasso regression, the default solutions
+from the library **Scikit-Learn** (when not shrinking $\beta_0$) for the unknown parameters
+$\boldsymbol{\beta}$, are derived under the assumption that both $\boldsymbol{y}$ and
+$\boldsymbol{X}$ are zero centered, that is we subtract the mean values.
+
+
+If our predictors represent different scales, then it is important to
+standardize the design matrix $\boldsymbol{X}$ by subtracting the mean of each
+column from the corresponding column and dividing the column with its
+standard deviation. Most machine learning libraries do this as a default. This means that if you compare your code with the results from a given library,
+the results may differ. 
+
+The
+[Standadscaler](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.StandardScaler.html)
+function in **Scikit-Learn** does this for us.  For the data sets we
+have been studying in our various examples, the data are in many cases
+already scaled and there is no need to scale them. You as a user of different machine learning algorithms, should always perform  a
+survey of your data, with a critical assessment of them in case you need to scale the data.
+
+If you need to scale the data, not doing so will give an *unfair*
+penalization of the parameters since their magnitude depends on the
+scale of their corresponding predictor.
+
+Suppose as an example that you 
+you have an input variable given by the heights of different persons.
+Human height might be measured in inches or meters or
+kilometers. If measured in kilometers, a standard linear regression
+model with this predictor would probably give a much bigger
+coefficient term, than if measured in millimeters.
+This can clearly lead to problems in evaluating the cost/loss functions.
+
+
+
+Keep in mind that when you transform your data set before training a model, the same transformation needs to be done
+on your eventual new data set  before making a prediction. If we translate this into a Python code, it would could be implemented as follows
+
+"""
+#Model training, we compute the mean value of y and X
+y_train_mean = np.mean(y_train)
+X_train_mean = np.mean(X_train,axis=0)
+X_train = X_train - X_train_mean
+y_train = y_train - y_train_mean
+
+# The we fit our model with the training data
+trained_model = some_model.fit(X_train,y_train)
+
+
+#Model prediction, we need also to transform our data set used for the prediction.
+X_test = X_test - X_train_mean #Use mean from training data
+y_pred = trained_model(X_test)
+y_pred = y_pred + y_train_mean
+"""
+
+Let us try to understand what this may imply mathematically when we
+subtract the mean values, also known as *zero centering*. For
+simplicity, we will focus on  ordinary regression, as done in the above example.
+
+The cost/loss function  for regression is
+
+$$
+C(\beta_0, \beta_1, ... , \beta_{p-1}) = \frac{1}{n}\sum_{i=0}^{n} \left(y_i - \beta_0 - \sum_{j=1}^{p-1} X_{ij}\beta_j\right)^2,.
+$$
+
+Recall also that we use the squared value since this leads to an increase of the penalty for higher differences between predicted and output/target values.
+
+What we have done is to single out the $\beta_0$ term in the definition of the mean squared error (MSE).
+The design matrix
+$X$ does in this case not contain any intercept column.
+When we take the derivative with respect to $\beta_0$, we want the derivative to obey
+
+$$
+\frac{\partial C}{\partial \beta_j} = 0,
+$$
+
+for all $j$. For $\beta_0$ we have
+
+$$
+\frac{\partial C}{\partial \beta_0} = -\frac{2}{n}\sum_{i=0}^{n-1} \left(y_i - \beta_0 - \sum_{j=1}^{p-1} X_{ij} \beta_j\right).
+$$
+
+Multiplying away the constant $2/n$, we obtain
+
+$$
+\sum_{i=0}^{n-1} \beta_0 = \sum_{i=0}^{n-1}y_i - \sum_{i=0}^{n-1} \sum_{j=1}^{p-1} X_{ij} \beta_j.
+$$
+
+We assume 
+that every column of $\boldsymbol{X}$ is centered, which we can do by subtracting the mean,
+
+X = X - np.mean(X,axis=0)
+
+This means that we need to rewrite $X_{ij}$ as $\tilde{X}_{ij}=X_{ij}-\mu_j$, where
+
+$$
+\mu_j = \frac{1}{n}\sum_{i=0}^{n-1}X_{ij}.
+$$
+
+Let us special first to the case where we have only two parameters $\beta_0$ and $\beta_1$.
+Our result for $\beta_0$ simplifies then to
+
+$$
+n\beta_0 = \sum_{i=0}^{n-1}y_i - \sum_{i=0}^{n-1} X_{i1} \beta_1.
+$$
+
+Assuming that the matrix elements $X_{i1}$ are centered, what we have is
+
+$$
+\beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \beta_1\frac{1}{n}\sum_{i=0}^{n-1} \left(X_{i1}-\mu_{1}\right),
+$$
+
+where
+
+$$
+\mu_1=\frac{1}{n}\sum_{i=0}^{n-1} (X_{i1},
+$$
+
+and if we define the mean value of the outputs as
+
+$$
+\mu_y=\frac{1}{n}\sum_{i=0}^{n-1}y_i,
+$$
+
+we have
+
+$$
+\beta_0 = \mu_y - \beta_1\frac{1}{n}\sum_{i=0}^{n-1} (X_{i1}-\mu_{1}),
+$$
+
+and it is easy to see that the last sum equals zero!  This means that we have
+
+$$
+\beta_0 = \mu_y,
+$$
+
+if the columns of the design matrix are centered. It is straight forward to generalize this results to more values of $\beta$.
+We have thus
+
+$$
+\beta_0 = \frac{1}{n}\sum_{i=0}^{n-1} y_i = \overline{\boldsymbol{y}},
+$$
+
+the average value of $\boldsymbol{y}$.
+
+Replacing $y_i$ with $y_i - \beta_0 = y_i - \overline{\boldsymbol{y}}$ and centering also our design matrix results in a cost function (in vector-matrix disguise)
+
+$$
+C(\boldsymbol{\beta}) = (\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\beta})^T(\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\beta}).
+$$
+
+If we minimize with respect to $\boldsymbol{\beta}$ we have then
+
+$$
+\hat{\boldsymbol{\beta}} = (\tilde{X}^T\tilde{X})^{-1}\tilde{X}^T\boldsymbol{\tilde{y}},
+$$
+
+where $\boldsymbol{\tilde{y}} = \boldsymbol{y} - \overline{\boldsymbol{y}}$
+and $\tilde{X}_{ij} = X_{ij} - \frac{1}{n}\sum_{k=0}^{n-1}X_{kj}$.
+
+For Ridge regression we need to add $\lambda \boldsymbol{\beta}^T\boldsymbol{\beta}$ to the cost function and get then
+
+$$
+\hat{\boldsymbol{\beta}} = (\tilde{X}^T\tilde{X} + \lambda I)^{-1}\tilde{X}^T\boldsymbol{\tilde{y}}.
+$$
+
+What does this mean? And why do we insist on all this? Let us look at some examples.
+
+
+This code shows a simple first-order fit to a data set using the above transformed data, where we consider the role of the intercept first, by either excluding it or including it (*code example thanks to  Øyvind Sigmundson Schøyen*). Here our scaling of the data is done by subtracting the mean values only.
+Note also that we do not split the data into training and test.
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from sklearn.linear_model import LinearRegression
+
+
+np.random.seed(2021)
+
+def MSE(y_data,y_model):
+    n = np.size(y_model)
+    return np.sum((y_data-y_model)**2)/n
+
+
+def fit_beta(X, y):
+    return np.linalg.pinv(X.T @ X) @ X.T @ y
+
+
+true_beta = [2, 0.5, 3.7]
+
+x = np.linspace(0, 1, 11)
+y = np.sum(
+    np.asarray([x ** p * b for p, b in enumerate(true_beta)]), axis=0
+) + 0.1 * np.random.normal(size=len(x))
+
+degree = 3
+X = np.zeros((len(x), degree))
+
+# Include the intercept in the design matrix
+for p in range(degree):
+    X[:, p] = x ** p
+
+beta = fit_beta(X, y)
+
+# Intercept is included in the design matrix
+skl = LinearRegression(fit_intercept=False).fit(X, y)
+
+print(f"True beta: {true_beta}")
+print(f"Fitted beta: {beta}")
+print(f"Sklearn fitted beta: {skl.coef_}")
+ypredictOwn = X @ beta
+ypredictSKL = skl.predict(X)
+print(f"MSE with intercept column")
+print(MSE(y,ypredictOwn))
+print(f"MSE with intercept column from SKL")
+print(MSE(y,ypredictSKL))
+
+
+plt.figure()
+plt.scatter(x, y, label="Data")
+plt.plot(x, X @ beta, label="Fit")
+plt.plot(x, skl.predict(X), label="Sklearn (fit_intercept=False)")
+
+
+# Do not include the intercept in the design matrix
+X = np.zeros((len(x), degree - 1))
+
+for p in range(degree - 1):
+    X[:, p] = x ** (p + 1)
+
+# Intercept is not included in the design matrix
+skl = LinearRegression(fit_intercept=True).fit(X, y)
+
+# Use centered values for X and y when computing coefficients
+y_offset = np.average(y, axis=0)
+X_offset = np.average(X, axis=0)
+
+beta = fit_beta(X - X_offset, y - y_offset)
+intercept = np.mean(y_offset - X_offset @ beta)
+
+print(f"Manual intercept: {intercept}")
+print(f"Fitted beta (wiothout intercept): {beta}")
+print(f"Sklearn intercept: {skl.intercept_}")
+print(f"Sklearn fitted beta (without intercept): {skl.coef_}")
+ypredictOwn = X @ beta
+ypredictSKL = skl.predict(X)
+print(f"MSE with Manual intercept")
+print(MSE(y,ypredictOwn+intercept))
+print(f"MSE with Sklearn intercept")
+print(MSE(y,ypredictSKL))
+
+plt.plot(x, X @ beta + intercept, "--", label="Fit (manual intercept)")
+plt.plot(x, skl.predict(X), "--", label="Sklearn (fit_intercept=True)")
+plt.grid()
+plt.legend()
+
+plt.show()
+
+The intercept is the value of our output/target variable
+when all our features are zero and our function crosses the $y$-axis (for a one-dimensional case). 
+
+Printing the MSE, we see first that both methods give the same MSE, as
+they should.  However, when we move to for example Ridge regression,
+the way we treat the intercept may give a larger or smaller MSE,
+meaning that the MSE can be penalized by the value of the
+intercept. Not including the intercept in the fit, means that the
+regularization term does not include $\beta_0$. For different values
+of $\lambda$, this may lead to differeing MSE values. 
+
+To remind the reader, the regularization term, with the intercept in Ridge regression, is given by
+
+$$
+\lambda \vert\vert \boldsymbol{\beta} \vert\vert_2^2 = \lambda \sum_{j=0}^{p-1}\beta_j^2,
+$$
+
+but when we take out the intercept, this equation becomes
+
+$$
+\lambda \vert\vert \boldsymbol{\beta} \vert\vert_2^2 = \lambda \sum_{j=1}^{p-1}\beta_j^2.
+$$
+
+For Lasso regression we have
+
+$$
+\lambda \vert\vert \boldsymbol{\beta} \vert\vert_1 = \lambda \sum_{j=1}^{p-1}\vert\beta_j\vert.
+$$
+
+It means that, when scaling the design matrix and the outputs/targets,
+by subtracting the mean values, we have an optimization problem which
+is not penalized by the intercept. The MSE value can then be smaller
+since it focuses only on the remaining quantities. If we however bring
+back the intercept, we will get a MSE which then contains the
+intercept.
+
+
+Armed with this wisdom, we attempt first to simply set the intercept equal to **False** in our implementation of Ridge regression for our well-known  vanilla data set.
+
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from sklearn.model_selection import train_test_split
+from sklearn import linear_model
+
+def MSE(y_data,y_model):
+    n = np.size(y_model)
+    return np.sum((y_data-y_model)**2)/n
+
+
+# A seed just to ensure that the random numbers are the same for every run.
+# Useful for eventual debugging.
+np.random.seed(3155)
+
+n = 100
+x = np.random.rand(n)
+y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
+
+Maxpolydegree = 20
+X = np.zeros((n,Maxpolydegree))
+#We include explicitely the intercept column
+for degree in range(Maxpolydegree):
+    X[:,degree] = x**degree
+# We split the data in test and training data
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
+
+p = Maxpolydegree
+I = np.eye(p,p)
+# Decide which values of lambda to use
+nlambdas = 6
+MSEOwnRidgePredict = np.zeros(nlambdas)
+MSERidgePredict = np.zeros(nlambdas)
+lambdas = np.logspace(-4, 2, nlambdas)
+for i in range(nlambdas):
+    lmb = lambdas[i]
+    OwnRidgeBeta = np.linalg.pinv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train
+    # Note: we include the intercept column and no scaling
+    RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
+    RegRidge.fit(X_train,y_train)
+    # and then make the prediction
+    ytildeOwnRidge = X_train @ OwnRidgeBeta
+    ypredictOwnRidge = X_test @ OwnRidgeBeta
+    ytildeRidge = RegRidge.predict(X_train)
+    ypredictRidge = RegRidge.predict(X_test)
+    MSEOwnRidgePredict[i] = MSE(y_test,ypredictOwnRidge)
+    MSERidgePredict[i] = MSE(y_test,ypredictRidge)
+    print("Beta values for own Ridge implementation")
+    print(OwnRidgeBeta)
+    print("Beta values for Scikit-Learn Ridge implementation")
+    print(RegRidge.coef_)
+    print("MSE values for own Ridge implementation")
+    print(MSEOwnRidgePredict[i])
+    print("MSE values for Scikit-Learn Ridge implementation")
+    print(MSERidgePredict[i])
+
+# Now plot the results
+plt.figure()
+plt.plot(np.log10(lambdas), MSEOwnRidgePredict, 'r', label = 'MSE own Ridge Test')
+plt.plot(np.log10(lambdas), MSERidgePredict, 'g', label = 'MSE Ridge Test')
+
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+
+The results here agree when we force **Scikit-Learn**'s Ridge function to include the first column in our design matrix.
+We see that the results agree very well. Here we have thus explicitely included the intercept column in the design matrix.
+What happens if we do not include the intercept in our fit?
+Let us see how we can change this code by zero centering.
+
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from sklearn.model_selection import train_test_split
+from sklearn import linear_model
+from sklearn.preprocessing import StandardScaler
+
+def MSE(y_data,y_model):
+    n = np.size(y_model)
+    return np.sum((y_data-y_model)**2)/n
+# A seed just to ensure that the random numbers are the same for every run.
+# Useful for eventual debugging.
+np.random.seed(315)
+
+n = 100
+x = np.random.rand(n)
+y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
+
+Maxpolydegree = 20
+X = np.zeros((n,Maxpolydegree-1))
+
+for degree in range(1,Maxpolydegree): #No intercept column
+    X[:,degree-1] = x**(degree)
+
+# We split the data in test and training data
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
+
+#For our own implementation, we will need to deal with the intercept by centering the design matrix and the target variable
+X_train_mean = np.mean(X_train,axis=0)
+#Center by removing mean from each feature
+X_train_scaled = X_train - X_train_mean 
+X_test_scaled = X_test - X_train_mean
+#The model intercept (called y_scaler) is given by the mean of the target variable (IF X is centered)
+#Remove the intercept from the training data.
+y_scaler = np.mean(y_train)           
+y_train_scaled = y_train - y_scaler   
+
+p = Maxpolydegree-1
+I = np.eye(p,p)
+# Decide which values of lambda to use
+nlambdas = 6
+MSEOwnRidgePredict = np.zeros(nlambdas)
+MSERidgePredict = np.zeros(nlambdas)
+
+lambdas = np.logspace(-4, 2, nlambdas)
+for i in range(nlambdas):
+    lmb = lambdas[i]
+    OwnRidgeBeta = np.linalg.pinv(X_train_scaled.T @ X_train_scaled+lmb*I) @ X_train_scaled.T @ (y_train_scaled)
+    intercept_ = y_scaler - X_train_mean@OwnRidgeBeta #The intercept can be shifted so the model can predict on uncentered data
+    #Add intercept to prediction
+    ypredictOwnRidge = X_test @ OwnRidgeBeta + intercept_ 
+    #Add intercept to prediction
+    ypredictOwnRidge = X_test_scaled @ OwnRidgeBeta + y_scaler 
+    RegRidge = linear_model.Ridge(lmb)
+    RegRidge.fit(X_train,y_train)
+    ypredictRidge = RegRidge.predict(X_test)
+    MSEOwnRidgePredict[i] = MSE(y_test,ypredictOwnRidge)
+    MSERidgePredict[i] = MSE(y_test,ypredictRidge)
+    print("Beta values for own Ridge implementation")
+    print(OwnRidgeBeta) #Intercept is given by mean of target variable
+    print("Beta values for Scikit-Learn Ridge implementation")
+    print(RegRidge.coef_)
+    print('Intercept from own implementation:')
+    print(intercept_)
+    print('Intercept from Scikit-Learn Ridge implementation')
+    print(RegRidge.intercept_)
+    print("MSE values for own Ridge implementation")
+    print(MSEOwnRidgePredict[i])
+    print("MSE values for Scikit-Learn Ridge implementation")
+    print(MSERidgePredict[i])
+
+
+# Now plot the results
+plt.figure()
+plt.plot(np.log10(lambdas), MSEOwnRidgePredict, 'b--', label = 'MSE own Ridge Test')
+plt.plot(np.log10(lambdas), MSERidgePredict, 'g--', label = 'MSE SL Ridge Test')
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+
+We see here, when compared to the code which includes explicitely the
+intercept column, that our MSE value is actually smaller. This is
+because the regularization term does not include the intercept value
+$\beta_0$ in the fitting.  This applies to Lasso regularization as
+well.  It means that our optimization is now done only with the
+centered matrix and/or vector that enter the fitting procedure. Note
+also that the problem with the intercept occurs mainly in these type
+of polynomial fitting problem.
+
+The next example is indeed an example where all these discussions about the role of intercept are not present.
+
+## More complicated Example: The Ising model
+
+The one-dimensional Ising model with nearest neighbor interaction, no
+external field and a constant coupling constant $J$ is given by
+
+<!-- Equation labels as ordinary links -->
+<div id="_auto1"></div>
+
+$$
+\begin{equation}
+    H = -J \sum_{k}^L s_k s_{k + 1},
+\label{_auto1} \tag{1}
+\end{equation}
+$$
+
+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 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))
+
+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.
+
+A more general form for the one-dimensional Ising model is
+
+<!-- Equation labels as ordinary links -->
+<div id="_auto2"></div>
+
+$$
+\begin{equation}
+    H = - \sum_j^L \sum_k^L s_j s_k J_{jk}.
+\label{_auto2} \tag{2}
+\end{equation}
+$$
+
+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
+
+<!-- Equation labels as ordinary links -->
+<div id="_auto3"></div>
+
+$$
+\begin{equation}
+    \boldsymbol{H} = \boldsymbol{X} J,
+\label{_auto3} \tag{3}
+\end{equation}
+$$
+
+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
+
+<!-- Equation labels as ordinary links -->
+<div id="_auto4"></div>
+
+$$
+\begin{equation}
+    \boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\epsilon},
+\label{_auto4} \tag{4}
+\end{equation}
+$$
+
+We split the data in training and test data as discussed in the previous example
+
+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)
+
+In the ordinary least squares method we choose the cost function
+
+<!-- Equation labels as ordinary links -->
+<div id="_auto5"></div>
+
+$$
+\begin{equation}
+    C(\boldsymbol{X}, \boldsymbol{\beta})= \frac{1}{n}\left\{(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})^T(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})\right\}.
+\label{_auto5} \tag{5}
+\end{equation}
+$$
+
+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
+)
+
+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
+
+<!-- Equation labels as ordinary links -->
+<div id="_auto6"></div>
+
+$$
+\begin{equation}
+    \boldsymbol{\beta} = \boldsymbol{V}\boldsymbol{\Sigma}^{+} \boldsymbol{U}^T \boldsymbol{y}.
+\label{_auto6} \tag{6}
+\end{equation}
+$$
+
+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?
+
+
+
+
+
+Let us now 
+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
+
+<!-- Equation labels as ordinary links -->
+<div id="_auto7"></div>
+
+$$
+\begin{equation}
+    H = -J \sum_{k}^L s_k s_{k + 1},
+\label{_auto7} \tag{7}
+\end{equation}
+$$
+
+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
+
+<!-- Equation labels as ordinary links -->
+<div id="_auto8"></div>
+
+$$
+\begin{equation}
+    H = - \sum_j^L \sum_k^L s_j s_k J_{jk}.
+\label{_auto8} \tag{8}
+\end{equation}
+$$
+
+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
+
+<!-- Equation labels as ordinary links -->
+<div id="_auto9"></div>
+
+$$
+\begin{equation}
+    H = X J,
+\label{_auto9} \tag{9}
+\end{equation}
+$$
+
+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.
+
+<!-- Equation labels as ordinary links -->
+<div id="_auto10"></div>
+
+$$
+\begin{equation}
+    \boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}.
+\label{_auto10} \tag{10}
+\end{equation}
+$$
+
+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 agree perfectly with our previous discussion where we used our own code.
+
+
+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
+
+6
+0
+ 
+<
+<
+<
+!
+!
+M
+A
+T
+H
+_
+B
+L
+O
+C
+K
+
+_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()
+
+In the **Least Absolute Shrinkage and Selection Operator** (LASSO)-method we get a third cost function.
+
+<!-- Equation labels as ordinary links -->
+<div id="_auto12"></div>
+
+$$
+\begin{equation}
+    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}}.
+\label{_auto12} \tag{12}
+\end{equation}
+$$
+
+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$.
+
+
+
+
+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.
+
+
+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$.
+
+
+
+
+
+
 ## Exercises and Projects
 
 
@@ -1439,6 +2402,7 @@ scipy.misc.imread
 Here is a simple part of a Python code which reads and plots the data
 from such files
 
+"""
 import numpy as np
 from imageio import imread
 import matplotlib.pyplot as plt
@@ -1454,6 +2418,7 @@ plt.imshow(terrain1, cmap='gray')
 plt.xlabel('X')
 plt.ylabel('Y')
 plt.show()
+"""
 
 If you should have problems in downloading the digital terrain data,
 we provide two examples under the data folder of project 1. One is