Try CusADi to use GPUs instead of CPUs ?

[EveCharbie]

I came across CusADi
DOI: 10.1109/LRA.2024.3512254
Repo: https://github.com/se-hwan/cusadi

I didn't try it, but I think it is worth spending some time some day to see how faster we could get with it :)

[mickaelbegon]

I put an agent to test it

[mickaelbegon]

Disclosure: this exploratory investigation and the text below were produced by OpenAI Codex while examining issue #1094 on behalf of Mickael Begon. The GPU benchmark has not yet been run on NVIDIA hardware; the current observations are code-generation and compilation feasibility checks, not measured GPU performance results.

Quick feasibility update, using an actual bioptim function rather than the CusADi toy example:

• Starting from current bioptim master (e8fbc729), I built the basic pendulum OCP with SX and extracted nlp[0].dynamics_func (ForwardDyn: 6 inputs, 1 output, 86 CasADi instructions).
• The current CusADi CUDA generator successfully emitted a CUDA kernel for that function (about 7.1 kB).
• Every CasADi operation used by this bioptim dynamics function is implemented by CusADi; no unsupported operation was encountered.
• A CPU-only Ubuntu GitHub runner installed the CUDA 12.6 compiler and successfully built the generated kernel as an x86-64 libForwardDyn.so for sm_75 (Tesla T4) and sm_86 (Ampere): https://github.com/mickaelbegon/BiorbdOptim/actions/runs/30477727487

This is encouraging for source-level and CUDA compilation compatibility. I have not claimed a speedup yet: the kernel has not been loaded or executed on an NVIDIA device. A real equivalence/performance benchmark therefore still needs a Linux NVIDIA machine or GPU runner.

Test branch, initially based exactly on the same commit as upstream master: https://github.com/mickaelbegon/BiorbdOptim/tree/codex/issue-1094-cusadi-gpu

Once a GPU runner is available, the next test will compare numerical error and batched evaluation time for CasADi CPU versus CusADi, both with and without CPU/GPU transfer costs.

[mickaelbegon]

Follow-up: measured NVIDIA GPU results

Disclosure: this exploratory benchmark, its implementation, and the analysis below were produced by OpenAI Codex while investigating issue #1094 on behalf of Mickael Begon. The measurements were run locally on Mickael's NVIDIA workstation.

This extends the earlier source-generation check with an actual GPU run.

Benchmark machine:

• Ubuntu 20.04.6 LTS, Linux 5.15, x86-64;
• AMD Ryzen 9 5950X: 16 physical cores, 32 hardware threads, one NUMA node, 64 MiB L3 cache;
• 125 GiB system RAM;
• NVIDIA GeForce RTX 3060: 12,288 MiB, 28 streaming multiprocessors, compute capability 8.6;
• NVIDIA driver 570.133.20 and CUDA runtime/toolkit 12.6;
• Python 3.11.15, NumPy 1.26.4, PyTorch 2.8.0+cu126, CasADi 3.7.2, and bioptim 3.5.0;
• bioptim revision 46041b4113ec4e090720e31d812f4ee3a47f7f85 and CusADi revision bffea25849a901277a8f8dc4db16748e629f9053.

The benchmark was expanded from the MadNLP solver-matrix examples to cover six complete OCPs plus one dynamics-only stress case:

• torque-driven pendulum (86 CasADi instructions);
• rigid cube (14);
• joint-acceleration-driven double pendulum (79);
• rigid contact with friction inequalities (958);
• muscle-driven arm (3,725);
• muscle-driven arm with rigid contact (3,781);
• whole-body gait with 34 DoFs, 20 muscles, measured external forces, and a 47,751-instruction dynamics graph.

The complete OCPs were solved with IPOPT on the CPU in fresh processes with Bioptim/CasADi and native-library thread limits set to 4, 8, 16, and 32. Times below are wall-clock seconds (one run per configuration); all 24 solves converged. Dynamics timings are medians of 10 evaluations after 3 warm-up evaluations.

OCP	4 threads	8 threads	16 threads	32 threads
Pendulum	1.887	1.911	1.950	2.812
Cube	0.092	0.107	0.115	0.250
Joint acceleration	1.585	1.712	1.926	4.147
Contact inequalities	8.290	8.499	8.579	10.322
Muscle arm	25.966	25.322	26.653	30.385
Muscle + rigid contact	35.239	34.552	35.779	39.354

Increasing the CPU thread count does not generally improve these IPOPT/MUMPS solves; 32 threads is consistently slower.

For the GPU portion, CUDA does not offer a portable way to reserve an arbitrary number of GPU cores. I therefore varied the number of independent CUDA evaluations (batch sizes 256, 2,048, and 16,384) and the block size (64, 128, and 256 threads). The following speedups include CPU-to-GPU and GPU-to-CPU transfers and compare against the fastest threaded CasADi CPU map among 4/8/16/32 threads:

Dynamics	batch 256	batch 2,048	batch 16,384
Pendulum	0.94x	4.46x	10.36x
Cube	0.69x	3.28x	9.44x
Joint acceleration	0.63x	4.43x	11.19x
Contact inequalities	2.87x	5.32x	6.60x
Muscle arm	3.41x	8.58x	8.56x
Muscle + rigid contact	3.78x	8.15x	8.81x
Whole-body gait	3.88x	8.45x	6.75x

The whole-body reference OCP uses 105 shooting intervals and documents a memory requirement of about 20 GB. Even a two-interval IPOPT trial took more than 90 seconds, so I kept it as a dynamics-only GPU stress test instead of publishing a misleading solver comparison. Its dynamics graph is about 13 times larger than the static muscle-arm graph.

As a separate CPU-only diagnostic, I nevertheless ran this whole-body example with a deliberately small three-interval transcription and fixed IPOPT budgets. These are partial solves, not converged solutions:

IPOPT budget	4 threads	8 threads	16 threads	32 threads
10 iterations	19.652 s	19.853 s	20.258 s	21.155 s
20 iterations	32.005 s	32.129 s	33.190 s	35.116 s

All runs stopped at the requested iteration limit with constraints still violated. They are useful only for comparing the per-iteration CPU cost; they do not establish solution quality.

Numerical validation passed for every generated kernel. The consolidated result contains 252 timing rows (63 GPU configurations x 4 CPU thread counts); maximum observed absolute and relative errors were approximately 2.18e-7 and 1.97e-9, respectively.

Follow-up: local dynamics Jacobians and Hessians are now included. CasADi first constructs the symbolic derivative functions on the CPU, then CusADi generates CUDA kernels that numerically evaluate, for every independent (x, u) sample:

• the continuous dynamics f(x, u);
• its sparse Jacobian J = d(f) / d(x, u);
• the sparse weighted Hessian H = d2(lambda^T f) / d(x, u)2, with one adjoint weight in lambda for each dynamics output.

The Hessian is therefore a local dynamics contribution to a Lagrangian Hessian, not the Hessian of the complete discretized OCP. The six derivative-capable cases were tested with batches 32, 256, and 2,048; blocks of 64, 128, and 256 CUDA threads; and 4, 8, 16, and 32 CPU threads. The table reports the best end-to-end GPU result (including transfers) against the fastest CPU thread configuration for each row:

Model	Jacobian: batch 32 / 256 / 2,048	Hessian: batch 32 / 256 / 2,048
Pendulum	0.43x / 1.66x / 4.75x	0.56x / 1.73x / 4.41x
Cube	0.60x / 0.67x / 3.47x	0.35x / 0.58x / 1.60x
Joint acceleration	0.42x / 1.24x / 4.49x	0.59x / 1.46x / 4.73x
Contact inequalities	1.16x / 2.56x / 4.20x	0.76x / 1.87x / 2.97x
Muscle arm	1.02x / 2.64x / 4.01x	0.90x / 2.04x / 3.63x
Muscle + rigid contact	1.18x / 2.22x / 3.08x	0.85x / 1.81x / 3.10x

All derivative kernels passed numpy.testing.assert_allclose against CasADi CPU values with rtol=1e-9 and atol=1e-9; the largest absolute difference was about 8.14e-10. CusADi's generated CUDA for CasADi OP_SIGN required a small compatibility helper because the current generator emits an undefined CUDA sign() call.

For clarity, what runs on the GPU at this stage is only the batched numerical evaluation of f, J, and H for one continuous dynamics call. Model/OCP construction, symbolic differentiation and CUDA compilation happen on the CPU. The integration scheme, shooting defects, path/boundary constraints, objective terms, assembled NLP Jacobian and Hessian, IPOPT iterations, and sparse KKT/linear solves also remain on the CPU. Consequently, the timings above are not complete OCP iterations or GPU OCP solves.

The full whole-body local Jacobian and weighted Hessian were also generated symbolically, but they contain 2,835,505 and 3,474,646 CasADi instructions and produce CUDA sources of about 224 MB and 275 MB. NVCC had not finished the Jacobian compilation after 10 minutes and was using about 15.5 GB of RAM, so I stopped it rather than report a GPU timing. Directional products (JVP/VJP/HVP) or derivative blocks are the sensible next step for this case.

The next useful layers to add are: GPU evaluation of an entire integration/defect function; GPU kernels for all objective and constraint blocks and their derivatives; solver callbacks that assemble the complete NLP Jacobian/Hessian without round trips; and finally a GPU-capable sparse linear/KKT solver. Only that last integrated stack would justify calling the result a GPU OCP solve.

The practical conclusion is unchanged but now covers derivatives too: CusADi is not generally advantageous for small batches of simple functions, but becomes useful when hundreds or thousands of independent dynamics/derivative evaluations can be grouped. This is promising for parallel simulations or multiple OCPs; exploiting it inside one standard multiple-shooting NLP still requires the integration steps described above.

[mickaelbegon]

Follow-up: complete NLP callbacks and end-to-end IPOPT solves

Disclosure: this implementation, benchmark, and analysis were produced by OpenAI Codex while investigating #1094 on behalf of Mickael Begon. The measurements below were run locally on the same RTX 3060 workstation described in the previous update.

The experiment now goes beyond local dynamics. Bioptim/CasADi constructs the complete discretized NLP, and CasADi nlpsol("ipopt", ...) is connected to interchangeable CPU CasADi or GPU CusADi callbacks for:

• the objective and all NLP constraints, f(x) and g(x);
• the objective gradient through CasADi's grad_f hook;
• the sparse complete constraint Jacobian through jac_g;
• the upper-triangular sparse Hessian of the complete Lagrangian through hess_lag.

IPOPT, globalization, MUMPS, and the sparse KKT solve remain on the CPU. The GPU is used only for the numerical callbacks. Thus this is a genuine hybrid IPOPT OCP solve, not an all-GPU solver.

A CasADi Callback bridge preserves the original input/output sparsities. Before solving, every GPU callback is checked against the corresponding CasADi function with rtol=1e-9 and atol=1e-9.

Exact-Hessian result: cube, 10 shooting intervals

The complete NLP has 97 decision variables, 60 constraints, 270 Jacobian nonzeros, and 94 upper-Hessian nonzeros.

Backend	IPOPT iterations	Solve wall time	Oracle time	Cost	Max constraint violation
CasADi CPU	12	0.048 s	0.0051 s	1117.997115269036	2.22e-15
CusADi GPU	12	0.055 s	0.0180 s	1117.997115269036	2.22e-15

The final decision vectors differ by only 1.78e-15. All four GPU callbacks, including the exact Lagrangian Hessian, were exercised. A cold build of the four CUDA libraries took 12.3 s; a repeated run reused byte-identical generated sources and took about 0.1 s of preparation.

Limited-memory result: pendulum, 10 shooting intervals

IPOPT was run with hessian_approximation=limited-memory, so the GPU supplies f/g, grad_f, and the 290-nonzero sparse Jacobian, while IPOPT maintains its L-BFGS Hessian approximation on the CPU.

Backend	IPOPT iterations	Solve wall time	Oracle time	Cost	Max constraint violation
CasADi CPU	36	0.127 s	0.0427 s	859.5990768962711	2.50e-12
CusADi GPU	36	0.324 s	0.2252 s	859.5990768962713	2.50e-12

The decision vectors differ by 4.94e-9. The callback counts are also identical: 182 primal calls, 76 gradient calls, and 38 Jacobian calls.

The first cold compilation of these three kernels took 190.9 s, dominated by the complete sparse Jacobian. With the new generated-source cache, a repeat reused all libraries and required 0.58 s of preparation.

For comparison, the exact pendulum Hessian contains 386,393 CasADi instructions and generates about 29 MB of CUDA. NVCC had still not completed that one kernel after more than ten minutes, so I stopped the build. This is a compilation-size limitation, not a numerical or API limitation; the smaller cube confirms that the exact-Hessian callback path works.

Current conclusion

The proposed architecture is feasible: the complete CasADi NLP can be solved by IPOPT while f, g, gradients, sparse Jacobians, and exact Hessians are evaluated with CusADi on the GPU, leaving KKT factorization on the CPU.

For these small single OCPs, however, the GPU is not faster. CusADi assigns one CUDA thread to one function instance, and an IPOPT callback for one OCP has batch size one. Kernel launches and host/device transfers dominate. The next performance-oriented step is therefore to expose shooting-node/penalty blocks as a batch within one OCP, or to solve multiple OCPs/multistarts concurrently, while retaining the same CPU KKT architecture.

The implementation and commands are documented on the existing test branch: https://github.com/mickaelbegon/BiorbdOptim/tree/codex/issue-1094-cusadi-gpu

[mickaelbegon]

Clarification on code availability: the new complete-NLP/IPOPT callback implementation and the latest README updates are currently in the local working tree on codex/issue-1094-cusadi-gpu and have not yet been committed or pushed. The linked remote branch is the existing benchmark baseline; the measured results above come from the local changes described in that comment.

[mickaelbegon]

Follow-up with the requested 100-node and richer-dynamics tests.

The implementation is now pushed in mickaelbegon/BiorbdOptim@b6552de5 on codex/issue-1094-cusadi-gpu.

Full hybrid IPOPT solve, 100 shooting intervals

The cube OCP uses 907 decision variables, 600 constraints, 2,700 sparse Jacobian nonzeros, and 904 upper-triangular Hessian nonzeros. The CPU CasADi callbacks and GPU CusADi callbacks both converged in 19 IPOPT iterations to cost 1121.046302678781, with a maximum constraint violation of 3.55e-15. The maximum CPU/GPU decision-vector difference was 1.03e-13.

• CPU callbacks: 0.125 s solve, 0.034 s in numerical oracles
• GPU callbacks including transfers: 0.176 s solve, 0.082 s in numerical oracles
• first CUDA generation/compilation: 274.6 s
• cached preparation on an identical rerun: 0.81 s

So the complete 100-node NLP path works, including exact sparse Hessians and a CPU KKT solve, but a single OCP evaluation is still too small to amortize GPU launch/transfer overhead.

More complex dynamics, OCPs built with 100 intervals

For batches of 100 node-local dynamics evaluations on an RTX 3060, compared with a four-thread CasADi map:

• contact + friction inequalities (958 instructions): 0.298 ms CPU vs 0.204 ms GPU including transfers, 1.46x speedup, max absolute error 1.32e-12
• six-muscle arm + rigid contact (3,781 instructions): 0.951 ms CPU vs 0.325 ms GPU including transfers, 2.93x speedup, max absolute error 2.46e-11

These are continuous-dynamics blocks, not complete IPOPT iterations.

Scaling limit of one monolithic full-NLP kernel

At 100 intervals, the contact-with-friction NLP has 909 variables and 1,200 constraints. Its sparse Jacobian has only 12,400 nonzeros, but the expanded CasADi Jacobian graph reaches 38,833,852 instructions; symbolic construction alone used about 4.15 GB RSS. Generating one giant CUDA translation unit is therefore not a scalable route for richer dynamics.

The practical direction is now clearer: keep IPOPT/KKT on the CPU, but generate reusable per-node integration-defect and derivative blocks (or matrix-free JVP/VJP/HVP kernels), then assemble their sparse outputs on the host. This preserves the multiple-shooting structure and offers enough parallel work for the GPU without materializing a huge monolithic graph.

[mickaelbegon]

Follow-up: the shooting-block Jacobian implementation and cold/warm benchmark are now pushed in mickaelbegon/BiorbdOptim@34808870.

DMS structure

The new benchmark uses OrderingStrategy.TIME_MAJOR. For the contact-with-friction OCP at 100 shooting intervals:

• 909 decision variables and 1,200 constraints
• 100 homogeneous shooting blocks
• 12 constraints, 18 local variables, and 124 Jacobian nonzeros per block
• 12,400 assembled global Jacobian nonzeros
• one 66,513-instruction JVP kernel evaluated for 100 * 18 = 1,800 node/direction pairs

The fixed phase duration at column 0 is the only column shared by every block. All other columns are contiguous: [1,17] for the first interval and [892,908] for the last interval. The assembled sparsity therefore has the expected block-banded DMS layout.

The exact sparse Jacobian is reconstructed on the CPU from the batched JVP outputs. It was checked against a global CasADi JVP: maximum global product error was 1.82e-10, and maximum CPU/GPU coefficient error was 2.95e-10.

Cold and warm timings

Both measurements used a fresh process, the RTX 3060, four CasADi CPU threads, CUDA block size 256, three warmups, and ten timed repetitions. The cold run forced --rebuild; the warm run regenerated the source for byte comparison and reused the compiled library.

Stage	Cold	Warm cache
OCP and block construction	3.220 s	3.084 s
CUDA source generation	0.211 s	0.214 s
CUDA compilation	105.005 s	0.002 s
CUDA library loading	0.119 s	0.115 s
Total GPU preparation	106.372 s	1.374 s
CasADi CPU Jacobian, 4 threads	125.188 ms	128.738 ms
CusADi kernel	11.377 ms	11.299 ms
GPU including transfers	11.568 ms	11.461 ms
Host sparse assembly	0.277 ms	0.314 ms
GPU transfers + assembly	11.845 ms	11.775 ms
Start to first assembled Jacobian	109.604 s	4.470 s

The steady-state Jacobian speedup, including transfers and sparse host assembly, is 10.57x in the cold run and 10.93x in the warm run. The 4.47 s warm-process startup still includes rebuilding the OCP graph, verifying the generated source, importing/initializing PyTorch CUDA, and loading the library; after initialization, the relevant per-NLP-iteration cost is about 11.8 ms.

This removes the 38.8-million-instruction monolithic Jacobian from the GPU path. The next integration step is to expose the assembled sparse values through the CasADi/IPOPT jac_g callback while keeping IPOPT and the KKT solve on the CPU.

[mickaelbegon]

Follow-up with complete OCP wall times, not only Jacobian evaluation times.

The shooting-block Jacobian was connected to CasADi/IPOPT through the jac_g callback. Both backends use the same 100-interval contact OCP, OrderingStrategy.TIME_MAJOR, 909 decision variables, 1,200 constraints, 12,400 Jacobian nonzeros, IPOPT with limited-memory Hessian, MUMPS on the CPU, four CasADi CPU threads, identical bounds/initial guess, and tolerance 1e-6. Only the block-JVP evaluator changes between CasADi CPU and CusADi GPU; the KKT solve remains on the CPU.

Both solves converged in 45 iterations:

• CPU cost: 0.15450797788175952
• GPU cost: 0.15450797788175913
• maximum constraint violation: 4.82e-7
• maximum CPU/GPU decision-vector difference: 3.70e-12
• jac_g calls: 47 for each solve

Complete OCP timings

The cold run forced CUDA recompilation. The warm run used a fresh process but reused the byte-identical compiled library.

Measurement	Cold CPU	Cold GPU	Warm CPU	Warm GPU
OCP/block construction	2.931 s	2.931 s	2.947 s	2.947 s
GPU preparation	—	103.694 s	—	1.349 s
Solver construction	4.074 s	3.892 s	3.871 s	3.881 s
IPOPT solve	9.560 s	2.014 s	9.108 s	1.924 s
Time inside jac_g	8.084 s	0.620 s	7.702 s	0.615 s
Solver construction + solve	13.635 s	5.906 s	12.979 s	5.804 s
Total OCP wall time	16.566 s	112.531 s	15.926 s	10.101 s

Gains

Scope	Cold gain	Warm gain
IPOPT solve only	4.75x	4.73x
Solver construction + solve	2.31x	2.24x
Total including graph build and GPU preparation	0.147x	1.58x

Thus:

• once initialized, the complete IPOPT solve is about 4.7x faster with the GPU block Jacobian;
• including solver construction, it is about 2.2–2.3x faster;
• the first cold GPU OCP is 6.79x slower end to end, because CUDA compilation alone takes about 102.34 s;
• with the compiled library cached, the full OCP from graph construction to solution is 1.58x faster end to end;
• at the measured solve times, roughly 15 solves of the same compiled graph are needed to amortize the one-time CUDA compilation.

These results show why both totals are needed: the block Jacobian materially accelerates the nonlinear solve, but JIT compilation must be cached or amortized before the complete workflow becomes faster.