lecture-12.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 12. Algebraic nature of the FMM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Previous lecture \n",
    "- Barnes-Hut algorithm, single tree\n",
    "- Barnes-Hut algorithm, double tree\n",
    "- Types of trees (octtree, KD-tree, recursive inertial bisection)\n",
    "- The bottom-top-left-right-top-bottom idea of the FMM\n",
    "- The concept of the multipole expansion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Todays lecture\n",
    "- Complexity analysis for the FMM\n",
    "- M2M, M2L, L2L operators via Taylor expansions\n",
    "- Separability criteria and low-rank approximation\n",
    "- What is the algebraic interpretation of the Barnes-Hut method\n",
    "- What is the algebraic interpretation of the FMM\n",
    "- The concept of equivalent charges\n",
    "- $\\mathcal{H}$, $\\mathcal{H}^2$-matrices\n",
    "- A brief recap on the uniform grid computations "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## N-body problem\n",
    "Write down again the basic N-body problem.\n",
    "\n",
    "$$V_i = \\sum_{j} q_j F(x_i, y_j).$$\n",
    "\n",
    "- We construct trees for the **sources** $y_j$ and **receivers** $x_i$.\n",
    "- We construct M2M, M2L, L2L operators\n",
    "- We fill the \"multipoles\" in the sources tree using M2M operators (**upward pass**)\n",
    "- We fill the \"local\" coefficients in the receivers tree using M2L operators (**translation operator**)\n",
    "- We fill the \"local\" coefficients using L2L operators (**downward pass**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Tree notation\n",
    "\n",
    "We will denote the tree as $\\mathcal{T}$, its nodes as $t$. \n",
    "\n",
    "The set **children(t)** is the set of all children of a given node (is empty for the leaf)  \n",
    "\n",
    "For each $t$ we will denote all pairs $(t, s)$ for all nodes $s$ that belong to the **interaction list**. \n",
    "\n",
    "Definition: For box $X$ the white squares are the **interaction list**: the fathers do not fullfil the separability criteria, but the children are separable)\n",
    "\n",
    "<img src='interaction2.png' width=60%>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## M2M operator\n",
    "\n",
    "M2M operator has a very simple form for the Barnes-Hut algorithm. \n",
    "\n",
    "Each node stores exactly **1 number** (total charge/mass),  \n",
    "\n",
    "Thus for a binary tree\n",
    "\n",
    "$$q_{father} = \\underbrace{\\begin{bmatrix} 1 & 1 \\end{bmatrix}}_{\\mbox{M2M operator}} \\begin{bmatrix} q_{leaf1} \\\\ q_{leaf2}\\end{bmatrix}$$\n",
    "\n",
    "In general, \n",
    "\n",
    "$$\n",
    "   M_{father} = M2M \\begin{bmatrix} M_{leaf1} \\\\ M_{leaf2}\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## M2L operator\n",
    "\n",
    "For a given node $s$ of the source cluster  in the **interaction list** \n",
    "we store a matrix $\\mbox{M2L(t, s)}$:\n",
    "\n",
    "$$L(s) = \\mbox{M2L}(t, s) M(t).$$\n",
    "\n",
    "For the Barnes-Hut example, $M2L$ is a $1 \\times 1$ matrix with the element $\\frac{1}{\\Vert R_t - R_s \\Vert}$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## L2L operator\n",
    "L2L operator is also a $1 \\times 1$ matrix for the BH!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Complexity estimate of the FMM\n",
    "A simple complexity estimate for the FMM in 2D gives  \n",
    "\n",
    "$$\\underbrace{N p}_{\\mbox{upward pass}} + \\underbrace{29 \\frac{N}{s} p^2}_{\\mbox{translation}} + \\underbrace{N p}_{\\mbox{downward pass}} + \\underbrace{9 N s}_{\\mbox{close interaction}}$$\n",
    "\n",
    "A good choice is $s \\approx p$ and that gives the complexity $40 N p$ ($p$ is the number of multipole coefficients involved)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The multipole method\n",
    "In the multipole method, we approximate the interaction between two separated boxes with the Taylor series, or in other words, the function $f(x, y)$ is approximated in the strip  \n",
    "\n",
    "$$\\Vert x - y \\Vert \\geq \\delta.$$\n",
    "\n",
    "\n",
    "<font color='red'> Taylor series is not optimal </font>\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Separated representations\n",
    "What we have from the Taylor series? \n",
    "\n",
    "We have the separated representation, where the separation is of \"sources\" and \"receivers\"\n",
    "\n",
    "$$f(x, y) \\approx \\sum_{\\alpha=1}^r u_{\\alpha}(x) v_{\\alpha}(y),$$\n",
    "\n",
    "where $x, y$ are in fact two-dimensional (or three-dimensional) vectors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Separated representation and fast methods\n",
    "If the interaction between two boxes is separable, then  \n",
    "\n",
    "$$V(x) = \\sum_{j} f(x, y_j) q_j = \\sum_{\\alpha=1}^r u_{\\alpha}(x) c_{\\alpha},$$\n",
    "\n",
    "where the coefficients $c_{\\alpha}$ are given by the formula\n",
    "\n",
    "$$c_{\\alpha} = \\sum_{j} v_{\\alpha}(y_j) q_j.$$\n",
    "\n",
    "\n",
    "Questions:\n",
    "\n",
    "- <font color='red'> Is this enough for BH-type complexity? </font>  \n",
    "\n",
    "- <font color='red'> Is this enough for the FMM-type complexity? </font>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Answers\n",
    "- It is enough for the BH-type complexity on the matvec stage (when the decomposition is known)\n",
    "- It is not enough for the FMM-type complexity, since we have to compute the coefficients and then add to all the points inside the box."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Algebraic interpretation \n",
    "Algebraic interpretation of the separated representation is simple.  \n",
    "\n",
    "Given two nodes $t, s$ of the cluster trees of sources and receivers, the interaction between boxes containing $t$ and $s$ is in fact the interaction between **all particles** in $t$ and **all particles** in $s$.\n",
    "\n",
    "This is then a (sub)matrix of the total interaction matrix $A$.  \n",
    "\n",
    "It is convenient to denote it by $A(t, s)$. The separability property is then reduced on the discrete level to \n",
    "\n",
    "$$A(t, s) \\approx U(t, s) V(t, s)^{\\top}, \\quad U(t, s) \\in \\mathbb{R}^{M_t \\times r}, \\quad V(t, s) \\in \\mathbb{R}^{N_s \\times r}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Algebraic interpretation (cont). \n",
    "The low-rank factorization of the separated blocks leads to the **hierarchical matrices** ($\\mathcal{H}$-matrices, which were actually introduced by Tyrtyshnikov in 1993, and later (re)introduced by Hackbusch, Khoromskij in 1999."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We split the matrix into blocks, and approximate blocks $A(t, s)$ that satisfy the **admissability condition** by low-rank matrices.\n",
    "\n",
    "<img src=\"h-matrix.jpg\" width=70%>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## This is not enough!\n",
    "The interpretation by \"block low-rank structure\" is viable, but not enough;\n",
    "\n",
    "A closer step to the FMM algorithm is to consider the decomposition\n",
    "\n",
    "\n",
    "$$A(t, s) = U(t) \\Psi(t, s) V^{\\top}(s),$$\n",
    "\n",
    "where we have **common factor matrices** $U(t)$ and $V(s)$ for each node, and coefficient matrices $\\Psi(t, s)$ which are **small**.\n",
    "\n",
    "Still not an FMM: will have to store $n \\times r$ matrices (there are nowhere in the FMM)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The final step: basis nestedness\n",
    "\n",
    "To get the matrix structure, that is algebraically equivalent to the FMM structure, we have to make one additional assumption: the bases $U(t)$ are <font color='red'> nested </font>\n",
    "\n",
    "In order to figure out what it means, we will have to define the **block row** of the matrix.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Block columns\n",
    "\n",
    "All the blocks from the **interaction list** form \"own far zone\". If we add the far zone of the father of the node, and of grandfather, then it will be the **block row** of the full **far zone** (i.e. everything except for the close zone)\n",
    "\n",
    "<img src='blockstring.png ' width=70%>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Algebraic characterization\n",
    "The algebraic characterization is that the **rank of the block row** is equal to $r$.\n",
    "\n",
    "Then the bases are nested, i.e. from $U(t_1)$ and $U(t_2)$ we can recompute the basis for the father:\n",
    "\n",
    "$$U(t) = \\begin{bmatrix} U_{t_1} & 0 \\\\\n",
    "                          0 & U_{t_2} \\end{bmatrix} \\Phi((t_1,t_2)\\rightarrow t).$$\n",
    "                          \n",
    "The transfer matrix $\\Phi$ is of size $(r_{t_1} + r_{t_2}) \\times r_t \\approx 2r \\times r$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## $H^2$-matrices\n",
    "- Thus, we only need the matrices to recompute the basises from the sons to the father.  \n",
    "- The same holds for the columns.\n",
    "- The transfer matrices $\\Phi$ are the M2M operators, the matrix $\\Psi(t, s)$ are the M2L operators, and the same column matrices are the $L2L$ operators!\n",
    "\n",
    "Let us do some whiteboard magic to show that it is indeed the case.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Speeding up translation operators\n",
    "Now when we know the algebraic interpretation of the FMM, we can easily describe how to speedup the computation of translation operators.\n",
    "\n",
    "The **translation** is just the multiplication by the matrices $\\Psi(t, s)$. \n",
    "\n",
    "Thus, if we select such basis $U'(t)$ and $V'(s)$ (maybe even larger, that the minimal one!) so $\\Psi(t, s)$ is **diagonal**, we are fine."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Diagonal translation operators\n",
    "Idea: the basis functions are created by charges located in the **uniform mesh** inside the box.\n",
    "\n",
    "Then, the interaction becomes the interaction between two uniform grids in two boxes, and \n",
    "it can be computed with the help of fast Toeplitz matrix-by-vector product."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Final interpretation of $H^2$-matrices\n",
    "\n",
    "If hierarchical matrices is the representation of a matrix as a set of **non-intersecting** low-rank blocks, the $\\mathcal{H}^2$-matrix structure is the representation as the set of **intesecting** blocks (block rows and columns) which are of low-rank. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Other types of splittings\n",
    "The splitting can be quite interesting, for example we can consider the simplest \"interaction\"\n",
    "\n",
    "$$F(x, y) = 1, \\quad \\Vert x - y \\Vert > \\delta, \\quad 0, \\quad \\mbox{otherwise}.$$\n",
    "\n",
    "Then in 1D this is the **excluded sum**.\n",
    "\n",
    "$$V_i = \\sum_{j \\ne i} q_j.$$\n",
    "\n",
    "The most efficient way is to compute the **suffix sum** and the **prefix sum**.\n",
    "\n",
    "This also can be interpreted as the splitting of the matrix into blocks!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Important questions\n",
    "- How you compute the factorization? \n",
    "- When such approximation exist?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## H-matrix software\n",
    "There is a well-known code HLIB for working with H-matrices. There is also **kernel-independent** variant. \n",
    "\n",
    "\n",
    "Several \"in-house\" codes, we also develop one (http://bitbucket.org/muxas/h2tools)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Overview of todays lecture\n",
    "\n",
    "- Algebraic interpretation of the FMM as block-low rank matrices\n",
    "- The nestedness property"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Next lecture\n",
    "- The computation of the H^2 factorization\n",
    "- Algebra on hierarchical matrices\n",
    "- Inversion of hierarchical matrices\n",
    "- Fast low-rank direct solvers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bonus track\n",
    "What happens, if the particles are on the uniform grid?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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: 100%;\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: 100%;\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: 100%;\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": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/alex.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "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
}