in medical image analysis to study deformations between images. Building on image registration approaches, image regression models [@ref:georegression; @hong2012simple; @hong2012metamorphic; @ref:nikhil; @fletcher2013geodesic; @hong2014time; @singh2014splines; @hong2014geodesic; @singh2015splines; @hong2016parametric] have been developed to analyze deformation trends in longitudinal imaging studies. One such approach is geodesic regression (GR) [@ref:georegression; @ref:nikhil; @fletcher2013geodesic] which (for images) build on the large displacement diffeomorphic metric mapping model (LDDMM) [@ref:lddmm]. In general, GR generalizes linear regression to Riemannian manifolds. When applied to longitudinal image data, it can compactly express spatial image transformations over time. However, the solution to the underlying optimization problem is computationally expensive. Hence, a simplified, approximate, GR approach has been proposed [@ref:hong] (SGR) to decouple the computation of the regression geodesic into pairwise image registrations. However, even such a simplified GR approach would require months of computation time on a single graphics processing unit (GPU) to process thousands of 3D image registrations for large-scale imaging studies such as ADNI [@jack2015magnetic]. The primary reason computational bottleneck for SGR are the optimization required to compute pair-wise registrations.

![image](image/Pipeline-Kwitt){width="\textwidth"}

1ex Recently, efficient approaches have been proposed for deformable image registration [@dinggang; @miao2016real; @nonrigid; @ref:yang2016; @quicksilver; @zhang2017frequency]. In particular, for LDDMM, which is the basis of GR approaches for 