magnetic resonance imaging ( mri ) is a noninvasive medical imaging modality commonly used to investigate soft tissues in the human body .
clinically , mri is attractive as it offers flexibility , superior contrast resolution , and use of only nonionizing radiation .
however , as the duration of a scan is directly proportional to the number of investigated spectral indices , obtaining high - resolution images under standard acquisition and reconstruction protocols can require a significant amount of time .
for example , during long examinations , patients often exhibit involuntary ( e.g. , respiration ) and/or voluntary motion ( e.g. , active response to discomfort ) , both of which can impart spatial blurring that may compromise diagnosis .
also , high temporal resolution is often needed to accurately depict physiological processes . under standard imaging protocols , spatial resolution
must unfortunately be sacrificed to permit quicker scan termination or more frequent temporal updates . rather than executing a low spatial resolution exam ,
contemporary mri protocols often acquire only a subset of the samples associated with a high - resolution exam and attempt to recover the image using alternative reconstruction methods such as homodyne detection   or compressive sensing ( cs ) .
cs theory asserts that the number of samples needed to form an accurate approximation of an image is largely determined by the image 's underlying complexity [ 2 , 3 ] .
thus , if there exists a means of transforming the image into a more efficient ( i.e. , sparse or compressible ) representation , less time may actually be required to collect the data set needed to form the high - resolution image .
background - suppressed contrast - enhanced mr angiography ( ce - mra ) is a very natural clinical target for cs methods .
as the diagnosis of many conditions like peripheral vascular disease are based on both vessel morphology and hemodynamics , high spatial and temporal resolution images are consequently needed .
although several authors ( e.g. , [ 58 ] ) have successfully demonstrated the application of cs methods to ce - mra , the computationally intensive nature of these applications has so far precluded their clinical viability ( e.g. , published cs reconstruction times for a single 3d volume ( ce - mra or not ) are often on the order of hours [ 4 , 7 , 911 ] ) . as the results of a ce - mra exam are often needed as soon as the acquisition completes ( either for immediate clinical intervention or to guide additional scans ) , it is not practical to wait for the result of any currently implemented cs reconstructions . instead , linear or other
noniterative reconstructions that can be executed online ( e.g. , ) must be used even if they provide suboptimal results . with the goal of reducing cs+mri reconstruction times to clinically practical levels , several authors have recently considered the use of advanced hardware environments for their reconstruction implementations .
most of these techniques have focussed on the reconstruction of mri data acquired using phased - array ( i.e. , multicoil ) receivers , as this is the dominant acquisition strategy used in clinical practice . chang and ji
[ 13 , 14 ] considered a coil - by - coil approach to reconstructing phased - array mri data .
although this strategy leads to natural task parallelization , with each element of a multicore processor independently handling the reconstruction on one coil image , they only demonstrated reconstruction times on the order of minutes , per 2d slice which is not clinically viable .
moreover , disjoint reconstruction of phased - array mri data is well known to exhibit suboptimal performance when compared against joint reconstruction strategies like sense   and grappa   and so , this strategy is of limited utility .
murphy et al .   later demonstrated that the spirit reconstruction algorithm , which is a generalization of grappa , can be significantly accelerated using graphics processors .
they generated high - quality parallel image reconstructions in on the order of minutes per 3d volume , representing a significant advance towards clinical feasibility .
[ 7 , 8 , 19 ] demonstrated cs reconstructions of time - resolved 3d ce - mra images acquired using parallel imaging and a state - of - the - art cartesian acquisition with a projection - reconstruction - like sampling ( capr ) strategy   also in a matter of only minutes per 3d volume using an advanced code implementation on a cluster system .
their algorithm , which is essentially a generalization of sense   that employs auxiliary sparsity penalties and an efficient inexact quasi - newton solver , was demonstrated to yield high quality reconstruction of 3d ce - mra data acquired at acceleration rates upwards of 50x . in this paper
, we focus on trzasko et al . 's cs+mri reconstruction strategy for 3d ce - mra and investigate the development , optimization , and performance analysis on several modern parallel architectures , including the latest quad- and six - core cpus , nvidia gpus , and intel many integrated core architecture ( intel mic ) .
our optimized implementation on a dual - socket six - core cpu is able to reconstruct a 256  160  80 volume of the neurovasculature from an 8-channel , 10x accelerated ( i.e. , 90% undersampled ) data set within 35 seconds , which is more than a 3x improvement over other conventional implementations [ 8 , 19 ] .
furthermore , we show that our cs implementation scales very well to the larger number of cores in today 's throughput - oriented architectures .
our nvidia tesla c2050 implementation reconstructs the same dataset in 19 seconds , while our intel 's knights ferry further reduces the reconstruction time to 12 seconds , which is considered clinically viable .
finally , our research simulator shows that the reconstruction can be done in 6 seconds on 128 cores , suggesting that many - core architectures are a promising platform for cs reconstruction .
, the targeted data acquisition process for one time frame can be modeled as 
 ( 1)[g1g2gc]=[12c]f+n , 
							 where f is a discrete approximation of the underlying image of interest , c is the cth coil sensitivity function ,  is the 3d discrete fourier transform ( dft ) operator ,  is a ( binary ) sampling operator that selects a prescribed subset of k - space values , n is complex additive white gaussian noise ( awgn ) , and gc is signal observed by the cth coil sensor . as described in ,
raw capr k - space data is background - subtracted and view - shared prior to execution of the cs reconstruction procedure .
it was demonstrated in [ 8 , 19 ] that background - suppressed ce - mra images acquired by systems of the form in ( 1 ) can be accurately recovered by ( approximately ) solving the following unconstrained optimization problem : 
 ( 2)v=argmin  v  j(v ) , 
							 where the cost functional 
 ( 3)j(v)=np(dnv)+c=1c||cvhc||22 , 
 d 
 n is the finite spatial difference operator for some offset direction n ( in the neighborhood  ) , and the penalty functional 
 ( 4)p(v)=x(v(x ) ) , 
							 for some concave metric functional ( ) . following , the nonconvex laplace functional 
 ( 5)()=1exp  ( 1|| ) ,  [0, ) 
							 is herein adopted .
although not considered here , the 1-norm ( ( ) = || ) could also be used if so desired . in
, an efficient inexact quasi - newton algorithm was proposed for ( approximately ) solving ( 2 ) to reconstruct capr ce - mra images . for completeness , this algorithm is briefly reviewed .
recall that complex quasi - newton iterations   are typically of the form 
 ( 6)vi+1=vib1(vi)l(vi ) , 
							 where the gradient of j(v ) ( taken with respect to v ) and b( ) is an approximation of the complex hessian of j(v ) .
the term  inexact  arises when i is only approximately determined , such as via truncated conjugate gradient ( cg ) iteration .
more specifically , given ( 3 ) , 
 ( 7)l(vi)=ndn(dnvi)dnvi  + c=1cc(cvikc(t ) ) , 
							 where the (  > 0 smoothed ) diagonal operator 
 ( 8)(dnvi)(x , x)=12|[dnvi](x)|i(|[dnvi](x)|i)|[dnvi](x)|i . 
							 in their work , trzasko et al .
[ 8 , 19 ] adopted the following analytical linear hessian approximation : 
 ( 9)b(vi)=2ndn(dnvi)dn+c=1ccc , 
							 which can be considered a generalization of vogel and oman 's  lagged diffusivity  model   for total variation ( tv ) denoising and deblurring . for improved convergence ,
decreasing continuation is also performed on the functional smoothing parameter ,  . in , an efficient c++ implementation of the above algorithm employing the templated class framework described by borisch et al .   and both the mpi ( http://www-unix.mcs.anl.gov/mpi/ ) and openmp ( http://www.openmp.org/ ) libraries was described and executed on an 8-node dedicated reconstruction cluster , where each node had two 3.4  ghz intel xeon processors and 16  gb memory . for a single 256  160  80 head volume reconstruction from 8-channel data and only 6 difference neighbors , reconstruction times of slightly less than 2 minutes
although these times represent a significant advancement over other existing works , they are still too long for routine clinical use .
we used five datasets ( two artificial and three clinical ) to analyze the cs performance on three platforms : intel cpus , nvidia gpus , and intel mic .
five datasets are used for our experiments , whose volume size and memory footprint are : ( 256  64  64 , 224  mb ) , ( 256  160  32 , 280  mb ) , ( 256  160  80 , 700  mb ) , ( 256  160  84 , 735  mb ) , and ( 256  160  88 , 770  mb ) , in the order of dataset 1 to 5 , respectively .
datasets 1 and 2 were artificially generated , dataset 3 represents a noncontrast - enhanced brain , and datasets 4 and 5 represent contrast - enhanced vasculature .
all mri data were acquired on a 3 t ge signa scanner ( v.20 ) using an 8-channel head array using the capr acquisition sequence . prior to reconstruction ,
view sharing was performed on datasets 35 , and background reference subtraction on dataset 4 and 5 as described in .
for all experiments , 5 outer and 15 inner ( cg ) iterations were executed under w = 1 ( corresponding to 26 finite difference neighbors ) .
figure 1 shows the current clinical   and cs - type reconstruction   results for dataset 5 .
the cs reconstruction results for all versions optimized for different architectures discussed below were visually identical to that shown here .
all architectures are compliant with ieee single - precision floating - piont arithmetic standard . 

intel core i7 processorthe intel core i7 processor is an 86-based multicore architecture which provides four / six cores ( 731  m/1.17  b transistors ) on the same die .
it features a superscalar out - of - order core supporting 2-way hyperthreading and 4-wide simd .
each core is backed by 32  kb l1 and 256  kb l2 caches , and all cores share an 8  mb/12  mb l3 cache .
quad and six - core cpus provide 100 gflops and 135 gflops of peak single - precision computation respectively , as well 32  gb / s of peak memory bandwidth . to optimize cs on core i7 processors , we took advantage of its sse4 instructions using the intel icc auto - vectorizing compiler as well as hand - vectorization intrinsics .
we parallelized the code with openmp and adopted a highly optimized fft implementation from intel 's math kernel library ( mkl ) 10.2 . 
 
the intel core i7 processor is an 86-based multicore architecture which provides four / six cores ( 731  m/1.17  b transistors ) on the same die .
it features a superscalar out - of - order core supporting 2-way hyperthreading and 4-wide simd .
each core is backed by 32  kb l1 and 256  kb l2 caches , and all cores share an 8  mb/12  mb l3 cache .
quad and six - core cpus provide 100 gflops and 135 gflops of peak single - precision computation respectively , as well 32  gb / s of peak memory bandwidth . to optimize cs on core i7 processors , we took advantage of its sse4 instructions using the intel icc auto - vectorizing compiler as well as hand - vectorization intrinsics .
we parallelized the code with openmp and adopted a highly optimized fft implementation from intel 's math kernel library ( mkl ) 10.2 . 

nvidia teslathe nvidia tesla c2050 [ 27 , 28 ] ( 3  b transistors ) provides 14 multiprocessors , each with 32 scalar processing units that share 128  kb of registers and a 64  kb on - chip memory .
32 scalar units are broken into two groups , where each group runs in lockstep .
it features a hardware multithreading , which allows hundreds of thread contexts running concurrently to hide memory latency .
the on - chip memory is software configurable , and it can be split into a 48  kb cache and a 16  kb shared memory space , or vice versa .
its peak single - precision computing performance is about 1.03 tflops and its on - board gddr memory provides up to 144  gb / s bandwidth .
cuda allows programmers to write a scalar program that is automatically organized into thread blocks to be run on multiprocessors .
cuda provides an open source fft library ( cufft 3.1 ) although more optimized fft implementations such as nukada and matsouka 's   have been published . 
 
the nvidia tesla c2050 [ 27 , 28 ] ( 3  b transistors ) provides 14 multiprocessors , each with 32 scalar processing units that share 128  kb of registers and a 64  kb on - chip memory .
32 scalar units are broken into two groups , where each group runs in lockstep .
it features a hardware multithreading , which allows hundreds of thread contexts running concurrently to hide memory latency .
all multiprocessors share a 768  kb l2 cache . the on - chip memory is software configurable , and it can be split into a 48  kb cache and a 16  kb shared memory space , or vice versa .
its peak single - precision computing performance is about 1.03 tflops and its on - board gddr memory provides up to 144  gb / s bandwidth .
cuda allows programmers to write a scalar program that is automatically organized into thread blocks to be run on multiprocessors .
cuda provides an open source fft library ( cufft 3.1 ) although more optimized fft implementations such as nukada and matsouka 's   have been published .
intel mic is an aubrey isle-   based platform and intel 's knights ferry   is its first silicon implementation with 32 cores running at 1.2  ghz .
it is an 86-based many - core processor based on small in - order cores that combines the full programmability of today 's general - purpose cpu architectures with the compute throughput and memory bandwidth capabilities of modern gpu architectures .
each core is a general - purpose processor , which has a scalar unit based on the pentium processor design , as well as a vector unit that supports 16 32-bit float or integer operations per clock .
it is equipped with two levels of cache : a low latency 32  kb l1 cache and a larger globally coherent total 8  mb l2 cache that is partitioned among the cores .
because intel mic is based on 86 , it provides a natural extension to the conventional 86 programming models .
thus , we could use similar data and thread level implementation as on core i7 processors . 

figure 2(a ) shows the overview of our cs implementation . the targeted cs reconstruction algorithm is composed of multiple iterations of 3d matrix arithmetics .
we divide the reconstruction model outlined in ( 6 ) into six stages based on the loop structure ( denoted as stage1 , stage2 , etc . ) .
more specifically , stage1 , stage2 , and stage3 correspond with computation of the left and right terms of ( 7 ) , respectively .
analogously , stage4 , stage5 , and stage6 correspond with computation of the left and right terms of ( 9 ) , respectively .
each stage performs a series of matrix computations such as elementwise additions and 3d ffts . the pie chart in figure 2(b ) shows the execution time breakdown of the key kernels .
fft3d ( performed in stage2 and stage5 ) is the most time - consuming and accounts for 46% of the total execution time . to achieve optimal performance
simple elementwise matrix arithmetics ( matrix ) are the second most time - consuming kernels . because they stream large amount of data from / to the main memory , our optimizations focus on hiding latency and utilizing bandwidth efficiently .
diff3d ( in stage1 and stage4 ) calculates the differences from the original matrix to its shifted copy . since the same data are used multiple times , we block the matrix to exploit data reuse in fast on - die memories .
geval ( in stage1 and stage3 ) performs transcendental operations such as division and exponentiation that are implicit within ( 8) . on gpus ,
we take advantage of fast math functions ; however , the performance gain due to the faster math is marginal because geval comprises only 7% of the total execution time .
naive implementation often misses a large amount of performance potential . in order to realize maximum performance potential
, we discuss a number of important architecture - aware optimizations for our cs implementation .
our optimization techniques are general , and thus can be applied to all three architectures ( cpu , gpu , and intel mic ) . in particular , since a cpu and mic share the same programming model , an optimized cpu code can be ported to intel mic without much modification . however , the cuda programming model is quite different from a cpu 's .
it requires significant effort to port a cpu code to gpus , and it may be easier to program it from scratch for gpus .  
because each element of the matrix is computed independently in most kernels , we exploit the element - level parallelism . in other words
, we can pack multiple elements into a vector computation and/or divide elements among multiple threads arbitrarily , without concern of data dependency .
a kernel is usually is composed of three nested loops for each of the dimensions , x , y , and z , as shown in figure 3(a ) . because there is no data dependency between elements , it is possible to vectorize within any of the iterations .
however , it is most efficient to vectorize the inner - most loop , since it exhibits sequential memory accesses along the x axis .
vectorization along the y or z axis requires gathering elements from nonsequential memory locations , resulting in poor performance . in addition , most loops in cs do not contain data - dependent control flow diversions ( i.e. ,  if  statements ) , which helps maintain high vector efficiency .
second , we perform 3d partitioning of the 3d matrix evenly among multiple threads , as shown in figure 3(b ) . for each partition ,
though fine - grain dynamic partitioning may provide better load balancing , it gives rise to interthread communication overhead .
in dynamic partitioning , a partition that is assigned to a thread at one stage may be assigned to another thread at another stage , which incurs data communication from the initial thread to the next . in our static partitioning ,
a partition is always assigned to the same thread throughout multiple stages , thus interthread communication is not required .
it is not trivial to vectorize and multithread fft due to its butterfly access patterns and bit - reversal element shuffling .
fft3d optimization in our cs implementation uses techniques discussed in , details of which are out of scope of this paper .
for example , figure 4 shows a high - level overview of stage3 which corresponds to assembly of the cost functional gradient defined in ( 7 ) following construction of its elements in stage1 and stage2 .
this stage is composed of five computation substages for two input matrices , a and b , and two output matrices , c and d. more specifically , a contains the set of all intermediary data generated using diff3d and b contains the set of all intermediary data generated using fft3d .
c then corresponds to the weighting matrix defined in ( 8) , whereas d is the entire composite variable in ( 7 ) .
one possible implementation is to execute each computation stage entirely before proceeding to the next state , as illustrated with the dotted arrows in the figure .
for example , we first multiply matrices a and b , then we perform an exponentiation of the result , and so forth . while this approach is easy to implement , its memory behavior is inefficient .
because each stage sweeps through the entire 3d matrix and the size of the matrix is usually larger than the last level cache size , temporary matrices between stages can not be retained within the cache , which results in cache thrashing , memory traffic increase , and , therefore , overall performance degradation .
a better implementation is to block the computation so that its temporary data is kept within the cache .
main memory accesses will occur only at the beginning to read the input matrix and at the end to write the output matrix .
we optimize even further to process the entire computation at the element level as shown in the solid arrow in the figure .
we read an element from each of matrix a and matrix b , perform the five computations , and write the result to matrix c and d. then , we move onto the next element , and so on .
this optimization is called loop fusion , because it fuses multiple small loops of individual computations into one big loop of a combined computation . because it handles one element at a time
, data can be kept in the registers , which eliminates the need for the temporary matrix and , therefore , removes intermediate memory loads / stores completely . also , because a fused loop performs more computation before it accesses the next element , it has more time to hide memory latency through data prefetches .
most matrix operations in cs read / write only one element from an input matrix to an output matrix .
however , diff3d requires 27 neighbor elements to compute an output element . in other words , an input element
is reused 27 times for 27 different output elements . to capture the data reuse , a cache - aware partitioning is required . for brevity , we explain our optimization with a 2d matrix shown in figure 5 .
it accesses elements from the beginning to the end of the current row before accessing elements in the the next row .
initially , it caches the first 3  3 inputs to calculate the first output .
but when it reaches the end of the row , the first 3  3 inputs are likely to be evicted from the cache if the number of elements touched during the row traversal exceeds the cache size . as a result , 3 
3 inputs are reloaded from memory to calculate the first output of the next row even though 2  3 of these inputs have been already read before .
it partitions the matrix in the middle of the row . instead of moving to the end of a row
, it stops at the end of the partition and moves down to the next row .
it can reuse the inputs from the previous row before they are replaced from the cache .
when we divide a matrix into four partitions , scheme  2 will show better cache behavior than scheme  1 .
we extend the same cache - aware data partitioning technique to a 3d matrix to implement diff3d .
though the main data structure in cs is a 3d matrix , the main memory access pattern is simple streaming that accesses data from the first to the last sequentially .
however , diff3d and fft3d also access data nonsequentially along the y and z axis .
memory accesses with a large power - of - two stride show poor cache behavior due to cache conflicts .
for example , this can occur when multiple data elements from adjacent matrix rows map to the same cache line . as the result , access to the second element results in the eviction of the first element from the cache . to solve the problem
each row along the x - axis is padded with one empty cache line at the end . without padding , accesses along the y - axis
we break the power - of - two stride by adding one extra cache line per row , which will reduce cache conflict misses .
note that , in addition to the padding at the end of row along the x - axis , we may also need to add another padding at the end of each xy 2d plane , if the size of y dimension is also power of two .
3d fft can be computed as multiple 1d ffts along the x- , y- and z - axis .
for a 256  160  80 matrix as our reference dataset , we can first perform 12800 256-point 1d ffts along the x - axis , followed by 20480 160-point 1d ffts along the y - axis , and finally 40960 80-point 1d ffts along the z - axis . however , performing 1d ffts , one axis at a time is not cache efficient , because it requires sweeping the entire matrix , which incurs a lot of cache misses due to the fact that the matrix does not fit into last - level cache .
instead , we perform 2d ffts for each 2d xy plane , then we perform 1d ffts for the z axis , shown in figure 7 . for a given z axis value
, we preload the corresponding xy plane to the cache , perform 1d ffts along the x - axis then the y axis ( a 2d fft for the entire xy plane ) , and then store the resulting xy plane to the memory . because last - level caches are usually larger than a single xy plane ( 320  kb for the reference dataset ) , our xy plane cache blocking is very effective in reducing memory bandwidth requirements . note that larger last - level caches in intel core i7 processor and intel 's knights ferry than tesla c2050 are beneficial for the 2d cache blocking . as each thread /
core works on a different xy plane , multiple 2d blocks should be kept in the cache .
in addition , as the size of datasets increases , the corresponding 2d blocks also get larger . therefore , to achieve good system - level and dataset - level scalability , large last - level caches are critical .
cs is composed of a large number of parallel stages separated by global barriers ( routines that synchronize threads in a parallel system ) .
when threads finish computation in the parallel region , they synchronize on a barrier before proceeding to the next parallel region .
efficient barrier synchronization is paramount for high scalability . on a system with a small number of cores such as intel core i7 processor ,
but as we increase the number of threads ( i.e. , 128 threads on intel 's knights ferry ) , we observe up to 10% overhead with our hand - optimized barrier implementation . without the optimized barrier
we implement a global barrier by launching a new kernel , that is , synchronizing with cpu , which costs one kernel launch overhead at minimum .
while there exist faster barrier implementations that run entirely on gpu , they require nonstandard memory consistency model assumption . another synchronization primitive used in cs is reduction . in stage3 and stage6 , it reduces a 3d matrix to one scalar value . to implement this reduction , we used a software privatization technique .
after all threads are done , the global reduction is performed , which aggregates all local values . to implement the global reduction
while atomics are generally slow on modern parallel architectures , their overhead in our cs implementation is small , due to the small fraction of time spent in the global reduction .
we compare performance of our cs implementation on three modern parallel architectures and provide in - depth performance analysis from a computer architecture perspective .
figure 8 shows the performance improvement of our cs implementation first as a single - threaded program on an intel core i7 processor as we incrementally applied our optimizations and subsequently as a multithreaded program .
the vertical bar represents the execution time in seconds for each optimization step , and the line shows the corresponding relative speedup over the baseline , base , which is the original single - threaded implementation of the algorithm compiled with the highest level of optimization , including autovectorization , function in - lining , and interprocedural optimization . as our first optimization , we replaced the fftw   used in the original implementation with the faster intel mkl .
second , we hand - vectorize the codes that can not be autovectorized by the compiler .
third , we apply cache blocking to exploit data reuse in the diff3d kernel , which shows another 1.10x speedup ( tile ) . through these three single - thread optimizations ,
we achieve an overall 1.51x speedup over the baseline implementation . to take advantage of multiple cores / threads , we parallelize the application . for fft3d
, we use the parallel implementation of the mkl library , and for the other kernels , we hand - parallelize using the openmp library .
parallelization achieves another 1.58x speedup on two cores over the single - core baseline and a 2.14x speedup on four cores .
overall , by combining the single - thread optimization and the multithread parallelization , we achieve a 3.21x performance improvement from the baseline implementation , which reduces the total execution time from 175 seconds ( base ) to 56 seconds ( 4 cores ) . 

figure 9 compares cs performance on three architectures : intel dual - socket six - core core i7 processor ( intel xeon processor x5670 at 2.93  ghz ) , nvidia tesla c2050 ( at 1.15  ghz attached to intel core i7 processor 960 at 3.2  ghz ) , and intel 's knights ferry ( at 1.2  ghz attached to intel core i7 processor 960 at 3.2  ghz ) . we normalize the speedups with respect to the optimized quad - core core i7 processor ( intel core i7 processor 975 at 3.33  ghz ) implementation ( 56 second runtime ) from the previous section and
since their dimensions are similar , the performance results for datasets 4 and 5 are correspondingly very similar to those for dataset 3 in terms of both execution time and relative performance across different architectures .
thus , for the sake of brevity , only the results for dataset 3 are shown in figure 9 .
for tesla c2050 and knights ferry , we show two speedup bars : one without data transfer overhead from the cpu host and the other with the overhead . the data transfer overhead results in small performance degradation , because cs spends significant time performing computation , and can hide most of the data transfer time .
the dual - socket six - core core i7 processor ( total 12 cores ) performs about 1.6x faster ( 35  s ) than the quad - core cpu ( 56  s ) , thanks to the increased core count and memory bandwidth .
the tesla c2050 gpu platform is 2.9x faster than the quad - core cpu for dataset 3 .
tesla c2050 shows about 7.5x speedup for dataset 1 but shows only 1.8x speedup for dataset 2 . for dataset 3 ( actual clinical data )
, knights ferry achieves 4.5x speedup over the quad - core cpu , which is about 1.57x faster than tesla c2050 .
knights ferry allows reconstructing the real anatomical dataset within a clinically feasible 12 seconds , which is a significant improvement over existing cs implementations .
note that core i7 processor and knights ferry are more efficient than tesla c2050 in terms of resource utilization . though tesla c2050 has about 4x peak flops and 3x peak bandwidth than core i7 processor , it only provides 2x performance .
also , while knights ferry delivers 10% flops / bandwidth of tesla c2050 , knights ferry shows about 55% speedup over tesla c2050 . finally , tesla c2050 's big performance variance across datasets is due to fft optimization .
we believe that cufft 3.1 is specially optimized for small power of two datasets like dataset 1 .
thus , its fft performance on dataset 1 is significantly better than on dataset 2 and 3 .
to elucidate the apparent impact architecture choice has on cs performance , we now provide an in - depth discussion on the performance of kernel - level operations .
we provide further insights into the achieved performance by breaking down the entire application into small microkernels and analyzing the microkernels individually .
we focus on the intel quad - core core i7 processor and intel 's knights ferry due to the lack of performance analysis tools in nvidia 's cuda environments .   table 1 shows the summary of our analysis .
the first column shows the fraction of execution time spent in the kernel inside the sequential code .
the second column shows the speedup achieved by knights ferry over core i7 processor on the kernel .
dataset 3 ( single - precision complex 256  160  80 ) is used for the analysis .
fft3d is the most important kernel occupying 46% of the total execution time . for the reference dataset ,
the intel mkl library achieves 30 gflops on core i7 processor and our in - house fft library achieves 175 gflops on knights ferry , which are the best performances that can be achieved on both architectures today .
therefore , for fft3d , we estimated 6x speedup of knights ferry over core i7 processor and actually obtained 4x speedup , which indicates that we might improve knights ferry performance further . 
 
it subtracts the original matrix and its shifted matrix by dx , dy , and dz : out = in - in_shifted ( dx , dy , dz ) . because the computation is a simple subtraction and the data size is large ( 200  mb ) , it is bandwidth bound . considering the memory bandwidth of core i7 processor and knights ferry , we expected 4x speedup of knights ferry over core i7 processor . in actual implementation
, we achieved 6x speedup , which indicates that core i7 processor may have room to improve . 
 
, exponentiation takes 37 cycles / element in core i7 processor and 0.88 cycles / element in knights ferry , and square - root computation takes 4.5 cycles / element in core i7 processor and 0.56 cycles / element in knights ferry .
based on the microbenchmark performance , we estimated 8x speedup of knights ferry over core i7 processor , and its actual performance ( 7x ) is close to our estimation . 
 
it is bandwidth bounded , because computations are simple addition , subtraction , or multiplication .
however , we also optimize it by fusing multiple computations together to exploit register and cache blocking . considering the peak memory bandwidth on both architectures ,
considering the peak floating point performance , knights ferry is expected to be 12x faster than core i7 processor . in reality
, we achieved 4x speedup of knights ferry over core i7 processor , which indicates that the kernel is currently bandwidth - bound on both architectures .  
chip multiprocessors ( cmps ) provide applications with an opportunity to achieve much higher performance , as the number of cores continues to grow over time in accordance with moore 's law .
for cmps to live up to their promise , it is important that as the number of cores continues to grow the performance of the application running on cmps also increases commensurately . to gain insights into how the cs algorithm will scale on future many - core architectures , we modeled a feasible but hypothetical future cmp processor ( running at 1  ghz with 256  kb / core cache ) on our cycle - accurate research simulator and figure 10 shows its cs performance scalability .
in particular , our 64-core configuration achieves 37x speedup over a single - core configuration , which is almost 60% parallel efficiency .
moving further , our 128-core configuration achieves 55x speedup over single - core configuration which is little less than 50% efficiency .  
table 2 compiles the predicted 128-core acceleration together with the previously discussed results for intel core i7 processor , nvidia tesla c2050 , and intel 's knights ferry .
although a 256  160  80 data volume may be considered large by numerical computing standards , it is relatively small by modern clinical standards and volumes many times larger are routinely encountered in practice .
moreover , many contemporary acquisition trends are migrating from single - phase to time - resolved paradigms , where a 3d  movie  of dynamic anatomy and physiology ( e.g. , of contrast flowing into and out of vessels ) is created .
thus , the scalability of cmps is paramount to seeing cs - type and other nonlinear reconstruction methods become practical for such imaging scenarios .
overall , existing and future many - core architectures are very promising platforms for accelerating the cs reconstruction algorithm to make it clinically viable .
in this work , we have shown that advanced computing architectures can facilitate significant improvements in the performance of cs mri reconstructions and particularly that optimized use of modern many - core architectures can significantly diminish the computational barrier associated with this class of techniques .
this suggests that as many - core architectures continue to evolve , cs methods can be employed in routine clinical mri practice .
although ce - mra was targeted in this work , the implication of the results apply to many other mri applications as well as other areas in medical imaging such as dose reduction in computed tomography .
 compressive sensing ( cs ) describes how sparse 
 signals can be accurately reconstructed from many fewer samples 
 than required by the nyquist criterion . 
 since mri scan duration 
 is proportional to the number of acquired samples , cs has been 
 gaining significant attention in mri . 
 however , the computationally 
 intensive nature of cs reconstructions has precluded their 
 use in routine clinical practice . in this work , 
 we investigate how 
 different throughput - oriented architectures can benefit one cs 
 algorithm and what levels of acceleration are feasible on different 
 modern platforms . 
 we demonstrate that a cuda - based code 
 running on an nvidia tesla c2050 gpu can reconstruct a 
 256  160  80 volume from an 8-channel acquisition in 19 seconds , 
 which is in itself a significant improvement over the state of the art . 
 we then 
 show that intel 's knights ferry can perform the same 3d mri 
 reconstruction in only 12 seconds , bringing cs methods even 
 closer to clinical viability . 
1. Introduction and Motivation
2. Methods
3. Experiments
4. Architecture-Aware Optimization
5. Results
6. Discussion
7. Conclusion