lecture-13.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{u'start_slideshow_at': 'selected', u'theme': 'sky', u'transition': 'zoom'}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.html.services.config import ConfigManager\n",
    "from IPython.utils.path import locate_profile\n",
    "cm = ConfigManager(profile_dir=locate_profile(get_ipython().profile))\n",
    "cm.update('livereveal', {\n",
    "              'theme': 'sky',\n",
    "              'transition': 'zoom',\n",
    "              'start_slideshow_at': 'selected',\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#Lecture 13. The algebra of hierarchical matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Previous lecture \n",
    "- Introduction to hierarhical matrices ($\\mathcal{H}, \\mathcal{H}^2$) as algebraic interpretation of the FMM\n",
    "- The concept of block rows and nested bases\n",
    "- The concept of splitting of the matrix into blocks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Todays lecture\n",
    "- How to construct the hierarhical approximation (both in H- and H-2 cases)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Book\n",
    "A good introductory book is S. Borm \"Efficient numerical methods for non-local operators: H2-matrix compression, algorithms and analysis\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hierarhical matrices\n",
    "\n",
    "- Split the matrix into blocks $A(t, s)$ corresponding to the mosaic partioning, approximate \"far\" blocks with low-rank matrices.\n",
    "- **$H^2$** matrix: using the block row (i.e. interaction of the box with everything outside) is of low-rank.\n",
    "- Computation of the factorization requires the treatment of block matrices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple case: H-matrices\n",
    "The case of H-matrices is simple: the matrix is represented as a collection of low-rank blocks, so we have to approximate each block independently, \n",
    "\n",
    "$$A(t, s) \\approx U(t, s) V(t, s)^{\\top}.$$\n",
    "\n",
    "\n",
    "How we can do that?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## NLA Flashback: approximation of low-rank matrices\n",
    "A rank-$r$ matrix can be represented as\n",
    "\n",
    "$$A = C \\widehat{A}^{-1} R, $$\n",
    "\n",
    "where $C$ are some **columns** of the matrix $A$, $R$ are some rows of the matrix $A$,  $\\widehat{A}$ is the submatrix on their intersection.\n",
    "\n",
    "Approximate case:\n",
    "\n",
    "If $\\widehat{A}$ is the submatrix with **maximal volume** (where volume is the absolute value of the determinant)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Cross-approximation\n",
    "Idea of the cross approximation: select the submatrix to maximize the determinant, i.e. in a **greedy** fashion.\n",
    "\n",
    "The term \"cross\" comes from the rank-$1$ update formula\n",
    "\n",
    "$$A_{ij} := A_{ij} - \\frac{A_{i j^*} A_{i^* j}}{A_{i^* j^*}},$$\n",
    "\n",
    "where the **pivots** $(i^*, j^*)$ has to be selected in such a way that $|A_{i^* j^*}|$ is as big as possible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Pivoting strategies\n",
    "- Row/column pivoting (select a column, find maximal in it).\n",
    "- Rook pivoting\n",
    "- Additionally, do some random sampling\n",
    "- There is a result by L. Demanet et. al on the class of matrices, where the approximation exist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Concept of approximation for H-matrices\n",
    "1. Create a list of blocks\n",
    "2. Sample rows/columns to get the low-rank factorization\n",
    "3. Profit!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## $H^2$-matrices\n",
    "For the $H^2$ matrices the situation is much more complicated.\n",
    "\n",
    "The standard way to go is to first compress into the $\\mathcal{H}$-form, and then do **recompression** into the $\\mathcal{H}^2$ form.  \n",
    "\n",
    "Such recompression can be done in $\\mathcal{O}(n \\log n)$ operations.\n",
    "\n",
    "Can we do it directly?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Nested cross approximation\n",
    "\n",
    "- A block row is of low-rank -> there exist **basis rows** that span the row space.\n",
    "- If we join the basis rows from the children, we get the basis rows for the father.\n",
    "- This requires the SVD of the matrix $2r \\times N$, and it has $\\mathcal{O}(N^2)$ complexity (although better, than $\\mathcal{O}(N^3)$ for direct SVD of big blocks)\n",
    "\n",
    "<font color='red'> Solution: </font> Approximate the \"far zone\" with few receivers (\"representor set\")."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Demo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.95954990387\n",
      "H2 1e-5\n",
      "('Function calls:', 36682)\n",
      "('Function values computed:', 33309887)\n",
      "('Function time:', 1.7329978942871094)\n",
      "('Average time per function value:', 5.202653177079554e-08)\n",
      "('Maxvol time:', 6.476972341537476)\n",
      "('Total MCBH time:', 9.493279933929443)\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "os.environ['OMP_NUM_THREADS'] = '1'\n",
    "os.environ['MKL_NUM_THREADS'] = '1'\n",
    "\n",
    "from h2py.tree.quadtree import ParticlesQuadTree as PQT\n",
    "from h2py.tree.octtree import ParticlesOctTree as POT\n",
    "from h2py.tree.inertial_tree import ParticlesInertialTree as PIT\n",
    "from h2py.data.particles_data import ParticlesData as PD\n",
    "from h2py.problem import Problem\n",
    "from h2py.data.particles import log_distance\n",
    "inv_distance = log_distance\n",
    "import numpy as np\n",
    "from time import time\n",
    "import sys\n",
    "\n",
    "N = 20000\n",
    "np.random.seed(0)\n",
    "particles = np.random.rand(2, N)\n",
    "data = PD(particles)\n",
    "tree = PIT(data, block_size = 20)\n",
    "problem = Problem(inv_distance, tree, tree)\n",
    "problem.build()\n",
    "problem.gen_queue(symmetric=0)\n",
    "\n",
    "print('H2 1e-5')\n",
    "problem.factorize('h2', tau=1e-5, onfly=0, verbose=1, iters=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Representor set\n",
    "- Way 1: Select it using apriori knowledge (geometrical approach)\n",
    "- Way 2: For \"good columns\" it is sufficient to know good columns of the father! \n",
    "For details, see http://arxiv.org/abs/1309.1773\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Inversion of the hierarchical matrices\n",
    "Recall our goal is often to solve the integral equation (i.e., compute the inverse).\n",
    "\n",
    "One of the possible ways is to use **recursive block elimination**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Block-LU (or Schur complement)\n",
    "Consider the matrix\n",
    "$$A = \\begin{bmatrix} A_{11} & A_{12} \\\\ A_{21} & A_{22}\\end{bmatrix},$$\n",
    "\n",
    "where the first group of variables correspond to the unknowns in the first node, the second group variables corresponds to the unknowns in the second node (binary tree implicitly assumed).\n",
    "\n",
    "Then, \n",
    "\n",
    "$$A \\begin{bmatrix} q_1 \\\\ q_2 \\end{bmatrix} = \\begin{bmatrix} f_1 \\\\ f_2 \\end{bmatrix}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "After the elimination we have the following equality for the **Schur complement**.  \n",
    "\n",
    "$$q_1 = A^{-1}_{11} f_1 - A^{-1}_{11} A_{12} q_2, \\quad \\underbrace{(A_{22} - A_{21} A^{-1}_{11} A_{12})}_{\\mbox{Schur complement}} q_2 = f_2 - A_{21} A^{-1}_{11} f_1.$$\n",
    "\n",
    "The core idea is **recursion**: if we know $A^{-1}_{11},$ then we can compute the matrix $S$, and invert it as well.\n",
    "\n",
    "The multiplication of H-matrices has $\\mathcal{O}(N \\log N)$ complexity, and is also (typically) implemented via recursion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Multiplication of H-matrices\n",
    "\n",
    "Consider 1D partioning, and multiplication of two matrices with H-structure (i.e., blocks(1,2) and (2,1) ) have low-rank\n",
    "\n",
    "$$\\begin{bmatrix} A_{11} & A_{12} \\\\ A_{21} & A_{22}\\end{bmatrix}\\begin{bmatrix} B_{11} & B_{12} \\\\ B_{21} & B_{22}\\end{bmatrix}$$\n",
    "\n",
    "The (1, 2) block of the result is \n",
    "\n",
    "$$\n",
    "   \\underbrace{A_{21} B_{11}}_{\\mbox{low rank}} + \\underbrace{A_{22} B_{21}}_{\\mbox{low rank}}\n",
    "$$\n",
    "\n",
    "(1, 1) and (2, 2) blocks are evaluated **recursively**, so it is a recursion inside recursion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary\n",
    "- Use Block elimination (by nodes of the tree)\n",
    "- Problem is reduced to the multiplication of $H$-matrices\n",
    "- Constant is high\n",
    "- Can compute $LU$-factorization instead."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Fast direct solvers for sparse matrices\n",
    "Sparse matrices coming from PDEs are in fact H-matrices as well!\n",
    "\n",
    "So, we can compute the inverse by the same block factorization technique.\n",
    "\n",
    "Works very well for the \"1D\" partioning (i.e., off-diagonal blocks are of low-rank), does not work for 2D/3D problems with optimal complexity (but the constants can be really good).\n",
    "\n",
    "We also have our own idea (an since it is unpublished, will use the whiteboard magic here :) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary \n",
    "- Nested cross approximation\n",
    "- Block Schur elimination idea for the inversion\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Next lecture (week)\n",
    "- We will talk about high-frequency problems (and there are problems there)\n",
    "- FFT on the non-uniform grid, butterfly algorithm\n",
    "- Precorrected FFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "    }\n",
       "    div.cell{\n",
       "        /*width:80%;*/\n",
       "        /*margin-left:auto !important;\n",
       "        margin-right:auto;*/\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\n",
       "    h2 {\n",
       "        font-family: 'Fenix', serif;\n",
       "    }\n",
       "    h3{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "        margin-top:12px;\n",
       "        margin-bottom: 3px;\n",
       "       }\n",
       "\th4{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "       }\n",
       "    h5 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\t   \n",
       "    div.text_cell_render{\n",
       "        font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 1.2;\n",
       "        font-size: 120%;\n",
       "        /*width:70%;*/\n",
       "        /*margin-left:auto;*/\n",
       "        margin-right:auto;\n",
       "    }\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 50pt;\n",
       "\t\tline-height: 110%;\n",
       "        color:#CD2305;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\t\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 300;\n",
       "        font-size: 16pt;\n",
       "        color: #CD2305;\n",
       "        font-style: italic;\n",
       "        margin-bottom: .5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\n",
       "    \n",
       "    li {\n",
       "        line-height: 110%;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}