Duplicate implementation of power iteration

[migrosser]

I just stumbled upon the fact that we have one method power_iterations! in MRIReco (which in turn is used by espirit) and a similar method power_iterations in RegularizedLeastSquares. The only difference between these methods seem to be the default values for rtoland maxiter and the fact that power_iterations! allows to pass a pre-allocated vector for b. I think it would be nice to adapt the implementation in RegularizedLeastSquares so that it can be used in espirit and then remove the implementation in MRIReco. I think it should not be a great deal of work and we would avoid this little bit of duplicate code. What do you think about this @JakobAsslaender?

[JakobAsslaender]

Hmm, I agree with your general intend, but it would further increase the dependencies? I would still prefer to move all coil estimation code out of MRIReco.jl into a dedicated, light weight package and in this case it would be nice to not introduce a dependency for this simple function. But I am interpreting @tknopp's silence on this topic as a preference of his to keep the code in MRIReco.jl (which is of course fair enough). Am I interpreting this correctly, @tknopp?

[tknopp]

I am actually not sure what the right modularization approach is here. Would you really see MRICoilSensitivities to be completely independent of any image reconstruction algorithm?

[JakobAsslaender]

Yes, I would. This way, you could easily use it any other package w/o heavy dependencies.

Just to give you one (egocentric) use example: I am currently setting up a MRFingerprintingRecon.jl package that will eventually implement a low rank ADMM algorithm. I think this is too specific to fit into something like MRIReco.jl, but, of course, it heavily builds on RegularizedLeastSquares.jl, NFFT.jl etc. And I include MRIReco.jl, purely for ESPIRIT. From my point of view, a package like MRICoilSensitivities.jl provides tools on the same abstraction level as, e.g. NFFT.jl and such a setup would streamline the implementation of specific image reconstruction packages.

[JeffFessler]

I too would like the convenience of a separate, smaller package for using ESPIRIT or such. A comment though is that estimateCoilSensitivities here works with input images, but some coil map estimation methods are joint with recon methods, like JSENSE: https://sigpy.readthedocs.io/en/latest/generated/sigpy.mri.app.JsenseRecon.html#sigpy.mri.app.JsenseRecon. So maybe the repo here should be MRICoilSensitivitiesFromImages. I'm kind of kidding there, but perhaps we should think about whether JSENSE would belong in MRICoilSensitivities or elsewhere like MRIreco.

[tknopp]

Jeff made my point much more clear. If we put JSENSE into MRICoilSensitivities, the package would be a dependent of MRIReco. Otherwise it would be a dependency. That makes a huge difference.

But I am all for proper modularization so it would be very helpful, if you two could give me some sketches on how your ideal package structure (as a tree) looks like. Then we can decide if we should name it MRICoilSensitivities or probably just ESPIRIT.jl.

@JakobAsslaender: It would be awesome if we could group all MR packages under the same Github organisation (MagneticResonanceImaging). You might consider moving MRFingerprintingRecon.jl. This would give us a much more clear picture what Julia packages we have and how we should group them properly. Same goes for BART.jl as I recently said. I feel that this gives MRIReco much more momentum. Please have a look at JuliaImages https://github.com/JuliaImages/Images.jl for my vision. They have joint documentation website where all the subpackages are explained.

[JeffFessler]

Thinking more, perhaps I muddied the water too much. JSENSE is a recon method that happens to also produce coil sensitivities, so arguably it belongs in MRIreco. I would prefer MRICoilSensitivities over ESPIRIT.jl both to better comply with julia package naming conventions and also to make room for future alternatives that also work in the image domain. It should suffice to have the readme explain that the repo is just for image domain methods and joint recon methods are elsewhere...

[JakobAsslaender]

I agree with @JeffFessler, JSENSE is foremost an image reconstruction and should be part of MRIReco.jl. The implementation of ESPIRIT, with a main function that takes an array of numbers as input. The wrapper function that takes the MRIReco.jl data format would stay in MRIReco.jl, and if someone were to introduce a new data format for whatever reason, they could write their own wrapper in their package.

And yes, I will move the packages. I just need to figure out how to do that w/o breaking the registry etc.

[tknopp]

I will have a look into this when I find some time. MRICoilSensitivities will certainly be based on MRIBase to keep all coil related functionality in a single package. If we do the modularization, I want to do this properly and don't want to write many wrappers in MRIReco.jl. Ideally MRIReco.jl is just a meta-package which reexports a certain package collection (of which MRFingerprintingRecon.jl for instance might be one).

[tknopp]

There is, by the way, yet another power iterations implementation :-)

https://github.com/JeffFessler/MIRT.jl/blob/26cfbc2a26b26814fa85739dbe01b5f1b8be5e21/src/algorithm/general/poweriter.jl

I would like to keep the duplicated code for the moment and first do the MRICoilSensitivities split off. But we keep this issue open as a reminder that there is room for improvement.

[tknopp]

I think I found a 4th implementation :-)

https://iterativesolvers.julialinearalgebra.org/stable/eigenproblems/power_method/

I would argue that we should just use that and remove all other implementations. IterativeSolvers is well established.

[tknopp]

I made some small tests and it seems the method in IterativeSolvers.jl is slower and allocating. Not sure if I tested properly though. In this commit one can see my (commented out) test:

https://github.com/MagneticResonanceImaging/MRIReco.jl/commit/a8407b2b1610285e244a3774a785f04ae7d89113

If there is a performance issue in IterativeSolvers, we should just fix that there and afterwards remove our implementation.