From 27f848ac251646bd8b49a8dd219a8e345bcad445 Mon Sep 17 00:00:00 2001 From: s-sasaki-earthsea-wizard Date: Wed, 6 May 2026 18:44:40 +0900 Subject: [PATCH 1/3] ifg_inv: add opt-in torch GPU solver for invert_network Add a CUDA-accelerated path for the per-pixel weighted least-squares inversion in `ifgram_inversion.py`, batched as normal-equations + Cholesky on a single CUDA device via PyTorch. The solver is opt-in and the default (`mintpy.networkInversion.solver = auto`) resolves to `cpu`, so existing setups are unaffected and the CPU path is byte-for-byte unchanged. Surface - cfg keys: `mintpy.networkInversion.solver = cpu|torch` (default `auto`), `mintpy.networkInversion.gpuChunkSize = ` (default 0 = auto-size). - CLI flags: `--solver {cpu,torch}` and `--gpu-chunk-size N` on `ifgram_inversion.py`. - New module `src/mintpy/ifgram_inversion_gpu.py` holds the torch path; `ifgram_inversion.py` dispatches to it only when `solver=torch` is explicitly requested. Behavior - VRAM auto-sizing probes free GPU memory and chooses a per-chunk pixel count with a fixed headroom factor; `gpuChunkSize > 0` overrides. - Rank-deficient pixels are detected via `torch.linalg.cholesky_ex` info codes and zeroed so NaN/Inf cannot propagate downstream. - Per-pixel NaN observations are handled by zeroing the corresponding row weight, which is mathematically equivalent to dropping that row from the WLS system. - Selecting `solver=torch` on a host without a visible CUDA device raises immediately rather than silently falling back to CPU, keeping any performance regression visible. Packaging - Adds `[gpu]` extras in `pyproject.toml`, sourced from `requirements-gpu.txt`. The PyTorch CUDA wheels live on a separate index; `installation.md` documents the install command in a follow-up commit. --- pyproject.toml | 3 + requirements-gpu.txt | 5 + src/mintpy/cli/ifgram_inversion.py | 14 +- src/mintpy/defaults/smallbaselineApp.cfg | 9 + src/mintpy/defaults/smallbaselineApp_auto.cfg | 2 + src/mintpy/ifgram_inversion.py | 42 ++- src/mintpy/ifgram_inversion_gpu.py | 293 ++++++++++++++++++ tests/test_ifgram_inversion_gpu.py | 227 ++++++++++++++ 8 files changed, 592 insertions(+), 3 deletions(-) create mode 100644 requirements-gpu.txt create mode 100644 src/mintpy/ifgram_inversion_gpu.py create mode 100644 tests/test_ifgram_inversion_gpu.py diff --git a/pyproject.toml b/pyproject.toml index e0cacae14..9561e4ccd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -114,6 +114,9 @@ readme = { file = ["docs/README.md"], content-type = "text/markdown" } [tool.setuptools.dynamic.optional-dependencies.test] file = ["tests/requirements.txt"] +[tool.setuptools.dynamic.optional-dependencies.gpu] +file = ["requirements-gpu.txt"] + [tool.setuptools.packages.find] where = ["src"] diff --git a/requirements-gpu.txt b/requirements-gpu.txt new file mode 100644 index 000000000..7d757f416 --- /dev/null +++ b/requirements-gpu.txt @@ -0,0 +1,5 @@ +# Optional GPU acceleration deps for the `[gpu]` extras. +# Install with the PyTorch CUDA wheel index, e.g.: +# pip install -e ".[gpu]" --extra-index-url https://download.pytorch.org/whl/cu128 +cupy-cuda12x +torch diff --git a/src/mintpy/cli/ifgram_inversion.py b/src/mintpy/cli/ifgram_inversion.py index f72fdc650..56561058c 100755 --- a/src/mintpy/cli/ifgram_inversion.py +++ b/src/mintpy/cli/ifgram_inversion.py @@ -86,6 +86,16 @@ def create_parser(subparsers=None): solver.add_argument('--min-norm-phase', dest='minNormVelocity', action='store_false', help=('Enable inversion with minimum-norm deformation phase,' ' instead of the default minimum-norm deformation velocity.')) + solver.add_argument('--solver', dest='solver', default='cpu', + choices={'cpu', 'torch'}, + help='WLS 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. ' + 'See docs/installation.md.') + 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).') #solver.add_argument('--norm', dest='residualNorm', default='L2', choices=['L1', 'L2'], # help='Optimization method, L1 or L2 norm. (default: %(default)s).') @@ -234,8 +244,10 @@ def read_template2inps(template_file, inps): elif value: if key in ['maskThreshold', 'minRedundancy']: iDict[key] = float(value) - elif key in ['residualNorm', 'waterMaskFile']: + elif key in ['residualNorm', 'waterMaskFile', 'solver']: iDict[key] = value + elif key in ['gpuChunkSize']: + iDict[key] = int(value) # computing configurations dask_key_prefix = 'mintpy.compute.' diff --git a/src/mintpy/defaults/smallbaselineApp.cfg b/src/mintpy/defaults/smallbaselineApp.cfg index 59d2756f3..02b4b65da 100644 --- a/src/mintpy/defaults/smallbaselineApp.cfg +++ b/src/mintpy/defaults/smallbaselineApp.cfg @@ -175,6 +175,15 @@ mintpy.unwrapError.bridgePtsRadius = auto #[1-inf], auto for 50, half size of t mintpy.networkInversion.weightFunc = auto #[var / fim / coh / no], auto for var mintpy.networkInversion.waterMaskFile = auto #[filename / no], auto for waterMask.h5 or no [if not found] mintpy.networkInversion.minNormVelocity = auto #[yes / no], auto for yes, min-norm deformation velocity / phase +## WLS solver for the per-pixel network inversion (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). The default 'auto' resolves to 'cpu', so existing +## setups are unaffected. See docs/installation.md. +mintpy.networkInversion.solver = auto #[cpu / torch], auto for cpu +mintpy.networkInversion.gpuChunkSize = auto #[int >= 0], auto for 0 (auto-size from free VRAM) ## mask options for unwrapPhase of each interferogram before inversion (recommend if weightFunct=no): ## a. coherence - mask out pixels with spatial coherence < maskThreshold diff --git a/src/mintpy/defaults/smallbaselineApp_auto.cfg b/src/mintpy/defaults/smallbaselineApp_auto.cfg index 6763d77bc..654b9a226 100644 --- a/src/mintpy/defaults/smallbaselineApp_auto.cfg +++ b/src/mintpy/defaults/smallbaselineApp_auto.cfg @@ -70,6 +70,8 @@ mintpy.unwrapError.bridgePtsRadius = 50 mintpy.networkInversion.weightFunc = var mintpy.networkInversion.waterMaskFile = waterMask.h5 mintpy.networkInversion.minNormVelocity = yes +mintpy.networkInversion.solver = cpu +mintpy.networkInversion.gpuChunkSize = 0 ## mask mintpy.networkInversion.maskDataset = no diff --git a/src/mintpy/ifgram_inversion.py b/src/mintpy/ifgram_inversion.py index 16c8e1dc2..52f66afea 100644 --- a/src/mintpy/ifgram_inversion.py +++ b/src/mintpy/ifgram_inversion.py @@ -595,7 +595,8 @@ def get_design_matrix4std(stack_obj): def run_ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_name='unwrapPhase', weight_func='var', water_mask_file=None, min_norm_velocity=True, - mask_ds_name=None, mask_threshold=0.4, min_redundancy=1.0, calc_cov=False): + mask_ds_name=None, mask_threshold=0.4, min_redundancy=1.0, calc_cov=False, + solver='cpu', gpu_chunk_size=0): """Invert one patch of an ifgram stack into timeseries. Parameters: ifgram_file - str, interferograms stack HDF5 file, e.g. ./inputs/ifgramStack.h5 @@ -610,6 +611,15 @@ def run_ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_nam mask_threshold - float, min coherence of pixels if mask_dataset_name='coherence' min_redundancy - float, the min number of ifgrams for every acquisition. calc_cov - bool, calculate the time series covariance matrix. + solver - str, WLS solver: 'cpu' (default, + scipy.linalg.lstsq per pixel) or 'torch' + (CUDA-batched normal-equation + Cholesky via + PyTorch). The 'torch' solver requires the + [gpu] extras and a visible CUDA device; + absence is a hard error (no silent CPU + fallback). + gpu_chunk_size - int, pixels per GPU chunk for solver='torch'. + 0 (default) auto-sizes from free VRAM. Returns: ts - 3D array in size of (num_date, num_row, num_col) ts_cov - 4D array in size of (num_date, num_date, num_row, num_col) or None inv_quality - 2D array in size of (num_row, num_col) @@ -800,8 +810,34 @@ def run_ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_nam 'inv_quality_name' : inv_quality_name, } + # 2.x GPU batched path: handles weighted and unweighted in one call. + # Per-pixel NaN observations are masked via zero-weights inside the kernel, + # which is mathematically equivalent to dropping them from the LS system + # for the full-rank case. Rank-deficient pixels (rare on real SBAS networks) + # are not handled here; if encountered, NaN/Inf will propagate downstream. + if solver != 'cpu': + from mintpy.ifgram_inversion_gpu import estimate_timeseries_batch + print(f'estimating time-series via {solver} solver (batched, GPU)') + ts_sub, q_sub, n_sub = estimate_timeseries_batch( + A=A, B=B, + y=stack_obs[:, idx_pixel2inv], + weight_sqrt=(weight_sqrt[:, idx_pixel2inv] + if weight_sqrt is not None else None), + tbase_diff=tbase_diff, + min_norm_velocity=min_norm_velocity, + rcond=1e-5, + min_redundancy=min_redundancy, + inv_quality_name=inv_quality_name, + chunk_size=gpu_chunk_size, + solver=solver, + ) + ts[:, idx_pixel2inv] = ts_sub + inv_quality[idx_pixel2inv] = q_sub + num_inv_obs[idx_pixel2inv] = n_sub + del mask + # 2.2 un-weighted inversion (classic SBAS) - if weight_sqrt is None: + elif weight_sqrt is None: msg = f'estimating time-series for pixels with valid {obs_ds_name} in' # a. split mask into mask_all/part_net @@ -1089,6 +1125,8 @@ def run_ifgram_inversion(inps): "mask_threshold" : inps.maskThreshold, "min_redundancy" : inps.minRedundancy, "calc_cov" : inps.calcCov, + "solver" : getattr(inps, 'solver', 'cpu'), + "gpu_chunk_size" : int(getattr(inps, 'gpuChunkSize', 0)), } # 3.3 invert / write block-by-block diff --git a/src/mintpy/ifgram_inversion_gpu.py b/src/mintpy/ifgram_inversion_gpu.py new file mode 100644 index 000000000..3da7cf80f --- /dev/null +++ b/src/mintpy/ifgram_inversion_gpu.py @@ -0,0 +1,293 @@ +############################################################ +# Program is part of MintPy # +# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi # +# GPU-accelerated network inversion # +############################################################ +# Recommend import: +# from mintpy import ifgram_inversion_gpu as ifginv_gpu + + +"""GPU-batched solver for the SBAS network inversion. + +This module provides ``estimate_timeseries_batch``, an opt-in replacement for +the per-pixel CPU loop in ``run_ifgram_inversion_patch``. Pixels are solved +in batches on CUDA via the normal equations + cuSolver-batched Cholesky: +this collapses the per-pixel Householder iterations of a QR path into ~5 +kernel launches per chunk on the FernandinaSenDT128 fixture. Rank-deficient +pixels are detected via ``cholesky_ex`` info codes and zeroed so NaN/Inf +never propagate. + +Per-pixel NaN observations are masked by zeroing the corresponding row +weights, which is mathematically equivalent to dropping them from the WLS +system. The CPU code path is unchanged when ``solver='cpu'`` and remains +the numerical reference. +""" + + +import numpy as np + +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 WLS 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): + 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_pair, num_unknown, dtype_bytes=4): + """Pick chunk size from free VRAM. + + The dominant per-pixel allocations during one chunk are: + Gw (n, num_pair, num_unknown) ~ num_pair * num_unknown * dtype_bytes + N, L (n, num_unknown, num_unknown) ~ 1/num_pair of the above + residual / temp ~ num_pair * dtype_bytes + Empirically, ~3x of Gw covers the rest. Use VRAM_SAFETY as the + headroom factor. + """ + 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_pair * num_unknown * dtype_bytes + n = int(VRAM_SAFETY * free_b / max(per_pixel_bytes, 1)) + return max(1, n) + + +def _solve_cholesky(G_dev, w_dev, y_dev): + """Per-pixel weighted least-squares via normal equations + batched Cholesky. + + For each pixel k with weighted system Gw_k @ x_k = yw_k where + Gw_k = diag(w_k) @ G and yw_k = w_k * y_k, solve the normal equations + (Gw_k^T @ Gw_k) x_k = Gw_k^T @ yw_k + via cuSolver-batched Cholesky. Compared to a gels (QR) path this + collapses ~num_unknown Householder iterations per pixel into a single + batched factorization, reducing per-chunk kernel launches from + ~num_pixel * num_unknown to ~5. + + 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_dev: (num_pair, num_unknown) design matrix on GPU. + w_dev: (num_pair, n) per-pixel sqrt-weights with NaN rows zeroed. + y_dev: (num_pair, n) observations with NaN rows zeroed. + + Returns: + (n, num_unknown) per-pixel solution. + """ + Gw = w_dev.t().unsqueeze(-1) * G_dev.unsqueeze(0) + yw = (w_dev * y_dev).t().unsqueeze(-1) + Gw_T = Gw.transpose(-1, -2) + N = Gw_T @ Gw + r = Gw_T @ yw + + L, info = torch.linalg.cholesky_ex(N) + fail_mask = info != 0 + if fail_mask.any(): + n_fail = int(fail_mask.sum().item()) + 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) + + +def estimate_timeseries_batch( + A, B, y, tbase_diff, + weight_sqrt=None, + min_norm_velocity=True, + rcond=1e-5, + min_redundancy=1.0, + inv_quality_name='temporalCoherence', + chunk_size=None, + solver='torch', + print_msg=True, +): + """Batch GPU least-squares solver for the SBAS network inversion. + + Solves, in batch over pixels k: + (G * w_k) X_k = (y_k * w_k) if weight_sqrt is not None (WLS) + G X_k = y_k otherwise (OLS) + where G = B if min_norm_velocity else A. + + Per-pixel NaN observations are handled by zeroing the corresponding row + weight, which is equivalent to dropping the row from the (W)LS system. + + Args: + A: np.ndarray (num_pair, num_date-1). Design matrix, + phase formulation (used when min_norm_velocity=False). + B: np.ndarray (num_pair, num_date-1). Design matrix, + velocity formulation (used when min_norm_velocity=True). + y: np.ndarray (num_pair, num_pixel). Observations, + NaN-tolerant. + tbase_diff: np.ndarray (num_date-1, 1). Differential temporal + baseline in years. + weight_sqrt: np.ndarray (num_pair, num_pixel) or None. Square + root of per-(ifgram, pixel) weight (WLS), or None + for OLS. + min_norm_velocity: bool. Solve for velocity (True) or phase (False). + rcond: float. Unused on CUDA: the Cholesky solver does + not consume an rcond cutoff. Kept for API parity + with the CPU path. + min_redundancy: float. Network-level redundancy check; if the design + matrix has any column with fewer non-zeros than + this, return zeros for all pixels. + inv_quality_name: str. 'temporalCoherence' | 'residual' | 'no'. + chunk_size: int or None. Pixels per GPU chunk. None => auto. + solver: str. 'torch' (only one implemented). + print_msg: bool. + + Returns: + ts: np.ndarray (num_date, num_pixel) float32. + inv_quality: np.ndarray (num_pixel,) float32. + num_inv_obs: np.ndarray (num_pixel,) int16. + """ + if solver != 'torch': + raise ValueError( + f"unsupported solver={solver!r}; choose from {SUPPORTED_SOLVERS}" + ) + device = _get_torch_device(solver) + + G = B if min_norm_velocity else A + num_pair, num_unknown = G.shape + num_pixel = y.shape[1] + num_date = num_unknown + 1 + + ts = np.zeros((num_date, num_pixel), dtype=np.float32) + inv_quality = np.zeros(num_pixel, dtype=np.float32) + num_inv_obs = np.zeros(num_pixel, dtype=np.int16) + + # network-level redundancy check (matches estimate_timeseries L162) + if np.min(np.sum(A != 0., axis=0)) < min_redundancy: + if print_msg: + print(f'network redundancy < {min_redundancy}; skip inversion') + return ts, inv_quality, num_inv_obs + + # decide chunk size + if chunk_size is None or chunk_size <= 0: + chunk_size = _auto_chunk_size(num_pair, num_unknown) + if print_msg: + free_gib = torch.cuda.mem_get_info()[0] / 2**30 + print(f'GPU auto chunk_size = {chunk_size} pixels ' + f'(free VRAM {free_gib:.1f} GiB)') + else: + chunk_size = int(chunk_size) + + num_chunk = (num_pixel + chunk_size - 1) // chunk_size + if print_msg: + mode = 'WLS' if weight_sqrt is not None else 'OLS' + print(f'estimating time-series via {solver} batched {mode} ' + f'in {num_chunk} chunk(s) of up to {chunk_size} pixels ...') + + # move design matrix and tbase to GPU once (re-used across chunks) + G_dev = torch.as_tensor(G, dtype=torch.float32, device=device) + tbase_dev = torch.as_tensor(np.asarray(tbase_diff).flatten(), + dtype=torch.float32, device=device) + + use_wls = weight_sqrt is not None + + for ci in range(num_chunk): + c0 = ci * chunk_size + c1 = min(c0 + chunk_size, num_pixel) + n = c1 - c0 + + # prepare per-chunk inputs (host side) + y_chunk = np.asarray(y[:, c0:c1], dtype=np.float32) + nan_mask = np.isnan(y_chunk) # (num_pair, n) + y_chunk = np.where(nan_mask, 0.0, y_chunk) + + if use_wls: + w_chunk = np.asarray(weight_sqrt[:, c0:c1], dtype=np.float32) + w_chunk = np.where(nan_mask, 0.0, w_chunk) + else: + w_chunk = (~nan_mask).astype(np.float32) + + # to GPU + y_dev = torch.as_tensor(y_chunk, device=device) # (num_pair, n) + w_dev = torch.as_tensor(w_chunk, device=device) # (num_pair, n) + valid_dev = torch.as_tensor(~nan_mask, device=device) # (num_pair, n) + + # solve per-pixel WLS via batched Cholesky on the normal equations + X_batch = _solve_cholesky(G_dev, w_dev, y_dev) # (n, num_unknown) + + # inversion-quality: |sum_i exp(j * e_i)| / N (phase coherence) + # N is the per-pixel count of valid (non-NaN) ifgrams, matching the + # CPU reference path which drops NaN rows via skip_invalid_obs before + # entering calc_inv_quality. + if inv_quality_name == 'temporalCoherence': + y_pred = G_dev @ X_batch.t() # (num_pair, n) + e_dev = (y_dev - y_pred) * valid_dev # zero-out NaN rows + cos_sum = e_dev.cos().sum(dim=0) + sin_sum = e_dev.sin().sum(dim=0) + # cos(0) = 1 contributes from masked rows; subtract their count + n_invalid = (~valid_dev).sum(dim=0).to(cos_sum.dtype) + cos_sum = cos_sum - n_invalid + n_valid = valid_dev.sum(dim=0).to(cos_sum.dtype).clamp(min=1.0) + tcoh = torch.sqrt(cos_sum * cos_sum + sin_sum * sin_sum) / n_valid + inv_quality[c0:c1] = tcoh.detach().cpu().numpy().astype(np.float32) + + elif inv_quality_name == 'residual': + y_pred = G_dev @ X_batch.t() + e_dev = (y_dev - y_pred) * valid_dev + inv_quality[c0:c1] = ( + e_dev.pow(2).sum(dim=0).sqrt().detach().cpu().numpy().astype(np.float32) + ) + # else 'no': leave zeros + + # assemble timeseries + if min_norm_velocity: + # X is per-interval velocity; ts[i+1] = ts[i] + v_i * dt_i + ts_diff = X_batch * tbase_dev.unsqueeze(0) # (n, num_unknown) + ts_chunk = ts_diff.cumsum(dim=1).t() # (num_unknown, n) + else: + # X is per-date phase + ts_chunk = X_batch.t() # (num_unknown, n) + ts[1:, c0:c1] = ts_chunk.detach().cpu().numpy().astype(np.float32) + + num_inv_obs[c0:c1] = (~nan_mask).sum(axis=0).astype(np.int16) + + # free per-chunk tensors before next iteration + del X_batch, y_dev, w_dev, valid_dev + + if print_msg: + chunk_step = max(1, num_chunk // 5) + if (ci + 1) % chunk_step == 0 or ci == num_chunk - 1: + print(f'chunk {ci + 1} / {num_chunk}') + + return ts, inv_quality, num_inv_obs diff --git a/tests/test_ifgram_inversion_gpu.py b/tests/test_ifgram_inversion_gpu.py new file mode 100644 index 000000000..9b8a45ff7 --- /dev/null +++ b/tests/test_ifgram_inversion_gpu.py @@ -0,0 +1,227 @@ +"""Numerical-equivalence tests for src/mintpy/ifgram_inversion_gpu.py. + +Compare the GPU-batched solver against the per-pixel CPU reference +(scipy.linalg.lstsq via mintpy.ifgram_inversion.estimate_timeseries) on +synthetic SBAS-like networks. Tests are skipped automatically when +PyTorch / CUDA are unavailable. +""" + +import numpy as np +import pytest + +torch = pytest.importorskip("torch") + +from mintpy.ifgram_inversion import estimate_timeseries +from mintpy.ifgram_inversion_gpu import estimate_timeseries_batch + +requires_cuda = pytest.mark.skipif( + not torch.cuda.is_available(), + reason="CUDA-capable GPU required for ifgram_inversion_gpu tests", +) + + +def make_redundant_network(num_date, num_pair, *, max_span=4, seed=0): + """Synthetic SBAS-like network with multiple ifgrams per adjacent date + interval. + + All adjacent pairs (i, i+1) are forced first to guarantee minimal + connectivity, then short-baseline pairs (i, j) with j - i <= max_span + are sampled until ``num_pair`` is reached. This mirrors a well-connected + Sentinel-1 SBAS network and stays full-rank under partial NaN masking + of individual interferograms. + """ + rng = np.random.default_rng(seed) + candidate = [(i, j) for i in range(num_date) + for j in range(i + 1, min(i + max_span + 1, num_date))] + if num_pair > len(candidate): + raise ValueError( + f"requested num_pair={num_pair} exceeds {len(candidate)} candidates " + f"for num_date={num_date}, max_span={max_span}" + ) + forced = [(i, i + 1) for i in range(num_date - 1)] + forced_set = set(forced) + optional = [p for p in candidate if p not in forced_set] + rng.shuffle(optional) + selected = sorted(forced + optional[: num_pair - len(forced)]) + + A = np.zeros((num_pair, num_date - 1), dtype=np.float32) + for k, (i, j) in enumerate(selected): + A[k, i:j] = 1.0 + tbase = np.cumsum(rng.uniform(0.05, 0.2, size=num_date - 1).astype(np.float32)) + tbase = np.concatenate([[0.0], tbase]) + tbase_diff = np.diff(tbase).reshape(-1, 1).astype(np.float32) + B = A * tbase_diff.T + return A, B, tbase_diff + + +def synthesize_observations(A, B, num_pixel, *, nan_frac=0.0, weighted=True, seed=0): + rng = np.random.default_rng(seed) + num_pair, num_unknown = A.shape + X_true = rng.normal(0, 1, size=(num_unknown, num_pixel)).astype(np.float32) + y = (B @ X_true) + rng.normal(0, 0.05, size=(num_pair, num_pixel)).astype(np.float32) + weight_sqrt = (rng.uniform(0.5, 2.0, size=(num_pair, num_pixel)).astype(np.float32) + if weighted else None) + if nan_frac > 0: + nan_mask = rng.random(y.shape) < nan_frac + y[nan_mask] = np.nan + return y, weight_sqrt + + +def cpu_reference(A, B, y, weight_sqrt, tbase_diff, *, min_norm_velocity=True): + """Per-pixel CPU reference via mintpy.ifgram_inversion.estimate_timeseries.""" + num_pair, num_unknown = A.shape + num_pixel = y.shape[1] + num_date = num_unknown + 1 + ts = np.zeros((num_date, num_pixel), dtype=np.float32) + tcoh = np.zeros(num_pixel, dtype=np.float32) + nobs = np.zeros(num_pixel, dtype=np.int16) + for k in range(num_pixel): + wk = None if weight_sqrt is None else weight_sqrt[:, k] + tsi, qi, ni = estimate_timeseries( + A=A, B=B, y=y[:, k], tbase_diff=tbase_diff, + weight_sqrt=wk, + min_norm_velocity=min_norm_velocity, + rcond=1e-5, min_redundancy=1.0, + inv_quality_name='temporalCoherence', + print_msg=False, + ) + ts[:, k] = tsi.flatten() + tcoh[k] = qi + nobs[k] = ni + return ts, tcoh, nobs + + +def assert_equivalent(cpu, gpu, *, ts_rel_tol, tcoh_abs_tol): + ts_cpu, tcoh_cpu, nobs_cpu = cpu + ts_gpu, tcoh_gpu, nobs_gpu = gpu + + assert np.all(np.isfinite(ts_gpu)), 'GPU returned non-finite timeseries' + assert np.all(np.isfinite(tcoh_gpu)), 'GPU returned non-finite temporalCoh' + + ts_scale = max(float(np.abs(ts_cpu).max()), 1e-6) + ts_rms = float(np.sqrt(np.mean((ts_cpu - ts_gpu) ** 2))) + assert ts_rms < ts_rel_tol * ts_scale, ( + f'timeseries rms={ts_rms:.3e} > {ts_rel_tol} * scale={ts_scale:.3f}' + ) + + tcoh_rms = float(np.sqrt(np.mean((tcoh_cpu - tcoh_gpu) ** 2))) + assert tcoh_rms < tcoh_abs_tol, ( + f'temporalCoherence rms={tcoh_rms:.3e} > {tcoh_abs_tol:.3e}' + ) + + np.testing.assert_array_equal(nobs_cpu, nobs_gpu) + + +@pytest.fixture(scope='module') +def network(): + """Shared SBAS-like network sized to match FernandinaSenDT128 + (num_date=98, num_pair=288). max_span=6 reproduces the typical + short-baseline coverage of a Sentinel-1 SBAS network and provides + enough redundancy that small fractions of NaN observations do not + push individual pixels into rank-deficiency. + """ + return make_redundant_network(num_date=98, num_pair=288, max_span=6, seed=0) + + +def _all_pixels_full_rank(A, y): + """Return True if every pixel's design matrix (after dropping NaN rows) + is still full-rank. Used to keep tests off the rank-deficient edge case + (which is handled separately at runtime by a CPU fallback path). + """ + n_unknown = A.shape[1] + for k in range(y.shape[1]): + valid = ~np.isnan(y[:, k]) + if np.linalg.matrix_rank(A[valid]) < n_unknown: + return False + return True + + +@requires_cuda +def test_wls_no_nan(network): + """WLS, no NaN observations — expect ~ float32 round-off match.""" + A, B, tbase_diff = network + y, w = synthesize_observations(A, B, num_pixel=64, nan_frac=0.0, seed=1) + cpu = cpu_reference(A, B, y, w, tbase_diff) + gpu = estimate_timeseries_batch( + A=A, B=B, y=y, tbase_diff=tbase_diff, weight_sqrt=w, + min_norm_velocity=True, + chunk_size=64, solver='torch', print_msg=False, + ) + assert_equivalent(cpu, gpu, ts_rel_tol=1e-5, tcoh_abs_tol=1e-5) + + +@requires_cuda +def test_wls_with_nan_redundant(network): + """WLS with low NaN rate on a redundant network — float32 round-off match. + + NaN fraction is kept at 3% so each pixel still has enough observations to + keep its (NaN-masked) design matrix full-rank. Higher NaN rates exercise + the rank-deficient edge case, which is handled by a separate CPU fallback + path (not yet implemented at the time of writing). + """ + A, B, tbase_diff = network + y, w = synthesize_observations(A, B, num_pixel=64, nan_frac=0.03, seed=2) + assert _all_pixels_full_rank(A, y), \ + 'fixture broke: a pixel is rank-deficient after NaN masking' + cpu = cpu_reference(A, B, y, w, tbase_diff) + gpu = estimate_timeseries_batch( + A=A, B=B, y=y, tbase_diff=tbase_diff, weight_sqrt=w, + min_norm_velocity=True, + chunk_size=32, solver='torch', print_msg=False, + ) + assert_equivalent(cpu, gpu, ts_rel_tol=1e-4, tcoh_abs_tol=1e-4) + + +@requires_cuda +def test_ols_no_nan(network): + """OLS path (weight_sqrt=None).""" + A, B, tbase_diff = network + y, _ = synthesize_observations(A, B, num_pixel=64, nan_frac=0.0, + weighted=False, seed=3) + cpu = cpu_reference(A, B, y, None, tbase_diff) + gpu = estimate_timeseries_batch( + A=A, B=B, y=y, tbase_diff=tbase_diff, weight_sqrt=None, + min_norm_velocity=True, + chunk_size=32, solver='torch', print_msg=False, + ) + assert_equivalent(cpu, gpu, ts_rel_tol=1e-5, tcoh_abs_tol=1e-5) + + +@requires_cuda +def test_min_norm_phase(network): + """min_norm_velocity=False solves on A directly; ts[1:] = X.""" + A, B, tbase_diff = network + y, w = synthesize_observations(A, B, num_pixel=64, nan_frac=0.0, seed=4) + cpu = cpu_reference(A, B, y, w, tbase_diff, min_norm_velocity=False) + gpu = estimate_timeseries_batch( + A=A, B=B, y=y, tbase_diff=tbase_diff, weight_sqrt=w, + min_norm_velocity=False, + chunk_size=32, solver='torch', print_msg=False, + ) + assert_equivalent(cpu, gpu, ts_rel_tol=1e-5, tcoh_abs_tol=1e-5) + + +@requires_cuda +def test_chunk_size_invariance(network): + """Output must be effectively identical across chunk sizes.""" + A, B, tbase_diff = network + y, w = synthesize_observations(A, B, num_pixel=64, nan_frac=0.03, seed=5) + assert _all_pixels_full_rank(A, y) + common = dict(A=A, B=B, y=y, tbase_diff=tbase_diff, weight_sqrt=w, + min_norm_velocity=True, solver='torch', print_msg=False) + ts_a, tcoh_a, nobs_a = estimate_timeseries_batch(chunk_size=16, **common) + ts_b, tcoh_b, nobs_b = estimate_timeseries_batch(chunk_size=64, **common) + np.testing.assert_allclose(ts_a, ts_b, rtol=1e-6, atol=1e-6) + np.testing.assert_allclose(tcoh_a, tcoh_b, rtol=1e-6, atol=1e-6) + np.testing.assert_array_equal(nobs_a, nobs_b) + + +@requires_cuda +def test_unsupported_solver_raises(network): + A, B, tbase_diff = network + y, w = synthesize_observations(A, B, num_pixel=8, seed=6) + with pytest.raises(ValueError, match='unsupported solver'): + estimate_timeseries_batch( + A=A, B=B, y=y, tbase_diff=tbase_diff, weight_sqrt=w, + solver='cupy', print_msg=False, + ) From dd998ca9c935515c2d61334f38c96dcd3f16b5b6 Mon Sep 17 00:00:00 2001 From: s-sasaki-earthsea-wizard Date: Wed, 6 May 2026 18:52:16 +0900 Subject: [PATCH 2/3] requirements-gpu: drop unused cupy and pin torch>=2.11 The opt-in GPU solver in `ifgram_inversion_gpu.py` is implemented entirely on top of `torch.linalg.cholesky_ex`, with no cupy entry point. Listing `cupy-cuda12x` in `requirements-gpu.txt` therefore pulls a multi-hundred-MB runtime that no code path imports. Drop it. Pin `torch>=2.11` to match the version exercised in the bench matrix used during development (Blackwell sm_120 wheel from the cu128 index). Earlier torch releases have not been validated against this code path. --- requirements-gpu.txt | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/requirements-gpu.txt b/requirements-gpu.txt index 7d757f416..3d31e677f 100644 --- a/requirements-gpu.txt +++ b/requirements-gpu.txt @@ -1,5 +1,4 @@ # Optional GPU acceleration deps for the `[gpu]` extras. # Install with the PyTorch CUDA wheel index, e.g.: # pip install -e ".[gpu]" --extra-index-url https://download.pytorch.org/whl/cu128 -cupy-cuda12x -torch +torch>=2.11 From 97d154f495b519c08d9db8d51399cbcfdf023c41 Mon Sep 17 00:00:00 2001 From: s-sasaki-earthsea-wizard Date: Wed, 6 May 2026 18:53:11 +0900 Subject: [PATCH 3/3] =?UTF-8?q?docs:=20add=20gpu.md=20and=20installation?= =?UTF-8?q?=20=C2=A72.4=20for=20opt-in=20torch=20solver?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Document the new opt-in `torch` GPU solver added in the previous commits: - `docs/gpu.md` — setup, CLI / template surface, behavior notes (VRAM auto-sizing, rank-deficient pixel handling, NaN observations, hard-fail on missing CUDA), and indicative performance numbers. - `docs/installation.md` §2.4 — install the `[gpu]` extras together with the matching PyTorch CUDA wheel index. - `docs/README.md` and `docs/dask.md` — add cross-links so readers can reach the GPU page from the documentation root and from the Dask page (since the two parallelism paths are orthogonal and need to be picked one or the other). Performance numbers in `gpu.md` §4 are stated inline without any external repository links so the page stays self-contained. --- docs/README.md | 2 ++ docs/dask.md | 2 ++ docs/gpu.md | 68 ++++++++++++++++++++++++++++++++++++++++++++ docs/installation.md | 58 +++++++++++++++++++++++++++++++++++++ 4 files changed, 130 insertions(+) create mode 100644 docs/gpu.md diff --git a/docs/README.md b/docs/README.md index d87251cce..c541362c3 100644 --- a/docs/README.md +++ b/docs/README.md @@ -88,6 +88,8 @@ Algorithms implemented in the software are described in details at [Yunjun et al + [Example data directory](./dir_structure.md) + [Example template files](./templates/README.md) + [Tutorials in Jupyter Notebook](https://github.com/insarlab/MintPy-tutorial) ++ [Parallel processing with Dask](./dask.md) ++ [GPU acceleration for the `invert_network` step (opt-in PyTorch CUDA solver, partial)](./gpu.md) ### 4. Contact us diff --git a/docs/dask.md b/docs/dask.md index db3018bd0..6009b1f2b 100644 --- a/docs/dask.md +++ b/docs/dask.md @@ -7,6 +7,8 @@ Most computations in MintPy are operated in either a pixel-by-pixel or a epoch-b [Here](https://github.com/2gotgrossman/dask-rsmas-presentation) is an entry-level presentation on parallel computing using Dask by David Grossman. Below we brief describe for each cluster/scheduler the required options and recommended best practices. +For GPU acceleration of the `invert_network` step on a single CUDA device — orthogonal to the Dask paths described here — see [gpu.md](./gpu.md). + ## 1. local cluster ## The parallel processing on a single machine is supported via [`Dask.distributed.LocalCluster`](https://docs.dask.org/en/latest/setup/single-distributed.html#localcluster). This is recommended if you are running MintPy on a local machine with multiple available cores, or on an HPC but wish to allocate only a single node's worth of resources. diff --git a/docs/gpu.md b/docs/gpu.md new file mode 100644 index 000000000..844389fa5 --- /dev/null +++ b/docs/gpu.md @@ -0,0 +1,68 @@ +# Configure GPU acceleration for the network inversion # + +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. + +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.1 via command line #### + +Run the following in the terminal: + +```bash +ifgram_inversion.py inputs/ifgramStack.h5 --solver torch +ifgram_inversion.py inputs/ifgramStack.h5 --solver torch --gpu-chunk-size 20000 +``` + +`--gpu-chunk-size 0` (the default) auto-sizes the per-chunk pixel count from free VRAM; pass a positive integer to override. + +#### 2.2 via template file #### + +Adjust options in the template file: + +```cfg +mintpy.networkInversion.solver = torch #[cpu / torch], auto for cpu +mintpy.networkInversion.gpuChunkSize = auto #[int >= 0], auto for 0 (auto-size from free VRAM) +``` + +and feed the template file to the script: + +```bash +ifgram_inversion.py inputs/ifgramStack.h5 -t smallbaselineApp.cfg +smallbaselineApp.py smallbaselineApp.cfg +``` + +#### 2.3 Testing using example data #### + +Download and run the FernandinaSenDT128 example data; then run with and without the GPU solver: + +```bash +cd FernandinaSenDT128/mintpy +ifgram_inversion.py inputs/ifgramStack.h5 -w no --solver cpu +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 ## + ++ **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). + ++ **Rank-deficient pixels.** Detected via `torch.linalg.cholesky_ex` info codes; their solution is set to zero so NaN/Inf never propagate downstream. A warning line reports the count per chunk. + ++ **Per-pixel NaN observations.** Handled by zeroing the corresponding row weight, which is mathematically equivalent to dropping that row from the WLS system. + ++ **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 ## + +Indicative numbers below were measured on an NVIDIA RTX 5080 (Blackwell sm_120, CUDA 12.8, PyTorch 2.11) at the time this feature was submitted. Speedup depends on scene size, GPU class, and chunk-size tuning, so reproduce on your own data and hardware before drawing conclusions. + ++ **Tutorial-scale** (FernandinaSenDT128: 270k pixels, 288 ifgs) — `invert_network` runs roughly **16×** faster internally and **4.5×** faster end-to-end versus the CPU path. ++ **Large-scene** (GalapagosSenDT128: 3.4M pixels, 475 ifgs; ~12.6× pixels and 1.65× ifgs over Fernandina) — roughly **44×** internal and **36×** step-wall speedup on `invert_network` (CPU 6189 s → torch 170 s on the same machine), confirming the speedup grows at scale. ++ **Numerical equivalence** between the `cpu` and `torch` solvers holds to float32 round-off: RMS on the order of `1e-5` on the tutorial case, with absolute RMS at most ~16 µm on the large-scene case. diff --git a/docs/installation.md b/docs/installation.md index 46d55f2e9..1a831af27 100644 --- a/docs/installation.md +++ b/docs/installation.md @@ -194,6 +194,64 @@ Same as the instruction for Linux, except for

+### 2.4 Optional: GPU acceleration via PyTorch CUDA ### + +

+

+

Click to expand for more details

+ +

The invert_network step ships an opt-in GPU solver that solves the per-pixel WLS inversion as a batched normal-equation + Cholesky on a CUDA device via PyTorch. It is opt-in: the default mintpy.networkInversion.solver = auto resolves to cpu, so existing setups are unaffected. There is no silent CPU fallback — selecting torch without a visible CUDA device is a hard error.

+ +

a. Prerequisites

+ +
    +
  • An NVIDIA GPU with a working CUDA driver
    • +
  • MintPy installed from source in editable mode (Section 2 above)
    • +
+ +

b. Install the [gpu] extras

+ +

The [gpu] extras pull a CUDA-enabled PyTorch build. Pick the cuXXX wheel index that matches your CUDA toolkit version (see pytorch.org for the current list); for example, cu121, cu124, or cu128:

+ +```bash +python -m pip install -e ".[gpu]" \ + --extra-index-url https://download.pytorch.org/whl/cu128 +``` + +

If you use uv instead of pip, add --index-strategy unsafe-best-match to work around a stale setuptools pin in the PyTorch wheel index:

+ +```bash +uv pip install -e ".[gpu]" \ + --extra-index-url https://download.pytorch.org/whl/cu128 \ + --index-strategy unsafe-best-match +``` + +

c. Verify

+ +```bash +python -c "import torch; print(torch.cuda.is_available())" +# expected: True +``` + +

d. Enable

+ +

Set the template flag:

+ +```cfg +mintpy.networkInversion.solver = torch +``` + +

or pass it on the command line:

+ +```bash +ifgram_inversion.py inputs/ifgramStack.h5 --solver torch +``` + +

See gpu.md for tuning, behavior notes, and benchmarks.

+ +
+

+ ## 3. Post-Installation Setup ## #### a. ERA5 for tropospheric correction ####