HierMat package

block_cluster_tree module

block_cluster_tree.py: BlockClusterTree, build_block_cluster_tree(), recursion_build_block_cluster_tree()

class HierMat.block_cluster_tree.BlockClusterTree(left_clustertree, right_clustertree, sons=None, level=0, is_admissible=False, plot_info=None)[source]

Bases: object

Tree structure containing two ClusterTree objects.

Parameters:
  • left_clustertree (ClusterTree) – “left side” cluster tree
  • right_clustertree (ClusterTree) – “right side” cluster tree
  • sons (list(BlockClusterTree)) – sons if the current tree is not a leaf
  • level (int) – level of the current tree
  • is_admissible (bool) – whether the current tree is admissible or not
  • plot_info (tuple(int, int)) – coordinates of the current tree
depth(root_level=None)[source]

Depth of the tree.

The root is at level 0 if not specified otherwise.

Parameters:root_level (int) – internal use.
Returns:depth
Return type:int

Hint

recursive function

draw(axes, admissible_color='#1e26bc', inadmissible_color='#bc1d38')[source]

Draw a patch into given axes.

Parameters:
  • axes (matplotlib.pyplot.axes) – axes instance to draw in
  • admissible_color (str) – color for admissible patch (see matplotlib for color specs)
  • inadmissible_color (str) – color for inadmissible patch
plot_recursion(axes, admissible_color='#1e26bc', inadmissible_color='#bc1d38')[source]
shape()[source]

Return length of left and right cluster tree.

Returns:x and y dimension
Return type:tuple(int, int)
to_dot()[source]

Return a string for .dot representation.

to_list()[source]

List representation for export.

Return a list containing the object and a list with its sons.

Hint

recursive function

to_xml()[source]

Return a string for xml representation.

HierMat.block_cluster_tree.build_block_cluster_tree(left_cluster_tree, right_cluster_tree=None, start_level=0, admissible_function=<function admissible>)[source]

Factory for building a block cluster tree.

Parameters:
  • left_cluster_tree (ClusterTree) – “left side” cluster tree
  • right_cluster_tree (ClusterTree) – “right side” cluster tree
  • start_level (int) – counter that keeps track of the level
  • admissible_function (function that returns a bool) – function that determines whether two cluster trees are admissible or not. This is crucial for the structure of the block cluster tree. (See utils.admissible for an example)
Returns:

block cluster tree

Return type:

BlockClusterTree

HierMat.block_cluster_tree.recursion_build_block_cluster_tree(current_tree, admissible_function)[source]

Recursion to build_block_cluster_tree().

cluster module

cluster.py: Cluster object and iterator

class HierMat.cluster.Cluster(grid, indices=None)[source]

Bases: object

Handles operations on a Grid object by manipulating an index list.

Parameters:
  • grid (Grid) – grid to build cluster around
  • indices (list) – index list (optional)
diameter()[source]

Return diameter.

Return the maximal Euclidean distance between two points. For big Clusters this is costly.

Returns:diameter
Return type:float
dim()[source]

Return dimension.

Returns:dim of Grid
Return type:int
distance(other)[source]

Compute distance to other cluster

Return minimal Euclidean distance between points of self and other

Parameters:other – another instance of Cluster
Returns:distance
Return type:float
get_grid_item(item)[source]

Return item from grid.

Parameters:item (int) – index
get_grid_item_support(item)[source]

Return support of item from grid.

Parameters:item (tuple(float)) – point
get_grid_item_support_by_index(item)[source]

Return support of i-th item from the current index list.

This enables the user to get the support of the “last” item by calling get_grid_item_support_by_index(-1)

Parameters:item (int) – index
get_index(item)[source]

Return index at item.

Parameters:item (int) – index
Returns:item-th index
Return type:int
get_patch_coordinates()[source]

Return min and max out of indices.

Returns:min and max
Return type:tuple(int, int)
class HierMat.cluster.ClusterIterator(cluster)[source]

Bases: object

Iterator to Cluster object

next()[source]

cluster_tree module

cluster_tree.py: ClusterTree, build_cluster_tree(), recursion_build_cluster_tree()

class HierMat.cluster_tree.ClusterTree(content, sons=None, max_leaf_size=1, level=0)[source]

Bases: object

Tree structure built according to the Splitable in use

Splits are done until max_leaf_size is reached

depth(root_level=None)[source]

Get depth of the tree

Parameters:root_level (int) – internal use.
Returns:depth
Return type:int

Hint

recursive function

diameter()[source]

Return the diameter of content

Returns:diameter
Return type:float
distance(other)[source]

Return distance to other

Parameters:other (ClusterTree) – other cluster tree
Returns:distance
Return type:float
get_grid_item(item)[source]

Get grid item from content

Parameters:item (int) – index of item to get
get_grid_item_support(item)[source]

Return supports ofitem from grid

Parameters:item (tuple(float)) – point
get_grid_item_support_by_index(item)[source]

Return supports of i-th item from grid

Parameters:item (int) – index
get_index(item)[source]

Get index from content

Parameters:item (int) – index to get
Returns:index
Return type:int
get_patch_coordinates()[source]

Return min and max out of indices

Returns:min and max
Return type:tuple(int, int)
to_dot()[source]

Return a string for .dot representation

to_list()[source]

Give list representation for export

Return a list containing the object and a list with its sons

Hint

recursive function

to_xml()[source]

Return a string for xml representation

HierMat.cluster_tree.admissible(left_clustertree, right_clustertree)[source]

Default admissible condition for BlockClusterTree

True if the smaller diameter of the input is smaller or equal to the distance between the two ClusterTrees

Parameters:
  • left_clustertree (ClusterTree) – “Left-side” ClusterTree
  • right_clustertree (ClusterTree) – “Right-side” ClusterTree
Returns:

admissible

Return type:

bool

HierMat.cluster_tree.build_cluster_tree(splitable, max_leaf_size=1, start_level=0)[source]

Factory for building a cluster tree

Parameters:
  • splitable (Splitable) – strategy that decides the structure
  • max_leaf_size (int) – cluster size at which recursion stops
  • start_level (int) – counter to identify levels
Returns:

cluster tree

Return type:

ClusterTree

HierMat.cluster_tree.recursion_build_cluster_tree(current_tree)[source]

Recursion to build_cluster_tree()

cuboid module

cuboid.py: Cuboid and minimal_cuboid()

class HierMat.cuboid.Cuboid(low_corner, high_corner)[source]

Bases: object

Cuboid that is parallel to the axis in Cartesian coordinates.

Characterized by two diagonal corners.

  • Attributes:

    low_corner: numpy.array with minimal values in each dimension

    high_corner: numpy.array with maximal values in each dimension

diameter()[source]

Return the diameter of the cuboid.

Diameter is the Euclidean distance between low_corner and high_corner.

Returns:diameter
Return type:float
distance(other)[source]

Return distance to other Cuboid.

Distance is the minimal Euclidean distance between points in self and points in other.

Parameters:other (Cuboid) – other cuboid
Returns:distance
Return type:float
split(axis=None)[source]

Split the cuboid in half

If axis is specified, the cuboid is restructure along the given axis, else the maximal axis is chosen.

Parameters:axis (int) – axis along which to restructure (optional)
Returns:cuboid1, cuboid2
Return type:tuple(Cuboid, Cuboid)

grid module

grid.py: Grid object and iterator

class HierMat.grid.Grid(points, supports)[source]

Bases: object

Discretized grid characterized by points and supports

Parameters:
  • points (list[tuple(float)]) –

    list of coordinates

    The points of the discretized Grid.

  • supports (dict{point: support}) – dictionary mapping points with their supports
Raises:

ValueError – if points and supports have different length

dim()[source]

Dimension of the Grid

Note

this is the dimension of the first point

if you store points of different dimensions, this will be misleading
Returns:dimension
Return type:int
get_point(item)[source]

Return point at position item

Parameters:item (int) – index
Returns:point
get_support(item)[source]

Return support for item

Parameters:item (float or tuple(float)) – point
Returns:support
get_support_by_index(index)[source]

Return support for the i-th item

Parameters:index (int) – index
Returns:support
class HierMat.grid.GridIterator(grid)[source]

Bases: object

Iterator to Grid object

next()[source]

hmat module

hmat.py: HMat, build_hmatrix(), recursion_build_hmatrix(), StructureWarning

class HierMat.hmat.HMat(blocks=(), content=None, shape=(), parent_index=())[source]

Bases: object

Implement a hierarchical Matrix

Parameters:
  • blocks (list(HMat)) – list of HMat instances (children) (optional)
  • content (RMat or numpy.matrix) – the content if the matrix has no children (optional)
  • shape (tuple(int, int)) – the shape of the matrix (same as for numpy matrices)
  • parent_index (tuple(int, int)) – the index of this matrix with respect to its containing parent-matrix, zero based

Supported operations

  • + (formatted addition)
  • * (formatted multiplication)
  • -a (unary minus)
  • - (formatted subtraction)
  • == (equal)
  • != (not equal)
  • abs (frobenius norm)
  • inv (inversion)
block_structure()[source]

return the block structure of self

produce a dict with parent_index: block paris

Return type:dict
Returns:dict with index: HMatrix pairs
column_sequence()[source]

Return the sequence of column groups as in thm. 1.9.4. in [Eve66]

only for consistent matrices

Returns:list of column-sizes in first row
Return type:list(int)
Raises:StructureWarning if the matrix is not consistent
inv(make_copy=True)[source]

Invert the matrix according to chapter 7.5.1 in [Hac15]

Note

If make_copy is True, a deepcopy is made at the beginning. This can be very slow for large matrices.

Parameters:make_copy (bool) – make a copy of the matrix to not change the original (Default True)
Returns:approximation of inverse
Return type:HMat
Raises:numpy.LinAlgError if matrix is singular
is_consistent()[source]

check if the blocks are aligned, i.e. we have consistent rows and columns as in def. 1.9.5 in [Eve66]

Returns:True on consistency, False otherwise
Return type:bool
norm(order=None)[source]

Norm of the matrix

Parameters:order – order of the norm (see in numpy.linalg.norm())
Returns:norm
Return type:float
restructure(structure)[source]

Restructure self into blocks to match structure

Parameters:structure (dict(index: size)) – a dict as returned by HMat.block_structure
Returns:HMat that has the specified structure
Return type:HMat
Raises:NotImplementedError if self has blocks or if self contains unknown content
row_sequence()[source]

Return the sequence of row groups as in thm. 1.9.4. in [Eve66]

defined only for consistent matrices

Returns:list of row-sizes in first column
Return type:list(int)
Raises:StructureWarning if the matrix is not consistent
split_hmat(start_x, start_y, end_x, end_y)[source]

Fetch the block specified by indices from a numpy.matrix

Parameters:
  • start_x (int) – vertical start index
  • start_y (int) – horizontal start index
  • end_x (int) – vertical end index
  • end_y (int) – horizontal end index
Returns:

numpy.matrix with corresponding block

Return type:

numpy.matrix

split_rmat(start_x, start_y, end_x, end_y)[source]

Fetch the block specified by indices from a rmat

Example

r = RMat(numpy.matrix([[2], [2], [2]]),
         numpy.matrix([[3], [3], [3]]))
h = HMat(content=r, shape=(3, 3), parent_index=(0, 0))
res = h.split_rmat(0, 2, 1, 3)
res
<RMat with left_mat: matrix([[2]]), right_mat: matrix([[3]]) and max_rank: None> 
Parameters:
  • start_x (int) – vertical start index
  • start_y (int) – horizontal start index
  • end_x (int) – vertical end index
  • end_y (int) – horizontal end index
Returns:

rmat with corresponding block

Return type:

RMat

to_matrix()[source]

Full matrix representation

Returns:full matrix
Return type:numpy.matrix
transpose()[source]

Return transposed copy of self

Return type:RMat
zero(make_copy=True)[source]

Return a copy of self with zero entries

i.e. an HMat instance with the same structure, but all content is replaced by numpy.zeros

Note

If make_copy is True, a deepcopy is made at the beginning. This can be very slow for large matrices.

Parameters:make_copy (bool) – make a copy of the matrix to not change the original (Default True)
Returns:zero matrix
Return type:HMat
exception HierMat.hmat.StructureWarning[source]

Bases: exceptions.Warning

Special Warning used in this module to indicate bad block structure

HierMat.hmat.build_hmatrix(block_cluster_tree=None, generate_rmat_function=None, generate_full_matrix_function=None)[source]

Factory to build an hierarchical matrix

Takes a block cluster tree and generating functions for full matrices and low rank matrices respectively

The generating functions take a block cluster tree as input and return a RMat or numpy.matrix respectively

Parameters:
  • block_cluster_tree (BlockClusterTree) – block cluster tree giving the structure
  • generate_rmat_function – function taking an admissible block cluster tree and returning a rank-k matrix
  • generate_full_matrix_function – function taking an inadmissible block cluster tree and returning a numpy.matrix
Returns:

hmatrix

Return type:

HMat

HierMat.hmat.recursion_build_hmatrix(current_hmat, block_cluster_tree, generate_rmat, generate_full_mat)[source]

Recursion to build_hmatrix()

rmat module

rmat.py: RMat

class HierMat.rmat.RMat(left_mat, right_mat=None, max_rank=None)[source]

Bases: object

Rank-k matrix

Implementation of the low rank matrix as described in [HB2015]

form_add(other, rank=None)[source]

Formatted addition of self and other, i.e. addition and reduction to rank:

(self + other).reduce(rank)

If rank is omitted, reduction to min(rank(self), rank(other))

Parameters:
  • other (RMat) – other rank-k matrix
  • rank (int) – rank after reduction
Returns:

reduced result

Return type:

RMat

static from_matrix(matrix, max_rank)[source]

Convert a full matrix to the rank k format

Parameters:
  • matrix (numpy.matrix) – Matrix to convert to RMat
  • max_rank (int) – maximum rank
Returns:

Best approximation of matrix in RMat format

Return type:

RMat

norm(order=None)[source]

Norm of the matrix

Parameters:order – order of the norm (see in numpy.linalg.norm())
Returns:norm
Return type:float
reduce(new_k)[source]

Perform a reduced QR decomposition to rank new_k

Parameters:new_k (int) – rank to reduce to
Returns:reduced matrix
Return type:RMat
split(start_x, start_y, end_x, end_y)[source]

Return the block of self that matches the supplied indices

Parameters:
  • start_x (int) – x start index
  • start_y (int) – y start index
  • end_x (int) – x end index
  • end_y (int) – y end index
Returns:

rmat that represents the corresponding block

Return type:

RMat

to_matrix()[source]

Full matrix representation:

left_mat * right_mat.T
Returns:full matrix
Return type:numpy.matrix
transpose()[source]

Return transposed copy of self

Return type:RMat

splitable module

splitable.py: Splitable (Interface) and iterator and different strategies

  • RegularCuboid

    Starting from a minimal cuboid that get’s split in half on every level, points are distributed according to their geometric location

  • MinimalCuboid

    Same as RegularCuboid, but after every split, the cuboids are reduced to minimal size

  • Balanced

    Distributes the points evenly on every level, depends solely on the order of the initial list

class HierMat.splitable.Balanced(cluster)[source]

Bases: HierMat.splitable.Splitable

Balanced strategy, where a cluster is always restructure in half only according to its size, without regarding the geometry

Gives a binary tree

diameter()[source]

Return the diameter

Diameter of the contained cluster

Returns:diameter
Return type:float
distance(other)[source]

Distance to other balanced

Parameters:other (Balanced) – other balanced
Returns:distance
Return type:float
get_grid_item(item)[source]

Get grid item from cluster

Parameters:item (int) – index of item to get
get_grid_item_support(item)[source]

Return supports of item from grid

Parameters:item (tuple(float)) – point
get_grid_item_support_by_index(item)[source]

Return supports of item from grid

Parameters:item (int) – index
get_index(item)[source]

Get index from cluster

Parameters:item (int) – index to get
Returns:index
Return type:int
get_patch_coordinates()[source]

Return min and max out of indices

Returns:min and max
Return type:tuple(int, int)
split()[source]

Split in two

distribute the first (smaller) half to the left and the rest to the right

return itself if it has only one point

Return type:list(Balanced)
class HierMat.splitable.MinimalCuboid(cluster)[source]

Bases: HierMat.splitable.RegularCuboid

Method of minimal cuboids

Split is implemented by splitting the surrounding cuboid in half along longest axis and distributing indices according to the cuboid they belong to. Afterwards all resulting cuboids are shrunk to minimal.

Gives a binary tree

split()[source]

Split the cuboid and distribute items in cluster according to the cuboid they belong to. Reduce every restructure to minimal size

Returns:list of minimal cuboids of smaller size
Return type:list(MinimalCuboid)
class HierMat.splitable.RegularCuboid(cluster, cuboid=None)[source]

Bases: HierMat.splitable.Splitable

Method of regular cuboids

If no cuboid is given, a minimal cuboid is built around initial list of indices. Split is then implemented by splitting the surrounding cuboid in half along longest axis and distributing indices according to the cuboid they belong to.

Gives a binary tree

diameter()[source]

Return the diameter of the cuboid

Diameter of the surrounding cuboid

Returns:diameter
Return type:float
distance(other)[source]

Return the distance between own cuboid and other cuboid

Parameters:other (RegularCuboid) – other regular cuboid
Returns:distance
Return type:float
split()[source]

Split the cuboid and distribute items in cluster according to the cuboid they belong to

Returns:list of smaller regular cuboids
Return type:list(RegularCuboid)
class HierMat.splitable.Splitable[source]

Bases: object

Interface to the different strategies that can be used to restructure a cluster in two or more.

Methods that need to be implemented by subclasses

  • __eq__: check for equality against other Splitable
  • split: restructure the object in two or more, return new instances
  • diameter: return the diameter of the object
  • distance: return the distance to other object
cluster = None
diameter()[source]
distance(other)[source]
get_grid_item(item)[source]

Get grid item from cluster

Parameters:item (int) – index of item to get
get_grid_item_support(item)[source]

Return supports of item from grid

Parameters:item (tuple(float)) – point
get_grid_item_support_by_index(item)[source]

Return supports of item from grid

Parameters:item (int) – index
get_index(item)[source]

Get index from cluster

Parameters:item (int) – index to get
Returns:index
Return type:int
get_patch_coordinates()[source]

Return min and max out of cluster indices

Returns:min and max
Return type:tuple(int, int)
split()[source]
class HierMat.splitable.SplitableIterator(obj)[source]

Bases: object

Iterator to the Splitable implementations.

next()[source]
HierMat.splitable.minimal_cuboid(cluster)[source]

Build minimal cuboid

Build minimal cuboid around cluster that is parallel to the axis in Cartesian coordinates

Parameters:cluster (Cluster) – cluster to build cuboid around
Returns:minimal cuboid
Return type:Cuboid

utils module

utils.py: Utilities for the HMatrix module

HierMat.utils.block_cluster_tree_plot(obj, filename=None, ticks=False, face_color='#133f52', admissible_color='#76f7a8', inadmissible_color='#ff234b')[source]

Plot the block cluster tree

Parameters:
  • obj (BlockClusterTree) – block cluster tree to plot
  • filename (str) – filename to save the plot. if omitted, the plot will be displayed
  • ticks (bool) – show ticks in the plot
  • face_color – background color (see matplotlib for color specs)
  • admissible_color (str) – color for admissible patch
  • inadmissible_color (str) – color for inadmissible patch

Note

depends on matplotlib.pyplot

HierMat.utils.divisor_generator(n)[source]

Return divisors of n

Parameters:n (int) – integer to find divisors of
Returns:divisors
Return type:list[int]

Warning

This is a generator! To get a list with all divisors call:

list(divisor_generator(n))

Note

found at StackOverflow on 2017.03.08

HierMat.utils.export(obj, form='xml', out_file='out')[source]

Export obj in specified format.

Parameters:
Raises:

NotImplementedError – if form is not supported

Note

implemented so far:

  • xml
  • dot
  • bin
HierMat.utils.load(filename)[source]

Load a ClusterTree or BlockClusterTree from file

Parameters:filename (String) – file to import
Returns:object
Return type:BlockClusterTree or ClusterTree

Note

Depends on pickle

HierMat.utils.plot(obj, filename=None, **kwargs)[source]

plot an object

Parameters:
  • obj (BlockClusterTree) – object to plot
  • filename (str) – filename to save the plot to (if omitted, the plot will be displayed)
  • kwargs – optional arguments to specific plot commands see the respective documentations

Examples

Examples to the HierMat package

In this module are example scripts to demonstrate the usage of the package

model_1d module

This is an implementation of the one-dimensional model example described in [thesis]. It shows a typical use-case of hierarchical matrices:

Integral Equation

\[\int_0^1 \]

rac{1}{2pi}log |x-y| u(y) dy = g

for \(x \in [0,1]\).

After Galerkin discretization, we end up with a linear system

\[\mathbf{A}^{\text{Gal}} \cdot \alpha = \mathbf{g}\]

where

\[\mathbf{A}^{\text{Gal}}_{t,\tau}:= \int_{t}\int_{\tau} \log |x-y| dy dx\]

To determine admissible blocks we start by building the geometric objects:

midpoints = [((i + 0.5)/n,) for i in xrange(n)]
intervals = {p: (p[0] - 0.5/n, p[0] + 0.5/n) for p in midpoints}
grid = HierMat.Grid(points=midpoints, supports=intervals)
cluster = HierMat.Cluster(grid=grid)
unit_cuboid = HierMat.Cuboid([0], [1])
strategy = HierMat.RegularCuboid(cluster=cluster, cuboid=unit_cuboid)
cluster_tree = HierMat.build_cluster_tree(splitable=strategy, max_leaf_size=n_min)
block_cluster_tree = HierMat.build_block_cluster_tree(left_cluster_tree=cluster_tree,
                                                  right_cluster_tree=cluster_tree,
                                                  admissible_function=HierMat.admissible
                                                  )

With the structure established, we can produce the hierarchical matrix:

hmat = HierMat.build_hmatrix(block_cluster_tree=block_cluster_tree,
                         generate_rmat_function=lambda bct: galerkin_1d_rank_k(bct, max_rank),
                         generate_full_matrix_function=galerkin_1d_full
                         )
HierMat.examples.model_1d.galerkin_1d_full(block_cluster_tree)[source]

Exact calculation of the integral

\[A_{i,j}=A_{\tau,t}^{Gal}=\int_t\int_\tau\log\Vert x-y\Vert \;dydx\]
Parameters:block_cluster_tree (HierMat.BlockClusterTree) – inadmissible block cluster tree
Returns:matrix with same shape as block_cluster_tree.shape()
Return type:numpy.matrix
HierMat.examples.model_1d.galerkin_1d_rank_k(block_cluster_tree, max_rank)[source]

Low-rank approximation of the kernel

\[ \begin{align}\begin{aligned}R: &= A \cdot B^T\\A_{\tau, k}: &= \int_\tau \log \Vert x-y_k\Vert dx\\B_{\tau, k}: &= \int_\tau \mathcal{L}_k(y) dy\end{aligned}\end{align} \]
Parameters:
  • block_cluster_tree (HierMat.BlockClusterTree) – admissible block cluster tree
  • max_rank (int) – separation rank
Returns:

HierMat.examples.model_1d.get_chebyshev_interpol_points(points, lower=0, upper=1)[source]

Get Chebyshev interpolation points on interval \((a, b)\).

\[\frac{1}{2} (a + b) + \frac{1}{2} (b - a) \cos\left(\frac{2k+1}{2n}\pi\right)\]

for \(k=0, ..., \text{points} - 1\)

Parameters:
  • points (int) – number of points
  • lower (float) – \(a\), (default 0)
  • upper (float) – \(b\), (default 0)
Returns:

\(\frac{1}{2} (a + b) + \frac{1}{2} (b - a) \cos\left(\frac{2k+1}{2n}\pi\right)\)

Return type:

list(floats)

HierMat.examples.model_1d.ker(x, y)[source]

Kernel to integrate

\[ker(x, y):= \log\left(\Vert x - y \Vert\right)\]
Parameters:
  • x (float) – real number
  • y (float) – real number
Returns:

\(\log\left(\Vert x - y \Vert\right)\)

Return type:

float

HierMat.examples.model_1d.kerlog(x)[source]

kerlog function as in [thesis].

\[kerlog(x):= x^2 \left( \log(\Vert x \Vert) - \frac{1}{2} \right)\]
Parameters:x (float) – real number
Returns:\(x^2 ( \log(\Vert x \Vert) - \frac{1}{2} )\)
Return type:float
HierMat.examples.model_1d.model_1d(n=8, max_rank=2, n_min=2)[source]

This is an implementation of the one-dimensional model example described in [thesis].

Parameters:
  • n (int) – number of discretization points
  • max_rank (int) – max rank of the low-rank approximation
  • n_min (int) – minimal leaf size for cluster trees
Returns:

error