New Algebraic Fast Algorithms for $N$-body Problems in Two and Three Dimensions

This paper introduces and evaluates two new algebraic multilevel hierarchical matrix algorithms, $\mathcal{H}^2_{*}$ and $(\mathcal{H}^2 + \mathcal{H})_{*}$, based on a weak admissibility condition and nested basis representations, which offer competitive memory and time efficiency for solving $N$-body problems in two and three dimensions compared to standard non-nested approaches.

The Gist

Imagine you are at a massive, crowded party with millions of people (let's call them "particles"). Everyone wants to talk to everyone else. In the world of physics and math, this is called an N-body problem.

If you tried to calculate exactly how much every single person influences every other person, you'd have to do a trillion calculations. That's like trying to read every book in the Library of Congress to find one specific sentence. It takes too long and uses too much memory. This is the "O(N²)" problem mentioned in the paper: as the party grows, the work grows explosively.

The authors of this paper, Ritesh Khan and Sivaram Ambikasaran, have built two new "super-fast" ways to handle this chaos. They call them algebraic fast algorithms. Here is the simple breakdown of what they did, using some everyday analogies.

The Problem: The "Too Many Neighbors" Issue

In these calculations, people (particles) are grouped into neighborhoods (clusters).

• Far-away neighbors: If you are in one neighborhood and someone is in a neighborhood far across the city, you don't need to know their exact face. You just need a general idea of their "group vibe." Mathematically, these interactions are "low-rank" (simple).
• Close neighbors: If you are next door, you need to know exactly who they are. These interactions are "high-rank" (complex).

Old methods had a rule: "If you are far away, we can simplify." But they were too strict. They said, "If you share a wall or a corner with a neighbor, you are too close to simplify." This forced the computer to do a lot of heavy lifting for people who were actually quite far apart in terms of influence.

The New Idea: The "Weak Admissibility" Rule

The authors introduced a new, smarter rule called Weak Admissibility.

• The Old Rule: "Only simplify if you are completely separated by a gap."
• The New Rule: "You can simplify even if you just share a corner (a vertex) with a neighbor."

Think of it like this: If you are in a house, and your neighbor's house shares a single corner with yours, you don't need to knock on their door to know they exist. You can just shout a general greeting. This tiny change allows the computer to treat many more interactions as "simple," saving a massive amount of time.

The Two New Algorithms

The paper presents two new ways to organize this party, both using a "Black Box" approach (meaning they work for any type of math problem without needing to know the specific physics rules beforehand).

1. The "Fully Nested" Algorithm (Efficient H2*)

The Analogy: Imagine a Russian Nesting Doll or a Family Tree.

• In the old way, every time you wanted to talk to a group, you had to introduce yourself to every single person in that group individually.
• In this new "Nested" way, you introduce yourself to the Head of the Family (the parent). The Head of the Family then passes the message down to their children, who pass it to their children.
• Why it's better: You don't need to carry a list of 1,000 names. You just carry the name of the Head of the Family, and they handle the rest. This saves huge amounts of memory and makes the calculation incredibly fast.
• The Trick: The authors realized that for "corner-sharing" neighbors, you can't just use the standard family tree method. So, they built a hybrid tree:
  • For people far away, they use a "Bottom-Up" tree (kids tell parents).
  • For people sharing a corner, they use a "Top-Down" tree (parents tell kids).
  • By mixing these two directions, they get the speed of the nesting dolls without the errors.

2. The "Semi-Nested" Algorithm (H2 + H)*

The Analogy: Imagine a Hybrid Car.

• This algorithm is a mix of two styles.
• For the "Far Away" groups, it uses the fancy, efficient Nesting Doll method (like the first algorithm).
• For the "Corner-Sharing" groups, it uses a simpler, older method (like a standard H-matrix).
• Why it's good: It's not as complex to build as the fully nested one, but it's still much faster than the old non-nested methods. It's a great "middle ground" that balances speed and setup time.

The Results: Why Should You Care?

The authors tested these algorithms on 2D (flat) and 3D (spatial) problems, simulating everything from heat flow to sound waves.

1. Speed: They can calculate the interactions of millions of particles in seconds, whereas the old methods would take hours or crash the computer.
2. Memory: They use much less computer memory. It's like packing a suitcase efficiently so you can fit a whole wardrobe in a backpack.
3. Versatility: Because they are "algebraic" (math-based) and not "analytic" (physics-rule-based), they work on any problem, even weird ones where the rules change from place to place.

The Bottom Line

The authors took a complex math problem that was like trying to count every grain of sand on a beach and turned it into a task that is like counting the buckets of sand.

They did this by realizing that sharing a corner is "far enough" to simplify, and by building a smarter, hybrid filing system (the nested algorithms) to organize the data. Their code is now open-source, meaning anyone can use these "super-speed" tools to solve difficult scientific problems faster than ever before.

Technical Summary

1. Problem Statement

The paper addresses the computational challenge of performing Matrix-Vector Products (MVP) for large, dense $N \times N$ kernel matrices arising from N-body problems in $d$ dimensions ($d=2, 3$).

• The Bottleneck: Direct evaluation of MVPs scales as $O(N^2)$, which is prohibitive for large $N$.
• Existing Solutions:
  • Fast Multipole Method (FMM) / H-matrices: Utilize "strong admissibility" (clusters must be well-separated by a distance greater than their diameter). While efficient, standard H-matrix algorithms often lack nested bases, leading to higher storage and computational costs.
  • Weak Admissibility (1D): In 1D, "weak admissibility" allows compressing interactions between adjacent (non-overlapping) intervals. However, a naive extension of this to higher dimensions ($d > 1$) fails because the rank of interactions between adjacent clusters (sharing edges/faces) grows with $N$, preventing quasi-linear complexity.
• The Gap: There is a need for purely algebraic (kernel-independent) hierarchical matrix algorithms that utilize a weak admissibility condition valid in higher dimensions (specifically including vertex-sharing clusters) while maintaining nested bases to achieve optimal $O(N)$ or quasi-linear $O(N \log N)$ complexity.

2. Methodology

The authors propose two new algebraic hierarchical matrix algorithms based on a specific weak admissibility condition introduced in their previous work [25].

Core Concept: Weak Admissibility in Higher Dimensions

Unlike standard admissibility (which excludes all neighbors), the proposed condition defines admissible clusters as:

1. Far-field clusters: Well-separated clusters (distance > 0).
2. Vertex-sharing clusters: Clusters that share a single vertex but do not share an edge or face.
   Crucially, clusters sharing an edge or face are treated as "near-field" and computed directly.

Proposed Algorithms

The paper introduces two algorithms utilizing Nested Cross Approximation (NCA) and Adaptive Cross Approximation (ACA):

A. Efficient $H2^*$ (Fully Nested Algorithm)

• Notation: $H2^*(b+t)$
• Strategy: Recognizing that vertex-sharing interactions have different rank properties than far-field interactions, the authors partition the interaction list ($IL^*$) into two sets:
  1. Far-field ($IL_{far}$): Compressed using Bottom-to-Top (B2T) NCA.
  2. Vertex-sharing ($IL_{ver}$): Compressed using Top-to-Bottom (T2B) NCA.
• Rationale:
  • B2T NCA works well for far-field interactions (standard H2 behavior).
  • T2B NCA is required for vertex-sharing interactions because the rank of these interactions grows polylogarithmically with $N$. A naive B2T approach on the entire list fails to capture the necessary pivots for vertex-sharing blocks, leading to poor approximation.
• Result: A fully nested basis representation that separates the two interaction types during initialization but combines them during the MVP.

B. Semi-Nested $(H2 + H)^*$

• Notation: $(H2 + H)^*$
• Strategy: A hybrid approach.
  • Far-field: Compressed using B2T NCA (nested).
  • Vertex-sharing: Compressed using standard ACA (non-nested).
• Result: This acts as a "cross" between H2 and H matrix algorithms. It sacrifices full nesting for the vertex-sharing part to reduce initialization complexity while still leveraging nesting for the far-field.

Implementation Details

• Kernel Independence: Both algorithms are purely algebraic, requiring only matrix entries, not analytic kernel expansions.
• Complexity:
  • Initialization: Scales quasi-linearly, $O(N \log^\alpha N)$.
  • MVP: Scales quasi-linearly, $O(N \log^\alpha N)$.
• Data Structures: Utilize uniform $2^d$ trees (Quad-trees in 2D, Oct-trees in 3D).

3. Key Contributions

1. Novel Admissibility Extension: Successfully extends the concept of weak admissibility to $d$-dimensions by identifying that only vertex-sharing interactions (not edge/face sharing) can be compressed with low rank.
2. Hybrid Pivot Selection: Demonstrates that a single pivot selection strategy (B2T or T2B) is insufficient for the mixed interaction list. The paper proposes a dual-strategy approach (B2T for far-field, T2B for vertex-sharing) to construct an efficient fully nested matrix ($H2^*$).
3. Semi-Nested Alternative: Introduces the $(H2 + H)^*$ algorithm, offering a trade-off between initialization speed and memory/MVP efficiency.
4. Comprehensive Benchmarking: Provides a rigorous comparison against standard algebraic H2 algorithms (using B2T and T2B NCA) and non-nested H/HODLR algorithms.
5. Open Source: The C++ implementation is released publicly, ensuring reproducibility.

4. Numerical Results

The algorithms were tested in 2D and 3D on various kernels (Laplacian, Matérn, Helmholtz, RBF) and applied to solve integral equations and RBF interpolation problems using GMRES.

• Memory Efficiency:
  • The proposed $H2^*(b+t)$ consistently uses less memory than the standard NCA-based H2 algorithms ($H2_{\sqrt{d}}$).
  • In 3D, for large $N$ (e.g., $N=1.7 \times 10^6$), the standard H2 algorithms ran out of memory, while $H2^*(b+t)$ and $(H2 + H)^*$ succeeded.
• Initialization Time:
  • $(H2 + H)^*$ has the fastest initialization time, often outperforming standard H2 algorithms.
  • $H2^*(b+t)$ has slightly higher initialization time than $(H2 + H)^*$ but is competitive with or faster than standard H2 algorithms in 3D due to smaller interaction lists.
• MVP and Solution Time:
  • $H2^*(b+t)$ achieves the fastest MVP and solution times among all nested algorithms in both 2D and 3D.
  • The semi-nested $(H2 + H)^*$ is competitive with standard H2 in 3D, often outperforming it in MVP time.
• Accuracy: All algorithms maintained relative errors within the specified tolerances ($10^{-8}$ to $10^{-12}$).

5. Significance

• Theoretical Advancement: The paper resolves the difficulty of applying weak admissibility in higher dimensions by proving that vertex-sharing clusters are compressible and by developing the specific algebraic machinery (mixed B2T/T2B NCA) to handle them efficiently.
• Practical Impact: The proposed algorithms offer a superior alternative to standard H2 matrices for kernel-independent applications. They provide:
  • Lower Memory Footprint: Critical for large-scale 3D simulations.
  • Faster Solves: Reduced MVP time accelerates iterative solvers like GMRES.
  • Robustness: Works effectively for non-translation invariant kernels and non-uniform particle distributions where analytic FMM methods struggle.
• Versatility: Being purely algebraic, these methods can be applied to any kernel function without deriving new translation operators, making them highly versatile for diverse scientific computing problems.

In conclusion, the authors successfully bridge the gap between weak admissibility and nested hierarchical matrices in higher dimensions, delivering algorithms that are theoretically sound, computationally efficient, and practically superior to existing standard H2 methods.