Fast Sparse Matrix Permutation for Mesh-Based Direct Solvers

Ahmed H. Mahmoud; Ana Dodik; Behrooz Zarebavami; Changcheng Yuan; John D. Owens; Justin Solomon; Maryam Mehri Dehnavi; Serban D. Porumbescu

REVIEW 1 major objections 6 minor 2 references

Fast Sparse Matrix Permutation for Mesh-Based Direct Solvers

T0 review · 1 major / 6 minor · reviewed 2026-08-03 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read A mesh-aware, patch-guided nested-dissection permutation cuts fill-reducing ordering time by up to 10.27x and end-to-end sparse-Cholesky solve time by up to 6.62x on large triangle-mesh systems without unacceptable fill-in.

desk verdict Solid and genuinely new systems work with transparent ablations, but the unverified elimination tree handed to cuDSS is a real gap that needs fixing before the end-to-end claims can be trusted. read the letter →

arxiv 2602.00898 v2 pith:KYIESK43 submitted 2026-01-31 cs.GR cs.DC

Behrooz Zarebavami , Ahmed H. Mahmoud , Ana Dodik , Changcheng Yuan , Serban D. Porumbescu , John D. Owens , Maryam Mehri Dehnavi , Justin Solomon This is my paper · ORCID

classification cs.GRcs.DC MSC 65F0565F50

keywords fill-reducingpermutationnesteddissectionsparseCholeskyeliminationtreetrianglemeshesquotientgraphdirectsolversGPU

verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper sets out to show that for linear systems coming from triangle meshes, the most expensive part of a sparse direct solver—the fill-reducing permutation—can be made dramatically faster by deliberately relaxing the two design goals that general-purpose nested-dissection tools chase: strict balance and minimal separator size. Instead of coarsening the whole graph level by level, the algorithm clusters mesh vertices into patches, builds one compact quotient graph, and runs nested dissection on that small graph at every recursion. It also constructs the elimination tree on the way, so the solver's symbolic phase can reuse that hierarchy. A sympathetic reader would see the claim as: for mesh-shaped matrices, preprocessing time can drop by an order of magnitude while the resulting factorization remains efficient enough to yield end-to-end speedups on big problems.

What carries the argument

The central object is the quotient graph of a mesh-patch partition: nodes are connected patches of mesh vertices, edges record adjacency across patches, and each patch is lifted to a group of vertices in the matrix graph. Separator search is done on this quotient graph at every recursion level; a bipartition of the quotient graph is lifted to a separator superset of actual graph vertices and then locally refined. Because the same patch structure is reused across all recursive calls, the algorithm avoids repeated multilevel coarsening. The elimination tree is built as the recursion produces separators, and the final permutation is a concatenation of local orderings (some minimum-degree style)

What would settle it

On a 1.5M-vertex triangle-mesh Laplacian, measure ordering time and the Cholesky fill ratio for the patch-guided permutation against a balanced nested-dissection permutation on the same hardware; the claim fails if the ordering-time advantage disappears or the fill ratio grows by more than the paper's own reported bounds (roughly +15%) at large patch counts.

Watch

Extended reading notes

Core claim

On the paper's own terms, the discovery is that a nested-dissection permutation for a mesh-derived matrix can be computed by first partitioning the mesh into connected patches, lifting that partition to the matrix graph as a group map, recursing only on the resulting quotient graph, and then refining the lifted separator sets with a local ordering step such as approximate minimum degree. This patch-guided recursion produces a near-nested-dissection ordering and its accompanying elimination tree in one pass. Because the quotient graph is much smaller than the input graph and is reused at every recursion level, separator computation and tree construction become cheap; the price is a modest inc

Load-bearing premise

The speedups depend on the matrix's sparsity graph matching the triangle-mesh adjacency graph, so that patches of mesh vertices correspond to groups of matrix rows/columns; for matrices without that correspondence (e.g., higher-order elements or constrained degrees of freedom), the patch-guided separator search loses its foundation.

Editorial extensions

If this is right

On meshes above roughly 50k–100k vertices, the patch-guided ordering beats general-purpose ordering tools in runtime, with up to a 10.27x ordering speedup and a 4.58x geometric mean.
End-to-end sparse-Cholesky solves on large meshes speed up by up to 6.62x on a GPU solver and 2.55x on a CPU solver, with the gains growing as permutation dominates preprocessing.
The elimination tree is produced together with the permutation, so symbolic analysis can reuse it, removing a separate tree-construction step.
The schedule used to concatenate local orderings can be chosen without changing fill-in, allowing the same permutation to adapt to different execution models—post-order for CPU locality, level-order/wavefront for GPU parallelism.
For repeated-factorization workloads, the speedup persists beyond a break-even iteration count (reported 148–160 iterations for mesh smoothing).

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

The same patch-guided separator mechanism could plausibly transfer to other mesh types (tetrahedral, hexahedral) if fast volumetric patching becomes available; the paper's preliminary tet experiment hints at this direction.
Because the schedule is a free parameter that does not affect fill-in, the approach opens an easy path to solver-specific scheduling optimization without recomputing the ordering.
The fill-in–runtime trade-off could be tuned per application by adjusting patch size; large patches favor speed, small patches favor factor quality, suggesting an auto-tuning knob for different solver pipelines.

Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Desk Editor's Note

read the letter

You should know two things about this paper. First, the core idea is genuinely new and the empirical work is transparent enough to be useful: patch-level nested dissection on a single reused quotient graph, with the elimination tree built as a byproduct and a schedule that can be swapped for CPU or GPU solvers. When each patch is a single vertex, it reduces to classical nested dissection. The authors integrate into cuDSS and MKL, leave the numerical factorization untouched, isolate the permutation effect, and even give break-even iteration counts in Table 1. That is the right way to report an ordering-quality versus runtime trade-off. Second, there is one gap that needs fixing before I'd trust the end-to-end speedups: they hand cuDSS an elimination tree whose leaves are single nodes containing thousands of vertices that were locally reordered with AMD. That collapsed tree is not the true elimination tree of the permuted matrix, and the paper never verifies that the solves are actually correct—no residuals, no comparison against a reference solver, no check that the separators are valid vertex separators. If cuDSS schedules internal elimination of those leaf blocks without respecting the AMD order, the factorization can be wrong. The applications would likely look terrible in that case, so the reported speedups are probably on real solutions, but "probably" isn't good enough for a paper whose central claim is an end-to-end speedup.

The smaller issues are easy to fix but real: the abstract says 6.27x, the caption and intro say 6.62x; the citation for AMD points to Tobler & Maierhofer (a mesh data structure paper) instead of Amestoy et al., and the NP-completeness citation is to Ashkiani et al. (a dynamic hash table paper); the normal-equation derivation in §6.2 is off by a factor—it should be (M + α L^T M^{-1} L)x = M y, not the printed form. Also, no run-to-run variance and no full input list. None of these breaks the core claim.

What holds up: the central claim is credible. The ablations isolate the contribution of each design choice, the baseline comparisons are against real tools, and the break-even analysis is a nice touch. The method is clearly scoped to triangle-mesh SPD systems, and the authors say so.

Who this is for: anyone working on sparse direct solvers for mesh-derived systems, especially on GPUs with cuDSS. Send it to peer review; require a correctness check on the etree and fix the mechanical issues before acceptance.

Referee Report

1 major / 6 minor

Summary. The paper proposes a patch-based nested-dissection permutation for sparse SPD systems whose sparsity graph is the adjacency graph of a triangle mesh. It computes a patch partition, lifts it to a group map, recursively bipartitions the quotient graph to obtain separators, builds a coarse elimination tree, and assembles the final permutation by concatenating local orderings (AMD) in a chosen schedule. The authors integrate the method into NVIDIA cuDSS and Intel MKL, replacing the default ordering and, for cuDSS, the elimination-tree construction. They report ordering-time speedups up to 10.27× vs METIS and end-to-end solve speedups up to 6.62×, with ablations showing the effect of patch size, separator computation, and local ordering.

Significance. If the claims are correct, this is a practically useful result for graphics workloads where the symbolic phase dominates sparse Cholesky. The paper is unusually transparent: it leaves the numerical factorization untouched when measuring integration effects, isolates the contributions of patch size, separator computation, and local ordering, reports break-even iteration counts, and releases code. The main gap is that the correctness of the cuDSS integration—specifically, the validity of the user-supplied elimination tree—is never verified, which is a prerequisite for any end-to-end timing claim.

major comments (1)

[§5, §4.2–4.3, §6] The central end-to-end speedup claims for cuDSS presuppose that the elimination tree handed to cuDSS is a valid schedule for the permuted matrix. The etree is built by making each separator and each leaf subgraph a single node, while leaves are locally reordered with AMD (§4.3). The internal elimination dependencies inside a leaf are therefore not represented in the etree. §5 warns that an incompatible schedule can cause runtime failure, but the evaluation reports no numerical residuals, no comparison against a reference direct solve, and no check that the separators used at each level are valid vertex separators of the current subgraph. Without such validation, the timings in Table 1 and Fig. 7 could reflect incorrect solves rather than genuine speedups. Please add residual/normwise-error checks for all benchmark configurations and verify (or construct) the exact etree of the final perm

minor comments (6)

[§1, abstract] The abstract reports 'up to 6.27×' end-to-end speedup, while §1 and Fig. 1 report 'up to 6.62×'. Please reconcile.
[§2, references] AMD is cited as [Tobler and Maierhofer 2006], which is a mesh data structure paper; the correct reference for AMD is Amestoy et al. 1996. Similarly, §3 cites [Ashkiani et al. 2018] for NP-completeness of fill-in minimization; the correct citation is Yannakakis 1981.
[§6.5] The sentence 'Hence, speedups here are smaller than for mesh smoothing (§6.4), permutation accounts for a larger share of runtime' contradicts the immediately preceding observation that increased density reduces the fraction of time spent in permutation. Presumably 'larger' should be 'smaller.'
[Algorithm 1 / §4] Algorithm 1 lists target_patch_size as a parameter only implicitly; the input line says 'A, M, nd_level' but GetPatches uses a patch-size parameter. Please include all tunable parameters explicitly and report the nd_level sensitivity. The paper sets nd_level to 9–10 with a citation but provides no ablation for it.
[§5] The text says the quotient graph is built once and then updated, but Algorithm 2 calls QuotientGraph(g, gmap) at every recursion level. Please clarify whether the initial construction is reused and only node/edge weights are updated.
[§2 / title] The title and abstract emphasize 'Sparse Matrix Permutation' broadly; the method actually requires the matrix sparsity graph to coincide with a triangle-mesh adjacency graph (or a block-wise contraction of it). A short formal statement of this precondition would make the contribution easier to evaluate and avoid overgeneralization.

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity; the central claims are empirical speedups measured against external solvers and ordering tools.

full rationale

The paper's central claim is an empirical performance improvement: replacing the fill-reducing permutation (and elimination-tree construction for cuDSS) with the proposed patch-guided nested-dissection ordering yields measured speedups against METIS, AMD, ParMETIS, Intel MKL, and NVIDIA cuDSS. No fitted parameter is renamed as a prediction, and no derived quantity reduces by construction to an input. The patch size and nd_level are hyperparameters tuned by ablation and explicitly disclosed as such (§6.6), which is standard empirical configuration rather than circular reasoning. Self-citations appear (RXMesh for GPU patch computation, Parth for elimination-tree reuse), but they are implementation and related-work references; they do not carry the load of the main performance claim, and no uniqueness theorem or ansatz is imported from them to force the algorithm's choices. The paper also transparently reports limitations (§5 disconnected groups degrade quality; §7 tet-mesh slowdowns), further indicating that the results are not manufactured by definition. Potential concerns about whether the supplied elimination tree is exactly valid for cuDSS are correctness/verification risks, not circularity, since no equation in the paper is equivalent to its own input.

Assumptions & free parameters 2 free parameters · 4 assumptions · 0 invented entities

The algorithm introduces no new physical or mathematical entity; it borrows patches, quotient graphs, METIS, AMD, and elimination trees. The main free choices are patch size and recursion depth, both tuned empirically. The load-bearing assumptions are the mesh-to-matrix correspondence and the usability of the patch partition as a coarse separator space; both are stated in the paper, but neither is formally guaranteed.

free parameters (2)

target_patch_size = 256 vertices (default)
Chosen via ablation in §6.6 (Fig. 11) to balance permutation runtime against fill-in; the reported end-to-end speedups use this default.
nd_level (nested-dissection recursion depth) = not specified for the experiments; text mentions moderate depths 9-10
Algorithm 1 takes nd_level as input; it controls parallelism versus separator quality, and the paper does not state the exact value used in the evaluation.

assumptions (4)

domain assumption The sparsity graph of A corresponds to the triangle-mesh adjacency graph (or to a block-wise version with 2x2 blocks merged), so mesh patches lift to graph groups.
Invoked in §2: 'there is a one-to-one correspondence between mesh vertices and unknowns... G matches the mesh adjacency graph.' This is the foundation of the patch-to-graph lifting in Algorithm 1.
domain assumption The patch partition is connected (or connected enough) that quotient-graph bipartitions yield useful separators when lifted.
Stated in §5: disconnected groups 'typically produce larger and less balanced separators.' The default GPU k-means patching is assumed to produce suitable patches, but connectivity is not proven.
standard math Standard nested-dissection and elimination-tree theory (Lipton-Rose-Tarjan), including fill behavior and Cholesky complexity O(sum d(i)^2), applies.
Used throughout §2 and §6.5 to motivate why separator hierarchy and fill-in matter; standard background for sparse direct solvers.
domain assumption METIS bipartitioning and node refinement on the quotient graph produce a valid, acceptable separator after lifting and refinement.
Algorithm 2 relies on METIS for BiPartition and RefineSeparator; no proof is given that the lifted separator always separates the subgraph with the claimed balance and size.

how reviews work

0 comments

Cite this review

Pith. "Pith review of Fast Sparse Matrix Permutation for Mesh-Based Direct Solvers." pith.science (2026). https://pith.science/paper/KYIESK43

@misc{pith2026260200898,
  author       = {Pith},
  title        = {Pith review of: Fast Sparse Matrix Permutation for Mesh-Based Direct Solvers},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/KYIESK43}},
  note         = {Machine review of arXiv:2602.00898}
}

read the original abstract

We present a fast sparse matrix permutation algorithm tailored to linear systems arising from triangle meshes. Our approach produces nested-dissection-style permutations while significantly reducing permutation runtime overhead. Rather than enforcing strict balance and separator optimality, the algorithm deliberately relaxes these design decisions to favor fast partitioning and efficient elimination-tree construction. Our method decomposes permutation into patch-level local orderings and a compact quotient-graph ordering of separators, preserving the essential structure required by sparse Cholesky factorization while avoiding its most expensive components. We integrate our algorithm into vendor-maintained sparse Cholesky solvers on both CPUs and GPUs. Across a range of graphics applications, including single factorizations and repeated factorizations, our method reduces permutation time and improves the sparse Cholesky solve performance by up to 6.27x. Our code is available at https://github.com/BehroozZare/fast-permute.

Figures

Figures reproduced from arXiv: 2602.00898 by the authors.

**Figure 1.** Scalable fill-reducing permutation computation on a triangular mesh. Our patch-based, nested-dissection-style method compresses the triangular mesh domain into a small quotient graph and reuses this compression across recursive calls, improving scalability over state-of-the-art methods. We achieve a 6.62× speedup over the highly optimized NVIDIA cuDSS by accelerating its permutation computation. Sparse Cholesky so… view at source ↗

**Figure 2.** Fraction of time spent on matrix permutation to end-to-end time for solving mean curvature flow [Desbrun et al. 1999] for different inputs. Sparse linear solvers are central computational kernels in computer graphics. Many mesh-centric workloads (e.g., parameterization, deformation, physical simulation, and geometric flows) solve large positive semidefinite systems, for which sparse Cholesky factorization is of… view at source ↗

**Figure 3.** Nested dissection example. Illustration of a single nested-dissection step. From left to right: the input matrix 𝐴, its graph 𝐺, the corresponding elimination tree 𝑇 , the permuted graph 𝐺 ′ , and the sparsity pattern of the Cholesky factor 𝐿. The separator vertices 𝑆 are ordered after the vertices in 𝑔1 and 𝑔2. The elimination tree encodes the dependencies induced by this order: computations for 𝑔1 and 𝑔2 are indep… view at source ↗

Figures from the paper (8 more)

**Figure 4.** Figure 4: Elimination tree construction. We compute the root separator {7, 2, 5}. For the resulting subgraphs gℓ = {1, 4, 6} and g𝑟 = {3, 8, 9}, we recurse to compute the next-level separators (e.g., {4} and {8}), continuing until reaching the desired depth. We store the full bi…

**Figure 5.** Figure 5: Separator computation. For the current subgraph g, we form a quotient graph q, bipartition q, lift the partition back to 𝑉g to obtain a separator superset, refine it to reduce its size while maintaining balance, and then extract the left and right subgraphs for the nex…

**Figure 6.** Figure 6: Permutation vectors from post-order (left) and level-order (right) tra [PITH_FULL_IMAGE:figures/full_fig_p005_6.png]

**Figure 9.** Figure 9: Separator computation: METIS vs. patch-based in our algorithm. We compare METIS vertex-separator computation against our algorithm’s patch-based separator stage at patch size 256, on a small mesh and a large mesh. Small-Patch and Small-METIS denote the corresponding …

**Figure 12.** Figure 12: Data smoothing via Laplacian regularization. We recover a smooth [PITH_FULL_IMAGE:figures/full_fig_p009_12.png]

**Figure 13.** Figure 13: Effect of patch size on separator quality. Decreasing the patch size results in a larger number of patches, allowing our algorithm to produce higher-quality separators, since the separator must follow patch boundaries. However, this comes with a trade-off where increa…

**Figure 15.** Figure 15: Results of six iterations of mesh smoothing applied on a mesh with [PITH_FULL_IMAGE:figures/full_fig_p010_15.png]

**Figure 16.** Figure 16: End-to-end solver speedup mesh smoothing. Using cuDSS, our fill-reducing permutation achieves an average 2.29× speedup across six smoothing iterations, with larger meshes benefiting the most. 2026-02-03 02:00. Page 10 of 1–10 [PITH_FULL_IMAGE:figures/full_fig_p010_16.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

2 extracted references

[2008]

Algorithm 887: CHOLMOD, Supernodal Sparse Cholesky Factorization and Update/Downdate.ACM Trans. Math. Softw.35, 3, Article 22 (Oct. 2008), 14 pages. doi:10.1145/1391989.1391995 Kazem Cheshmi, Danny M. Kaufman, Shoaib Kamil, and Maryam Mehri Dehnavi. 2020. NASOQ: Numerically Accurate Sparsity-Oriented QP Solver.ACM Transactions on Graphics39, 4 (Aug. 2020)...

arXiv 2008
[2025]

doi:10.1145/3731179 2026-02-03 02:00

Adaptive Algebraic Reuse of Reordering in Cholesky Factorizations with Dynamic Sparsity Patterns.ACM Transactions on Graphics44, 4 (July 2025), 119:1– 119:17. doi:10.1145/3731179 2026-02-03 02:00. Page 8 of 1–10. Fast Sparse Matrix Permutation for Mesh-Based Direct Solvers•111:9 104 105 106 Number of mesh vertices 0 2 4Speedup vs. AMD Mean: 1.71x Min: 0.1...

doi:10.1145/3731179 2025

Pith tools

Reviewed August 3, 2026 · model on record in the stance chip above.

[1] [2008]

Algorithm 887: CHOLMOD, Supernodal Sparse Cholesky Factorization and Update/Downdate.ACM Trans. Math. Softw.35, 3, Article 22 (Oct. 2008), 14 pages. doi:10.1145/1391989.1391995 Kazem Cheshmi, Danny M. Kaufman, Shoaib Kamil, and Maryam Mehri Dehnavi. 2020. NASOQ: Numerically Accurate Sparsity-Oriented QP Solver.ACM Transactions on Graphics39, 4 (Aug. 2020)...

arXiv 2008

[2] [2025]

doi:10.1145/3731179 2026-02-03 02:00

Adaptive Algebraic Reuse of Reordering in Cholesky Factorizations with Dynamic Sparsity Patterns.ACM Transactions on Graphics44, 4 (July 2025), 119:1– 119:17. doi:10.1145/3731179 2026-02-03 02:00. Page 8 of 1–10. Fast Sparse Matrix Permutation for Mesh-Based Direct Solvers•111:9 104 105 106 Number of mesh vertices 0 2 4Speedup vs. AMD Mean: 1.71x Min: 0.1...

doi:10.1145/3731179 2025