CG may be numerically suboptimal and seems to be allocating

[JeffFessler]

The conjugate gradient (CG) routine here looks to be written to "solve" the system of equations H*x=b
where H is positive (semi) definite, rather than minimizing the least-squares cost function | A x - y |_2^2.

https://github.com/guanhuaw/MIRTorch/blob/master/mirtorch/alg/cg.py

On paper the two are the same by letting H = A'A and b = A'y, but computing the gradient like
H x - b = (A'Ax) - A'y
is less numerically accurate than working with the "residual" form A' (A x - y)
because A'A has the square of the condition number of A.

Of course, for the "Toeplitz" version in MRI, we have to use T=A'A.
But this toolbox is more general than that special case.

Also the cg.py code looks like it has a lot of (presumably) allocating operations like r0 = b - A * x0.

The Astra toolbox implementation seems to use the residual form A' (A x - y)
and also does lots of in-place operations like -= which presumably would save time.
https://github.com/astra-toolbox/astra-toolbox/blob/master/python/astra/plugins/cgls.py

FWIW, the stable Julia version in MIRT.jl uses the residual form but does lots of allocations.
https://github.com/JeffFessler/MIRT.jl/blob/main/src/algorithm/general/ncg.jl
There is a WIP to make an in-place version:
https://github.com/JeffFessler/MIRT.jl/blob/ncg/src/algorithm/general/ncg2.jl

[guanhuaw]

Thanks @JeffFessler! You made a good point about the numerical precision, especially for fp16/bf16 scenarios. I am thinking about how to make the 'residual form' backward compatible.

I realized that there are unnecessary calculations of A*x in the code, as well as memory allocations. Will fix this first.