PaperlayerView exact revision on arXiv

Source-linked AI summary

Accelerating Branch MPC with Two-Level Parallel Direct Solves on GPUs

Fenglong Song, Luyao Zhang, Liang Wu, Ján Drgoňa, Colin N. Jones

arXiv:2609.29977v1eess.SYcs.SYmath.OC

TL;DR

Branch MPC plans for multiple future trajectories, but computation grows with scenario count and horizon length. This paper presents a GPU direct solver that parallelizes across scenarios and along the horizon for root-coupled branch systems. It reports speedups over cuDSS and PARDISO, with the strongest stated maxima reaching 27.6× for factorization and 15.8× for triangular solves over PARDISO.

  • Problem

    Branch MPC optimizes multiple possible futures, but its computational cost grows with both the number of trajectories and prediction horizon.

    Source: §I · p. 1

  • Method

    A GPU direct linear-system solver combines scenario-level and horizon-level parallelism for branch MPC systems with a shared root and independent tails.

    Source: §I · pp. 1–2

  • Results

    The solver was faster than cuDSS at every tested point (2.0–6.0× factorization; 1.7–3.5× triangular solves), while PARDISO speedups reached 27.6× and 15.8×, respectively.

    Source: B. Results · pp. 6–7

  • Takeaways & Limitations

    The solver can serve as a reusable backend for multiple optimization methods whose reduced systems have the required symmetric positive-definite structure.

    Source: §I · p. 2

    The solver targets branch MPC systems with one shared root and independent tails, and applies to reduced systems with the required symmetric positive-definite structure.

    Source: §I · p. 1

Abstract

from arXiv · show

Branch model predictive control optimizes multiple future trajectories coupled through shared decisions, with computational demands increasing as the number of scenarios and prediction horizon grow. We present a GPU-accelerated direct linear solver for branch MPC formulations in which all trajectories share a single root decision node and evolve independently thereafter. By operating at the linear-algebra level, the solver provides a reusable backend for multiple optimization algorithms whose reduced systems have the required symmetric positive-definite structure. The solver exploits two levels of parallelism: across scenarios and along each prediction horizon. A tailored variable ordering enables horizon-parallel Cholesky factorization while preserving a single root-tail coupling block per scenario in the factor. Numerical experiments demonstrate substantial speedups over state-of-the-art sparse direct solvers, achieving factorization speedups of up to 6.0$\times$ over cuDSS and 27.6$\times$ over eight-thread PARDISO, with triangular solve speedups of up to 3.5$\times$ and 15.8$\times$, respectively.

I. INTRODUCTION

Branch MPC plans for multiple possible futures, but its computational cost grows with both the number of trajectories and the prediction horizon. This paper develops a GPU linear-algebra solver for root-coupled branch systems and a recursive ordering that exposes parallelism across scenarios and time.

I. INTRODUCTION

  • Branch MPC addresses uncertainty by optimizing trajectories that can diverge as information becomes available, but computational cost grows with trajectory count and prediction horizon.The paper considers trajectories that share one root decision node and evolve independently afterward.Source: §I · p. 1
  • The approach targets differently named MPC formulations, including stochastic, robust, and risk-averse MPC, that produce linear systems with the same tree-induced sparsity.The paper also lists contingency MPC for automated driving and branch MPC for interactive, multi-modal motion planning.Source: §I · p. 1
  • Existing branch-parallel methods retain sequential temporal recursions, while associative scans and hierarchical factorizations offer horizon parallelism.The cited temporal methods can reduce time complexity from O(N) to O(log N) with sufficient resources.Source: §I · p. 1
  • The proposed direct factorization and triangular solves exploit scenario- and horizon-level parallelism at the linear-algebra level, serving optimization methods with the required system structure.The paper names IPMs and ADMM as examples of methods that can use the backend when their linear systems have the target structure.Source: §I · p. 2
  • Numerical experiments across scenario counts and prediction horizons show substantial speedups over PARDISO and cuDSS.The introduction identifies these as state-of-the-art general sparse direct solvers.Source: §I · p. 2

II. PRELIMINARIES

  • Symmetric permutation reorders the variables equivalently, while recursively selecting every other block creates uncoupled groups for parallel factorization and solves.Each recursive step reduces the remaining blocks by approximately half, enabling O(log N) time complexity with sufficient parallel resources.Source: §II · p. 2

A. Problem Formulation

The formulation separates one shared root block from scenario-dependent tails. Its resulting linear system has block-tridiagonal tails coupled through the root.

A. Problem Formulation

  • The formulation separates the shared root block x0 from the scenario-dependent tail variables x^i_1,...,x^i_N, which contain states and inputs at each stage.At the last stage, inputs do not appear.Source: A. Problem Formulation · p. 3, B. Structure of the Linear System · p. 3
  • The shared root and tail variables are assumed to have unified size; the smaller last-stage block can be padded to satisfy this assumption.The last stage is smaller because it has no input.Source: A. Problem Formulation · p. 3

B. Structure of the Linear System

  • The optimization methods repeatedly solve KΔx = r, where K is SPD and Δx is the unknown.The stated formulation introduces repeated linear-system solutions for these methods.Source: B. Structure of the Linear System · p. 3
  • The system groups scenario tails together and places the root last, with the right-hand side partitioned correspondingly.This ordering is introduced before specifying the root and tail blocks.Source: B. Structure of the Linear System · p. 3
  • Each tail has a block-tridiagonal matrix, and different tails have no direct coupling; they interact only through the root block.The root block is associated with the shared root variables, while each coupling matrix links a tail to the root.Source: B. Structure of the Linear System · p. 3
  • The root couples directly only to the first temporal block of each tail, so each root-tail coupling matrix has just one nonzero block.That block is at the beginning of the tail.Source: B. Structure of the Linear System · p. 3
  • Together, the block-tridiagonal tails and their individual root couplings give the coefficient matrix a block-diagonal-arrow structure.The paper develops a direct solver tailored to this structure.Source: B. Structure of the Linear System · p. 3

IV. TWO-LEVEL PARALLEL DIRECT SOLUTION

The direct solution uses scenario-level and horizon-level parallelism, with a block permutation that exposes parallel elimination while keeping each tail’s root coupling confined to one block.

IV. TWO-LEVEL PARALLEL DIRECT SOLUTION

  • Scenario tails are processed concurrently, while a tailored permutation exposes parallel factorization and triangular solves along each horizon; tails synchronize through the root block.The two levels are scenario-level parallelism and horizon-level parallelism.Source: §IV · p. 3

A. Block Permutation for Two-Level Parallelism

  • The complete-system permutation reorders the variables, unknown, and right-hand side together, making factorization equivalent to using that variable order.Source: A. Block Permutation for Two-Level… · p. 3
  • Within each reversed tail, the recursive permutation is applied to all but the root-coupled block, which is kept last.The same permutation is used for tails with equal horizon lengths and block dimensions.Source: A. Block Permutation for Two-Level… · p. 3
  • Applying the permutation independently to each tail retains scenario independence and allows horizon parallelism while leaving a single nonzero block in each permuted root coupling.Source: A. Block Permutation for Two-Level… · p. 4
  • In the Figure 2 example, reversing variables avoids fill but leaves sequential horizon factorization, while the proposed ordering enables horizon-parallel elimination and confines root coupling to one factor block.The illustration uses M = 2 scenarios, N = 9, and block sizes of 3.Source: A. Block Permutation for Two-Level… · p. 4, Figure 2 · p. 4

B. Parallel Cholesky Factorization

The Cholesky factorization handles each tail independently before combining their Schur contributions at the root.

1) Tail factorization:

Tail factorization uses a horizon-parallel method, and the coupling solve exploits the sparsity of the coupling matrix.

  • Each permuted tail is factorized using the method in [18, Algorithm 4] to achieve horizon-level parallelization.Source: 1) Tail factorization · p. 5
  • The coupling solve uses the fact that the transpose of the permuted coupling matrix has nonzeros only in its last block row.Source: 2) Tail coupling solve · p. 5

2) Tail coupling solve:

The tail coupling solve reduces its computation by exploiting the structure of the coupling factor and the Cholesky factor.

3) Local Schur update:

The section describes a two-level direct solve whose tail work is parallelized across scenarios and whose root operations combine the tail contributions. It also reports the overall complexity of the factorization and triangular solve.

3) Local Schur update:

  • Each tail matrix F_i has one nonzero block, allowing the local update to be simplified.The supplied passage does not include the resulting simplified expression.Source: 3) Local Schur update · p. 5
  • Tail forward substitutions run in parallel, then their results are reduced at the root before root forward substitution.The algorithm labels the root reduction as O(n^2 b + n b log M).Source: 3) Local Schur update · p. 5
  • After the root backward substitution, the tail solutions are assembled and returned.The algorithm labels root backward substitution as O(n^2) and identifies solution assembly.Source: 3) Local Schur update · p. 5

C. Parallel Triangular Solve

The triangular solve combines independent tail solves with root-level reduction and propagation. Tail operations run concurrently across scenarios and use horizon-parallel solves within each tail.

C. Parallel Triangular Solve

  • The solve begins with parallel tail forward substitutions, reduces their contributions at the root, and then solves the root system.These stages use the permuted right-hand side and tail factors.Source: C. Parallel Triangular Solve · p. 5
  • After root backward substitution, the root solution is propagated to the tails, followed by parallel tail backward substitutions.Root reduction and propagation use only the last block of each tail solution because F_i has only one nonzero block.Source: C. Parallel Triangular Solve · p. 5
  • Independent tail stages run concurrently across scenarios, while each tail's forward and backward solves use a horizon-parallel algorithm.The stated overall complexity is O(n^2 b log N + n b log M).Source: C. Parallel Triangular Solve · p. 5

V. NUMERICAL EXPERIMENTS

The experiments benchmark a GPU solver for root-coupled branch MPC against PARDISO and cuDSS across scenario counts and horizons. The proposed solver is faster across the tested cuDSS cases, with substantial factorization and triangular-solve speedups over both baselines; the conclusion identifies shared-input formulations as future work.

V. NUMERICAL EXPERIMENTS

  • The implementation uses NVIDIA Warp and socu, with persistent workspaces and CUDA graphs to reduce kernel-launch overhead.Factorization and solve pipelines are captured separately as CUDA graphs.Source: §V · pp. 5–6

B. Results

  • 2.0–6.0× factorization and 1.7–3.5× triangular-solve speedups over cuDSS; factorization is 2.4–27.6× faster than PARDISO, with triangular-solve speedups ranging from approximately 0.4× to 15.8×.The cuDSS comparison remains faster at every tested point; the PARDISO triangular solve is faster on the smallest problem.Source: Figure 3 · p. 6, B. Results · pp. 6–7

A. Experiment Settings

  • The solver was evaluated on a quadrotor scenario grid spanning 20–200 scenarios and 20–160 stages, using double-precision arithmetic.Reported times exclude one-time setup costs and are medians of 50 executions after 10 warm-up runs.Source: A. Experiment Settings · p. 6, §V · p. 6

B. Results

  • With 100 scenarios, increasing the horizon from 40 to 160 stages raises proposed factorization and solve times by 3.1× and 2.3×, versus 6.0× and 16.6× for PARDISO.At 100 scenarios and 160 stages, factorization remains below 1 ms and triangular solve below 0.35 ms.Source: Figure 4 · p. 7, B. Results · p. 7

VI. CONCLUSION AND FUTURE WORK

  • The solver targets systems with one shared root node; formulations sharing the entire input sequence across scenarios require methods for their different sparsity pattern.The paper identifies this formulation as future work because its horizon-wide input coupling changes the sparsity structure.Source: §VI · p. 7

APPENDIX I LINEAR SYSTEM SPARSITY IN IPM AND ADMM

The appendix compares IPM and ADMM linear-system structures, then proves that the tail-coupling solve has support only in its terminal block.

APPENDIX I LINEAR SYSTEM SPARSITY IN IPM AND ADMM

  • IPM formulations differ in numerical weights, but their reduced systems share the same structural sparsity as ADMM.The appendix considers QP subproblems preserving the coupling structure of the original problem.Source: Appendix I · p. 8
  • Eliminating the second and third block rows of the ADMM system yields a reduced matrix K that is guaranteed symmetric positive definite.K combines P + τI with weighted contributions from A and G.Source: Appendix I · p. 8

APPENDIX II BLOCK-SPARSITY IN TAIL COUPLING SOLVE

  • Theorem 1 shows that the coupling solution F_i has the same terminal-block support as the permuted coupling matrix C_i.The proof uses the nonsingularity of L_i: the first block equation implies F_i's earlier blocks are zero.Source: APPENDIX II BLOCK-SPARSITY IN TAIL… · p. 8
Loading 2609.29977v1…