PaperlayerView exact revision on arXiv

Source-linked AI summary

Fast construction of hierarchical matrix representation from matrix-vector multiplication

Lin Lin, Jianfeng Lu, Lexing Ying

arXiv:1001.0149v2math.NAcs.NA

TL;DR

Constructing an H-matrix from a black-box matrix-vector routine is useful when an elliptic Green’s-function matrix is needed explicitly. The paper combines randomized low-rank approximation with top-down peeling; its numerical tests report efficient and accurate constructions, with uniform H1 outperforming H1 and H2 further improving results for large N.

  • Problem

    The task is to construct an unknown symmetric H-matrix from black-box matrix-vector products, motivated by obtaining elliptic Green’s-function matrices explicitly.

    Source: §1 · p. 1, §1.1 · p. 2

  • Method

    The algorithm uses randomized SVD to approximate low-rank blocks and peels them off level by level using matrix-vector multiplications.

    Source: §1.3 · p. 4, §4 · p. 30

  • Results

    Numerical tests report efficient and accurate implementations; uniform H1 is significantly more advantageous than H1 in time and memory, and H2 further improves results for large N.

    Source: §4 · p. 31

  • Takeaways & Limitations

    The resulting H-matrix can support faster matrix-vector multiplication and further compression for numerical homogenization or structured data.

    Source: §4 · p. 30

    Although the algorithm uses O(log n) matvecs asymptotically, examples can require thousands, so it may not suit many applications.

    Source: §4 · p. 31

Abstract

from arXiv · show

We develop a hierarchical matrix construction algorithm using matrix-vector multiplications, based on the randomized singular value decomposition of low-rank matrices. The algorithm uses $\mathcal{O}(\log n)$ applications of the matrix on structured random test vectors and $\mathcal{O}(n \log n)$ extra computational cost, where $n$ is the dimension of the unknown matrix. Numerical examples on constructing Green's functions for elliptic operators in two dimensions show efficiency and accuracy of the proposed algorithm.

1. Introduction

The paper asks how to construct a symmetric H-matrix from black-box matrix-vector multiplications, motivated especially by explicit representations of elliptic Green’s functions. It combines randomized low-rank approximation with a top-down peeling strategy and relates the approach to prior H-matrix and homogenization work.

1. Introduction

  • The target is an unknown symmetric matrix with low-rank or approximately low-rank off-diagonal blocks, accessible only through matrix-vector multiplication.The authors focus on symmetric matrices and say the method extends straightforwardly to the nonsymmetric case.Source: §1 · pp. 1–2

1.1. Motivation and applications

  • For elliptic equations, constructing the Green’s-function matrix explicitly can support data-dependent compression, while the construction can be parallelized.The paper notes that direct fast solvers may be preferable for only a few right-hand sides.Source: §1.1 · p. 2

1.2. Randomized singular value decomposition algorithm

  • Randomized SVD approximates numerically low-rank blocks using matrix-vector products; its stated cost is O(r+c) matvecs plus O(n(r+c)^2) extra steps.The method adopted from [7] forms products with random matrices and uses SVD to obtain low-rank factors.Source: §1.2 · pp. 3–4

1.3. Top-down construction of H-matrix

  • The construction peels off off-diagonal blocks level by level, combining random tests to reduce matrix-vector products and yielding logarithmically many applications overall.The introduction presents a one-dimensional symmetric example and notes its simplifying low-rank assumption may not hold for general higher-dimensional elliptic operators.Source: §1.3 · pp. 4, 7

1.4. Related works

  • The approach requires a fast matrix-vector multiplication routine, which may come from elliptic solvers or boundary-integral methods, rather than direct access to matrix entries.Unlike Martinsson’s method [29], it does not assume arbitrary entry access, which the authors note can be difficult for inverse stiffness matrices.Source: §1.4 · pp. 8–9
  • The authors also connect explicit Green’s-function construction to homogenization, contrasting O(log n) solves with an approach requiring d solves under the Cordes condition.They state their approach makes no analytic assumptions about the operator, while being more costly in solves.Source: §1.4 · p. 9

2. Algorithm

The algorithm section develops the construction in two dimensions, using top-down peeling and the H2 property to reduce complexity. It also outlines variants that construct H1 and uniform H1 representations.

2. Algorithm

  • The two-dimensional algorithm combines top-down peeling with the H2 property to further reduce complexity.The authors describe extension to three dimensions as straightforward.Source: §2 · p. 9
  • The section first introduces a more conservative two-dimensional domain partition intended to guarantee low-rank matrix blocks.The geometric setup is followed by hierarchical-matrix definitions and construction procedures.Source: §2 · p. 9
  • Variants for constructing H1 and uniform H1 representations are also described.These variants are presented after the H2 construction.Source: §2 · p. 10

2.1. Geometric setup and notations

The geometric setup partitions a periodic two-dimensional grid into dyadic boxes and defines neighbor and interaction lists at each level. The matrix is decomposed into levelwise interactions and a final adjacent-box term.

2.1. Geometric setup and notations

  • The periodic square domain is discretized on an N × N uniform grid and partitioned into dyadic levels, with four child boxes per non-leaf box.The final level can be chosen so each leaf contains a constant number of grid points as N grows.Source: §2.1 · p. 10
  • Each box’s neighbor list contains itself and adjacent boxes, while its interaction list contains non-neighbors inherited from its parent’s neighbor set.At level 3 the interaction list has 55 boxes; at later levels it has 27.Source: §2.1 · p. 10, Figure 1 · p. 11, Figure 2 · p. 11
  • The matrix is decomposed into interactions associated with levels 3 through L and a final term for adjacent boxes at level L.For non-neighbor boxes, the relevant interaction is represented at a level where their containing boxes belong to each other’s interaction lists.Source: §2.1 · p. 12

2.2. Hierarchical matrix

The section distinguishes H1, uniform H1, and H2 representations by how their off-diagonal low-rank factors are organized and reused. Their application costs differ, with H2 reducing the logarithmic factor associated with H1.

2.2.1. H1 matrices

  • An H1 matrix approximates interaction-list blocks with low-rank factors, giving matrix-vector application complexity O(N^2 log N).The stated approximation error is O(ε).Source: §2.2.1 · pp. 12–13

2.2.2. Uniform H1 matrix

  • A uniform H1 matrix reuses one orthogonal basis per box, improving the application prefactor and providing a stated speedup of roughly 27r(ε)/rU(ε).Its asymptotic application complexity remains O(N^2 log N).Source: §2.2.2 · p. 13

2.2.3. H2 matrices

  • An H2 matrix adds nested bases across parent-child boxes, giving application complexity O(N^2), one logarithmic factor below H1.Its application corresponds to the M2M, M2L, and L2L translations of the fast multipole method.Source: §2.2.3 · p. 14
  • The algorithm also requires applying a partial representation containing the levelwise terms from G(3) through G(L′), for some L′ ≤ L.This operation is handled by a variant of the H2 application algorithm.Source: §2.2.3 · p. 14

2.3. Peeling algorithm: outline and preparation

The algorithm constructs an H2 representation from black-box matrix-vector products through a top-down sequence of peeling, uniformization, and projection. Randomized SVD supplies low-rank bases, while projection updates the nested bases and transfer matrices.

2.3. Peeling algorithm: outline and preparation

  • At each level, peeling extracts a new H1 representation, uniformization makes it uniform, and projection incorporates it into the H2 representation built so far.The steps repeat from level 3 upward.Source: §2.3 · pp. 15–16
  • Uniformization finds a shared column space by applying SVD to weighted neighboring-block bases, with randomized SVD offered as the more effective option.The bases are weighted by their coupling matrices so directions associated with larger singular values matter more.Source: §2.3 · pp. 16–17
  • Uniformization and projection each cost O(N^2) per level across all boxes.The stated costs follow from the SVD work summed over the boxes on a level.Source: §2.3 · pp. 17–18
  • Projection combines a child’s new basis with its existing basis, updates the child basis by SVD, and refreshes transfer and coupling matrices.The existing basis may not span the new basis, so it must be updated as well.Source: §2.3 · p. 18

2.4. Peeling algorithm: details

The detailed construction applies peeling level by level, using grouped test vectors to recover off-diagonal blocks, then extracts the finest-level diagonal part. Grouping exploits the limited support of each residual product, and the full peeling procedure costs O(log N) matrix-vector products plus O(N^2 log N) extra steps.

2.4.1. Level

  • At level 3, randomized SVD approximates neighboring off-diagonal blocks from matvecs with localized random probes, followed by uniformization and projection.The construction uses r(ε) + c matvecs per box, or 64(r(ε) + c) total at that level.Source: §2.4.1 · p. 19

2.4.2. Level

  • At level 4, 64 groups of four boxes allow residual products to be computed together when their nonzero regions do not overlap.The residual is G − G(3); the level costs 64(r(ε) + c) matvecs and O(N^2) extra steps.Source: §2.4.2 · pp. 20–21

2.4.3. Level l

  • At each later level, the algorithm subtracts the previously constructed partial H2 matrix, groups boxes into 64 sets, and repeats peeling, uniformization, and projection.The per-level cost remains 64(r(ε) + c) matvecs and O(N^2) extra steps.Source: §2.4.3 · p. 22

2.4.4. Computation of D(L)

  • The finest-level diagonal part is recovered from the remaining matrix using 16 groups whose boxes have nonoverlapping neighbor lists.Because N/2^L is constant, this final step costs O(1) matvecs and O(N^2) extra steps.Source: §2.4.4 · p. 23

2.4. Peeling algorithm: details

  • O(log N) matvecs and O(N^2 log N) extra steps suffice for the full peeling algorithm; error accumulation across levels has no proof, though tested examples show no apparent problem.At each level the algorithm uses 64(r(ε) + c) matvecs and O(N^2) extra steps.Source: §2.4.4 · pp. 23–24

2.5. Peeling algorithm: variants

The paper describes H1-only and uniform H1 variants that omit some steps of the H2 construction. Although they have higher extra-step complexity, the uniform H1 version is reported to be faster and more memory-efficient than the H1 version.

2.5. Peeling algorithm: variants

  • The H1 variant performs only peeling, while the uniform H1 variant adds uniformization; each uses its respective partial representation in the next level’s peeling calculation.The variants construct H1 and uniform H1 representations, respectively.Source: §2.5 · p. 24
  • Both simplified variants use O(log N) matvecs but O(N^2 log^2 N) extra steps, and apply to matrices in H1 or uniform H1 that are not H2.Their extra-step cost is log N times that of the H2 version.Source: §2.5 · pp. 24–25
  • Despite uniformization, uniform H1 constructs much faster than H1, partly by reusing representations at future levels; it stores one UI matrix per box versus about 27, improving memory efficiency, especially in higher dimensions.Source: §2.5 · p. 25

3. Numerical results

The numerical tests compare H1, uniform H1, and H2 constructions for discretized elliptic operators, measuring cost and accuracy. With optimal level selection, uniform H1 improves on H1, while H2 brings further gains, especially for large N.

Cost and accuracy

  • The tests use the inverse of a discretized two-dimensional elliptic operator with a strongly random potential, and report both wall-clock time and matvec calls.The matvec routine is nested-dissection reordered block Gauss elimination, with O(N^2 log N) cost for an N × N grid.Source: §3 · p. 25
  • Matvec calls are equal across the three representations and grow logarithmically with N when the domain size at level L is fixed; the optimal L is shared across methods.The reported optimal levels are 4 for N = 32, 64, 128; 5 for N = 256; and 6 for N = 512.Source: §3 · p. 26, Table 1 · p. 27
  • Memory cost per degree of freedom generally decreases as L increases, because more of the original dense matrix is represented by low-rank approximations.Table 2 reports memory costs for H1, uniform H1, and H2 across grid sizes and levels.Source: §3 · p. 26, Table 2 · p. 28
  • Errors remain well controlled as N and L increase despite more aggressive compression; uniform H1 ranks are only slightly above H1, while H2 ranks are comparable or lower.The reported 2-norm errors are estimated by the power method, and Table 4 compares average ranks at N = 256.Source: §3 · pp. 26–27, Table 3 · p. 29, Table 4 · p. 29
  • With optimal L, uniform H1 is significantly more advantageous than H1, and H2 further improves over uniform H1, especially for large N.Figure 4 compares time and memory costs per degree of freedom for N = 32, 64, 128, 256, 512.Source: §3 · p. 26, Figure 4 · p. 28
  • For high-contrast divergence-form operators with V = 10^-3W or 10^-6W, the tested algorithms work well despite slower off-diagonal decay.The N = 64, L = 4 H2 tests report absolute and relative 2-norm errors; the small relative error for V = 10^-6W is attributed to the large 2-norm of H^-1.Source: §3 · p. 30, Table 5 · p. 30

4. Conclusions and future work

The proposed algorithm constructs hierarchical matrices from matrix-vector products, with numerical results showing efficiency and accuracy and advantages for uniform H1 and H2 representations. Its practical suitability depends on the application, and the current implementation and extensions have stated scope boundaries.

4. Conclusions and future work

  • The algorithm uses O(log n) matrix-vector multiplications and O(n log n) additional steps, computes off-diagonal blocks by top-down peeling, and is friendly to parallelization.The resulting H-matrix can support faster matrix-vector multiplication and numerical homogenization or upscaling.Source: §4 · p. 30
  • Uniform H1 is significantly more advantageous than H1 in time and memory, while H2 leads to further improvement for large N.The implementations are reported as efficient and accurate.Source: §4 · p. 31
  • Despite requiring O(log n) matvecs, the algorithms may use several thousand in the Section 3 example and may not suit many applications.The authors identify inverting the same very large system and homogenization without applicable analytic approaches as settings where the algorithm provides an effective alternative.Source: §4 · p. 31
  • The current implementation depends on a rectangular-domain geometric partition, while unstructured-grid and 3D extensions are described as possible.For 3D cases, the authors expect the constant to grow significantly.Source: §4 · p. 31
Loading 1001.0149v2…