1 Introduction↩︎

Inverse Kinematics (IK) algorithms compute (collision-free) joint configurations that achieve a desired end-effector pose, a fundamental requirement for numerous robotics applications [@kucuk2004inverse; @isaacs1987controlling; @aristidou2018inverse; @tevatia2000inverse; @angeles1985numerical; @cohn2024constrained]. While analytical solvers (e.g., IKFast [@diankov2010automated], IKBT [@zhang2017ikbt]) are computationally efficient, the inherent redundancy of most robots means general closed-form solutions rarely exist, restricting these solvers to low-degree-of-freedom (DoF) systems with specific topologies.

For general or high-DoF manipulators, numerical approaches provide necessary flexibility but incur higher computational costs. Iterative methods popular in computer graphics, such as FABRIK [@aristidou2011fabrik] and Cyclic Coordinate Descent (CCD) [@wang1991combined; @welman1993inverse; @lander1998making; @kulpa2005fast; @canutescu2003cyclic; @kenwright2012inverse; @chen2019general; @ojha2023singularity], yield high-quality solutions but often lack the strict orientation constraint support required for robotics. Consequently, robotics typically relies on first- [@balestrino1984robust; @wolovich1984computational; @nakamura1986inverse; @wampler1986manipulator; @buss2005selectively; @beeson2015trac] and second-order [@xie2022speedup; @goldenberg1985complete] Jacobian-based optimization, which can readily incorporate additional constraints [@beeson2015trac; @pink; @rakita2018relaxedik; @dantam2021robust; @wang2023rangedik]. While generally effective, these methods are all local methods, and do not guarantee that a solution will be found given an arbitrary initialization [@angeles1985numerical].

To overcome local minima, “global” IK solvers [@huang2012particle; @rokbani2013inverse; @farzan2014parallel; @dai2019global; @starke2020bio] have been proposed. However, they remain computationally prohibitive for real-time applications or restricted to specific robot topologies. Similarly, learning-based neural network approaches [@demers1991learning; @hasan2006adaptive; @ren2020learning; @bensadoun2022neural; @ames2022ikflow] offer fast, parallelizable inference but struggle with precision, often yielding unacceptable errors on the order of centimeters.

Given these challenges and opportunities, parallel computing provides promise through the parallel exploration of multiple candidate configurations that can escape local minima, particularly in cluttered or obstacle-rich environments [@126022; @594541; @gosselin1993parallel; @farzan2014parallel; @sun2024optimization]. In particular GPU-accelerated IK approaches specifically aimed at robotics [@6696709; @126022; @ames2022ikflow; @sundaralingam2023curobo; @danaci2023inverse; @pyroki2025] have been used to produce state-of-the-art results.

Inspired by these recent advances, in this work, we introduce Hybrid Jacobian Coordinate Descent IK (HJCD-IK), a GPU-accelerated hybrid three-phase sampling-based IK solver that supports both self and environment collision avoidance. HJCD-IK combines a novel orientation-aware, greedy coordinate descent initialization scheme, that provides fast, diverse seeds, with a parallel Jacobian-based polishing routine, and a GPU-accelerated collision detection filter.

As previewed in Fig. 1, and detailed in Section 5, this design improves both convergence speed and accuracy over the state-of-the-art methods, both learning- and non-learning-based, consistently finding solutions along the accuracy-latency Pareto frontier and often achieving order-of-magnitude gains. In addition, our method produces not only one, but a broad distribution of high-quality samples, yielding the lowest maximum mean discrepancy compared to baseline approaches. We deploy our solver onto a physical Franka manipulator for real-world demonstrations of real-time, collision-free inverse kinematics (9, 5.4) and release our solver open-source at:

https://github.com/a2r-lab/HJCD-IK.

2 Background↩︎

Forward kinematics (FK) determines the pose of a robot manipulator’s end-effector given a joint configuration. For an \(n\)-DoF manipulator, the FK function can be expressed as follows, where \(\boldsymbol{\theta} = [\theta_1, \theta_2, ..., \theta_n]\) is the joint configuration, and \(\boldsymbol{P}_{ee} \in \text{SE(3)}\) is the end-effector workspace pose: \[\label{eq:fk} \boldsymbol{P}_{ee} = f(\theta_1, \theta_2, ..., \theta_n).\tag{1}\] Inverse kinematics (IK) seeks to determine the inverse, the joint angles \(\boldsymbol{\theta}^*\) that result in a target pose \(\boldsymbol{P}_t \in \text{SE(3)}\): \[\label{eq:ik} \boldsymbol{\theta}^* = f^{-1}(\boldsymbol{P}_t)\tag{2}\] Numerical IK solvers generally turn the IK problem into a constrained optimization problem, minimizing the distance between the end-effector pose \(\boldsymbol{P}_{ee}\) and the target pose \(\boldsymbol{P}_t\), subject to joint limit constraints defined by \(\boldsymbol{\theta}_{min}\), \(\boldsymbol{\theta}_{max}\), and other constraints (e.g., obstacles, task goals), \(g(\boldsymbol{\theta})\):\[\label{eq:opt} \begin{align} \min_{\boldsymbol{\theta}} \quad &||\boldsymbol{P}_{ee}(\boldsymbol{\theta}) - \boldsymbol{P}_t||^2 \\ s.t. & \quad \boldsymbol{\theta}_{min} \leq \boldsymbol{\theta} \leq \boldsymbol{\theta}_{max}, \quad g(\boldsymbol{\theta}) \geq 0. \end{align}\tag{3}\]

Popular methods for solving 3 , e.g., Jacobian-based methods [@10.5555/533662; @balestrino1984robust; @wolovich1984computational; @nakamura1986inverse; @wampler1986manipulator; @buss2005selectively; @beeson2015trac], operate by iteratively linearizing and minimizing the task residual \(r(\boldsymbol{\theta})\), which is the position and orientation error from the target end-effector pose: \[\label{eq:task95residual} r(\boldsymbol{\theta}) = \begin{bmatrix} r_p(\boldsymbol{\theta}) \\ r_o(\boldsymbol{\theta})\end{bmatrix} = \begin{bmatrix} \boldsymbol{P}_t - \boldsymbol{P}_{ee}(\boldsymbol{\theta}) \\ \omega(\boldsymbol{\theta}) \end{bmatrix},\tag{4}\] where \(\omega(\boldsymbol{\theta})\) encodes the quaternion difference: \[\label{eq:quat95diff} \begin{align} q_{err} &= q_t \otimes q_e^{-1} = [w, \mathbf{v}], \\ \omega(\boldsymbol{\theta}) &= \frac{2 \arctan2(\|\mathbf{v}\|, |w|)}{\|\mathbf{v}\|}\mathbf{v}. \end{align}\tag{5}\]

This results in a least-squares update to \(\boldsymbol{\theta}\) [@deo1995overview]. To improve the robustness of this approach near singularities and joint limits, this system is usually solved with a Levenberg–Marquardt (LM) formulation [@sugihara2011solvability]: \[\label{eq:levenberg-marquardt} [J(\boldsymbol{\theta})^\intercal J(\boldsymbol{\theta}) + \lambda D]\Delta \boldsymbol{\theta} = J(\boldsymbol{\theta})^\intercal r(\boldsymbol{\theta}),\tag{6}\] where \(\lambda > 0\) is a damping factor, \(D\) is a positive diagonal matrix that scales joint updates to stabilize the step, and \(J(\boldsymbol{\theta})\) is the manipulator Jacobian at the current joint configuration. We note that the columns of \(J(\boldsymbol{\theta})\) are as follows (for joint \(i\) with world-frame axis \(z_i(\boldsymbol{\theta})\) and position \(P_i(\boldsymbol{\theta})\)): \[\label{eq:jacobian} J(\boldsymbol{\theta}) = \begin{bmatrix} z_i(\boldsymbol{\theta}) \times (\boldsymbol{P}_{ee}(\boldsymbol{\theta}) - \boldsymbol{P}_i(\boldsymbol{\theta}))\\ z_i(\boldsymbol{\theta}) \end{bmatrix}_{i=1...N}.\tag{7}\]

Finally, since multiple solutions often exist for a given target pose, it would be ideal to find the global optimal solution. While this is often computationally infeasible in practice, by leveraging multiple local optima, we can enhance overall IK solver performance for real-time robotic systems. This strategy has led to the current state-of-the-art GPU-accelerated solvers [@sundaralingam2023curobo; @pyroki2025]. For ease of notation later, we refer to such a batch of \(M\) solutions \([\boldsymbol{\theta}_1 \ldots \boldsymbol{\theta}_M]\) as \([\Theta]_M\).

2.1 Cyclic Coordinate Descent (CCD)↩︎

The CCD algorithm (3), is an iterative heuristic search technique designed for solving the position-only IK problem [@wang1991combined], that is, solving for the joint angles \(\boldsymbol{\theta}\) that produce the desired \(x,y,z\) position of the end effector.

In CCD, joints are numbered \(i = 1 \text{ to } n\) starting at the base. As shown in 2, at each step of the algorithm a joint, \(j\) is selected from \(j = n \text{ to } 1\) (from the tip to the root of the kinematic tree). Using joint \(j\)’s rotational axis, \(\vec{r}_j(\boldsymbol{\theta})\), referred to as \(\vec{r}_j\) in the following, the desired rotation angle update, \(\Delta \theta_j\), is computed by projecting the vectors between the position of the current rotational axis and both the target and end-effector positions onto rotational axis plane: \[\label{eq:vec95proj} \begin{align} &\vec{u}_j = \frac{\boldsymbol{P}_{ee}(\boldsymbol{\theta}) - \boldsymbol{P}_{j}(\boldsymbol{\theta})}{||\boldsymbol{P}_{ee}(\boldsymbol{\theta}) - \boldsymbol{P}_{j}(\boldsymbol{\theta})||} &\vec{v}_j = \frac{\boldsymbol{P}_t - \boldsymbol{P}_{j}(\boldsymbol{\theta})}{||\boldsymbol{P}_t - \boldsymbol{P}_{j}(\boldsymbol{\theta})||} \\ &\vec{u}^r_j = \left(\frac{\vec{u}_j \cdot \vec{r}_j}{|\vec{r}_j|^2}\right)\vec{r}_j \quad &\vec{v}^r_j = \left(\frac{\vec{v}_j \cdot \vec{r}_j}{|\vec{r}_j|^2}\right)\vec{r}_j \\ &\vec{u}^\text{~proj}_j = \vec{u}_j - \vec{u}^r_j \quad &\vec{v}^\text{~proj}_j = \vec{v}_j - \vec{v}^r_j. \\ \end{align}\tag{8}\] These projected vectors can then be used to compute \(\Delta \theta_j\): \[\label{eq:theta95proj} \Delta \theta_j = \cos^{-1}\left(\frac{\vec{v}^\text{~proj}_j}{||\vec{v}^\text{~proj}_j||} \cdot \frac{\vec{u}^\text{~proj}_j}{||\vec{u}^\text{~proj}_j||}\right)\tag{9}\] Once all \(n\) joints have been updated, the algorithm measures the distance between the end-effector and target position, often through the \(L_2\) norm, \(||\boldsymbol{P}_{ee}(\boldsymbol{\theta}) - \boldsymbol{P}_t||_2^2\), and exits upon \(\epsilon\)-convergence. However, if the target is out of reach or CCD becomes locked in a singularity [@ojha2023singularity], the algorithm will continue until it reaches a preset iteration limit.

3 Orientation-Aware CCD↩︎

In this section, we develop an orientation-aware CCD scheme leveraging screw theory and manifold optimization to project not only in position space but also in orientation space. To the best of the authors’ knowledge, this is the first time that orientation-based projections have been integrated into a CCD framework for robotics (e.g., [@canutescu2003cyclic] derived custom update steps specifically for protein loop closure).

Given the end-effector quaternion \(q_{ee}(\boldsymbol{\theta})\) and the target quaternion \(q_t\), the orientation error \(q_{err}(\boldsymbol{\theta})\) is defined as the quaternion difference 5 . We decompose the resulting quaternion error \(q_{err} = [w, \boldsymbol{v}]\) into angle-axis form. The rotation angle \(\phi\) is obtained from the scalar part of the quaternion error, while the corresponding axis of rotation is obtained by solving for its vector component [@lechuga2022iterative]: \[\label{eq:angle95axis} \begin{align} &\phi = 2 \times \arccos (w)\\ &[w, \boldsymbol{v}] = \cos\left(\frac{\phi}{2}\right) + \hat{a}\sin\left(\frac{\phi}{2}\right) \;\; \Rightarrow \;\; \hat{a} = \frac{\boldsymbol{v}}{\sin\left(\frac{\phi}{2}\right)} \end{align}\tag{10}\] The angle-axis pair \((\phi, \hat{a})\) defines the orientation update.

For each joint \(j\) with axis \(\vec{r}_j(\boldsymbol{\theta})\), again referred to as \(\vec{r}_j\) for notational simplicity, the update is projected onto the component of \(\hat{a}\) aligned with \(\vec{r}_j\) and scaled with an annealed factor \(\delta(k)\) for each algorithmic iteration \(k\): \[\label{eq:ori95update} \Delta \theta_j = \delta(k) \operatorname{sgn}(\hat{a}\cdot \vec{r}_j)\phi \quad \text{with } \operatorname{sgn}(0)=0.\tag{11}\] By considering this orientation update with the standard positional CCD step, our method performs greedy joint selection in both position and orientation space, enabling rapid convergence with respect to Euclidean and angular error. As such, this scheme enables a CCD-style algorithm to solve the full IK problem as shown in 2 .

4 The HJCD-IK Algorithm↩︎

In this section we present HJCD-IK (4, 5), a GPU-accelerated, sampling-based hybrid solver that couples an orientation-aware greedy coordinate-descent (PO-CCD) initializer with a parallel Jacobian-based polishing stage (PJ-IK), and leverages recent advances in GPU-accelerated collision detection [@huang2026prrtc] to ensure that the returned solution is feasible. We present the full algorithm in 5 and detail its sub-components in later sections and algorithms.

At the highest level, the first stage evaluates hundreds to low-thousands of orientation-aware CCD seeds in parallel (4.1) and returns \(M\) seeds (5 Line 1). These are collected and ranked by their residuals, of which the top \(K\) candidates are then duplicated with small perturbations to reduce the risk of converging to local minima, producing a batch of \(B\) seeds (5 Lines 3-7). These are further refined by a batched Jacobian solver (4.2) in stage 2 to return an optimized batch of \(B^*\) configurations (5 Line 8). These are then filtered for collision avoidance (5 Line 9) to produce the batch \(CF\), and then scored, returning the final optimized configuration \(\boldsymbol{\theta}^*\) (5 Line 10) in stage 3.

We note that alone, orientation-aware CCD cannot converge quickly to high-precision solutions, and batched Jacobian solvers are prone to getting trapped in spurious local minima. However, when combined, orientation-aware CCD generates a diverse set of coarse solutions from many randomized initializations, providing excellent seeds for the Jacobian-based polishing stage, which can then be filtered for collision avoidance. This design enables fast high-quality batch generation and fast post-hoc feasibility checks. We also note that HJCD-IK is explicitly co-designed for GPUs, with kernels that exploit block-, warp-, and thread-level parallelism to maximize throughput while minimizing synchronization. This algorithm-hardware-software co-design approach enables real-time performance and, as shown in 5, produces IK solutions that are both faster and more accurate than state-of-the-art alternatives.

4.1 Massively-Parallel Orientation-Aware CCD (PO-CCD)↩︎

Our co-designed, orientation-aware, greedy CCD stage (6) efficiently harnesses GPU parallelism by distributing orientation-aware CCD computations across blocks, warps, and threads, enabling the simultaneous processing of hundreds to low-thousands of sampled configuration seeds. At the highest level this is done by launching \(M\) parallel blocks that run asynchronously, one for each sample. Within these samples synchronizations are kept to a minimum, and all underlying linear algebra is parallelized across threads. As shown in 6 Lines 2-4, the algorithm begins with uniform sampling of \(M\) initial joint configuration seeds within the joint limits of each joint.

We then enter the main orientation-aware CCD loop (6 Lines 5-14) where each block uses two warps of threads for each joint. This ensures that each pair of warps is responsible for performing an update for a single joint of a single IK problem, one for orientation and one for position. At every iteration, these warps calculate coordinate descent updates while holding all other values constant, producing for each joint, in parallel, the orientation and position task residual (6 Lines 6-8).

To determine the best update for each step, these residuals are stored in shared memory, allowing for fast, low-overhead comparisons across each block in which we greedily select both the best orientation and position update and apply it (6 Lines 9-10). If both updates target the same joint, we choose the larger of the two. To reduce the chance of this process jumping to a local minimum, we require a minimum threshold for improvement, \(\gamma\), along either the position or orientation task space. If this update fails, we instead randomly perturb the joint configuration. This greedy joint selection process allows us to skip joint updates that will result in no progress (e.g., when the normal axis of the rotation plane is close to that of the vector made between the joint and the target position) and reduce total computation time (6 Lines 11-13). Once a seed satisfies the position and orientation error thresholds (6 Line 14), the parallel loop is broken and all samples are returned.

Because this initializer is very fast at finding reasonably good solutions for large scales of hundreds to thousands of seeds, but can struggle to satisfy extremely tight tolerances, it terminates at a coarse tolerance with only a modest number of iterations, deferring fine accuracy to the downstream Jacobian-based refinement.

4.2 Parallel Jacobian IK↩︎

Stage 2 of our algorithm polishes the returned high-quality batch of refined seeds \(B\) into an optimized batch for collision filtering \(B^*\) through parallel Jacobian IK (PJ-IK) as shown in 7. Here, each block handles one seed and parallelizes all underlying linear algebra across the warps of threads within that block. As Jacobian calculations are more expensive than CCD calculations, this reduced batch size of tens to low-hundreds of refined seeds \(B\) also scales well on modern GPU hardware. To ensure that these calculations are done maximally efficiently, we leverage and extend the GRiD library [@plancher2021accelerating; @plancher2022grid] for all Jacobian calculations, which is also designed for overall block-level parallelism and underlying thread-level linear algebra parallelism. In particular, we add support for arbitrary end-effectors to enable real-world robot hardware validation.

For each candidate \(b\), we solve 2 by first minimizing a weighted task residual, where \(W(\boldsymbol{\theta})\) is a diagonal weighting matrix that normalizes the rows of the Jacobian and adaptively scales the translational and rotational components of the residual. Adapting this to the LM formulation, 6 , and dropping the dependence of all variables on \(\boldsymbol{\theta}\) for readability, we solve the following, where \(D= \text{diag}(J^\intercal J)\): \[\label{eq:normal-eqs} (J^\intercal W^\intercal J + \lambda D) \Delta \boldsymbol{\theta} = -J^\intercal W^\intercal W r,\tag{12}\] The interpolation between the gradient and Gauss-Newton step via the \(\lambda D\) damping term provides additional robustness near singularities and joint limits. To further stabilize convergence, each update \(\Delta\boldsymbol{\theta}\) is constrained by a trust region of radius \(R\), preventing aggressive steps far from the solution. Candidate steps are also validated through a backtracking line-search for the scaling factor \(\alpha\), accepting the first iterate that reduces the task residual. This LM step is shown in 7 Lines 3-8. We note that in general, such a line search will solve for \(\alpha^*\) when given an update, \(\Delta \boldsymbol{\theta}\), original value, \(\boldsymbol{\theta}\), and line-search range \(\mathcal{A} = [1, \frac{1}{\beta}, \ldots, \frac{1}{\beta^A}]\), under a cost function, \(r(\cdot)\): \[\label{eq:line-search} \begin{align} \alpha^* = \min \{ \alpha_i \mid r(\boldsymbol{\theta} + \alpha_i \Delta\boldsymbol{\theta}) < r(\boldsymbol{\theta}), \alpha \in \mathcal{A} \} \text{ or } 0\\ \end{align}\tag{13}\]

If no scaled LM step can be found, we resort to two fallback strategies [@martinez2013towards]. The first, a dogleg method [@mizutani1999powell], shown in 7 lines 9–11, constructs a step \(\Delta\boldsymbol{\theta}(\tau)\) that interpolates between the steepest gradient descent step and the Gauss-Newton step by selecting a value \(\tau \in [0, 1]\) such that the update lies within the trust-region \(R\) [@rosen2014rise]: \[\label{eq:dogleg} \begin{align} &\Delta\boldsymbol{\theta}_{GD} = -\alpha J^\intercal r, \quad \quad \Delta\boldsymbol{\theta}_{GN} = -(J^\intercal J)^{-1} J^\intercal r,\\ &\Delta\boldsymbol{\theta}(\tau) = \tau\Delta \boldsymbol{\theta}_{GD} + (1-\tau) \Delta\boldsymbol{\theta}_{GN},\\ &||\Delta\boldsymbol{\theta}(\tau)||_2^2 \leq R^2. \end{align}\tag{14}\] If the dogleg also fails, a single-coordinate line search, shown in 7 Lines 12-14, is attempted by updating the single joint with largest gradient magnitude according to the norm of the weighted residual [@wright2015coordinate]: \[\label{eq:single-coord} \Delta\theta = \max\limits_i \frac{\partial}{\partial \theta_i} \frac{1}{2}||W(\boldsymbol{\theta})r(\boldsymbol{\theta})||^2.\tag{15}\] If all of these fail, a random perturbation is made to avoid local minima (7 Line 15). This is repeated until convergence or maximum iterations (7 Line 16). These combined refinement stages are able to consistently return precise solutions from the coarse candidates, achieving sub-millimeter positional and sub-degree rotational accuracy.

4.3 GPU-Parallel Collision Checking↩︎

The final algorithm stage is a GPU-parallel collision avoidance filter. Here we leverage similar block-, warp-, and thread-optimized approaches as with our prior two stages, adapted from [@huang2026prrtc]. Offline, we use the foam tool [@coumar2025foam] to generate bounding spheres for the robot geometry. Online, we check for collisions against those spheres through a two-stage hierarchical check: a low-resolution pass flags potential collisions, followed by a high-resolution pass on flagged links for precision. To maximize throughput, groups of four threads execute 4x4 matrix operations for forward kinematics, while warp primitives enable low-latency “early exits” upon collision detection. This architecture enables fast queries that scale efficiently across the tens to low hundreds batch sizes common for \(B^*\). We refer the reader to [@huang2026prrtc] for more details on the underlying kernel design.

5 Results↩︎

We benchmark HJCD-IK against three other GPU-parallel IK solvers: cuRobo [@sundaralingam2023curobo], IKFlow [@ames2022ikflow], and PyRoki [@pyroki2025]. These solvers represent different state-of-the-art approaches for solving IK problems based on: GPU-accelerated optimization, generative modeling with normalizing flows, and differentiable JAX-based computation, respectively. All results were collected using a laptop with an Intel Core i7-14700HX CPU (20 core, 2.1 GHz base), an NVIDIA GeForce RTX 4060, WSL Ubuntu 24.04, and CUDA 12.5. We benchmark in simulation the 7-DoF Franka Panda [@haddadin2022franka] and Fetch [@wise2016fetch] arms. We evaluate open-world poses using joint configuration samples from a Halton Sequence to generate 100 feasible target poses. We evaluate collision avoidance using end-effector poses from the box_panda scene in MotionBenchMaker [@chamzas2021motionbenchmaker]. Hardware experiments use a 7-DoF Franka Research 3 [@FrankaResearch3].

5.1 Batch-Size Scalability↩︎

We first present a scalability study computing open-world (1) and collision-free ([fig:meta_colfree], 2) motions with sample batch sizes of \(B\in \{1, 10, 100, 1000, 2000\}\). Across both robots and scenarios, HJCD-IK benefits from increased batch size, yielding rapid error reductions with only modest per-target latency increases. This results in our approach (shown in orange) outperforming state-of-the-art baselines in terms of latency across all batch sizes, while also generally surpassing all baselines in accuracy, remaining on or near the accuracy-latency Pareto frontier.

In open-world benchmarks, increasing batch size reduces Panda position and orientation errors by three to four orders of magnitude, respectively, with Fetch showing even greater improvements. Crucially, HJCD-IK incurs minimal latency penalties with solve times increasing only slightly (4.04ms to 4.37ms for Panda; 2.59ms to 2.73ms for Fetch), yielding order-of-magnitude speedups over IKFlow and \({>3.1}\)x over PyRoki and cuRobo. Collision-free benchmarks on Panda follow this trend, with error reductions of five to six orders of magnitude and latency varying only from 4.19ms to 5.44ms, with order-of-magnitude speedups over PyRoki and cuRobo. Overall, once batch size is 100 or grater, HJCD-IK is the fastest and most accurate solver, and achieves 100% success rate on the challenging box_panda task.

4pt

4pt

4pt

5.2 DoF Scalability↩︎

We assess scalability with respect to manipulator complexity by fixing the batch size to \(B = 1000\) and evaluating 7-, 12-, 18-, and 24-DoF Panda Arm variants⁴ using HJCD-IK, cuRobo, and PyRoki in an open-world environment. As shown in 8 and 3, HJCD-IK maintains the lowest pose error at every DoF, with position confined to \(1.71\times 10^{-5}\) to \(3.84 \times 10^{-5}\) mm and orientation error to \(4.11\times 10^{-8}\) to \(7.32 \times 10^{-8}\) rad, indicating added redundancy does not inflate error. Notably, as DoF increases, this error improvement grows to multiple orders-of-magnitude over baselines. HJCD-IK also demonstrates immense scalability, incurring minimal latency penalties over increases in DoFs. HJCD-IK is consistently faster while maintaining better accuracy compared to both solvers, presenting speedups ranging from \(3.55\)x–\(5.22\)x over PyRoki and \(2.14\)x–\(4.19\)x over cuRobo. Overall, HJCD-IK remains on the accuracy-latency Pareto frontier across all DoFs, showcasing faster IK solve times while converging to lower errors.

5.3 Solution Space Distribution↩︎

Returning batches of diverse solutions indicates that a batched local solver is doing a good job at sampling across the space of local optima to find good globalizing solutions. To evaluate this, we computed the Maximum Mean Discrepancy (MMD) between solver joint configurations and an approximate ground-truth distribution generated via TRAC-IK [@beeson2015trac]. For each of 100 target poses, we compared the 50 best configurations from a batch of 2000 solves from each solver against 50 randomly seeded ground-truth samples. As shown in Table 4, HJCD-IK achieves the lowest MMD (0.02261) and squared MMD (0.00051), indicating that its returned solutions most closely approximate the ground-truth manifold. While IKFlow is competitive in MMD (0.03670), both PyRoki and cuRobo exhibit significantly higher scores, suggesting reduced diversity. The superior performance of HJCD-IK, particularly in \(\text{MMD}^2\), demonstrates its ability to consistently preserve solution redundancy and maintain high-variance coverage even with a minimized return batch. We diagram the impact of this in 1, demonstrating the larger breadth of solutions returned by HJCD-IK over baselines.

5.4 Hardware Deployment↩︎

We demonstrate HJCD-IK’s real-time capabilities, solution diversity, and collision avoidance on a 7-DoF Franka Research 3 manipulator. With every camera frame, we obtain both the target and obstacle poses and give them to HJCD-IK, using a batch size of 1000 to obtain a diverse set of collision-free configurations on the self-motion manifold. To ensure sufficiently smooth motion, we choose to move towards the solution that is nearest to the current configuration.

Through this demonstration, we show the robustness and real-world applicability of HJCD-IK for accurate tracking in the presence of moving obstacles. Additionally, the diversity of solutions offered by HJCD-IK ensures that we are able to generate continuous, collision-free motions, which will be useful in future applications to pathwise IK problems.

6 Conclusion and Future Work↩︎

HJCD-IK provides fast, accurate, collision-free IK solutions for any kinematically redundant manipulator operating in SE(3) by combining a massively parallel, orientation-aware greedy coordinate descent initialization, with a parallel Jacobian-based polishing scheme, and a parallel collision filter. Our experiments show that HJCD-IK outperforms the state-of-the-art methods, remaining on or near the accuracy-latency Pareto frontier across batch sizes and robot DoFs, resulting in order-of-magnitude improvements in latency or accuracy. HJCD-IK also provides lower MMD and MMD\(^2\) scores for batches of solutions, indicating greater solution diversity and coverage. Finally, we deploy our solver for real-world demonstrations of real-time, collision-free inverse kinematics and release our solver open-source.

Looking ahead, there are several promising directions for future work. One of particular note is that while our empirical results demonstrate HJCD-IK’s strong performance, formal analysis of its convergence and optimality remains an important open area of research.

---

1. This material is based upon work supported by the National Science Foundation (under Award 2411369). Any opinions, findings, conclusions, or recommendations expressed in this material are those of the authors and do not necessarily reflect those of the funding organizations.↩︎

2. \(^{1}\)Columbia University, NY, USA. \(^{2}\)Purdue University, IN, USA.↩︎

3. \(^{3,4}\)Barnard College, Columbia University, NY, USA, & Dartmouth College, NH, USA. Correspondence to: plancher@dartmouth.edu↩︎

4. We create the 12-, 18-, and 24-DoF variants using chained copies of the first 6 joints and links of the Panda.↩︎