correct_topography: add batched-Cholesky GPU solver (opt-in via mintpy.topographicResidual.solver)

docs/gpu.md

@@ -1,14 +1,19 @@
-# Configure GPU acceleration for the network inversion #
+# Configure GPU acceleration #
 
-The `invert_network` step (in `ifgram_inversion.py`) ships an opt-in GPU solver that batches the per-pixel weighted least-squares inversion as normal-equations + Cholesky on a CUDA device via PyTorch. This is a partial GPU implementation: only `invert_network` is offloaded to the GPU; every other step in `smallbaselineApp.py` continues to run on the CPU. The solver is opt-in — the default `mintpy.networkInversion.solver = auto` resolves to `cpu`, so existing setups are unaffected.
+Two `smallbaselineApp.py` steps ship an opt-in GPU solver that batches their per-pixel inversion as normal-equations + Cholesky on a CUDA device via PyTorch:
+
+- **`invert_network`** (`ifgram_inversion.py`) — weighted least-squares network inversion. Toggle: `mintpy.networkInversion.solver = torch` / `--solver torch`.
+- **`correct_topography`** (`dem_error.py`) — pixelwise DEM-error fit. Toggle: `mintpy.topographicResidual.solver = torch` / `--solver torch`. Only the pixelwise-geometry branch (`pixelwiseGeometry = yes`, the production default) is GPU-dispatched; the mean-geometry branch is already pixel-batched on CPU and stays there.
+
+Every other step in `smallbaselineApp.py` continues to run on the CPU. Both solvers are opt-in — `mintpy.*.solver = auto` resolves to `cpu`, so existing setups are unaffected.
 
 The `torch` solver is orthogonal to Dask parallel processing (see [dask.md](./dask.md)): the former replaces the per-pixel CPU loop with a single batched Cholesky on one CUDA device, the latter distributes that same per-pixel loop across multiple worker processes. The two paths are not currently combined; pick one.
 
 ## 1. Setup ##
 
 See [installation.md](./installation.md) section 2.4 for installing the `[gpu]` extras with the matching CUDA wheel index. Selecting `solver = torch` on a host without a visible CUDA device is a hard error (no silent CPU fallback).
 
-## 2. Enable ##
+## 2. Enable on `invert_network` ##
 
 #### 2.1 via command line ####
 
@@ -49,7 +54,41 @@ ifgram_inversion.py inputs/ifgramStack.h5 -w no --solver torch
 
 The two outputs should agree to float32 round-off (RMS on the order of 1e-5).
 
-## 3. Behavior notes ##
+## 3. Enable on `correct_topography` ##
+
+The same `--solver torch` opt-in is exposed on `dem_error.py`, controlled by `mintpy.topographicResidual.solver`. Only the pixelwise-geometry branch (`pixelwiseGeometry = yes`) is GPU-dispatched.
+
+#### 3.1 via command line ####
+
+```bash
+dem_error.py timeseries_ERA5_ramp.h5 -g inputs/geometryRadar.h5 --solver torch
+dem_error.py timeseries_ERA5_ramp.h5 -g inputs/geometryRadar.h5 --solver torch --gpu-chunk-size 200000
+```
+
+#### 3.2 via template file ####
+
+```cfg
+mintpy.topographicResidual.solver       = torch  #[cpu / torch], auto for cpu
+mintpy.topographicResidual.gpuChunkSize = auto   #[int >= 0], auto for 0 (auto-size from free VRAM)
+```
+
+then:
+
+```bash
+smallbaselineApp.py smallbaselineApp.cfg
+```
+
+#### 3.3 Testing using example data ####
+
+```bash
+cd FernandinaSenDT128/mintpy
+dem_error.py timeseries_ERA5_ramp.h5 -g inputs/geometryRadar.h5 --solver cpu -o ts_demErr_cpu.h5
+dem_error.py timeseries_ERA5_ramp.h5 -g inputs/geometryRadar.h5 --solver torch -o ts_demErr_gpu.h5
+```
+
+CPU and GPU outputs should agree to float32 round-off (rms / |cpu|.max on the order of 1e-6 for `delta_z`, 1e-8 for `ts_cor`).
+
+## 4. Behavior notes ##
 
 + **VRAM auto-sizing.** `gpuChunkSize = 0` (auto) probes free VRAM at runtime and chooses a per-chunk pixel count with a fixed headroom factor. Set an explicit integer to override (e.g. for reproducible chunking across hosts with different VRAM).
 
@@ -59,10 +98,12 @@ The two outputs should agree to float32 round-off (RMS on the order of 1e-5).
 
 + **No silent CPU fallback.** Selecting `solver = torch` on a host without a visible CUDA device raises immediately rather than silently falling back to CPU; this keeps performance regressions visible.
 
-## 4. Performance ##
+## 5. Performance ##
 
 Benchmarks on RTX 5080 (Blackwell sm_120, CUDA 12.8, PyTorch 2.11) are tracked in the sibling [`mintpy-benchmark`](https://github.com/s-sasaki-earthsea-wizard/mintpy-benchmark) repository.
 
+### `invert_network`
+
 Tutorial-scale on FernandinaSenDT128 (270k pixels, 288 ifgs):
 
 + [report_torch.md](https://github.com/s-sasaki-earthsea-wizard/mintpy-benchmark/blob/main/reports/report_torch.md) — `cpu` vs `torch` end-to-end
@@ -73,3 +114,14 @@ Tutorial-scale on FernandinaSenDT128 (270k pixels, 288 ifgs):
 Large-scene on GalapagosSenDT128 (3.4M pixels, 475 kept ifgs; ~12.6× pixels and 1.65× ifgs vs Fernandina):
 
 + [report_large_scene.md](https://github.com/s-sasaki-earthsea-wizard/mintpy-benchmark/blob/main/reports/report_large_scene.md) — `solver=torch` reaches **36.4× step wall** / **44.4× internal** speedup on `invert_network` (cpu 6189 s → torch 170 s on RTX 5080 / SSD), confirming the speedup grows at scale; output equivalence preserved at float32 round-off (abs RMS max ~16 µm)
+
+### `correct_topography`
+
+Five-scene survey (Fernandina + Galapagos + three Zenodo Tier 1 scenes covering ISCE2 / GMTSAR / ARIA / ROI_PAC ingest pipelines, two wavelengths, and D ∈ [24, 333]):
+
++ [report_bench_survey.md](https://github.com/s-sasaki-earthsea-wizard/mintpy-benchmark/blob/main/reports/dem_error/report_bench_survey.md) — combined speedup-vs-K curve, numeric gate validation, axes coverage analysis
++ Per-scene: [report_fernandina.md](https://github.com/s-sasaki-earthsea-wizard/mintpy-benchmark/blob/main/reports/dem_error/report_fernandina.md), [report_galapagos.md](https://github.com/s-sasaki-earthsea-wizard/mintpy-benchmark/blob/main/reports/dem_error/report_galapagos.md), [report_sanfranbay.md](https://github.com/s-sasaki-earthsea-wizard/mintpy-benchmark/blob/main/reports/dem_error/report_sanfranbay.md), [report_sanfran_aria.md](https://github.com/s-sasaki-earthsea-wizard/mintpy-benchmark/blob/main/reports/dem_error/report_sanfran_aria.md), [report_kuju.md](https://github.com/s-sasaki-earthsea-wizard/mintpy-benchmark/blob/main/reports/dem_error/report_kuju.md)
+
+Headline: on local SSD, tutorial-scale scenes (~250–325k pixels) typically see **2–2.5×** speedup; production-scale (K = 3.4M Galapagos) reaches **6.15×**. CPU/GPU numeric agreement at `rms / |cpu|.max < 1e-5` on all 5 scenes. Note that the GPU path is **more I/O sensitive** than CPU because its compute portion is sub-second — networked storage (CIFS / NAS) can erase the speedup at small scenes (Fernandina drops from 2.56× on SSD to 0.97× on NAS).
+
+The closed-form speedup-vs-(K,D,P) model fitted against the survey lives in the wiki: [GPU Speedup Scaling Model](https://github.com/s-sasaki-earthsea-wizard/MintPy/wiki/GPU-Speedup-Scaling-Model). Short version: K is the GPU-parallelisable batch dimension; D and P are per-pixel internal dimensions that scale CPU and GPU walls similarly. Expect ≳ 4× speedup once K ≳ 1M.

src/mintpy/cli/dem_error.py

@@ -73,6 +73,21 @@ def create_parser(subparsers=None):
                              '1) output time-series file already exists, readable '+
                              'and newer than input interferograms file\n' +
                              '2) all configuration parameters are the same.')
+
+    # solver (pixelwise-geometry branch only; mean-geometry stays CPU)
+    solver = parser.add_argument_group('solver', 'per-pixel solver for the pixelwise-geometry branch')
+    solver.add_argument('--solver', dest='solver', default='cpu',
+                        choices={'cpu', 'torch'},
+                        help='per-pixel solver: cpu (scipy.linalg.lstsq, default) '
+                             'or torch (CUDA-batched normal-equation + Cholesky via '
+                             'PyTorch). torch requires the [gpu] extras and a visible '
+                             'CUDA device; absence is a hard error. Only the '
+                             'pixelwise-geometry branch is affected; the '
+                             'mean-geometry branch is already pixel-batched on CPU.')
+    solver.add_argument('--gpu-chunk-size', dest='gpuChunkSize', type=int, default=0,
+                        help='pixels per GPU chunk for --solver=torch '
+                             '(0=auto-size from free VRAM; default).')
+
     # computing
     parser = arg_utils.add_memory_argument(parser)
     parser = arg_utils.add_parallel_argument(parser)
@@ -145,6 +160,10 @@ def read_template2inps(template_file, inps):
                 iDict[key] = int(value)
             elif key in ['excludeDate','stepDate']:
                 iDict[key] = ptime.yyyymmdd(value.split(','))
+            elif key in ['solver']:
+                iDict[key] = value
+            elif key in ['gpuChunkSize']:
+                iDict[key] = int(value)
 
     # computing configurations
     dask_key_prefix = 'mintpy.compute.'

src/mintpy/defaults/smallbaselineApp.cfg

@@ -283,6 +283,16 @@ mintpy.topographicResidual.phaseVelocity     = auto  #[yes / no], auto for no -
 mintpy.topographicResidual.stepDate          = auto  #[20080529,20190704T1733 / no], auto for no, date of step jump
 mintpy.topographicResidual.excludeDate       = auto  #[20070321 / txtFile / no], auto for exclude_date.txt
 mintpy.topographicResidual.pixelwiseGeometry = auto  #[yes / no], auto for yes, use pixel-wise geometry info
+## Per-pixel solver for the pixelwise-geometry branch (GPU path is opt-in):
+## a. cpu   - scipy.linalg.lstsq, per-pixel (default, original behavior)
+## b. torch - batched normal-equation + Cholesky on CUDA via PyTorch.
+##            Requires the [gpu] extras (CUDA-enabled torch build) and a
+##            visible CUDA device; absence is a hard error (no silent CPU
+##            fallback). Only the pixelwise-geometry branch is affected;
+##            the mean-geometry branch is already pixel-batched on CPU.
+##            See docs/installation.md and docs/gpu.md.
+mintpy.topographicResidual.solver            = auto  #[cpu / torch], auto for cpu
+mintpy.topographicResidual.gpuChunkSize      = auto  #[int >= 0], auto for 0 (auto-size from free VRAM)
 
 
 ########## 11.1 residual_RMS (root mean squares for noise evaluation)

src/mintpy/defaults/smallbaselineApp_auto.cfg

@@ -122,6 +122,8 @@ mintpy.topographicResidual.phaseVelocity      = no
 mintpy.topographicResidual.stepDate           = no
 mintpy.topographicResidual.excludeDate        = exclude_date.txt
 mintpy.topographicResidual.pixelwiseGeometry  = yes
+mintpy.topographicResidual.solver             = cpu
+mintpy.topographicResidual.gpuChunkSize       = 0
 
 
 ########## residual_RMS

src/mintpy/dem_error.py

@@ -281,7 +281,8 @@ def estimate_dem_error(ts0, G0, tbase, date_flag=None, phase_velocity=False,
 
 
 def correct_dem_error_patch(G_defo, ts_file, geom_file=None, box=None,
-                            date_flag=None, phase_velocity=False):
+                            date_flag=None, phase_velocity=False,
+                            solver='cpu', chunk_size=None):
     """
     Correct one path of a time-series for DEM error.
 
@@ -291,6 +292,13 @@ def correct_dem_error_patch(G_defo, ts_file, geom_file=None, box=None,
                 box            - tuple of 4 int in (x0, y0, x1, y1) for the area of interest
                 date_flag      - 1D np.ndarray in bool in size of (num_date), dates used for the estimation
                 phase_velocity - bool, minimize the resdiual phase or phase velocity
+                solver         - str, 'cpu' (default, scipy.linalg.lstsq per pixel) or
+                                 'torch' (batched Cholesky on CUDA via mintpy.gpu.dem_error).
+                                 Only affects the pixelwise-geometry branch; the
+                                 mean-geometry branch is already pixel-batched on CPU.
+                chunk_size     - int or None. Pixels per GPU chunk for solver='torch'.
+                                 None / <=0 selects auto-sizing from free VRAM.
+                                 Ignored for solver='cpu'.
     Returns:    delta_z        - 2D np.ndarray in size of (num_row, num_col)
                 ts_cor         - 3D np.ndarray in size of (num_date, num_row, num_col)
                 ts_res         - 3D np.ndarray in size of (num_date, num_row, num_col)
@@ -386,6 +394,40 @@ def correct_dem_error_patch(G_defo, ts_file, geom_file=None, box=None,
         ts_cor[:, mask] = ts_cor_i
         ts_res[:, mask] = ts_res_i
 
+    elif solver != 'cpu':
+        print(f'estimating DEM error pixel-wisely via solver={solver!r} ...')
+        from mintpy.gpu.dem_error import (
+            estimate_dem_error_pixelwise_batch,
+            is_solver_available,
+        )
+        if not is_solver_available(solver):
+            raise RuntimeError(
+                f"solver={solver!r} requested but not available "
+                "(missing torch / CUDA)"
+            )
+
+        ts_valid = ts_data[:, mask]
+        if pbase.shape[1] == 1:
+            pbase_valid = pbase
+        else:
+            pbase_valid = pbase[:, mask]
+
+        delta_z_valid, ts_cor_valid, ts_res_valid = estimate_dem_error_pixelwise_batch(
+            ts_valid=ts_valid,
+            pbase_valid=pbase_valid,
+            range_dist_valid=range_dist[mask],
+            sin_inc_angle_valid=sin_inc_angle[mask],
+            G_defo=G_defo,
+            tbase=tbase,
+            date_flag=date_flag,
+            phase_velocity=phase_velocity,
+            solver=solver,
+            chunk_size=chunk_size,
+        )
+        delta_z[mask] = delta_z_valid
+        ts_cor[:, mask] = ts_cor_valid
+        ts_res[:, mask] = ts_res_valid
+
     else:
         print('estimating DEM error pixel-wisely ...')
         prog_bar = ptime.progressBar(maxValue=num_pixel2inv)
@@ -499,12 +541,17 @@ def correct_dem_error(inps):
     )
 
     # 3.2 prepare the input arguments for *_patch()
+    # chunk_size <=0 means "auto from free VRAM" inside the GPU dispatch;
+    # the patch function passes it through to mintpy.gpu.dem_error.
+    chunk_size = getattr(inps, 'gpuChunkSize', 0) or 0
     data_kwargs = {
         'G_defo'         : G_defo,
         'ts_file'        : inps.ts_file,
         'geom_file'      : inps.geom_file,
         'date_flag'      : date_flag,
         'phase_velocity' : inps.phaseVelocity,
+        'solver'         : getattr(inps, 'solver', 'cpu'),
+        'chunk_size'     : chunk_size if chunk_size > 0 else None,
     }
 
     # 3.3 invert / write block-by-block

src/mintpy/gpu/_common.py

@@ -0,0 +1,129 @@
+############################################################
+# Program is part of MintPy                                #
+# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi         #
+# Shared helpers for the mintpy.gpu batched solvers        #
+############################################################
+# Recommend import:
+#     from mintpy.gpu import _common as gpu_common
+
+
+"""Shared device/solver helpers for the mintpy.gpu batched solvers.
+
+Both ``mintpy.gpu.ifgram_inversion`` and ``mintpy.gpu.dem_error`` reduce
+a per-pixel weighted (or unweighted) least-squares problem to batched
+normal equations + cuSolver-batched Cholesky. This module factors out
+the pieces that do not depend on the specific solver:
+
+* PyTorch / CUDA availability probing (``HAS_TORCH``,
+  ``is_solver_available``, ``get_torch_device``).
+* VRAM-aware chunk sizing (``auto_chunk_size``).
+* The normal-equations + rank-deficient fallback core
+  (``solve_normal_equations_batched``).
+
+Callers are still responsible for assembling the per-pixel design
+matrix and right-hand-side (e.g. applying sqrt-weights, broadcasting a
+shared geometric column, ...) before passing them in.
+"""
+
+
+try:
+    import torch
+    HAS_TORCH = True
+except ImportError:
+    HAS_TORCH = False
+
+
+SUPPORTED_SOLVERS = ('cpu', 'torch')
+
+# default chunk size when caller does not provide one and VRAM probing fails
+DEFAULT_CHUNK_SIZE = 20000
+
+# safety factor applied to free VRAM when auto-sizing the chunk
+VRAM_SAFETY = 0.4
+
+
+def is_solver_available(solver):
+    """Return True if the named solver is importable and usable."""
+    if solver == 'cpu':
+        return True
+    if solver == 'torch':
+        return HAS_TORCH and torch.cuda.is_available()
+    return False
+
+
+def get_torch_device(solver):
+    """Return ``torch.device('cuda')`` for a GPU solver, else raise."""
+    if not HAS_TORCH:
+        raise ImportError(
+            f"solver='{solver}' requires PyTorch. "
+            "Install with `pip install -e \".[gpu]\" "
+            "--extra-index-url https://download.pytorch.org/whl/cu128 "
+            "--index-strategy unsafe-best-match`."
+        )
+    if not torch.cuda.is_available():
+        raise RuntimeError(
+            f"solver='{solver}' requires CUDA, but torch.cuda.is_available() is False."
+        )
+    return torch.device('cuda')
+
+
+def auto_chunk_size(num_rows, num_cols, dtype_bytes=4):
+    """Pick chunk size from free VRAM.
+
+    The dominant per-pixel allocation in both batched solvers is the
+    design matrix tensor of shape ``(n, num_rows, num_cols)``; ~3x of
+    that empirically covers temporaries (normal matrix, Cholesky
+    factor, residual). Returns ``DEFAULT_CHUNK_SIZE`` when CUDA is
+    unavailable.
+    """
+    if not (HAS_TORCH and torch.cuda.is_available()):
+        return DEFAULT_CHUNK_SIZE
+    free_b, _ = torch.cuda.mem_get_info()
+    per_pixel_bytes = 3 * num_rows * num_cols * dtype_bytes
+    n = int(VRAM_SAFETY * free_b / max(per_pixel_bytes, 1))
+    return max(1, n)
+
+
+def solve_normal_equations_batched(G_batch, y_batch, print_msg=True):
+    """Per-pixel least-squares via normal equations + batched Cholesky.
+
+    For each pixel k with system ``G_k @ x_k = y_k`` (shapes
+    ``(num_row, num_col)`` and ``(num_row,)``), solve the normal
+    equations
+        ``(G_k^T @ G_k) x_k = G_k^T @ y_k``
+    via cuSolver-batched Cholesky. Compared to a per-pixel QR path this
+    collapses ~num_col Householder iterations into a single batched
+    factorization.
+
+    Rank-deficient pixels are detected via ``cholesky_ex`` info codes;
+    their factor is replaced with identity and right-hand-side with
+    zeros so the downstream ``cholesky_solve`` produces an all-zero
+    solution for those pixels and never propagates NaN/Inf.
+
+    Args:
+        G_batch:   (n, num_row, num_col) per-pixel design matrix.
+                   Callers solving a weighted system must pre-apply
+                   sqrt-weights to ``G_batch`` and ``y_batch``.
+        y_batch:   (n, num_row) per-pixel observation vector.
+        print_msg: bool. Warn on rank-deficient pixels.
+
+    Returns:
+        (n, num_col) per-pixel solution.
+    """
+    G_T = G_batch.transpose(-1, -2)
+    N = G_T @ G_batch
+    r = G_T @ y_batch.unsqueeze(-1)
+
+    L, info = torch.linalg.cholesky_ex(N)
+    fail_mask = info != 0
+    if fail_mask.any():
+        n_fail = int(fail_mask.sum().item())
+        if print_msg:
+            print(f'WARNING: {n_fail} rank-deficient pixel(s) in chunk; '
+                  'setting solution to zero')
+        eye = torch.eye(N.shape[-1], device=L.device, dtype=L.dtype)
+        L = torch.where(fail_mask.view(-1, 1, 1), eye, L)
+        r = torch.where(fail_mask.view(-1, 1, 1), torch.zeros_like(r), r)
+
+    X = torch.cholesky_solve(r, L)
+    return X.squeeze(-1)