update
This commit is contained in:
Morten Hjorth-Jensen
Binary file not shown.
Binary file not shown.
File diff suppressed because it is too large
Load Diff
File diff suppressed because one or more lines are too long
@@ -390,7 +390,6 @@ document.write(`
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#example-3">Example 3</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#example-4">Example 4</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-mean-squared-error-and-its-derivative">The mean squared error and its derivative</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#other-useful-relations">Other useful relations</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#meet-the-hessian-matrix">Meet the Hessian Matrix</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#id2">Interpretations and optimizing our parameters</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#example-relevant-for-the-exercises">Example relevant for the exercises</a></li>
|
||||
@@ -413,6 +412,8 @@ document.write(`
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-singular-value-decomposition">The singular value decomposition</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#linear-regression-problems">Linear Regression Problems</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#fixing-the-singularity">Fixing the singularity</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#ridge-and-lasso-regression">Ridge and LASSO Regression</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#deriving-the-ridge-regression-equations">Deriving the Ridge Regression Equations</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#basic-math-of-the-svd">Basic math of the SVD</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-svd-a-fantastic-algorithm">The SVD, a Fantastic Algorithm</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#economy-size-svd">Economy-size SVD</a></li>
|
||||
@@ -429,13 +430,11 @@ document.write(`
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#covariance-matrix-examples">Covariance Matrix Examples</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#correlation-matrix">Correlation Matrix</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#correlation-matrix-with-pandas">Correlation Matrix with Pandas</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#correlation-matrix-with-pandas-and-the-franke-function">Correlation Matrix with Pandas and the Franke function</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#rewriting-the-covariance-and-or-correlation-matrix">Rewriting the Covariance and/or Correlation Matrix</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#linking-with-the-svd">Linking with the SVD</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#what-does-it-mean">What does it mean?</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#and-finally-boldsymbol-x-boldsymbol-x-t">And finally <span class="math notranslate nohighlight">\(\boldsymbol{X}\boldsymbol{X}^T\)</span></a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#ridge-and-lasso-regression">Ridge and LASSO Regression</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#deriving-the-ridge-regression-equations">Deriving the Ridge Regression Equations</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#back-to-ridge-and-lasso-regression">Back to Ridge and LASSO Regression</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#interpreting-the-ridge-results">Interpreting the Ridge results</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#more-interpretations">More interpretations</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#deriving-the-lasso-regression-equations">Deriving the Lasso Regression Equations</a></li>
|
||||
@@ -461,10 +460,10 @@ doconce format html week35.do.txt --no_mako -->
|
||||
<p>The main topics are:</p>
|
||||
<ol class="arabic simple">
|
||||
<li><p>Brief repetition from last week</p></li>
|
||||
<li><p>Discussions of the equations for ordinary least squares</p></li>
|
||||
<li><p>Discussions of the equations for ordinary least squares (<strong>OLS</strong>)</p></li>
|
||||
<li><p>Discussion on how to prepare data and examples of applications of linear regression</p></li>
|
||||
<li><p>Material for the lecture on Monday: Mathematical interpretations of linear regression</p></li>
|
||||
<li><p>Monday: Ridge and Lasso regression and Singular Value Decomposition</p></li>
|
||||
<li><p>Mathematical interpretations of OLS</p></li>
|
||||
<li><p>Introduction of Ridge and Lasso regression</p></li>
|
||||
</ol>
|
||||
<section id="reading-recommendations">
|
||||
<h3>Reading recommendations:<a class="headerlink" href="#reading-recommendations" title="Link to this heading">#</a></h3>
|
||||
@@ -474,7 +473,7 @@ doconce format html week35.do.txt --no_mako -->
|
||||
<!-- o [Video of lecture](https://youtu.be/VKakN-e4aUA) -->
|
||||
<!-- o [Video for exercises week 35](https://youtu.be/yiY0OltU1s8) -->
|
||||
<ol class="arabic simple" start="2">
|
||||
<li><p>Goodfellow, Bengio and Courville, Deep Learning, chapter 2 on linear algebra and sections 3.1-3.10 on elements of statistics (background)</p></li>
|
||||
<li><p>Goodfellow, Bengio and Courville, Deep Learning, chapter 2 on linear algebra</p></li>
|
||||
<li><p>Raschka et al on preprocessing of data, relevant for exercise 3 this week, see chapter 4.</p></li>
|
||||
<li><p>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”<a class="reference external" href="https://mml-book.github.io/">https://mml-book.github.io/</a>” (section 5.5 on derivatives) is very useful for exercise 1 this coming week.</p></li>
|
||||
</ol>
|
||||
@@ -568,16 +567,16 @@ y_{i}=\langle y_i \rangle = \theta_0x_{i,0}+\theta_1x_{i,1}+\theta_2x_{i,2}+\dot
|
||||
\]</div>
|
||||
<p>where <span class="math notranslate nohighlight">\(\langle y_i \rangle\)</span> is the mean value. Keep in mind also that
|
||||
till now we have treated <span class="math notranslate nohighlight">\(y_i\)</span> as the exact value. Normally, the
|
||||
response (dependent or outcome) variable <span class="math notranslate nohighlight">\(y_i\)</span> is the outcome of a
|
||||
output (response, target, dependent or outcome) variable <span class="math notranslate nohighlight">\(y_i\)</span> 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 <span class="math notranslate nohighlight">\(y_i\)</span> as our exact value for the response variable.</p>
|
||||
will treat <span class="math notranslate nohighlight">\(y_i\)</span> as our exact value for the output variable.</p>
|
||||
<p>In order to find the parameters <span class="math notranslate nohighlight">\(\theta_i\)</span> we will then minimize the spread of <span class="math notranslate nohighlight">\(C(\boldsymbol{\theta})\)</span>, that is we are going to solve the problem</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
{\displaystyle \min_{\boldsymbol{\theta}\in
|
||||
\hat{\boldsymbol{\theta}}={\displaystyle \min_{\boldsymbol{\theta}\in
|
||||
{\mathbb{R}}^{p}}}\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)\right\}.
|
||||
\]</div>
|
||||
<p>In practical terms it means we will require</p>
|
||||
@@ -611,7 +610,7 @@ will treat <span class="math notranslate nohighlight">\(y_i\)</span> as our exac
|
||||
<p>and if the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span> is invertible we have the solution</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\boldsymbol{\theta} =\left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}.
|
||||
\hat{\boldsymbol{\theta}} =\left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}.
|
||||
\]</div>
|
||||
<p>We note also that since our design matrix is defined as <span class="math notranslate nohighlight">\(\boldsymbol{X}\in
|
||||
{\mathbb{R}}^{n\times p}\)</span>, the product <span class="math notranslate nohighlight">\(\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 <strong>LU</strong> decomposition or <strong>Singular Value Decomposition</strong> (SVD) for finding the inverse of the matrix
|
||||
<span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span>. This is discussed on Thursday this week.</p>
|
||||
<p><strong>Small question</strong>: Do you think the example we have at hand here (the nuclear binding energies) can lead to problems in inverting the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span>? What kind of problems can we expect?</p>
|
||||
<span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span>.</p>
|
||||
<p><strong>Small question</strong>: When inverting the matrix $\boldsymbol{X}^T\boldsymbol{X}, what kind of problems can we expect?</p>
|
||||
</section>
|
||||
<section id="some-useful-matrix-and-vector-expressions">
|
||||
<h2>Some useful matrix and vector expressions<a class="headerlink" href="#some-useful-matrix-and-vector-expressions" title="Link to this heading">#</a></h2>
|
||||
@@ -821,26 +820,6 @@ C(\boldsymbol{\theta})=\frac{1}{n}\boldsymbol{w}^T\boldsymbol{w},
|
||||
\frac{\partial C(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}^T}=-\frac{2}{n}\boldsymbol{X}^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right).
|
||||
\]</div>
|
||||
</section>
|
||||
<section id="other-useful-relations">
|
||||
<h2>Other useful relations<a class="headerlink" href="#other-useful-relations" title="Link to this heading">#</a></h2>
|
||||
<p>We list here some other useful relations we may encounter (recall that vectors are defined by boldfaced low-key letters)</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\frac{\partial (\boldsymbol{x}^T\boldsymbol{a})}{\partial \boldsymbol{x}} = \boldsymbol{a}^T,
|
||||
\]</div>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\frac{\partial (\boldsymbol{a}^T\boldsymbol{x})}{\partial \boldsymbol{x}} = \boldsymbol{a}^T,
|
||||
\]</div>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\frac{\partial tr(\boldsymbol{B}\boldsymbol{A})}{\partial \boldsymbol{A}} = \boldsymbol{B}^T,
|
||||
\]</div>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\frac{\partial \log{\vert\boldsymbol{A}\vert}}{\partial \boldsymbol{A}} = (\boldsymbol{A}^{-1})^T.
|
||||
\]</div>
|
||||
</section>
|
||||
<section id="meet-the-hessian-matrix">
|
||||
<h2>Meet the Hessian Matrix<a class="headerlink" href="#meet-the-hessian-matrix" title="Link to this heading">#</a></h2>
|
||||
<p>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,</p>
|
||||
\]</div>
|
||||
<p>For ordinary least squares, it is inversely proportional (derivation
|
||||
next week) with the variance of the optimal parameters
|
||||
<span class="math notranslate nohighlight">\(\hat{\boldsymbol{\theta}}\)</span>. Furthermore, we will see later this week that it is
|
||||
<span class="math notranslate nohighlight">\(\hat{\boldsymbol{\theta}}\)</span>. Furthermore, we will see next week that it is
|
||||
(aside the factor <span class="math notranslate nohighlight">\(1/n\)</span>) 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
|
||||
<strong>Linear algebra question:</strong> 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!).</p>
|
||||
problem. We will discuss this in greater detail next week when we introduce gradient methods.</p>
|
||||
<p><strong>Linear algebra question:</strong> 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!).</p>
|
||||
</section>
|
||||
<section id="id2">
|
||||
<h2>Interpretations and optimizing our parameters<a class="headerlink" href="#id2" title="Link to this heading">#</a></h2>
|
||||
@@ -903,7 +881,7 @@ We assume our data can represented by a fourth-order polynomial. For the <span c
|
||||
<p>It is rather straightforward to implement the matrix inversion and obtain the parameters <span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span>. After having defined the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> and the outputs <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> we have</p>
|
||||
<div class="cell docutils container">
|
||||
<div class="cell_input docutils container">
|
||||
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span># matrix inversion to find beta
|
||||
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span># 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
|
||||
</pre></div>
|
||||
</div>
|
||||
</div>
|
||||
@@ -1021,16 +999,16 @@ X[:,3] = x**3
|
||||
X[:,4] = x**4
|
||||
# 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)
|
||||
# matrix inversion to find beta
|
||||
beta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train
|
||||
print(beta)
|
||||
# matrix inversion to find theta
|
||||
theta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train
|
||||
print(theta)
|
||||
# and then make the prediction
|
||||
ytilde = X_train @ beta
|
||||
ytilde = X_train @ theta
|
||||
print("Training R2")
|
||||
print(R2(y_train,ytilde))
|
||||
print("Training MSE")
|
||||
print(MSE(y_train,ytilde))
|
||||
ypredict = X_test @ beta
|
||||
ypredict = X_test @ theta
|
||||
print("Test R2")
|
||||
print(R2(y_test,ypredict))
|
||||
print("Test MSE")
|
||||
@@ -1319,7 +1297,7 @@ decomposition.</p>
|
||||
<p>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.</p>
|
||||
algorithm based on say LU, QR or Cholesky decomposition may lead to singularities. We will see examples of this below and in other examples.</p>
|
||||
<p>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.</p>
|
||||
@@ -1404,9 +1382,175 @@ the regression parameters <span class="math notranslate nohighlight">\(\theta_i\
|
||||
\]</div>
|
||||
<p>where <span class="math notranslate nohighlight">\(\boldsymbol{I}\)</span> is the identity matrix. When we discuss <strong>Ridge</strong> regression this is actually what we end up evaluating. The parameter <span class="math notranslate nohighlight">\(\lambda\)</span> is called a hyperparameter. More about this later.</p>
|
||||
</section>
|
||||
<section id="ridge-and-lasso-regression">
|
||||
<h2>Ridge and LASSO Regression<a class="headerlink" href="#ridge-and-lasso-regression" title="Link to this heading">#</a></h2>
|
||||
<p>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</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
{\displaystyle \min_{\boldsymbol{\theta}\in {\mathbb{R}}^{p}}}\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)\right\}.
|
||||
\]</div>
|
||||
<p>or we can state it as</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
{\displaystyle \min_{\boldsymbol{\theta}\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 \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\vert\vert_2^2,
|
||||
\]</div>
|
||||
<p>where we have used the definition of a norm-2 vector, that is</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\vert\vert \boldsymbol{x}\vert\vert_2 = \sqrt{\sum_i x_i^2}.
|
||||
\]</div>
|
||||
<p>By minimizing the above equation with respect to the parameters
|
||||
<span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span> we could then obtain an analytical expression for the
|
||||
parameters <span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span>. We can add a regularization parameter <span class="math notranslate nohighlight">\(\lambda\)</span> by
|
||||
defining a new cost function to be optimized, that is</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
{\displaystyle \min_{\boldsymbol{\theta}\in
|
||||
{\mathbb{R}}^{p}}}\frac{1}{n}\vert\vert \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\vert\vert_2^2+\lambda\vert\vert \boldsymbol{\theta}\vert\vert_2^2
|
||||
\]</div>
|
||||
<p>which leads to the Ridge regression minimization problem where we
|
||||
require that <span class="math notranslate nohighlight">\(\vert\vert \boldsymbol{\theta}\vert\vert_2^2\le t\)</span>, where <span class="math notranslate nohighlight">\(t\)</span> is
|
||||
a finite number larger than zero. By defining</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
C(\boldsymbol{X},\boldsymbol{\theta})=\frac{1}{n}\vert\vert \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\vert\vert_2^2+\lambda\vert\vert \boldsymbol{\theta}\vert\vert_1,
|
||||
\]</div>
|
||||
<p>we have a new optimization equation</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
{\displaystyle \min_{\boldsymbol{\theta}\in
|
||||
{\mathbb{R}}^{p}}}\frac{1}{n}\vert\vert \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\vert\vert_2^2+\lambda\vert\vert \boldsymbol{\theta}\vert\vert_1
|
||||
\]</div>
|
||||
<p>which leads to Lasso regression. Lasso stands for least absolute shrinkage and selection operator.</p>
|
||||
<p>Here we have defined the norm-1 as</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\vert\vert \boldsymbol{x}\vert\vert_1 = \sum_i \vert x_i\vert.
|
||||
\]</div>
|
||||
</section>
|
||||
<section id="deriving-the-ridge-regression-equations">
|
||||
<h2>Deriving the Ridge Regression Equations<a class="headerlink" href="#deriving-the-ridge-regression-equations" title="Link to this heading">#</a></h2>
|
||||
<p>Using the matrix-vector expression for Ridge regression and dropping the parameter <span class="math notranslate nohighlight">\(1/n\)</span> in front of the standard means squared error equation, we have</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
C(\boldsymbol{X},\boldsymbol{\theta})=\left\{(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta})^T(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta})\right\}+\lambda\boldsymbol{\theta}^T\boldsymbol{\theta},
|
||||
\]</div>
|
||||
<p>and
|
||||
taking the derivatives with respect to <span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span> we obtain then
|
||||
a slightly modified matrix inversion problem which for finite values
|
||||
of <span class="math notranslate nohighlight">\(\lambda\)</span> does not suffer from singularity problems. We obtain
|
||||
the optimal parameters</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\hat{\boldsymbol{\theta}}_{\mathrm{Ridge}} = \left(\boldsymbol{X}^T\boldsymbol{X}+\lambda\boldsymbol{I}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y},
|
||||
\]</div>
|
||||
<p>with <span class="math notranslate nohighlight">\(\boldsymbol{I}\)</span> being a <span class="math notranslate nohighlight">\(p\times p\)</span> identity matrix with the constraint that</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\sum_{i=0}^{p-1} \theta_i^2 \leq t,
|
||||
\]</div>
|
||||
<p>with <span class="math notranslate nohighlight">\(t\)</span> a finite positive number.</p>
|
||||
<p>If we keep the <span class="math notranslate nohighlight">\(1/n\)</span> factor, the equation for the optimal <span class="math notranslate nohighlight">\(\theta\)</span> changes to</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\hat{\boldsymbol{\theta}}_{\mathrm{Ridge}} = \left(\boldsymbol{X}^T\boldsymbol{X}+n\lambda\boldsymbol{I}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}.
|
||||
\]</div>
|
||||
<p>In many textbooks the <span class="math notranslate nohighlight">\(1/n\)</span> term is often omitted. Note that a library like <strong>Scikit-Learn</strong> does not include the <span class="math notranslate nohighlight">\(1/n\)</span> factor in the setup of the cost function.</p>
|
||||
<p>When we compare this with the ordinary least squares result we have</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\hat{\boldsymbol{\theta}}_{\mathrm{OLS}} = \left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y},
|
||||
\]</div>
|
||||
<p>which can lead to singular matrices. However, with the SVD, we can always compute the inverse of the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span>.</p>
|
||||
<p>We see that Ridge regression is nothing but the standard OLS with a
|
||||
modified diagonal term added to <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span>. The consequences, in
|
||||
particular for our discussion of the bias-variance tradeoff are rather
|
||||
interesting. We will see that for specific values of <span class="math notranslate nohighlight">\(\lambda\)</span>, we may
|
||||
even reduce the variance of the optimal parameters <span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span>. These topics and other related ones, will be discussed after the more linear algebra oriented analysis here.</p>
|
||||
<p>When we have discussed the singular value decomposition of the design
|
||||
matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>, we will in turn perform a more rigorous mathematical
|
||||
discussion of Ridge regression.</p>
|
||||
<p>The code here is a simple demonstration of how to implement Ridge regression with our own code and compare this with scikit-learn.</p>
|
||||
<div class="cell docutils container">
|
||||
<div class="cell_input docutils container">
|
||||
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span>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()
|
||||
</pre></div>
|
||||
</div>
|
||||
</div>
|
||||
</div>
|
||||
<p>The results here agree when we force <strong>Scikit-Learn</strong>’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.</p>
|
||||
</section>
|
||||
<section id="basic-math-of-the-svd">
|
||||
<h2>Basic math of the SVD<a class="headerlink" href="#basic-math-of-the-svd" title="Link to this heading">#</a></h2>
|
||||
<p>From standard linear algebra we know that a square matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> can be diagonalized if and only it is
|
||||
<p>From standard linear algebra we know that a square matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> can be diagonalized if and only if it is
|
||||
a so-called <a class="reference external" href="https://en.wikipedia.org/wiki/Normal_matrix">normal matrix</a>, that is if <span class="math notranslate nohighlight">\(\boldsymbol{X}\in {\mathbb{R}}^{n\times n}\)</span>
|
||||
we have <span class="math notranslate nohighlight">\(\boldsymbol{X}\boldsymbol{X}^T=\boldsymbol{X}^T\boldsymbol{X}\)</span> or if <span class="math notranslate nohighlight">\(\boldsymbol{X}\in {\mathbb{C}}^{n\times n}\)</span> we have <span class="math notranslate nohighlight">\(\boldsymbol{X}\boldsymbol{X}^{\dagger}=\boldsymbol{X}^{\dagger}\boldsymbol{X}\)</span>.
|
||||
The matrix has then a set of eigenpairs</p>
|
||||
@@ -1946,66 +2090,6 @@ print(correlation_matrix)
|
||||
</div>
|
||||
</div>
|
||||
</div>
|
||||
<p>We expand this model to the Franke function discussed above.</p>
|
||||
</section>
|
||||
<section id="correlation-matrix-with-pandas-and-the-franke-function">
|
||||
<h2>Correlation Matrix with Pandas and the Franke function<a class="headerlink" href="#correlation-matrix-with-pandas-and-the-franke-function" title="Link to this heading">#</a></h2>
|
||||
<div class="cell docutils container">
|
||||
<div class="cell_input docutils container">
|
||||
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span># 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)
|
||||
</pre></div>
|
||||
</div>
|
||||
</div>
|
||||
</div>
|
||||
<p>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 <span class="math notranslate nohighlight">\(n\)</span>).</p>
|
||||
<p>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.</p>
|
||||
</section>
|
||||
<section id="rewriting-the-covariance-and-or-correlation-matrix">
|
||||
<h2>Rewriting the Covariance and/or Correlation Matrix<a class="headerlink" href="#rewriting-the-covariance-and-or-correlation-matrix" title="Link to this heading">#</a></h2>
|
||||
@@ -2136,8 +2220,8 @@ measure how much correlations are contained in the rows of <span class="math not
|
||||
of our data (the columns of <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>, the quantity of interest for us are the non-zero singular
|
||||
values and the column vectors of <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span>.</p>
|
||||
</section>
|
||||
<section id="ridge-and-lasso-regression">
|
||||
<h2>Ridge and LASSO Regression<a class="headerlink" href="#ridge-and-lasso-regression" title="Link to this heading">#</a></h2>
|
||||
<section id="back-to-ridge-and-lasso-regression">
|
||||
<h2>Back to Ridge and LASSO Regression<a class="headerlink" href="#back-to-ridge-and-lasso-regression" title="Link to this heading">#</a></h2>
|
||||
<p>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</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
@@ -2183,42 +2267,7 @@ C(\boldsymbol{X},\boldsymbol{\theta})=\frac{1}{n}\vert\vert \boldsymbol{y}-\bold
|
||||
\[
|
||||
\vert\vert \boldsymbol{x}\vert\vert_1 = \sum_i \vert x_i\vert.
|
||||
\]</div>
|
||||
</section>
|
||||
<section id="deriving-the-ridge-regression-equations">
|
||||
<h2>Deriving the Ridge Regression Equations<a class="headerlink" href="#deriving-the-ridge-regression-equations" title="Link to this heading">#</a></h2>
|
||||
<p>Using the matrix-vector expression for Ridge regression and dropping the parameter <span class="math notranslate nohighlight">\(1/n\)</span> in front of the standard means squared error equation, we have</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
C(\boldsymbol{X},\boldsymbol{\theta})=\left\{(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta})^T(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta})\right\}+\lambda\boldsymbol{\theta}^T\boldsymbol{\theta},
|
||||
\]</div>
|
||||
<p>and
|
||||
taking the derivatives with respect to <span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span> we obtain then
|
||||
a slightly modified matrix inversion problem which for finite values
|
||||
of <span class="math notranslate nohighlight">\(\lambda\)</span> does not suffer from singularity problems. We obtain
|
||||
the optimal parameters</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\hat{\boldsymbol{\theta}}_{\mathrm{Ridge}} = \left(\boldsymbol{X}^T\boldsymbol{X}+\lambda\boldsymbol{I}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y},
|
||||
\]</div>
|
||||
<p>with <span class="math notranslate nohighlight">\(\boldsymbol{I}\)</span> being a <span class="math notranslate nohighlight">\(p\times p\)</span> identity matrix with the constraint that</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\sum_{i=0}^{p-1} \theta_i^2 \leq t,
|
||||
\]</div>
|
||||
<p>with <span class="math notranslate nohighlight">\(t\)</span> a finite positive number.</p>
|
||||
<p>If we keep the <span class="math notranslate nohighlight">\(1/n\)</span> factor, the equation for the optimal <span class="math notranslate nohighlight">\(\theta\)</span> changes to</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\hat{\boldsymbol{\theta}}_{\mathrm{Ridge}} = \left(\boldsymbol{X}^T\boldsymbol{X}+n\lambda\boldsymbol{I}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}.
|
||||
\]</div>
|
||||
<p>In many textbooks the <span class="math notranslate nohighlight">\(1/n\)</span> term is often omitted. Note that a library like <strong>Scikit-Learn</strong> does not include the <span class="math notranslate nohighlight">\(1/n\)</span> factor in the setup of the cost function.</p>
|
||||
<p>When we compare this with the ordinary least squares result we have</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\hat{\boldsymbol{\theta}}_{\mathrm{OLS}} = \left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y},
|
||||
\]</div>
|
||||
<p>which can lead to singular matrices. However, with the SVD, we can always compute the inverse of the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span>.</p>
|
||||
<p>We see that Ridge regression is nothing but the standard OLS with a
|
||||
<p>Ridge regression, as discussed above, is nothing but the standard OLS with a
|
||||
modified diagonal term added to <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span>. The consequences, in
|
||||
particular for our discussion of the bias-variance tradeoff are rather
|
||||
interesting. We will see that for specific values of <span class="math notranslate nohighlight">\(\lambda\)</span>, 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}.
|
||||
\]</div>
|
||||
<p>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 <a class="reference external" href="https://cvxopt.org/">CVXOPT</a>. We will discuss this later.</p>
|
||||
<p>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 <a class="reference external" href="https://cvxopt.org/">CVXOPT</a>. We will discuss how to code LASSO regression next week, when we have introduced gradient methods.</p>
|
||||
</section>
|
||||
</section>
|
||||
|
||||
@@ -2376,7 +2425,6 @@ C(\boldsymbol{X},\boldsymbol{\theta})=\frac{1}{n}\left\{(\boldsymbol{y}-\boldsym
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#example-3">Example 3</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#example-4">Example 4</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-mean-squared-error-and-its-derivative">The mean squared error and its derivative</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#other-useful-relations">Other useful relations</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#meet-the-hessian-matrix">Meet the Hessian Matrix</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#id2">Interpretations and optimizing our parameters</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#example-relevant-for-the-exercises">Example relevant for the exercises</a></li>
|
||||
@@ -2399,6 +2447,8 @@ C(\boldsymbol{X},\boldsymbol{\theta})=\frac{1}{n}\left\{(\boldsymbol{y}-\boldsym
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-singular-value-decomposition">The singular value decomposition</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#linear-regression-problems">Linear Regression Problems</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#fixing-the-singularity">Fixing the singularity</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#ridge-and-lasso-regression">Ridge and LASSO Regression</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#deriving-the-ridge-regression-equations">Deriving the Ridge Regression Equations</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#basic-math-of-the-svd">Basic math of the SVD</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-svd-a-fantastic-algorithm">The SVD, a Fantastic Algorithm</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#economy-size-svd">Economy-size SVD</a></li>
|
||||
@@ -2415,13 +2465,11 @@ C(\boldsymbol{X},\boldsymbol{\theta})=\frac{1}{n}\left\{(\boldsymbol{y}-\boldsym
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#covariance-matrix-examples">Covariance Matrix Examples</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#correlation-matrix">Correlation Matrix</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#correlation-matrix-with-pandas">Correlation Matrix with Pandas</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#correlation-matrix-with-pandas-and-the-franke-function">Correlation Matrix with Pandas and the Franke function</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#rewriting-the-covariance-and-or-correlation-matrix">Rewriting the Covariance and/or Correlation Matrix</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#linking-with-the-svd">Linking with the SVD</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#what-does-it-mean">What does it mean?</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#and-finally-boldsymbol-x-boldsymbol-x-t">And finally <span class="math notranslate nohighlight">\(\boldsymbol{X}\boldsymbol{X}^T\)</span></a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#ridge-and-lasso-regression">Ridge and LASSO Regression</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#deriving-the-ridge-regression-equations">Deriving the Ridge Regression Equations</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#back-to-ridge-and-lasso-regression">Back to Ridge and LASSO Regression</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#interpreting-the-ridge-results">Interpreting the Ridge results</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#more-interpretations">More interpretations</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#deriving-the-lasso-regression-equations">Deriving the Lasso Regression Equations</a></li>
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
+773
-649
File diff suppressed because it is too large
Load Diff
Reference in New Issue
Block a user