🔢 mathematics

Sensor Placement for Tsunami Early Warning via Large-Scale Bayesian Optimal Experimental Design

This paper presents a scalable, multi-GPU Bayesian optimal experimental design framework that overcomes the computational challenges of hyperbolic PDEs by reformulating the problem as a dense matrix subset selection task, successfully optimizing a 175-sensor network for tsunami early warning in the Cascadia Subduction Zone with over one billion degrees of freedom.

Original authors: Sreeram Venkat, Stefan Henneking, Omar Ghattas

Published 2026-04-13

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Sreeram Venkat, Stefan Henneking, Omar Ghattas

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine the ocean floor off the coast of the Pacific Northwest is a giant, ticking time bomb. If it explodes (a massive earthquake), it could trigger a tsunami that hits the shore in minutes, leaving no time for warning.

To stop this, scientists want to build a giant "ear" underwater—a network of sensors that listens to the ocean floor the moment it moves. But here's the catch: Sensors are incredibly expensive. You can't just buy 1,000 of them and drop them everywhere. You have a limited budget, maybe enough for 175.

The big question is: Where do you put those 175 sensors so they give you the best possible warning?

This paper is about solving that puzzle using a super-smart, super-fast computer algorithm. Here is how they did it, explained simply:

1. The Problem: The "Needle in a Haystack"

Imagine you have a giant map with 600 possible spots to drop a sensor. You need to pick the best 175.

The Math Nightmare: If you tried to test every possible combination of 175 sensors out of 600, you would have more combinations than there are atoms in the universe. Even the world's fastest supercomputers would take thousands of years to check them all.
The Physics Problem: Tsunamis are governed by "hyperbolic" physics (like sound waves). Unlike heat spreading out slowly (which is easy to predict), these waves are chaotic and don't smooth out. This makes them very hard to model with standard math shortcuts.

2. The Solution: A "Data Space" Shortcut

The researchers realized they didn't need to look at the entire ocean floor (which has over 1 billion variables) to figure out where to put the sensors.

The Analogy: Imagine trying to figure out where a storm is by looking at the entire sky (1 billion clouds) versus just looking at the barometer readings from a few weather stations.
The Trick: They used a mathematical magic trick (called the Sherman-Morrison-Woodbury identity) to flip the problem. Instead of trying to solve the puzzle with 1 billion pieces, they shrank it down to a manageable puzzle with only about 100,000 pieces (the data from the sensors).
The Result: Suddenly, a problem that was impossible became solvable.

3. The Algorithm: The "Smart Greedy" Picker

Even with the smaller puzzle, checking every combination is too slow. So, they used a Greedy Algorithm.

How it works: Imagine you are building a team of 175 players.
1. You pick the single best player available.
2. Then, you look at the remaining players and pick the one who works best with the first player you chose.
3. You keep doing this, one by one, until you have 175.
Why it's good: While this doesn't guarantee the perfect theoretical team, math proves it will get you at least 63% of the way to perfection. In the real world, that's a massive improvement over random guessing.

4. The Speed: The "Conveyor Belt" Supercomputer

The biggest hurdle was that the computer had to read huge files from the hard drive while doing complex math. Usually, the computer waits for the hard drive, then does math, then waits again. This is like a chef waiting for the dishwasher to finish before chopping the next vegetable.

The authors built a Pipeline System:

The Analogy: Imagine a factory assembly line. While one robot is chopping vegetables (doing math), a second robot is already grabbing the next batch of vegetables from the fridge (reading data). A third robot is washing the dishes (I/O) in the background.
The Tech: They used Multi-GPU computers (like the ones at the world's fastest supercomputers, Perlmutter and Frontier). They made sure the computer never stopped working. While the graphics cards were crunching numbers, the hard drives were feeding them the next batch of data.
The Result: They achieved "near-perfect" speed. If they doubled the number of computers, they cut the time in half.

5. The Outcome: A Real-World Blueprint

They tested this on a "Digital Twin" of the Cascadia Subduction Zone (the real earthquake zone).

The Task: Pick the best 175 sensors out of 600 candidates.
The Time: It took 1.5 hours on 16 powerful GPUs.
The Comparison: A random guess would have been terrible. Their algorithm found a setup that drastically reduced the uncertainty of the tsunami forecast.

Why This Matters

This isn't just about math; it's about saving lives.

Before: We had to guess where to put expensive sensors, or we couldn't afford a network big enough to be accurate.
Now: We have a mathematical tool that tells us exactly where to spend our money to get the best possible warning system.

In a nutshell: The authors turned a math problem that would take a supercomputer 50,000 years to solve, into a 90-minute task. They did this by shrinking the problem, using a smart step-by-step selection process, and building a computer pipeline that never stops moving. This paves the way for a real, life-saving tsunami warning network in the Pacific Northwest.

1. Problem Statement

The paper addresses the critical challenge of optimizing sensor placement for Tsunami Early Warning Systems (TEWS), specifically for the Cascadia Subduction Zone (CSZ).

Context: Traditional warning systems rely on seismic data and simplified fault assumptions, which often fail for near-field events or complex rupture dynamics. Modern approaches use real-time offshore ocean bottom pressure and GNSS data to infer seafloor motion (the tsunami source) via Bayesian Inverse Problems.
The Challenge:
- Scale: The parameter space (spatiotemporal evolution of seafloor motion) is extreme, exceeding $10^9$ degrees of freedom.
- Computational Intractability: Standard Bayesian Optimal Experimental Design (OED) requires solving many inverse problems. For hyperbolic Partial Differential Equations (PDEs) governing tsunami propagation, the lack of spectral decay prevents the use of standard low-rank approximations.
- Cost: Offshore sensor deployment is expensive. Determining the optimal subset of sensors (e.g., selecting 175 out of 600 candidates) to maximize information gain is an NP-hard combinatorial problem.
- Previous Limitations: Directly solving this OED problem was considered impossible; previous estimates suggested it would take 50 years on hundreds of GPUs to find the posterior mean, let alone optimize the sensor network.

2. Methodology

The authors propose a scalable framework that transforms the intractable parameter-space problem into a tractable data-space problem, solved via a highly optimized greedy algorithm.

A. Mathematical Reformulation (Data Space)

Linear Time-Invariant (LTI) Systems: The tsunami propagation model is treated as an LTI system.
Sherman-Morrison-Woodbury Identity: Instead of working in the billion-dimensional parameter space, the authors reformulate the Bayesian inversion in the lower-dimensional data space.
Data-Space Hessian ( $K$ ): They construct a dense Hessian matrix $K$ of size $(N_d N_t) \times (N_d N_t)$ , where $N_d$ is the number of sensors and $N_t$ is the number of timesteps. This matrix can be explicitly formed using FFT-based algorithms, avoiding direct PDE solves during the optimization phase.
Objective Function: The problem is cast as D-optimal design, minimizing the log-determinant of the posterior covariance. This is equivalent to maximizing the Expected Information Gain (EIG). The objective function is shown to be submodular, allowing for greedy approximation with theoretical guarantees.

B. Algorithmic Innovation: Block Schur Complement Update

To solve the greedy sensor selection efficiently:

Naive Approach: Refactorizing the test matrix for every candidate sensor at every step has $O(k^3)$ complexity per candidate, leading to $O(|C|B^4)$ total complexity.
Proposed Approach: The authors develop a parallel block Schur complement update algorithm.
- It reuses the Cholesky factor from the previous iteration.
- It avoids redundant dense matrix factorizations.
- Complexity Reduction: Reduces per-candidate cost to $O(k^2)$ , lowering total complexity to $O(|C|B^3)$ . This is a massive reduction, making the problem tractable for large budgets ( $B$ ).

C. High-Performance Computing (HPC) Implementation

The algorithm is implemented using PyTorch and MPI for multi-GPU execution, featuring:

Pipelined I/O: A double-buffered architecture where CPU I/O (reading matrix blocks from disk) is fully overlapped with GPU computations.
Zero-Allocation Memory Management: Pre-allocating memory buffers to prevent fragmentation and dynamic allocation latency during the inner loop.
In-Place Operations: Strict in-place linear algebra to minimize GPU memory footprint, allowing for larger sensor budgets.
Precision Strategy: Uses single-precision (float32) for the greedy selection loop (sufficient for ranking candidates) while storing the master matrix in double-precision (float64) for stability.
Scalability: Designed to run on leadership-class supercomputers (Perlmutter and Frontier) with varying architectures (NVIDIA and AMD GPUs).

3. Key Contributions

First Extreme-Scale PDE-Constrained OED: The first framework to solve a Bayesian OED problem with a parameter space of $O(10^9)$ dimensions without using reduced-order models or surrogates.
Algorithmic Acceleration: Introduction of a Schur-complement-based greedy algorithm that eliminates redundant factorizations, significantly reducing memory and compute requirements.
Extreme-Scale Parallelization: A co-designed, multi-GPU implementation that achieves near-perfect weak and strong scaling across hundreds of GPUs, utilizing pipelined I/O to hide data transfer latency.
Real-World Application: Successful application to the Cascadia Subduction Zone digital twin, optimizing a 175-sensor network from 600 candidates to minimize uncertainty in seafloor motion inference.

4. Results

Performance Scaling:
- Tested on Perlmutter (NVIDIA A100) and Frontier (AMD MI250X).
- Achieved near-perfect weak and strong scaling over a 128x increase in GPU count (up to 512 GPUs / 1,024 GCDs).
- The Schur formulation was orders of magnitude faster than the naive refactorization approach and allowed for larger sensor budgets due to lower memory usage.
Application to Tsunami Warning:
- Setup: 600 candidate sensor locations; 512 GPUs used to precompute the adjoint PDE solutions (offline phase, ~500 hours).
- Optimization: The greedy selection of 175 sensors from 600 candidates was completed in 1.5 hours using only 16 A100 GPUs.
- Quality: The greedy-optimal configuration significantly outperformed 100 random sensor configurations in terms of the objective function value (information gain).
- Uncertainty Reduction: The selected sensors effectively reduced the pointwise standard deviation of the inferred seafloor displacement field across the domain.

5. Significance

Feasibility of Extreme-Scale Digital Twins: This work proves that high-fidelity, billion-parameter Bayesian OED is computationally feasible, removing the need for inaccurate model reductions.
Cost-Effective Infrastructure Planning: By providing a rigorous mathematical tool to optimize sensor placement, the framework helps justify the massive investment required for offshore sensor networks (like the proposed CSZ network), ensuring maximum predictive power for a constrained budget.
Generalizability: While demonstrated for tsunamis, the framework applies to any Bayesian inverse problem governed by LTI systems, including other wave propagation phenomena (acoustics, electromagnetics) where hyperbolic PDEs prevent low-rank approximations.
Operational Impact: Enables rapid "what-if" analysis for sensor network design, allowing stakeholders to evaluate different budget constraints, cost models, and noise characteristics in an agile manner.