Positivity Constraint Convergence and Relaxation Factor

[XiaoWang-Github]

The current C code reduces the number of parallel threads to stabilize convergence when the positivity constraint is turned off. If we implement the relaxation factor for voxel update, then there will be no need to reduce the parallel threads. I tested the C code with relaxation factor set to be 0.8 on NERSC cori cluster. That stabilizes the algorithm convergence.

[cabouman]

Xiao, I'm concerned that a single fixed relaxation factor of 0.8 could slow convergence on regular computers. Have you you done any testing across different platforms?

[XiaoWang-Github]

Hi Professor Bouman, we can automate how we choose the relaxation factor. I sent you, Jordan and Professor Buzzard an email this morning with presentation slides from 2016. We can approximate Jun Zheng's 1999 paper equation for the relaxation factor, and it can be approximated as Nv/(Nv+nc-1), where Nv is the number of view angles, nc is the number of parallel threads.

[sjkisner]

Practically, the thread limitation that happens when the positivity constraint is off will only have an effect if the image volume is small (on any given node).

I want to be very cautious about turning off any of the fixes described toward the bottom of Issue #212. Since then we haven't had any reports of divergent behavior.

But I'm strongly in favor of adding over/under relaxation. In a lot of problems this can accelerate convergence.

[sjkisner]

Also a technical point, a source of the instability we saw with positivity=False is a time lag between the pixel update and the corresponding update to the error sinogram due to the double buffering. I have an idea on how to tighten that up but it's not implemented yet.

[cabouman]

The formula I get is:

[XiaoWang-Github]

Fantastic!! It's simple to compute and makes a lot of sense to me.

[XiaoWang-Github]

Intuitive too.

[sjkisner]

Ran an experiment with the relaxation formula above with unfavorable results. I don't think we should set a default relaxation to anything but 1.0.

Based on demo_2D_shepp_logan.py :

threads=1:
relax_factor: 1.0
iterations to converge: 25
recon time: 0.813 sec

threads=2:
relax_factor: 0.41
iterations to converge: 52
recon time: 1.39 sec

threads=20:
relax_factor: 0.04
iterations to converge: 163
recon time: 4.41 sec

[sjkisner]

Relaxation factor added to recon() arguments in PR #243

[XiaoWang-Github]

Hi Jordan, did you use the equation that Professor Bouman and I proposed to calculate the relaxation factor? The relaxation factor is usually like 0.8 or 0.9, and it shouldn't go below 0.5. And it's most useful for datasets that have phase contrast features and negative voxels. Then you don't need to restrict the number of parallel threads like what the current SV-MBIR implementation does.

[sjkisner]

Yes in demo_2D_shepp_logan.py, Np=256x256, Nsv=19x19, so for nc=20 you get
256/(256+19^3) = 0.036.

[cabouman]

Jordan and Xiao, that formula is probably much too conservative. The problem is that it is based on an old formula which assumes individual greedy update of every voxel. But the voxels in an SV are not full individually greedy because they depend on each other. It might be possible to derive a less conservative bound.

Perhaps the three of us can meet to discuss this?

Intuitively, something like this might work:

[XiaoWang-Github]

Yes. That looks better. We can meet.
In the previous version of the formula, we had Nsv *Nc by assuming that all voxels in every super-voxel is updated in parallel. That's probably not correct. The super-voxels are updated in parallel, but the voxels inside the super-voxels are updated in sequence. This version of the equation acknowledges this hybrid parallel-sequential updates by using square root, and that makes a lot of sense to me.

[cabouman]

Yes, exactly. Maybe we can derive a theory to explain this based on block updates of pixels.