Pith. sign in

REVIEW 3 major objections 5 minor 4 references

MGPBD: A Multigrid Accelerated Global XPBD Solver

T0 review · 3 major / 5 minor · reviewed 2026-08-15 · deepseek-v4-flash

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

Pith's one-line read Solving XPBD's dual-space system with an AMG-preconditioned CG solver removes the low-frequency stalling that breaks high-resolution deformable simulation.

desk verdict A solid dual-space AMG solver for XPBD; the lazy-setup trick is new and useful, though its validation is thinner than the performance claims. read the letter →

arxiv 2505.13390 v1 pith:DLYUWNSM submitted 2025-05-19 cs.GR

Chunlei Li , Peng Yu , Tiantian Liu , Siyuan Yu , Yuting Xiao , Shuai Li , Aimin Hao , Yang Gao
show 1 more author
Qinping Zhao
This is my paper · ORCID
classification cs.GR
keywords MGPBDXPBDalgebraicmultigriddual-spaceoptimizationdeformableobjectsimulationpreconditionedconjugategradientposition-baseddynamicsGPU
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

Extended Position-Based Dynamics (XPBD) is widely used because it is fast and simple, but its per-constraint Gauss-Seidel updates cannot propagate low-frequency errors, so stiff or finely meshed simulations stall or blow up. This paper claims that solving XPBD's global dual-space linear system with an algebraic multigrid (AMG) preconditioned conjugate-gradient solver removes those low-frequency errors, and that two such iterations eliminate residual spikes that 300 XPBD iterations cannot. To keep the multigrid cost practical, the paper introduces an unsmoothed aggregation step for sparser coarse grids, a lazy setup that reuses prolongators for up to 100 frames, and a near-kernel built from a few Gauss-Seidel sweeps on the homogeneous system. It also concedes that the global solve needs an explicit sparse matrix and introduces extra oscillation that requires damping. If the claims hold, production PBD/XPBD pipelines gain a drop-in solver for stiff, high-resolution cloth, muscle, and soft-body scenes that currently stall or crash.

What carries the argument

The load-bearing object is the global dual-space linear system whose unknowns are the Lagrange multipliers of the constraints, assembled each frame as $A = \nabla C M^{-1}\nabla C^T + \tilde{\alpha}$. The argument runs through three construction choices: unsmoothed aggregation (UA) builds the interpolation operator without the smoothing pass of smoothed aggregation, keeping coarse-grid matrices sparse; the near-kernel components that seed the hierarchy are generated by roughly 20 Gauss-Seidel sweeps on the homogeneous equation $Ax = 0$, repeated six times; and a lazy setup strategy reuses the prolongators for many frames because the matrix structure is fixed when the mesh topology is static and the entry values change slowly. The solving phase wraps one V-cycle inside a preconditioned conjugate-gradient loop, with $\omega$-Jacobi as the smoother for soft bodies and Chebyshev for cloth.

What would settle it

Run MGPBD with setup interval 20 on a scene with large inter-frame rotations or violent collisions, such as a stiff beam twisted 180 degrees within a few frames or cloth smashing into a sphere, and compare the relative dual-residual curve against a run that rebuilds the AMG hierarchy every frame; if the lazy-setup curve flattens or diverges while the fresh-setup curve converges, the reusability premise fails.

Watch

Extended reading notes

Core claim

The paper's central claim is that the stalling of XPBD comes from its diagonal-only local update, which discards the off-diagonal coupling between constraints, and that solving the full Schur-complement system $A = \nabla C M^{-1}\nabla C^T + \tilde{\alpha}$ in the dual space with a multigrid-preconditioned conjugate-gradient loop removes both low- and high-frequency residual errors. The proposed MGPBD solver assembles $A$ each iteration, runs a V-cycle of unsmoothed-aggregation AMG as the preconditioner, and updates positions via $\Delta x = M^{-1}\nabla C^T \Delta \lambda$. Its central evidence is the residual power-spectral-density plot showing that two MGPBD iterations flatten the low- and high-frequency spikes that remain after 300 XPBD iterations, together with convergence curves in which XPBD stalls or crashes at 850K tetrahedra while MGPBD keeps decreasing the relative dual residual.

Load-bearing premise

The load-bearing premise is the paper's Section 4.2 claim that prolongators stay accurate across many frames because the dual-system entries change little unless constraint gradients or the angles between them change substantially.

Editorial extensions

If this is right

  • Under the same time budget, MGPBD keeps decreasing the relative dual residual in high-stiffness, high-resolution scenarios where XPBD stalls; at 850K tetrahedra XPBD diverges while MGPBD stays stable.
  • High stiffness and large time steps become usable together: MGPBD keeps a beam with stiffness $10^{12}$ stiff at time steps of 10, 20, and 30 ms, whereas XPBD shows resolution-dependent softness.
  • The AMG hierarchy amortizes: lazy setup cuts the setup share of runtime from about two-thirds to about 2%, so per-frame cost is dominated by the V-cycle solve rather than hierarchy construction.
  • Time per iteration scales roughly linearly with resolution ($R^2 = 0.9978$), unlike direct solvers, which makes million-element models feasible, including a 1.667M-tet human muscle at 40.6 seconds per frame with 20 iterations per frame.
  • MGPBD can be blended into existing XPBD pipelines: it takes the same inputs and outputs as XPBD and can mix with off-the-shelf implementations for the parts that still need classic XPBD sweeps.

Reading between the lines

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

  • Editorial inference: the lazy-setup result implies the AMG hierarchy could in principle be built once per scene or even once per stiffness regime when the topology is static and motion is mild; the practical limit is the rate of change of constraint gradients, not the frame count itself.
  • Editorial inference: because the system lives in the dual constraint space, the method inherits XPBD's insensitivity to stiffness-to-weight ratio; the same UA-AMG preconditioning could be tested on other dual-space or primal-dual dynamics solvers to see whether the convergence advantage transfers.
  • Editorial inference: the near-kernel technique suggests a general cheap bootstrap for aggregation AMG: run any simple relaxation on the homogeneous equation and use the results as near-nullspace vectors; how few sweeps suffice probably depends on the smoother and material, and that trade-off is not characterized in the paper.
  • The paper's own limitation note, that explicit matrix assembly blocks dynamic topology and the global system adds oscillation that needs damping, points to the clearest next test: whether the convergence gain survives tearing, fracturing, or severe collisions that force frequent hierarchy rebuilds.
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.

Referee Report

3 major / 5 minor

Summary. The paper presents MGPBD, a global dual-space XPBD solver that replaces the per-constraint Gauss-Seidel update of XPBD with a PCG solve preconditioned by an unsmoothed aggregation algebraic multigrid (UA-AMG) hierarchy. The global system is the Schur-complement system in Eq. (3). The contributions are (i) a lazy setup strategy that reuses the AMG prolongators for many frames with only occasional rebuilds, based on the claim that the system matrix changes slowly; (ii) a cheap near-kernel construction obtained by a few Gauss-Seidel sweeps on the homogeneous equation, repeated six times; and (iii) a GPU implementation in CUDA/Taichi. Experiments compare convergence and runtime against XPBD, AMGX, AMGCL, and PARDISO on tetrahedral soft bodies, cloth, muscle, beam, ball, and collision scenes, with reported speedups of two to three orders of magnitude over the third-party solvers.

Significance. If the results hold, MGPBD would be a practical remedy for the well-known stalling of XPBD on high-resolution and high-stiffness problems, and it would be the first AMG-based solver applied in the dual space of XPBD. The claim that two MGPBD iterations remove both low- and high-frequency residual spikes that 300 XPBD iterations cannot (Fig. 2) is strong and falsifiable. The paper ships code, uses standard AMG machinery, and includes comparisons against external libraries, which are real strengths. However, the efficiency claims rest on several empirically validated assumptions, most notably the lazy setup premise, and on comparisons whose reporting is incomplete. The near-kernel simplification is an interesting practical contribution, but its robustness and the fairness of the baseline comparisons need to be established before the headline numbers can be taken at face value.

major comments (3)
  1. [Sec. 4.2, Fig. 5] The lazy setup premise is the load-bearing component of the efficiency claim (setup cost drops from roughly two-thirds of runtime to 2%), yet it is validated only on the 270K-tet bunny squash scenario. The cloth (Fig. 14) and muscle (Figs. 1 and 11) experiments use setup interval 20 without any sensitivity test. For cloth distance constraints, the off-diagonal entry A_ij is proportional to the cosine of the angle between two edges sharing a vertex, so draping and stretching can rotate those angles substantially; for ARAP muscle constraints, matrix entries scale with the deformation gradient, which can change by large factors. Since stale prolongators degrade the Galerkin coarse operators and increase V-cycle counts, the paper should report residual-versus-frame or total-frame-time data for setup intervals 1, 20, and 100 on cloth and muscle cases. Without that, the generality of the 2% setup overhead claim is unsupported.
  2. [Sec. 5.2, Fig. 10] The comparison against AMGX, AMGCL, and PARDISO in Fig. 10(b) reports only per-iteration time (Ours: 0.54s; AMGX: 3.85s; AMGCL: 1.63s; PARDISO: 91.98s). Unless the number of iterations to reach the same residual tolerance is reported, these numbers do not support the 'two or three orders of magnitude improvements' claim. A direct solver such as PARDISO performs a factorization whose cost is amortized over one solve, so comparing its single-call time to one MGPBD iteration is not an equal-footing comparison. The authors should provide total wall-clock time (including setup and all iterations) to reach the same relative dual residual for all solvers.
  3. [Sec. 4.5 and Sec. 4.3] The smoother's omega_opt depends on lambda_min(D^{-1}A), approximated by a user-defined estimate (e.g., 0.1), and the paper states that this approximation 'boosts performance by 24%' in the 850K-tet bunny case. The near-kernel construction uses 20 GS sweeps repeated six times from a random initial vector. These are free parameters, and the smoother is the dominant cost of the V-cycle (85% of an iteration). The paper should provide a sensitivity analysis for the lambda_min estimate and for the near-kernel sweep count/random seed, or explicitly state which values were used in each experiment. Without this, the reproducibility and robustness of the convergence results are not fully established.
minor comments (5)
  1. [General] There are several typos and formatting errors: 'Comparision' in Sec. 5.1, 'interleved' in Sec. 4.1, 'aross' in the Fig. 5 caption, and 'comapred' in the Fig. 8 caption.
  2. [Sec. 4.2] The formulas for A_ii and A_ij contain garbled symbols (e.g., '˝') and undefined notation for m^{-1}_{s_v}; the mathematical typesetting should be fixed.
  3. [Fig. 2] The caption and surrounding text say the red subgraph shows dual residuals before solving and the green subgraphs after solving, but the figure layout is described as 'from left to right' in a way that is hard to follow; please add clear labels directly in the figure.
  4. [Sec. 5.2, Fig. 9] The text says 'Fig. 9 shows the time taken to achieve a 10^{-2} relative dual residual,' but the y-axis is labeled 'Time Per Iteration (ms).' This is misleading: the linear scaling claim (R^2 = 0.9978) applies to per-iteration cost, not end-to-end time, unless the iteration count is shown to be resolution-independent. Please clarify the caption and, if appropriate, also plot total time to the tolerance.
  5. [Sec. 4.3] The justification that six near-kernel vectors are used 'because the solid motion consists of six rigid body modes' is heuristic in the dual space, where the unknowns are Lagrange multipliers rather than positions. Please either provide a dual-space derivation or soften the claim.

Circularity Check

0 steps flagged · score 2.0 of 10

No circular derivation: the AMG hierarchy and near-kernel are built from the system matrix itself, and the central claim is tested against external baselines.

full rationale

The paper's derivation chain is self-contained. The global dual system (Eq. 3) is the standard Schur complement of the XPBD Newton system; the AMG setup (Sec. 4.1) uses standard strength-of-connection filtering, aggregation, and Galerkin projection computed from the matrix A itself, not from a target result. The near-kernel components (Sec. 4.3) are generated by a few sweeps of GS on the homogeneous equation Ax=0, a conventional AMG bootstrap heuristic; no fitted constant is fed back into the derivation as a 'prediction.' The lazy setup strategy (Sec. 4.2) is an empirical approximation validated in Fig. 5, not a circular definition. Performance claims are compared against external baselines (XPBD, AMGX, AMGCL, PARDISO) and multiple scenes. Self-citations to Xian et al. 2019 and Ruan et al. 2024, which include co-author Tiantian Liu, appear only for background positioning and to motivate unsmoothed aggregation; the sparsity advantage is independently demonstrated in Fig. 4. No step reduces by construction to its own input.

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

The central claim rests on the standard XPBD dual-space formulation, the SPD assumption for A, and the temporal-stationarity heuristic for the system matrix; five hand-tuned parameters are introduced. No new physical entities are postulated.

free parameters (5)
  • strength threshold theta_s = 0.1
    User-defined threshold for weak-connection filtering in AMG setup (Sec 4.1); stated as optimal for most softbody cases, but no sensitivity analysis.
  • lazy setup interval = 20 frames
    Chosen as conservative in all experiments; Fig 5 shows little difference up to 100 frames.
  • near-kernel GS sweeps = 20 sweeps, 6 random starts
    Sec 4.3: 20 sweeps of GS on Ax=0 with random initial vector in (0, max|A_ij|), repeated six times; no sensitivity study.
  • lambda_min estimate for omega-Jacobi = 0.1
    Sec 4.5: user-defined estimate of smallest eigenvalue of D^{-1}A; boosts performance by 24% in one case; value is heuristic.
  • relaxation factor omega (line 11) = 0.1 softbody / 0.25 cloth or backtracking
    User-specified or halved when residual increases (Algorithm 1, line 11).
assumptions (4)
  • domain assumption The dual-space XPBD linear system A Delta_lambda = b (Eq. 3) is the correct root system for solving the constraint dynamics.
    Taken from Macklin et al. 2016/2020; MGPBD relies on this formulation.
  • domain assumption A is symmetric positive definite so PCG and the power method apply.
    Needed for PCG; for compliant constraints with positive masses this generally holds, but the paper does not prove it for ARAP constraint gradients.
  • domain assumption Rigid body modes are the relevant near-nullspace for the dual system.
    Sec 4.3 uses six vectors assuming six rigid body modes; for tet-based ARAP this is plausible but not derived for the constraint dual matrix.
  • ad hoc to paper The system matrix A changes slowly enough that prolongators remain effective for up to 100 frames.
    Sec 4.2 lazy setup premise; argued by gradient stability, validated only on the 270K bunny case.

how reviews work

0 comments
Cite this review

Pith. "Pith review of MGPBD: A Multigrid Accelerated Global XPBD Solver." pith.science (2026). https://pith.science/paper/DLYUWNSM

@misc{pith2026250513390,
  author       = {Pith},
  title        = {Pith review of: MGPBD: A Multigrid Accelerated Global XPBD Solver},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/DLYUWNSM}},
  note         = {Machine review of arXiv:2505.13390}
}
read the original abstract

We introduce a novel Unsmoothed Aggregation (UA) Algebraic Multigrid (AMG) method combined with Preconditioned Conjugate Gradient (PCG) to overcome the limitations of Extended Position-Based Dynamics (XPBD) in high-resolution and high-stiffness simulations. While XPBD excels in simulating deformable objects due to its speed and simplicity, its nonlinear Gauss-Seidel (GS) solver often struggles with low-frequency errors, leading to instability and stalling issues, especially in high-resolution, high-stiffness simulations. Our multigrid approach addresses these issues efficiently by leveraging AMG. To reduce the computational overhead of traditional AMG, where prolongator construction can consume up to two-thirds of the runtime, we propose a lazy setup strategy that reuses prolongators across iterations based on matrix structure and physical significance. Furthermore, we introduce a simplified method for constructing near-kernel components by applying a few sweeps of iterative methods to the homogeneous equation, achieving convergence rates comparable to adaptive smoothed aggregation (adaptive-SA) at a lower computational cost. Experimental results demonstrate that our method significantly improves convergence rates and numerical stability, enabling efficient and stable high-resolution simulations of deformable objects.

Figures

Figures reproduced from arXiv: 2505.13390 by the authors.

Figure 1
Figure 1. Human muscle case (1.667M tetrahedrons) simulated via MGPBD. 20 iterations per frame. Performance: 40.6 s/frame. [PITH_FULL_IMAGE:figures/full_fig_p001_1.png] view at source ↗
Figure 2
Figure 2. Power density graph after 2 iterations of MGPBD and 300 iterations [PITH_FULL_IMAGE:figures/full_fig_p003_2.png] view at source ↗
Figure 3
Figure 3. Pipeline of the setup phase of UAAMG. Input: [PITH_FULL_IMAGE:figures/full_fig_p004_3.png] view at source ↗
Figures from the paper (12 more)
Figure 4
Figure 4. Figure 4: Sparsity patterns of SA and UA. Number of nonzeros ( [PITH_FULL_IMAGE:figures/full_fig_p004_4.png]
Figure 5
Figure 5. Figure 5: Residual curves for different setup intervals. There is no signif [PITH_FULL_IMAGE:figures/full_fig_p005_5.png]
Figure 7
Figure 7. Figure 7: Comparison between our improved near nullspace method and the [PITH_FULL_IMAGE:figures/full_fig_p005_7.png]
Figure 8
Figure 8. Figure 8: Relative energy of Bar Twist under various stiffness comapred with primal space method. Relative energy is normalized by dividing it by the energy before the loop. 104 105 106 102 103 Resolution (Tet) Time Per Iteration (ms) [PITH_FULL_IMAGE:figures/full_fig_p006_8.png]
Figure 9
Figure 9. Figure 9: Time per iteration to achieve 10−2 relative dual residual with respect to the resolution. It shows that the time per iteration is linearly scaled up with resolution, where 𝑅 2 = 0.9978 for linear regression. Data are from bunny-squash cases. same time budget for both M…
Figure 10
Figure 10. Figure 10: Performance of AMG Implementation. a) Time distribution of each [PITH_FULL_IMAGE:figures/full_fig_p007_10.png]
Figure 11
Figure 11. Figure 11: Monster. A monster animated with bones. The model has 72K tets and 14K external constraints connecting bones and muscles. The average speed of MGPBD is 2.3 s/frame. The XPBD runs under time budget 2.3 s/frame and failed. The stiffness is 1𝑒9 [PITH_FULL_IMAGE:figures/…
Figure 12
Figure 12. Figure 12: Ball. Different resolutions(1K/22K/99K tetrahedrons from left to right) for a ball falling cases under the same time budget. MGPBD shows the stiff results in high-res cases [PITH_FULL_IMAGE:figures/full_fig_p010_12.png]
Figure 13
Figure 13. Figure 13: Bunny Squash. A bunny recovers from a squashed state under ARAP constraints. The MGPBD and XPBD are compared under the same time budget with different resolutions (0.1s for 5K, 5s for 270K, and 10s for 850K). The 850K XPBD crashes. 𝑑𝑡 = 10ms, stiffness=1𝑒9. XPBD runs …
Figure 14
Figure 14. Figure 14: Cloth. For mass-spring clothes with maxiter = 1e5, the simulation runs until the dual residual converges to 10−4 at various resolutions: #vertices=(𝑁 +1) 2 (N = 64/128/256/512). This case shows that MGPBD (red) has better inextensibility because it is free of stalling…
Figure 16
Figure 16. Figure 16: Collision. Three armadillos collided with 20 cylinders defined in static SDF. A position based collision response is adopted after the MGPBD elasticity loop. 0 100 200 300 400 500 10−6 10−4 10−2 100 Iteration Number Relative Residual 5K MGPBD 12K MGPBD 270K MGPBD 850K…
Figure 17
Figure 17. Figure 17: Residual of Bunny Squash at various resolutions. Relative dual residual versus iteration number (the left subfigure) and time (the right subfigure) in the bunny squash case. Our MGPBD solver avoids stalling, whereas XPBD stalls at all resolutions, and even crashes whe…

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

4 extracted references · 2 canonical work pages

  1. [4]

    In SIGGRAPH Asia

    Subspace-Preconditioned GPU Projective Dynamics with Contact for Cloth Simulation. In SIGGRAPH Asia. Association for Computing Machinery, New York, NY, USA, 1:1–1:12. https://doi.org/10.1145/3610548.3618157 Xuan Li, Minchen Li, and Chenfanfu Jiang. 2022. Energetically consistent inelasticity for optimization time integration. ACM Trans. Graph. 41, 4, Arti...

    arXiv 2020

  2. [2007]

    ACM Trans

    Efficient simulation of inextensible cloth. ACM Trans. Graph. 26, 3 (July 2007), 49–es. doi:10.1145/1276377.1276438 Yuanming Hu. 2020. The Taichi programming language. InACM SIGGRAPH 2020 Courses (Virtual Event, USA) (SIGGRAPH ’20). Association for Computing Machinery, New York, NY, USA, Article 21, 50 pages. doi:10.1145/3388769.3407493 Krishna Murthy Jat...

    arXiv 2007

  3. [2021]

    International Conference on Learning Representations (ICLR) (2021)

    gradSim: Differentiable simulation for system identification and visuomotor control. International Conference on Learning Representations (ICLR) (2021). https: //openreview.net/forum?id=c_E8kFWfhp0 Theodore Kim. 2020. A Finite Element Formulation of Baraff-Witkin Cloth. Computer Graphics Forum 39, 8 (2020), 171–179. doi:10.1111/cgf.14111 arXiv:https://onl...

    work page arXiv 2021

  4. [2023]

    In Proceedings of the 16th ACM SIGGRAPH Conference on Motion, Interaction and Games (Rennes, France) (MIG ’23)

    Primal Extended Position Based Dynamics for Hyperelasticity. In Proceedings of the 16th ACM SIGGRAPH Conference on Motion, Interaction and Games (Rennes, France) (MIG ’23) . Association for Computing Machinery, New York, NY, USA, Article 21, 10 pages. doi:10.1145/3623264.3624437 Yizhou Chen, Yushan Han, Jingyu Chen, Zhan Zhang, Alex Mcadams, and Joseph Te...

    work page arXiv 2020

Pith tools

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