Research map / Architectures and efficiency / RNNs and normalization
Schwarz–Schur Involution: Lightspeed Differentiable Sparse Linear Solvers
Read the paperCheck your idea against 108,967 papers
AI-extracted summary. Written by AI from the paper's full text. Numbers and claims can be wrong or missing context, so check the paper before you cite or build on it.
Proposes Schwarz–Schur involution, a batched dense-tensor implementation of nested-dissection Gaussian elimination that makes exact direct sparse solves on image-grid Laplacian-like systems up to ~1000x faster than SciPy and ~40–170x faster than cuDSS, usable as a differentiable module.
Problem
Direct sparse linear solvers are accurate and reliable but considered too slow for neural architectures and interactive vision/graphics, while iterative solvers are problem-dependent, parameter-sensitive, and have unpredictable runtime.
Approach
Condense a sparse grid matrix into a batched 4-D tensor of Dirichlet-to-Neumann (Schur complement) blocks, then recursively merge adjacent subdomains via batched dense inversions (Schwarz step plus Schur steps), a parallel nested-dissection Gaussian elimination that exploits dense GPU BLAS kernels.
Key results
- On a 2561x2561 Laplacian system, the method takes 220 ms versus 25612 ms for CUDA and 36926 ms for SciPy, with 168x and 1151x speedups per the reported table (the ~1000x claim is up to three orders of magnitude).
- On 513x513 images the method solves in 10.9 ms, a 95x speedup over CUDA and 188x over SciPy.
- On Darcy flow, the zero-shot solver achieves 0.869% ± 0.0927% relative error versus 1.56% reported for the learned neural operator, in 5-7 ms at 257x257 resolution.
Measured effects reported in the paper
| What was changed | Where | Measured | Effect |
|---|---|---|---|
| Use Schwarz–Schur Involution instead of SciPy | Laplacian Dirichlet solver · Image grids from 257×257 to 2561×2561 | runtime (ms): 253318 -> 220.1; 143100 -> 158.5; 16512 -> 36.45; 2051 -> 10.90; 355 -> 5.82 | Helped |
| Use Schwarz–Schur Involution instead of CUDA (cuDSS) | Laplacian Dirichlet solver · Image grids from 257×257 to 2561×2561 | runtime (ms): 36926 -> 220.1; 21354 -> 158.5; 4710 -> 36.45; 1036 -> 10.90; 234 -> 5.82 | Helped |
| Use Schwarz–Schur Involution instead of SciPy-LU | Sparse solver · Laplacian Dirichlet 513×513 | runtime (ms): 2208 -> 12.3 | Helped |
| Use Schwarz–Schur Involution instead of CUDSS-LDU | Sparse solver · Laplacian Dirichlet 513×513 | runtime (ms): 1011 -> 12.3 | Helped |
| Use Schwarz–Schur Involution instead of SciPy-LU | Sparse solver · Helmholtz κ²=1/64×0.1, 513×513 | runtime (ms): 2687 -> 35.5 | Helped |
| Use Schwarz–Schur Involution instead of CUDSS-LDU | Sparse solver · Helmholtz κ²=1/64×0.1, 513×513 | runtime (ms): 1037 -> 35.5 | Helped |
| Use Schwarz–Schur Involution without training instead of the learned solver | Darcy flow coefficient-to-solution mapping · Darcy flow, 257×257 | relative error: 1.56% -> 0.869% ± 0.0927% | Helped |
| Use Schwarz–Schur Involution instead of SciPy | Generalized deconvolution · 257×257 images | runtime (ms): 494.4ms -> 7.51ms | Helped |
| Use Schwarz–Schur Involution instead of CG | Laplacian solver · 513×513 problems, Lap 1e-2 | runtime (ms): 65.1 -> 12.3 | Helped |
| Use Schwarz–Schur Involution instead of CG | Helmholtz solver · 513×513 problems, Hel 1.0 | runtime (ms): 6442 -> 35.5 | Helped |
| Use Schwarz–Schur Involution instead of biCGstab | Helmholtz solver · 513×513 problems, Hel 0.1 | runtime (ms): NaN -> 35.5 | Helped |
| Increase anisotropy in the Laplacian kernel | Iterative solvers · Anisotropic versus isotropic Laplacian | number of iterations: 6× | Hurt |
Surprising findings
- The direct solver's error in float-64 is consistently smaller than SciPy's SuperLU on elliptic PDEs, despite not optimizing for accuracy.
- A naive PyTorch prototype with torch.linalg.inv on batched dense blocks outperforms optimized sparse direct solvers and iterative solvers by orders of magnitude.
Limitations
- Higher memory consumption than general-purpose sparse solvers, and designed for single-GPU problems around 1M variables rather than billions.
- Currently supports one variable per pixel and 2D image-grid stencils (9-point); extensions to 3D, larger stencils, and arbitrary meshes are future work.
- Gradients via adjoint method require transposed solve; speedups depend on reuse of factorization.
Experimental setup
Datasets: neuraloperator Darcy flow examples
Benchmarks: Laplacian (Dirichlet) systems on 513x513 to 2561x2561 images Helmholtz systems random sparse matrices Darcy flow (neuraloperator) image matting Laplacian (PyMatting) generalized deconvolution with spatially varying kernels
Metrics: runtime (ms) relative residual error ||Ax-b||/||b|| speedup
Baselines: SciPy spsolve (SuperLU) CUDA cuDSS LDU PyAMG Ruge-Stuben AMG CG MINRES GMRES LSQR biCGstab PCG+AMG (AMGX) FNO (Li et al., 2024)
Builds on: nested dissection multifrontal solvers domain decomposition (Schwarz alternating method) adjoint method for differentiable linear solvers
Closest papers by meaning
- Learning Sparse Approximate Inverse Preconditioners for Conjugate Gradient Solvers on GPUs (NeurIPS 2025 · similarity 0.85)
- Provable Acceleration of Nesterov's Accelerated Gradient for Asymmetric Matrix Factorization and Linear Neural Networks (NeurIPS 2024 · similarity 0.85)
- Toward Efficient Kernel-Based Solvers for Nonlinear PDEs (ICML 2025 · similarity 0.85)
- Dynamic Sparse Training of Diagonally Sparse Networks (ICML 2025 · similarity 0.85)
- Searching for Efficient Linear Layers over a Continuous Space of Structured Matrices (NeurIPS 2024 · similarity 0.84)
- Optimal deep learning of holomorphic operators between Banach spaces (NeurIPS 2024 · similarity 0.84)
Writing a paper in this area? scime draft is a free LaTeX editor with AI included, real-time co-authors and a check of your references and venue rules before you submit. It also imports OpenAI Prism projects.