The main topics are:
Brief repetition from last week
Discussions of the equations for ordinary least squares
Discussions of the equations for ordinary least squares (OLS)
Discussion on how to prepare data and examples of applications of linear regression
Material for the lecture on Monday: Mathematical interpretations of linear regression
Monday: Ridge and Lasso regression and Singular Value Decomposition
Mathematical interpretations of OLS
Introduction of Ridge and Lasso regression
Goodfellow, Bengio and Courville, Deep Learning, chapter 2 on linear algebra and sections 3.1-3.10 on elements of statistics (background)
Goodfellow, Bengio and Courville, Deep Learning, chapter 2 on linear algebra
Raschka et al on preprocessing of data, relevant for exercise 3 this week, see chapter 4.
For exercise 1 of week 35, the book by A. Aldo Faisal, Cheng Soon Ong, and Marc Peter Deisenroth on the Mathematics of Machine Learning, may be very relevant. In particular chapter 5 at URL”https://mml-book.github.io/” (section 5.5 on derivatives) is very useful for exercise 1 this coming week.
where \(\langle y_i \rangle\) is the mean value. Keep in mind also that till now we have treated \(y_i\) as the exact value. Normally, the -response (dependent or outcome) variable \(y_i\) is the outcome of a +output (response, target, dependent or outcome) variable \(y_i\) is the outcome of a numerical experiment or another type of experiment and could thus be treated itself as an approximation to the true value. It is then always accompanied by an error estimate, often limited to a statistical error estimate given by the standard deviation discussed earlier. In the discussion here we -will treat \(y_i\) as our exact value for the response variable.
+will treat \(y_i\) as our exact value for the output variable.In order to find the parameters \(\theta_i\) we will then minimize the spread of \(C(\boldsymbol{\theta})\), that is we are going to solve the problem
In practical terms it means we will require
@@ -611,7 +610,7 @@ will treat \(y_i\) as our exacand if the matrix \(\boldsymbol{X}^T\boldsymbol{X}\) is invertible we have the solution
We note also that since our design matrix is defined as \(\boldsymbol{X}\in {\mathbb{R}}^{n\times p}\), the product \(\boldsymbol{X}^T\boldsymbol{X} \in @@ -621,8 +620,8 @@ matrices to invert. The methods discussed here and for many other supervised learning algorithms like classification with logistic regression or support vector machines, exhibit dimensionalities which allow for the usage of direct linear algebra methods such as LU decomposition or Singular Value Decomposition (SVD) for finding the inverse of the matrix -\(\boldsymbol{X}^T\boldsymbol{X}\). This is discussed on Thursday this week.
-Small question: Do you think the example we have at hand here (the nuclear binding energies) can lead to problems in inverting the matrix \(\boldsymbol{X}^T\boldsymbol{X}\)? What kind of problems can we expect?
+\(\boldsymbol{X}^T\boldsymbol{X}\). +Small question: When inverting the matrix $\boldsymbol{X}^T\boldsymbol{X}, what kind of problems can we expect?
We list here some other useful relations we may encounter (recall that vectors are defined by boldfaced low-key letters)
-A very important matrix we will meet again and again in machine @@ -859,13 +838,12 @@ second derivative of the mean squared error as cost function is,
\]For ordinary least squares, it is inversely proportional (derivation next week) with the variance of the optimal parameters -\(\hat{\boldsymbol{\theta}}\). Furthermore, we will see later this week that it is +\(\hat{\boldsymbol{\theta}}\). Furthermore, we will see next week that it is (aside the factor \(1/n\)) equal to the covariance matrix. It plays also a very important role in optmization algorithms and Principal Component Analysis as a way to reduce the dimensionality of a machine learning/data analysis -problem. -v -Linear algebra question: Can we use the Hessian matrix to say something about properties of the cost function (our optmization problem)? (hint: think about convex or concave problems and how to relate these to a matrix!).
+problem. We will discuss this in greater detail next week when we introduce gradient methods. +Linear algebra question: Can we use the Hessian matrix to say something about properties of the cost function (our optmization problem)? (hint: think about convex or concave problems and how to relate these to a matrix!).
# matrix inversion to find beta
+# matrix inversion to find theta
# First we set up the data
import numpy as np
x = np.random.rand(100)
@@ -916,9 +894,9 @@ X[:,1] = x
X[:,2] = x**2
X[:,3] = x**3
X[:,4] = x**4
-beta = (np.linalg.inv(X.T @ X) @ X.T ) @ y
+theta = (np.linalg.inv(X.T @ X) @ X.T ) @ y
# and then make the prediction
-ytilde = X @ beta
+ytilde = X @ theta
As we will also see in the first project, this may however not the be case in general and a standard matrix inversion -algorithm based on say LU, QR or Cholesky decomposition may lead to singularities. We will see examples of this below.
+algorithm based on say LU, QR or Cholesky decomposition may lead to singularities. We will see examples of this below and in other examples.There is however a way to circumvent this problem and also gain some insights about the ordinary least squares approach, and later shrinkage methods like Ridge and Lasso regressions.
@@ -1404,9 +1382,175 @@ the regression parameters \(\theta_i\ \]where \(\boldsymbol{I}\) is the identity matrix. When we discuss Ridge regression this is actually what we end up evaluating. The parameter \(\lambda\) is called a hyperparameter. More about this later.
+Let us remind ourselves about 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
+or we can state it as
+where we have used the definition of a norm-2 vector, that is
+By minimizing the above equation with respect to the parameters +\(\boldsymbol{\theta}\) we could then obtain an analytical expression for the +parameters \(\boldsymbol{\theta}\). We can add a regularization parameter \(\lambda\) by +defining a new cost function to be optimized, that is
+which leads to the Ridge regression minimization problem where we +require that \(\vert\vert \boldsymbol{\theta}\vert\vert_2^2\le t\), where \(t\) is +a finite number larger than zero. By defining
+we have a new optimization equation
+which leads to Lasso regression. Lasso stands for least absolute shrinkage and selection operator.
+Here we have defined the norm-1 as
+Using the matrix-vector expression for Ridge regression and dropping the parameter \(1/n\) in front of the standard means squared error equation, we have
+and +taking the derivatives with respect to \(\boldsymbol{\theta}\) we obtain then +a slightly modified matrix inversion problem which for finite values +of \(\lambda\) does not suffer from singularity problems. We obtain +the optimal parameters
+with \(\boldsymbol{I}\) being a \(p\times p\) identity matrix with the constraint that
+with \(t\) a finite positive number.
+If we keep the \(1/n\) factor, the equation for the optimal \(\theta\) changes to
+In many textbooks the \(1/n\) term is often omitted. Note that a library like Scikit-Learn does not include the \(1/n\) factor in the setup of the cost function.
+When we compare this with the ordinary least squares result we have
+which can lead to singular matrices. However, with the SVD, we can always compute the inverse of the matrix \(\boldsymbol{X}^T\boldsymbol{X}\).
+We see that Ridge regression is nothing but the standard OLS with a +modified diagonal term added to \(\boldsymbol{X}^T\boldsymbol{X}\). The consequences, in +particular for our discussion of the bias-variance tradeoff are rather +interesting. We will see that for specific values of \(\lambda\), we may +even reduce the variance of the optimal parameters \(\boldsymbol{\theta}\). These topics and other related ones, will be discussed after the more linear algebra oriented analysis here.
+When we have discussed the singular value decomposition of the design +matrix \(\boldsymbol{X}\), we will in turn perform a more rigorous mathematical +discussion of Ridge regression.
+The code here is a simple demonstration of how to implement Ridge regression with our own code and compare this with scikit-learn.
+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]
+ OwnRidgeTheta = 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 @ OwnRidgeTheta
+ ypredictOwnRidge = X_test @ OwnRidgeTheta
+ ytildeRidge = RegRidge.predict(X_train)
+ ypredictRidge = RegRidge.predict(X_test)
+ MSEOwnRidgePredict[i] = MSE(y_test,ypredictOwnRidge)
+ MSERidgePredict[i] = MSE(y_test,ypredictRidge)
+ print("Theta values for own Ridge implementation")
+ print(OwnRidgeTheta)
+ print("Theta 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? We will discuss this in more detail next week.
+From standard linear algebra we know that a square matrix \(\boldsymbol{X}\) can be diagonalized if and only it is +
From standard linear algebra we know that a square matrix \(\boldsymbol{X}\) can be diagonalized if and only if it is a so-called normal matrix, that is if \(\boldsymbol{X}\in {\mathbb{R}}^{n\times n}\) we have \(\boldsymbol{X}\boldsymbol{X}^T=\boldsymbol{X}^T\boldsymbol{X}\) or if \(\boldsymbol{X}\in {\mathbb{C}}^{n\times n}\) we have \(\boldsymbol{X}\boldsymbol{X}^{\dagger}=\boldsymbol{X}^{\dagger}\boldsymbol{X}\). The matrix has then a set of eigenpairs
@@ -1946,66 +2090,6 @@ print(correlation_matrix)We expand this model to the Franke function discussed above.
-# Common imports
-import numpy as np
-import pandas as pd
-
-
-def FrankeFunction(x,y):
- term1 = 0.75*np.exp(-(0.25*(9*x-2)**2) - 0.25*((9*y-2)**2))
- term2 = 0.75*np.exp(-((9*x+1)**2)/49.0 - 0.1*(9*y+1))
- term3 = 0.5*np.exp(-(9*x-7)**2/4.0 - 0.25*((9*y-3)**2))
- term4 = -0.2*np.exp(-(9*x-4)**2 - (9*y-7)**2)
- return term1 + term2 + term3 + term4
-
-
-def create_X(x, y, n ):
- if len(x.shape) > 1:
- x = np.ravel(x)
- y = np.ravel(y)
-
- N = len(x)
- l = int((n+1)*(n+2)/2) # Number of elements in beta
- X = np.ones((N,l))
-
- for i in range(1,n+1):
- q = int((i)*(i+1)/2)
- for k in range(i+1):
- X[:,q+k] = (x**(i-k))*(y**k)
-
- return X
-
-
-# Making meshgrid of datapoints and compute Franke's function
-n = 4
-N = 100
-x = np.sort(np.random.uniform(0, 1, N))
-y = np.sort(np.random.uniform(0, 1, N))
-z = FrankeFunction(x, y)
-X = create_X(x, y, n=n)
-
-Xpd = pd.DataFrame(X)
-# subtract the mean values and set up the covariance matrix
-Xpd = Xpd - Xpd.mean()
-covariance_matrix = Xpd.cov()
-print(covariance_matrix)
-We note here that the covariance is zero for the first rows and -columns since all matrix elements in the design matrix were set to one -(we are fitting the function in terms of a polynomial of degree \(n\)).
-This means that the variance for these elements will be zero and will -cause problems when we set up the correlation matrix. We can simply -drop these elements and construct a correlation -matrix without these elements.
Let us remind ourselves about 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
Using the matrix-vector expression for Ridge regression and dropping the parameter \(1/n\) in front of the standard means squared error equation, we have
-and -taking the derivatives with respect to \(\boldsymbol{\theta}\) we obtain then -a slightly modified matrix inversion problem which for finite values -of \(\lambda\) does not suffer from singularity problems. We obtain -the optimal parameters
-with \(\boldsymbol{I}\) being a \(p\times p\) identity matrix with the constraint that
-with \(t\) a finite positive number.
-If we keep the \(1/n\) factor, the equation for the optimal \(\theta\) changes to
-In many textbooks the \(1/n\) term is often omitted. Note that a library like Scikit-Learn does not include the \(1/n\) factor in the setup of the cost function.
-When we compare this with the ordinary least squares result we have
-which can lead to singular matrices. However, with the SVD, we can always compute the inverse of the matrix \(\boldsymbol{X}^T\boldsymbol{X}\).
-We see that Ridge regression is nothing but the standard OLS with a +
Ridge regression, as discussed above, is nothing but the standard OLS with a modified diagonal term added to \(\boldsymbol{X}^T\boldsymbol{X}\). The consequences, in particular for our discussion of the bias-variance tradeoff are rather interesting. We will see that for specific values of \(\lambda\), we may @@ -2301,7 +2350,7 @@ C(\boldsymbol{X},\boldsymbol{\theta})=\frac{1}{n}\left\{(\boldsymbol{y}-\boldsym \[ \boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\theta}+\lambda sgn(\boldsymbol{\theta})=2\boldsymbol{X}^T\boldsymbol{y}. \] -
This equation does not lead to a nice analytical equation as in either Ridge regression or ordinary least squares. This equation can however be solved by using standard convex optimization algorithms using for example the Python package CVXOPT. We will discuss this later.
+This equation does not lead to a nice analytical equation as in either Ridge regression or ordinary least squares. This equation can however be solved by using standard convex optimization algorithms using for example the Python package CVXOPT. We will discuss how to code LASSO regression next week, when we have introduced gradient methods.