{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Sparse Matrix Factorizations and Fill-In\n",
        "Copyright (C) 2020 Andreas Kloeckner\n",
        "\n",
        "<details>\n",
        "<summary>MIT License</summary>\n",
        "Permission is hereby granted, free of charge, to any person obtaining a copy\n",
        "of this software and associated documentation files (the \"Software\"), to deal\n",
        "in the Software without restriction, including without limitation the rights\n",
        "to use, copy, modify, merge, publish, distribute, sublicense, and/or sell\n",
        "copies of the Software, and to permit persons to whom the Software is\n",
        "furnished to do so, subject to the following conditions:\n",
        "\n",
        "The above copyright notice and this permission notice shall be included in\n",
        "all copies or substantial portions of the Software.\n",
        "\n",
        "THE SOFTWARE IS PROVIDED \"AS IS\", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR\n",
        "IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,\n",
        "FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE\n",
        "AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER\n",
        "LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,\n",
        "OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN\n",
        "THE SOFTWARE.\n",
        "</details>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": []
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "collapsed": false,
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import scipy.linalg as la\n",
        "\n",
        "import matplotlib.pyplot as pt\n",
        "\n",
        "import random"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here's a helper routine to make a random **symmetric** sparse matrix:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "collapsed": false,
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [],
      "source": [
        "def make_random_sparse_matrix(n, row_fill):\n",
        "    nentries = (n*row_fill) // 2 # because of symmetry\n",
        "    data = np.random.randn(nentries)\n",
        "    rows = np.random.randint(0, n-1, nentries)\n",
        "    cols = np.random.randint(0, n-1, nentries)\n",
        "    \n",
        "    import scipy.sparse as sps\n",
        "    \n",
        "    coo = sps.coo_matrix((data, (rows, cols)), shape=(n,n))\n",
        "    \n",
        "    # NOTE: Cuthill-McKee applies only to symmetric matrices!\n",
        "    return (100*np.eye(n) + np.array(coo.todense() +  coo.todense().T))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we will take a look at that matrix from a \"birds eye view\". Every entry with absolute value greater that $10^{-10}$ will show up as a 'dot':"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "collapsed": false,
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [],
      "source": [
        "prec = 1e-10\n",
        "\n",
        "np.random.seed(15)\n",
        "random.seed(15)\n",
        "\n",
        "A = make_random_sparse_matrix(200, 3)\n",
        "print(\"%d non-zeros\" % len(np.where(np.abs(A)>prec)[0]))\n",
        "pt.figure()\n",
        "pt.spy(A, marker=\",\", precision=prec)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, let's apply the same visualization to the inverse:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "collapsed": false,
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [],
      "source": [
        "Ainv = la.inv(A)\n",
        "print(\"%d non-zeros\" % len(np.where(np.abs(Ainv) > prec)[0]))\n",
        "pt.spy(Ainv, marker=\",\", precision=prec)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And the Cholesky factorization:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {
        "collapsed": false,
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [],
      "source": [
        "L = la.cholesky(A)\n",
        "print(\"%d non-zeros\" % len(np.where(np.abs(L) > prec)[0]))\n",
        "pt.spy(L, marker=\",\", precision=prec)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Cholesky is often less bad, but in principle affected the same way."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Reducing the fill-in"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Define the *degree* of a row as the number of non-zeros in it."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {
        "collapsed": false,
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [],
      "source": [
        "def degree(mat, row):\n",
        "    return len(np.where(mat[row])[0])\n",
        "\n",
        "print(degree(A, 3))\n",
        "print(degree(A, 4))\n",
        "print(degree(A, 5))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Then find an ordering so that all the low degrees come first.\n",
        "\n",
        "The [Cuthill-McKee algorithm](https://en.wikipedia.org/wiki/Cuthill%E2%80%93McKee_algorithm) is a greedy algorithm to find such an ordering:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 29,
      "metadata": {
        "collapsed": false,
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [],
      "source": [
        "def argmin2(iterable, return_value=False):\n",
        "    it = iter(iterable)\n",
        "    try:\n",
        "        current_argmin, current_min = next(it)\n",
        "    except StopIteration:\n",
        "        raise ValueError(\"argmin of empty iterable\")\n",
        "\n",
        "    for arg, item in it:\n",
        "        if item < current_min:\n",
        "            current_argmin = arg\n",
        "            current_min = item\n",
        "\n",
        "    if return_value:\n",
        "        return current_argmin, current_min\n",
        "    else:\n",
        "        return current_argmin\n",
        "\n",
        "def argmin(iterable):\n",
        "    return argmin2(enumerate(iterable))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 30,
      "metadata": {
        "collapsed": false,
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [],
      "source": [
        "def cuthill_mckee(mat):\n",
        "    \"\"\"Return a Cuthill-McKee ordering for the given matrix.\n",
        "\n",
        "    See (for example)\n",
        "    Y. Saad, Iterative Methods for Sparse Linear System,\n",
        "    2nd edition, p. 76.\n",
        "    \"\"\"\n",
        "\n",
        "    # this list is called \"old_numbers\" because it maps a\n",
        "    # \"new number to its \"old number\"\n",
        "    old_numbers = []\n",
        "    visited_nodes = set()\n",
        "    levelset = []\n",
        "\n",
        "    all_nodes = set(range(len(mat)))\n",
        "\n",
        "    while len(old_numbers) < len(mat):\n",
        "        if not levelset:\n",
        "            unvisited = list(all_nodes - visited_nodes)\n",
        "\n",
        "            if not unvisited:\n",
        "                break\n",
        "\n",
        "            start_node = unvisited[\n",
        "                    argmin(degree(mat, node) for node in unvisited)]\n",
        "            visited_nodes.add(start_node)\n",
        "            old_numbers.append(start_node)\n",
        "            levelset = [start_node]\n",
        "\n",
        "        next_levelset = set()\n",
        "        levelset.sort(key=lambda row: degree(mat, row))\n",
        "        #print(levelset)\n",
        "        \n",
        "        for node in levelset:\n",
        "            row = mat[node]\n",
        "            neighbors, = np.where(row)\n",
        "            \n",
        "            for neighbor in neighbors:\n",
        "                if neighbor in visited_nodes:\n",
        "                    continue\n",
        "\n",
        "                visited_nodes.add(neighbor)\n",
        "                next_levelset.add(neighbor)\n",
        "                old_numbers.append(neighbor)\n",
        "\n",
        "        levelset = list(next_levelset)\n",
        "\n",
        "    return np.array(old_numbers, dtype=np.intp)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "metadata": {
        "collapsed": false,
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [],
      "source": [
        "cmk = cuthill_mckee(A)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Someone (empirically) observed that the *reverse* of the Cuthill-McKee ordering often does better than forward Cuthill-McKee.\n",
        "\n",
        "So construct a permutation matrix corresponding to that:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 31,
      "metadata": {
        "collapsed": true,
        "jupyter": {
          "outputs_hidden": true
        }
      },
      "outputs": [],
      "source": [
        "P = np.eye(len(A))[cmk[::-1]]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And then reorder both rows and columns according to that--a similarity transform:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 32,
      "metadata": {
        "collapsed": false,
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [],
      "source": [
        "A_reordered = P @ A @ P.T\n",
        "\n",
        "pt.spy(A_reordered, marker=\",\", precision=prec)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, let's try Cholesky again:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 33,
      "metadata": {
        "collapsed": false,
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [],
      "source": [
        "L = la.cholesky(A_reordered)\n",
        "print(\"%d non-zeros\" % len(np.where(np.abs(L) > prec)[0]))\n",
        "pt.spy(L, marker=\",\", precision=prec)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false,
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [],
      "source": []
    }
  ],
  "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.8.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}