{ "cells": [ { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "***********************************\n", "Numerical Linear Algebra Background\n", "***********************************\n", "\n", "Solving Systems of Linear Equations — Matrix Factorizations\n", "===========================================================\n", "\n", "- The typical :math:`\\mathbf{A} \\mathbf{x} = \\mathbf{b}`, where\n", " :math:`\\mathbf{A} \\in \\mathbb{R}^{n \\times n}` and\n", " :math:`\\mathbf{b}, \\mathbf{x} \\in \\mathbb{R}^n`, can be solved using QR or\n", " PLU.\n", "\n", " - QR's orthogonal transformations (e.g. Householder) is numerically stable\n", " because they do not affect the magnitude of a matrix's norm.\n", " - PLU is accurate and half as expensive as QR, but could magnify the elements\n", " of :math:`\\mathbf{A}`.\n", "\n", "- Matrices with trivial analytical structure can be solved directly.\n", "\n", " - Permutation\n", " - Orthogonal\n", " - Nonsingular diagonal\n", " - Nonsingular diagonal block\n", " - Nonsingular lower triangular\n", " - Nonsingular upper triangular\n", "\n", "- Symmetric Positive Definite matrices can be factorized using Cholesky\n", " decomposition: :math:`\\mathbf{A} = \\mathbf{L} \\mathbf{L}^\\top` or\n", " :math:`\\mathbf{A} = \\mathbf{L} \\mathbf{D} \\mathbf{L}^\\top`.\n", "- Symmetric Indefinite matrices can be factorized using the Bunch-Parlett\n", " algorithm: :math:`\\mathbf{A} = \\mathbf{P} \\mathbf{L} \\mathbf{D}_B\n", " \\mathbf{L}^\\top \\mathbf{P}^\\top` where :math:`\\mathbf{P}` is a permutation\n", " matrix.\n", "\n", "Errors in Solving Linear Systems\n", "================================\n", "\n", "The condition number :math:`\\kappa_p(\\mathbf{A}) =\n", "\\lVert \\mathbf{A} \\rVert_p \\lVert \\mathbf{A}^{-1} \\rVert_p` of a matrix measures\n", "its sensitivity to perturbations. It could be used to indicate the matrix's\n", "nearness to singularity. Since computing the inverse is an expensive and\n", "unstable operation, the condition number is often estimated as\n", ":math:`\\hat{\\kappa}_p(\\mathbf{A}) = \\lVert \\mathbf{M} \\rVert_p \\alpha`.\n", ":math:`\\mathbf{M}` could be :math:`\\mathbf{A}` or a factor of :math:`\\mathbf{A}`\n", "with similar conditioning, and :math:`\\alpha` is an estimate of the induced\n", "matrix norm :math:`\\lVert \\mathbf{M}^{-1} \\rVert_p`.\n", "\n", "If ill-conditioned systems occur far away from the solution, then perturbing\n", "the system into a better-conditioned one is acceptable. If ill-conditioned\n", "systems occur near the solution, then the underlying problem is not well posed.\n", "\n", "Updating Matrix Factorizations\n", "==============================\n", "\n", "Jacobi rotation can be used to efficiently update a QR decomposition. Each\n", "Jacobi rotation zeros out one element of a matrix while changing only two rows\n", "of the matrix.\n", "\n", "Eigenvalues and Positive Definiteness\n", "=====================================\n", "\n", "Given a symmetric matrix, the Gerschgorin Theorem (3.5.9) provides a lower bound\n", "of the minimum eigenvalue and an upper bound for the maximum eigenvalue. If the\n", "matrix is strictly diagonally dominant, it is positive definite.\n", "\n", "Linear Least Squares\n", "====================\n", "\n", "The geometric interpretation of\n", ":math:`\\min_{\\mathbf{x}} \\lVert \\mathbf{A} \\mathbf{x} - \\mathbf{b} \\rVert_2` is\n", "given in Figure 3.6.2. The closest point to :math:`\\mathbf{b}` in subspace\n", ":math:`C(\\mathbf{A})` will be :math:`\\mathbf{A} \\mathbf{x}_* \\in C(\\mathbf{A})`\n", "such that the residual :math:`\\mathbf{A} \\mathbf{x}_* - \\mathbf{b}` is\n", "perpendicular to the entire subspace of :math:`C(\\mathbf{A})`. The solution\n", ":math:`\\mathbf{A}^\\top (\\mathbf{A} \\mathbf{x}_* - \\mathbf{b}) = \\boldsymbol{0}`\n", "is known as the normal equation because the residual vector is perpendicular to\n", "the column space of :math:`\\mathbf{A}` i.e. the residual is in the null space of\n", ":math:`\\mathbf{A}^\\top`. Note that SVD is more expensive than the other matrix\n", "factorization techniques, so only use it when :math:`\\mathbf{A}` is not known to\n", "have full row or column rank." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ ".. _dennis1996numerical-ex-3.1:\n", "\n", "Exercise 1\n", "==========\n", "\n", ".. math::\n", "\n", " \\lim_{p \\rightarrow \\infty} \\lVert x \\rVert_p\n", " &= \\lim_{p \\rightarrow \\infty} \\left(\n", " \\sum_{i = 1}^n \\left\\vert x_i \\right\\vert^p\n", " \\right)^{1 / p}\\\\\n", " &= \\lim_{p \\rightarrow \\infty} \\left(\n", " \\sum_i \\left\\vert \\frac{x_i}{M} \\right\\vert^p M^p\n", " \\right)^{1 / p}\n", " & \\quad & M = \\max_i \\left\\{ \\left\\vert x_i \\right\\vert \\right\\} =\n", " \\left\\vert x_m \\right\\vert\\\\\n", " &= \\lim_{p \\rightarrow \\infty} \\left(\n", " \\left\\vert \\frac{x_m}{M} \\right\\vert^p M^p +\n", " \\sum_{i \\neq m} \\left\\vert \\frac{x_i}{M} \\right\\vert^p M^p\n", " \\right)^{1 / p}\n", " & \\quad & \\left\\vert \\frac{x_i}{M} \\right\\vert^p < 1 \\forall\\, i \\neq m\n", " \\text{ and }\n", " \\left\\vert \\frac{x_m}{M} \\right\\vert^p = 1\\\\\n", " &= \\max_i \\left\\{ \\left\\vert x_i \\right\\vert \\right\\}\n", " & \\quad & \\lim_{p \\rightarrow \\infty}\n", " \\left\\vert \\frac{x_i}{M} \\right\\vert^p = \\delta(i - m)\\\\\n", " &= \\lVert x \\rVert_{\\infty}." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 2\n", "==========\n", "\n", "The following proof is from :cite:`olverlnn`.\n", "\n", "Recall that a norm defines a continuous real-valued function\n", ":math:`f(\\mathbf{v}) = \\lVert \\mathbf{v} \\rVert` on :math:`\\mathbb{R}^n`. Let\n", ":math:`\\left\\Vert \\cdot \\right\\Vert_1` and\n", ":math:`\\left\\Vert \\cdot \\right\\Vert_2` be any two norms on :math:`\\mathbb{R}^n`.\n", "Let :math:`S_1 = \\left\\{ \\left\\Vert \\mathbf{u} \\right\\Vert_1 = 1 \\right\\}`\n", "denote the unit sphere of the first norm. Any continuous function defined on a\n", "compact set achieves both a maximum and a minimum value. Thus, restricting the\n", "second norm function to the unit sphere :math:`S_1` of the first norm yields\n", "\n", ".. math::\n", "\n", " c^\\star = \\min \\left\\{\n", " \\left\\Vert \\mathbf{u} \\right\\Vert_2 \\mid \\mathbf{u} \\in S_1\n", " \\right\\}\n", " \\quad \\text{and} \\quad\n", " C^\\star = \\max \\left\\{\n", " \\left\\Vert \\mathbf{u} \\right\\Vert_2 \\mid \\mathbf{u} \\in S_1\n", " \\right\\}.\n", "\n", "By inspection, :math:`0 < c^\\star \\leq C^\\star < \\infty` with equality holding\n", "if and only if the norms are the same. Hence\n", "\n", ".. math::\n", "\n", " c^\\star \\leq \\left\\Vert \\mathbf{u} \\right\\Vert_2 \\leq C^\\star\n", " \\quad \\text{when} \\quad\n", " \\left\\Vert \\mathbf{u} \\right\\Vert_1 = 1.\n", "\n", "To generalize this result, define\n", ":math:`\\mathbf{u} = \\mathbf{v} / \\left\\Vert \\mathbf{v} \\right\\Vert_1 \\in S_1`.\n", "Consequently,\n", "\n", ".. math::\n", "\n", " \\left\\Vert \\mathbf{u} \\right\\Vert_2 =\n", " \\left\\Vert\n", " \\frac{\\mathbf{v}}{\\left\\Vert \\mathbf{v} \\right\\Vert_1}\n", " \\right\\Vert_2 =\n", " \\frac{\n", " \\left\\Vert \\mathbf{v} \\right\\Vert_2\n", " }{\n", " \\left\\Vert \\mathbf{v} \\right\\Vert_1\n", " }.\n", "\n", "Therefore,\n", "\n", ".. math::\n", "\n", " c^\\star \\left\\Vert \\mathbf{v} \\right\\Vert_1 \\leq\n", " \\left\\Vert \\mathbf{v} \\right\\Vert_2 \\leq\n", " C^\\star \\left\\Vert \\mathbf{v} \\right\\Vert_1." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 3\n", "==========\n", "\n", "(1)\n", "---\n", "\n", ".. math::\n", "\n", " \\lVert v \\rVert_1 \\geq \\lVert v \\rVert_2 \\geq \\lVert v \\rVert_\\infty\n", "\n", "because\n", "\n", ".. math::\n", "\n", " \\lVert v \\rVert_2^2 &= \\sum_{i = 1}^n \\lvert v_i \\rvert^2\\\\\n", " &\\leq \\sum_{i = 1}^n \\lvert v_i \\rvert \\sum_{j = 1}^n \\lvert v_j \\rvert\\\\\n", " &= \\lVert v \\rVert_1 \\lVert v \\rVert_1\\\\\n", " &= \\lVert v \\rVert_1^2\\\\\n", " \\lVert v \\rVert_2 &\\leq \\lVert v \\rVert_1\n", "\n", "and\n", "\n", ".. math::\n", "\n", " \\lVert v \\rVert_\\infty^2 &= \\left( \\max_{1 \\leq i \\leq n} v_i \\right)^2\\\\\n", " &\\leq v_\\max^2 + \\sum_{i \\neq \\max} v_i^2\\\\\n", " &= \\lVert v \\rVert_2^2\\\\\n", " \\lVert v \\rVert_\\infty &\\leq \\lVert v \\rVert_2.\n", "\n", "(2)\n", "---\n", "\n", ".. math::\n", "\n", " \\lVert v \\rVert_1^2 &= \\left( \\sum_{i = 1}^n \\lvert v_i \\rvert \\right)^2\\\\\n", " &= \\left( \\sum_{i = 1}^n \\lvert v_i \\rvert (1) \\right)^2\\\\\n", " &\\leq \\lVert v \\rVert_2^2 \\cdot \\lVert \\boldsymbol{1} \\rVert_2^2\n", " & \\quad & \\text{Cauchy–Schwarz inequality}\\\\\n", " &= n \\lVert v \\rVert_2^2\\\\\n", " \\lVert v \\rVert_1 &\\leq \\sqrt{n} \\lVert v \\rVert_2\n", "\n", "(3)\n", "---\n", "\n", ".. math::\n", "\n", " \\lVert v \\rVert_2^2 &= \\sum_{i = 1}^n \\lvert v_i \\rvert^2\\\\\n", " &\\leq \\sum_{i = 1}^n v_\\max^2\\\\\n", " &= n \\lVert v \\rVert_\\infty^2\\\\\n", " \\lVert v \\rVert_2 &\\leq \\sqrt{n} \\lVert v \\rVert_\\infty" ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ ".. _dennis1996numerical-ex-3.4:\n", "\n", "Exercise 4\n", "==========\n", "\n", "The :math:`p`-norm of a matrix :math:`A` induced by a given vector norm (3.1.7)\n", "is\n", "\n", ".. math::\n", "\n", " \\lVert A \\rVert_p\n", " &= \\max_{v \\in \\mathbb{R}^n \\setminus 0}\n", " \\left\\{\\frac{\\lVert A v \\rVert_p}{\\lVert v \\rVert_p} \\right\\}\\\\\n", " &= \\frac{\\lVert A v_* \\rVert_p}{\\lVert v_* \\rVert_p}\\\\\n", " &= \\frac{1}{\\lVert v_* \\rVert_p}\n", " \\left( \\sum_{i} \\lvert A_{i.} v_* \\rvert^p \\right)^{1 / p}\\\\\n", " &\\geq 0\n", "\n", "for every :math:`A \\in \\mathbb{R}^{n \\times n}` and\n", ":math:`\\lVert A \\rVert_p = 0` only if\n", ":math:`\\lVert A v_* \\rVert_p = 0`. An induced matrix norm is said to be\n", "compatible or consistent if\n", ":math:`\\lVert A v \\rVert_p \\leq \\lVert A \\rVert_p \\lVert v \\rVert_p`.\n", "\n", "(1)\n", "---\n", "\n", "Given that :math:`A \\in \\mathbb{R}^{n \\times n}` is invertible, the following\n", "property holds\n", "\n", ".. math::\n", "\n", " \\lVert A^{-1} \\rVert_p\n", " &= \\left(\n", " \\min_{v \\in \\mathbb{R}^n \\setminus 0}\n", " \\frac{\\lVert A v \\rVert_p}{\\lVert v \\rVert_p}\n", " \\right)^{-1}\n", " & \\quad & \\text{(3.1.18)}\\\\\n", " \\frac{1}{\\lVert A^{-1} \\rVert_p}\n", " &= \\min_{v \\in \\mathbb{R}^n \\setminus 0}\n", " \\frac{\\lVert A v \\rVert_p}{\\lVert v \\rVert_p}\\\\\n", " &\\leq \\max_{v \\in \\mathbb{R}^n \\setminus 0}\n", " \\frac{\\lVert A v \\rVert_p}{\\lVert v \\rVert_p}\\\\\n", " \\frac{\\lVert v_* \\rVert_p}{\\lVert A^{-1} \\rVert_p}\n", " &\\leq \\lVert A v_* \\rVert_p.\n", "\n", "(2)\n", "---\n", "\n", "For every :math:`\\alpha \\in \\mathbb{R}` and\n", ":math:`A \\in \\mathbb{R}^{n \\times n}`,\n", "\n", ".. math::\n", "\n", " \\lVert \\alpha A \\rVert_p\n", " &= \\max_{v \\in \\mathbb{R}^n \\setminus 0}\n", " \\left\\{\\frac{\\lVert \\alpha A v \\rVert_p}{\\lVert v \\rVert_p} \\right\\}\\\\\n", " &= \\frac{\\lVert \\alpha A v_* \\rVert_p}{\\lVert v_* \\rVert_p}\\\\\n", " &= \\frac{1}{\\lVert v_* \\rVert_p}\n", " \\left( \\sum_{i} \\lvert \\alpha A_{i.} v_* \\rvert^p \\right)^{1 / p}\\\\\n", " &= \\frac{\\lvert \\alpha \\rvert}{\\lVert v_* \\rVert_p}\n", " \\left( \\sum_{i} \\lvert A_{i.} v_* \\rvert^p \\right)^{1 / p}\\\\\n", " &= \\lvert \\alpha \\rvert \\lVert A \\rVert_p.\n", "\n", "(3)\n", "---\n", "\n", "For every :math:`A, B \\in \\mathbb{R}^{n \\times n}`,\n", "\n", ".. math::\n", "\n", " \\lVert A + B \\rVert_p\n", " &= \\max_{v \\in \\mathbb{R}^n \\setminus 0}\n", " \\left\\{\\frac{\\lVert (A + B) v \\rVert_p}{\\lVert v \\rVert_p} \\right\\}\\\\\n", " &= \\frac{\\lVert (A + B) v_* \\rVert_p}{\\lVert v_* \\rVert_p}\\\\\n", " &= \\frac{1}{\\lVert v_* \\rVert_p}\n", " \\left( \\sum_{i} \\lvert (A_{i.} + B_{i.}) v_* \\rvert^p \\right)^{1 / p}\\\\\n", " &\\leq \\frac{1}{\\lVert v_* \\rVert_p}\n", " \\left(\n", " \\left( \\sum_{i} \\lvert A_{i.} v_* \\rvert^p \\right)^{1 / p} +\n", " \\left( \\sum_{i} \\lvert B_{i.} v_* \\rvert^p \\right)^{1 / p}\n", " \\right)\n", " & \\quad & \\text{Minkowski inequality}\\\\\n", " &= \\lVert A \\rVert_p + \\lVert B \\rVert_p.\n", "\n", "(4)\n", "---\n", "\n", "For every :math:`A, B \\in \\mathbb{R}^{n \\times n}`,\n", "\n", ".. math::\n", "\n", " \\lVert A B \\rVert_p\n", " &= \\max_{v \\in \\mathbb{R}^n \\setminus 0}\n", " \\left\\{\\frac{\\lVert A B v \\rVert_p}{\\lVert v \\rVert_p} \\right\\}\\\\\n", " &= \\frac{\\lVert A B v^* \\rVert_p}{\\lVert v^* \\rVert_p}\\\\\n", " &\\leq \\frac{\\lVert A \\rVert_p \\lVert B v \\rVert_p}{\\lVert v \\rVert_p}\\\\\n", " &\\leq \\frac{\n", " \\lVert A \\rVert_p \\lVert B \\rVert_p \\lVert v \\rVert_p\n", " }{\n", " \\lVert v \\rVert_p\n", " }\\\\\n", " &= \\lVert A \\rVert_p \\lVert B \\rVert_p." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 5\n", "==========\n", "\n", "(1)\n", "---\n", "\n", "For every :math:`A \\in \\mathbb{R}^{n \\times n}` and :math:`v \\in \\mathbb{R}^n`,\n", "\n", ".. math::\n", "\n", " \\lVert A v \\rVert_1 &= \\sum_i \\lvert A_i v \\rvert\\\\\n", " &= \\sum_i \\lvert \\sum_j A_{ij} v_j \\rvert\\\\\n", " &\\leq \\sum_i \\sum_j \\lvert A_{ij} \\rvert \\lvert v_j \\rvert\\\\\n", " &= \\sum_j \\left( \\sum_i \\lvert A_{ij} \\rvert \\right) \\lvert v_j \\rvert\\\\\n", " &\\leq \\max_{1 \\leq j \\leq n}\n", " \\left( \\sum_i \\lvert A_{ij} \\rvert \\right) \\lVert v \\rVert_1\\\\\n", " \\frac{\\lVert A v \\rVert_1}{\\lVert v \\rVert_1}\n", " &\\leq \\max_{1 \\leq j \\leq n} \\sum_i \\lvert A_{ij} \\rvert.\n", "\n", "Hence the induced matrix norm\n", ":math:`\\lVert A \\rVert_1 \\leq \\max_j \\left\\{ \\lVert A_{.j} \\rVert_1 \\right\\}`.\n", "\n", "(2)\n", "---\n", "\n", "For every :math:`A \\in \\mathbb{R}^{n \\times n}` and :math:`v \\in \\mathbb{R}^n`,\n", "\n", ".. math::\n", "\n", " \\lVert A v \\rVert_\\infty &= \\max_{1 \\leq i \\leq n} \\lvert A_i v \\rvert\\\\\n", " &= \\max_{1 \\leq i \\leq n} \\lvert \\sum_j A_{ij} v_j \\rvert\\\\\n", " &\\leq \\max_{1 \\leq i \\leq n}\n", " \\sum_j \\lvert A_{ij} \\rvert \\lVert v \\rVert_\\infty\\\\\n", " \\frac{\\lVert A v \\rVert_\\infty}{\\lVert v \\rVert_\\infty}\n", " &\\leq \\max_{1 \\leq i \\leq n} \\sum_j \\lvert A_{ij} \\rvert.\n", "\n", "Hence the induced matrix norm\n", ":math:`\\lVert A \\rVert_\\infty \\leq\n", "\\max_i \\left\\{ \\lVert A_{i.} \\rVert_1 \\right\\}`." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ ".. _dennis1996numerical-ex-3.6:\n", "\n", "Exercise 6\n", "==========\n", "\n", "For every :math:`A \\in \\mathbb{R}^{n \\times n}` and\n", ":math:`w \\in \\mathbb{R}^n`,\n", "\n", ".. math::\n", "\n", " \\lVert A w \\rVert_2^2 &= \\sum_i \\lvert A_i w \\rvert^2\\\\\n", " &= w^\\top A^\\top A w\\\\\n", " &= \\sum_{i = 1}^n \\sum_{j = 1}^n \\alpha_i v_i^\\top A^\\top A \\alpha_j v_j\n", " & \\quad & \\text{by definition of orthonormal basis and Theorem 3.5.3}\\\\\n", " &= \\sum_{i = 1}^n \\sum_{j = 1}^n \\alpha_i v_i^\\top \\alpha_j \\lambda_j v_j\n", " & \\quad & \\text{by definition of eigenvectors and eigenvalues}\\\\\n", " &= \\sum_{i = 1}^n \\sum_{j = 1}^n \\lambda_j \\alpha_i \\alpha_j v_i^\\top v_j\\\\\n", " &= \\sum_{i = 1}^n \\lambda_i \\alpha_i^2\\\\\n", " &\\leq \\max_\\lambda \\left\\{ A^\\top A \\right\\} w^\\top w\n", " & \\quad & \\text{where}\\ w = \\sum_{i = 1}^n \\alpha_i v_i\\\\\n", " &= \\max_\\lambda \\left\\{ A^\\top A \\right\\} \\lVert w \\rVert_2^2\\\\\n", " \\frac{\\lVert A w \\rVert_2}{\\lVert w \\rVert_2}\n", " &\\leq \\left( \\max_\\lambda \\left\\{ A^\\top A \\right\\} \\right)^{1/2}.\n", "\n", "Hence the induced matrix norm :math:`\\lVert A \\rVert_2 \\leq\n", "\\left( \\max_{\\lambda} \\left\\{ A^\\top A \\right\\} \\right)^{1/2}`.\n", "\n", "By inspection, :math:`\\min_w \\frac{\\lVert A w \\rVert_2}{\\lVert w \\rVert_2} \\geq\n", "\\left( \\min_{\\lambda} \\left\\{ A^\\top A \\right\\} \\right)^{1/2}`.\n", "\n", "Thus, the maximum and minimum value can be attained via picking :math:`w` to be\n", "the eigenvector corresponding to :math:`\\lambda_\\max` and\n", ":math:`\\lambda_\\min` respectively." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 7\n", "==========\n", "\n", "Recall that\n", "\n", ".. math::\n", "\n", " \\text{tr}(A^\\top B) = \\sum_i A^\\top_i B_{.i} = a^\\top b\n", "\n", "where :math:`a =\n", "\\begin{bmatrix} {A_{.0}}^\\top & \\cdots & {A_{.n}}^\\top \\end{bmatrix}^\\top` and\n", ":math:`b =\n", "\\begin{bmatrix} {B_{.0}}^\\top & \\cdots & {B_{.n}}^\\top \\end{bmatrix}^\\top`.\n", "\n", "The Frobenius norm of :math:`A` can be derived as\n", "\n", ".. math::\n", "\n", " \\text{tr}\\left( A^\\top A \\right) &= \\sum_i A^\\top_i A_{.i}\\\\\n", " &= \\sum_i \\sum_j A^\\top_{ij} A_{ji}\\\\\n", " &= \\sum_j \\sum_i \\lvert A_{ji} \\rvert^2\\\\\n", " \\text{tr}\\left( A^\\top A \\right)^{1 / 2}\n", " &= \\left( \\sum_j \\sum_i \\lvert A_{ij} \\rvert^2 \\right)^{1 / 2}\\\\\n", " &= \\lVert A \\rVert_F." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ ".. _dennis1996numerical-ex-3.8:\n", "\n", "Exercise 8\n", "==========\n", "\n", "The following are using matrix norms.\n", "\n", "(1)\n", "---\n", "\n", "For every :math:`A \\in \\mathbb{R}^{n \\times n}`,\n", "\n", ".. math::\n", "\n", " \\lVert A \\rVert_{p}\n", " &= \\left( \\sum_j \\sum_i \\lvert A_{ij} \\rvert^p \\right)^{1 / p}\\\\\n", " &= \\left( \\sum_j \\sum_i \\lvert A_{ji} \\rvert^p \\right)^{1 / p}\\\\\n", " &= \\left( \\sum_j \\sum_i \\lvert A^\\top_{ij} \\rvert^p \\right)^{1 / p}\\\\\n", " &= \\lVert A^\\top \\rVert_{p}.\n", "\n", "(2)\n", "---\n", "\n", "For every :math:`A, B \\in \\mathbb{R}^{n \\times n}`,\n", "\n", ".. math::\n", "\n", " \\lVert AB \\rVert_F^2 &= \\sum_i \\sum_j \\lvert (AB)_{ij} \\rvert^2\\\\\n", " &= \\sum_i \\sum_j \\lvert A_{i.} B_{.j} \\rvert^2\\\\\n", " &= \\sum_i \\sum_j \\left\\vert \\sum_k A_{ik} B_{kj} \\right\\vert^2\\\\\n", " &\\leq \\sum_i \\sum_j \\sum_k \\lvert A_{ik} \\rvert^2 \\lvert B_{kj} \\rvert^2\\\\\n", " \\lVert AB \\rVert_F\n", " &\\leq \\min\\left\\{\n", " \\lVert A \\rVert_F \\lVert B \\rVert_2,\n", " \\lVert A \\rVert_2 \\lVert B \\rVert_F\n", " \\right\\}.\n", "\n", "(3)\n", "---\n", "\n", "For every :math:`v \\in \\mathbb{R}^n` and :math:`A \\in \\mathbb{R}^{n \\times n}`,\n", "\n", ".. math::\n", "\n", " \\lVert A v \\rVert_2^2 &= \\sum_i \\lvert A_i v \\rvert^2\\\\\n", " &= \\sum_i \\left\\vert \\sum_j A_{ij} v_j \\right\\vert^2\\\\\n", " &\\leq \\sum_i \\sum_j \\lvert A_{ij} \\rvert^2 \\lvert v_j \\rvert^2\\\\\n", " &\\leq \\sum_i \\sum_j \\lvert A_{ij} \\rvert^2 \\cdot\n", " \\left( \\sum_k \\lvert v_k \\rvert^2 \\right)\\\\\n", " \\lVert A v \\rVert_2 &\\leq \\lVert A \\rVert_F \\lVert v \\rVert_2.\n", "\n", "(4)\n", "---\n", "\n", "For every :math:`A, B \\in \\mathbb{R}^{n \\times n}`,\n", "\n", ".. math::\n", "\n", " \\lVert A + B \\rVert_F^2 &= \\sum_k \\lvert A_k + B_k \\rvert^2\\\\\n", " &= \\lVert A \\rVert_F^2 + \\lVert B \\rVert_F^2 +\n", " \\sum_k \\lvert 2 A_k B_k \\rvert\\\\\n", " &\\leq \\lVert A \\rVert_F^2 + \\lVert B \\rVert_F^2 +\n", " 2 \\lVert A \\rVert_F \\lVert B \\rVert_F\n", " & \\quad & \\text{Cauchy–Schwarz inequality (3.1.10)}\\\\\n", " \\lVert A + B \\rVert_F &\\leq \\lVert A \\rVert_F + \\lVert B \\rVert_F." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 9\n", "==========\n", "\n", "Notice that :math:`\\lVert v w^\\top \\rVert_F = \\lVert v w^\\top \\rVert_2` is true\n", "by definition of matrix norm.\n", "\n", "For any :math:`v, w \\in \\mathbb{R}^n`,\n", "\n", ".. math::\n", "\n", " \\lVert v w^\\top \\rVert_2^2\n", " &= \\sum_i \\sum_j \\left\\vert \\left( v w^\\top \\right)_{ij} \\right\\vert^2\\\\\n", " &= \\sum_i \\sum_j \\lvert v_i \\rvert^2 \\lvert w_j \\rvert^2\\\\\n", " &= \\sum_i \\lvert v_i \\rvert^2 \\sum_j \\lvert w_j \\rvert^2\\\\\n", " \\lVert v w^\\top \\rVert_2 &= \\lVert v \\rVert_2 \\lVert w \\rVert_2." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 10\n", "===========\n", "\n", "(1)\n", "---\n", "\n", "This solution was adapted from :cite:`solverfn`.\n", "\n", "Let :math:`\\lVert \\cdot \\rVert` be any norm on :math:`\\mathbb{R}^{n \\times n}`\n", "that obeys (3.1.10). Define :math:`E, I \\in \\mathbb{R}^{n \\times n}` such that\n", ":math:`\\lVert E \\rVert < 1` and :math:`\\lVert I \\rVert = 1`.\n", "\n", "Given :math:`m, n \\in \\mathbb{Z}^*` and :math:`m > n`,\n", ":math:`S_k = I + E + E^2 + \\cdots + E^k` for :math:`k \\in \\mathbb{Z}^*` is a\n", "Cauchy sequence because\n", "\n", ".. math::\n", "\n", " \\lVert S_m - S_n \\rVert &= \\lVert E^{n + 1} + \\cdots + E^m \\rVert\\\\\n", " &\\leq \\lVert E \\rVert^{n + 1} + \\cdots + \\lVert E \\rVert^m\\\\\n", " &= \\epsilon\n", "\n", "where :math:`\\epsilon \\in \\mathbb{R}^+`.\n", "\n", "Since :math:`E` is a convergent matrix, notice that\n", "\n", ".. math::\n", "\n", " (I - E) S_k &= (I - E) (I + E + E^2 + \\cdots + E^k)\\\\\n", " &= I - E + E - E^2 + E^2 - E^3 + \\cdots - E^{k + 1}\\\\\n", " &= I - E^{k + 1}\n", "\n", "and\n", "\n", ".. math::\n", "\n", " (I - E) \\lim_{k \\rightarrow \\infty} S_k\n", " &= \\lim_{k \\rightarrow \\infty} I - E^{k + 1}\\\\\n", " &= I - \\lim_{k \\rightarrow \\infty} E^{k + 1}\\\\\n", " &= I\n", "\n", "implies :math:`(I - E)^{-1} = \\lim_{k \\rightarrow \\infty} S_k`.\n", "\n", "Combining the foregoing observations with the fact that\n", ":math:`\\lVert S_k \\rVert` is a convergent geometric series yield\n", "\n", ".. math::\n", "\n", " \\lVert (I - E)^{-1} \\rVert &= \\lim_{k \\rightarrow \\infty} \\lVert S_k \\rVert\\\\\n", " &\\leq \\sum_{i = 0}^\\infty \\lVert E \\rVert^i\\\\\n", " &= \\frac{1}{1 - \\lVert E \\rVert}.\n", "\n", "(2)\n", "---\n", "\n", "Let :math:`A, B \\in \\mathbb{R}^{n \\times n}`, :math:`A` is nonsingular, and\n", ":math:`\\lVert A^{-1} (B - A) \\rVert < 1`. Observe that\n", "\n", ".. math::\n", "\n", " \\lVert A^{-1} (B - A) \\rVert &= \\lVert A^{-1} B - I \\rVert\\\\\n", " &= \\lVert -(I - A^{-1} B) \\rVert\\\\\n", " &= \\lvert -1 \\rvert \\lVert I - A^{-1} B \\rVert\\\\\n", " &< 1.\n", "\n", "Hence :math:`\\left( I - \\left( I - A^{-1} B \\right) \\right)^{-1} = B^{-1} A`\n", "exists. Thus\n", "\n", ".. math::\n", "\n", " \\left\\Vert \\left( I - A^{-1} (B - A) \\right)^{-1} \\right\\Vert\n", " &\\leq \\frac{1}{1 - \\lVert A^{-1} (B - A) \\rVert}\n", " & \\quad & \\text{(3.1.19)}\\\\\n", " \\lVert B^{-1} A \\rVert &\\leq \\frac{1}{1 - \\lVert A^{-1} (B - A) \\rVert}\\\\\n", " \\lVert B^{-1} \\rVert \\lVert A \\rVert\n", " &\\leq \\frac{1}{1 - \\lVert A^{-1} (B - A) \\rVert}\n", " & \\quad & \\text{(3.1.10)}\\\\\n", " \\lVert B^{-1} \\rVert\n", " &\\leq \\frac{1}{1 - \\lVert A^{-1} (B - A) \\rVert}\n", " \\frac{1}{\n", " \\max_{v \\in \\mathbb{R}^n, v \\neq 0}\n", " \\frac{\\lVert A v \\rVert}{\\lVert v \\rVert}\n", " }\\\\\n", " &\\leq \\frac{1}{1 - \\lVert A^{-1} (B - A) \\rVert}\n", " \\frac{1}{\n", " \\min_{v \\in \\mathbb{R}^n, v \\neq 0}\n", " \\frac{\\lVert A v \\rVert}{\\lVert v \\rVert}\n", " }\\\\\n", " \\lVert B^{-1} \\rVert\n", " &\\leq \\frac{\\lVert A^{-1} \\rVert}{1 - \\lVert A^{-1} (B - A) \\rVert}\n", " & \\quad & \\text{(3.1.18).}" ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ ".. _dennis1996numerical-ex-3.11:\n", "\n", "Exercise 11\n", "===========\n", "\n", "For any :math:`A \\in \\mathbb{R}^{n \\times n}` and any orthonormal set\n", ":math:`\\left\\{ v_i \\in \\mathbb{R}^n \\right\\}_{i = 1}^n`,\n", "\n", ".. math::\n", "\n", " \\sum_{i = 1}^n \\lVert A v_i \\rVert_2^2\n", " &= \\lVert A \\begin{bmatrix} v_0 & \\cdots & v_n \\end{bmatrix} \\rVert_F^2\\\\\n", " &= \\lVert A Q \\rVert_F^2\\\\\n", " &= \\text{tr}\\left( Q^\\top A^\\top A Q \\right)\\\\\n", " &= \\text{tr}\\left( A^\\top A Q Q^\\top \\right)\\\\\n", " &= \\text{tr}\\left( A^\\top A \\right)\\\\\n", " &= \\lVert A \\rVert_F^2." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 12\n", "===========\n", "\n", "(1)\n", "---\n", "\n", "For any :math:`v \\in \\mathbb{R}^n` and orthogonal matrix\n", ":math:`Q \\in \\mathbb{R}^{n \\times n}`,\n", "\n", ".. math::\n", "\n", " \\lVert Q v \\rVert_2^2 &= v^\\top Q^\\top Q v\\\\\n", " &= v^\\top v\\\\\n", " \\lVert Q v \\rVert_2 &= \\lVert v \\rVert_2.\n", "\n", "(2)\n", "---\n", "\n", "By definition of element-wise matrix norm (3.1.6) and the result of\n", ":ref:`Exercise 3.11 `, the following holds for\n", "every :math:`A \\in \\mathbb{R}^{n \\times n}`:\n", "\n", ".. math::\n", "\n", " \\begin{aligned}\n", " \\lVert Q A \\rVert_2 &= \\lVert Q A \\rVert_F\\\\\n", " &= \\lVert A \\rVert_F\\\\\n", " &= \\lVert A \\rVert_2\n", " \\end{aligned}\n", " \\quad \\text{and} \\quad\n", " \\begin{aligned}\n", " \\lVert A Q \\rVert_2 &= \\lVert A Q \\rVert_F\\\\\n", " &= \\lVert A \\rVert_F\\\\\n", " &= \\lVert A \\rVert_2.\n", " \\end{aligned}\n", "\n", "(3)\n", "---\n", "\n", "For any :math:`v \\in \\mathbb{R}^n` and orthogonal matrix\n", ":math:`Q \\in \\mathbb{R}^{n \\times n}`, the induced matrix norm gives\n", "\n", ".. math::\n", "\n", " \\lVert Q \\rVert_2\n", " &= \\max_{v \\in \\mathbb{R}^n, v \\neq 0}\n", " \\frac{\\lVert Q v \\rVert_2}{\\lVert v \\rVert_2}\\\\\n", " &= \\frac{\\left( v_*^\\top Q^\\top Q v_* \\right)^{1/2}}{\\lVert v_* \\rVert_2}\\\\\n", " &= \\frac{\\left( v_*^\\top v_* \\right)^{1/2}}{\\lVert v_* \\rVert_2}\\\\\n", " &= \\frac{\\lVert v_* \\rVert_2}{\\lVert v_* \\rVert_2}\\\\\n", " &= 1\n", "\n", "while the element-wise matrix norm yields\n", "\n", ".. math::\n", "\n", " \\lVert Q \\rVert_2 &= \\lVert Q \\rVert_F\\\\\n", " &= \\text{tr}(Q^\\top Q)^{1/2}\\\\\n", " &= \\text{tr}(I)^{1/2}\\\\\n", " &= \\sqrt{n}." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 13\n", "===========\n", "\n", ".. math::\n", "\n", " 2 x_2 = 6 \\Rightarrow x_2 = 3.\n", "\n", ".. math::\n", "\n", " \\begin{bmatrix}\n", " 1 & 2\\\\\n", " 3 & 1\n", " \\end{bmatrix} \\begin{bmatrix} x_0\\\\ x_1 \\end{bmatrix} =\n", " \\begin{bmatrix} 5\\\\ 5 \\end{bmatrix}\n", " \\begin{matrix} \\\\ R_2 - 3 R_1 \\end{matrix}\n", " \\Rightarrow\n", " \\begin{bmatrix}\n", " 1 & 2\\\\\n", " 0 & -5\n", " \\end{bmatrix} \\begin{bmatrix} x_0\\\\ x_1 \\end{bmatrix} =\n", " \\begin{bmatrix} 5\\\\ -10 \\end{bmatrix}\n", " \\Rightarrow\n", " \\begin{aligned}\n", " -5 x_1 &= -10\\\\\n", " x_1 &= 2\n", " \\end{aligned}\n", " \\Rightarrow\n", " \\begin{aligned}\n", " x_0 + 2 x_1 &= 5\\\\\n", " x_0 &= 1\n", " \\end{aligned}\n", "\n", ".. math::\n", "\n", " \\begin{bmatrix}\n", " -1 & 4\\\\\n", " 2 & -1\n", " \\end{bmatrix} \\begin{bmatrix} x_3\\\\ x_4 \\end{bmatrix} =\n", " \\begin{bmatrix} 16\\\\ 3 \\end{bmatrix}\n", " \\begin{matrix} \\\\ R_2 - (-2) R_1 \\end{matrix}\n", " \\Rightarrow\n", " \\begin{bmatrix}\n", " -1 & 4\\\\\n", " 0 & 7\n", " \\end{bmatrix} \\begin{bmatrix} x_3\\\\ x_4 \\end{bmatrix} =\n", " \\begin{bmatrix} 16\\\\ 35 \\end{bmatrix}\n", " \\Rightarrow\n", " \\begin{aligned}\n", " 7 x_4 &= 35\\\\\n", " x_4 &= 5\n", " \\end{aligned}\n", " \\Rightarrow\n", " \\begin{aligned}\n", " 2 x_3 - x_4 &= 3\\\\\n", " x_3 &= 4\n", " \\end{aligned}\n", "\n", "Hence :math:`x = \\begin{bmatrix} 1 & 2 & 3 & 4 & 5 \\end{bmatrix}^\\top`." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ ".. _dennis1996numerical-ex-3.14:\n", "\n", "Exercise 14\n", "===========\n", "\n", "Let :math:`u \\in \\mathbb{R}^n`. Observe that :math:`I - u u^\\top` is a\n", "symmetric matrix. In order for that matrix to be orthogonal, the following must\n", "hold\n", "\n", ".. math::\n", "\n", " (I - u u^\\top)^\\top (I - u u^\\top) &= I - 2 u u^\\top + u u^\\top u u^\\top\\\\\n", " &= I - 2 u u^\\top + u (u^\\top u) u^\\top\\\\\n", " &= I.\n", "\n", "The cancellation of terms will occur when :math:`u^\\top u = 2`.\n", "\n", "Let :math:`v \\in \\mathbb{R}^n \\setminus \\boldsymbol{0}` and choose\n", ":math:`u = \\frac{\\sqrt{2} v}{\\lVert v \\rVert_2}` so that\n", "\n", ".. math::\n", "\n", " u^\\top u &= \\left( \\frac{\\sqrt{2}}{\\lVert v \\rVert_2} v^\\top \\right)\n", " \\left( \\frac{\\sqrt{2}}{\\lVert v \\rVert_2} v \\right)\\\\\n", " &= 2 \\frac{v^\\top v}{\\lVert v \\rVert_2^2}\\\\\n", " &= 2.\n", "\n", "Thus :math:`Q = I - 2 \\frac{u u^\\top}{u^\\top u}` will be orthogonal for any\n", ":math:`u \\in \\mathbb{R}^n`.\n", "\n", "The final property that the Householder transformation must satisfy is\n", "\n", ".. math::\n", "\n", " \\left( I - 2 \\frac{u u^\\top}{u^\\top u} \\right) v\n", " &= v - 2 \\frac{u u^\\top}{u^\\top u} v\\\\\n", " &= v - u \\frac{2 u^\\top v}{u^\\top u}\\\\\n", " &= \\alpha e_1\n", "\n", "where :math:`e_1 =\n", "\\begin{bmatrix} 1 & 0 & \\cdots & 0 \\end{bmatrix}^\\top \\in \\mathbb{R}^n` and\n", ":math:`\\alpha \\in \\mathbb{R} \\setminus 0`.\n", "\n", "Notice that for :math:`2 \\leq i \\leq n`,\n", ":math:`v_i - u_i \\frac{2 u^\\top v}{u^\\top u} = 0` constrains\n", ":math:`\\frac{2 u^\\top v}{u^\\top u} = 1` with :math:`u_i = v_i` such that\n", "\n", ".. math::\n", "\n", " 2 u^\\top v &= u^\\top u\\\\\n", " 2 \\left( u_1 v_1 + v_2^2 + \\ldots + v_n^2 \\right)\n", " &= \\left( u_1^2 + v_2^2 + \\ldots + v_n^2 \\right)\\\\\n", " 0 &= u_1^2 - 2 u_1 v_1 - v_2^2 - v_3^2 - \\ldots - v_n^2\\\\\n", " u_1 &= \\frac{-(-2 v_1) \\pm\n", " \\sqrt{(-2 v_1)^2 - 4 (1) (- v_2^2 - v_3^2 - \\ldots - v_n^2)}}{2 (1)}\\\\\n", " &= \\frac{2 v_1 \\pm \\sqrt{4 (v_1^2 + v_2^2 + \\ldots + v_n^2)}}{2}\\\\\n", " &= v_1 \\pm \\sqrt{v_1^2 + v_2^2 + \\ldots + v_n^2}.\n", "\n", "In the more general case of\n", "\n", ".. math::\n", "\n", " Q v = \\begin{bmatrix}\n", " \\alpha_1\\\\\n", " \\alpha_2\\\\\n", " \\vdots\\\\\n", " \\alpha_k\\\\\n", " 0\\\\\n", " 0\\\\\n", " \\vdots\\\\\n", " \\end{bmatrix}\n", "\n", "where :math:`\\alpha_1, \\alpha_2, \\ldots, \\alpha_k \\in \\mathbb{R} \\setminus 0`,\n", "the same analysis\n", "\n", ".. math::\n", "\n", " 2 u^\\top v &= u^\\top u\\\\\n", " 2 \\left( u_1 v_1 + v_2^2 + \\ldots + v_n^2 \\right)\n", " &= \\left( u_1^2 + v_2^2 + \\ldots + v_n^2 \\right)\\\\\n", " 0 &= u_1^2 - 2 u_1 v_1 + u_2^2 - 2 u_2 v_2 + \\ldots + u_k^2 - 2 u_k v_k -\n", " v_{k + 1}^2 - v_{k + 2}^2 - \\ldots - v_n^2\n", "\n", "still holds. Hence the solution must satisfy\n", "\n", ".. math::\n", "\n", " u_1^2 - 2 u_1 v_1 + u_2^2 - 2 u_2 v_2 + \\ldots + u_k^2 - 2 u_k v_k =\n", " \\lVert v \\rVert^2 - v_1^2 - v_2^2 - \\ldots - v_k^2.\n", "\n", "One solution would be to set :math:`u_i = 0` for :math:`i < k` and\n", ":math:`u_j = v_j` for :math:`k < j` resulting in\n", "\n", ".. math::\n", "\n", " u_k^2 - 2 u_k v_k &= \\lVert v \\rVert^2 - v_1^2 - v_2^2 - \\ldots - v_k^2\\\\\n", " u_k^2 - 2 u_k v_k + v_1^2 + v_2^2 + \\ldots + v_k^2 - \\lVert v \\rVert^2 &= 0\\\\\n", " u_k &= \\frac{\n", " -(-2 v_k) \\pm\n", " \\sqrt{\n", " (-2 v_k)^2 -\n", " 4 (1) (v_1^2 + v_2^2 + \\ldots + v_k^2 - \\lVert v \\rVert^2)\n", " }\n", " }{2 (1)}\\\\\n", " u_k &= \\frac{\n", " 2 v_k \\pm\n", " \\sqrt{\n", " 4 \\lVert v \\rVert^2 -\n", " 4 \\left( v_1^2 + v_2^2 + \\ldots + v_{k - 1}^2 \\right)\n", " }\n", " }{2}\\\\\n", " u_k &= \\frac{\n", " 2 v_k \\pm\n", " \\sqrt{\n", " 4 \\left( v_k^2 + v_{k + 1}^2 + \\ldots + v_n^2 \\right)\n", " }\n", " }{2}\\\\\n", " u_k &= v_k \\pm \\sqrt{v_k^2 + v_{k + 1}^2 + \\ldots + v_n^2}." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 15\n", "===========" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy\n", "\n", "def QR_Householder(A):\n", " m, n = A.shape\n", " Q = numpy.eye(m)\n", " R = A\n", " for i in range(min(m - 1, n)):\n", " x = R[i:, i]\n", " #use operations of the same sign to avoid cancellation\n", " if x[0] < 0:\n", " u = numpy.hstack((x[0] - numpy.linalg.norm(x), x[1:])).T\n", " else:\n", " u = numpy.hstack((x[0] + numpy.linalg.norm(x), x[1:])).T\n", " Q_i = numpy.eye(m)\n", " I = numpy.eye(len(u))\n", " Q_i[i:,i:] = I - 2 * numpy.outer(u, u.T) / numpy.dot(u.T, u)\n", " Q = numpy.dot(Q, Q_i)\n", " R = numpy.dot(Q.T, A)\n", " return Q, R\n", "\n", "A = numpy.asfarray([[1, -1],\n", " [1, 0],\n", " [1, 1],\n", " [1, 2]])\n", "Q, R = QR_Householder(A)\n", "print(Q)\n", "print(R)\n", "print(numpy.round(numpy.dot(Q, R)))" ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 16\n", "===========\n", "\n", "Given :math:`x = \\begin{bmatrix} -1 & 1 \\end{bmatrix}^\\top`,\n", "\n", ".. math::\n", "\n", " x^\\top \\begin{bmatrix} 1 & 2 \\\\ 2 & 3 \\end{bmatrix} x\n", " &= x^\\top \\begin{bmatrix} -1 + 2 \\\\ -2 + 3 \\end{bmatrix}\\\\\n", " &= 1 - 2 - 2 + 3\\\\\n", " &= 0\n", "\n", "illustrates that the positive definite criterion is not satisfied. Running\n", "Cholesky decomposition gives some negative diagonal matrix values." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy\n", "\n", "def Cholesky_LDL(A):\n", " n = A.shape[0]\n", " D = numpy.zeros(n)\n", " L = numpy.eye(n)\n", "\n", " for i in range(n):\n", " #compute lower triangular\n", " for j in range(i):\n", " _ = A[i,j] - sum([L[i,k] * L[j,k] * D[k] for k in range(j)])\n", " L[i,j] = _ / D[j]\n", "\n", " #compute diagonal\n", " D[i] = A[i,i] - sum([L[i,k]**2 * D[k] for k in range(i)])\n", " return L, D\n", "\n", "A = numpy.asfarray([[1, 2],\n", " [2, 3]\n", " ])\n", "L, D = Cholesky_LDL(A)\n", "print(L, D)\n", "print(numpy.dot(L, numpy.dot(numpy.diag(D), L.T)))" ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 17\n", "===========\n", "\n", "The matrix will no longer be positive definite when :math:`10` is replaced by\n", ":math:`7`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "A = numpy.asfarray([[4, 6, -2],\n", " [6, 10, 1],\n", " [-2, 1, 22]\n", " ])\n", "L, D = Cholesky_LDL(A)\n", "print(L, D)\n", "print(numpy.dot(L, numpy.dot(numpy.diag(D), L.T)))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "A = numpy.asfarray([[4, 6, -2],\n", " [6, 7, 1],\n", " [-2, 1, 22]\n", " ])\n", "L, D = Cholesky_LDL(A)\n", "print(L, D)\n", "print(numpy.dot(L, numpy.dot(numpy.diag(D), L.T)))" ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 18\n", "===========\n", "\n", "(1)\n", "---\n", "\n", "Since the rows and columns of an orthogonal matrix :math:`Q` are orthonormal\n", "vectors, :math:`\\left\\{ \\lVert Q_{i.} \\rVert_1, \\lVert Q_{.j} \\rVert_1 \\right\\}\n", "\\leq \\sqrt{n}`.\n", "\n", "Given :math:`A = QR`,\n", ":math:`\\frac{1}{n} \\kappa_1(A) \\leq \\kappa_1(R) \\leq n \\kappa_1(A)` because\n", "\n", ".. math::\n", "\n", " \\kappa_1(A) &= \\lVert A \\rVert_1 \\lVert A^{-1} \\rVert_1\\\\\n", " &= \\lVert Q R \\rVert_1 \\lVert R^{-1} Q^\\top \\rVert_1\\\\\n", " &\\leq \\lVert Q \\rVert_1 \\lVert R \\rVert_1\n", " \\lVert R^{-1} \\rVert_1 \\lVert Q^\\top \\rVert_1\n", " & \\quad & \\text{(3.1.10)}\\\\\n", " &\\leq n \\lVert R \\rVert_1 \\lVert R^{-1} \\rVert_1\n", " & \\quad & \\text{(3.1.8a)}\\\\\n", " \\frac{1}{n} \\kappa_1(A) &\\leq \\kappa_1(R)\n", "\n", "and\n", "\n", ".. math::\n", "\n", " \\kappa_1(A)\n", " &= \\left(\n", " \\max_{v \\in \\mathbb{R}^n, v \\neq 0}\n", " \\frac{\\lVert Q R v \\rVert_1}{\\lVert v \\rVert_1}\n", " \\right)\n", " \\left(\n", " \\min_{v \\in \\mathbb{R}^n, v \\neq 0}\n", " \\frac{\\lVert Q R v \\rVert_1}{\\lVert v \\rVert_1}\n", " \\right)^{-1}\n", " & \\quad & \\text{(3.1.18)}\\\\\n", " &= \\left(\n", " \\max_{v \\in \\mathbb{R}^n, v \\neq 0}\n", " \\frac{\n", " \\sum_i \\left\\vert Q_{i \\cdot} \\cdot R v \\right\\vert\n", " }{\\lVert v \\rVert_1}\n", " \\right)\n", " \\left(\n", " \\min_{v \\in \\mathbb{R}^n, v \\neq 0}\n", " \\frac{\n", " \\sum_i \\left\\vert Q_{i \\cdot} \\cdot R v \\right\\vert\n", " }{\\lVert v \\rVert_1}\n", " \\right)^{-1}\\\\\n", " &= \\left(\n", " \\max_{v \\in \\mathbb{R}^n, v \\neq 0}\n", " \\frac{\n", " \\sum_i\n", " \\left\\Vert Q_{i \\cdot} \\right\\Vert_2\n", " \\left\\Vert R v \\right\\Vert_2\n", " \\left\\vert \\cos \\theta_i \\right\\vert\n", " }{\\lVert v \\rVert_1}\n", " \\right)\n", " \\left(\n", " \\min_{v \\in \\mathbb{R}^n, v \\neq 0}\n", " \\frac{\n", " \\sum_i\n", " \\left\\Vert Q_{i \\cdot} \\right\\Vert_2\n", " \\left\\Vert R v \\right\\Vert_2\n", " \\left\\vert \\cos \\theta_i \\right\\vert\n", " }{\\lVert v \\rVert_1}\n", " \\right)^{-1}\n", " & \\quad & \\text{(3.1.21)}\\\\\n", " &= \\left(\n", " \\max_{v \\in \\mathbb{R}^n, v \\neq 0}\n", " \\frac{\\left\\Vert R v \\right\\Vert_2}{\\lVert v \\rVert_1}\n", " \\right)\n", " \\left(\n", " \\min_{v \\in \\mathbb{R}^n, v \\neq 0}\n", " \\frac{\\left\\Vert R v \\right\\Vert_2}{\\lVert v \\rVert_1}\n", " \\right)^{-1}\\\\\n", " &\\geq\n", " \\left(\n", " \\max_{v \\in \\mathbb{R}^n, v \\neq 0}\n", " \\frac{\\left\\Vert R v \\right\\Vert_1}{\\sqrt{n} \\lVert v \\rVert_1}\n", " \\right)\n", " \\left(\n", " \\min_{v \\in \\mathbb{R}^n, v \\neq 0}\n", " \\frac{\\sqrt{n} \\left\\Vert R v \\right\\Vert_1}{\\lVert v \\rVert_1}\n", " \\right)^{-1}\n", " & \\quad & \\text{(3.1.5)}\\\\\n", " &= \\frac{1}{n} \\lVert R \\rVert_1 \\lVert R^{-1} \\rVert_1\\\\\n", " n \\kappa_1(A) &\\geq \\kappa_1(R).\n", "\n", "(2)\n", "---\n", "\n", ".. math::\n", "\n", " \\kappa_2(A) &= \\lVert A \\rVert_2 \\lVert A^{-1} \\rVert_2\\\\\n", " &= \\left(\n", " \\max_{v \\in \\mathbb{R}^n, v \\neq 0}\n", " \\frac{\\lVert Q R v \\rVert_2}{\\lVert v \\rVert_2}\n", " \\right)\n", " \\left(\n", " \\min_{v \\in \\mathbb{R}^n, v \\neq 0}\n", " \\frac{\\lVert Q R v \\rVert_2}{\\lVert v \\rVert_2}\n", " \\right)^{-1}\n", " & \\quad & \\text{(3.1.18)}\\\\\n", " &= \\left(\n", " \\max_{v \\in \\mathbb{R}^n, v \\neq 0}\n", " \\frac{\\lVert R v \\rVert_2}{\\lVert v \\rVert_2}\n", " \\right)\n", " \\left(\n", " \\min_{v \\in \\mathbb{R}^n, v \\neq 0}\n", " \\frac{\\lVert R v \\rVert_2}{\\lVert v \\rVert_2}\n", " \\right)^{-1}\\\\\n", " &= \\lVert R \\rVert_2 \\lVert R^{-1} \\rVert_2\\\\\n", " &= \\kappa_2(R)." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 19\n", "===========\n", "\n", "Given\n", "\n", ".. math::\n", "\n", " R = \\begin{bmatrix} 1 & 1\\\\ 0 & 10^{-7} \\end{bmatrix}\n", " \\quad \\text{and} \\quad\n", " R^{-1} = \\begin{bmatrix} 1 & -10^{7}\\\\ 0 & 10^7 \\end{bmatrix},\n", "\n", ".. math::\n", "\n", " \\kappa_1(R) &= \\left\\Vert R \\right\\Vert_1 \\left\\Vert R^{-1} \\right\\Vert_1\\\\\n", " &= \\left(\n", " \\max_j \\left\\Vert R_{\\cdot j} \\right\\Vert_1\n", " \\right)\n", " \\left(\n", " \\max_j \\left\\Vert R^{-1}_{\\cdot j} \\right\\Vert_1\n", " \\right)\\\\\n", " &= \\left\\Vert R_{\\cdot 2} \\right\\Vert_1\n", " \\left\\Vert R^{-1}_{\\cdot 2} \\right\\Vert_1\\\\\n", " &= 20000002." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 20\n", "===========\n", "\n", "See :cite:`cline1979estimate` for an elegant and concise exposition on\n", "estimating the upper triangular condition number." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy\n", "\n", "def R_solve(R, b):\n", " x = numpy.copy(b)\n", " n = b.shape[0]\n", " d = numpy.diag(R)\n", "\n", " x[-1] = x[-1] / d[-1]\n", " for i in range(n - 2, -1, -1):\n", " _ = x[i] - sum(x[i+1:n] * R[i,i+1:n])\n", " x[i] = _ / d[i]\n", " return x" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def cond_est(R):\n", " #estimate ||R||_1\n", " R_1 = max(numpy.linalg.norm(R, ord=1, axis=0))\n", "\n", " #solve R^T z = e\n", " d = numpy.diag(R)\n", " n = d.shape[0]\n", " z = numpy.zeros(d.shape)\n", "\n", " z[0] = 1 / d[0]\n", " p = R[0,:] * z[0]\n", " m = numpy.zeros(p.shape)\n", " for j in range(1, n):\n", " zp = (1 - p[j]) / d[j]\n", " zm = (-1 - p[j]) / d[j]\n", " tmp_p = abs(zp)\n", " tmp_m = abs(zm)\n", "\n", " m[j+1:n] = p[j+1:n] + R[j,j+1:n] * zm\n", " tmp_m += sum(numpy.abs(m[j+1:n]) / numpy.abs(d[j+1:n]))\n", "\n", " p[j+1:n] = p[j+1:n] + R[j,j+1:n] * zp\n", " tmp_p += sum(numpy.abs(p[j+1:n]) / numpy.abs(d[j+1:n]))\n", "\n", " if tmp_p >= tmp_m:\n", " z[j] = zp\n", " else:\n", " z[j] = zm\n", " p[j+1:n] = m[j+1:n]\n", " z_1 = numpy.linalg.norm(z)\n", "\n", " #solve R y = z\n", " y_1 = numpy.linalg.norm(R_solve(R, z), ord=1)\n", "\n", " return R_1 * y_1 / z_1" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "R = numpy.asfarray([[1, 1],\n", " [0, 1e-7]\n", " ])\n", "cond_est(R)" ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 21\n", "===========\n", "\n", "See Lemma 3.4.2 and :cite:`golub2012matrix` for an algorithm to robustly\n", "compute the Jacobi matrix.\n", "\n", "Given :math:`\\hat{R}(\\alpha, \\beta) \\in \\mathbb{R}^{n \\times n}` where\n", ":math:`n = 2`, the Schur decomposition states that :math:`\\hat{R} = Q U Q^*`\n", "where :math:`Q` is a unitary matrix. Let\n", ":math:`A_{ij} = \\mathbb{\\hat{R}}(\\alpha, \\beta)_{ij}` where\n", ":math:`1 \\leq i,j \\leq 2`.\n", "\n", ".. math::\n", "\n", " \\det(A) &= A_{11} A_{22} - A_{12} A_{21}\\\\\n", " &= \\frac{\\alpha^2}{\\alpha^2 + \\beta^2} -\n", " \\frac{-\\beta^2}{\\alpha^2 + \\beta^2}\\\\\n", " &= \\frac{\\alpha^2 + \\beta^2}{\\alpha^2 + \\beta^2}\\\\\n", " &= 1\n", "\n", "and\n", "\n", ".. math::\n", "\n", " \\text{tr}(A) &= A_{11} + A_{22}\\\\\n", " &= \\frac{2 \\alpha}{\\sqrt{\\alpha^2 + \\beta^2}}.\n", "\n", "The eigenvalues are\n", "\n", ".. math::\n", "\n", " \\lambda_1 &= \\frac{\\text{Tr}(A) + \\sqrt{\\text{Tr}(A)^2 - 4 \\det(A)}}{2}\\\\\n", " &= \\frac{\\alpha}{\\sqrt{\\alpha^2 + \\beta^2}} +\n", " \\sqrt{\\frac{\\alpha^2}{\\alpha^2 + \\beta^2} - 1}\\\\\n", " &= \\frac{\\alpha}{\\sqrt{\\alpha^2 + \\beta^2}} +\n", " \\sqrt{\\frac{-\\beta^2}{\\alpha^2 + \\beta^2}}\\\\\n", " &= \\frac{\\alpha + i \\beta}{\\sqrt{\\alpha^2 + \\beta^2}}\n", "\n", "and\n", "\n", ".. math::\n", "\n", " \\lambda_2 &= \\frac{\\text{Tr}(A) - \\sqrt{\\text{Tr}(A)^2 - 4 \\det(A)}}{2}\\\\\n", " &= \\frac{\\alpha - i \\beta}{\\sqrt{\\alpha^2 + \\beta^2}}.\n", "\n", "Solving for the eigenvectors yields\n", "\n", ".. math::\n", "\n", " A v_1 &= \\lambda_1 v_1\\\\\n", " \\left( A - \\lambda_1 I \\right) v_1 &= 0\\\\\n", " \\frac{1}{\\sqrt{\\alpha^2 + \\beta^2}}\n", " \\begin{bmatrix} - i \\beta & -\\beta\\\\ \\beta & - i \\beta \\end{bmatrix} v_1\n", " &= 0\n", "\n", "and\n", "\n", ".. math::\n", "\n", " A v_2 &= \\lambda_2 v_2\\\\\n", " \\left( A - \\lambda_2 I \\right) v_2 &= 0\\\\\n", " \\frac{1}{\\sqrt{\\alpha^2 + \\beta^2}}\n", " \\begin{bmatrix} i \\beta & -\\beta\\\\ \\beta & i \\beta \\end{bmatrix} v_2\n", " &= 0\n", "\n", "which shows that :math:`v_1 = \\begin{bmatrix} t_1 \\\\ -i t_1 \\end{bmatrix}` and\n", ":math:`v_2 = \\begin{bmatrix} t_2 \\\\ i t_2 \\end{bmatrix}` where\n", ":math:`t_1, t_2 \\in \\mathbb{R}` because each system of linear equations has a\n", "free variable.\n", "\n", "Defining\n", "\n", ".. math::\n", "\n", " Q =\n", " \\begin{bmatrix} v_1, v_2 \\end{bmatrix} =\n", " \\begin{bmatrix} t_1 & t_2 \\\\ -i t_1 & i t_2 \\end{bmatrix}\n", " \\quad \\text{and} \\quad\n", " Q^* = \\begin{bmatrix} t_1 & i t_1 \\\\ t_2 & -i t_2 \\end{bmatrix}\n", "\n", "gives\n", "\n", ".. math::\n", "\n", " Q Q^* &= \\begin{bmatrix} t_1 & t_2 \\\\ -i t_1 & i t_2 \\end{bmatrix}\n", " \\begin{bmatrix} t_1 & i t_1 \\\\ t_2 & -i t_2 \\end{bmatrix}\\\\\n", " &= \\begin{bmatrix}\n", " t_1^2 + t_2^2 & i t_1^2 - i t_2^2 \\\\\n", " -i t_1^2 + i t_2^2 & t_1^2 + t_2^2\n", " \\end{bmatrix}\n", "\n", "and\n", "\n", ".. math::\n", "\n", " Q^* Q &= \\begin{bmatrix} t_1 & i t_1 \\\\ t_2 & -i t_2 \\end{bmatrix}\n", " \\begin{bmatrix} t_1 & t_2 \\\\ -i t_1 & i t_2 \\end{bmatrix}\\\\\n", " &= \\begin{bmatrix}\n", " t_1^2 + t_1^2 & t_1 t_2 - t_1 t_2 \\\\\n", " t_1 t_2 - t_1 t_2 & t_2^2 + t_2^2\n", " \\end{bmatrix}\\\\\n", " &= \\begin{bmatrix}\n", " 2 t_1^2 & 0 \\\\\n", " 0 & 2 t_2^2\n", " \\end{bmatrix}.\n", "\n", "In order to satisfy the orthogonality constraint\n", ":math:`Q Q^* = Q^* Q = I`, set :math:`t_1 = t_2 = \\frac{1}{\\sqrt{2}}`.\n", "\n", ".. math::\n", "\n", " U &= Q^* A Q\\\\\n", " &= \\frac{1}{\\sqrt{\\alpha^2 + \\beta^2}} Q^*\n", " \\begin{bmatrix} \\alpha & -\\beta \\\\ \\beta & \\alpha \\end{bmatrix}\n", " \\begin{bmatrix} t_1 & t_2 \\\\ -i t_1 & i t_2 \\end{bmatrix}\\\\\n", " &= \\frac{1}{\\sqrt{\\alpha^2 + \\beta^2}}\n", " \\begin{bmatrix} t_1 & i t_1 \\\\ t_2 & -i t_2 \\end{bmatrix}\n", " \\begin{bmatrix}\n", " \\alpha t_1 + i \\beta t_1 & \\alpha t_2 - i \\beta t_2\\\\\n", " \\beta t_1 - i \\alpha t_1 & \\beta t_2 + i \\alpha t_2\n", " \\end{bmatrix}\\\\\n", " &= \\frac{1}{\\sqrt{\\alpha^2 + \\beta^2}}\n", " \\begin{bmatrix}\n", " \\alpha t_1^2 + i \\beta t_1^2 + i \\beta t_1^2 + \\alpha t_1^2 &\n", " \\alpha t_1 t_2 - i \\beta t_1 t_2 + i \\beta t_1 t_2 - \\alpha t_1 t_2\\\\\n", " \\alpha t_1 t_2 + i \\beta t_1 t_2 - i \\beta t_1 t_2 - \\alpha t_1 t_2 &\n", " \\alpha t_2^2 - i \\beta t_2^2 - i \\beta t_2^2 + \\alpha t_2^2\n", " \\end{bmatrix}\\\\\n", " &= \\frac{1}{\\sqrt{\\alpha^2 + \\beta^2}}\n", " \\begin{bmatrix}\n", " 2 \\alpha t_1^2 + 2 i \\beta t_1^2 & 0\\\\\n", " 0 & 2 \\alpha t_2^2 - 2 i \\beta t_2^2\n", " \\end{bmatrix}\\\\\n", " &= \\frac{\\sqrt{\\alpha^2 + \\beta^2}}{\\alpha^2 + \\beta^2}\n", " \\begin{bmatrix}\n", " \\alpha + i \\beta & 0 \\\\ 0 & \\alpha - i \\beta\n", " \\end{bmatrix}\n", "\n", "The Schur decomposition allowed :math:`\\hat{R}(\\alpha, \\beta)` to be computed\n", "more accurately because the values in :math:`U` are bounded between two\n", "conjugate values and :math:`Q, Q^*` can be computed to any decimal place via\n", "Newton's Method." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy\n", "\n", "def jacobi_rotation(n, i, j, alpha, beta):\n", " _ = numpy.sqrt(alpha**2 + beta**2)\n", " alpha_bar = alpha / _\n", " beta_bar = beta / _\n", "\n", " J = numpy.zeros((n,n))\n", " J[i,i] = J[j,j] = alpha_bar\n", " J[j,i] = beta_bar\n", " J[i,j] = -J[j,i]\n", "\n", " for k in range(n):\n", " if k != i and k != j:\n", " J[k,k] = 1\n", "\n", " return J" ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 22\n", "===========" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "M = numpy.asfarray([[1, 3, 5],\n", " [2, 4, 6],\n", " [3, 5, 9]])\n", "n = M.shape[0]\n", "\n", "#zero out M[2,0]\n", "J_il = jacobi_rotation(n, 2, 0, M[0,0], M[2,0])\n", "J_jl = jacobi_rotation(n, 2, 0, M[2,0], -M[0,0])\n", "print(numpy.dot(J_il, M))\n", "print(numpy.dot(J_jl, M))" ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 23\n", "===========\n", "\n", "If :math:`A` has full column rank, then the nullspace of :math:`A` is zero.\n", "---------------------------------------------------------------------------\n", "\n", ".. math::\n", "\n", " \\lVert Ax - b \\rVert^2_2 &= \\lVert A(x - x_*) + (A x_* - b) \\rVert^2_2\\\\\n", " &= \\left[ A(x - x_*) + (A x_* - b) \\right]^\\top\n", " \\left[ A(x - x_*) + (A x_* - b) \\right]\\\\\n", " &= \\lVert A (x - x_*) \\rVert^2_2 +\n", " \\lVert A x_* - b \\rVert^2_2 +\n", " 2 (x - x_*)^\\top A^\\top (A x_* - b)\\\\\n", " &= \\lVert A (x - x_*) \\rVert^2_2 + \\lVert A x_* - b \\rVert^2_2\n", " & \\quad & \\text{since}\\ A^\\top (A x_* - b) = 0\\\\\n", " \\lVert Ax - b \\rVert_2 &\\geq \\lVert A x_* - b \\rVert_2.\n", "\n", "This means for :math:`A v = 0` to be true, :math:`v = 0`; since\n", ":math:`x \\neq x_*`,\n", ":math:`\\lVert Ax - b \\rVert_2 > \\lVert A x_* - b \\rVert_2`.\n", "\n", ":math:`A^\\top A` is nonsingular if and only if :math:`A` has full column rank.\n", "------------------------------------------------------------------------------\n", "\n", "One way to prove this is to show that :math:`A^\\top A` has the same nullspace as\n", ":math:`A`. See :cite:`strang2009ila` for background on the fundamental theorem\n", "of linear algebra, subspaces, and orthogonality.\n", "\n", "If :math:`x` is in the nullspace of :math:`A`, then :math:`Ax = 0`. Multiplying\n", "by :math:`A^\\top` yields :math:`A^\\top A x = 0`, which shows that :math:`x` is\n", "also in the nullspace of :math:`A^\\top A`.\n", "\n", "If :math:`x` is in the nullspace of :math:`A^\\top A`, then\n", ":math:`A^\\top A x = 0`. Multiplying by :math:`x^\\top` yields\n", ":math:`x^\\top A^\\top A x = \\lVert A x \\rVert_2^2 = 0`, which means every\n", ":math:`x` in the nullspace of :math:`A^\\top A` also needs to be in the nullspace\n", "of :math:`A` in order to satisfy :math:`A x = 0`." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 24\n", "===========\n", "\n", "Since :math:`A \\in \\mathbb{R}^{m \\times n}` is full rank, its nullspace is zero\n", "and its columns are linearly independent. A set of vectors spans a space if\n", "their linear combinations fill the space i.e. the columns of :math:`A` span its\n", "column space. A basis for a vector space is a sequence of vectors that are\n", "linearly independent and span that space. Hence the column space of :math:`A`,\n", "whose dimension is :math:`n`, could span the entire space\n", ":math:`\\mathbb{R}^m` or only a subspace.\n", "\n", "Notice that the QR decomposition of\n", "\n", ".. math::\n", "\n", " A =\n", " \\begin{bmatrix} Q_1 & Q_2 \\end{bmatrix}\n", " \\begin{bmatrix} R_u \\\\ 0 \\end{bmatrix} =\n", " Q_1 R_u\n", "\n", "where :math:`Q_1 \\in \\mathbb{R}^{m \\times n}` and\n", ":math:`R_u \\in \\mathbb{R}^{n \\times n}`. Since :math:`Q` is an orthogonal\n", "matrix, the columns of :math:`Q_1` must be orthonormal. Hence the column space\n", "of :math:`Q_1`, whose dimension is :math:`n`, could span the entire space\n", ":math:`\\mathbb{R}^m` or only a subspace.\n", "\n", "By the theorem of Equal Dimensions Yields Equal Subspaces\n", ":cite:`beezer2008first`, the first :math:`n` columns of :math:`Q` form an\n", "orthonormal basis for the space spanned by the columns of :math:`A`." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 25\n", "===========" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "A = numpy.asfarray([[1, -1],\n", " [1, 0],\n", " [1, 1],\n", " [1, 2]])\n", "b = numpy.asfarray([[3],\n", " [2],\n", " [0],\n", " [4]])\n", "\n", "print(numpy.dot(numpy.dot(numpy.linalg.inv(numpy.dot(A.T, A)), A.T), b))\n", "\n", "n = A.shape[1]\n", "Q, R = QR_Householder(A)\n", "print(numpy.dot(numpy.linalg.inv(R[0:n,0:n]), numpy.dot(Q.T, b)[0:n,0:n]))" ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 26\n", "===========\n", "\n", "Given :math:`A = LQ`, the equivalent minimization is\n", "\n", ".. math::\n", "\n", " \\min_{x \\in \\mathbb{R}^n} \\left\\{ \\lVert x \\rVert_2 \\mid LQx = b\\right\\}.\n", "\n", "Since :math:`\\lVert Q x \\rVert_2 = \\lVert x \\rVert_2` (3.1.23), defining\n", ":math:`z = Q x` gives\n", "\n", ".. math::\n", "\n", " \\min_{z \\in \\mathbb{R}^n} \\left\\{ \\lVert z \\rVert_2 \\mid Lz = b\\right\\}.\n", "\n", "According to the analysis done in\n", ":ref:`Exercise 3.27 `,\n", ":math:`z = \\begin{bmatrix} L_l^{-1} b \\\\ 0 \\end{bmatrix}`. Thus\n", ":math:`x_* = Q^\\top z = Q^\\top \\begin{bmatrix} L_l^{-1} b \\\\ 0 \\end{bmatrix}`." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ ".. _dennis1996numerical-ex-3.27:\n", "\n", "Exercise 27\n", "===========\n", "\n", "Observe that the LQ decomposition of\n", "\n", ".. math::\n", "\n", " A =\n", " \\begin{bmatrix} L_l & 0 \\end{bmatrix}\n", " \\begin{bmatrix} Q_1 \\\\ Q_2 \\end{bmatrix} =\n", " L_l Q_1\n", "\n", "where :math:`L_l \\in \\mathbb{R}^{m \\times m}`,\n", ":math:`Q_1 \\in \\mathbb{R}^{m \\times n}`, and\n", ":math:`Q_2 \\in \\mathbb{R}^{(n - m) \\times n}`.\n", "\n", "To carry out LQ decomposition, apply QR decomposition to the full column rank\n", "matrix\n", "\n", ".. math::\n", "\n", " A^\\top =\n", " \\begin{bmatrix} Q_1^\\top & Q_2^\\top \\end{bmatrix}\n", " \\begin{bmatrix} L_l^\\top \\\\ 0 \\end{bmatrix}." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 28\n", "===========\n", "\n", ".. math::\n", "\n", " \\kappa_2(A) &= \\lVert A \\rVert_2 \\lVert A^{-1} \\rVert_2\\\\\n", " &= \\left(\n", " \\max_{v \\in \\mathbb{R}^n, v \\neq 0}\n", " \\frac{\\lVert A v \\rVert_2}{\\lVert v \\rVert_2}\n", " \\right)\n", " \\left(\n", " \\min_{v \\in \\mathbb{R}^n, v \\neq 0}\n", " \\frac{\\lVert A v \\rVert_2}{\\lVert v \\rVert_2}\n", " \\right)^{-1}\n", " & \\quad & \\text{(3.1.18)}\\\\\n", " &= \\frac{\n", " \\left( \\max_\\lambda \\left\\{ A^\\top A \\right\\} \\right)^{1/2}\n", " }{\n", " \\left( \\min_\\lambda \\left\\{ A^\\top A \\right\\} \\right)^{1/2}\n", " }\n", " & \\quad & \\text{(3.1.8b)}\\\\\n", " &= \\frac{\n", " \\left( \\max_\\lambda \\left\\{ V D^2 V^\\top \\right\\} \\right)^{1/2}\n", " }{\n", " \\left( \\min_\\lambda \\left\\{ V D^2 V^\\top \\right\\} \\right)^{1/2}\n", " }\n", " & \\quad & A = U D V^\\top\\\\\n", " &= \\frac{\\sigma_1}{\\sigma_n}\n", " & \\quad & \\text{Theorem 3.6.4}\n", "\n", "where :math:`\\sigma_1` is :math:`A`'s largest singular value and\n", ":math:`\\sigma_n` is :math:`A`'s smallest singular value.\n", "\n", "Let :math:`u_{\\cdot j}` and :math:`v_{\\cdot j}` be the columns of :math:`U` and\n", ":math:`V`. The following are the results of perturbing :math:`b`:\n", "\n", ".. math::\n", "\n", " A \\Delta x &= \\Delta b\\\\\n", " \\Delta x &= V D^{-1} U^\\top \\alpha u_{\\cdot j}\n", " & \\quad & \\Delta b = \\alpha u_{\\cdot j}\\\\\n", " &= V D^{-1}\n", " \\begin{bmatrix} 0 & \\cdots & \\alpha & \\cdots & 0 \\end{bmatrix}^\\top\\\\\n", " &= \\frac{\\alpha}{\\sigma_j} V\n", " \\begin{bmatrix} 0 & \\cdots & 1 & \\cdots & 0 \\end{bmatrix}^\\top\\\\\n", " &= \\frac{\\alpha}{\\sigma_j} v_{\\cdot j}." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ ".. rubric:: References\n", "\n", ".. bibliography:: chapter-03.bib" ] } ], "metadata": { "anaconda-cloud": {}, "celltoolbar": "Raw Cell Format", "kernelspec": { "display_name": "Python [default]", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.5.2" } }, "nbformat": 4, "nbformat_minor": 0 }