642 KiB
Week 35: From Ordinary Linear Regression to Ridge and Lasso Regression
Morten Hjorth-Jensen, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University
Date: August 28-September 1
Plans for week 35
The main topics are:
-
Brief repetition from last week
-
Derivation of the equations for ordinary least squares
-
Discussion on how to prepare data and examples of applications of linear regression
-
Material for the lecture on Thursday: Mathematical interpretations of linear regression
-
Thursday: Ridge and Lasso regression and Singular Value Decomposition
Reading recommendations:
-
See lecture notes for week 35 at https://compphysics.github.io/MachineLearning/doc/web/course.html
-
Goodfellow, Bengio and Courville, Deep Learning, chapter 2 on linear algebra and sections 3.1-3.10 on elements of statistics (background)
-
Hastie, Tibshirani and Friedman, The elements of statistical learning, sections 3.1-3.4 (on relevance for the discussion of linear regression).
Why Linear Regression (aka Ordinary Least Squares and family), repeat from last week
We need first a reminder from last week about linear regression.
Fitting a continuous function with linear parameterization in terms of the parameters \boldsymbol{\beta}.
-
Method of choice for fitting a continuous function!
-
Gives an excellent introduction to central Machine Learning features with understandable pedagogical links to other methods like Neural Networks, Support Vector Machines etc
-
Analytical expression for the fitting parameters
\boldsymbol{\beta} -
Analytical expressions for statistical propertiers like mean values, variances, confidence intervals and more
-
Analytical relation with probabilistic interpretations
-
Easy to introduce basic concepts like bias-variance tradeoff, cross-validation, resampling and regularization techniques and many other ML topics
-
Easy to code! And links well with classification problems and logistic regression and neural networks
-
Allows for easy hands-on understanding of gradient descent methods
-
and many more features
For more discussions of Ridge and Lasso regression, Wessel van Wieringen's article is highly recommended. Similarly, Mehta et al's article is also recommended.
The equations for ordinary least squares
Our data which we want to apply a machine learning method on, consist
of a set of inputs \boldsymbol{x}^T=[x_0,x_1,x_2,\dots,x_{n-1}] and the
outputs we want to model \boldsymbol{y}^T=[y_0,y_1,y_2,\dots,y_{n-1}].
We assume that the output data can be represented (for a regression case) by a continuous function f
through
y_i=f(x_i)+\epsilon_i,
or in general
\boldsymbol{y}=f(\boldsymbol{x})+\boldsymbol{\epsilon},
where \boldsymbol{\epsilon} represents some noise which is normally assumed to
be distributed via a normal probability distribution with zero mean
value and a variance \sigma^2.
In linear regression we approximate the unknown function with another
continuous function \tilde{\boldsymbol{y}}(\boldsymbol{x}) which depends linearly on
some unknown parameters
\boldsymbol{\beta}^T=[\beta_0,\beta_1,\beta_2,\dots,\beta_{p-1}].
Last week we introduced the so-called design matrix in order to define
the approximation \boldsymbol{\tilde{y}} via the unknown quantity
\boldsymbol{\beta} as
\boldsymbol{\tilde{y}}= \boldsymbol{X}\boldsymbol{\beta},
and in order to find the optimal parameters \beta_i we defined a function which
gives a measure of the spread between the values y_i (which
represent the output values we want to reproduce) and the parametrized
values \tilde{y}_i, namely the so-called cost/loss function.
The cost/loss function
We used the mean squared error to define the way we measure the quality of our model
C(\boldsymbol{\beta})=\frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)^T\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)\right\},
or using the matrix \boldsymbol{X} and in a more compact matrix-vector notation as
C(\boldsymbol{\beta})=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right\}.
This function represents one of many possible ways to define the so-called cost function.
It is also common to define
the function C as
C(\boldsymbol{\beta})=\frac{1}{2n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2,
since when taking the first derivative with respect to the unknown parameters \beta, the factor of 2 cancels out.
Interpretations and optimizing our parameters
The function
C(\boldsymbol{\beta})=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right\},
can be linked to the variance of the quantity y_i if we interpret the latter as the mean value.
When linking (see the discussions next week) with the maximum likelihood approach below, we will indeed interpret y_i as a mean value
y_{i}=\langle y_i \rangle = \beta_0x_{i,0}+\beta_1x_{i,1}+\beta_2x_{i,2}+\dots+\beta_{n-1}x_{i,n-1}+\epsilon_i,
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
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.
In order to find the parameters \beta_i we will then minimize the spread of C(\boldsymbol{\beta}), that is we are going to solve the problem
{\displaystyle \min_{\boldsymbol{\beta}\in
{\mathbb{R}}^{p}}}\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right\}.
In practical terms it means we will require
\frac{\partial C(\boldsymbol{\beta})}{\partial \beta_j} = \frac{\partial }{\partial \beta_j}\left[ \frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}\right)^2\right]=0,
which results in
\frac{\partial C(\boldsymbol{\beta})}{\partial \beta_j} = -\frac{2}{n}\left[ \sum_{i=0}^{n-1}x_{ij}\left(y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}\right)\right]=0,
or in a matrix-vector form as (multiplying away the factor -2/n, see derivation below)
\frac{\partial C(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}^T} = 0 = \boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right).
Interpretations and optimizing our parameters
We can rewrite, see the derivations below,
\frac{\partial C(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}^T} = 0 = \boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right),
as
\boldsymbol{X}^T\boldsymbol{y} = \boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\beta},
and if the matrix \boldsymbol{X}^T\boldsymbol{X} is invertible we have the solution
\boldsymbol{\beta} =\left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}.
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
{\mathbb{R}}^{p\times p}$. In most cases we have that p \ll n. In our example case below we have p=5 meaning. We end up with inverting a small
5\times 5 matrix. This is a rather common situation, in many cases we end up with low-dimensional
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?
Some useful matrix and vector expressions
The following matrix and vector relation will be useful here and for the rest of the course. Vectors are always written as boldfaced lower case letters and matrices as upper case boldfaced letters. In the following we will discuss how to calculate derivatives of various matrices relevant for machine learning. We will often represent our data in terms of matrices and vectors.
Let us introduce first some conventions. We assume that \boldsymbol{y} is a
vector of length m, that is it has m elements $y_0,y_1,\dots,
y_{m-1}$. By convention we start labeling vectors with the zeroth
element, as are arrays in Python and C++/C, for example. Similarly, we
have a vector \boldsymbol{x} of length n, that is
\boldsymbol{x}^T=[x_0,x_1,\dots, x_{n-1}].
We assume also that \boldsymbol{y} is a function of \boldsymbol{x} through some
given function f
\boldsymbol{y}=f(\boldsymbol{x}).
The Jacobian
We define the partial derivatives of the various components of \boldsymbol{y} as functions of x_i in terms of the so-called Jacobian matrix
\boldsymbol{J}=\frac{\partial \boldsymbol{y}}{\partial \boldsymbol{x}}=\begin{bmatrix} \frac{\partial y_0}{\partial x_0} & \frac{\partial y_0}{\partial x_1} & \frac{\partial y_0}{\partial x_2} & \dots & \dots & \frac{\partial y_0}{\partial x_{n-1}} \\ \frac{\partial y_1}{\partial x_0} & \frac{\partial y_1}{\partial x_1} & \frac{\partial y_1}{\partial x_2} & \dots & \dots & \frac{\partial y_1}{\partial x_{n-1}} \\
\frac{\partial y_2}{\partial x_0} & \frac{\partial y_2}{\partial x_1} & \frac{\partial y_2}{\partial x_2} & \dots & \dots & \frac{\partial y_2}{\partial x_{n-1}} \\
\dots & \dots & \dots & \dots & \dots & \dots \\
\dots & \dots & \dots & \dots & \dots & \dots \\
\frac{\partial y_{m-1}}{\partial x_0} & \frac{\partial y_{m-1}}{\partial x_1} & \frac{\partial y_{m-1}}{\partial x_2} & \dots & \dots & \frac{\partial y_{m-1}}{\partial x_{n-1}} \end{bmatrix},
which is an m\times n matrix. If \boldsymbol{x} is a scalar, then the
Jacobian is only a single-column vector, or an m\times 1 matrix. If
on the other hand \boldsymbol{y} is a scalar, the Jacobian becomes a
1\times n matrix.
When this matrix is a square matrix m=n, its determinant is often referred to as the Jacobian
determinant. Both the matrix and (if m=n) the determinant are
often referred to simply as the Jacobian. The Jacobian matrix represents the differential of \boldsymbol{y} at every point where the
vector is differentiable.
Derivatives, example 1
Let now \boldsymbol{y}=\boldsymbol{A}\boldsymbol{x}, where \boldsymbol{A} is an m\times n matrix and the matrix does not depend on \boldsymbol{x}. If we write out the vector \boldsymbol{y} compoment by component we have
y_i = \sum_{j=0}^{n-1}a_{ij}x_j,
with \forall i=0,1,2,\dots,m-1. The individual matrix elements of \boldsymbol{A} are given by the symbol a_{ij}.
It follows that the partial derivatives of y_i with respect to x_k
\frac{\partial y_i }{\partial x_k}= a_{ik} \forall i=0,1,2,\dots,m-1.
From this we have, using the definition of the Jacobian
\frac{\partial \boldsymbol{y} }{\partial \boldsymbol{x}}= \boldsymbol{A}.
Example 2
We define a scalar (our cost/loss functions are in general also scalars, just think of the mean squared error) as the result of some matrix vector multiplications
\alpha = \boldsymbol{y}^T\boldsymbol{A}\boldsymbol{x},
with \boldsymbol{y} a vector of length m, \boldsymbol{A} an m\times n matrix and \boldsymbol{x} a vector of length n. We assume also that \boldsymbol{A} does not depend on any of the two vectors.
In order to find the derivative of \alpha with respect to the two vectors, we define an intermediate vector \boldsymbol{z}. We define first
\boldsymbol{z}^T=\boldsymbol{y}^T\boldsymbol{A}, a vector of length n. We have then, using the definition of the Jacobian,
\alpha = \boldsymbol{z}^T\boldsymbol{x},
which means that (using our previous example) we have
\frac{\partial \alpha}{\partial \boldsymbol{x}} = \boldsymbol{z}=bm{A}^T\boldsymbol{y}.
Note that the resulting vector elements are the same for \boldsymbol{z}^T and \boldsymbol{z}, the only difference is that one if just the transpose of the other.
Since \alpha is a scalar we have \alpha =\alpha^T=\boldsymbol{x}^T\boldsymbol{A}^T\boldsymbol{y}. Defining now \boldsymbol{z}=\boldsymbol{x}^T\boldsymbol{A}^T we find that
\frac{\partial \alpha}{\partial \boldsymbol{y}} = \boldsymbol{z}^T=\boldsymbol{x}^T\boldsymbol{A}^T.
Example 3
We start with a new scalar but where now the vector \boldsymbol{y} is
replaced by a vector \boldsymbol{x} and the matrix \boldsymbol{A} is a square
matrix with dimension n\times n.
\alpha = \boldsymbol{x}^T\boldsymbol{A}\boldsymbol{x},
with \boldsymbol{x} a vector of length n.
We write out the specific sums involved in the calculation of \alpha
\alpha = \sum_{i=0}^{n-1}\sum_{j=0}^{n-1}x_i a_{ij}x_j,
taking the derivative of \alpha with respect to a given component x_k we get the two sums
\frac{\partial \alpha}{\partial x_k} = \sum_{i=0}^{n-1}a_{ik}x_i+\sum_{j=0}^{n-1}a_{kj}x_j,
for \forall k =0,1,2,\dots,n-1. We identify these sums as
\frac{\partial \alpha}{\partial \boldsymbol{x}} = \boldsymbol{x}^T\left(\boldsymbol{A}^T+\boldsymbol{A}\right).
If the matrix \boldsymbol{A} is symmetric, that is \boldsymbol{A}=\boldsymbol{A}^T, we have
\frac{\partial \alpha}{\partial \boldsymbol{x}} = 2\boldsymbol{x}^T\boldsymbol{A}.
Example 4
We let the scalar \alpha be defined by
\alpha = \boldsymbol{y}^T\boldsymbol{x},
where both \boldsymbol{y} and \boldsymbol{x} have the same length n, or if we
wish to think of them as column vectors, they have dimensions $n\times
1$. We assume that both \boldsymbol{y} and \boldsymbol{x} depend on a vector
\boldsymbol{z} of the same length. To calculate the derivative of \alpha
with respect to a given component z_k we need first to write out the
inner product that defines \alpha as
\alpha = \sum_{i=0}^{n-1}y_ix_i,
and the partial derivative
\frac{\partial \alpha}{\partial z_k} = \sum_{i=0}^{n-1}\left(x_i\frac{\partial y_i}{\partial z_k}+y_i\frac{\partial x_i}{\partial z_k}\right),
for \forall k =0,1,2,\dots,n-1. We can rewrite the partial derivative in a more compact form as
\frac{\partial \alpha}{\partial \boldsymbol{z}} = \boldsymbol{x}^T\frac{\partial \boldsymbol{y}}{\partial \boldsymbol{z}}+\boldsymbol{y}^T\frac{\partial \boldsymbol{x}}{\partial \boldsymbol{z}},
and if \boldsymbol{y}=\boldsymbol{x} we have
\frac{\partial \alpha}{\partial \boldsymbol{z}} = 2\boldsymbol{x}^T\frac{\partial \boldsymbol{x}}{\partial \boldsymbol{z}}.
The mean squared error and its derivative
We defined earlier a possible cost function using the mean squared error
C(\boldsymbol{\beta})=\frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)^T\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)\right\},
or using the design/feature matrix \boldsymbol{X} we have the more compact matrix-vector
C(\boldsymbol{\beta})=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right\}.
We note that the design matrix \boldsymbol{X} does not depend on the unknown parameters defined by the vector \boldsymbol{\beta}.
We are now interested in minimizing the cost function with respect to the unknown parameters \boldsymbol{\beta}.
The mean squared error is a scalar and if we use the results from example three above, we can define a new vector
\boldsymbol{w}=\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta},
which depends on \boldsymbol{\beta}. We rewrite the cost function as
C(\boldsymbol{\beta})=\frac{1}{n}\boldsymbol{w}^T\boldsymbol{w},
with partial derivative
\frac{\partial C(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}=\frac{2}{n}\boldsymbol{w}^T\frac{\partial \boldsymbol{w}}{\partial \boldsymbol{\beta}},
and using that
\frac{\partial \boldsymbol{w}}{\partial \boldsymbol{\beta}}=-\boldsymbol{X},
where we used the result from example two above. Inserting the last expression we obtain
\frac{\partial C(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}=-\frac{2}{n}\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\boldsymbol{X},
or as
\frac{\partial C(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}^T}=-\frac{2}{n}\boldsymbol{X}^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right).
Other useful relations
We list here some other useful relations we may encounter (recall that vectors are defined by boldfaced low-key letters)
\frac{\partial (\boldsymbol{b}^T\boldsymbol{a})}{\partial \boldsymbol{a}} = \boldsymbol{b},
\frac{\partial tr(\boldsymbol{B}\boldsymbol{A})}{\partial \boldsymbol{A}} = \boldsymbol{B}^T,
\frac{\partial \log{\vert\boldsymbol{A}\vert}}{\partial \boldsymbol{A}} = (\boldsymbol{A}^{-1})^T.
Meet the Hessian Matrix
A very important matrix we will meet again and again in machine
learning is the Hessian. It is given by the second derivative of the
cost function with respect to the parameters \boldsymbol{\beta}. Using the above
expression for derivatives of vectors and matrices, we find that the
second derivative of the mean squared error as cost function is,
\frac{\partial}{\partial \boldsymbol{\beta}}\frac{\partial C(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}^T} =\frac{\partial}{\partial \boldsymbol{\beta}}\left[-\frac{2}{n}\boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right]=\frac{2}{n}\boldsymbol{X}^T\boldsymbol{X}.
The Hessian matrix plays an important role and is defined here as