{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Sherman-Morrison"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import scipy.linalg as la"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's set up some matrices and data for the *rank-one modification*:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n = 5\n",
        "A = np.random.randn(n, n)\n",
        "u = np.random.randn(n)\n",
        "v = np.random.randn(n)\n",
        "\n",
        "b = np.random.randn(n)\n",
        "\n",
        "Ahat = A + np.outer(u, v)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's start by computing the \"base\" factorization.\n",
        "\n",
        "We'll use `lu_factor` from `scipy`, which stuffs both `L` and `U` into a single matrix (why can it do that?) and also returns pivoting information:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[[ 1.08364401 -2.18363086 -1.07875426 -0.17119634  0.40477559]\n",
            " [ 0.94001453  2.11377723  0.49515679  0.04580586 -0.48416682]\n",
            " [ 0.11060844 -0.01756305 -1.72250622 -0.49623917 -1.42697035]\n",
            " [-0.2463909  -0.41778398  0.78182516  1.19130607  1.03773055]\n",
            " [-0.06177908 -0.79834934 -0.1559638   0.71449418 -0.04804558]]\n",
            "[3 3 2 4 4]\n"
          ]
        }
      ],
      "source": [
        "LU, piv = la.lu_factor(A)\n",
        "print(LU)\n",
        "print(piv)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we set up a subroutine to solve using that factorization and check that it works:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "3.074360713233696e-15"
            ]
          },
          "execution_count": 9,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "def solveA(b):\n",
        "    return la.lu_solve((LU, piv), b)\n",
        "\n",
        "la.norm(np.dot(A, solveA(b)) - b)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As a last step, we try the Sherman-Morrison formula:\n",
        "\n",
        "$$(A+uv^T)^{-1} = A^{-1} - {A^{-1}uv^T A^{-1} \\over 1 + v^T A^{-1}u}$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To see that we got the right answer, we first compute the right solution of the modified system:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "xhat = la.solve(Ahat, b)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, apply Sherman-Morrison to find `xhat2`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "xhat2 = solveA(b) - solveA(u)*np.dot(v, solveA(b))/(1+np.dot(v, solveA(u)))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "5.010784718882614e-15"
            ]
          },
          "execution_count": 14,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "la.norm(xhat - xhat2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "* What's the cost of the Sherman-Morrison procedure?"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "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.1+"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}