{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Julia Version 1.6.2\n", "Commit 1b93d53fc4 (2021-07-14 15:36 UTC)\n", "Platform Info:\n", " OS: macOS (x86_64-apple-darwin18.7.0)\n", " CPU: Intel(R) Core(TM) i5-8279U CPU @ 2.40GHz\n", " WORD_SIZE: 64\n", " LIBM: libopenlibm\n", " LLVM: libLLVM-11.0.1 (ORCJIT, skylake)\n" ] } ], "source": [ "versioninfo()" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[32m\u001b[1m Activating\u001b[22m\u001b[39m environment at `~/Dropbox/class/M1399.000200/2021/M1399_000200-2021fall/Project.toml`\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\u001b[32m\u001b[1m Status\u001b[22m\u001b[39m `~/Dropbox/class/M1399.000200/2021/M1399_000200-2021fall/Project.toml`\n", " \u001b[90m [7d9fca2a] \u001b[39mArpack v0.5.3\n", " \u001b[90m [6e4b80f9] \u001b[39mBenchmarkTools v1.1.4\n", " \u001b[90m [1e616198] \u001b[39mCOSMO v0.8.1\n", " \u001b[90m [f65535da] \u001b[39mConvex v0.14.13\n", " \u001b[90m [a93c6f00] \u001b[39mDataFrames v1.2.2\n", " \u001b[90m [31a5f54b] \u001b[39mDebugger v0.6.8\n", " \u001b[90m [31c24e10] \u001b[39mDistributions v0.24.18\n", " \u001b[90m [e2685f51] \u001b[39mECOS v0.12.3\n", " \u001b[90m [f6369f11] \u001b[39mForwardDiff v0.10.19\n", " \u001b[90m [28b8d3ca] \u001b[39mGR v0.58.1\n", " \u001b[90m [c91e804a] \u001b[39mGadfly v1.3.3\n", " \u001b[90m [bd48cda9] \u001b[39mGraphRecipes v0.5.7\n", " \u001b[90m [82e4d734] \u001b[39mImageIO v0.5.8\n", " \u001b[90m [6218d12a] \u001b[39mImageMagick v1.2.1\n", " \u001b[90m [916415d5] \u001b[39mImages v0.24.1\n", " \u001b[90m [b6b21f68] \u001b[39mIpopt v0.7.0\n", " \u001b[90m [42fd0dbc] \u001b[39mIterativeSolvers v0.9.1\n", " \u001b[90m [4076af6c] \u001b[39mJuMP v0.21.9\n", " \u001b[90m [b51810bb] \u001b[39mMatrixDepot v1.0.4\n", " \u001b[90m [1ec41992] \u001b[39mMosekTools v0.9.4\n", " \u001b[90m [76087f3c] \u001b[39mNLopt v0.6.3\n", " \u001b[90m [47be7bcc] \u001b[39mORCA v0.5.0\n", " \u001b[90m [a03496cd] \u001b[39mPlotlyBase v0.4.3\n", " \u001b[90m [f0f68f2c] \u001b[39mPlotlyJS v0.15.0\n", " \u001b[90m [91a5bcdd] \u001b[39mPlots v1.21.2\n", " \u001b[90m [438e738f] \u001b[39mPyCall v1.92.3\n", " \u001b[90m [d330b81b] \u001b[39mPyPlot v2.9.0\n", " \u001b[90m [dca85d43] \u001b[39mQuartzImageIO v0.7.3\n", " \u001b[90m [6f49c342] \u001b[39mRCall v0.13.12\n", " \u001b[90m [ce6b1742] \u001b[39mRDatasets v0.7.5\n", " \u001b[90m [c946c3f1] \u001b[39mSCS v0.7.1\n", " \u001b[90m [276daf66] \u001b[39mSpecialFunctions v1.6.1\n", " \u001b[90m [2913bbd2] \u001b[39mStatsBase v0.33.10\n", " \u001b[90m [b8865327] \u001b[39mUnicodePlots v2.0.1\n", " \u001b[90m [0f1e0344] \u001b[39mWebIO v0.8.15\n", " \u001b[90m [8f399da3] \u001b[39mLibdl\n", " \u001b[90m [2f01184e] \u001b[39mSparseArrays\n", " \u001b[90m [10745b16] \u001b[39mStatistics\n" ] } ], "source": [ "using Pkg\n", "Pkg.activate(\"../..\")\n", "Pkg.status()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# QR Decomposition\n", "\n", "* We learned Cholesky decomposition as **one** approach for solving linear regression.\n", "\n", "* Another approach for linear regression uses the QR decomposition. \n", " **This is how the `lm()` function in R does linear regression.**\n", " \\\n", " **This is also how Julia's (and MATLAB's) `\\` works for rectangular matrices.**" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Vector{Float64}:\n", " 0.3795466676698624\n", " 0.6508866456093487\n", " 0.39225041956535506" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Random\n", "\n", "Random.seed!(280) # seed\n", "\n", "n, p = 5, 3\n", "X = randn(n, p) # predictor matrix\n", "y = randn(n) # response vector\n", "\n", "# backslash finds the (minimum L2 norm) least squares solution\n", "X \\ y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We want to understand what is QR and how it is used for solving least squares problem." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Definitions\n", "\n", "* Assume $\\mathbf{X} \\in \\mathbb{R}^{n \\times p}$ has full column rank. Necessarilly $n \\ge p$.\n", "\n", "* **Full QR decomposition**: \n", "$$\n", " \\mathbf{X} = \\mathbf{Q} \\mathbf{R}, \n", "$$\n", "where \n", " - $\\mathbf{Q} \\in \\mathbb{R}^{n \\times n}$, $\\mathbf{Q}^T \\mathbf{Q} = \\mathbf{Q}\\mathbf{Q}^T = \\mathbf{I}_n$. In other words, $\\mathbf{Q}$ is an orthogonal matrix. \n", " - First $p$ columns of $\\mathbf{Q}$ form an orthonormal basis of ${\\cal R}(\\mathbf{X})$ (**range** or column space of $\\mathbf{X}$) \n", " - Last $n-p$ columns of $\\mathbf{Q}$ form an orthonormal basis of ${\\cal N}(\\mathbf{X}^T)$ (**null space** of $\\mathbf{X}^T$)\n", " - Recall that $\\mathcal{N}(\\mathbf{X}^T)=\\mathcal{R}(\\mathbf{X})^{\\perp}$ and $\\mathcal{R}(\\mathbf{X}) \\oplus \\mathcal{N}(\\mathbf{X}^T) = \\mathbb{R}^n$.\n", " - $\\mathbf{R} \\in \\mathbb{R}^{n \\times p}$ is upper triangular with positive diagonal entries. \n", " - The lower $(n-p)\\times p$ block of $\\mathbf{R}$ is $\\mathbf{0}$ (why?)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " * **Reduced QR decomposition**:\n", "$$\n", " \\mathbf{X} = \\mathbf{Q}_1 \\mathbf{R}_1,\n", "$$\n", "where\n", " - $\\mathbf{Q}_1 \\in \\mathbb{R}^{n \\times p}$, $\\mathbf{Q}_1^T \\mathbf{Q}_1 = \\mathbf{I}_p$. In other words, $\\mathbf{Q}_1$ is a partially orthogonal matrix. Note $\\mathbf{Q}_1\\mathbf{Q}_1^T \\neq \\mathbf{I}_n$.\n", " - $\\mathbf{R}_1 \\in \\mathbb{R}^{p \\times p}$ is an upper triangular matrix with positive diagonal entries. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Given QR decomposition $\\mathbf{X} = \\mathbf{Q} \\mathbf{R}$,\n", " $$\n", " \\mathbf{X}^T \\mathbf{X} = \\mathbf{R}^T \\mathbf{Q}^T \\mathbf{Q} \\mathbf{R} = \\mathbf{R}^T \\mathbf{R} = \\mathbf{R}_1^T \\mathbf{R}_1.\n", " $$\n", " - Once we have a (reduced) QR decomposition of $\\mathbf{X}$, we automatically have the Cholesky decomposition of the *Gram matrix* $\\mathbf{X}^T \\mathbf{X}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Application: least squares\n", "\n", "* Normal equation\n", "$$\n", " \\mathbf{X}^T\\mathbf{X}\\beta = \\mathbf{X}^T\\mathbf{y}\n", "$$\n", "is equivalently written with reduced QR as\n", "$$\n", " \\mathbf{R}_1^T\\mathbf{R}_1\\beta = \\mathbf{R}_1^T\\mathbf{Q}_1^T\\mathbf{y}\n", "$$\n", "\n", "* Since $\\mathbf{R}_1$ is invertible, we only need to solve the triangluar system\n", "$$\n", " \\mathbf{R}_1\\beta = \\mathbf{Q}_1^T\\mathbf{y}\n", "$$\n", "Multiplication $\\mathbf{Q}_1^T \\mathbf{y}$ is done implicitly (see below).\n", "\n", "* This method is numerically more stable than directly solving the normal equation, since $\\kappa(\\mathbf{X}^T\\mathbf{X}) = \\kappa(\\mathbf{X})^2$!\n", "\n", "* In case we need standard errors, compute inverse of $\\mathbf{R}_1^T \\mathbf{R}_1$. This involves triangular solves." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Gram-Schmidt procedure\n", "\n", "* Wait! Does $\\mathbf{X}$ always have a QR decomposition?\n", " - Yes. It is equivalent to the Gram-Schmidt procedure for basis orthonormalization.\n", "\n", "\n", "\n", "[Jørgen Pedersen Gram, 1850-1916](https://en.wikipedia.org/wiki/Jørgen_Pedersen_Gram)\n", "\n", "\n", "\n", "[Erhard Schmidt, 1876-1959](https://en.wikipedia.org/wiki/Erhard_Schmidt)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Assume $\\mathbf{X} = [\\mathbf{x}_1 | \\dotsb | \\mathbf{x}_p] \\in \\mathbb{R}^{n \\times p}$ has full column rank. That is, $\\mathbf{x}_1,\\ldots,\\mathbf{x}_p$ are *linearly independent*.\n", "\n", "* Gram-Schmidt (GS) procedure produces nested orthonormal basis vectors $\\{\\mathbf{q}_1, \\dotsc, \\mathbf{q}_p\\}$ that spans $\\mathcal{R}(\\mathbf{X})$, i.e.,\n", "$$\n", "\\begin{split}\n", " \\text{span}(\\{\\mathbf{x}_1\\}) &= \\text{span}(\\{\\mathbf{q}_1\\}) \\\\\n", " \\text{span}(\\{\\mathbf{x}_1, \\mathbf{x}_2\\}) &= \\text{span}(\\{\\mathbf{q}_1, \\mathbf{q}_2\\}) \\\\\n", " & \\vdots \\\\\n", " \\text{span}(\\{\\mathbf{x}_1, \\mathbf{x}_2, \\dotsc, \\mathbf{x}_p\\}) &= \\text{span}(\\{\\mathbf{q}_1, \\mathbf{q}_2, \\dotsc, \\mathbf{q}_p\\}) \n", "\\end{split}\n", "$$\n", "and $\\langle \\mathbf{q}_i, \\mathbf{q}_j \\rangle = \\delta_{ij}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* The algorithm:\n", " 0. Initialize $\\mathbf{q}_1 = \\mathbf{x}_1 / \\|\\mathbf{x}_1\\|_2$\n", " 0. For $k=2, \\ldots, p$, \n", "$$\n", "\\begin{eqnarray*}\n", "\t\\mathbf{v}_k &=& \\mathbf{x}_k - P_{\\text{span}(\\{\\mathbf{q}_1,\\ldots,\\mathbf{q}_{k-1}\\})}(\\mathbf{x}_k) = \\mathbf{x}_k - \\sum_{j=1}^{k-1} \\langle \\mathbf{q}_j, \\mathbf{x}_k \\rangle \\cdot \\mathbf{q}_j \\\\\n", "\t\\mathbf{q}_k &=& \\mathbf{v}_k / \\|\\mathbf{v}_k\\|_2\n", "\\end{eqnarray*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### GS conducts reduced QR\n", "\n", "* $\\mathbf{Q} = [\\mathbf{q}_1 | \\dotsb | \\mathbf{q}_p]$. Obviously $\\mathbf{Q}^T \\mathbf{Q} = \\mathbf{I}_p$.\n", "\n", "* Where is $\\mathbf{R}$? \n", " - Let $r_{jk} = \\langle \\mathbf{q}_j, \\mathbf{x}_k \\rangle$ for $j < k$, and $r_{kk} = \\|\\mathbf{v}_k\\|_2$.\n", " - Re-write the above expression:\n", "$$\n", " r_{kk} \\mathbf{q}_k = \\mathbf{x}_k - \\sum_{j=1}^{k-1} r_{jk} \\cdot \\mathbf{q}_j\n", "$$\n", " or\n", "$$\n", " \\mathbf{x}_k = r_{kk} \\mathbf{q}_k + \\sum_{j=1}^{k-1} r_{jk} \\cdot \\mathbf{q}_j\n", " .\n", "$$ \n", " - If we let $r_{jk} = 0$ for $j > k$, then $\\mathbf{R}=(r_{jk})$ is upper triangular and\n", "$$\n", " \\mathbf{X} = \\mathbf{Q}\\mathbf{R}\n", " .\n", "$$\n", " \n", " \n", " \n", "Source: \n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Classical Gram-Schmidt" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "cgs (generic function with 1 method)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using LinearAlgebra\n", "function cgs(X::Matrix{T}) where T<:AbstractFloat\n", " n, p = size(X)\n", " Q = Matrix{T}(undef, n, p)\n", " R = zeros(T, p, p)\n", " for j=1:p\n", " Q[:, j] .= X[:, j]\n", " for i=1:j-1\n", " R[i, j] = dot(Q[:, i], X[:, j])\n", " Q[:, j] .-= R[i, j] * Q[:, i]\n", " end\n", " R[j, j] = norm(Q[:, j])\n", " Q[:, j] /= R[j, j]\n", " end\n", " Q, R\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* CGS is *unstable* (we lose orthogonality due to roundoff errors) when columns of $\\mathbf{X}$ are almost collinear." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4×3 Matrix{Float32}:\n", " 1.0 1.0 1.0\n", " 1.19209f-7 0.0 0.0\n", " 0.0 1.19209f-7 0.0\n", " 0.0 0.0 1.19209f-7" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "e = eps(Float32)\n", "A = [1f0 1f0 1f0; e 0 0; 0 e 0; 0 0 e]" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4×3 Matrix{Float32}:\n", " 1.0 0.0 0.0\n", " 1.19209f-7 -0.707107 -0.707107\n", " 0.0 0.707107 0.0\n", " 0.0 0.0 0.707107" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Q, R = cgs(A)\n", "Q" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Matrix{Float32}:\n", " 1.0 -8.42937f-8 -8.42937f-8\n", " -8.42937f-8 1.0 0.5\n", " -8.42937f-8 0.5 1.0" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "transpose(Q)*Q" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* `Q` is hardly orthogonal.\n", "* Where exactly does the problem occur? (HW)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Modified Gram-Schmidt\n", "\n", "* The algorithm:\n", " 0. Initialize $\\mathbf{q}_1 = \\mathbf{x}_1 / \\|\\mathbf{x}_1\\|_2$\n", " 0. For $k=2, \\ldots, p$, \n", "$$\n", "\\begin{eqnarray*}\n", "\t\\mathbf{v}_k &=& \\mathbf{x}_k - P_{\\text{span}(\\{\\mathbf{q}_1,\\ldots,\\mathbf{q}_{k-1}\\})}(\\mathbf{x}_k) = \\mathbf{x}_k - \\sum_{j=1}^{k-1} \\langle \\mathbf{q}_j, \\mathbf{x}_k \\rangle \\cdot \\mathbf{q}_j \\\\\n", " &=& \\mathbf{x}_k - \\sum_{j=1}^{k-1} \\left\\langle \\mathbf{q}_j, \\mathbf{x}_k - \\sum_{l=1}^{j-1}\\langle \\mathbf{q}_l, \\mathbf{x}_k \\rangle \\mathbf{q}_l \\right\\rangle \\cdot \\mathbf{q}_j \\\\\n", "\t\\mathbf{q}_k &=& \\mathbf{v}_k / \\|\\mathbf{v}_k\\|_2\n", "\\end{eqnarray*}\n", "$$" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "mgs! (generic function with 1 method)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function mgs!(X::Matrix{T}) where T<:AbstractFloat\n", " n, p = size(X)\n", " R = zeros(T, p, p)\n", " for j=1:p\n", " for i=1:j-1\n", " R[i, j] = dot(X[:, i], X[:, j])\n", " X[:, j] -= R[i, j] * X[:, i]\n", " end\n", " R[j, j] = norm(X[:, j])\n", " X[:, j] /= R[j, j]\n", " end\n", " X, R\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* $\\mathbf{X}$ is overwritten by $\\mathbf{Q}$ and $\\mathbf{R}$ is stored in a separate array." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4×3 Matrix{Float32}:\n", " 1.0 0.0 0.0\n", " 1.19209f-7 -0.707107 -0.408248\n", " 0.0 0.707107 -0.408248\n", " 0.0 0.0 0.816497" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Q, R = mgs!(copy(A))\n", "Q" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Matrix{Float32}:\n", " 1.0 -8.42937f-8 -4.8667f-8\n", " -8.42937f-8 1.0 3.14007f-8\n", " -4.8667f-8 3.14007f-8 1.0" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "transpose(Q)*Q" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* So MGS is more stable than CGS. However, even MGS is not completely immune to instability." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2×2 Matrix{Float32}:\n", " 0.7 0.707107\n", " 0.7 0.707107" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "B = [0.7f0 0.7071068f0; 0.7000001f0 0.7071068f0]" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2×2 Matrix{Float32}:\n", " 0.707107 1.0\n", " 0.707107 0.0" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Q, R = mgs!(copy(B))\n", "Q" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2×2 Matrix{Float32}:\n", " 1.0 0.707107\n", " 0.707107 1.0" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "transpose(Q)*Q" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* `Q` is hardly orthogonal.\n", "* Where exactly the problem occurs? (HW)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Computational cost of CGS and MGS is $\\sum_{k=1}^p 4n(k-1) \\approx 2np^2$.\n", "\n", "* There are 3 algorithms to compute QR: (modified) Gram-Schmidt, Householder transform, (fast) Givens transform.\n", "\n", " In particular, the **Householder transform** for QR is implemented in LAPACK and thus used in R and Julia." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## QR by Householder transform\n", "\n", "\n", "\n", "[Alston Scott Householder (1904-1993)](https://en.wikipedia.org/wiki/Alston_Scott_Householder)\n", "\n", "* **This is the algorithm for solving linear regression in R**.\n", "\n", "* Assume again $\\mathbf{X} = [\\mathbf{x}_1 | \\dotsb | \\mathbf{x}_p] \\in \\mathbb{R}^{n \\times p}$ has full column rank.\n", "\n", "* Gram-Schmidt can be understood as:\n", "$$\n", " \\mathbf{X}\\mathbf{R}_{1} \\mathbf{R}_2 \\cdots \\mathbf{R}_n = \\mathbf{Q}_1\n", "$$\n", "where $\\mathbf{R}_j$ are a sequence of upper triangular matrices." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Householder QR does\n", "$$\n", " \\mathbf{H}_{p} \\cdots \\mathbf{H}_2 \\mathbf{H}_1 \\mathbf{X} = \\begin{pmatrix} \\mathbf{R}_1 \\\\ \\mathbf{0} \\end{pmatrix},\n", "$$\n", "where $\\mathbf{H}_j \\in \\mathbf{R}^{n \\times n}$ are a sequence of Householder transformation matrices.\n", "\n", " It yields the **full QR** where $\\mathbf{Q} = \\mathbf{H}_1 \\cdots \\mathbf{H}_p \\in \\mathbb{R}^{n \\times n}$. Recall that CGS/MGS only produces the **reduced QR** decomposition." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* For arbitrary $\\mathbf{v}, \\mathbf{w} \\in \\mathbb{R}^{n}$ with $\\|\\mathbf{v}\\|_2 = \\|\\mathbf{w}\\|_2$, we can construct a **Householder matrix** (or **Householder reflector**)\n", "$$\n", " \\mathbf{H} = \\mathbf{I}_n - 2 \\mathbf{u} \\mathbf{u}^T, \\quad \\mathbf{u} = - \\frac{1}{\\|\\mathbf{v} - \\mathbf{w}\\|_2} (\\mathbf{v} - \\mathbf{w}),\n", "$$\n", "that transforms $\\mathbf{v}$ to $\\mathbf{w}$:\n", "$$\n", "\t\\mathbf{H} \\mathbf{v} = \\mathbf{w}.\n", "$$\n", "$\\mathbf{H}$ is symmetric and orthogonal. Calculation of Householder vector $\\mathbf{u}$ costs $4n$ flops for general $\\mathbf{u}$ and $\\mathbf{w}$.\n", "\n", "\n", "\n", "Source: " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Now choose $\\mathbf{H}_1$ so that\n", "$$\n", "\t\\mathbf{H}_1 \\mathbf{x}_1 = \\begin{pmatrix} \\|\\mathbf{x}_{1}\\|_2 \\\\ 0 \\\\ \\vdots \\\\ 0 \\end{pmatrix}.\n", "$$\n", "That is, $\\mathbf{v} = \\mathbf{x}_1$ and $\\mathbf{w} = \\|\\mathbf{x}_1\\|_2\\mathbf{e}_1$.\n", "\n", "* Left-multiplying $\\mathbf{H}_1$ zeros out the first column of $\\mathbf{X}$ below (1, 1)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Take $\\mathbf{H}_2$ to zero the second column below diagonal:\n", "$$\n", "\\mathbf{H}_2\\mathbf{H}_1\\mathbf{X} = \n", "\\begin{bmatrix} \n", "\\times & \\times & \\times & \\times \\\\ \n", "0 & \\boldsymbol{\\times} & \\boldsymbol{\\times} & \\boldsymbol{\\times} \\\\\n", "0 & \\mathbf{0} & \\boldsymbol{\\times} & \\boldsymbol{\\times} \\\\\n", "0 & \\mathbf{0} & \\boldsymbol{\\times} & \\boldsymbol{\\times} \\\\\n", "0 & \\mathbf{0} & \\boldsymbol{\\times} & \\boldsymbol{\\times} \n", "\\end{bmatrix} \n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* In general, choose the $j$-th Householder transform $\\mathbf{H}_j = \\mathbf{I}_n - 2 \\mathbf{u}_j \\mathbf{u}_j^T$, where \n", "$$\n", " \\mathbf{u}_j = \\begin{bmatrix} \\mathbf{0}_{j-1} \\\\ {\\tilde u}_j \\end{bmatrix}, \\quad {\\tilde u}_j \\in \\mathbb{R}^{n-j+1},\n", "$$\n", "to zero the $j$-th column below diagonal. $\\mathbf{H}_j$ takes the form\n", "$$\n", "\t\\mathbf{H}_j = \\begin{bmatrix}\n", "\t\\mathbf{I}_{j-1} & \\\\\n", "\t& \\mathbf{I}_{n-j+1} - 2 {\\tilde u}_j {\\tilde u}_j^T\n", "\t\\end{bmatrix} = \\begin{bmatrix}\n", "\t\\mathbf{I}_{j-1} & \\\\\n", "\t& {\\tilde H}_{j}\n", "\t\\end{bmatrix}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Applying a Householder transform $\\mathbf{H} = \\mathbf{I} - 2 \\mathbf{u} \\mathbf{u}^T$ to a matrix $\\mathbf{X} \\in \\mathbb{R}^{n \\times p}$\n", "$$\n", "\t\\mathbf{H} \\mathbf{X} = \\mathbf{X} - 2 \\mathbf{u} (\\mathbf{u}^T \\mathbf{X})\n", "$$\n", "costs $4np$ flops. **Householder updates never entails explicit formation of the Householder matrices.**\n", "\n", "* Note applying ${\\tilde H}_j$ to $\\mathbf{X}$ only needs $4(n-j+1)(p-j+1)$ flops." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Algorithm\n", "\n", "```Julia\n", "for j=1:p\n", " u = House!(X[j:n, j])\n", " for i=j:p\n", " X[j:n, j:p] .-= 2u*(u'X[j:n, j:p])\n", " end\n", "end\n", "```\n", "\n", "* The process is done in place. Upper triangular part of $\\mathbf{X}$ is overwritten by $\\mathbf{R}_1$ and the essential Householder vectors ($\\tilde u_{j1}$ is normalized to 1) are stored in $\\mathbf{X}[j:n,j]$.\n", "\n", "* At $j$-th stage\n", " 0. computing the Householder vector ${\\tilde u}_j$ costs $3(n-j+1)$ flops\n", " 0. applying the Householder transform ${\\tilde H}_j$ to the $\\mathbf{X}[j:n, j:p]$ block costs $4(n-j+1)(p-j+1)$ flops \n", " \n", "* In total we need $\\sum_{j=1}^p [3(n-j+1) + 4(n-j+1)(p-j+1)] \\approx 2np^2 - \\frac 23 p^3$ flops." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Where is $\\mathbf{Q}$? \n", " - $\\mathbf{Q} = \\mathbf{H}_1 \\cdots \\mathbf{H}_p$. In some applications, it's necessary to form the orthogonal matrix $\\mathbf{Q}$. \n", "\n", " Accumulating $\\mathbf{Q}$ costs another $2np^2 - \\frac 23 p^3$ flops.\n", "\n", "* When computing $\\mathbf{Q}^T \\mathbf{v}$ or $\\mathbf{Q} \\mathbf{v}$ as in some applications (e.g., solve linear equation using QR), no need to form $\\mathbf{Q}$. Simply apply Householder transforms successively to the vector $\\mathbf{v}$. (HW)\n", "\n", "* Computational cost of Householder QR for linear regression: $2n p^2 - \\frac 23 p^3$ (regression coefficients and $\\hat \\sigma^2$) or more (fitted values, s.e., ...)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Householder QR with column pivoting\n", "\n", "Consider rank deficient $\\mathbf{X}$.\n", "\n", "* At the $j$-th stage, swap the column `X[:, j]` with `X[:, k]` where `k` is the column number in `X[j:n,j:p]` with maximum $\\ell_2$ norm to be the pivot column. If the maximum $\\ell_2$ norm is 0, it stops, ending with\n", "$$\n", "\\mathbf{X} \\mathbf{P} = \\mathbf{Q} \\begin{bmatrix} \\mathbf{R}_{11} & \\mathbf{R}_{12} \\\\ \\mathbf{0}_{(n-r) \\times r} & \\mathbf{0}_{(n-r) \\times (p-r)} \\end{bmatrix},\n", "$$\n", "where $\\mathbf{P} \\in \\mathbb{R}^{p \\times p}$ is a permutation matrix and $r$ is the rank of $\\mathbf{X}$. QR with column pivoting is rank revealing." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Implementation\n", "\n", "* Julia functions: [`qr`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.qr), [`qrfact!`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.qr!), or call LAPACK wrapper functions [`geqrf!`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.LAPACK.geqrf!) and [`geqp3!`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.LAPACK.geqp3!)\n", "\n", "* R function: [`qr`](https://stat.ethz.ch/R-manual/R-devel/library/base/html/qr.html). Wraps LAPACK routine [`dgeqp3`](http://www.netlib.org/lapack/explore-html/dd/d9a/group__double_g_ecomputational_ga1b0500f49e03d2771b797c6e88adabbb.html) (with `LAPACK=TRUE`; default uses LINPACK, an ancient version of LAPACK)." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "([0.12623784496408763 -1.1278278196026448 -0.8826696457022515; -2.3468794813834966 1.1478628422667982 1.7138424693341203; … ; -0.23920888512829233 -0.23706547458394342 1.0818212935057896; -0.5784508270929073 -0.6809935026697 -0.17040614729185025], [-1.794259674143309, 1.0913793110025305, 0.4266277597108536, -0.6244337204329091, 0.03204861737738283])" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X, y" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Vector{Float64}:\n", " 0.3795466676698624\n", " 0.6508866456093487\n", " 0.39225041956535506" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X \\ y # least squares solution by QR" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Vector{Float64}:\n", " 0.37954666766986195\n", " 0.6508866456093481\n", " 0.3922504195653549" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# same as\n", "qr(X) \\ y" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Vector{Float64}:\n", " 0.37954666766986256\n", " 0.6508866456093485\n", " 0.3922504195653555" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cholesky(X'X) \\ (X'y) # least squares solution by Cholesky" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "QRPivoted{Float64, Matrix{Float64}}\n", "Q factor:\n", "5×5 LinearAlgebra.QRPackedQ{Float64, Matrix{Float64}}:\n", " -0.0407665 -0.692007 0.318693 0.185526 -0.619257\n", " 0.757887 -0.0465712 -0.260086 -0.522259 -0.288166\n", " -0.618938 -0.233814 -0.404293 -0.624592 -0.0931611\n", " 0.0772486 -0.235405 -0.808135 0.53435 0.00216687\n", " 0.186801 -0.639431 0.119392 -0.131075 0.724429\n", "R factor:\n", "3×3 Matrix{Float64}:\n", " -3.09661 1.60888 1.84089\n", " 0.0 1.53501 0.556903\n", " 0.0 0.0 -1.32492\n", "permutation:\n", "3-element Vector{Int64}:\n", " 1\n", " 2\n", " 3" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# QR factorization with column pivoting\n", "xqr = qr(X, Val(true))" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Vector{Float64}:\n", " 0.3795466676698624\n", " 0.6508866456093487\n", " 0.39225041956535506" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "xqr \\ y # least squares solution" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.071020016095422e-15" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# thin Q matrix multiplication (a sequence of Householder transforms)\n", "norm(xqr.Q * xqr.R - X[:, xqr.p]) # recovers X (with columns permuted)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## QR by Givens rotation\n", "\n", "* Householder transform $\\mathbf{H}_j$ introduces batch of zeros into a vector.\n", "\n", "* **Givens transform** (aka **Givens rotation**, **Jacobi rotation**, **plane rotation**) selectively zeros one element of a vector.\n", "\n", "* Overall QR by Givens rotation is less efficient than the Householder method, but is better suited for matrices with structured patterns of nonzero elements.\n", "\n", "* **Givens/Jacobi rotations**: \n", "$$\n", "\t\\mathbf{G}(i,k,\\theta) = \\begin{bmatrix} \n", "\t1 & & 0 & & 0 & & 0 \\\\\n", "\t\\vdots & \\ddots & \\vdots & & \\vdots & & \\vdots \\\\\n", "\t0 & & c & & s & & 0 \\\\ \n", "\t\\vdots & & \\vdots & \\ddots & \\vdots & & \\vdots \\\\\n", "\t0 & & - s & & c & & 0 \\\\\n", "\t\\vdots & & \\vdots & & \\vdots & \\ddots & \\vdots \\\\\n", "\t0 & & 0 & & 0 & & 1 \\end{bmatrix},\n", "$$\n", "where $c = \\cos(\\theta)$ and $s = \\sin(\\theta)$. $\\mathbf{G}(i,k,\\theta)$ is orthogonal." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Pre-multiplication by $\\mathbf{G}(i,k,\\theta)^T$ rotates counterclockwise $\\theta$ radians in the $(i,k)$ coordinate plane. If $\\mathbf{x} \\in \\mathbb{R}^n$ and $\\mathbf{y} = \\mathbf{G}(i,k,\\theta)^T \\mathbf{x}$, then\n", "$$\n", "\ty_j = \\begin{cases}\n", "\tcx_i - s x_k & j = i \\\\\n", "\tsx_i + cx_k & j = k \\\\\n", "\tx_j & j \\ne i, k\n", "\t\\end{cases}.\n", "$$\n", "Apparently if we choose $\\tan(\\theta) = -x_k / x_i$, or equivalently,\n", "$$\n", "\\begin{eqnarray*}\n", "\tc = \\frac{x_i}{\\sqrt{x_i^2 + x_k^2}}, \\quad s = \\frac{-x_k}{\\sqrt{x_i^2 + x_k^2}},\n", "\\end{eqnarray*}\n", "$$\n", "then $y_k=0$.\n", "\n", "* Pre-applying Givens transform $\\mathbf{G}(i,k,\\theta)^T \\in \\mathbb{R}^{n \\times n}$ to a matrix $\\mathbf{A} \\in \\mathbb{R}^{n \\times m}$ only effects two rows of $\\mathbf{\n", "A}$:\n", "$$\n", "\t\\mathbf{A}([i, k], :) \\gets \\begin{bmatrix} c & s \\\\ -s & c \\end{bmatrix}^T \\mathbf{A}([i, k], :),\n", "$$\n", "costing $6m$ flops.\n", "\n", "* Post-applying Givens transform $\\mathbf{G}(i,k,\\theta) \\in \\mathbb{R}^{m \\times m}$ to a matrix $\\mathbf{A} \\in \\mathbb{R}^{n \\times m}$ only effects two columns of $\\mathbf{A}$:\n", "$$\n", "\t\\mathbf{A}(:, [i,k]) \\gets \\mathbf{A}(:, [i,k]) \\begin{bmatrix} c & s \\\\ -s & c \\end{bmatrix},\n", "$$\n", "costing $6n$ flops.\n", "\n", "* QR by Givens: $\\mathbf{G}_t^T \\cdots \\mathbf{G}_1^T \\mathbf{X} = \\begin{bmatrix} \\mathbf{R}_1 \\\\ \\mathbf{0} \\end{bmatrix}$.\n", "\n", "\n", "\n", "* Zeros in $\\mathbf{X}$ can also be introduced row-by-row.\n", "\n", "* If $\\mathbf{X} \\in \\mathbb{R}^{n \\times p}$, the total cost is $3np^2 - p^3$ flops and $O(np)$ square roots.\n", "\n", "* Note each Givens transform can be summarized by a single number, which is stored in the zeroed entry of $\\mathbf{X}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Applications\n", "\n", "### Linear regression\n", "\n", "* QR decomposition of $\\mathbf{X}$: $2np^2 - \\frac 23 p^3$ flops.\n", "\n", "* Solve $\\mathbf{R}^T \\mathbf{R} \\beta = \\mathbf{R}^T \\mathbf{Q}^T \\mathbf{y}$ for $\\beta$.\n", "\n", "* If $\\mathbf{X}$ is full rank, then $\\mathbf{R}$ is invertible, so we only need to solve the triangular system\n", "$$\n", " \\mathbf{R} \\beta = \\mathbf{Q}^T \\mathbf{y}\n", " .\n", "$$\n", "Multiplication $\\mathbf{Q}^T \\mathbf{y}$ is done implicitly.\n", "\n", "* If need standard errors, compute inverse of $\\mathbf{R}^T \\mathbf{R}$. This involves triangular solves." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Further reading\n", "\n", "* Section II.5.3 of [Computational Statistics](https://link.springer.com/book/10.1007%2F978-0-387-98144-4) by James Gentle (2010).\n", "\n", "* Chapter 5 of [Matrix Computation](https://www.amazon.com/Computations-Hopkins-Studies-Mathematical-Sciences/dp/1421407949/ref=sr_1_1?keywords=matrix+computation+golub&qid=1567157884&s=gateway&sr=8-1) by Gene Golub and Charles Van Loan (2013)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Acknowledgment\n", "\n", "Many parts of this lecture note is based on [Dr. Hua Zhou](http://hua-zhou.github.io)'s 2019 Spring Statistical Computing course notes available at ." ] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "kernelspec": { "display_name": "Julia 1.6.2", "language": "julia", "name": "julia-1.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.6.2" }, "toc": { "colors": { "hover_highlight": "#DAA520", "running_highlight": "#FF0000", "selected_highlight": "#FFD700" }, "moveMenuLeft": true, "nav_menu": { "height": "65px", "width": "252px" }, "navigate_menu": true, "number_sections": true, "sideBar": true, "skip_h1_title": true, "threshold": 4, "toc_cell": true, "toc_section_display": "block", "toc_window_display": true, "widenNotebook": false } }, "nbformat": 4, "nbformat_minor": 4 }