Problem setting or how to use h2tools

Main feature of h2tools is a special method of approximation of given matrix by \(\mathcal{H}^2\)-matrix directly from some of its entries. So, problem setting here is about constructing approximation and working with approximant. To build any hierarchical approximation we need hierarchical (mosaic) partitioning, which is usually based on block cluster tree.

Assume we have matrix \(A\), which corresponds to our initial problem. Since size of the initial problem can be incredibly large, we cannot store whole matrix \(A\) in memory and number of operations for a simple matrix-vector operation can be too big (linear in size of matrix \(A\)). So, we need a function, so-called block function, that will return submatrix of \(A\).

Since matrix \(A\) can be assumed as linear operator, we denote rows and columns of \(A\) as range of values and domain of corresponding linear operator (discretized initial problem). In accordance to N-body problems, we also name rows and columns of \(A\) as destinations and sources of corresponding interaction.

To make it short, there are following virtual stages of working with h2tools:

  1. Build discretization of initial problem, so that it has corresponding matrix \(A\) and block function.
  2. Build cluster tree for destination discrete elements (rows of \(A\)) and and for source discrete elements (columns of \(A\)). Resulting cluster trees are also called row and column cluster trees.
  3. Link row and column cluster trees into block cluster tree by acquiring near-field and far-field relations between row and column cluster trees.
  4. Build hierarchical approximation by passing block function and block cluster tree as arguments to approximation procedure.
  5. Use resulting approximant in any way you like.

Easiest way to understand following documentation is to read minimal_example in examples directory with ipython notebook.

1a. Building discretization

There is no special procedure in h2tools, which will build discretization. However, Python objects, containing discretizations (range of values and domain) or destinations and sources, should have special methods. These methods are used to build hierarchical partitioning.

Assume we have data object data (i.e. corresponding to destinations/range of values), consisting of some discrete elements or particles. Then, data should have method __len__, which returns number of discrete elements in data. If you plan to build hierarchical partitioning automatically by h2tools, data object should additionally have following methods:

  • divide: divides cluster of discrete elements into subclusters of discrete elements,
  • check_far: returns True if two clusters of discrete elements are geometrically far from each other,
  • compute_aux: returns auxiliary data for given cluster (i.e. bounding box).

Both sources and destinations (domain and range of values) should have corresponding data object, except for case of symmetric problem (one object is enough). To understand minimal requirements on data, see this example:

1b. Providing block function

Block function must have following definition:

def block_func(row_data, rows, col_data, cols):
    # row_data is a data, corresponding to destinations (range of values)
    # col_data is a data, corresponding to sources (domain)
    # rows and cols are arrays of indexes of rows and columns,
    # corresponding to required submatrix
    return submatrix

h2tools supports vector/matrix/tensor kernels, but to use such kernels block function must return submatrix of a special shape. If element_shape is a shape of kernel element, then entire submatrix must have shape (n_rows, element_shape[0], ..., element_shape[-1], n_columns). Here is simple example:

def block_func(row_data, rows, col_data, cols):
    # Function, that returns radius-vectors from sources (col_data) to
    #   destinations (row_data)
    shape = (rows.size, 3, cols.size)
    submatrix = np.zeros(shape)
    for i in range(rows.size):
        for j in range(cols.size):
            submatrix[i, :, j] = row_data[rows[i]]-col_data[cols[j]]
    return submatrix

2. Building cluster trees

Current implementation of cluster trees generation is built into procedure of linking cluster trees into block cluster tree. This stage requires only initialization of row and column cluster trees.

Assume row_data corresponds to destinations and col_data corresponds to sources. Then, following simple example shows cluster tree initialization:

from h2tools import ClusterTree
# row_tree_block_size and col_tree_block_size stand for maximum size
# of leaf nodes of row and column cluster trees
row_tree_block_size = 50
col_tree_block_size = 50
row_tree = ClusterTree(row_data, row_tree_block_size)
col_tree = ClusterTree(col_data, col_tree_block_size)

More information on ClusterTree class is here:

3. Linking block cluster tree

If cluster trees are already initialized, hierarchical partitioning can be easily computed with following command:

from h2tools import Problem
# func is a block function, returning submatrix of initial matrix
# symmetric shows if initial matrix is symmetric
# in symmetric case row_tree and col_tree must be the same Python object
# (not two different instances)
# verbose is a flag of additional output
problem = Problem(func, row_tree, col_tree, symmetric, verbose)

More information on Problem class is here:

4. Building approximation

After Problem object, containing source data, destination data, block function and block cluster tree, is ready, MCBH-approximation by \(\mathcal{H}^2\)-matrix can be computed by following command:

from h2tools.mcbh import mcbh
# tau is a parameter of relative spectral tolerance of approximation
# iters is a number of iterations of MCBH-method
# onfly is a boolean flag, showing if basis submatrices should be computed
# on fly or saved in memory (onfly=1 dramatically reduces memory
# consumption by approximant)
# verbose is a boolean flag, showing if additional information on
# approximation should go to output
matrix = mcbh(problem, tau, iters, onfly, verbose)
# matrix is an instance of H2matrix

More information on multicharge Barnes-Hut (MCBH) method is here:

5. Working with provided approximant

MCBH approximation procedure returns instance of H2matrix class. Full description of H2matrix class is here:

Gaining Performance boost

As was mentioned previously, approximation requires block function and block cluster tree. Block function is a numerical bottleneck of entire approximation problem and can be accelerated by means of Cython or Numba (becomes up to 100 times faster). Bottleneck of block cluster tree generation is check_far method of source and destination data objects. More information on accelerating tree generation is here: