Add theory sections

docs/source/glossary.rst

@@ -0,0 +1,93 @@
+Glossary
+========
+
+.. glossary::
+
+    HODLR
+        Hierarchical Off-Diagonal Low-Rank matrix. See 
+        :doc:`theory/hodlr_theory` for more information.
+
+    low-rank matrix
+    low-rank format
+        Matrix representation obtained by approximating a matrix using a 
+        truncated singular value decomposition. For more information, see
+        :ref:`low-rank-explanation`.
+
+    tree
+        Tree-shaped structure used to represent the hierarchical nature of a 
+        :term:`HODLR` matrix. Composed of :term:`nodes` connected in a 
+        tree-like shape. For more information, see 
+        :doc:`theory/library/hodlr`.
+
+    node
+    nodes
+        Basic unit of a :term:`tree` data structure - multiple nodes connected
+        together in a tree shape make up a :term:`tree`. Depending on its 
+        type, a node may either hold data or connect to children nodes.
+
+    height
+        The number of edges on the longest path from the :term:`root node` to
+        the bottommost :term:`leaf node`. Also, the number of :term:`levels` 
+        composing a :term:`tree`. E.g., a tree with one :term:`root node`
+        and four :term:`children` will have a height of 1.
+
+    level
+    levels
+        The number of edges on the longest path from the :term:`root node` to
+        a particular :term:`node`. Nodes that are the same number of edges 
+        away from the :term:`root` are on the same level. E.g., the 
+        :term:`root node` is always at ``level==0``, its children are at
+        ``level==1``, etc.
+
+    root
+    root node
+        The topmost :term:`node` of a :term:`tree`, i.e. its beginning. Has no 
+        :term:`parent`. In ``hmat_lib``, a root node is always an 
+        :term:`internal node`.
+
+    leaf
+    leaves
+    leaf node
+    leaf nodes
+        A terminal :term:`node`, i.e. its end. Has no :term:`children`.
+
+    internal node
+    internal nodes
+        A :term:`node` that has one or more :term:`children`. In ``hmat_lib``,
+        an internal node does not store any data. For more information, see
+        :ref:`HODLR structure explanation<internal-node-explanation>`
+
+    diagonal node
+    diagonal nodes
+    diagonal leaf node
+    diagonal leaf nodes
+        A :term:`leaf node` which represents a dense block on the diagonal
+        of the :term:`HODLR` matrix. For more information, see
+        :ref:`HODLR structure explanation<diagonal-node-explanation>`
+
+    off-diagonal node
+    off-diagonal nodes
+    off-diagonal leaf node
+    off-diagonal leaf nodes
+        A :term:`leaf node` which represents a low-rank block off the diagonal
+        of the :term:`HODLR` matrix. For more information, see
+        :ref:`HODLR structure explanation<offdiagonal-node-explanation>`
+
+    parent
+    parent node
+        A :term:`node` that is the ancestor of another :term:`node` (its 
+        :term:`child`). In ``hmat_lib``, a parent node is always an 
+        :term:`internal node`.
+
+    child
+    children
+    child node
+    child nodes
+        A :term:`node` that descends from another :term:`node` (its 
+        :term:`parent`)
+
+    subtree
+        A :term:`tree` formed by a :term:`node` and all its descendants.
+
+
+

docs/source/index.rst

@@ -10,6 +10,8 @@ matrices, with a focus on HODLR-HODLR matrix-matrix multiplication.
    :caption: Contents:
 
    installation
+   theory
+   glossary
    api
    dev
 

docs/source/theory.rst

@@ -0,0 +1,12 @@
+Theory
+======
+
+These pages explain various background aspects of ``hmat_lib``.
+
+
+.. toctree::
+   :maxdepth: 1
+   :glob:
+
+   theory/*
+

docs/source/theory/hmat_lib.rst

@@ -0,0 +1,95 @@
+hmat_lib Library
+================
+
+How does the library work?
+--------------------------
+
+The below pages explain *how* ``hmat_lib`` works, describing the data 
+structures and how they relate to :term:`HODLR`. For, instead, *why* the 
+library works the way it works, read on in the 
+:ref:`next section<design-decisions>`.
+
+.. toctree::
+   :maxdepth: 1
+   :glob:
+
+   library/*
+
+
+.. _design-decisions:
+
+Design decisions
+----------------
+
+``hmat_lib`` was designed with HODLR-HODLR matrix-matrix multiplication in 
+mind. In the process of its implementation, several design choices have been
+made that influenced the process. These and their reasons are detailed here:
+
+HODLR as a perfectly balanced tree
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The HODLR matrix in ``hmat_lib`` is represented via a perfectly balanced 
+:term:`tree`. For more details on the data structure and how it relates to 
+HODLR matrices, see :doc:`library/hodlr` - here we discuss the reasons for 
+and consequences of this decision. In this case, the reasons for this 
+decision are quite straightforward:
+
+1. Simplicity
+
+   * It is simpler to conceptualise and work with a :term:`tree` that is 
+     perfectly balanced (+ fewer things to unit test!) - this way there are no 
+     extraneous cases and it is simpler and more efficient to 
+     :ref:`iterate over<tree-iteration>`.
+
+2. Practicality
+
+   * In many :doc:`applications<why_hodlr>` when a matrix is converted into
+     the HODLR format, the HODLR ends up fairly balanced - the ranks on all
+     off-diagonal nodes, both between different levels and within a level,
+     tend to be similar. A perfectly balanced :term:`tree` represents such
+     arrangement well.
+
+The consequence of this decision is that only uniformly deep :term:`HODLR` 
+matrices can be represented in ``hmat_lib``. Fortunately, this is the default,
+and if a more complex arrangement ever becomes of interest, it should be 
+possible to extend the current framework.
+
+
+.. _hodlr-always-square:
+
+HODLR matrix is always square
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The :c:struct:`TreeHODLR` data structure has been designed to always represent
+a square :term:`HODLR` matrix - not only are there no routines for converting
+a rectangular matrix into the :term:`HODLR` format, with the current 
+``struct``\ s, it is impossible to generate one at all. This is because, 
+again, there is limited interest in such an arrangement and would require
+:ref:`additional work<diagonal-always-square>`.
+
+
+.. _diagonal-always-square:
+
+Diagonal blocks are square matrices
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The :c:struct:`TreeHODLR` data structure has also been designed so that all 
+the :term:`diagonal leaf nodes` store *square* blocks. This way, the 
+:term:`HODLR` always captures the diagonal - which is typically the densest 
+region of the kind of matrix well represented by a :term:`HODLR` - using dense
+data. As a consequence, however, a rectangular :term:`HODLR` may be difficult 
+to represent, though there are currently 
+:ref:`no plans to do so<hodlr-always-square>`.
+
+
+.. _tree-iteration:
+
+HODLR tree traversal uses iteration
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+For all operations on :c:struct:`TreeHODLR`, it has been decided to iterate
+over the :term:`tree` rather than use recursion. This was mostly done because
+of concerns for how well recursion would be handled (i.e. potential to run
+into stack overflows etc.), but no rigorous tests on the topic were preformed
+and so is more of a preference choice.
+

docs/source/theory/hodlr_theory.rst

@@ -0,0 +1,187 @@
+What is HODLR?
+==============
+
+Hierarchical Off-Diagonal Low-Rank (HODLR) matrix is a matrix representation
+in which a matrix is recursively partitioned into quarters, with the 
+off-diagonal blocks stored as low-rank matrices. 
+
+Off-Diagonal
+------------
+
+In other words, given a square dense matrix :math:`D`, we construct a HODLR 
+matrix :math:`H` by first dividing :math:`D` into 4 blocks:
+
+.. image:: img/convert.svg
+   :alt: Visual representation of the conversion of a dense matrix into the 
+         HODLR format. Shows a dense matrix as a filled square being turned
+         into HODLR matrix, a square with the top left and bottom right
+         quarters filled.
+
+The two blocks on the diagonal are then stored as dense matrices while the two
+*off-diagonal* blocks are compressed and stored as *low-rank* matrices:
+
+.. math::
+
+   H=
+   \left[ {\begin{array}{cc}
+   {}^{0,0}D & {}^{0,1}U {}^{0,1}V^T \\
+   {}^{1,0}U {}^{1,0}V^T & {}^{1,1}D \\
+   \end{array} } \right]
+
+where :math:`{}^{i,i}D` is a dense block and :math:`{}^{i,j}U {}^{i,j}V^T` is 
+a low-rank block. Visually:
+
+.. image:: img/partition.svg
+
+
+.. _low-rank-explanation:
+   :alt: Visual representation of a HODLR matrix. Shows a square with the top
+         left and bottom right quarters filled, both of which are labelled as
+         dense. The top right and bottom left quarters are labelled as 
+         off-diagonal.
+
+
+Low-Rank
+--------
+
+The dense blocks are simple copies of the blocks of the dense matrix 
+:math:`D`, but the low-rank blocks are obtained by compressing the blocks of
+:math:`D` using `Singular Value Decomposition`_ (SVD). Using SVD, any real
+matrix can be decomposed into:
+
+.. math::
+
+   D = U \Sigma V^T
+
+
+where :math:`U` and :math:`V^T` are square orthogonal matrices and 
+:math:`\Sigma` is a rectangular diagonal matrix containing the singular 
+values. The singular values (:math:`\sigma_k = \Sigma_{k,k}`) are real
+non-negative numbers and, when obtained computationally, they usually come 
+sorted in descending order (:math:`\Sigma_{k,k} < \Sigma_{k+1,k+1}`). 
+Therefore, when performing an SVD of an off-diagonal (``i!=j``) block:
+
+.. math::
+
+   {}^{i,j}D = {}^{i,j}U {}^{i,j}\Sigma {}^{i,j}V^T
+
+where :math:`{}^{i,j}D` is a dense diagonal block of size ``m×n``, 
+:math:`{}^{i,j}U` is the ``m×m`` left singular matrix, :math:`{}^{i,j}\Sigma` 
+is the ``m×n`` matrix of singular values, and :math:`{}^{i,j}V^T` is the 
+``n×n`` right-singular matrix. Visually:
+
+.. image:: img/svd.svg
+   :alt: Diagram illustrating SVD using a filled square for D, U, and V 
+         transpose matrices and an empty square with a diagonal for sigma.
+
+
+Truncating zeroes
+^^^^^^^^^^^^^^^^^
+
+If :math:`D` is indeed structured correctly and suitable for conversion
+to HODLR, it will be the case that the singular values, of this off-diagonal
+block will decay rapidly (:math:`{}^{i,j}\Sigma_{r,r} \approx 0` for a 
+:math:`r << \min(m, n)`). In that case, the decomposition can be written as:
+
+.. math::
+
+   {}^{i,j}D = \sum_{k=0}^{k<r}{{}^{i,j}\sigma_k {}^{i,j}\mathbf{u}_k {}^{i,j}\mathbf{v^T}_k}
+
+where :math:`r` is the number of non-zero singular values, also called 
+the **rank** of the matrix. :math:`{}^{i,j}\mathbf{u}_k` is a vector 
+representing a column of the :math:`{}^{i,j}U` and 
+:math:`{}^{i,j}\mathbf{v^T}_k` is a vector representing a row of the 
+:math:`{}^{i,j}V^T` matrix:
+
+.. image:: img/svd2.svg
+   :alt: Diagram illustrating truncated SVD using the previous diagram, but
+         with the diagonal of the sigma square going only halfway from the
+         top right corner. A vertical line halfway through U and a horizontal
+         line halfway through V transpose also indicate that only half the
+         matrices will be kept.
+
+
+Truncating insignificant singular values
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Furthermore, it is expected that, due to this decay, the 
+singular values quickly become small enough in comparison to the first 
+singular value that they can be considered insignificant (
+:math:`{}^{i,j}\Sigma_{r',r'} << {}^{i,j}\Sigma_{0,0}` for :math:`r' < r`):
+
+.. image:: img/svd3.svg
+   :alt: Diagram illustrating an even more truncated SVD. Uses the previous
+         diagram, but the diagonal of sigma only extends a tenth of the way,
+         after which it appears as a very faint line. The lines through the
+         U and V transpose are now placed at a tenth of the matrices.
+
+Then, we can discard all the insignificant singular values and truncate the 
+:math:`{}^{i,j}U` and :math:`{}^{i,j}V^T` matrices, keeping only a small 
+number of columns/rows. This way, the off-diagonal block :math:`{}^{i,j}D` 
+is approximated as:
+
+.. math::
+
+   {}^{i,j}D \approx {}^{i,j}\hat{U} {}^{i,j}\hat{\Sigma} {}^{i,j}\hat{V}^T
+
+or
+
+.. math::
+
+   {}^{i,j}D \approx \sum_{k=0}^{k<r'}{{}^{i,j}\sigma_k {}^{i,j}\mathbf{u}_k {}^{i,j}\mathbf{v^T}_k}
+
+
+.. _u-scaling:
+
+In practice
+^^^^^^^^^^^
+
+In practice, to save memory, the :math:`{}^{i,j}U` matrix can be scaled by the
+singular values (i.e. :math:`{}^{i,j}U' = {}^{i,j}U {}^{i,j}\Sigma`) with
+only the result and the :math:`{}^{i,j}V^T` matrices stored:
+
+.. math::
+
+   {}^{i,j}D \approx {}^{i,j}\hat{U'} {}^{i,j}\hat{V}^T
+
+or
+
+.. math::
+
+   {}^{i,j}D \approx \sum_{k=0}^{k<r'}{ {}^{i,j}\mathbf{u'}_k {}^{i,j}\mathbf{v^T}_k}
+
+or visually:
+
+.. image:: img/svd4.svg
+   :alt: Diagram illustrating the fully truncated SVD. Uses a similar diagram
+         to the previous ones, but with the middle square representing sigma
+         completely gone, the U becoming a tall and narrow rectangle, and V 
+         becoming a wide and short rectangle.
+
+For more information, see the 
+:doc:`explanation from the code perspective<library/hodlr>`.
+
+
+Hierarchical
+------------
+
+Lastly, the above procedure can be repeated for the two diagonal blocks:
+
+.. image:: img/conversion2.svg
+   :alt: Diagram showing the conversion from a height 1 HODLR to a height 2
+         HODLR.
+
+splitting each block into quarters and compressing the two off-diagonal 
+sub-blocks. This can be repeated any number of times, yielding a HODLR matrix:
+
+.. math::
+
+   H = 
+   \left[ {\begin{array}{cc}
+   {}^{0,0}H & {}^{0,1}U {}^{0,1}V^T \\
+   {}^{1,0}U {}^{1,0}V^T & {}^{1,1}H \\
+   \end{array} } \right]
+
+
+.. _Singular Value Decomposition: https://en.wikipedia.org/wiki/Singular_value_decomposition
+