lars/FYS-STK4155
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
9e2f8c701a2af6e145ba0488017dd198c6883aa9
FYS-STK4155/doc/LectureNotes/_build/html/week36.html
T
Morten Hjorth-Jensen 9e2f8c701a typo in title
2025-09-07 22:10:57 +02:00

2564 lines
165 KiB
HTML
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
<!DOCTYPE html>
<html lang="en" data-content_root="./" >
<head>
<meta charset="utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1.0" /><meta name="viewport" content="width=device-width, initial-scale=1" />
<title>Week 36: Linear Regression and Gradient descent &#8212; Applied Data Analysis and Machine Learning</title>
<script data-cfasync="false">
document.documentElement.dataset.mode = localStorage.getItem("mode") || "";
document.documentElement.dataset.theme = localStorage.getItem("theme") || "";
</script>
<!-- Loaded before other Sphinx assets -->
<link href="_static/styles/theme.css?digest=dfe6caa3a7d634c4db9b" rel="stylesheet" />
<link href="_static/styles/bootstrap.css?digest=dfe6caa3a7d634c4db9b" rel="stylesheet" />
<link href="_static/styles/pydata-sphinx-theme.css?digest=dfe6caa3a7d634c4db9b" rel="stylesheet" />
<link href="_static/vendor/fontawesome/6.5.2/css/all.min.css?digest=dfe6caa3a7d634c4db9b" rel="stylesheet" />
<link rel="preload" as="font" type="font/woff2" crossorigin href="_static/vendor/fontawesome/6.5.2/webfonts/fa-solid-900.woff2" />
<link rel="preload" as="font" type="font/woff2" crossorigin href="_static/vendor/fontawesome/6.5.2/webfonts/fa-brands-400.woff2" />
<link rel="preload" as="font" type="font/woff2" crossorigin href="_static/vendor/fontawesome/6.5.2/webfonts/fa-regular-400.woff2" />
<link rel="stylesheet" type="text/css" href="_static/pygments.css?v=fa44fd50" />
<link rel="stylesheet" type="text/css" href="_static/styles/sphinx-book-theme.css?v=eba8b062" />
<link rel="stylesheet" type="text/css" href="_static/togglebutton.css?v=13237357" />
<link rel="stylesheet" type="text/css" href="_static/copybutton.css?v=76b2166b" />
<link rel="stylesheet" type="text/css" href="_static/mystnb.8ecb98da25f57f5357bf6f572d296f466b2cfe2517ffebfabe82451661e28f02.css?v=6644e6bb" />
<link rel="stylesheet" type="text/css" href="_static/sphinx-thebe.css?v=4fa983c6" />
<link rel="stylesheet" type="text/css" href="_static/sphinx-design.min.css?v=95c83b7e" />
<!-- Pre-loaded scripts that we'll load fully later -->
<link rel="preload" as="script" href="_static/scripts/bootstrap.js?digest=dfe6caa3a7d634c4db9b" />
<link rel="preload" as="script" href="_static/scripts/pydata-sphinx-theme.js?digest=dfe6caa3a7d634c4db9b" />
<script src="_static/vendor/fontawesome/6.5.2/js/all.min.js?digest=dfe6caa3a7d634c4db9b"></script>
<script src="_static/documentation_options.js?v=9eb32ce0"></script>
<script src="_static/doctools.js?v=9a2dae69"></script>
<script src="_static/sphinx_highlight.js?v=dc90522c"></script>
<script src="_static/clipboard.min.js?v=a7894cd8"></script>
<script src="_static/copybutton.js?v=f281be69"></script>
<script src="_static/scripts/sphinx-book-theme.js?v=887ef09a"></script>
<script>let toggleHintShow = 'Click to show';</script>
<script>let toggleHintHide = 'Click to hide';</script>
<script>let toggleOpenOnPrint = 'true';</script>
<script src="_static/togglebutton.js?v=4a39c7ea"></script>
<script>var togglebuttonSelector = '.toggle, .admonition.dropdown';</script>
<script src="_static/design-tabs.js?v=f930bc37"></script>
<script>const THEBE_JS_URL = "https://unpkg.com/thebe@0.8.2/lib/index.js"; const thebe_selector = ".thebe,.cell"; const thebe_selector_input = "pre"; const thebe_selector_output = ".output, .cell_output"</script>
<script async="async" src="_static/sphinx-thebe.js?v=c100c467"></script>
<script>var togglebuttonSelector = '.toggle, .admonition.dropdown';</script>
<script>const THEBE_JS_URL = "https://unpkg.com/thebe@0.8.2/lib/index.js"; const thebe_selector = ".thebe,.cell"; const thebe_selector_input = "pre"; const thebe_selector_output = ".output, .cell_output"</script>
<script>window.MathJax = {"options": {"processHtmlClass": "tex2jax_process|mathjax_process|math|output_area"}}</script>
<script defer="defer" src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"></script>
<script>DOCUMENTATION_OPTIONS.pagename = 'week36';</script>
<link rel="index" title="Index" href="genindex.html" />
<link rel="search" title="Search" href="search.html" />
<link rel="next" title="Exercises week 37" href="exercisesweek37.html" />
<link rel="prev" title="Exercises week 36" href="exercisesweek36.html" />
<meta name="viewport" content="width=device-width, initial-scale=1"/>
<meta name="docsearch:language" content="en"/>
</head>
<body data-bs-spy="scroll" data-bs-target=".bd-toc-nav" data-offset="180" data-bs-root-margin="0px 0px -60%" data-default-mode="">
<div id="pst-skip-link" class="skip-link d-print-none"><a href="#main-content">Skip to main content</a></div>
<div id="pst-scroll-pixel-helper"></div>
<button type="button" class="btn rounded-pill" id="pst-back-to-top">
<i class="fa-solid fa-arrow-up"></i>Back to top</button>
<input type="checkbox"
class="sidebar-toggle"
id="pst-primary-sidebar-checkbox"/>
<label class="overlay overlay-primary" for="pst-primary-sidebar-checkbox"></label>
<input type="checkbox"
class="sidebar-toggle"
id="pst-secondary-sidebar-checkbox"/>
<label class="overlay overlay-secondary" for="pst-secondary-sidebar-checkbox"></label>
<div class="search-button__wrapper">
<div class="search-button__overlay"></div>
<div class="search-button__search-container">
<form class="bd-search d-flex align-items-center"
action="search.html"
method="get">
<i class="fa-solid fa-magnifying-glass"></i>
<input type="search"
class="form-control"
name="q"
id="search-input"
placeholder="Search this book..."
aria-label="Search this book..."
autocomplete="off"
autocorrect="off"
autocapitalize="off"
spellcheck="false"/>
<span class="search-button__kbd-shortcut"><kbd class="kbd-shortcut__modifier">Ctrl</kbd>+<kbd>K</kbd></span>
</form></div>
</div>
<div class="pst-async-banner-revealer d-none">
<aside id="bd-header-version-warning" class="d-none d-print-none" aria-label="Version warning"></aside>
</div>
<header class="bd-header navbar navbar-expand-lg bd-navbar d-print-none">
</header>
<div class="bd-container">
<div class="bd-container__inner bd-page-width">
<div class="bd-sidebar-primary bd-sidebar">
<div class="sidebar-header-items sidebar-primary__section">
</div>
<div class="sidebar-primary-items__start sidebar-primary__section">
<div class="sidebar-primary-item">
<a class="navbar-brand logo" href="intro.html">
<img src="_static/logo.png" class="logo__image only-light" alt="Applied Data Analysis and Machine Learning - Home"/>
<script>document.write(`<img src="_static/logo.png" class="logo__image only-dark" alt="Applied Data Analysis and Machine Learning - Home"/>`);</script>
</a></div>
<div class="sidebar-primary-item">
<script>
document.write(`
<button class="btn search-button-field search-button__button" title="Search" aria-label="Search" data-bs-placement="bottom" data-bs-toggle="tooltip">
<i class="fa-solid fa-magnifying-glass"></i>
<span class="search-button__default-text">Search</span>
<span class="search-button__kbd-shortcut"><kbd class="kbd-shortcut__modifier">Ctrl</kbd>+<kbd class="kbd-shortcut__modifier">K</kbd></span>
</button>
`);
</script></div>
<div class="sidebar-primary-item"><nav class="bd-links bd-docs-nav" aria-label="Main">
<div class="bd-toc-item navbar-nav active">
<ul class="nav bd-sidenav bd-sidenav__home-link">
<li class="toctree-l1">
<a class="reference internal" href="intro.html">
Applied Data Analysis and Machine Learning
</a>
</li>
</ul>
<p aria-level="2" class="caption" role="heading"><span class="caption-text">About the course</span></p>
<ul class="nav bd-sidenav">
<li class="toctree-l1"><a class="reference internal" href="schedule.html">Course setting</a></li>
<li class="toctree-l1"><a class="reference internal" href="teachers.html">Teachers and Grading</a></li>
<li class="toctree-l1"><a class="reference internal" href="textbooks.html">Textbooks</a></li>
</ul>
<p aria-level="2" class="caption" role="heading"><span class="caption-text">Review of Statistics with Resampling Techniques and Linear Algebra</span></p>
<ul class="nav bd-sidenav">
<li class="toctree-l1"><a class="reference internal" href="statistics.html">1. Elements of Probability Theory and Statistical Data Analysis</a></li>
<li class="toctree-l1"><a class="reference internal" href="linalg.html">2. Linear Algebra, Handling of Arrays and more Python Features</a></li>
</ul>
<p aria-level="2" class="caption" role="heading"><span class="caption-text">From Regression to Support Vector Machines</span></p>
<ul class="nav bd-sidenav">
<li class="toctree-l1"><a class="reference internal" href="chapter1.html">3. Linear Regression</a></li>
<li class="toctree-l1"><a class="reference internal" href="chapter2.html">4. Ridge and Lasso Regression</a></li>
<li class="toctree-l1"><a class="reference internal" href="chapter3.html">5. Resampling Methods</a></li>
<li class="toctree-l1"><a class="reference internal" href="chapter4.html">6. Logistic Regression</a></li>
<li class="toctree-l1"><a class="reference internal" href="chapteroptimization.html">7. Optimization, the central part of any Machine Learning algortithm</a></li>
<li class="toctree-l1"><a class="reference internal" href="chapter5.html">8. Support Vector Machines, overarching aims</a></li>
</ul>
<p aria-level="2" class="caption" role="heading"><span class="caption-text">Decision Trees, Ensemble Methods and Boosting</span></p>
<ul class="nav bd-sidenav">
<li class="toctree-l1"><a class="reference internal" href="chapter6.html">9. Decision trees, overarching aims</a></li>
<li class="toctree-l1"><a class="reference internal" href="chapter7.html">10. Ensemble Methods: From a Single Tree to Many Trees and Extreme Boosting, Meet the Jungle of Methods</a></li>
</ul>
<p aria-level="2" class="caption" role="heading"><span class="caption-text">Dimensionality Reduction</span></p>
<ul class="nav bd-sidenav">
<li class="toctree-l1"><a class="reference internal" href="chapter8.html">11. Basic ideas of the Principal Component Analysis (PCA)</a></li>
<li class="toctree-l1"><a class="reference internal" href="clustering.html">12. Clustering and Unsupervised Learning</a></li>
</ul>
<p aria-level="2" class="caption" role="heading"><span class="caption-text">Deep Learning Methods</span></p>
<ul class="nav bd-sidenav">
<li class="toctree-l1"><a class="reference internal" href="chapter9.html">13. Neural networks</a></li>
<li class="toctree-l1"><a class="reference internal" href="chapter10.html">14. Building a Feed Forward Neural Network</a></li>
<li class="toctree-l1"><a class="reference internal" href="chapter11.html">15. Solving Differential Equations with Deep Learning</a></li>
<li class="toctree-l1"><a class="reference internal" href="chapter12.html">16. Convolutional Neural Networks</a></li>
<li class="toctree-l1"><a class="reference internal" href="chapter13.html">17. Recurrent neural networks: Overarching view</a></li>
</ul>
<p aria-level="2" class="caption" role="heading"><span class="caption-text">Weekly material, notes and exercises</span></p>
<ul class="current nav bd-sidenav">
<li class="toctree-l1"><a class="reference internal" href="exercisesweek34.html">Exercises week 34</a></li>
<li class="toctree-l1"><a class="reference internal" href="week34.html">Week 34: Introduction to the course, Logistics and Practicalities</a></li>
<li class="toctree-l1"><a class="reference internal" href="exercisesweek35.html">Exercises week 35</a></li>
<li class="toctree-l1"><a class="reference internal" href="week35.html">Week 35: From Ordinary Linear Regression to Ridge and Lasso Regression</a></li>
<li class="toctree-l1"><a class="reference internal" href="exercisesweek36.html">Exercises week 36</a></li>
<li class="toctree-l1 current active"><a class="current reference internal" href="#">Week 36: Linear Regression and Gradient descent</a></li>
<li class="toctree-l1"><a class="reference internal" href="exercisesweek37.html">Exercises week 37</a></li>
<li class="toctree-l1"><a class="reference internal" href="week37.html">Week 37: Gradient descent methods</a></li>
<li class="toctree-l1"><a class="reference internal" href="exercisesweek38.html">Exercises week 37</a></li>
</ul>
<p aria-level="2" class="caption" role="heading"><span class="caption-text">Projects</span></p>
<ul class="nav bd-sidenav">
<li class="toctree-l1"><a class="reference internal" href="project1.html">Project 1 on Machine Learning, deadline October 6 (midnight), 2025</a></li>
</ul>
</div>
</nav></div>
</div>
<div class="sidebar-primary-items__end sidebar-primary__section">
</div>
<div id="rtd-footer-container"></div>
</div>
<main id="main-content" class="bd-main" role="main">
<div class="sbt-scroll-pixel-helper"></div>
<div class="bd-content">
<div class="bd-article-container">
<div class="bd-header-article d-print-none">
<div class="header-article-items header-article__inner">
<div class="header-article-items__start">
<div class="header-article-item"><button class="sidebar-toggle primary-toggle btn btn-sm" title="Toggle primary sidebar" data-bs-placement="bottom" data-bs-toggle="tooltip">
<span class="fa-solid fa-bars"></span>
</button></div>
</div>
<div class="header-article-items__end">
<div class="header-article-item">
<div class="article-header-buttons">
<div class="dropdown dropdown-download-buttons">
<button class="btn dropdown-toggle" type="button" data-bs-toggle="dropdown" aria-expanded="false" aria-label="Download this page">
<i class="fas fa-download"></i>
</button>
<ul class="dropdown-menu">
<li><a href="_sources/week36.ipynb" target="_blank"
class="btn btn-sm btn-download-source-button dropdown-item"
title="Download source file"
data-bs-placement="left" data-bs-toggle="tooltip"
>
<span class="btn__icon-container">
<i class="fas fa-file"></i>
</span>
<span class="btn__text-container">.ipynb</span>
</a>
</li>
<li>
<button onclick="window.print()"
class="btn btn-sm btn-download-pdf-button dropdown-item"
title="Print to PDF"
data-bs-placement="left" data-bs-toggle="tooltip"
>
<span class="btn__icon-container">
<i class="fas fa-file-pdf"></i>
</span>
<span class="btn__text-container">.pdf</span>
</button>
</li>
</ul>
</div>
<button onclick="toggleFullScreen()"
class="btn btn-sm btn-fullscreen-button"
title="Fullscreen mode"
data-bs-placement="bottom" data-bs-toggle="tooltip"
>
<span class="btn__icon-container">
<i class="fas fa-expand"></i>
</span>
</button>
<script>
document.write(`
<button class="btn btn-sm nav-link pst-navbar-icon theme-switch-button" title="light/dark" aria-label="light/dark" data-bs-placement="bottom" data-bs-toggle="tooltip">
<i class="theme-switch fa-solid fa-sun fa-lg" data-mode="light"></i>
<i class="theme-switch fa-solid fa-moon fa-lg" data-mode="dark"></i>
<i class="theme-switch fa-solid fa-circle-half-stroke fa-lg" data-mode="auto"></i>
</button>
`);
</script>
<script>
document.write(`
<button class="btn btn-sm pst-navbar-icon search-button search-button__button" title="Search" aria-label="Search" data-bs-placement="bottom" data-bs-toggle="tooltip">
<i class="fa-solid fa-magnifying-glass fa-lg"></i>
</button>
`);
</script>
<button class="sidebar-toggle secondary-toggle btn btn-sm" title="Toggle secondary sidebar" data-bs-placement="bottom" data-bs-toggle="tooltip">
<span class="fa-solid fa-list"></span>
</button>
</div></div>
</div>
</div>
</div>
<div id="jb-print-docs-body" class="onlyprint">
<h1>Week 36: Linear Regression and Gradient descent</h1>
<!-- Table of contents -->
<div id="print-main-content">
<div id="jb-print-toc">
<div>
<h2> Contents </h2>
</div>
<nav aria-label="Page">
<ul class="visible nav section-nav flex-column">
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#plans-for-week-36">Plans for week 36</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#material-for-lecture-monday-september-2">Material for lecture Monday September 2</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#mathematical-interpretation-of-ordinary-least-squares">Mathematical Interpretation of Ordinary Least Squares</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#residual-error">Residual Error</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#simple-case">Simple case</a></li>
<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>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#codes-for-the-svd">Codes for the SVD</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#note-about-svd-calculations">Note about SVD Calculations</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#mathematics-of-the-svd-and-implications">Mathematics of the SVD and implications</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#example-matrix">Example Matrix</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#setting-up-the-matrix-to-be-inverted">Setting up the Matrix to be inverted</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#further-properties-important-for-our-analyses-later">Further properties (important for our analyses later)</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>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#optimization-and-gradient-descent-the-central-part-of-any-machine-learning-algortithm">Optimization and gradient descent, the central part of any Machine Learning algortithm</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#reminder-on-newton-raphson-s-method">Reminder on Newton-Raphson’s method</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-equations">The equations</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#simple-geometric-interpretation">Simple geometric interpretation</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#extending-to-more-than-one-variable">Extending to more than one variable</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#steepest-descent">Steepest descent</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#more-on-steepest-descent">More on Steepest descent</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-ideal">The ideal</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-sensitiveness-of-the-gradient-descent">The sensitiveness of the gradient descent</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#convex-functions">Convex functions</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#convex-function">Convex function</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#conditions-on-convex-functions">Conditions on convex functions</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#more-on-convex-functions">More on convex functions</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#some-simple-problems">Some simple problems</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#revisiting-ordinary-least-squares">Revisiting Ordinary Least Squares</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#gradient-descent-example">Gradient descent example</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-derivative-of-the-cost-loss-function">The derivative of the cost/loss function</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-hessian-matrix">The Hessian matrix</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#simple-program">Simple program</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#id1">Gradient Descent Example</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#gradient-descent-and-ridge">Gradient descent and Ridge</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-hessian-matrix-for-ridge-regression">The Hessian matrix for Ridge Regression</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#program-example-for-gradient-descent-with-ridge-regression">Program example for gradient descent with Ridge Regression</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#using-gradient-descent-methods-limitations">Using gradient descent methods, limitations</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#material-for-lab-sessions-sessions-tuesday-and-wednesday">Material for lab sessions sessions Tuesday and Wednesday</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#linear-regression-and-the-svd">Linear Regression and 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="#id2">Ridge and LASSO Regression</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#from-ols-to-ridge-and-lasso">From OLS to Ridge and Lasso</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#id3">Deriving the Ridge Regression Equations</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#note-on-scikit-learn">Note on Scikit-Learn</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#comparison-with-ols">Comparison with OLS</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#svd-analysis">SVD analysis</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#id4">Interpreting the Ridge results</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#id5">More interpretations</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#id6">Deriving the Lasso Regression Equations</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#simple-example-to-illustrate-ordinary-least-squares-ridge-and-lasso-regression">Simple example to illustrate Ordinary Least Squares, Ridge and Lasso Regression</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#ridge-regression">Ridge Regression</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#lasso-regression">Lasso Regression</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#yet-another-example">Yet another Example</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-ols-case">The OLS case</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-ridge-case">The Ridge case</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#writing-the-cost-function">Writing the Cost Function</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#lasso-case">Lasso case</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-first-case">The first Case</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#simple-code-for-solving-the-above-problem">Simple code for solving the above problem</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#with-lasso-regression">With Lasso Regression</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#another-example-now-with-a-polynomial-fit">Another Example, now with a polynomial fit</a></li>
</ul>
</nav>
</div>
</div>
</div>
<div id="searchbox"></div>
<article class="bd-article">
<!-- HTML file automatically generated from DocOnce source (https://github.com/doconce/doconce/)
doconce format html week36.do.txt --no_mako -->
<!-- dom:TITLE: Week 36: Linear Regression and Gradient descent --><section class="tex2jax_ignore mathjax_ignore" id="week-36-linear-regression-and-gradient-descent">
<h1>Week 36: Linear Regression and Gradient descent<a class="headerlink" href="#week-36-linear-regression-and-gradient-descent" title="Link to this heading">#</a></h1>
<p><strong>Morten Hjorth-Jensen</strong>, Department of Physics, University of Oslo, Norway</p>
<p>Date: <strong>September 1-5, 2025</strong></p>
<section id="plans-for-week-36">
<h2>Plans for week 36<a class="headerlink" href="#plans-for-week-36" title="Link to this heading">#</a></h2>
<p><strong>Material for the lecture on Monday September 1:</strong></p>
<ol class="arabic simple">
<li><p>Linear Regression, ordinary least squares (OLS), Ridge and Lasso and mathematical analysis</p></li>
<li><p>Derivation of Gradient descent and discussion of implementations for</p></li>
<li><p>Video of lecture at <a class="reference external" href="https://youtu.be/nVE_FRnGAHw">https://youtu.be/nVE_FRnGAHw</a></p></li>
<li><p>Whiteboard notes at <a class="github reference external" href="https://github.com/CompPhysics/MachineLearning/blob/master/doc/HandWrittenNotes/2025/FYSSTKweek36.pdf">CompPhysics/MachineLearning</a></p></li>
</ol>
<p><strong>Material for the lab sessions on Tuesday and Wednesday (see at the end of these slides):</strong></p>
<ol class="arabic simple">
<li><p>Technicalities concerning Ridge and Lasso linear regression.</p></li>
<li><p>Presentation and discussion of the first project</p></li>
</ol>
<!-- * [Video of lab session](https://youtu.be/ZrIdZdZtHe0) -->
<p><strong>Reading suggestion:</strong></p>
<ol class="arabic simple">
<li><p>Goodfellow et al, Deep Learning, introduction to gradient descent, see chapter 4.3 at <a class="reference external" href="https://www.deeplearningbook.org/contents/numerical.html">https://www.deeplearningbook.org/contents/numerical.html</a></p></li>
<li><p>Rashcka et al, pages 37-44 and pages 278-283 with focus on linear regression.</p></li>
<li><p>Video on gradient descent at <a class="reference external" href="https://www.youtube.com/watch?v=sDv4f4s2SB8">https://www.youtube.com/watch?v=sDv4f4s2SB8</a></p></li>
</ol>
</section>
<section id="material-for-lecture-monday-september-2">
<h2>Material for lecture Monday September 2<a class="headerlink" href="#material-for-lecture-monday-september-2" title="Link to this heading">#</a></h2>
</section>
<section id="mathematical-interpretation-of-ordinary-least-squares">
<h2>Mathematical Interpretation of Ordinary Least Squares<a class="headerlink" href="#mathematical-interpretation-of-ordinary-least-squares" title="Link to this heading">#</a></h2>
<p>What is presented here is a mathematical analysis of various regression algorithms (ordinary least squares, Ridge and Lasso Regression). The analysis is based on an important algorithm in linear algebra, the so-called Singular Value Decomposition (SVD).</p>
<p>We have shown that in ordinary least squares the optimal parameters <span class="math notranslate nohighlight">\(\theta\)</span> are given by</p>
<div class="math notranslate nohighlight">
\[
\hat{\boldsymbol{\theta}} = \left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}.
\]</div>
<p>The <strong>hat</strong> over <span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span> means we have the optimal parameters after minimization of the cost function.</p>
<p>This means that our best model is defined as</p>
<div class="math notranslate nohighlight">
\[
\tilde{\boldsymbol{y}}=\boldsymbol{X}\hat{\boldsymbol{\theta}} = \boldsymbol{X}\left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}.
\]</div>
<p>We now define a matrix</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{A}=\boldsymbol{X}\left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T.
\]</div>
<p>We can rewrite</p>
<div class="math notranslate nohighlight">
\[
\tilde{\boldsymbol{y}}=\boldsymbol{X}\hat{\boldsymbol{\theta}} = \boldsymbol{A}\boldsymbol{y}.
\]</div>
<p>The matrix <span class="math notranslate nohighlight">\(\boldsymbol{A}\)</span> has the important property that <span class="math notranslate nohighlight">\(\boldsymbol{A}^2=\boldsymbol{A}\)</span>. This is the definition of a projection matrix.
We can then interpret our optimal model <span class="math notranslate nohighlight">\(\tilde{\boldsymbol{y}}\)</span> as being represented by an orthogonal projection of <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> onto a space defined by the column vectors of <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>. In our case here the matrix <span class="math notranslate nohighlight">\(\boldsymbol{A}\)</span> is a square matrix. If it is a general rectangular matrix we have an oblique projection matrix.</p>
</section>
<section id="residual-error">
<h2>Residual Error<a class="headerlink" href="#residual-error" title="Link to this heading">#</a></h2>
<p>We have defined the residual error as</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{\epsilon}=\boldsymbol{y}-\tilde{\boldsymbol{y}}=\left[\boldsymbol{I}-\boldsymbol{X}\left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\right]\boldsymbol{y}.
\]</div>
<p>The residual errors are then the projections of <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> onto the orthogonal component of the space defined by the column vectors of <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>.</p>
</section>
<section id="simple-case">
<h2>Simple case<a class="headerlink" href="#simple-case" title="Link to this heading">#</a></h2>
<p>If the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> is an orthogonal (or unitary in case of complex values) matrix, we have</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{X}^T\boldsymbol{X}=\boldsymbol{X}\boldsymbol{X}^T = \boldsymbol{I}.
\]</div>
<p>In this case the matrix <span class="math notranslate nohighlight">\(\boldsymbol{A}\)</span> becomes</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{A}=\boldsymbol{X}\left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T)=\boldsymbol{I},
\]</div>
<p>and we have the obvious case</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{\epsilon}=\boldsymbol{y}-\tilde{\boldsymbol{y}}=0.
\]</div>
<p>This serves also as a useful test of our codes.</p>
</section>
<section id="the-singular-value-decomposition">
<h2>The singular value decomposition<a class="headerlink" href="#the-singular-value-decomposition" title="Link to this heading">#</a></h2>
<p>The examples we have looked at so far are cases where we normally can
invert the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span>. Using a polynomial expansion where we fit of various functions leads to
row vectors of the design matrix which are essentially orthogonal due
to the polynomial character of our model. Obtaining the inverse of the
design matrix is then often done via a so-called LU, QR or Cholesky
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 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>
<p>This is given by the <strong>Singular Value Decomposition</strong> (SVD) algorithm,
perhaps the most powerful linear algebra algorithm. The SVD provides
a numerically stable matrix decomposition that is used in a large
swath oc applications and the decomposition is always stable
numerically.</p>
<p>In machine learning it plays a central role in dealing with for
example design matrices that may be near singular or singular.
Furthermore, as we will see here, the singular values can be related
to the covariance matrix (and thereby the correlation matrix) and in
turn the variance of a given quantity. It plays also an important role
in the principal component analysis where high-dimensional data can be
reduced to the statistically relevant features.</p>
</section>
<section id="linear-regression-problems">
<h2>Linear Regression Problems<a class="headerlink" href="#linear-regression-problems" title="Link to this heading">#</a></h2>
<p>One of the typical problems we encounter with linear regression, in particular
when the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> (our so-called design matrix) is high-dimensional,
are problems with near singular or singular matrices. The column vectors of <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>
may be linearly dependent, normally referred to as super-collinearity.<br />
This means that the matrix may be rank deficient and it is basically impossible to
to model the data using linear regression. As an example, consider the matrix</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{align*}
\mathbf{X} &amp; = \left[
\begin{array}{rrr}
1 &amp; -1 &amp; 2
\\
1 &amp; 0 &amp; 1
\\
1 &amp; 2 &amp; -1
\\
1 &amp; 1 &amp; 0
\end{array} \right]
\end{align*}
\end{split}\]</div>
<p>The columns of <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> are linearly dependent. We see this easily since the
the first column is the row-wise sum of the other two columns. The rank (more correct,
the column rank) of a matrix is the dimension of the space spanned by the
column vectors. Hence, the rank of <span class="math notranslate nohighlight">\(\mathbf{X}\)</span> is equal to the number
of linearly independent columns. In this particular case the matrix has rank 2.</p>
<p>Super-collinearity of an <span class="math notranslate nohighlight">\((n \times p)\)</span>-dimensional design matrix <span class="math notranslate nohighlight">\(\mathbf{X}\)</span> implies
that the inverse of the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span> (the matrix we need to invert to solve the linear regression equations) is non-invertible. If we have a square matrix that does not have an inverse, we say this matrix singular. The example here demonstrates this</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{align*}
\boldsymbol{X} &amp; = \left[
\begin{array}{rr}
1 &amp; -1
\\
1 &amp; -1
\end{array} \right].
\end{align*}
\end{split}\]</div>
<p>We see easily that <span class="math notranslate nohighlight">\(\mbox{det}(\boldsymbol{X}) = x_{11} x_{22} - x_{12} x_{21} = 1 \times (-1) - 1 \times (-1) = 0\)</span>. Hence, <span class="math notranslate nohighlight">\(\mathbf{X}\)</span> is singular and its inverse is undefined.
This is equivalent to saying that the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> has at least an eigenvalue which is zero.</p>
</section>
<section id="fixing-the-singularity">
<h2>Fixing the singularity<a class="headerlink" href="#fixing-the-singularity" title="Link to this heading">#</a></h2>
<p>If our design matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> which enters the linear regression problem</p>
<!-- Equation labels as ordinary links -->
<div id="_auto1"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\boldsymbol{\theta} = (\boldsymbol{X}^{T} \boldsymbol{X})^{-1} \boldsymbol{X}^{T} \boldsymbol{y},
\label{_auto1} \tag{1}
\end{equation}
\]</div>
<p>has linearly dependent column vectors, we will not be able to compute the inverse
of <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span> and we cannot find the parameters (estimators) <span class="math notranslate nohighlight">\(\theta_i\)</span>.
The estimators are only well-defined if <span class="math notranslate nohighlight">\((\boldsymbol{X}^{T}\boldsymbol{X})^{-1}\)</span> exists.
This is more likely to happen when the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> is high-dimensional. In this case it is likely to encounter a situation where
the regression parameters <span class="math notranslate nohighlight">\(\theta_i\)</span> cannot be estimated.</p>
<p>A cheap <em>ad hoc</em> approach is simply to add a small diagonal component to the matrix to invert, that is we change</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{X}^{T} \boldsymbol{X} \rightarrow \boldsymbol{X}^{T} \boldsymbol{X}+\lambda \boldsymbol{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>%matplotlib inline
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(&quot;Theta values for own Ridge implementation&quot;)
print(OwnRidgeTheta)
print(&quot;Theta values for Scikit-Learn Ridge implementation&quot;)
print(RegRidge.coef_)
print(&quot;MSE values for own Ridge implementation&quot;)
print(MSEOwnRidgePredict[i])
print(&quot;MSE values for Scikit-Learn Ridge implementation&quot;)
print(MSERidgePredict[i])
# Now plot the results
plt.figure()
plt.plot(np.log10(lambdas), MSEOwnRidgePredict, &#39;r&#39;, label = &#39;MSE own Ridge Test&#39;)
plt.plot(np.log10(lambdas), MSERidgePredict, &#39;g&#39;, label = &#39;MSE Ridge Test&#39;)
plt.xlabel(&#39;log10(lambda)&#39;)
plt.ylabel(&#39;MSE&#39;)
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 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>
<div class="math notranslate nohighlight">
\[
(\lambda_1,\boldsymbol{u}_1),\dots, (\lambda_n,\boldsymbol{u}_n),
\]</div>
<p>and the eigenvalues are given by the diagonal matrix</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{\Sigma}=\mathrm{Diag}(\lambda_1, \dots,\lambda_n).
\]</div>
<p>The matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> can be written in terms of an orthogonal/unitary transformation <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span></p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{X} = \boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^T,
\]</div>
<p>with <span class="math notranslate nohighlight">\(\boldsymbol{U}\boldsymbol{U}^T=\boldsymbol{I}\)</span> or <span class="math notranslate nohighlight">\(\boldsymbol{U}\boldsymbol{U}^{\dagger}=\boldsymbol{I}\)</span>.</p>
<p>Not all square matrices are diagonalizable. A matrix like the one discussed above</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\boldsymbol{X} = \begin{bmatrix}
1&amp; -1 \\
1&amp; -1\\
\end{bmatrix}
\end{split}\]</div>
<p>is not diagonalizable, it is a so-called <a class="reference external" href="https://en.wikipedia.org/wiki/Defective_matrix">defective matrix</a>. It is easy to see that the condition
<span class="math notranslate nohighlight">\(\boldsymbol{X}\boldsymbol{X}^T=\boldsymbol{X}^T\boldsymbol{X}\)</span> is not fulfilled.</p>
</section>
<section id="the-svd-a-fantastic-algorithm">
<h2>The SVD, a Fantastic Algorithm<a class="headerlink" href="#the-svd-a-fantastic-algorithm" title="Link to this heading">#</a></h2>
<p>However, and this is the strength of the SVD algorithm, any general
matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> can be decomposed in terms of a diagonal matrix and
two orthogonal/unitary matrices. The <a class="reference external" href="https://en.wikipedia.org/wiki/Singular_value_decomposition">Singular Value Decompostion
(SVD) theorem</a>
states that a general <span class="math notranslate nohighlight">\(m\times n\)</span> matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> can be written in
terms of a diagonal matrix <span class="math notranslate nohighlight">\(\boldsymbol{\Sigma}\)</span> of dimensionality <span class="math notranslate nohighlight">\(m\times n\)</span>
and two orthognal matrices <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> and <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span>, where the first has
dimensionality <span class="math notranslate nohighlight">\(m \times m\)</span> and the last dimensionality <span class="math notranslate nohighlight">\(n\times n\)</span>.
We have then</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{X} = \boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^T
\]</div>
<p>As an example, the above defective matrix can be decomposed as</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\boldsymbol{X} = \frac{1}{\sqrt{2}}\begin{bmatrix} 1&amp; 1 \\ 1&amp; -1\\ \end{bmatrix} \begin{bmatrix} 2&amp; 0 \\ 0&amp; 0\\ \end{bmatrix} \frac{1}{\sqrt{2}}\begin{bmatrix} 1&amp; -1 \\ 1&amp; 1\\ \end{bmatrix}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^T,
\end{split}\]</div>
<p>with eigenvalues <span class="math notranslate nohighlight">\(\sigma_1=2\)</span> and <span class="math notranslate nohighlight">\(\sigma_2=0\)</span>.
The SVD exits always!</p>
<p>The SVD
decomposition (singular values) gives eigenvalues
<span class="math notranslate nohighlight">\(\sigma_i\geq\sigma_{i+1}\)</span> for all <span class="math notranslate nohighlight">\(i\)</span> and for dimensions larger than <span class="math notranslate nohighlight">\(i=p\)</span>, the
eigenvalues (singular values) are zero.</p>
<p>In the general case, where our design matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> has dimension
<span class="math notranslate nohighlight">\(n\times p\)</span>, the matrix is thus decomposed into an <span class="math notranslate nohighlight">\(n\times n\)</span>
orthogonal matrix <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span>, a <span class="math notranslate nohighlight">\(p\times p\)</span> orthogonal matrix <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span>
and a diagonal matrix <span class="math notranslate nohighlight">\(\boldsymbol{\Sigma}\)</span> with <span class="math notranslate nohighlight">\(r=\mathrm{min}(n,p)\)</span>
singular values <span class="math notranslate nohighlight">\(\sigma_i\geq 0\)</span> on the main diagonal and zeros filling
the rest of the matrix. There are at most <span class="math notranslate nohighlight">\(p\)</span> singular values
assuming that <span class="math notranslate nohighlight">\(n &gt; p\)</span>. In our regression examples for the nuclear
masses and the equation of state this is indeed the case, while for
the Ising model we have <span class="math notranslate nohighlight">\(p &gt; n\)</span>. These are often cases that lead to
near singular or singular matrices.</p>
<p>The columns of <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> are called the left singular vectors while the columns of <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span> are the right singular vectors.</p>
</section>
<section id="economy-size-svd">
<h2>Economy-size SVD<a class="headerlink" href="#economy-size-svd" title="Link to this heading">#</a></h2>
<p>If we assume that <span class="math notranslate nohighlight">\(n &gt; p\)</span>, then our matrix <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> has dimension <span class="math notranslate nohighlight">\(n
\times n\)</span>. The last <span class="math notranslate nohighlight">\(n-p\)</span> columns of <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> become however
irrelevant in our calculations since they are multiplied with the
zeros in <span class="math notranslate nohighlight">\(\boldsymbol{\Sigma}\)</span>.</p>
<p>The economy-size decomposition removes extra rows or columns of zeros
from the diagonal matrix of singular values, <span class="math notranslate nohighlight">\(\boldsymbol{\Sigma}\)</span>, along with the columns
in either <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> or <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span> that multiply those zeros in the expression.
Removing these zeros and columns can improve execution time
and reduce storage requirements without compromising the accuracy of
the decomposition.</p>
<p>If <span class="math notranslate nohighlight">\(n &gt; p\)</span>, we keep only the first <span class="math notranslate nohighlight">\(p\)</span> columns of <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> and <span class="math notranslate nohighlight">\(\boldsymbol{\Sigma}\)</span> has dimension <span class="math notranslate nohighlight">\(p\times p\)</span>.
If <span class="math notranslate nohighlight">\(p &gt; n\)</span>, then only the first <span class="math notranslate nohighlight">\(n\)</span> columns of <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span> are computed and <span class="math notranslate nohighlight">\(\boldsymbol{\Sigma}\)</span> has dimension <span class="math notranslate nohighlight">\(n\times n\)</span>.
The <span class="math notranslate nohighlight">\(n=p\)</span> case is obvious, we retain the full SVD.
In general the economy-size SVD leads to less FLOPS and still conserving the desired accuracy.</p>
</section>
<section id="codes-for-the-svd">
<h2>Codes for the SVD<a class="headerlink" href="#codes-for-the-svd" 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>import numpy as np
# SVD inversion
def SVD(A):
&#39;&#39;&#39; Takes as input a numpy matrix A and returns inv(A) based on singular value decomposition (SVD).
SVD is numerically more stable than the inversion algorithms provided by
numpy and scipy.linalg at the cost of being slower.
&#39;&#39;&#39;
U, S, VT = np.linalg.svd(A,full_matrices=True)
print(&#39;test U&#39;)
print( (np.transpose(U) @ U - U @np.transpose(U)))
print(&#39;test VT&#39;)
print( (np.transpose(VT) @ VT - VT @np.transpose(VT)))
print(U)
print(S)
print(VT)
D = np.zeros((len(U),len(VT)))
for i in range(0,len(VT)):
D[i,i]=S[i]
return U @ D @ VT
X = np.array([ [1.0,-1.0], [1.0,-1.0]])
#X = np.array([[1, 2], [3, 4], [5, 6]])
print(X)
C = SVD(X)
# Print the difference between the original matrix and the SVD one
print(C-X)
</pre></div>
</div>
</div>
</div>
<p>The matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> has columns that are linearly dependent. The first
column is the row-wise sum of the other two columns. The rank of a
matrix (the column rank) is the dimension of space spanned by the
column vectors. The rank of the matrix is the number of linearly
independent columns, in this case just <span class="math notranslate nohighlight">\(2\)</span>. We see this from the
singular values when running the above code. Running the standard
inversion algorithm for matrix inversion with <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span> results
in the program terminating due to a singular matrix.</p>
</section>
<section id="note-about-svd-calculations">
<h2>Note about SVD Calculations<a class="headerlink" href="#note-about-svd-calculations" title="Link to this heading">#</a></h2>
<p>The <span class="math notranslate nohighlight">\(U\)</span>, <span class="math notranslate nohighlight">\(S\)</span>, and <span class="math notranslate nohighlight">\(V\)</span> matrices returned from the <strong>svd()</strong> function
cannot be multiplied directly.</p>
<p>As you can see from the code, the <span class="math notranslate nohighlight">\(S\)</span> vector must be converted into a
diagonal matrix. This may cause a problem as the size of the matrices
do not fit the rules of matrix multiplication, where the number of
columns in a matrix must match the number of rows in the subsequent
matrix.</p>
<p>If you wish to include the zero singular values, you will need to
resize the matrices and set up a diagonal matrix as done in the above
example</p>
</section>
<section id="mathematics-of-the-svd-and-implications">
<h2>Mathematics of the SVD and implications<a class="headerlink" href="#mathematics-of-the-svd-and-implications" title="Link to this heading">#</a></h2>
<p>Let us take a closer look at the mathematics of the SVD and the various implications for machine learning studies.</p>
<p>Our starting point is our design matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> of dimension <span class="math notranslate nohighlight">\(n\times p\)</span></p>
<div class="math notranslate nohighlight">
\[\begin{split}
\boldsymbol{X}=\begin{bmatrix}
x_{0,0} &amp; x_{0,1} &amp; x_{0,2}&amp; \dots &amp; \dots x_{0,p-1}\\
x_{1,0} &amp; x_{1,1} &amp; x_{1,2}&amp; \dots &amp; \dots x_{1,p-1}\\
x_{2,0} &amp; x_{2,1} &amp; x_{2,2}&amp; \dots &amp; \dots x_{2,p-1}\\
\dots &amp; \dots &amp; \dots &amp; \dots \dots &amp; \dots \\
x_{n-2,0} &amp; x_{n-2,1} &amp; x_{n-2,2}&amp; \dots &amp; \dots x_{n-2,p-1}\\
x_{n-1,0} &amp; x_{n-1,1} &amp; x_{n-1,2}&amp; \dots &amp; \dots x_{n-1,p-1}\\
\end{bmatrix}.
\end{split}\]</div>
<p>We can SVD decompose our matrix as</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{X}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^T,
\]</div>
<p>where <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> is an orthogonal matrix of dimension <span class="math notranslate nohighlight">\(n\times n\)</span>, meaning that <span class="math notranslate nohighlight">\(\boldsymbol{U}\boldsymbol{U}^T=\boldsymbol{U}^T\boldsymbol{U}=\boldsymbol{I}_n\)</span>. Here <span class="math notranslate nohighlight">\(\boldsymbol{I}_n\)</span> is the unit matrix of dimension <span class="math notranslate nohighlight">\(n \times n\)</span>.</p>
<p>Similarly, <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span> is an orthogonal matrix of dimension <span class="math notranslate nohighlight">\(p\times p\)</span>, meaning that <span class="math notranslate nohighlight">\(\boldsymbol{V}\boldsymbol{V}^T=\boldsymbol{V}^T\boldsymbol{V}=\boldsymbol{I}_p\)</span>. Here <span class="math notranslate nohighlight">\(\boldsymbol{I}_n\)</span> is the unit matrix of dimension <span class="math notranslate nohighlight">\(p \times p\)</span>.</p>
<p>Finally <span class="math notranslate nohighlight">\(\boldsymbol{\Sigma}\)</span> contains the singular values <span class="math notranslate nohighlight">\(\sigma_i\)</span>. This matrix has dimension <span class="math notranslate nohighlight">\(n\times p\)</span> and the singular values <span class="math notranslate nohighlight">\(\sigma_i\)</span> are all positive. The non-zero values are ordered in descending order, that is</p>
<div class="math notranslate nohighlight">
\[
\sigma_0 &gt; \sigma_1 &gt; \sigma_2 &gt; \dots &gt; \sigma_{p-1} &gt; 0.
\]</div>
<p>All values beyond <span class="math notranslate nohighlight">\(p-1\)</span> are all zero.</p>
</section>
<section id="example-matrix">
<h2>Example Matrix<a class="headerlink" href="#example-matrix" title="Link to this heading">#</a></h2>
<p>As an example, consider the following <span class="math notranslate nohighlight">\(3\times 2\)</span> example for the matrix <span class="math notranslate nohighlight">\(\boldsymbol{\Sigma}\)</span></p>
<div class="math notranslate nohighlight">
\[\begin{split}
\boldsymbol{\Sigma}=
\begin{bmatrix}
2&amp; 0 \\
0 &amp; 1 \\
0 &amp; 0 \\
\end{bmatrix}
\end{split}\]</div>
<p>The singular values are <span class="math notranslate nohighlight">\(\sigma_0=2\)</span> and <span class="math notranslate nohighlight">\(\sigma_1=1\)</span>. It is common to rewrite the matrix <span class="math notranslate nohighlight">\(\boldsymbol{\Sigma}\)</span> as</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\boldsymbol{\Sigma}=
\begin{bmatrix}
\boldsymbol{\tilde{\Sigma}}\\
\boldsymbol{0}\\
\end{bmatrix},
\end{split}\]</div>
<p>where</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\boldsymbol{\tilde{\Sigma}}=
\begin{bmatrix}
2&amp; 0 \\
0 &amp; 1 \\
\end{bmatrix},
\end{split}\]</div>
<p>contains only the singular values. Note also (and we will use this below) that</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\boldsymbol{\Sigma}^T\boldsymbol{\Sigma}=
\begin{bmatrix}
4&amp; 0 \\
0 &amp; 1 \\
\end{bmatrix},
\end{split}\]</div>
<p>which is a <span class="math notranslate nohighlight">\(2\times 2 \)</span> matrix while</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\boldsymbol{\Sigma}\boldsymbol{\Sigma}^T=
\begin{bmatrix}
4&amp; 0 &amp; 0\\
0 &amp; 1 &amp; 0\\
0 &amp; 0 &amp; 0\\
\end{bmatrix},
\end{split}\]</div>
<p>is a <span class="math notranslate nohighlight">\(3\times 3 \)</span> matrix. The last row and column of this last matrix
contain only zeros. This will have important consequences for our SVD
decomposition of the design matrix.</p>
</section>
<section id="setting-up-the-matrix-to-be-inverted">
<h2>Setting up the Matrix to be inverted<a class="headerlink" href="#setting-up-the-matrix-to-be-inverted" title="Link to this heading">#</a></h2>
<p>The matrix that may cause problems for us is <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span>. Using the SVD we can rewrite this matrix as</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{X}^T\boldsymbol{X}=\boldsymbol{V}\boldsymbol{\Sigma}^T\boldsymbol{U}^T\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^T,
\]</div>
<p>and using the orthogonality of the matrix <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> we have</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{X}^T\boldsymbol{X}=\boldsymbol{V}\boldsymbol{\Sigma}^T\boldsymbol{\Sigma}\boldsymbol{V}^T.
\]</div>
<p>We define <span class="math notranslate nohighlight">\(\boldsymbol{\Sigma}^T\boldsymbol{\Sigma}=\tilde{\boldsymbol{\Sigma}}^2\)</span> which is a diagonal matrix containing only the singular values squared. It has dimensionality <span class="math notranslate nohighlight">\(p \times p\)</span>.</p>
<p>We can now insert the result for the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span> into our equation for ordinary least squares where</p>
<div class="math notranslate nohighlight">
\[
\tilde{y}_{\mathrm{OLS}}=\boldsymbol{X}\left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y},
\]</div>
<p>and using our SVD decomposition of <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> we have</p>
<div class="math notranslate nohighlight">
\[
\tilde{y}_{\mathrm{OLS}}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^T\left(\boldsymbol{V}\tilde{\boldsymbol{\Sigma}}^{2}(\boldsymbol{V}^T\right)^{-1}\boldsymbol{V}\boldsymbol{\Sigma}^T\boldsymbol{U}^T\boldsymbol{y},
\]</div>
<p>which gives us, using the orthogonality of the matrix <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span>,</p>
<div class="math notranslate nohighlight">
\[
\tilde{y}_{\mathrm{OLS}}=\sum_{i=0}^{p-1}\boldsymbol{u}_i\boldsymbol{u}^T_i\boldsymbol{y},
\]</div>
<p>which is not the same as <span class="math notranslate nohighlight">\(\tilde{y}_{\mathrm{OLS}}=\boldsymbol{U}\boldsymbol{U}^T\boldsymbol{y}\)</span>, which due to the orthogonality of <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> would have given us that the model equals the output.</p>
<p>It means that the ordinary least square model (with the optimal
parameters) <span class="math notranslate nohighlight">\(\boldsymbol{\tilde{y}}\)</span>, corresponds to an orthogonal
transformation of the output (or target) vector <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> by the
vectors of the matrix <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span>. <strong>Note that the summation ends at</strong>
<span class="math notranslate nohighlight">\(p-1\)</span>, that is <span class="math notranslate nohighlight">\(\boldsymbol{\tilde{y}}\ne \boldsymbol{y}\)</span>. We can thus not use the
orthogonality relation for the matrix <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span>.</p>
</section>
<section id="further-properties-important-for-our-analyses-later">
<h2>Further properties (important for our analyses later)<a class="headerlink" href="#further-properties-important-for-our-analyses-later" title="Link to this heading">#</a></h2>
<p>Let us study again <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span> in terms of our SVD,</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{X}^T\boldsymbol{X}=\boldsymbol{V}\boldsymbol{\Sigma}^T\boldsymbol{U}^T\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^T=\boldsymbol{V}\boldsymbol{\Sigma}^T\boldsymbol{\Sigma}\boldsymbol{V}^T.
\]</div>
<p>If we now multiply from the right with <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span> (using the orthogonality of <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span>) we get</p>
<div class="math notranslate nohighlight">
\[
\left(\boldsymbol{X}^T\boldsymbol{X}\right)\boldsymbol{V}=\boldsymbol{V}\boldsymbol{\Sigma}^T\boldsymbol{\Sigma}.
\]</div>
<p>This means the vectors <span class="math notranslate nohighlight">\(\boldsymbol{v}_i\)</span> of the orthogonal matrix <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span> are the eigenvectors of the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span>
with eigenvalues given by the singular values squared, that is</p>
<div class="math notranslate nohighlight">
\[
\left(\boldsymbol{X}^T\boldsymbol{X}\right)\boldsymbol{v}_i=\boldsymbol{v}_i\sigma_i^2.
\]</div>
<p>Similarly, if we use the SVD decomposition for the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\boldsymbol{X}^T\)</span>, we have</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{X}\boldsymbol{X}^T=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^T\boldsymbol{V}\boldsymbol{\Sigma}^T\boldsymbol{U}^T=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{\Sigma}^T\boldsymbol{U}^T.
\]</div>
<p>If we now multiply from the right with <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> (using the orthogonality of <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span>) we get</p>
<div class="math notranslate nohighlight">
\[
\left(\boldsymbol{X}\boldsymbol{X}^T\right)\boldsymbol{U}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{\Sigma}^T.
\]</div>
<p>This means the vectors <span class="math notranslate nohighlight">\(\boldsymbol{u}_i\)</span> of the orthogonal matrix <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> are the eigenvectors of the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\boldsymbol{X}^T\)</span>
with eigenvalues given by the singular values squared, that is</p>
<div class="math notranslate nohighlight">
\[
\left(\boldsymbol{X}\boldsymbol{X}^T\right)\boldsymbol{u}_i=\boldsymbol{u}_i\sigma_i^2.
\]</div>
<p><strong>Important note</strong>: we have defined our design matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> to be an
<span class="math notranslate nohighlight">\(n\times p\)</span> matrix. In most supervised learning cases we have that <span class="math notranslate nohighlight">\(n
\ge p\)</span>, and quite often we have <span class="math notranslate nohighlight">\(n &gt;&gt; p\)</span>. For linear algebra based methods like ordinary least squares or Ridge regression, this leads to a matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span> which is small and thereby easier to handle from a computational point of view (in terms of number of floating point operations).</p>
<p>In our lectures, the number of columns will
always refer to the number of features in our data set, while the
number of rows represents the number of data inputs. Note that in
other texts you may find the opposite notation. This has consequences
for the definition of for example the covariance matrix and its relation to the SVD.</p>
</section>
<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">
\[
{\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>
<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
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>Using our insights about the SVD of the design matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>
We have already analyzed the OLS solutions in terms of the eigenvectors (the columns) of the right singular value matrix <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> as</p>
<div class="math notranslate nohighlight">
\[
\tilde{\boldsymbol{y}}_{\mathrm{OLS}}=\boldsymbol{X}\boldsymbol{\theta} =\boldsymbol{U}\boldsymbol{U}^T\boldsymbol{y}.
\]</div>
<p>For Ridge regression this becomes</p>
<div class="math notranslate nohighlight">
\[
\tilde{\boldsymbol{y}}_{\mathrm{Ridge}}=\boldsymbol{X}\boldsymbol{\theta}_{\mathrm{Ridge}} = \boldsymbol{U\Sigma V^T}\left(\boldsymbol{V}\boldsymbol{\Sigma}^2\boldsymbol{V}^T+\lambda\boldsymbol{I} \right)^{-1}(\boldsymbol{U\Sigma V^T})^T\boldsymbol{y}=\sum_{j=0}^{p-1}\boldsymbol{u}_j\boldsymbol{u}_j^T\frac{\sigma_j^2}{\sigma_j^2+\lambda}\boldsymbol{y},
\]</div>
<p>with the vectors <span class="math notranslate nohighlight">\(\boldsymbol{u}_j\)</span> being the columns of <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> from the SVD of the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>.</p>
</section>
<section id="interpreting-the-ridge-results">
<h2>Interpreting the Ridge results<a class="headerlink" href="#interpreting-the-ridge-results" title="Link to this heading">#</a></h2>
<p>Since <span class="math notranslate nohighlight">\(\lambda \geq 0\)</span>, it means that compared to OLS, we have</p>
<div class="math notranslate nohighlight">
\[
\frac{\sigma_j^2}{\sigma_j^2+\lambda} \leq 1.
\]</div>
<p>Ridge regression finds the coordinates of <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> with respect to the
orthonormal basis <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span>, it then shrinks the coordinates by
<span class="math notranslate nohighlight">\(\frac{\sigma_j^2}{\sigma_j^2+\lambda}\)</span>. Recall that the SVD has
eigenvalues ordered in a descending way, that is <span class="math notranslate nohighlight">\(\sigma_i \geq
\sigma_{i+1}\)</span>.</p>
<p>For small eigenvalues <span class="math notranslate nohighlight">\(\sigma_i\)</span> it means that their contributions become less important, a fact which can be used to reduce the number of degrees of freedom. More about this when we have covered the material on a statistical interpretation of various linear regression methods.</p>
</section>
<section id="more-interpretations">
<h2>More interpretations<a class="headerlink" href="#more-interpretations" title="Link to this heading">#</a></h2>
<p>For the sake of simplicity, let us assume that the design matrix is orthonormal, that is</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{X}^T\boldsymbol{X}=(\boldsymbol{X}^T\boldsymbol{X})^{-1} =\boldsymbol{I}.
\]</div>
<p>In this case the standard OLS results in</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{\theta}^{\mathrm{OLS}} = \boldsymbol{X}^T\boldsymbol{y},
\]</div>
<p>and</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{\theta}^{\mathrm{Ridge}} = \left(\boldsymbol{I}+\lambda\boldsymbol{I}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}=\left(1+\lambda\right)^{-1}\boldsymbol{\theta}^{\mathrm{OLS}},
\]</div>
<p>that is the Ridge estimator scales the OLS estimator by the inverse of a factor <span class="math notranslate nohighlight">\(1+\lambda\)</span>, and
the Ridge estimator converges to zero when the hyperparameter goes to
infinity.</p>
<p>We will come back to more interpretions after we have gone through some of the statistical analysis part.</p>
<p>For more discussions of Ridge and Lasso regression, <a class="reference external" href="https://arxiv.org/abs/1509.09169">Wessel van Wieringen’s</a> article is highly recommended.
Similarly, <a class="reference external" href="https://arxiv.org/abs/1803.08823">Mehta et al’s article</a> is also recommended.</p>
</section>
<section id="deriving-the-lasso-regression-equations">
<h2>Deriving the Lasso Regression Equations<a class="headerlink" href="#deriving-the-lasso-regression-equations" title="Link to this heading">#</a></h2>
<p>Using the matrix-vector expression for Lasso regression, we have the following <strong>cost</strong> function</p>
<div class="math notranslate nohighlight">
\[
C(\boldsymbol{X},\boldsymbol{\theta})=\frac{1}{n}\left\{(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta})^T(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta})\right\}+\lambda\vert\vert\boldsymbol{\theta}\vert\vert_1,
\]</div>
<p>Taking the derivative with respect to <span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span> and recalling that the derivative of the absolute value is (we drop the boldfaced vector symbol for simplicty)</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\frac{d \vert \theta\vert}{d \theta}=\mathrm{sgn}(\theta)=\left\{\begin{array}{cc} 1 &amp; \theta &gt; 0 \\-1 &amp; \theta &lt; 0, \end{array}\right.
\end{split}\]</div>
<p>we have that the derivative of the cost function is</p>
<div class="math notranslate nohighlight">
\[
\frac{\partial C(\boldsymbol{X},\boldsymbol{\theta})}{\partial \boldsymbol{\theta}}=-\frac{2}{n}\boldsymbol{X}^T(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta})+\lambda sgn(\boldsymbol{\theta})=0,
\]</div>
<p>and reordering we have</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\theta}+\frac{n}{2}\lambda sgn(\boldsymbol{\theta})=\boldsymbol{X}^T\boldsymbol{y}.
\]</div>
<p>We can redefine <span class="math notranslate nohighlight">\(\lambda\)</span> to absorb the constant <span class="math notranslate nohighlight">\(n/2\)</span> and we rewrite the last equation as</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\theta}+\lambda sgn(\boldsymbol{\theta})=\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.We will discuss how to code the above methods using gradient descent methods.</p>
</section>
<section id="optimization-and-gradient-descent-the-central-part-of-any-machine-learning-algortithm">
<h2>Optimization and gradient descent, the central part of any Machine Learning algortithm<a class="headerlink" href="#optimization-and-gradient-descent-the-central-part-of-any-machine-learning-algortithm" title="Link to this heading">#</a></h2>
<p>Almost every problem in machine learning and data science starts with
a dataset <span class="math notranslate nohighlight">\(X\)</span>, a model <span class="math notranslate nohighlight">\(g(\theta)\)</span>, which is a function of the
parameters <span class="math notranslate nohighlight">\(\theta\)</span> and a cost function <span class="math notranslate nohighlight">\(C(X, g(\theta))\)</span> that allows
us to judge how well the model <span class="math notranslate nohighlight">\(g(\theta)\)</span> explains the observations
<span class="math notranslate nohighlight">\(X\)</span>. The model is fit by finding the values of <span class="math notranslate nohighlight">\(\theta\)</span> that minimize
the cost function. Ideally we would be able to solve for <span class="math notranslate nohighlight">\(\theta\)</span>
analytically, however this is not possible in general and we must use
some approximative/numerical method to compute the minimum.</p>
</section>
<section id="reminder-on-newton-raphson-s-method">
<h2>Reminder on Newton-Raphson’s method<a class="headerlink" href="#reminder-on-newton-raphson-s-method" title="Link to this heading">#</a></h2>
<p>Let us quickly remind ourselves how we derive the above method.</p>
<p>Perhaps the most celebrated of all one-dimensional root-finding
routines is Newton’s method, also called the Newton-Raphson
method. This method requires the evaluation of both the
function <span class="math notranslate nohighlight">\(f\)</span> and its derivative <span class="math notranslate nohighlight">\(f'\)</span> at arbitrary points.
If you can only calculate the derivative
numerically and/or your function is not of the smooth type, we
normally discourage the use of this method.</p>
</section>
<section id="the-equations">
<h2>The equations<a class="headerlink" href="#the-equations" title="Link to this heading">#</a></h2>
<p>The Newton-Raphson formula consists geometrically of extending the
tangent line at a current point until it crosses zero, then setting
the next guess to the abscissa of that zero-crossing. The mathematics
behind this method is rather simple. Employing a Taylor expansion for
<span class="math notranslate nohighlight">\(x\)</span> sufficiently close to the solution <span class="math notranslate nohighlight">\(s\)</span>, we have</p>
<!-- Equation labels as ordinary links -->
<div id="eq:taylornr"></div>
<div class="math notranslate nohighlight">
\[
f(s)=0=f(x)+(s-x)f'(x)+\frac{(s-x)^2}{2}f''(x) +\dots.
\label{eq:taylornr} \tag{2}
\]</div>
<p>For small enough values of the function and for well-behaved
functions, the terms beyond linear are unimportant, hence we obtain</p>
<div class="math notranslate nohighlight">
\[
f(x)+(s-x)f'(x)\approx 0,
\]</div>
<p>yielding</p>
<div class="math notranslate nohighlight">
\[
s\approx x-\frac{f(x)}{f'(x)}.
\]</div>
<p>Having in mind an iterative procedure, it is natural to start iterating with</p>
<div class="math notranslate nohighlight">
\[
x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}.
\]</div>
</section>
<section id="simple-geometric-interpretation">
<h2>Simple geometric interpretation<a class="headerlink" href="#simple-geometric-interpretation" title="Link to this heading">#</a></h2>
<p>The above is Newton-Raphson’s method. It has a simple geometric
interpretation, namely <span class="math notranslate nohighlight">\(x_{n+1}\)</span> is the point where the tangent from
<span class="math notranslate nohighlight">\((x_n,f(x_n))\)</span> crosses the <span class="math notranslate nohighlight">\(x\)</span>-axis. Close to the solution,
Newton-Raphson converges fast to the desired result. However, if we
are far from a root, where the higher-order terms in the series are
important, the Newton-Raphson formula can give grossly inaccurate
results. For instance, the initial guess for the root might be so far
from the true root as to let the search interval include a local
maximum or minimum of the function. If an iteration places a trial
guess near such a local extremum, so that the first derivative nearly
vanishes, then Newton-Raphson may fail totally</p>
</section>
<section id="extending-to-more-than-one-variable">
<h2>Extending to more than one variable<a class="headerlink" href="#extending-to-more-than-one-variable" title="Link to this heading">#</a></h2>
<p>Newton’s method can be generalized to systems of several non-linear equations
and variables. Consider the case with two equations</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{array}{cc} f_1(x_1,x_2) &amp;=0\\
f_2(x_1,x_2) &amp;=0,\end{array}
\end{split}\]</div>
<p>which we Taylor expand to obtain</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{array}{cc} 0=f_1(x_1+h_1,x_2+h_2)=&amp;f_1(x_1,x_2)+h_1
\partial f_1/\partial x_1+h_2
\partial f_1/\partial x_2+\dots\\
0=f_2(x_1+h_1,x_2+h_2)=&amp;f_2(x_1,x_2)+h_1
\partial f_2/\partial x_1+h_2
\partial f_2/\partial x_2+\dots
\end{array}.
\end{split}\]</div>
<p>Defining the Jacobian matrix <span class="math notranslate nohighlight">\({\bf \boldsymbol{J}}\)</span> we have</p>
<div class="math notranslate nohighlight">
\[\begin{split}
{\bf \boldsymbol{J}}=\left( \begin{array}{cc}
\partial f_1/\partial x_1 &amp; \partial f_1/\partial x_2 \\
\partial f_2/\partial x_1 &amp;\partial f_2/\partial x_2
\end{array} \right),
\end{split}\]</div>
<p>we can rephrase Newton’s method as</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\left(\begin{array}{c} x_1^{n+1} \\ x_2^{n+1} \end{array} \right)=
\left(\begin{array}{c} x_1^{n} \\ x_2^{n} \end{array} \right)+
\left(\begin{array}{c} h_1^{n} \\ h_2^{n} \end{array} \right),
\end{split}\]</div>
<p>where we have defined</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\left(\begin{array}{c} h_1^{n} \\ h_2^{n} \end{array} \right)=
-{\bf \boldsymbol{J}}^{-1}
\left(\begin{array}{c} f_1(x_1^{n},x_2^{n}) \\ f_2(x_1^{n},x_2^{n}) \end{array} \right).
\end{split}\]</div>
<p>We need thus to compute the inverse of the Jacobian matrix and it
is to understand that difficulties may
arise in case <span class="math notranslate nohighlight">\({\bf \boldsymbol{J}}\)</span> is nearly singular.</p>
<p>It is rather straightforward to extend the above scheme to systems of
more than two non-linear equations. In our case, the Jacobian matrix is given by the Hessian that represents the second derivative of cost function.</p>
</section>
<section id="steepest-descent">
<h2>Steepest descent<a class="headerlink" href="#steepest-descent" title="Link to this heading">#</a></h2>
<p>The basic idea of gradient descent is
that a function <span class="math notranslate nohighlight">\(F(\mathbf{x})\)</span>,
<span class="math notranslate nohighlight">\(\mathbf{x} \equiv (x_1,\cdots,x_n)\)</span>, decreases fastest if one goes from <span class="math notranslate nohighlight">\(\bf {x}\)</span> in the
direction of the negative gradient <span class="math notranslate nohighlight">\(-\nabla F(\mathbf{x})\)</span>.</p>
<p>It can be shown that if</p>
<div class="math notranslate nohighlight">
\[
\mathbf{x}_{k+1} = \mathbf{x}_k - \gamma_k \nabla F(\mathbf{x}_k),
\]</div>
<p>with <span class="math notranslate nohighlight">\(\gamma_k &gt; 0\)</span>.</p>
<p>For <span class="math notranslate nohighlight">\(\gamma_k\)</span> small enough, then <span class="math notranslate nohighlight">\(F(\mathbf{x}_{k+1}) \leq
F(\mathbf{x}_k)\)</span>. This means that for a sufficiently small <span class="math notranslate nohighlight">\(\gamma_k\)</span>
we are always moving towards smaller function values, i.e a minimum.</p>
</section>
<section id="more-on-steepest-descent">
<h2>More on Steepest descent<a class="headerlink" href="#more-on-steepest-descent" title="Link to this heading">#</a></h2>
<p>The previous observation is the basis of the method of steepest
descent, which is also referred to as just gradient descent (GD). One
starts with an initial guess <span class="math notranslate nohighlight">\(\mathbf{x}_0\)</span> for a minimum of <span class="math notranslate nohighlight">\(F\)</span> and
computes new approximations according to</p>
<div class="math notranslate nohighlight">
\[
\mathbf{x}_{k+1} = \mathbf{x}_k - \gamma_k \nabla F(\mathbf{x}_k), \ \ k \geq 0.
\]</div>
<p>The parameter <span class="math notranslate nohighlight">\(\gamma_k\)</span> is often referred to as the step length or
the learning rate within the context of Machine Learning.</p>
</section>
<section id="the-ideal">
<h2>The ideal<a class="headerlink" href="#the-ideal" title="Link to this heading">#</a></h2>
<p>Ideally the sequence <span class="math notranslate nohighlight">\(\{\mathbf{x}_k \}_{k=0}\)</span> converges to a global
minimum of the function <span class="math notranslate nohighlight">\(F\)</span>. In general we do not know if we are in a
global or local minimum. In the special case when <span class="math notranslate nohighlight">\(F\)</span> is a convex
function, all local minima are also global minima, so in this case
gradient descent can converge to the global solution. The advantage of
this scheme is that it is conceptually simple and straightforward to
implement. However the method in this form has some severe
limitations:</p>
<p>In machine learing we are often faced with non-convex high dimensional
cost functions with many local minima. Since GD is deterministic we
will get stuck in a local minimum, if the method converges, unless we
have a very good intial guess. This also implies that the scheme is
sensitive to the chosen initial condition.</p>
<p>Note that the gradient is a function of <span class="math notranslate nohighlight">\(\mathbf{x} =
(x_1,\cdots,x_n)\)</span> which makes it expensive to compute numerically.</p>
</section>
<section id="the-sensitiveness-of-the-gradient-descent">
<h2>The sensitiveness of the gradient descent<a class="headerlink" href="#the-sensitiveness-of-the-gradient-descent" title="Link to this heading">#</a></h2>
<p>The gradient descent method
is sensitive to the choice of learning rate <span class="math notranslate nohighlight">\(\gamma_k\)</span>. This is due
to the fact that we are only guaranteed that <span class="math notranslate nohighlight">\(F(\mathbf{x}_{k+1}) \leq
F(\mathbf{x}_k)\)</span> for sufficiently small <span class="math notranslate nohighlight">\(\gamma_k\)</span>. The problem is to
determine an optimal learning rate. If the learning rate is chosen too
small the method will take a long time to converge and if it is too
large we can experience erratic behavior.</p>
<p>Many of these shortcomings can be alleviated by introducing
randomness. One such method is that of Stochastic Gradient Descent
(SGD), to be discussed next week.</p>
</section>
<section id="convex-functions">
<h2>Convex functions<a class="headerlink" href="#convex-functions" title="Link to this heading">#</a></h2>
<p>Ideally we want our cost/loss function to be convex(concave).</p>
<p>First we give the definition of a convex set: A set <span class="math notranslate nohighlight">\(C\)</span> in
<span class="math notranslate nohighlight">\(\mathbb{R}^n\)</span> is said to be convex if, for all <span class="math notranslate nohighlight">\(x\)</span> and <span class="math notranslate nohighlight">\(y\)</span> in <span class="math notranslate nohighlight">\(C\)</span> and
all <span class="math notranslate nohighlight">\(t \in (0,1)\)</span> , the point <span class="math notranslate nohighlight">\((1 − t)x + ty\)</span> also belongs to
C. Geometrically this means that every point on the line segment
connecting <span class="math notranslate nohighlight">\(x\)</span> and <span class="math notranslate nohighlight">\(y\)</span> is in <span class="math notranslate nohighlight">\(C\)</span> as discussed below.</p>
<p>The convex subsets of <span class="math notranslate nohighlight">\(\mathbb{R}\)</span> are the intervals of
<span class="math notranslate nohighlight">\(\mathbb{R}\)</span>. Examples of convex sets of <span class="math notranslate nohighlight">\(\mathbb{R}^2\)</span> are the
regular polygons (triangles, rectangles, pentagons, etc…).</p>
</section>
<section id="convex-function">
<h2>Convex function<a class="headerlink" href="#convex-function" title="Link to this heading">#</a></h2>
<p><strong>Convex function</strong>: Let <span class="math notranslate nohighlight">\(X \subset \mathbb{R}^n\)</span> be a convex set. Assume that the function <span class="math notranslate nohighlight">\(f: X \rightarrow \mathbb{R}\)</span> is continuous, then <span class="math notranslate nohighlight">\(f\)</span> is said to be convex if $<span class="math notranslate nohighlight">\(f(tx_1 + (1-t)x_2) \leq tf(x_1) + (1-t)f(x_2) \)</span><span class="math notranslate nohighlight">\( for all \)</span>x_1, x_2 \in X<span class="math notranslate nohighlight">\( and for all \)</span>t \in [0,1]<span class="math notranslate nohighlight">\(. If \)</span>\leq<span class="math notranslate nohighlight">\( is replaced with a strict inequaltiy in the definition, we demand \)</span>x_1 \neq x_2<span class="math notranslate nohighlight">\( and \)</span>t\in(0,1)<span class="math notranslate nohighlight">\( then \)</span>f<span class="math notranslate nohighlight">\( is said to be strictly convex. For a single variable function, convexity means that if you draw a straight line connecting \)</span>f(x_1)<span class="math notranslate nohighlight">\( and \)</span>f(x_2)<span class="math notranslate nohighlight">\(, the value of the function on the interval \)</span>[x_1,x_2]$ is always below the line as illustrated below.</p>
</section>
<section id="conditions-on-convex-functions">
<h2>Conditions on convex functions<a class="headerlink" href="#conditions-on-convex-functions" title="Link to this heading">#</a></h2>
<p>In the following we state first and second-order conditions which
ensures convexity of a function <span class="math notranslate nohighlight">\(f\)</span>. We write <span class="math notranslate nohighlight">\(D_f\)</span> to denote the
domain of <span class="math notranslate nohighlight">\(f\)</span>, i.e the subset of <span class="math notranslate nohighlight">\(R^n\)</span> where <span class="math notranslate nohighlight">\(f\)</span> is defined. For more
details and proofs we refer to: [S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press](<a class="reference external" href="http://stanford.edu/boyd/cvxbook/">http://stanford.edu/boyd/cvxbook/</a>, 2004).</p>
<p><strong>First order condition.</strong></p>
<p>Suppose <span class="math notranslate nohighlight">\(f\)</span> is differentiable (i.e <span class="math notranslate nohighlight">\(\nabla f(x)\)</span> is well defined for
all <span class="math notranslate nohighlight">\(x\)</span> in the domain of <span class="math notranslate nohighlight">\(f\)</span>). Then <span class="math notranslate nohighlight">\(f\)</span> is convex if and only if <span class="math notranslate nohighlight">\(D_f\)</span>
is a convex set and $<span class="math notranslate nohighlight">\(f(y) \geq f(x) + \nabla f(x)^T (y-x) \)</span><span class="math notranslate nohighlight">\( holds
for all \)</span>x,y \in D_f<span class="math notranslate nohighlight">\(. This condition means that for a convex function
the first order Taylor expansion (right hand side above) at any point
a global under estimator of the function. To convince yourself you can
make a drawing of \)</span>f(x) = x^2+1<span class="math notranslate nohighlight">\( and draw the tangent line to \)</span>f(x)$ and
note that it is always below the graph.</p>
<p><strong>Second order condition.</strong></p>
<p>Assume that <span class="math notranslate nohighlight">\(f\)</span> is twice
differentiable, i.e the Hessian matrix exists at each point in
<span class="math notranslate nohighlight">\(D_f\)</span>. Then <span class="math notranslate nohighlight">\(f\)</span> is convex if and only if <span class="math notranslate nohighlight">\(D_f\)</span> is a convex set and its
Hessian is positive semi-definite for all <span class="math notranslate nohighlight">\(x\in D_f\)</span>. For a
single-variable function this reduces to <span class="math notranslate nohighlight">\(f''(x) \geq 0\)</span>. Geometrically this means that <span class="math notranslate nohighlight">\(f\)</span> has nonnegative curvature
everywhere.</p>
<p>This condition is particularly useful since it gives us an procedure for determining if the function under consideration is convex, apart from using the definition.</p>
</section>
<section id="more-on-convex-functions">
<h2>More on convex functions<a class="headerlink" href="#more-on-convex-functions" title="Link to this heading">#</a></h2>
<p>The next result is of great importance to us and the reason why we are
going on about convex functions. In machine learning we frequently
have to minimize a loss/cost function in order to find the best
parameters for the model we are considering.</p>
<p>Ideally we want the
global minimum (for high-dimensional models it is hard to know
if we have local or global minimum). However, if the cost/loss function
is convex the following result provides invaluable information:</p>
<p><strong>Any minimum is global for convex functions.</strong></p>
<p>Consider the problem of finding <span class="math notranslate nohighlight">\(x \in \mathbb{R}^n\)</span> such that <span class="math notranslate nohighlight">\(f(x)\)</span>
is minimal, where <span class="math notranslate nohighlight">\(f\)</span> is convex and differentiable. Then, any point
<span class="math notranslate nohighlight">\(x^*\)</span> that satisfies <span class="math notranslate nohighlight">\(\nabla f(x^*) = 0\)</span> is a global minimum.</p>
<p>This result means that if we know that the cost/loss function is convex and we are able to find a minimum, we are guaranteed that it is a global minimum.</p>
</section>
<section id="some-simple-problems">
<h2>Some simple problems<a class="headerlink" href="#some-simple-problems" title="Link to this heading">#</a></h2>
<ol class="arabic simple">
<li><p>Show that <span class="math notranslate nohighlight">\(f(x)=x^2\)</span> is convex for <span class="math notranslate nohighlight">\(x \in \mathbb{R}\)</span> using the definition of convexity. Hint: If you re-write the definition, <span class="math notranslate nohighlight">\(f\)</span> is convex if the following holds for all <span class="math notranslate nohighlight">\(x,y \in D_f\)</span> and any <span class="math notranslate nohighlight">\(\lambda \in [0,1]\)</span> <span class="math notranslate nohighlight">\(\lambda f(x)+(1-\lambda)f(y)-f(\lambda x + (1-\lambda) y ) \geq 0\)</span>.</p></li>
<li><p>Using the second order condition show that the following functions are convex on the specified domain.</p></li>
</ol>
<ul class="simple">
<li><p><span class="math notranslate nohighlight">\(f(x) = e^x\)</span> is convex for <span class="math notranslate nohighlight">\(x \in \mathbb{R}\)</span>.</p></li>
<li><p><span class="math notranslate nohighlight">\(g(x) = -\ln(x)\)</span> is convex for <span class="math notranslate nohighlight">\(x \in (0,\infty)\)</span>.</p></li>
</ul>
<ol class="arabic simple" start="3">
<li><p>Let <span class="math notranslate nohighlight">\(f(x) = x^2\)</span> and <span class="math notranslate nohighlight">\(g(x) = e^x\)</span>. Show that <span class="math notranslate nohighlight">\(f(g(x))\)</span> and <span class="math notranslate nohighlight">\(g(f(x))\)</span> is convex for <span class="math notranslate nohighlight">\(x \in \mathbb{R}\)</span>. Also show that if <span class="math notranslate nohighlight">\(f(x)\)</span> is any convex function than <span class="math notranslate nohighlight">\(h(x) = e^{f(x)}\)</span> is convex.</p></li>
<li><p>A norm is any function that satisfy the following properties</p></li>
</ol>
<ul class="simple">
<li><p><span class="math notranslate nohighlight">\(f(\alpha x) = |\alpha| f(x)\)</span> for all <span class="math notranslate nohighlight">\(\alpha \in \mathbb{R}\)</span>.</p></li>
<li><p><span class="math notranslate nohighlight">\(f(x+y) \leq f(x) + f(y)\)</span></p></li>
<li><p><span class="math notranslate nohighlight">\(f(x) \leq 0\)</span> for all <span class="math notranslate nohighlight">\(x \in \mathbb{R}^n\)</span> with equality if and only if <span class="math notranslate nohighlight">\(x = 0\)</span></p></li>
</ul>
<p>Using the definition of convexity, try to show that a function satisfying the properties above is convex (the third condition is not needed to show this).</p>
</section>
<section id="revisiting-ordinary-least-squares">
<h2>Revisiting Ordinary Least Squares<a class="headerlink" href="#revisiting-ordinary-least-squares" title="Link to this heading">#</a></h2>
<p>We will use linear regression as a case study for the gradient descent
methods. Linear regression is a great test case for the gradient
descent methods discussed in the lectures since it has several
desirable properties such as:</p>
<ol class="arabic simple">
<li><p>An analytical solution (recall homework sets for week 35).</p></li>
<li><p>The gradient can be computed analytically.</p></li>
<li><p>The cost function is convex which guarantees that gradient descent converges for small enough learning rates</p></li>
</ol>
<p>We revisit an example similar to what we had in the first homework set. We had a function of the type</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span>x = 2*np.random.rand(m,1)
y = 4+3*x+np.random.randn(m,1)
</pre></div>
</div>
</div>
</div>
<p>with <span class="math notranslate nohighlight">\(x_i \in [0,1] \)</span> is chosen randomly using a uniform distribution. Additionally we have a stochastic noise chosen according to a normal distribution <span class="math notranslate nohighlight">\(\cal {N}(0,1)\)</span>.
The linear regression model is given by</p>
<div class="math notranslate nohighlight">
\[
h_\theta(x) = \boldsymbol{y} = \theta_0 + \theta_1 x,
\]</div>
<p>such that</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{y}_i = \theta_0 + \theta_1 x_i.
\]</div>
</section>
<section id="gradient-descent-example">
<h2>Gradient descent example<a class="headerlink" href="#gradient-descent-example" title="Link to this heading">#</a></h2>
<p>Let <span class="math notranslate nohighlight">\(\mathbf{y} = (y_1,\cdots,y_n)^T\)</span>, <span class="math notranslate nohighlight">\(\mathbf{\boldsymbol{y}} = (\boldsymbol{y}_1,\cdots,\boldsymbol{y}_n)^T\)</span> and <span class="math notranslate nohighlight">\(\theta = (\theta_0, \theta_1)^T\)</span></p>
<p>It is convenient to write <span class="math notranslate nohighlight">\(\mathbf{\boldsymbol{y}} = X\theta\)</span> where <span class="math notranslate nohighlight">\(X \in \mathbb{R}^{100 \times 2} \)</span> is the design matrix given by (we keep the intercept here)</p>
<div class="math notranslate nohighlight">
\[\begin{split}
X \equiv \begin{bmatrix}
1 &amp; x_1 \\
\vdots &amp; \vdots \\
1 &amp; x_{100} &amp; \\
\end{bmatrix}.
\end{split}\]</div>
<p>The cost/loss/risk function is given by (</p>
<div class="math notranslate nohighlight">
\[
C(\theta) = \frac{1}{n}||X\theta-\mathbf{y}||_{2}^{2} = \frac{1}{n}\sum_{i=1}^{100}\left[ (\theta_0 + \theta_1 x_i)^2 - 2 y_i (\theta_0 + \theta_1 x_i) + y_i^2\right]
\]</div>
<p>and we want to find <span class="math notranslate nohighlight">\(\theta\)</span> such that <span class="math notranslate nohighlight">\(C(\theta)\)</span> is minimized.</p>
</section>
<section id="the-derivative-of-the-cost-loss-function">
<h2>The derivative of the cost/loss function<a class="headerlink" href="#the-derivative-of-the-cost-loss-function" title="Link to this heading">#</a></h2>
<p>Computing <span class="math notranslate nohighlight">\(\partial C(\theta) / \partial \theta_0\)</span> and <span class="math notranslate nohighlight">\(\partial C(\theta) / \partial \theta_1\)</span> we can show that the gradient can be written as</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\nabla_{\theta} C(\theta) = \frac{2}{n}\begin{bmatrix} \sum_{i=1}^{100} \left(\theta_0+\theta_1x_i-y_i\right) \\
\sum_{i=1}^{100}\left( x_i (\theta_0+\theta_1x_i)-y_ix_i\right) \\
\end{bmatrix} = \frac{2}{n}X^T(X\theta - \mathbf{y}),
\end{split}\]</div>
<p>where <span class="math notranslate nohighlight">\(X\)</span> is the design matrix defined above.</p>
</section>
<section id="the-hessian-matrix">
<h2>The Hessian matrix<a class="headerlink" href="#the-hessian-matrix" title="Link to this heading">#</a></h2>
<p>The Hessian matrix of <span class="math notranslate nohighlight">\(C(\theta)\)</span> is given by</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\boldsymbol{H} \equiv \begin{bmatrix}
\frac{\partial^2 C(\theta)}{\partial \theta_0^2} &amp; \frac{\partial^2 C(\theta)}{\partial \theta_0 \partial \theta_1} \\
\frac{\partial^2 C(\theta)}{\partial \theta_0 \partial \theta_1} &amp; \frac{\partial^2 C(\theta)}{\partial \theta_1^2} &amp; \\
\end{bmatrix} = \frac{2}{n}X^T X.
\end{split}\]</div>
<p>This result implies that <span class="math notranslate nohighlight">\(C(\theta)\)</span> is a convex function since the matrix <span class="math notranslate nohighlight">\(X^T X\)</span> always is positive semi-definite.</p>
</section>
<section id="simple-program">
<h2>Simple program<a class="headerlink" href="#simple-program" title="Link to this heading">#</a></h2>
<p>We can now write a program that minimizes <span class="math notranslate nohighlight">\(C(\theta)\)</span> using the gradient descent method with a constant learning rate <span class="math notranslate nohighlight">\(\gamma\)</span> according to</p>
<div class="math notranslate nohighlight">
\[
\theta_{k+1} = \theta_k - \gamma \nabla_\theta C(\theta_k), \ k=0,1,\cdots
\]</div>
<p>We can use the expression we computed for the gradient and let use a
<span class="math notranslate nohighlight">\(\theta_0\)</span> be chosen randomly and let <span class="math notranslate nohighlight">\(\gamma = 0.001\)</span>. Stop iterating
when <span class="math notranslate nohighlight">\(||\nabla_\theta C(\theta_k) || \leq \epsilon = 10^{-8}\)</span>. <strong>Note that the code below does not include the latter stop criterion</strong>.</p>
<p>And finally we can compare our solution for <span class="math notranslate nohighlight">\(\theta\)</span> with the analytic result given by
<span class="math notranslate nohighlight">\(\theta= (X^TX)^{-1} X^T \mathbf{y}\)</span>.</p>
</section>
<section id="id1">
<h2>Gradient Descent Example<a class="headerlink" href="#id1" title="Link to this heading">#</a></h2>
<p>Here our simple example</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span>
# Importing various packages
from random import random, seed
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import sys
# the number of datapoints
n = 100
x = 2*np.random.rand(n,1)
y = 4+3*x+np.random.randn(n,1)
X = np.c_[np.ones((n,1)), x]
# Hessian matrix
H = (2.0/n)* X.T @ X
# Get the eigenvalues
EigValues, EigVectors = np.linalg.eig(H)
print(f&quot;Eigenvalues of Hessian Matrix:{EigValues}&quot;)
theta_linreg = np.linalg.inv(X.T @ X) @ X.T @ y
print(theta_linreg)
theta = np.random.randn(2,1)
eta = 1.0/np.max(EigValues)
Niterations = 1000
for iter in range(Niterations):
gradient = (2.0/n)*X.T @ (X @ theta-y)
theta -= eta*gradient
print(theta)
xnew = np.array([[0],[2]])
xbnew = np.c_[np.ones((2,1)), xnew]
ypredict = xbnew.dot(theta)
ypredict2 = xbnew.dot(theta_linreg)
plt.plot(xnew, ypredict, &quot;r-&quot;)
plt.plot(xnew, ypredict2, &quot;b-&quot;)
plt.plot(x, y ,&#39;ro&#39;)
plt.axis([0,2.0,0, 15.0])
plt.xlabel(r&#39;$x$&#39;)
plt.ylabel(r&#39;$y$&#39;)
plt.title(r&#39;Gradient descent example&#39;)
plt.show()
</pre></div>
</div>
</div>
</div>
</section>
<section id="gradient-descent-and-ridge">
<h2>Gradient descent and Ridge<a class="headerlink" href="#gradient-descent-and-ridge" title="Link to this heading">#</a></h2>
<p>We have also discussed Ridge regression where the loss function contains a regularized term given by the <span class="math notranslate nohighlight">\(L_2\)</span> norm of <span class="math notranslate nohighlight">\(\theta\)</span>,</p>
<div class="math notranslate nohighlight">
\[
C_{\text{ridge}}(\theta) = \frac{1}{n}||X\theta -\mathbf{y}||^2 + \lambda ||\theta||^2, \ \lambda \geq 0.
\]</div>
<p>In order to minimize <span class="math notranslate nohighlight">\(C_{\text{ridge}}(\theta)\)</span> using GD we adjust the gradient as follows</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\nabla_\theta C_{\text{ridge}}(\theta) = \frac{2}{n}\begin{bmatrix} \sum_{i=1}^{100} \left(\theta_0+\theta_1x_i-y_i\right) \\
\sum_{i=1}^{100}\left( x_i (\theta_0+\theta_1x_i)-y_ix_i\right) \\
\end{bmatrix} + 2\lambda\begin{bmatrix} \theta_0 \\ \theta_1\end{bmatrix} = 2 (\frac{1}{n}X^T(X\theta - \mathbf{y})+\lambda \theta).
\end{split}\]</div>
<p>We can easily extend our program to minimize <span class="math notranslate nohighlight">\(C_{\text{ridge}}(\theta)\)</span> using gradient descent and compare with the analytical solution given by</p>
<div class="math notranslate nohighlight">
\[
\theta_{\text{ridge}} = \left(X^T X + n\lambda I_{2 \times 2} \right)^{-1} X^T \mathbf{y}.
\]</div>
</section>
<section id="the-hessian-matrix-for-ridge-regression">
<h2>The Hessian matrix for Ridge Regression<a class="headerlink" href="#the-hessian-matrix-for-ridge-regression" title="Link to this heading">#</a></h2>
<p>The Hessian matrix of Ridge Regression for our simple example is given by</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\boldsymbol{H} \equiv \begin{bmatrix}
\frac{\partial^2 C(\theta)}{\partial \theta_0^2} &amp; \frac{\partial^2 C(\theta)}{\partial \theta_0 \partial \theta_1} \\
\frac{\partial^2 C(\theta)}{\partial \theta_0 \partial \theta_1} &amp; \frac{\partial^2 C(\theta)}{\partial \theta_1^2} &amp; \\
\end{bmatrix} = \frac{2}{n}X^T X+2\lambda\boldsymbol{I}.
\end{split}\]</div>
<p>This implies that the Hessian matrix is positive definite, hence the stationary point is a
minimum.
Note that the Ridge cost function is convex being a sum of two convex
functions. Therefore, the stationary point is a global
minimum of this function.</p>
</section>
<section id="program-example-for-gradient-descent-with-ridge-regression">
<h2>Program example for gradient descent with Ridge Regression<a class="headerlink" href="#program-example-for-gradient-descent-with-ridge-regression" 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>from random import random, seed
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import sys
# the number of datapoints
n = 100
x = 2*np.random.rand(n,1)
y = 4+3*x+np.random.randn(n,1)
X = np.c_[np.ones((n,1)), x]
XT_X = X.T @ X
#Ridge parameter lambda
lmbda = 0.001
Id = n*lmbda* np.eye(XT_X.shape[0])
# Hessian matrix
H = (2.0/n)* XT_X+2*lmbda* np.eye(XT_X.shape[0])
# Get the eigenvalues
EigValues, EigVectors = np.linalg.eig(H)
print(f&quot;Eigenvalues of Hessian Matrix:{EigValues}&quot;)
theta_linreg = np.linalg.inv(XT_X+Id) @ X.T @ y
print(theta_linreg)
# Start plain gradient descent
theta = np.random.randn(2,1)
eta = 1.0/np.max(EigValues)
Niterations = 100
for iter in range(Niterations):
gradients = 2.0/n*X.T @ (X @ (theta)-y)+2*lmbda*theta
theta -= eta*gradients
print(theta)
ypredict = X @ theta
ypredict2 = X @ theta_linreg
plt.plot(x, ypredict, &quot;r-&quot;)
plt.plot(x, ypredict2, &quot;b-&quot;)
plt.plot(x, y ,&#39;ro&#39;)
plt.axis([0,2.0,0, 15.0])
plt.xlabel(r&#39;$x$&#39;)
plt.ylabel(r&#39;$y$&#39;)
plt.title(r&#39;Gradient descent example for Ridge&#39;)
plt.show()
</pre></div>
</div>
</div>
</div>
</section>
<section id="using-gradient-descent-methods-limitations">
<h2>Using gradient descent methods, limitations<a class="headerlink" href="#using-gradient-descent-methods-limitations" title="Link to this heading">#</a></h2>
<ul class="simple">
<li><p><strong>Gradient descent (GD) finds local minima of our function</strong>. Since the GD algorithm is deterministic, if it converges, it will converge to a local minimum of our cost/loss/risk function. Because in ML we are often dealing with extremely rugged landscapes with many local minima, this can lead to poor performance.</p></li>
<li><p><strong>GD is sensitive to initial conditions</strong>. One consequence of the local nature of GD is that initial conditions matter. Depending on where one starts, one will end up at a different local minima. Therefore, it is very important to think about how one initializes the training process. This is true for GD as well as more complicated variants of GD.</p></li>
<li><p><strong>Gradients are computationally expensive to calculate for large datasets</strong>. In many cases in statistics and ML, the cost/loss/risk function is a sum of terms, with one term for each data point. For example, in linear regression, <span class="math notranslate nohighlight">\(E \propto \sum_{i=1}^n (y_i - \mathbf{w}^T\cdot\mathbf{x}_i)^2\)</span>; for logistic regression, the square error is replaced by the cross entropy. To calculate the gradient we have to sum over <em>all</em> <span class="math notranslate nohighlight">\(n\)</span> data points. Doing this at every GD step becomes extremely computationally expensive. An ingenious solution to this, is to calculate the gradients using small subsets of the data called “mini batches”. This has the added benefit of introducing stochasticity into our algorithm.</p></li>
<li><p><strong>GD is very sensitive to choices of learning rates</strong>. GD is extremely sensitive to the choice of learning rates. If the learning rate is very small, the training process take an extremely long time. For larger learning rates, GD can diverge and give poor results. Furthermore, depending on what the local landscape looks like, we have to modify the learning rates to ensure convergence. Ideally, we would <em>adaptively</em> choose the learning rates to match the landscape.</p></li>
<li><p><strong>GD treats all directions in parameter space uniformly.</strong> Another major drawback of GD is that unlike Newton’s method, the learning rate for GD is the same in all directions in parameter space. For this reason, the maximum learning rate is set by the behavior of the steepest direction and this can significantly slow down training. Ideally, we would like to take large steps in flat directions and small steps in steep directions. Since we are exploring rugged landscapes where curvatures change, this requires us to keep track of not only the gradient but second derivatives. The ideal scenario would be to calculate the Hessian but this proves to be too computationally expensive.</p></li>
<li><p>GD can take exponential time to escape saddle points, even with random initialization. As we mentioned, GD is extremely sensitive to initial condition since it determines the particular local minimum GD would eventually reach. However, even with a good initialization scheme, through the introduction of randomness, GD can still take exponential time to escape saddle points.</p></li>
</ul>
</section>
<section id="material-for-lab-sessions-sessions-tuesday-and-wednesday">
<h2>Material for lab sessions sessions Tuesday and Wednesday<a class="headerlink" href="#material-for-lab-sessions-sessions-tuesday-and-wednesday" title="Link to this heading">#</a></h2>
<p>The material here contains a summary of the lecture on Monday and discussion of SVD, Ridge and Lasso regression with examples</p>
</section>
<section id="linear-regression-and-the-svd">
<h2>Linear Regression and the SVD<a class="headerlink" href="#linear-regression-and-the-svd" title="Link to this heading">#</a></h2>
<p>We used the SVD to analyse the matrix to invert in ordinary lineat regression</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{X}^T\boldsymbol{X}=\boldsymbol{V}\boldsymbol{\Sigma}^T\boldsymbol{U}^T\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^T=\boldsymbol{V}\boldsymbol{\Sigma}^T\boldsymbol{\Sigma}\boldsymbol{V}^T.
\]</div>
<p>Since the matrices here have dimension <span class="math notranslate nohighlight">\(p\times p\)</span>, with <span class="math notranslate nohighlight">\(p\)</span> corresponding to the singular values, we defined last week the matrix</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\boldsymbol{\Sigma}^T\boldsymbol{\Sigma} = \begin{bmatrix} \tilde{\boldsymbol{\Sigma}} &amp; \boldsymbol{0}\\ \end{bmatrix}\begin{bmatrix} \tilde{\boldsymbol{\Sigma}} \\ \boldsymbol{0}\end{bmatrix},
\end{split}\]</div>
<p>where the tilde-matrix <span class="math notranslate nohighlight">\(\tilde{\boldsymbol{\Sigma}}\)</span> is a matrix of dimension <span class="math notranslate nohighlight">\(p\times p\)</span> containing only the singular values <span class="math notranslate nohighlight">\(\sigma_i\)</span>, that is</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\tilde{\boldsymbol{\Sigma}}=\begin{bmatrix} \sigma_0 &amp; 0 &amp; 0 &amp; \dots &amp; 0 &amp; 0 \\
0 &amp; \sigma_1 &amp; 0 &amp; \dots &amp; 0 &amp; 0 \\
0 &amp; 0 &amp; \sigma_2 &amp; \dots &amp; 0 &amp; 0 \\
0 &amp; 0 &amp; 0 &amp; \dots &amp; \sigma_{p-2} &amp; 0 \\
0 &amp; 0 &amp; 0 &amp; \dots &amp; 0 &amp; \sigma_{p-1} \\
\end{bmatrix},
\end{split}\]</div>
<p>meaning we can write</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{X}^T\boldsymbol{X}=\boldsymbol{V}\tilde{\boldsymbol{\Sigma}}^2\boldsymbol{V}^T.
\]</div>
<p>Multiplying from the right with <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span> (using the orthogonality of <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span>) we get</p>
<div class="math notranslate nohighlight">
\[
\left(\boldsymbol{X}^T\boldsymbol{X}\right)\boldsymbol{V}=\boldsymbol{V}\tilde{\boldsymbol{\Sigma}}^2.
\]</div>
</section>
<section id="what-does-it-mean">
<h2>What does it mean?<a class="headerlink" href="#what-does-it-mean" title="Link to this heading">#</a></h2>
<p>This means the vectors <span class="math notranslate nohighlight">\(\boldsymbol{v}_i\)</span> of the orthogonal matrix <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span>
are the eigenvectors of the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span> with eigenvalues
given by the singular values squared, that is</p>
<div class="math notranslate nohighlight">
\[
\left(\boldsymbol{X}^T\boldsymbol{X}\right)\boldsymbol{v}_i=\boldsymbol{v}_i\sigma_i^2.
\]</div>
<p>In other words, each non-zero singular value of <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> is a positive
square root of an eigenvalue of <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span>. It means also that
the columns of <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span> are the eigenvectors of
<span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span>. Since we have ordered the singular values of
<span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> in a descending order, it means that the column vectors
<span class="math notranslate nohighlight">\(\boldsymbol{v}_i\)</span> are hierarchically ordered by how much correlation they
encode from the columns of <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>.</p>
<p>Note that these are also the eigenvectors and eigenvalues of the
Hessian matrix.</p>
</section>
<section id="id2">
<h2>Ridge and LASSO Regression<a class="headerlink" href="#id2" 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>
</section>
<section id="from-ols-to-ridge-and-lasso">
<h2>From OLS to Ridge and Lasso<a class="headerlink" href="#from-ols-to-ridge-and-lasso" title="Link to this heading">#</a></h2>
<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. We do not include such a constraints in the discussions here.</p>
<p>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="id3">
<h2>Deriving the Ridge Regression Equations<a class="headerlink" href="#id3" 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>
</section>
<section id="note-on-scikit-learn">
<h2>Note on Scikit-Learn<a class="headerlink" href="#note-on-scikit-learn" title="Link to this heading">#</a></h2>
<p>Note well that a library like <strong>Scikit-Learn</strong> does not include the <span class="math notranslate nohighlight">\(1/n\)</span> factor in the expression for the mean-squared error. If you include it, the optimal parameter <span class="math notranslate nohighlight">\(\theta\)</span> becomes</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 our codes where we compare our own codes with <strong>Scikit-Learn</strong>, we do thus not include the <span class="math notranslate nohighlight">\(1/n\)</span> factor in the cost function.</p>
</section>
<section id="comparison-with-ols">
<h2>Comparison with OLS<a class="headerlink" href="#comparison-with-ols" title="Link to this heading">#</a></h2>
<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>
</section>
<section id="svd-analysis">
<h2>SVD analysis<a class="headerlink" href="#svd-analysis" title="Link to this heading">#</a></h2>
<p>Using our insights about the SVD of the design matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>
We have already analyzed the OLS solutions in terms of the eigenvectors (the columns) of the right singular value matrix <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> as</p>
<div class="math notranslate nohighlight">
\[
\tilde{\boldsymbol{y}}_{\mathrm{OLS}}=\boldsymbol{X}\boldsymbol{\theta} =\boldsymbol{U}\boldsymbol{U}^T\boldsymbol{y}.
\]</div>
<p>For Ridge regression this becomes</p>
<div class="math notranslate nohighlight">
\[
\tilde{\boldsymbol{y}}_{\mathrm{Ridge}}=\boldsymbol{X}\boldsymbol{\theta}_{\mathrm{Ridge}} = \boldsymbol{U\Sigma V^T}\left(\boldsymbol{V}\boldsymbol{\Sigma}^2\boldsymbol{V}^T+\lambda\boldsymbol{I} \right)^{-1}(\boldsymbol{U\Sigma V^T})^T\boldsymbol{y}=\sum_{j=0}^{p-1}\boldsymbol{u}_j\boldsymbol{u}_j^T\frac{\sigma_j^2}{\sigma_j^2+\lambda}\boldsymbol{y},
\]</div>
<p>with the vectors <span class="math notranslate nohighlight">\(\boldsymbol{u}_j\)</span> being the columns of <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> from the SVD of the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>.</p>
</section>
<section id="id4">
<h2>Interpreting the Ridge results<a class="headerlink" href="#id4" title="Link to this heading">#</a></h2>
<p>Since <span class="math notranslate nohighlight">\(\lambda \geq 0\)</span>, it means that compared to OLS, we have</p>
<div class="math notranslate nohighlight">
\[
\frac{\sigma_j^2}{\sigma_j^2+\lambda} \leq 1.
\]</div>
<p>Ridge regression finds the coordinates of <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> with respect to the
orthonormal basis <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span>, it then shrinks the coordinates by
<span class="math notranslate nohighlight">\(\frac{\sigma_j^2}{\sigma_j^2+\lambda}\)</span>. Recall that the SVD has
eigenvalues ordered in a descending way, that is <span class="math notranslate nohighlight">\(\sigma_i \geq
\sigma_{i+1}\)</span>.</p>
<p>For small eigenvalues <span class="math notranslate nohighlight">\(\sigma_i\)</span> it means that their contributions become less important, a fact which can be used to reduce the number of degrees of freedom.</p>
</section>
<section id="id5">
<h2>More interpretations<a class="headerlink" href="#id5" title="Link to this heading">#</a></h2>
<p>For the sake of simplicity, let us assume that the design matrix is orthonormal, that is</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{X}^T\boldsymbol{X}=(\boldsymbol{X}^T\boldsymbol{X})^{-1} =\boldsymbol{I}.
\]</div>
<p>In this case the standard OLS results in</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{\theta}^{\mathrm{OLS}} = \boldsymbol{X}^T\boldsymbol{y}=\sum_{i=0}^{n-1}\boldsymbol{u}_i\boldsymbol{u}_i^T\boldsymbol{y},
\]</div>
<p>and</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{\theta}^{\mathrm{Ridge}} = \left(\boldsymbol{I}+\lambda\boldsymbol{I}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}=\left(1+\lambda\right)^{-1}\boldsymbol{\theta}^{\mathrm{OLS}},
\]</div>
<p>that is the Ridge estimator scales the OLS estimator by the inverse of a factor <span class="math notranslate nohighlight">\(1+\lambda\)</span>, and
the Ridge estimator converges to zero when the hyperparameter goes to
infinity.</p>
<p>We will come back to more interpreations after we have gone through some of the statistical analysis part.</p>
<p>For more discussions of Ridge and Lasso regression, <a class="reference external" href="https://arxiv.org/abs/1509.09169">Wessel van Wieringen’s</a> article is highly recommended.
Similarly, <a class="reference external" href="https://arxiv.org/abs/1803.08823">Mehta et al’s article</a> is also recommended.</p>
</section>
<section id="id6">
<h2>Deriving the Lasso Regression Equations<a class="headerlink" href="#id6" title="Link to this heading">#</a></h2>
<p>Using the matrix-vector expression for Lasso regression, we have the following <strong>cost</strong> function</p>
<div class="math notranslate nohighlight">
\[
C(\boldsymbol{X},\boldsymbol{\theta})=\frac{1}{n}\left\{(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta})^T(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta})\right\}+\lambda\vert\vert\boldsymbol{\theta}\vert\vert_1,
\]</div>
<p>Taking the derivative with respect to <span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span> and recalling that the derivative of the absolute value is (we drop the boldfaced vector symbol for simplicity)</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\frac{d \vert \theta\vert}{d \theta}=\mathrm{sgn}(\theta)=\left\{\begin{array}{cc} 1 &amp; \theta &gt; 0 \\-1 &amp; \theta &lt; 0, \end{array}\right.
\end{split}\]</div>
<p>we have that the derivative of the cost function is</p>
<div class="math notranslate nohighlight">
\[
\frac{\partial C(\boldsymbol{X},\boldsymbol{\theta})}{\partial \boldsymbol{\theta}}=-\frac{2}{n}\boldsymbol{X}^T(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta})+\lambda sgn(\boldsymbol{\theta})=0,
\]</div>
<p>and reordering we have</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\theta}+\lambda sgn(\boldsymbol{\theta})=\boldsymbol{X}^T\boldsymbol{y}.
\]</div>
<p>This equation does not lead to a nice analytical equation as in Ridge regression or ordinary least squares. We have absorbed the factor <span class="math notranslate nohighlight">\(2/n\)</span> in a redefinition of the parameter <span class="math notranslate nohighlight">\(\lambda\)</span>. We will solve this type of problems using libraries like <strong>scikit-learn</strong> and using our own gradient descent code in project 1.</p>
</section>
<section id="simple-example-to-illustrate-ordinary-least-squares-ridge-and-lasso-regression">
<h2>Simple example to illustrate Ordinary Least Squares, Ridge and Lasso Regression<a class="headerlink" href="#simple-example-to-illustrate-ordinary-least-squares-ridge-and-lasso-regression" title="Link to this heading">#</a></h2>
<p>Let us assume that our design matrix is given by unit (identity) matrix, that is a square diagonal matrix with ones only along the
diagonal. In this case we have an equal number of rows and columns <span class="math notranslate nohighlight">\(n=p\)</span>.</p>
<p>Our model approximation is just <span class="math notranslate nohighlight">\(\tilde{\boldsymbol{y}}=\boldsymbol{\theta}\)</span> and the mean squared error and thereby the cost function for ordinary least sqquares (OLS) is then (we drop the term <span class="math notranslate nohighlight">\(1/n\)</span>)</p>
<div class="math notranslate nohighlight">
\[
C(\boldsymbol{\theta})=\sum_{i=0}^{p-1}(y_i-\theta_i)^2,
\]</div>
<p>and minimizing we have that</p>
<div class="math notranslate nohighlight">
\[
\hat{\theta}_i^{\mathrm{OLS}} = y_i.
\]</div>
</section>
<section id="ridge-regression">
<h2>Ridge Regression<a class="headerlink" href="#ridge-regression" title="Link to this heading">#</a></h2>
<p>For Ridge regression our cost function is</p>
<div class="math notranslate nohighlight">
\[
C(\boldsymbol{\theta})=\sum_{i=0}^{p-1}(y_i-\theta_i)^2+\lambda\sum_{i=0}^{p-1}\theta_i^2,
\]</div>
<p>and minimizing we have that</p>
<div class="math notranslate nohighlight">
\[
\hat{\theta}_i^{\mathrm{Ridge}} = \frac{y_i}{1+\lambda}.
\]</div>
</section>
<section id="lasso-regression">
<h2>Lasso Regression<a class="headerlink" href="#lasso-regression" title="Link to this heading">#</a></h2>
<p>For Lasso regression our cost function is</p>
<div class="math notranslate nohighlight">
\[
C(\boldsymbol{\theta})=\sum_{i=0}^{p-1}(y_i-\theta_i)^2+\lambda\sum_{i=0}^{p-1}\vert\theta_i\vert=\sum_{i=0}^{p-1}(y_i-\theta_i)^2+\lambda\sum_{i=0}^{p-1}\sqrt{\theta_i^2},
\]</div>
<p>and minimizing we have that</p>
<div class="math notranslate nohighlight">
\[
-2\sum_{i=0}^{p-1}(y_i-\theta_i)+\lambda \sum_{i=0}^{p-1}\frac{(\theta_i)}{\vert\theta_i\vert}=0,
\]</div>
<p>which leads to</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\hat{\boldsymbol{\theta}}_i^{\mathrm{Lasso}} = \left\{\begin{array}{ccc}y_i-\frac{\lambda}{2} &amp;\mathrm{if} &amp; y_i&gt; \frac{\lambda}{2}\\
y_i+\frac{\lambda}{2} &amp;\mathrm{if} &amp; y_i&lt; -\frac{\lambda}{2}\\
0 &amp;\mathrm{if} &amp; \vert y_i\vert\le \frac{\lambda}{2}\end{array}\right.\\.
\end{split}\]</div>
<p>Plotting these results shows clearly that Lasso regression suppresses (sets to zero) values of <span class="math notranslate nohighlight">\(\theta_i\)</span> for specific values of <span class="math notranslate nohighlight">\(\lambda\)</span>. Ridge regression reduces on the other hand the values of <span class="math notranslate nohighlight">\(\theta_i\)</span> as function of <span class="math notranslate nohighlight">\(\lambda\)</span>.</p>
</section>
<section id="yet-another-example">
<h2>Yet another Example<a class="headerlink" href="#yet-another-example" title="Link to this heading">#</a></h2>
<p>Let us assume we have a data set with outputs/targets given by the vector</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\boldsymbol{y}=\begin{bmatrix}4 \\ 2 \\3\end{bmatrix},
\end{split}\]</div>
<p>and our inputs as a <span class="math notranslate nohighlight">\(3\times 2\)</span> design matrix</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\boldsymbol{X}=\begin{bmatrix}2 &amp; 0\\ 0 &amp; 1 \\ 0 &amp; 0\end{bmatrix},
\end{split}\]</div>
<p>meaning that we have two features and two unknown parameters <span class="math notranslate nohighlight">\(\theta_0\)</span> and <span class="math notranslate nohighlight">\(\theta_1\)</span> to be determined either by ordinary least squares, Ridge or Lasso regression.</p>
</section>
<section id="the-ols-case">
<h2>The OLS case<a class="headerlink" href="#the-ols-case" title="Link to this heading">#</a></h2>
<p>For ordinary least squares (OLS) we know that the optimal solution is</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>Inserting the above values we obtain that</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\hat{\boldsymbol{\theta}}^{\mathrm{OLS}}=\begin{bmatrix}2 \\ 2\end{bmatrix},
\end{split}\]</div>
<p>The code which implements this simpler case is presented after the discussion of Ridge and Lasso.</p>
</section>
<section id="the-ridge-case">
<h2>The Ridge case<a class="headerlink" href="#the-ridge-case" title="Link to this heading">#</a></h2>
<p>For Ridge regression we have</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>Inserting the above values we obtain that</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\hat{\boldsymbol{\theta}}^{\mathrm{Ridge}}=\begin{bmatrix}\frac{8}{4+\lambda} \\ \frac{2}{1+\lambda}\end{bmatrix},
\end{split}\]</div>
<p>There is normally a constraint on the value of <span class="math notranslate nohighlight">\(\vert\vert \boldsymbol{\theta}\vert\vert_2\)</span> via the parameter <span class="math notranslate nohighlight">\(\lambda\)</span>.
Let us for simplicity assume that <span class="math notranslate nohighlight">\(\theta_0^2+\theta_1^2=1\)</span> as constraint. This will allow us to find an expression for the optimal values of <span class="math notranslate nohighlight">\(\theta\)</span> and <span class="math notranslate nohighlight">\(\lambda\)</span>.</p>
<p>To see this, let us write the cost function for Ridge regression.</p>
</section>
<section id="writing-the-cost-function">
<h2>Writing the Cost Function<a class="headerlink" href="#writing-the-cost-function" title="Link to this heading">#</a></h2>
<p>We define the MSE without the <span class="math notranslate nohighlight">\(1/n\)</span> factor and have then, using that</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\boldsymbol{X}\boldsymbol{\theta}=\begin{bmatrix} 2\theta_0 \\ \theta_1 \\0 \end{bmatrix},
\end{split}\]</div>
<div class="math notranslate nohighlight">
\[
C(\boldsymbol{\theta})=(4-2\theta_0)^2+(2-\theta_1)^2+\lambda(\theta_0^2+\theta_1^2),
\]</div>
<p>and taking the derivative with respect to <span class="math notranslate nohighlight">\(\theta_0\)</span> we get</p>
<div class="math notranslate nohighlight">
\[
\theta_0=\frac{8}{4+\lambda},
\]</div>
<p>and for <span class="math notranslate nohighlight">\(\theta_1\)</span> we obtain</p>
<div class="math notranslate nohighlight">
\[
\theta_1=\frac{2}{1+\lambda},
\]</div>
<p>Using the constraint for <span class="math notranslate nohighlight">\(\theta_0^2+\theta_1^2=1\)</span> we can constrain <span class="math notranslate nohighlight">\(\lambda\)</span> by solving</p>
<div class="math notranslate nohighlight">
\[
\left(\frac{8}{4+\lambda}\right)^2+\left(\frac{2}{1+\lambda}\right)^2=1,
\]</div>
<p>which gives <span class="math notranslate nohighlight">\(\lambda=4.571\)</span> and <span class="math notranslate nohighlight">\(\theta_0=0.933\)</span> and <span class="math notranslate nohighlight">\(\theta_1=0.359\)</span>.</p>
</section>
<section id="lasso-case">
<h2>Lasso case<a class="headerlink" href="#lasso-case" title="Link to this heading">#</a></h2>
<p>For Lasso we need now, keeping a constraint on <span class="math notranslate nohighlight">\(\vert\theta_0\vert+\vert\theta_1\vert=1\)</span>, to take the derivative of the absolute values of <span class="math notranslate nohighlight">\(\theta_0\)</span>
and <span class="math notranslate nohighlight">\(\theta_1\)</span>. This gives us the following derivatives of the cost function</p>
<div class="math notranslate nohighlight">
\[
C(\boldsymbol{\theta})=(4-2\theta_0)^2+(2-\theta_1)^2+\lambda(\vert\theta_0\vert+\vert\theta_1\vert),
\]</div>
<div class="math notranslate nohighlight">
\[
\frac{\partial C(\boldsymbol{\theta})}{\partial \theta_0}=-4(4-2\theta_0)+\lambda\mathrm{sgn}(\theta_0)=0,
\]</div>
<p>and</p>
<div class="math notranslate nohighlight">
\[
\frac{\partial C(\boldsymbol{\theta})}{\partial \theta_1}=-2(2-\theta_1)+\lambda\mathrm{sgn}(\theta_1)=0.
\]</div>
<p>We have now four cases to solve besides the trivial cases <span class="math notranslate nohighlight">\(\theta_0\)</span> and/or <span class="math notranslate nohighlight">\(\theta_1\)</span> are zero, namely</p>
<ol class="arabic simple">
<li><p><span class="math notranslate nohighlight">\(\theta_0 &gt; 0\)</span> and <span class="math notranslate nohighlight">\(\theta_1 &gt; 0\)</span>,</p></li>
<li><p><span class="math notranslate nohighlight">\(\theta_0 &gt; 0\)</span> and <span class="math notranslate nohighlight">\(\theta_1 &lt; 0\)</span>,</p></li>
<li><p><span class="math notranslate nohighlight">\(\theta_0 &lt; 0\)</span> and <span class="math notranslate nohighlight">\(\theta_1 &gt; 0\)</span>,</p></li>
<li><p><span class="math notranslate nohighlight">\(\theta_0 &lt; 0\)</span> and <span class="math notranslate nohighlight">\(\theta_1 &lt; 0\)</span>.</p></li>
</ol>
</section>
<section id="the-first-case">
<h2>The first Case<a class="headerlink" href="#the-first-case" title="Link to this heading">#</a></h2>
<p>If we consider the first case, we have then</p>
<div class="math notranslate nohighlight">
\[
-4(4-2\theta_0)+\lambda=0,
\]</div>
<p>and</p>
<div class="math notranslate nohighlight">
\[
-2(2-\theta_1)+\lambda=0.
\]</div>
<p>which yields</p>
<div class="math notranslate nohighlight">
\[
\theta_0=\frac{16+\lambda}{8},
\]</div>
<p>and</p>
<div class="math notranslate nohighlight">
\[
\theta_1=\frac{4+\lambda}{2}.
\]</div>
<p>Using the constraint on <span class="math notranslate nohighlight">\(\theta_0\)</span> and <span class="math notranslate nohighlight">\(\theta_1\)</span> we can then find the optimal value of <span class="math notranslate nohighlight">\(\lambda\)</span> for the different cases. We leave this as an exercise to you.</p>
</section>
<section id="simple-code-for-solving-the-above-problem">
<h2>Simple code for solving the above problem<a class="headerlink" href="#simple-code-for-solving-the-above-problem" title="Link to this heading">#</a></h2>
<p>Here we set up the OLS, Ridge and Lasso functionality in order to study the above example. Note that here we have opted for a set of values of <span class="math notranslate nohighlight">\(\lambda\)</span>, meaning that we need to perform a search in order to find the optimal values.</p>
<p>First we study and compare the OLS and Ridge results. The next code compares all three methods.</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span>import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
def R2(y_data, y_model):
return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
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.
X = np.array( [ [ 2, 0], [0, 1], [0,0]])
y = np.array( [4, 2, 3])
# matrix inversion to find beta
OLSbeta = np.linalg.inv(X.T @ X) @ X.T @ y
print(OLSbeta)
# and then make the prediction
ytildeOLS = X @ OLSbeta
print(&quot;Training MSE for OLS&quot;)
print(MSE(y,ytildeOLS))
ypredictOLS = X @ OLSbeta
# Repeat now for Ridge regression and various values of the regularization parameter
I = np.eye(2,2)
# Decide which values of lambda to use
nlambdas = 100
MSEPredict = np.zeros(nlambdas)
lambdas = np.logspace(-4, 4, nlambdas)
for i in range(nlambdas):
lmb = lambdas[i]
Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y
# print(Ridgebeta)
# and then make the prediction
ypredictRidge = X @ Ridgebeta
MSEPredict[i] = MSE(y,ypredictRidge)
# print(MSEPredict[i])
# Now plot the results
plt.figure()
plt.plot(np.log10(lambdas), MSEPredict, &#39;r--&#39;, label = &#39;MSE Ridge Train&#39;)
plt.xlabel(&#39;log10(lambda)&#39;)
plt.ylabel(&#39;MSE&#39;)
plt.legend()
plt.show()
</pre></div>
</div>
</div>
</div>
<p>We see here that we reach a plateau. What is actually happening?</p>
</section>
<section id="with-lasso-regression">
<h2>With Lasso Regression<a class="headerlink" href="#with-lasso-regression" 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>import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn import linear_model
def R2(y_data, y_model):
return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
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.
X = np.array( [ [ 2, 0], [0, 1], [0,0]])
y = np.array( [4, 2, 3])
# matrix inversion to find beta
OLSbeta = np.linalg.inv(X.T @ X) @ X.T @ y
print(OLSbeta)
# and then make the prediction
ytildeOLS = X @ OLSbeta
print(&quot;Training MSE for OLS&quot;)
print(MSE(y,ytildeOLS))
ypredictOLS = X @ OLSbeta
# Repeat now for Ridge regression and various values of the regularization parameter
I = np.eye(2,2)
# Decide which values of lambda to use
nlambdas = 100
MSERidgePredict = np.zeros(nlambdas)
MSELassoPredict = np.zeros(nlambdas)
lambdas = np.logspace(-4, 4, nlambdas)
for i in range(nlambdas):
lmb = lambdas[i]
Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y
print(Ridgebeta)
# and then make the prediction
ypredictRidge = X @ Ridgebeta
MSERidgePredict[i] = MSE(y,ypredictRidge)
RegLasso = linear_model.Lasso(lmb,fit_intercept=False)
RegLasso.fit(X,y)
ypredictLasso = RegLasso.predict(X)
print(RegLasso.coef_)
MSELassoPredict[i] = MSE(y,ypredictLasso)
# Now plot the results
plt.figure()
plt.plot(np.log10(lambdas), MSERidgePredict, &#39;r--&#39;, label = &#39;MSE Ridge Train&#39;)
plt.plot(np.log10(lambdas), MSELassoPredict, &#39;r--&#39;, label = &#39;MSE Lasso Train&#39;)
plt.xlabel(&#39;log10(lambda)&#39;)
plt.ylabel(&#39;MSE&#39;)
plt.legend()
plt.show()
</pre></div>
</div>
</div>
</div>
</section>
<section id="another-example-now-with-a-polynomial-fit">
<h2>Another Example, now with a polynomial fit<a class="headerlink" href="#another-example-now-with-a-polynomial-fit" 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>import os
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 R2(y_data, y_model):
return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
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)
x = np.random.rand(100)
y = 2.0+5*x*x+0.1*np.random.randn(100)
# number of features p (here degree of polynomial
p = 3
# The design matrix now as function of a given polynomial
X = np.zeros((len(x),p))
X[:,0] = 1.0
X[:,1] = x
X[:,2] = x*x
# 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
OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train
print(OLSbeta)
# and then make the prediction
ytildeOLS = X_train @ OLSbeta
print(&quot;Training MSE for OLS&quot;)
print(MSE(y_train,ytildeOLS))
ypredictOLS = X_test @ OLSbeta
print(&quot;Test MSE OLS&quot;)
print(MSE(y_test,ypredictOLS))
# Repeat now for Lasso and Ridge regression and various values of the regularization parameter
I = np.eye(p,p)
# Decide which values of lambda to use
nlambdas = 100
MSEPredict = np.zeros(nlambdas)
MSETrain = np.zeros(nlambdas)
MSELassoPredict = np.zeros(nlambdas)
MSELassoTrain = np.zeros(nlambdas)
lambdas = np.logspace(-4, 4, nlambdas)
for i in range(nlambdas):
lmb = lambdas[i]
Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train
# include lasso using Scikit-Learn
RegLasso = linear_model.Lasso(lmb,fit_intercept=False)
RegLasso.fit(X_train,y_train)
# and then make the prediction
ytildeRidge = X_train @ Ridgebeta
ypredictRidge = X_test @ Ridgebeta
ytildeLasso = RegLasso.predict(X_train)
ypredictLasso = RegLasso.predict(X_test)
MSEPredict[i] = MSE(y_test,ypredictRidge)
MSETrain[i] = MSE(y_train,ytildeRidge)
MSELassoPredict[i] = MSE(y_test,ypredictLasso)
MSELassoTrain[i] = MSE(y_train,ytildeLasso)
# Now plot the results
plt.figure()
plt.plot(np.log10(lambdas), MSETrain, label = &#39;MSE Ridge train&#39;)
plt.plot(np.log10(lambdas), MSEPredict, &#39;r--&#39;, label = &#39;MSE Ridge Test&#39;)
plt.plot(np.log10(lambdas), MSELassoTrain, label = &#39;MSE Lasso train&#39;)
plt.plot(np.log10(lambdas), MSELassoPredict, &#39;r--&#39;, label = &#39;MSE Lasso Test&#39;)
plt.xlabel(&#39;log10(lambda)&#39;)
plt.ylabel(&#39;MSE&#39;)
plt.legend()
plt.show()
</pre></div>
</div>
</div>
</div>
</section>
</section>
<script type="text/x-thebe-config">
{
requestKernel: true,
binderOptions: {
repo: "binder-examples/jupyter-stacks-datascience",
ref: "master",
},
codeMirrorConfig: {
theme: "abcdef",
mode: "python"
},
kernelOptions: {
name: "python3",
path: "./."
},
predefinedOutput: true
}
</script>
<script>kernelName = 'python3'</script>
</article>
<footer class="prev-next-footer d-print-none">
<div class="prev-next-area">
<a class="left-prev"
href="exercisesweek36.html"
title="previous page">
<i class="fa-solid fa-angle-left"></i>
<div class="prev-next-info">
<p class="prev-next-subtitle">previous</p>
<p class="prev-next-title">Exercises week 36</p>
</div>
</a>
<a class="right-next"
href="exercisesweek37.html"
title="next page">
<div class="prev-next-info">
<p class="prev-next-subtitle">next</p>
<p class="prev-next-title">Exercises week 37</p>
</div>
<i class="fa-solid fa-angle-right"></i>
</a>
</div>
</footer>
</div>
<div class="bd-sidebar-secondary bd-toc"><div class="sidebar-secondary-items sidebar-secondary__inner">
<div class="sidebar-secondary-item">
<div class="page-toc tocsection onthispage">
<i class="fa-solid fa-list"></i> Contents
</div>
<nav class="bd-toc-nav page-toc">
<ul class="visible nav section-nav flex-column">
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#plans-for-week-36">Plans for week 36</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#material-for-lecture-monday-september-2">Material for lecture Monday September 2</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#mathematical-interpretation-of-ordinary-least-squares">Mathematical Interpretation of Ordinary Least Squares</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#residual-error">Residual Error</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#simple-case">Simple case</a></li>
<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>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#codes-for-the-svd">Codes for the SVD</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#note-about-svd-calculations">Note about SVD Calculations</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#mathematics-of-the-svd-and-implications">Mathematics of the SVD and implications</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#example-matrix">Example Matrix</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#setting-up-the-matrix-to-be-inverted">Setting up the Matrix to be inverted</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#further-properties-important-for-our-analyses-later">Further properties (important for our analyses later)</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>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#optimization-and-gradient-descent-the-central-part-of-any-machine-learning-algortithm">Optimization and gradient descent, the central part of any Machine Learning algortithm</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#reminder-on-newton-raphson-s-method">Reminder on Newton-Raphson’s method</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-equations">The equations</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#simple-geometric-interpretation">Simple geometric interpretation</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#extending-to-more-than-one-variable">Extending to more than one variable</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#steepest-descent">Steepest descent</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#more-on-steepest-descent">More on Steepest descent</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-ideal">The ideal</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-sensitiveness-of-the-gradient-descent">The sensitiveness of the gradient descent</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#convex-functions">Convex functions</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#convex-function">Convex function</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#conditions-on-convex-functions">Conditions on convex functions</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#more-on-convex-functions">More on convex functions</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#some-simple-problems">Some simple problems</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#revisiting-ordinary-least-squares">Revisiting Ordinary Least Squares</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#gradient-descent-example">Gradient descent example</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-derivative-of-the-cost-loss-function">The derivative of the cost/loss function</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-hessian-matrix">The Hessian matrix</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#simple-program">Simple program</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#id1">Gradient Descent Example</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#gradient-descent-and-ridge">Gradient descent and Ridge</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-hessian-matrix-for-ridge-regression">The Hessian matrix for Ridge Regression</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#program-example-for-gradient-descent-with-ridge-regression">Program example for gradient descent with Ridge Regression</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#using-gradient-descent-methods-limitations">Using gradient descent methods, limitations</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#material-for-lab-sessions-sessions-tuesday-and-wednesday">Material for lab sessions sessions Tuesday and Wednesday</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#linear-regression-and-the-svd">Linear Regression and 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="#id2">Ridge and LASSO Regression</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#from-ols-to-ridge-and-lasso">From OLS to Ridge and Lasso</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#id3">Deriving the Ridge Regression Equations</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#note-on-scikit-learn">Note on Scikit-Learn</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#comparison-with-ols">Comparison with OLS</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#svd-analysis">SVD analysis</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#id4">Interpreting the Ridge results</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#id5">More interpretations</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#id6">Deriving the Lasso Regression Equations</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#simple-example-to-illustrate-ordinary-least-squares-ridge-and-lasso-regression">Simple example to illustrate Ordinary Least Squares, Ridge and Lasso Regression</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#ridge-regression">Ridge Regression</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#lasso-regression">Lasso Regression</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#yet-another-example">Yet another Example</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-ols-case">The OLS case</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-ridge-case">The Ridge case</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#writing-the-cost-function">Writing the Cost Function</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#lasso-case">Lasso case</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-first-case">The first Case</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#simple-code-for-solving-the-above-problem">Simple code for solving the above problem</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#with-lasso-regression">With Lasso Regression</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#another-example-now-with-a-polynomial-fit">Another Example, now with a polynomial fit</a></li>
</ul>
</nav></div>
</div></div>
</div>
<footer class="bd-footer-content">
<div class="bd-footer-content__inner container">
<div class="footer-item">
<p class="component-author">
By Morten Hjorth-Jensen
</p>
</div>
<div class="footer-item">
<p class="copyright">
© Copyright 2023.
<br/>
</p>
</div>
<div class="footer-item">
</div>
<div class="footer-item">
</div>
</div>
</footer>
</main>
</div>
</div>
<!-- Scripts loaded after <body> so the DOM is not blocked -->
<script src="_static/scripts/bootstrap.js?digest=dfe6caa3a7d634c4db9b"></script>
<script src="_static/scripts/pydata-sphinx-theme.js?digest=dfe6caa3a7d634c4db9b"></script>
<footer class="bd-footer">
</footer>
</body>
</html>
Reference in New Issue View Git Blame Copy Permalink