magnetic resonance imaging ( mri ) is an important application of compressive sensing ( cs ) [ 14 ] . cs in mri   has the potential to significantly improve both the speed of acquisition and quality of mr images , but requires an iterative reconstruction that is more computationally intensive than traditional inverse fourier reconstruction . compressive sensing accelerates mr acquisitions by reducing the amount of data that must be acquired . reconstruction of this partial data set is then accomplished by iteratively constraining the resulting image to be sparse in some domain while enforcing consistency of the measured subset of fourier data . one practical barrier to the routine adoption of cs mri is the delay between acquisition and reconstruction of images . compressed sensing solvers work almost entirely with vector and image arithmetic , making them an excellent candidate for acceleration through using graphics processing units ( gpus ) for parallelization . here we illustrate how gpus can be used to achieve significant increases in speed of cs reconstructions of large mri data sets . high - end video cards can contain hundreds of separate floating - point units , allowing for massively parallel computations at a fraction of the cost of cpu - based supercomputers , as measured on a per gigaflop basis ( one gigaflop is equivalent to one billion floating point operations per second ) . gpu computing works best in single instruction , multiple data ( simd ) situations , such as the solution of large systems of linear equations . the power of gpu computing is already being realized in several advanced medical image reconstruction applications [ 69 ] . the reconstruction platform tested was a high - performance gpu server designed to serve an mri scanner as a dedicated reconstructor . this system contained a six - core , 2.67  ghz xeon x5650 with 12  gb of 1333  mhz ddr3 ram , 32  kb of l1 cache , 256  kb of l2 cache , and 12 mb of l3 cache . the tesla c2050 contains 448 cuda cores arrayed as 14 streaming multiprocessors ( sms ) with 32 cores per sm . it has 3  gb of vram and 64  kb of memory per sm . all sms share 768  kb of l2 cache . the optimization problem we are choosing to explore is the reconstruction of partial mri data sets using cs . mri is proving to be a fertile application for cs because many mr images can be highly compressed by transform coding with little loss of information . conventional mri data is acquired according to the fourier sampling pattern required to satisfy the nyquist criterion while producing an image of a given field of view and spatial resolution . to sample the fourier domain compressively ,  only a random subset of the full nyquist fourier sampling scheme is acquired . reconstruction of the conventional fully sampled mri data requires simply an inverse fourier transform , 
 ( 1)u = f1b , 
 but the inverse problem becomes underdetermined when part of the fourier data is omitted , so approximate methods must be used . one highly successful method in particular has been to formulate the cs mri reconstruction as a sparse recovery problem , in which an image is found that is consistent with the acquired fourier data while having the sparsest representation in a chosen basis ( e.g. , gradient and wavelet ) . the typical formulation of the reconstruction of a complex image u from a partial fourier data set b is then 
 ( 2)u = arg  minx||x||1+||axb||22 , 
 where a is a measurement operator that transforms the sparse representation x to the image domain then performs the subsampled fourier measurement , and the l1 and l2 norms are , respectively , 
 ( 3)||x||1=i|xi|,||x||22=ixixi , 
 where the bar denotes complex conjugation . with the addition of the l1 norm , the problem is more difficult to solve , and iterative techniques such as interior point methods [ 10 , 11 ] , iterative soft thresholding [ 12 , 13 ] , and gradient projection [ 1416 ] are typically employed . the open - source split bregman code of goldstein and osher   was chosen as the starting point for the gpu - based cs solver . this solver was chosen for its rapid convergence and lack of array reduction steps , which hinders parallelization . we modified the original code to work with jacket 1.8.0 ( accelereyes , atlanta , ga ) and matlab r2010b ( mathworks , natick , ma ) . algorithm 1 briefly outlines the procedure for running the split bregman reconstruction on the gpu . note that the split bregman algorithm runs for a fixed number of iterations , so there is no variation in run time due to different descent trajectories as with a tolerance - based stopping criterion . furthermore , the choice of image reconstructed has no bearing on the results , since the fixed number of iterations ensure that the same number of operations are performed on any input data set .   at the beginning of the reconstruction , the fourier data in main memory must be transferred to the gpu with jacket 's gsingle and gdouble matlab commands . next , temporary storage is allocated on the gpu using jacket 's gzero and gones commands . all subsequent arithmetic operations , including the fourier transform , are carried out on the gpu using function overloading . function overloading simplifies code syntax by enabling a single function to encapsulate different functionality for different types of arguments . the specific behavior is typically chosen by the compiler or at run time . in our case , matlab automatically calls the jacket library if an operation is requested on a matrix that lies in gpu memory , while identical operations on a matrix in main memory are carried out with matlab 's built - in functions . after the last loop iteration on the gpu , the solution is transferred back to main memory with overloaded versions of matlab 's double and single commands ; all temporary storage is automatically freed . the first was a pure matrix multiplication , designed to measure practical peak floating - point performance of the cpu and gpu as realized by jacket 1.8.0 . to remove dependencies on the multithreading performance of matlab r2010b and provide easier comparison of the cs reconstructions the second experiment , and the focus of this work , was a cs mri reconstruction of a t1-weighted breast image subjected to a 50% undersampling in fourier ( spatial frequency ) space . total variation was used as the sparsity constraint and was defined as the sum of the magnitudes of pixels in the gradient image . cs mri reconstructions were performed for powers - of - two image sizes ranging from 32 to 8192 ( up to 4096 only for double precision , due to memory limitations ) . this range covers the range of realistic mr acquisition matrix sizes for 2d scans with allowance for specialty techniques at very low or very high resolutions or future developments in imaging capabilities . the largest matrix sizes can also be indicative of the performance of three - dimensional reconstruction problems . for example , a 256 data set is the same size as a 4096 one . all reconstructions were timed eleven times , with the first iteration discarded and the following ten iterations averaged . jacket uses just - in - time ( jit ) compilation to improve performance of repeated function calls , so the first call to the jacket library is slowed by this compilation step . the cuda driver , upon initial invocation , optimizes the low - level gpu code for the particular hardware being used . these two processes increase code performance across multiple runs but reduce it for the initial function calls .   for timing experiments , warmed up  with a similar problem before timing the full - scale computation . here we used the simplest approach of discarding the first iteration . in principle , though , the warmup problem can be much smaller in data size as long as it uses the same set of jacket functions needed in the full - size reconstruction . table 1 shows the result of the matrix multiplication experiment with and without cpu multithreading enabled . in terms of cpu performance , matlab r2010b accelerated the cpu - based matrix multiplication by ~6 on this six - core processor using multithreading . using the gpu then yielded an additional factor of ~5 beyond this , with a combined speedup of ~30 over a single cpu core . based on the measured gpu single- and double - precision performance of 650 and 311 gflops , respectively , we can see that matlab r2010b combined with jacket 1.8.0 reached 60% , and 63% of the theoretical maximum single- and double - precision performance , respectively , of the gpu . the results of the cs mr image reconstructions are shown in tables 2 and 3 and figure 2 . the maximum speedup was 27 for a single - precision image of size 2048 . for a typical double - precision mri acquisition matrix of 256 to 512 figure 2 shows the speed advantage of the gpu - based code over using both one and multiple cpu cores . the difference between enabling multithreading or not was large , and should be considered to fairly evaluate the speed improvement of using a gpu . as can be seen in figure 2 , the speedup was less than one for images smaller than about 64 , followed by a rapid gain in speedup factor for images of 2048 and a decline in speedup factor for the largest matrix sizes . despite the performance falloff , the single - precision gpu code was able to reconstruct images up to 8192 before running out of gpu memory since single - precision matrices require half the storage per element of double - precision matrices . two notable features are evident in tables 2 and 3 . first , the gpu reconstruction times are very similar for images below 256 , regardless of numerical precision . this is due to the communication overhead in transferring the data from the cpu memory to the gpu memory . for the smallest images , this cost dominated the computation time and , for images below 64 , even caused the gpu reconstruction to take longer than the cpu reconstruction this suggests that for very small acquisition matrices an efficient gpu reconstruction should combine multiple 2d data sets , such as different slices or echoes , into one reconstruction problem . the second interesting feature of tables 2 and 3 is that gpu reconstruction times for images of size 1024 and smaller were effectively instantaneous . a typical rapid gradient echo sequence may employ a repetition time of 5  ms , so a 50% undersampling of an n fourier matrix would require roughly 2.5  n  ms to acquire ( ignoring other acceleration methods for simplicity ) . the smallest image tested here would thus take 40  ms to acquire , while the cs reconstruction would take only 60  ms to complete , even with the severe communication penalty . the largest double - precision image we tested ( 4096 ) would take 15  s to acquire and 10  s to reconstruct . thus our gpu - based platform has the capacity to produce iteratively reconstructed cs - accelerated gradient echo images in real time for some mri applications . the decline in acceleration at 4096 was likely due to memory limitations ( yalamanchili , private communication ) . the nvidia tesla c2050 card is designed such that each streaming multiprocessor ( sm ; comprised of 32 floating point units ) shares a single 64  kb block of local memory . for a double - precision complex matrix multiplication , each matrix element requires 16 bytes of storage , so a maximum of 4096 elements can be stored in the sm 's shared memory . matlab stores arrays in column major order , so a matrix column must fit entirely into the shared memory of the sm in order to minimize memory access overhead . thus memory access patterns will become inefficient above matrix sizes of 2048 for double - precision complex matrix multiplication . we do in fact see a leveling off of performance at 1024 and a dramatic decline in performance above 2048 , which is consistent with this prediction .   many alternate methods of accelerating the cs reconstruction on a gpu platform exist , including writing custom low - level code in cuda c or opencl , using free low - level libraries like cula and cufft , using matlab 's built - in gpu library through the parallel computing toolbox ( r2010b later only ) , or using other free high - level matlab - cuda interfaces , such as gpumat . the pace of algorithmic development is accelerating in compressed sensing , as demonstrated in figure 3 , and one advantage of using a high - level interface to the gpu is rapid prototyping and implementation . debugging time is shorter with high - level languages , and coding effort is reduced , allowing the newest algorithms to be implemented quickly while still retaining a majority of the theoretical computational benefit of the gpu . matlab 's built - in gpu library does not support the array indexing needed for the gradient operation , so we could not compare it here . gpumat is free and open source and a viable alternative to jacket in theory , but we were unable to use it with our reconstruction code due to an unresolvable memory access error . jacket is the only high - level software package to support sparse matrix operations on the gpu , which allows classes of compressed sensing algorithms that use sparse matrix operators to be used . also , the overloaded gpu functions in jacket are implemented as mex files , which are precompiled c / c++ functions callable from within matlab ; so writing custom cuda subroutines could only eliminate the function call overhead and not speed up the individual simd operations . we have shown that gpu computation significantly accelerated cs mri reconstruction of all but the smallest of the tested image sizes . the combination of matlab and jacket yields a processing package that is able to realize over half of the theoretical maximum performance of the gpu , while requiring minimal code development . the speedup realized by the gpu for the smallest images was progressively hampered by communication overhead , while the largest images suffered from the limited gpu ram . the optimal image dimensions , however , seem to be serendipitously close to that of high - resolution mri data ; so gpu computing , coupled with the goldstein and osher split bregman algorithm , appears to be a well - suited platform for rapid cs mri reconstruction . future improvements to these methods include algorithm modifications to allow unified reconstruction of multiple two - dimensional or a single three - dimensional fourier data set on the gpu with a single call to the reconstructor , thus reducing the communication penalty . additional cores on the gpu card could allow higher acceleration , since we found that the split bregman algorithm parallelizes extremely well . and finally , more ram onboard the gpu would allow larger data sets to be reconstructed more efficiently .