Enable TV, L1-wavelet and other sparse regularisation in LSQR and CGLS, via a new IRLS algorithm and Tikhonov operators by bryceshirley · Pull Request #2369 · TomographicImaging/CIL

bryceshirley · 2026-09-07T18:42:21Z

Description

Adds an IRLS (Iteratively Reweighted Least Squares) outer algorithm that solves the sparsity-promoting problem

$$\min_u |Au-b|_2^2 + \alpha^2|Lu|_1$$

as a sequence of reweighted least-squares problems, each handled by an inner LSQR or CGLS. To support it, LSQR and CGLS gain structured Tikhonov regularisation in block and standard form.

See the notebook regularisation-for-krylov-methods.ipynb‎ in CIL-Demos/misc pull request 292 for example usage.

Contribution Notes

The content of this Pull Request (the Contribution) is intentionally submitted for inclusion in CIL (the Work) under the terms and conditions of theApache-2.0 License
I confirm that the contribution does not violate any intellectual property rights of third parties

Checklist

I have performed a self-review of my code
I have added docstrings in line with the guidance in the developer guide
I have updated the relevant documentation
I have implemented unit tests that cover any new or modified functionality
CHANGELOG.md has been updated with any functionality change
Request review from all relevant developers

gfardell · 2026-09-09T09:00:36Z

This looks fantastic @bryceshirley, thank you! Would you be able to join a developer meeting to walk us through it (I think the 22nd is our next one)? I know you've already spoken to/worked with the team, but it would be good to go through the current status together.

bryceshirley · 2026-09-09T14:03:08Z

This looks fantastic @bryceshirley, thank you! Would you be able to join a developer meeting to walk us through it (I think the 22nd is our next one)? I know you've already spoken to/worked with the team, but it would be good to go through the current status together.

Sounds good! I can present that pull request and explain the decisions in terms of interface and memory management, as well as alternatives. How much time should it be?

MargaretDuff

I have left a bunch of comments, none of them very substantial. A lot of it was deleting unnecessary comments.

What I suggest, and would like to discuss at the developer meeting next week, is a change to the user interaction with IRLS. Currently, to use IRLS the user must pass a set up CGLS or LSQR algorithm, with the right operators, the right choice of the weighted flag in CGLS or LSQR. This confuses the standalone CGLS or LSQR algorithm with parameters that might not be relevant and makes things like the form parameter difficult to explain. It would also allow us to hide most of the new create_tikhonov_operator, resolve_form, WeightedStructOpeartor, BlockTikhonovOperator - they become implementation details rather than things a user must understand to call IRLS.

On the other hand, it looses some flexibility - it would be harder to pass a custom inner solver. Though possibly later on we could have a IRLS.from_solver(inner_solver,..) alternative set up later?

For IRLS it would allow us to have a safe default inner solver - LSQR is better conditions.

MargaretDuff

I have left a bunch of comments, none of them very substantial. A lot of it was deleting unnecessary comments.

What I suggest, and would like to discuss at the developer meeting next week, is a change to the user interaction with IRLS. Currently, to use IRLS the user must pass a set up CGLS or LSQR algorithm, with the right operators, the right choice of the weighted flag in CGLS or LSQR. This confuses the standalone CGLS or LSQR algorithm with parameters that might not be relevant and makes things like the form parameter difficult to explain. It would also allow us to hide most of the new create_tikhonov_operator, resolve_form, WeightedStructOpeartor, BlockTikhonovOperator - they become implementation details rather than things a user must understand to call IRLS.

On the other hand, it looses some flexibility - it would be harder to pass a custom inner solver. Though possibly later on we could have a IRLS.from_solver(inner_solver,..) alternative set up later?

For IRLS it would allow us to have a safe default inner solver - LSQR is better conditions.

Co-authored-by: Margaret Duff <43645617+MargaretDuff@users.noreply.github.com>

Signed-off-by: Bryce Shirley <114733577+bryceshirley@users.noreply.github.com>

bryceshirley · 2026-10-14T11:01:13Z

I have updated the PR to now make IRLS a stand alone algorithm with default inner loop solver LSQR in block form and moved the demo notebook into a CIL-Demos pull request.

MargaretDuff

A quick review of the CGLS algorithm. Will try and keep looking this afternoon

MargaretDuff · 2026-10-15T12:30:54Z

+        # Catch self.regalpha and the operator's regalpha disagreeing.
+        operator_alpha = getattr(self.operator, 'regalpha', None)
+        if operator_alpha is not None and operator_alpha != self.regalpha:
+            raise ValueError(


I am not sure a user would understand what this means if they see this error! When would you expect it to come up?

MargaretDuff · 2026-10-15T12:33:31Z

There are a lot of if block and if standard statements in this file. I wonder if instead we use a factory method? Could the user call CGLS and then inside that file it either sets up and returns classes _standardCGLS or _blockCGLS ? Or perhaps a standard CGLS that factories off to _blockCGLS when needed but that _blockCGLS also uses things from CGLS?

Realise we don't do this elsewhere in our algorithm class - @gfardell or @paskino would appreciate if you had any thoughts?

Co-authored-by: Margaret Duff <43645617+MargaretDuff@users.noreply.github.com> Signed-off-by: Bryce Shirley <114733577+bryceshirley@users.noreply.github.com>

MargaretDuff

A pass at IRLS - Thanks Bryce!

MargaretDuff · 2026-10-15T14:35:24Z

+
+    Allowing the inner solver to solve
+    .. math::
+        \min_x \|Ax - b\|_2^2 + \alpha^2 \|Lx\|_1


Can we put this as the first equation on the algorithm? So,

This equation solves

It does this by

Any notes

Parameters

MargaretDuff · 2026-10-15T14:38:10Z

+    ----------
+    initial : DataContainer, optional
+        Starting point of the first inner solve, and the point the first
+        weights are computed from. Defaults to zeros.


Maybe not a code question but for when we recommend this to users - are the first weights defaulted zeros? Would initialising the weights with a FBP recon be a good idea? If we put in bad initial weights, could the algorithm recover?

MargaretDuff · 2026-10-15T14:41:36Z

+    struct_operator : Operator, optional
+        The structural operator :math:`L`. ``None`` means the identity,
+        i.e. a plain L1 penalty on the solution.
+    inner_solver : class, default LSQR


Do you pass the class or pass a string? Could we pass both? I.e. a string either LSQR or CGLS or a class with these properties...

I think our basic users would prefer just a string but our mathematical users might appreciate being able to pass their own custom solver.

MargaretDuff · 2026-10-15T14:53:30Z

+
+        # Solution buffer for standard form's (WL)^-1 mapping; block form uses the live iterate.
+        self.tmp_solution = (
+            self.inner_solver.solution_geometry().allocate(0)


Note that this is probably not SIRF compatible

MargaretDuff · 2026-10-15T14:58:28Z

+        u = self.inner_solver.get_output()
+
+        residual = self.inner_solver.operator.operator.direct(u)
+        residual.sapyb(1.0, self.inner_solver.data, -1.0, out=residual)
+        data_term = residual.squared_norm()


Any way to reuse memory for this and not reallocate?

MargaretDuff · 2026-10-16T18:19:58Z

@lauramurgatroyd @gfardell @paskino - Just to say that I ran out of time looking more at this. It would be great if you have time to give Bryce some feedback! Thank-you!

bryceshirley marked this pull request as draft September 7, 2026 18:42

lauramurgatroyd added this to CIL work Sep 9, 2026

github-project-automation Bot moved this to Todo in CIL work Sep 9, 2026

MargaretDuff self-assigned this Sep 16, 2026

MargaretDuff reviewed Sep 16, 2026

View reviewed changes

lauramurgatroyd requested a review from gfardell September 22, 2026 13:43

Add BlockDataContainer diagonal support to DiagonalOperator

01c38d2

Co-authored-by: Margaret Duff <43645617+MargaretDuff@users.noreply.github.com>

bryceshirley force-pushed the pull-request-1-refactor branch from 7fd6acb to aebf28e Compare October 13, 2026 19:28

bryceshirley changed the title ~~DRAFT: Enable TV, L1-wavelet and other sparse regularisation in LSQR and CGLS, via a new IRLS algorithm and Tikhonov operators~~ Enable TV, L1-wavelet and other sparse regularisation in LSQR and CGLS, via a new IRLS algorithm and Tikhonov operators Oct 14, 2026

bryceshirley mentioned this pull request Oct 14, 2026

Demonstration for CIL pull request #2369 on regularisation in Krylov algorithms (IRLS, LSQR and CGLS) TomographicImaging/CIL-Demos#292

Open

bryceshirley force-pushed the pull-request-1-refactor branch 2 times, most recently from 68e9097 to 172c647 Compare October 14, 2026 10:47

bryceshirley marked this pull request as ready for review October 14, 2026 10:50

bryceshirley and others added 3 commits October 14, 2026 10:57

Add generalised Tikhonov regularisation to LSQR and CGLS

8cd86e3

Add IRLS algorithm for L1-regularised least squares

0653849

Merge branch 'master' into pull-request-1-refactor

050128e

Signed-off-by: Bryce Shirley <114733577+bryceshirley@users.noreply.github.com>

bryceshirley force-pushed the pull-request-1-refactor branch from 172c647 to 050128e Compare October 14, 2026 10:58

MargaretDuff reviewed Oct 15, 2026

View reviewed changes

Apply batched suggestions from code review

00a848f

Co-authored-by: Margaret Duff <43645617+MargaretDuff@users.noreply.github.com> Signed-off-by: Bryce Shirley <114733577+bryceshirley@users.noreply.github.com>

MargaretDuff reviewed Oct 15, 2026

View reviewed changes

fixing comments

e7f2821

bryceshirley force-pushed the pull-request-1-refactor branch from 61b70d8 to e7f2821 Compare October 28, 2026 20:21

Merge branch 'TomographicImaging:master' into pull-request-1-refactor

a0ef82b

Conversation

bryceshirley commented Sep 7, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Description

Contribution Notes

Checklist

Uh oh!

gfardell commented Sep 9, 2026

Uh oh!

bryceshirley commented Sep 9, 2026

Uh oh!

MargaretDuff left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

MargaretDuff left a comment

Choose a reason for hiding this comment

Uh oh!

bryceshirley commented Oct 14, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

MargaretDuff left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

MargaretDuff Oct 15, 2026

Choose a reason for hiding this comment

Uh oh!

MargaretDuff Oct 15, 2026

Choose a reason for hiding this comment

Uh oh!

MargaretDuff left a comment

Choose a reason for hiding this comment

Uh oh!

MargaretDuff Oct 15, 2026

Choose a reason for hiding this comment

Uh oh!

MargaretDuff Oct 15, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

MargaretDuff Oct 15, 2026

Choose a reason for hiding this comment

Uh oh!

MargaretDuff Oct 15, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

MargaretDuff Oct 15, 2026

Choose a reason for hiding this comment

Uh oh!

MargaretDuff commented Oct 16, 2026

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

bryceshirley commented Sep 7, 2026 •

edited

Loading

bryceshirley commented Oct 14, 2026 •

edited

Loading