SPDHG compute norms

[epapoutsellis]

In SPDHG, if we do not precompute norms outside the algorithm, we need to iterate over all operators and append the norm for each operator into a list.

https://github.com/TomographicImaging/CIL/blob/45e0b235ac79efc03dc5cb46b8ddda9d16f6331b/Wrappers/Python/cil/optimisation/algorithms/SPDHG.py#L153-L161

If the number of operators, i.e., number of subsets is small, this is ok. However, it will take some time to setup SPDHG if we have a large number of subsets and large datasets.

In the case where the number of subsets divides the total number of projection angles, we can check if the corresponding geometries for each of these operators is the same and compute only one norm of these operators.

all_geoms = [i.range_geometry().shape for i in operator]

if all(x == all_geoms[0] for x in all_geoms):
    norms = [operator[0].norm()]*self.ndual_subsets

I can create a PR if you are interested.

[MargaretDuff]

@epapoutsellis

Thanks for suggesting this! We had a chat about maybe incorporating this in the block operator norm, which currently operator, in turn, calculates the norm and squares and sums the result. If instead, for each new operator, there could be some check, which perhaps compares the geometry, to see if the norm had been calculated before it could save the norms being recalculated. The question is then how could we tell, without calculating the norm, that the norm should be the same?

I loaded in some example data:

Number of dimensions: 3 Shape: (128, 300, 128) Axis labels: ('vertical', 'angle', 'horizontal')

Partitioned the data into 20 subsets:

Produced a blockoperator from the newly partitioned data

For each block calculated the norm:

As you can see, the norms are almost equal but not exactly:

So even with the same number of angles in a partition, we get slightly different norms. Is this what we expect? Is it an error with our norm calculation? Is it a small enough difference that, in practice, we can treat the norms to be equal?

[MargaretDuff]

So I have had a play with a possible (messy) implementation that can be found here in the BlockOperator object: https://github.com/MargaretDuff/CIL-margaret/blob/257ff1112430360ca0115d79cfe1f4bdc121e9e3/Wrappers/Python/cil/optimisation/operators/BlockOperator.py#L147-L182 Here if the flag, "block_norms=True" is called when you set up the block projector or block operator the norm will be calculated based on the geometries. If the geometry of each block is the same, the norms will be calculated based on the first block. If they differ, for each new block the code will check it has not seen the geometry already before using a pre-saved value or recalculating. The flag "block_norms=False", results in the same behaviour as before.

For the (128,300,128) data above, split into 10 equal subsets, on Edo's VM calculating the norm takes 1.25seconds using the previous defaults and this reduces to about 0.15 seconds. For 11 non-equal subsets it goes from 1.3 to 0.28 seconds. The calculated norms are slightly different, but close, e.g. for the 11 subsets calculating each subset individually (the default) gives 3066.8028451348337 and for the new implementation 3066.338020925357. Code for testing can be found https://github.com/MargaretDuff/cil_test_notebooks/blob/master/block_operator_norm/test_norms.py