Model Predictive Control for Fluid Human-to-Robot Handovers

Human-robot handover is a fundamental yet challenging task in human-robot interaction and collaboration. Recently, remarkable progressions have been made in human-to-robot handovers of unknown objects by using learning-based grasp generators. However, how to responsively generate smooth motions to take an object from a human is still an open question. Specifically, planning motions that take human comfort into account is not a part of the human-robot handover process in most prior works. In this paper, we propose to generate smooth motions via an efficient model-predictive control (MPC) framework that integrates perception and complex domain-specific constraints into the optimization problem. We introduce a learning-based grasp reachability model to select candidate grasps which maximize the robot's manipulability, giving it more freedom to satisfy these constraints. Finally, we integrate a neural net force/torque classifier that detects contact events from noisy data. We conducted human-to-robot handover experiments on a diverse set of objects with several users (N=4) and performed a systematic evaluation of each module. The study shows that the users preferred our MPC approach over the baseline system by a large margin. More results and videos are available at https://sites.google.com/nvidia.com/mpc-for-handover.

I. INTRODUCTION
Recently, remarkable advances have been made in generic human-to-robot handovers of arbitrary graspable objects [1], [2] thanks to development of learning-based grasp planners for unknown objects [3], [4]. In these works, grasps are selected based on grasp stability or task constraints, and the robot is driven towards the end-effector pose to complete the handover process.
While promising results have been demonstrated using this pipeline, how to select the best grasp and associated motion to execute is still an open challenge. Prior works select grasps based on grasp quality [1], [2], which neither guarantees the reachability and manipulability of the grasp nor the smoothness of the motion generated to reach the grasp. In this paper, we propose to learn a neural network to predict the manipulability of a grasp pose for motion-aware grasp selection. Together with the grasp stability, it enables us to select a stable grasp which is also reachable and manipulable. Additionally, compared with using inverse kinematics to check the reachability of each grasp sequentially, our model can evaluate the reachability and the manipulability of a batch of grasps in parallel, which is more efficient for realtime robotics tasks such as human-to-robot handovers.
Another challenge is how to plan smooth and natural motion for handovers. Given a selected grasp, the robot is Initial State Fig. 1: How to responsively generate smooth motions for handovers is an open question. The pipeline of planning motion after selecting a grasp, such as in [2], does not always guarantee smooth motions (top). We propose an approach to use Model-Predictive Control (MPC) for motion generation with either a single selected grasp or a set of grasps (bottom), which results in a better, more fluid user experience. More results and videos are available at https://sites.google.com/nvidia. com/mpc-for-handover.
usually driven towards the end-effector pose by local (nonglobal) policies such as Riemannian Motion Policies [5] as in [2], [6] or visual servoing [1]. However, natural motions are not just accelerations towards a particular end-effector pose; they consider how to reach there while obeying complex task constraints. In addition, without planned motions, the robot might choose a grasp that leads to slower or less efficient motions. In this work, we employ a highly-parallelized Stochastic Model Predictive Control (MPC) framework [7] for real-time motion planning. This framework also allows us to incorporate a range of complex task-specific constraints proposed in [6] such as preventing hand occlusions.
To bridge the gap between grasp selection and motion planning, we further extend the MPC framework by incorporating the grasp selection procedure into the optimization given a set of candidate grasps. This unified MPC framework enables the system to select the grasp and plan motion simultaneously, while taking grasp reachability and domain specific constraints into account.
Finally, during the physical handover phase, i.e., from the first contact of the receiver's hand on the object till the release of the object, force feedback can be used to decide when to release the object (robot-to-human) or when to grasp the object (human-to-robot) [8]. However, the raw force/torque sensor readings are often noisy and hard to interpret directly [9]. To handle this, we collected a dataset of humans interacting with a moving robot, and trained a neural net to detect when the contact between hand and gripper happens. The detection is then coupled with vision during the physical handover phase to trigger the robot grasp action.

Drop Grasp Wait
Generated grasps Standoffs target grasp target standoff Approach Fig. 2: Different stages of human-to-robot handovers. The robot is idle until it detects a human holding an object (Wait). Then it approaches to a standoff position using MPC (Approach). Once the robot reaches the standoff pose, the robot grasps the object with a blocking policy (Grasp). Finally, the robot places the object to a preset position (Drop).
Our contributions are summarized as follows: 1) We integrate a stochastic MPC framework for efficient motion planning with complex task-specific objectives, and proposed a unified framework for grasp selection and motion planning. 2) We introduce a learned reachability/manipulability model to facilitate motion-aware grasp selection. 3) We collected a dataset and trained a neural net to detect the contact events for improving action triggering in the physical handover phase. 4) We evaluate handovers with the proposed systems on three objects from different handover locations and orientations, as well as dynamically changing the object orientation after the robot starts to move. We further conducted a user evaluation (N = 4) with a diverse set of household objects comparing the proposed system with the baseline.

II. RELATED WORK
Human-robot object handover is essential for assistive and collaborative robots in environments such as retirement homes and factories. Many efforts have been made in humanrobot communication, grasp planning, perception, motion planning, and control to enable natural and fluent handovers [10]. Among them, grasp planning is a crucial component in pre-handover phase to ensure a success physical handover. For example, human gaze was used to determine the contact point on the object to direct a grasp in [11]. Cini et al. [12] chose grasps based on the object constraints and the receiver's task. For human-to-robot handovers, object properties such as shape, function and safety were considered when planning grasp in [13]. Yang et al. [6] decided robot grasps based on how the object is holding by a human giver. Recent advances [1], [2] used deep-learning based grasp planners to generate grasps for unknown objects in human hand. These works selected the grasp to be executed based on the grasp quality or the stability associated with the grasps, which could not guarantee the selected grasps are all manipulable. Our reachability model directly evaluates the robot reachability/manipulability given a grasp and is complementary to the grasp stability and quality metrics.
Model predictive control (MPC) is being increasingly used for reactive motion generation in manipulators [7], [14], [15], [16], [17]. Kuindersma et al. [17] leveraged MPC to generate dynamic motions for humanoid robots, that can navigate varying terrains. Faroni et al. [14] explored leveraging MPC to dynamically slow down the robot when a human is in the workspace of the robot. Minniti et al. [15] used MPC in combination with an observer to handle contact with the environment such as contact during opening of doors. Kramer et al. [16] showed their approach avoiding dynamic obstacles leveraging non-linear programming to solve MPC. Bhardwaj et al. [7] took a different approach from the above methods for MPC, specifically they explore leveraging sampling-based MPC for real-time reactive manipulator control. By leveraging sampling, their approach does not require cost terms to be smooth or differentiable. We leverage [7] in this paper to solve MPC for human-robot handovers as it does not require 1) extensive engineering effort to tune the many heuristics and 2) the cost terms to be smooth/differentiable.

III. TASK MODEL
Our high-level task model largely follows the one in [2], [6]. As illustrated in Fig. 2, there are four states the robot can be in: (1) waiting for human or object, (2) approaching the object standoff position, (3) grasping the object, and (4) placing the object when it is done.
(1) Wait. During this stage, the robot is idle. It will move to a home position, and not otherwise move around. This step will continue until the robot detects a human holding an object in the workspace.
(2) Approach. During this phase, the robot approaches the object and reaches a standoff pose from the object. This standoff pose is 15 cm from the final grasp on the object. Getting to this standoff position is the role of the MPC system discussed in Sec.V.
(3) Grasp. Once the robot reaches the standoff pose, we begin our grasping policy. Unlike the Approach or Wait policies, this is a blocking policy because as the robot gets closer to the object, body tracking, segmentation, and grasp estimation all become less reliable, as discussed in [2]. After moving the end-effector forward to the grasp position, the robot closes its gripper and retreats to the standoff pose.
We detect missed grasps by observing the distance between the gripper fingers after closing. After a failure grasp, the robot goes back to the Approach phase, and attempts to grasp the object again.
We trained a classifier given raw force/torque signal to detect the contact with the gripper. If the human pushes the object into the robot's hand, the grasp motion will terminate early and the gripper will be closed. The force classifier is described in more detail in Sec. IV-C.
(4) Drop. If the robot has the object, it will place it on the table at a preset position.

IV. PERCEPTION PIPELINE
In this section, we first briefly review the modules we used for grasp generation in Section IV-A. Then we present our reachability model which measures the manipulability of a grasp pose in IV-B. Finally, a contact event detector is introduced in IV-C to improve the robustness of the system when the contact happens (e.g., when humans push/pull the object into/from the robot gripper or when the robot gripper touches the object) or when vision modules fail.

A. Grasp Generation
Our system generates a set of candidate grasps for humanto-robot handovers by following [2]. Specifically, we track the robot pose by DART [18] and the human body by the Azure body tracking SDK at 15 Hz, then segment the hand and the object in hand by a pretrained hand segmentation model at 9 Hz. Given the segmented object point cloud, we use a temporally consistent version of 6-DoF GraspNet [2], [4] to generate grasps at 5 Hz.

B. Grasp Reachability Model
While GraspNet [4] estimates a score for each predicted grasp to measure the grasp stability, not all the generated stable grasps can be realized by the robot. To reduce the set of grasps being considered, we introduce a learned ranking/cost function that predicts the manipulability of grasp pose without an intermediate Inverse Kinematics (IK) step. The manipulability of a grasp pose is a measure of how far the robot is to a singular configuration [19]. This metric is particularly valuable when trying to reach a moving target; a highly manipulable grasp is likely to stay reachable as the object moves. While previous works have developed different methods for encoding reachability of grasps, their approaches use simple heuristics [20], proxy representations [21] or require interpolation of precomputed data [22], [23]. Our approach directly learns the manipulability metric from data to quickly predict scores.
Concretely, our manipulability metric M is defined as |J T J| where J is the Jacobian of the IK of the grasp g. If no IK solution exists, we use the negative twist-distance between the closest achievable end-effector pose and the target: where | · | is the determinant operator, dist twist (g, e) is the twist distance between grasp g and the closest achievable pose of the end-effector pose e.
We further incorporate joint limits of the robot into the manipulability score, we use a joint-limit performance metric derived in previous work [24] to form a weighting matrix. The metric for each joint θ = {θ i : i = 1, 2, · · · , dof}: Note that this metric is minimum (0) in the middle for the joint ranges and maximum (infinity) at the extremes of the joint limits. The weighting matrix is a diagonal matrix with entry on the diagonal corresponding to each joint given as The joint-aware manipulability is then computed as |J T WJ| for grasp poses that have IK solutions. This is used to generate the dataset that is used to train the joint-limit-aware reachability model. To facilitate training, we generated a dataset of grasp poses and the corresponding manipulability scores, and trained an multi-layer perception with the mean squared error loss. To account for IK redundancies during data generation, the manipulability of a given grasp pose is the maximum of the manipulability scores over the IK solutions that realize that grasp. The input of the multilayer perceptron (MLP) is the 6-DOF grasp pose (x, y, z, roll, pitch, yaw) ∈ R 6 and the output is the manipulability score. The median inference time of the trained model is 1.3 milliseconds in experiments.

C. Physical Contact Detection
To better coordinate the timing for the physical handover phase, when the visual perception becomes noisier due to the close proximity between the hand and the robot gripper, we trained a feed-forward neural network to detect the contact event between the hand/object and the robot gripper. The network takes the raw sensor data including joint velocities, efforts, force, and torque for the last T = 5 steps as input, and predicts a probability of a contact event.
Our model is a temporal convolutional neural net which encodes each individual timestep through a two-layer MLP, then performs two 1D convolutions over history, before another MLP (with dropout) predicts a single output indicating if a force event was detected within the window.
We collected a physical contact dataset mimicking the handover procedures: we moved the robot gripper to a random position in the workspace. When the robot was about to reach the target pose, we pushed/pulled on the robot gripper, and pressed a key to label the physical contact. We recorded joint positions, velocities, efforts, forces, and torques during the whole procedure. With T = 5, our dataset consists of 20K samples with about 8% positive samples.
Our model achieved 93.6% accuracy on a held-out nonoverlapping test set of 5k examples collected.

V. MODEL PREDICTIVE CONTROL FOR GRASPING
For human-to-robot handover, the robot needs to move to one of the grasp poses X g∈G , and then grasp the object from the human, as shown in Fig. 2. During this robot motion, we also want to encode heuristics that promote fluid human-robot handovers, specifically: 1) encouraging the robot's gripper to move in straight lines; 2) avoiding collisions between robot and human; 3) having the human hand always in view; and 4) reducing jerk of the robot during motion. As some heuristics are in the joint space and some in Cartesian space of the links, we formulate the problem as a kinematic joint space trajectory optimization problem [25], [26] (Sec. V-A). We then solve the formulated problem in real time leveraging stochastic MPC [7] (Sec. V-B).

A. Trajectory Optimization Formulation
Given the robot's joint position θ 0 and velocityθ 0 at timestep 0, we want to compute the joint accelerationsθ t∈[0,H−1] across the horizon H timesteps that minimizes cost terms C(·) while also satisfying constraints formulated below, where (3) contains the cost terms formulated for humanrobot handover, and (4) lists the manipulability and stop cost terms from Bhardwaj et al. [7] that aid MPC in avoiding local minima and overshooting respectively (refer [7] for more details). Constraints (5-6) are collision avoidance constraints that prevent the robot from colliding with the environment, human, and itself. Eq. (7-8) are euler integration equations to obtain joint position, and joint velocity from joint acceleration. We also have box constraints on the robot's joint position, velocity, and acceleration to satisfy joint limits.
1) Reaching Grasp Pose: Given a target pose X g and the current gripper pose X t = F K(θ t ) computed using the forward kinematics of the robot at the current joint configuration θ t , we define a pose distance metric dist(X t , X g ), where w R g and w R t are the rotation matrices of pose X g and X t respectively. The translation vectors of X g and X t are represented by w d g and w d t respectively. We can reach a given gripper pose by using the above distance metric as the goal cost C g (·) = dist(·). We can also optimize for reaching one pose from a set of grasp poses G by writing the goal cost as the distance between the current gripper pose X t and the closest pose in the goal set, The above formulation of reaching a goal set enables MPC to reason about all the grasp poses inside the optimization while also considering all the other heuristics and constraints. This would enable MPC to move fluidly to another grasp from the current grasp given new constraints (e.g., human hand moves to current grasp region). We discuss results from both these formulations in Sec. VI.
2) Straight Line Cost: We found the robot to take very circular trajectories in the Cartesian space due to the presence of many revolute joints in the manipulator. These circular trajectories are hard to predict by users, especially those with minimal domain knowledge. We hence formulate a cost term that penalizes gripper's linear velocity in directions that

Hand Collision Avoidance Hand Occlusion Avoidance Formulation
Human hand location Robot hand location are not parallel to the vector connecting the current gripper position to the goal positiond g = dg−dt ||dg−dt|| . We compute the gripper's current velocity leveraging the kinematic Jacobian J(θ t ) at the current joint configuration θ t and the current joint velocityθ t asḋ t = J(θ t )θ t and then normalize this vector to get the current gripper motion directiond t . The cost is then written as, This cost becomes zero when both vectors are parallel.
3) Collision Avoidance: During gripper motion, the robot needs to avoid colliding with the table, and also the human hand. In addition, we want to keep the hand from being occluded by the robot during the robot motion. To do so, we represent the table as a cuboid, the human hand as a sphere, and use the line connecting the camera's origin with the human hand position to build a capsule as shown in Fig. 3 (left). We also represent the robot's links with spheres similar to [7] and then compute the signed distance between the robot and the environment via an analytic function S e (θ t ). We additionally avoid robot self collisions by using a neural network trained by [7]. The neural network outputs the signed distance 1 between the two closest links given a joint configuration θ t . We then use this as S r (θ t ) in Eq. (6). As demonstrated in Fig. 3, the robot avoids colliding with human hand and keeps the hand from occlusion.

B. Solving as Stochastic MPC
We solve the optimization problem by leveraging stochastic MPC (also known as sampling-based MPC) as it has shown to work with many cost terms [7]. Stochastic MPC optimizes by sampling action sequences for many particles from a distribution, rolling out these actions, computing the cost incurred by each particle, and then using these costs to update the distribution. By iterating through this process, given large number of particles, stochastic MPC can generate motions at 50-100 Hz which is comparable to gradient based MPC methods while having the benefit of not requiring differentiability of the cost terms or the dynamics [7]. Stochastic MPC cannot handle constraints, hence we setup the constraints in our optimization problem as cost terms with large weights which has shown to be sufficient for satisfying constraints on realistic settings [7], [27]. We also leverage this optimization problem for motion generation in other stages of the handover, e.g., "Drop", where we change the goal cost to be a L2 loss on the joint configuration. Our MPC uses a horizon of 2 seconds while optimizing the motion. All constraints have a weight of 5000. We set α 1 and α 2 in Eq. (9) as 70 and 220 respectively and 30 for the straight line cost weight. by starting with a small value and slowly increasing until the motion reached sufficient precision (a total of 30 minutes was spent to tune these parameters).

VI. SYSTEMATIC EVALUATION
System Setup. We conducted experiments with a Franka Emika Panda arm, and an externally mounted Azure Kinect RGBD camera with 1280 × 720 resolution. The system is distributed on 3 desktop computers with an additional realtime desktop for Franka control. Four NVIDIA RTX 2080 Ti and one 3090 GPUs were used. Protocols and Metrics. We evaluated the system performance with three objects with different shapes and sizes, including a banana from the YCB object [28], a cracker box, and a pepper bottle, as shown in Fig. 5. For all the experiments in this section, we repeated the handover trials until three successful handovers have been recorded. Then we measured the system performance by computing the success rate, approach time, and the total time for a successful handover. We also recorded the velocity, acceleration and the jerk during the handover. Baselines. We considered Yang et al. [2] as the baseline method. In addition to the grasp selection criterion proposed in [2], we use the learned reachability metric R(X g ) described in Sec. IV-B to determine which grasps are reachable and restrict the set to ones that are highly manipulable: where S g and S min is the score for grasp g and the minimum acceptable score, X g and X prev denote the pose of grasp g and the pose of previous chosen grasp respectively. X home denotes the end-effector pose at the home position. ω s , ω prev , ω home , ω R are weights. d(·, ·) is a distance metric with both position and rotation components.
We compared the above baselines with the following variations of the proposed approach: 1) MPC/MPC-R: using MPC for motion planning towards one selected grasp (without/with reachability metric); 2) MPC-GoalSet: using MPC to generate motion given a set of grasps. Noted that reachability metric is inherently embedded in the cost function Eq.(4).

A. Handover Reactivity
We first investigate the performance of our approach by rotating the object along it's standing axis 45 degrees after the robots starts moving. As shown in Fig. 4 (Rotation), all MPC-based systems have lower approach/total time for successful handovers since they can react quickly to the changing of the object orientation. The reachability model promotes lower approach time while maintaining a similar success rate for the baseline method [2], and improves the success rate for the MPC variant. We observe that while the MPC-GoalSet achieves the lowest approach time for successful handovers, it has a reduced success rate especially when handling changes in orientation. We suspect this to be because the number of grasps sufficiently reduces after the filtering by the reachability model ( Sec. IV-B), and at least one grasp is consistent across time.

B. Handover Location
We further evaluate the efficiency of our approach in response to H2R handovers at different locations and different ways of holding objects. We handed over three objects at three locations: left, center, and right. For each location, we handed over the object three times with three ways of holding the object, as shown in Fig. 5. Results are reported in Fig. 4.
The MPC is comparable to Yang et al. [2], where we used MPC for motion planning while keeping the robot velocity (Fig. 6 left) and all the other modules the same. MPC reduces the approach time (11.5s to 7.1s) and total time (13.3s to 7.9s) and improve the success rate from 85.7% to 88.3%.
We also investigate integrating the reachability model with both Yang et al. [2] and MPC (denoted as *-R). In this setting, the reachability metric allows for better grasp selection. It improves the approach time of [2] by two seconds and also improves the success rate of MPC by 7.5%. This is because, with the reachability model, the robot tends to choose more reachable/manipulable grasp poses. Our MPC-GoalSet variant has slightly longer approach time and slightly lower success rate compared to MPC/MPC-R, but it achieves smoother motion with less jerk (see Fig. 6 and Sec. VI-C).

C. Motion Metrics
To better understand the motion pattern, we recorded the robot position, velocity, acceleration and jerk during the above experiments. The statistics are reported in Fig. 6. In general, all our MPC-based approaches have more consistent velocity and accelerations compared with the baseline approach [2]. Our approach has less sudden accelerations, resulting in less jerk. In particular, our MPC-GoalSet achieves the least jerk, as in this formulation MPC optimizes to reach a single grasp from the grasp set, while considering the robot's current position, velocity, and acceleration.

VII. USER EVALUATION
We also conducted a user evaluation to compare the proposed system (MPC-R as Ours) with the prior system Yang et al. [2] to validate that our system enables fluid H2R handovers. We recruited four participants from the lab 2 . Each participant did two rounds of handovers to interact with two systems without knowing the details of the systems. Participants were instructed to hand over ten objects from Household-A objects used in [2] to the robot one at a time. They were asked to rate the systems with a set of Likert scale questions that are commonly asked in prior works [10] immediately after each round of interaction with the system. We additionally asked the users which system is less jerky to evaluate the subjective motion smoothness. We counterbalanced the order in which participants interacted with the two systems to avoid the bias. After completing all interactions, participants were asked to rank the two systems in different dimensions and share their opinions and comments through open-ended questions. Objective Evaluation. We report the approach time and the success rate in Table. I. Our approach was able to take over the object with less time except for scissors and toothpaste, and the overall approach time is shorter than the prior system [2]. The success rate of our system is on par with prior system [2]. 2 We were not able to recruit outside participants due to the ongoing COVID-19 pandemic.
Which system was more natural? Which system had more predictable motion?
Which system was more human-like? Which system was LESS jerky? Which system was LESS aggressive? Which system did you feel safer while working with? Which system did you feel more comfortable while working Overall, which system did you prefer? Subjective Evaluation. Fig. 7 reports user ratings over the two systems. 3 out of 4 users preferred our proposed system. Particularly, they thought our system was more predictable, less jerky, less aggressive, safer, and more comfortable to work with. The main reasons for choosing our systems are "I was able to predict the robot motion better", "it felt less aggressive", and "it did not crash" and "required less motion before picking the object". One user preferred [2], commenting "([2]) had a better approach for some of the items and changed trajectories for less items, thus providing a better experience for me". But this same user also commented "some of these questions are really A=B", which proves the user experience of our system is better or on par with the prior system [2] in these aspects. For questions "which system was more natural" and "human-like", half of the users chose [2] and the other half preferred our system.

VIII. DISCUSSION AND CONCLUSION
We proposed a method for fluid human-to-robot handovers by incorporating more knowledge about the robot motion planning, specifically using MPC to find smooth, consistent trajectories. To facilitate grasp selection, we learned a reachability model to prioritize robot grasps with higher manipulability scores. Our approach was demonstrated more efficient through a systematic evaluation and more favorable by users through a user evaluation compared to prior work. Nevertheless, the reachability model and the grasp generation are independent, which sometimes leads to robot oscillations when the reachability scores and the candidate grasps are in disagreement, especially when using MPC given a set of grasps (MPC-GoalSet). In the future, we plan to better incorporate grasp and motion planning as well as reachability modeling to improve the robustness of the system. Fig. 8 illustrates the structure of our system. We used a Franka Emika Panda arm for all experiments. The Franka control box was connected to a realtime Ubuntu desktop. An Azure Kinect RGBD camera was mounted next to the robot to capture the full scene. The perception, planning and control modules are distributed on three desktop computers. Specifically, the Azure Kinect software, including depth image processing and body tracking, and the camera-robot calibration using DART [18] ran on two NVIDIA RTX 2080 Ti GPUs from the same desktop. The hand-object segmentation, GraspNet [4], and the physical contact detection net ran on another two NVIDIA RTX 2080 Ti GPUs from a different desktop. Finally, the stochastic MPC is computed from an NVIDIA RTX 3090 GPU and the task planning is computed without GPUs. All the desktops are connected through a local area network (LAN) and communicated through the Robot Operating System (ROS).

A. System Diagram
B. Neural Network Architectures 1) Physical Contact Detection Model: Our contact detection model takes a concatenation of joint velocities (7 dim), efforts (7 dim), force (3 dim) and torque (3 dim) for the past T = 5 steps as input, which is a T × 20 tensor. Each individual timestep is first encoded through a twolayer perception. Each layer consists of a 1D convolutions with kernel size one and a ReLU activation (the number of channels is given by 20-256,256-512). Then the model performs two 1D convolutions with kernel size one and 512 channels with ReLU activation and Batch Normalization over history. Finally, another two-layer perception (1D convolutions with kernel size 1 and number of channels as 512-256, 256-1, ReLU activation and dropout) predicts a single output indicating if a force event was detected within the window. The model is trained with the binary cross-entropy loss. Fig. 9 visualizes the recorded force, the ground truth label and our prediction of a test sequence.
2) Reachability Model: The reachability prediction model takes in the 6-dim grasp pose as input, which is passed through five fully-connected layers, each with 64 neurons and ReLU activation. The final output is the predicted reachability value and the model is trained with the meansquared-error loss.

3) Hand Segmentation Model:
The backbone of the hand segmentation model is a Feature Pyramid Network [29] based on ResNet-50 [30] pretrained on the ImageNet [31]. The four convolutions with upsampling operations are used to recover the feature map with the original resolution and a binary segmentation mask is predicted to indicate whether a pixel is hand or background. Please see [2] for details of the training dataset.

C. Failure cases
Graspnet [4] and reachability together do not guarantee consistent grasps. As a result, we do see a number of cases where no grasps are possible -for example, GraspNet might predict only grasps that the reachability model says are very low quality. One problematic case sees a number of generated grasps that are right on the threshold. Due to the noise from one step to the next, the method might be trying to reach very different results, leading to oscillations that are unnerving to a human user, as shown in Fig. 10.


Fast and memory-efficient reconstruction of sparse Poisson data in listmode with non-smooth priors with application to time-of-flight PET

Complete time of flight (TOF) sinograms of state-of-the-art TOF PET scanners have a large memory footprint. Currently, they contain ~4e9 data bins which amount to ~17GB in 32bit floating point precision. Using iterative algorithms to reconstruct such enormous TOF sinograms becomes increasingly challenging due to the memory requirements and the computation time needed to evaluate the forward model for every data bin. This is especially true for more advanced optimization algorithms such as the SPDHG algorithm which allows for the use of non-smooth priors using subsets with guaranteed convergence. SPDHG requires the storage of additional sinograms in memory, which severely limits its application to data sets from state-of-the-art TOF PET systems. Motivated by the generally sparse nature of the TOF sinograms, we propose and analyze a new listmode (LM) extension of the SPDHG algorithm for reconstruction of sparse data following a Poisson distribution. The new algorithm is evaluated based on 2D and 3D simulations, and a real dataset acquired on a recent TOF PET/CT system. The performance of the newly proposed LM SPDHG algorithm is compared against the conventional sinogram SPDHG and the listmode EM-TV algorithm. We show that the speed of convergence of LM-SPDHG is equivalent the original SPDHG using binned data. However, we find that for a TOF PET system with 400ps TOF resolution, LM-SPDHG reduces the required memory from ~56GB to 0.7GB for a short dynamic frame with 1e7 counts and to 12.4GB for a long static acquisition with 5e8 counts. In contrast to SPDHG, the reduced memory requirements of LM-SPDHG enable a pure GPU implementation on state-of-the-art GPUs which will substantially accelerate reconstruction times. This in turn will allow the application of LM-SPDHG in routine clinical practice where short reconstruction times are crucial.

Introduction
A major challenge of image reconstruction in positron emission tomography (PET) is noise suppression since the acquired emission data suffer from high levels of Poisson noise due to limitations in acquisition time, injectable dose and scanner sensitivity. To limit the transfer of the data noise into the image during model-based iterative reconstruction (MBIR), different strategies exist. One possibility is to add a "smoothing" prior to the data fidelity term in the cost function that is being optimized. In general, we can formulate the resulting optimization problem for any imaging system where the acquired data follow a Poisson distribution as where x is the image to be reconstructed, P is the linear forward model, d are the acquired data and s are additive contaminations. m i=1 D i ((P x + s) i ) is the negative Poisson log-likelihood, i is the index of the data bin and m is the total number of data bins. In the specific case of time of flight (TOF) PET, P is the time of flight (TOF) PET forward model including the effects of attenuation, normalization and limited spatial resolution, d are the acquired prompt TOF PET coincidences (the emission sinogram), and s are the estimated random and scattered coincidences. R(K·) is a "smoothing prior" consisting of a generic linear operator K that calculates local differences and a proper, convex, lower-semicontinous function R. The level of regularization is controlled by the nonnegative scalar factor β. A specific example for K would be the gradient operator ∇, e.g. approximated by finite forward differences in the discretized setting. Combining the gradient operator for K with the mixed L2-L1 norm for R leads to the well-known Total Variation (TV) prior [1]. The TV prior, as well as many other advanced smoothing priors aiming for edge-preservation such as e.g. Total Generalized Variation (TGV) [2], Joint T(G)V [3,4] Parallel Level Sets [5,6] or directional Total Variation (DTV) [7], require the use of non-smooth functions for R as instrumental building block. This prevents the use of simple and efficient purely gradient-based optimization algorithms to solve (1).

PDHG and SPDHG for PET reconstruction with non-smooth priors
Using the fact that z → D(z) := m i=1 D i (z i ) and z → βR(z) are convex, lower-semicontinuous functions and thus are equal to their convex biconjugates D * * (z) = sup y z, y − m i=1 D * i (y i ) and (βR) * * (z) = sup w z, w − (βR) * (w), respectively, where D * i and (βR) * are the convex conjugates, and that (βR) * (w) = βR * (w/β), we can rewrite (1) as the saddle point problem introducing the dual variables y and w, and the convex dual of the Poisson log-likelihood given as Under mild assumptions on R, which hold for all of the above-mentioned smoothing priors, Problem (2) is equivalent to (1) and can be solved even for non-smooth priors using the generic primal-dual hybrid gradient (PDHG) algorithm by Chambolle and Pock [8]. PDHG is an iterative algorithm that requires the evaluation of the complete forward and adjoint operator in every update. The usage of the original PDHG algorithm to solve (2) for real-world state-of-the-art TOF PET systems, however, usually results in extremely long computation times, because the evaluation of P and P T for state-of-the-art TOF PET systems is computationally very demanding, and because several hundreds to thousands of updates are needed to obtain reasonable convergence. To overcome this limitation, Chambolle et al. published a stochastic extension of PDHG called SPDHG for saddle point problems that are separable in the dual variable in 2018 [9]. In contrast to PDHG, SPDHG has the advantage that the complete forward and adjoint operator are split into n subsets, and that, in every update, only a random subset of the forward and adjoint operator chosen according to a probability p k has to be evaluated. In [10], Ehrhardt et al. applied SPDHG to 3D non-TOF PET reconstruction with TV-like priors and showed that around 10 complete projections and back projections can be sufficient to obtain reasonable convergence using SPDHG with 252 subsets. Moreover, the authors also demonstrated that preconditioning further accelerates convergence. The resulting SPDHG algorithm to solve (2) is summarized in Algorithm 1, where the proximal operator for a convex function f using the weighted norm · 2 M = M −1 ·, · induced by a symmetric and positive-definite matrix M is defined as For the convex conjugate of the negative Poisson log-likehood D * i , this proximal operator can be calculated point-wise and is given by Algorithm 1 SPDHG for PET reconstruction [10] 1: Initialize x(= 0), y(= 0), (w = 0), x ← proj ≥0 (x − T z) 5: Select i ∈ {1, . . . , n + 1} randomly according to (p i ) i 6: if i ≤ n then 7: : 12: δz ← K T (w + − w) 13: w ← w +

14:
end if 15: z ← z + δz 16: z ← z + (δz/p i ) 17: until stopping criterion fulfilled 18: return x The proximal operator for R * obviously depends on the choice of R but can be also efficiently computed using point-wise operations for many popular choices of R. As mentioned in [10], Algorithm 1 converges if we use the preconditioned step sizes setting T = min k=1,...,n+1 T k pointwise, and choosing 0 < ρ < 1 and γ > 0.

Limitations of PDHG and SPDHG for PET reconstruction
Whilst SPDHG is a big step forward for an efficient solution of (2) in terms of computational speed, SPDHG also comes with two main limitations. First, as discussed in Remark 2 of [10], a potential drawback of SPDHG is that it requires keeping at least one more complete (TOF) sinogram in memory (the dual variable y). Moreover, if the proposed preconditioning is used, a second complete (TOF) sinogram (the sequence of step sizes (S k ) n k=1 ) needs to be stored in memory. Storing one or two extra sinograms in memory of modern computers is not a major problem for static single-bed non-TOF PET data, where sinogram sizes are relatively small. However, for simultaneous multi-bed, dynamic or TOF PET data, the size of the complete data sinograms can be become problematic. This issue gets even more severe when aiming for a complete GPU implementation of Algorithm 1 since the available memory on state-of-the-art GPUs is nowadays much smaller compared to CPUs. As an example, for a TOF PET scanner with 25 cm axial FOV and a TOF resolution of ca. 400 ps, a complete unmashed static TOF sinogram for one bed position has approximately 4.4 · 10 9 data bins, requiring ca. 17 GB of memory in 32 bit floating point precision. Note that with improved TOF resolution and increasing axial field of view (e.g. total body PET scanners), the memory required to store a complete TOF sinogram will continue to substantially increase in the future. Second, PDHG and SPDHG only work with "binned" data. In the case of TOF PET reconstruction, that means that the acquired raw listmode data first need to be binned into TOF sinograms and that P and P T need to be evaluated using sinogram projectors. For most acquisitions with modern TOF PET scanners, this is inefficient both in terms of memory and computational time since the TOF emission data are extremely sparse. In contrast, storage and processing of the data in listmode format (event by event) is usually more efficient.

Sparsity of TOF PET data
Compared to non-TOF PET emission sinograms, TOF PET emission sinograms of most acquisitions with state-of-the-art TOF PET scanners are extremely sparse. This is because every geometrical line of response (LOR) has to be subdivided into several small TOF bins. To achieve sufficient sampling of the TOF information, the number of TOF bins has to be inversely proportional to the TOF resolution of the scanner. Consequently, for a fixed number of acquired prompt coincidences, the sparsity of the sinogram is proportional to the TOF resolution. As an example, for a typical 80 s acquisition of a liver bed position in an FDG scan with an injected dose of 323 MBq acquired 60 min p.i. on a state-of-the-art TOF PET/CT scanner with 20 cm axial FOV, more than 94% of the data (TOF sinogram) bins are empty. An example that demonstrates this extreme sparsity of TOF sinograms is shown in Fig. 1. For shorter frames, as present, e.g., in the early phase of dynamic scans or when respiratory or cardiac gating is used, the fraction of empty bins can even higher. Note that we expect that the sparsity of the TOF PET emission data of future PET systems will increase faster than linear compared to the improvement of the TOF resolution. This is because with better TOF resolution, every detected event carries more information such that fewer detected events are needed to reconstruct images with the similar variance [11]. Moreover, the method presented in this work is applicable not only to TOF PET image reconstruction, but to all imaging reconstrution problems with sparse data following a Poisson distribution such as (low count) SPECT or photon counting CT.

Contributions and Aim
To improve the efficiency of SPDHG in terms of required memory and computation time when reconstructing sparse TOF PET data, we propose and analyze a listmode extension of the SPDHG algorithm called LM-SPDHG that allows event-by-event processing using dedicated listmode forward and back projectors. We first derive LM-SPDHG from SPDHG and show that the convergence of LM-SPDHG is as fast as the convergence of SPDHG based on dedicated numeric examples in 2D. Moreover, we analyze the memory requirements for LM-SPDHG compared to SPDHG for typical scans acquired on state-of-the-art TOF PET scanners. We emphasize that the focus of this work is not on finding the optimal (non-smooth) prior or prior strength for a given clinical task in PET imaging. Instead, we aim to provide a framework enabling fast and memory efficient reconstruction of sparse TOF PET listmode data which will finally facilitate future research on the use of (non-smooth) priors in PET reconstruction.

Theory and Algorithms
Before deriving LM-SPDHG, the next subsection first of all shows how to reduce the memory requirements of SPDHG when reconstructing sparse TOF sinograms.

Memory-efficient SPDHG for sparse TOF sinograms
As shown in [12], the memory requirements for SPDHG can be substantially reduced by choosing a better initialization of the dual variable y. From Eq. (5) we can observe that for data bins i where d i = 0 (bins in the TOF emission sinogram with zero counts), since all other quantities are non-negative. Consequently, if we initialize all bins of y where the data sinogram d equals zero with one, these bins of y remain 1 during all iterations. This in turn means that these bins do not contribute to the update of δz, z,z, and x, since only the difference between y and y + is backprojected in line 8 of Algorithm 1 meaning that all bins without data do not have to be kept in memory during the iteration loop (lines 3 until 17). The only place where the empty data bins contribute is the initialization of z andz in line 2. However, this single back projection can be split into smaller chunks to also reduce the required memory of this step. While at a first glance, the initialization of y proposed above might seem artificial, in fact it directly corresponds to choosing the optimal value for those y i where d i = 0. To see this, note that an optimal solution (x,ŷ,ŵ) of (2) in particular needs to satisfy the optimality condition where ∂D * is the subdifferential of D * . Since this is equivalent toŷ ∈ ∂D(Px), and since D (the negative Poisson log likelihood) is differentiable this leads to the condition where the division is to be understood point-wise. With this we see that the proposed initialization directly corresponds to choosing the optimal value for y i where d i = 0, which is known explicitly in this case. As argued above, this improved initialization of y naturally reduces the memory requirements of SPDHG and also improves the speed of convergence when a "warm start" for x 0 is chosen. The latter can, e.g., be achieved by applying one iteration and a reasonable amount of subsets of the EM-TV algorithm [13,14]. Since in EM-TV every update is split into a classical EM step followed by a weighted denoising optimization problem in image space, EM-TV can be used in listmode as well by simply modifying the EM step.

Listmode SPDHG
As indicated by the name, emission data in listmode format are a chronological list N of detected events e ∈ N , where each event is characterized by a small set of (integer) numbers (e.g. the number of the two detectors, the discretized TOF difference, and a time stamp). To process the listmode data during reconstruction without binning it into a sinogram, we introduce the listmode forward operator P LM N mapping the image data x to a data-vector of dimension |N | via where i e is the sinogram bin in which event e was detected. Also, we denote by s LM with s LM e = s ie the listmode-based scatter and random estimate. Re-writing the gradient of the negative Poisson log likelihood using listmode data and the listmode operators is straightforward, such that any gradient-based PET reconstruction algorithm can be easily adapted to listmode data. For a Bayesian approach with prior R(K·) as in (1) this is less immediate, but we can show that indeed (1) can be equivalently re-written into a minimization problem involving only the listmode forward operator Algorithm 2 LM-SPDHG for PET reconstruction 1: Input event list N , contamination list s N 2: Calculate event counts µ e for each e in N (see text) 3: 6: Split lists N , s N and y N into n sublists N i , y Ni and s Ni 7: repeat 8: x ← proj ≥0 (x − T z) 9: Select i ∈ {1, . . . , n + 1} randomly according to (p i ) i 10: if i ≤ n then 11: y Ni ← y + Ni 14: end if 19: z ← z + δz 20: z ← z + (δz/p i ) 21: until stopping criterion fulfilled 22: return x that is of the same form as (1) as shown in appendix A. This allows us to extend the SPDHG algorithm for listmode data and yields the algorithm as shown in Algorithm 2. In contrast to the original SPDHG using binned data (sinograms), the forward and adjoint PET operators (lines 11 and 12) have been replaced by their listmode equivalents as in (8). Moreover, the dual variable for the data fidelity is replaced by the list y N , which has the same length as the measured event list N . If an event with a fixed sinogram bin i e occurs more than once in the event list N , it is also forward and back-projected multiple times in steps 11 and 12. To compensate for this fact, the algorithm divides by the event count 1 µ e before back projection in line 12. Note that (i) as shown in Fig. 1 in most standard acquisitions the event count of most events is 1 and that (ii) calculating the event count µ e which is a prerequisit creates a small pre-processing overhead (step 2). However, when implemented on a modern GPU this overhead is small compared to the computation time needed to calculate all iterations. Moreover, it is in the similar order of the time needed to unlist the native listmode data into a sinogram which is a prerequisit and pre-processing overhead for SPDHG. Another difference of LM-SPDHG compared to SPDHG is the fact that we split the data into n subsets by assigning every n-th event of the complete event list N to the n-th sub list -as commonly done when defining subsets in listmode OSEM. In that way, we can think of the subset listmode forward operator P LM N i as the full forward operator P with a sensitivity reduced by a factor of n. Accordingly, we set the step sizes associated with the subset listmode PET operators to At a first glance, the initialization of z andz in step 5 of Algorithm 2 might look odd. However, the first part of the expression is equivalent to applying the adjoint sinogram operator P T to a sinogram initialized according to (7). Note that step 5 still requires to calculate a sinogram back projection of a unity sinogram (P T 1) which we would like to avoid in listmode data processing. Unfortunately, avoiding this single sinogram backprojection is not possible. However, note that this "sensitivity image" is also needed in  Table 1: Estimation of required memory for SPDHG and LM-SPDHG assuming an image size of (300,300,125) and a TOF sinogram with 357 radial elemnts, 224 views, 1981 direct and oblique planes, and 27 TOF bins corresponding to a TOF PET scanner with 400 ps TOF resolution and 25 cm axial FOV for three counts levels. 5e8 counts approximately correspond to a high count 20 min late static FDG brain scan, 7e7 counts to an 80 s static FDG body bed position 1 h p.i., and 1e7 counts to a short early frame in a dynamic acquisition. In this estimation, we assume that every TOF PET listmode even can be encoded using 10 bytes.
the calculation of the step sizes T k and also in the listmode EM-TV algorithm that we use to initialize x. Hence, we recommend pre-computing the sensitivity image (P T 1) and storing it in memory.
In summary, compared to SPDHG, LM-SPDHG as shown in Algorithm 2 has the advantage that (i) only lists (N , y N , s N , µ N ) instead of complete sinograms (y and d) have to be stored in memory and (ii) all projections and back projections can be performed using listmode projectors. This means that for all acquisitions where the number of detected events is substantially smaller than the number of data bins, the required memory and computation time is reduced. The latter depends on the actual implementation of the sinogram and listmode projectors and on the computational hardware. According to our experience using a state-of-the art GPU implementation of Joseph projectors [15] for a TOF PET scanner with 400 ps TOF resolution and 25 cm axial field of view, a forward and back projection is faster in listmode if approximately less than 3e8 events have to be processed. For 7e7 and 1e7 counts, the projections are approximately faster by a factor of 3 and 5, respectively 2 . Note that the difference in the computation time of the projection for sinogram and listmode depends on the specific implementation -especially on the optimization of the memory access -and the computational hardware. A detailed comparison of the required memory for SPDHG and LM-SPDHG for different typical PET acquisitions is listed in Table 1.

Numerical Experiments
Recalling the fact that for most common TOF PET acquisitions, LM-SPDHG requires less memory and is also faster, one would probably prefer LM-SPDHG over SPDHG if the speed of convergence is the same. In the following, we describe a set of numerical experiments that we conducted to show that this is indeed the case. In the absence of an analytical solution to the optimization problems (1) and (2), we analyzed the convergence of LM-SPDHG and SPDHG with respect to a reference solution x * as done in [10]. Convergence was monitored by tracking the relative cost function where c(x) is the cost function to be optimized in (1) and x 0 is the initialization used for x. Moreover, convergence was also montiored in image space by tracking the peak signal-to-noise ratio with respect to where MSE is the mean squared error. The reference solution was obtained by running the deterministic PDHG (SPDHG without subsets) for 20000 iterations. Since running PDHG with 20000 iterations using realistic 3D TOF PET data takes a very long time (approx. 250 h), all numerical convergence experiments were performed using simulated 2D TOF PET data. A 2D software brain phantom with a gray to white w ← g/(βx) 8: x + ← argmin u≥0 j w j 2 (u j − z j ) 2 + R(Ku) 9: until stopping criterion fulfilled 10: return x matter contrast of 4:1 was created based on the brainweb phantom [16] and used to generate simulated 2D TOF data including the effects of limited spatial resolution, attenuation and a flat contamination mimicking random and scattered coincidences with a contamination fraction of 42%. The geometry and TOF resolution of the simulated 2D PET system was chosen to mimic one direct plane of a state-of-the art TOF PET scanner with a ring diameter of 650 mm and a 400 ps TOF resolution (sinogram dimension: 357 radial bins, 224 projection angles, 27 TOF bins). Noisy simulated prompt emission TOF sinograms and corresponding listmode data were generated for two different count levels (5e5 and 5e6 prompt counts). Unless stated otherwise, we always applied the EM-TV algorithm, summarized in Algorithm 3, using 1 iteration and 28 subsets to initialize x 0 and y 0 according to (7). To solve the weighted denoising problem in step 8 of the EM-TV algorithm, we applied the accelerated PDHG algorithm [8] using 20 iterations. In all SPDHG and LM-SPDHG reconstructios, the step size ratio γ was set to 3/ x 0 ∞ and ρ was set to 0.999. When splitting the data into n subsets, the vector of probabilities determining whether an update with respect to a subset of the data or with respect to the prior is done was set to such that on average an update with respect to the prior was done on every second update as suggested in [10]. In this work, we use the term iteration for 2n updates such that on average in every iteration the complete data are forward and backprojected once and n updates with respect to the prior are performed.
As benchmark examples for non-smooth priors, we consider two "Total Variation like" priors in this work. These priors have the form where Kx 2,1 is the sum over all entries of the pointwise Euclidean norm of Kx (mixed L2-L1 norm). The proximal operator for the convex dual of this prior is given by For the linear operator K, we first used the finite forward difference operator resulting in the classical Total Variation (TV) prior [1]. Moreover, we also implemented the Directional Total Variation (DTV) prior incorporating structural information by only considering the component of the finite difference vector in every voxel that is perpendicular to the spatial gradient of a prior image in that voxel [7]. In our convergence experiments using DTV, we assumed perfect structural prior information meaning that the prior image used for DTV, was the image itself with a flipped contrast. Note that in real acquisitions the quality of the available structural prior information is of course inferior. However, we do not expect that this fact strongly affects the convergence properties of SPDHG and LM-SPDHG.
To demonstrate that LM-SPDHG also works for realistic 3D TOF PET listmode data, we simulated and reconstructed 3D TOF PET data based on the 3D XCAT [17] phantom and a state-of-the-art TOF PET Note that in the calculation of the relative cost the same x 0 was used and that in both cases the same step size ratio γ was used. Moreover, in the very early iterations the relative cost of SPDHG is greater 1 meaning that the cost is worse compared the the initialization. This is due to the fact that (S)PDHG in general is a non-monotonic algorithm.
scanner with 20 cm axial field of view, a ring diameter of 650 mm and a 400 ps TOF resolution (sinogram dimension: 357 radial bins, 224 projection angles, 1296 direct and oblique planes, 27 TOF bins). As before, he data simulation included the effects of attenuation, finite resolution and a flat contamination with a contamination fraction of 42%. In total, 7e7 prompt counts were simulated corresponding approximately to a standard 80 s 1 h p.i. FDG liver bed position acquisition. Simulated listmode data were reconstructed using LM-SPDHG with 50 iterations and 224 subsets. Last but not least, we also reconstructed a real TOF PET data set of the NEMA image quality phantom acquired on a GE Discovery MI 4 ring TOF PET/CT [18] (20 cm axial FOV, sinogram dimension: 357 radial bins, 272 projection angles, 1261 direct and oblique planes, 29 TOF bins, 400 ps TOF resolution). The activity contrast between the sphere inserts and the background was 10:1 and the acquisition contained 4.7e7 prompt (2.2e7 true) coincidences. Reconstructions using a TV prior with β = 6 were performed with EM-TV using 100 iterations with 1 and 28 subsets, and with LM-SPDHG using 100 iterations and 224 subsets. We also performed a reconstruction of the same data set with LM-PDHG using 20000 iterations and 1 subset. This reconstruction was used as a reference to calculate the relative cost and PSNR. All reconstruction algorithms used in this work were implemented in python3 using the open-source "parallelproj" CUDA projector libraries available at https://github.com/gschramm/parallelproj. An Nvidia V100 GPU with 16 GB RAM was used to perform all reconstructions. Upon publication of this article, an open-source implementation of the LM-SPDHG algorithm will be made availalbe in the associated "pyparallelproj" python package. Preprocessing of the NEMA image quality phantom data set was performed using GE's duetto PET reconstruction toolbox v2.17.

SPDHG using a warm start vs a cold start
Before investigating the convergence behavior of LM-SPDHG, we first tested whether a warm start could already help to lead to faster convergence of SPDHG using sinograms. Figure 2 shows the results of SPDHG reconstructions with a cold (x 0 = 0 and y 0 = 0) and warm start as described above for a data set with 3e5 true counts and a TV prior with β = 0.03. It can be seen that SPDHG with the warm start performs better in terms of relative cost and PSNR with respect to the reference reconstruction. A similar trend was observed for higher count levels, with different β values and using the DTV prior indicating that as expected the warm start helps to accelerate convergence in the early iterations. Figure 3 summarizes the convergence comparison between sinogram SPDHG and LM-SPDHG for the same data set and prior as described in the previous subsection using the same warm start for both algorithms. As demonstrated by the convergence metrics and the reconstructed images, the convergence of sinogram SPDHG and LM-SPDHG is almost identical. This also holds for different count levels and for the structural DTV priors as shown in subfigures (b) and (c) and for different levels of regularization as shown in supplementary Figs. 1 and 2. For the example with the DTV prior shown in subfigure (c), the convergence in terms of PSNR seems to be even slightly faster with LM-SPDHG, which can be also seen in the difference images with respect to the reference reconstruction.

Convergence of LM-SPDHG compared to listmode EM-TV
A comparison between the convergence of LM-SPDHG using 224 subsets and listmode EM-TV using 1, 7 and 28 subsets is shown in Fig. 4. In all sub figures, LM-SPDHG converges much faster than EM-TV. Interestingly, when using listmode EM-TV with more than one subset, the convergence metrics saturate meaning that the algorithm remains on a limited cycle and the optimal solution is not reached. When using one subset, the convergence of EM-TV is much slower but does not seem to saturate after 100 iterations. Figure 5 shows the convergence of LM-SPDHG as a function of the number of data subsets used for low (3e5) and high (3e6) counts and the TV and DTV prior. In general, using 224 subsets leads to faster convergence compared to 56 and 112 subsets. However, for the TV prior, especially at low counts, there is almost no difference between using 56, 112, and 224 subsets. In contrast, for the DTV prior, the difference between using 56, 112, and 224 subsets is more pronounced. For a fixed number of subsets, convergence is faster for high count compared to low count data sets. Figure 6 shows the results of the reconstruction of the listmode data simulated from the 3D XCAT phantom using a TV prior with β = 0.03. From the plot of the cost function and the reconstructions shown in the middle, we can see that reasonable convergence is reached after approximately 25 iterations with 224 subsets. Note that in contrast to our 2D experiments, the cost function initially increases until the 3rd iteration before it decreases and stabilizes after around 25 iterations. As mentioned above, calculating the projections in listmode for 7e7 counts is roughly a factor of 3 faster compared to sinogram-based processing which substantially speeds up every iteration. Note that in our proof-of-concept implementation, the effective speed up is less than a factor of 3 since we did not implement the gradient operators and proximal mappings on a GPU yet which leads to overhead due to gradient-based updates. However, we expect that, once properly implemented on GPUs, this overhead will be small compared to the time needed to calculate the actual TOF projections. Figure 7 shows the results of the reconstructions of the listmode data from the acquisition of the NEMA image quality phantom on the GE Discovey MI 4 ring TOF PET system. After 100 iterations, LM-SPDHG using 224 subsets has lower cost compared to both EM-TV reconstructions and is closest to the referene reconstruction in terms of PSNR. From the differene images it can be seen that the stochastic LM-SDPHG mostly differs from the reference in the background outside the phantom. In contrast, both EM-TV    reconstructions show bigger residual differences within the phantom which are most pronounced in the "warm" background for EM-TV using 28 subsets and in the "cold" central insert for EM-TV using 1 subset.

Discussion
The results of our numerical experiments presented in this work demonstrate that the speed of convergence of LM-SPDHG is essentially the same as the one of the original SPDHG using sinograms. For clinical acquisitions with state-of-the-art TOF PET scanners resulting in extremely sparse data, LM-SPDHG has two distinct advantages compared to SPDHG. First, during the iterations all forward and back projections can be performed in listmode which is faster compared to sinogram projectors for sparse TOF data. Second, as shown in Table 1, the memory requirements are substantially reduced such that even for an aquisition with 5e8 prompt coincidences only 12.5 GB of memory are needed. This actually enables a pure GPU-implementation of LM-SPDHG avoiding intermediate memory transfer between the host and a state-of-the-art GPU with a memory of approximately 16 GB. Note that in our proof of concept implementation of LM-SPDHG we used a hybrid computing approach where only the TOF PET forward and back projections where computed using a GPU. We expect that once properly implemented purely on a GPU, the computation time required for LM-SPDHG will further decrease since memory transfter between host and GPU is a bottleneck in our implementation. Note that a pure GPU implementation of the conventional SPDHG algorithm for modern TOF PET data is more complicated since more than 50 GB of GPU is required. In this proof of concept work, we only used two non-smooth TV-like priors to benchmark LM-SPDHG. Note that in general, SPDHG and LM-SPDHG can be also used for smooth priors as long as the proximal operator of the convex dual of the priors can be efficiently calculated which is e.g. the case for prior Figure 7: Convergence of LM-SPDHG for reconstruction of real 3D TOF data from an acqusition of the NEMA image quality phantom on a GE Discovey MI 4 ring TOF PET system with a TV prior and β = 6 with 2.2e7 true (4.7e7 prompt) counts. First column, top two rows: transaxial and coronal slice of reference LM-PDHG reconstruction using 20000 iterations and 1 subset. Second column, top two rows: LM-SPDHG reconstruction using 100 iterations and 224 subsets. Third column, top two rows: EM-TV reconstruction using 100 iterations and 1 subset. Fourth column, top two rows: EM-TV reconstruction using 100 iterations and 28 subsets. Fist column, bottom two rows: relative cost and PSNR to reference reconstruction. Second, third, fourth column, bottom two rows: absolute difference compared to reference reconstruction.
penalizing the squared L2 norm or Huber TV. However, for smooth priors the speed of convergence of (LM)-SPDHG should be benchmarked against other stochastic gradient-based optimization techniques such as the stochastic variance reduced gradient (SVRG) [19] or SAGA [20] which is beyond the scope of this work. We also note that in this work we used a scalar spatially-invariant prior strength β which leads to a spatially-variant local pertubation response (LPR). For pratical applications where a spatiallyinvariant reponse is desirable, a spatially-varient prior strength should be used as described in [21,22]. We emphasize that the listmode EM-TV algorithm is still a very practical and useful algorithm to approximate the solution of the optimization problem (1) for listmode data and non-smooth priors. However, as shown in Fig. 4, it seems that when using ordered subsets for acceleration, EM-TV does not reach the optimal solution but rather remains on a limit cycle similar to the behaviour of OSEM. Whether the difference between the limit cycle and the true optimal solution is of importance for a given count level and clinical task, should be investigated in the future.
As shown in Fig. 5, choosing the "optimal" number of subsets for LM-SPDHG seems to depend on the prior and the number of acquired counts. For a reconstruction of a high-count scan with DTV, it is possible that LM-SPDHG with more than 224 subsets converges even faster. A detailed investigation of the choice of the number of subsets for different acquistion and reconstruction scenarios is left for future research. Last but not least, we would like to note that according to our experience, SPDHG and LM-SPDHG for TOF PET reconstruction also converge for ρ slightly greater than one which might be of interest for practical applications where stopping as early as possible can be important. Suppl. Fig. 3 shows that in our 3D XCAT example, convergence can be further accelerated by using ρ = 8 and γ = 30/ x 0 ∞ . A detailed investigation of the convergence speed as a function of ρ and γ potentially involving an iteration-dependent choice of ρ [23] is, however, beyond the scope of this manuscript and left for future research.

Conclusion
For sparse TOF data, the proposed LM-SPDHG algorithm severely reduces the memory requirements and computation times compared to the original SPDHG enabling the application of LM-SPDHG in routine clinical practice where short reconstruction times are crucial. respectively. Now the last line of the above reformulation provides an equivalent data term which (after removing the constant e∈N log(d ie )) can be used to obtain an equivalent list-mode-based reformulation of (1).
Applying the SPDHG for this reformulation results in the proposed LM-SPDHG Algorithm 2. Note that there, the non-rescaled operator P LM and estimate s LM appear since i) on line 11 of the algorithm, the rescaling of P LM and s LM cancels out with the rescaling of P LM in the step size S k of (9), and ii) the rescaling of P LM in line 12 is incorporated explicitly. Furthermore, note that we use the full forward operator P instead of P LM N k for the stepsize T k , for reasons explained in the paragraph above the respective equation (9).