Repository navigation
Title: Krylov methods: float64 norm accumulation, working restarts, and results that were being discarded - #774
Merged
Conversation
Two defects that only bite at production sizes, found while explaining why
720-view CGLS at 1024^3 always exited after one iteration reporting
divergence, while the same code at 45-180 views ran to completion.
1. np.linalg.norm sums squares in the array's OWN dtype. On a
projection-sized float32 array the accumulator dwarfs the increments still
being added to it and the norm comes back LOW - measured against a chunked
float64 reference on 3072^2 detector data:
45 views (4.2e8 elements) 0.4 % low
180 views (1.7e9 elements) 3 % low
360 views (3.4e9 elements) 16 % low
720 views (6.8e9 elements) 26 % low
This is not a reporting detail. alpha = gamma / q_norm**2, so a q_norm 26 %
low inflates the CGLS step by ~1.8x; the algorithm overshoots, the residual
rises, and the loss-of-orthogonality test - itself a comparison of two of
these norms - reads the overshoot as divergence. A 1024^3 volume (1.1e9
elements) is already in the bad regime. Every norm in the module now routes
through a float64 accumulator, returning float32 so no caller's dtype
changes (Ax/Atb reject float64).
The MATLAB reference anticipates exactly this in a comment: "may be caused
either by the mismatch of the backprojection w.r.t the real adjoint, or
numerical issues related to doing several order of magnitude difference
operations on single precission numbers." Measurement picks the second: the
pair is already Ax 'Siddon' / Atb 'matched', which an adjoint inner-product
test puts at ~4e-5 agreement.
2. When orthogonality IS lost, MATLAB's CGLS.m breaks out of an inner loop
into an outer one that rebuilds r/p/gamma and carries on. The Python port
kept only the inner loop, so that break left run_main_iter altogether -
there was no restart, only an exit. Its give-up test also read
`re_init_at_iteration + 1 == i` against an initial 0, true at i == 1, where
MATLAB compares `remember == iter` against an initial 0 its 1-based counter
cannot hit on the first check. So the very first lost step ended the run.
Restored the two-loop structure and the sentinel (-1).
The sibling Krylov algorithms (LSQR, LSMR, hybrid_LSQR, AB/BA-GMRES,
hybrid_fLSQR_TV) share both the flattened loop and the sentinel. They get the
norm fix, which is mechanical; the loop restructure is deliberately left for a
change that can test each one, because fixing the sentinel alone would turn
their loud early exit into a silent one.
Tests: Python/tests/test_krylov_restart_and_norms.py (5, no GPU - the numerics
are pure numpy and the control-flow tests drive run_main_iter with scripted
residual norms).
The float32 accumulation fixed for the Krylov family in 5aac944 is not
confined to it. im3DNORM(x, 2) is a bare np.linalg.norm, and ASD-POCS drives
its control flow with it:
dd = im3DNORM(g - self.proj, 2) # projection-sized: 26 % low at 720v
dp = im3DNORM(dp_vec, 2) # volume-sized: ~1 % low at 1024^3
dg = im3DNORM(dg_vec, 2)
`dd` gates the data-fit constraint against epsilon and the dp/dg ratio gates
the TV step, so at full detector size and high view counts the algorithm is
steering on numbers that are wrong by tens of percent - not merely reporting
them wrongly. iterative_recon_alg's error_measurement uses it as well.
l2norm() now lives in im3Dnorm.py as the single implementation, im3DNORM
routes its L2 case through it (other norms untouched), and the Krylov module's
_norm delegates rather than keeping a second copy.
Note for anyone comparing results across this commit: POCS runs at 360 and 720
views were measured on the old norms. At 45-90 views the projection arrays are
small enough (0.4 % error) that the difference is immaterial.
result-discarding bugs found doing it
Follow-on from 5aac944 / 5f0a845, which fixed CGLS's norms and im3DNORM.
Auditing the rest of the module turned up more of the same class, plus
defects that made some algorithms return the wrong array entirely.
WHICH REDUCTIONS ARE ACTUALLY DANGEROUS (measured, not assumed)
np.sum is fine - numpy's pairwise summation gives 4e-9 relative error at
1024^3. Only BLAS-dispatched reductions degrade: np.linalg.norm and np.dot.
The mechanism is an accumulator growing against terms that stay small, so a
sum of SQUARES (all positive, nothing cancels) is far worse than a dot of two
different vectors. At 1024^3 on the Linux box: norm 1.4e-2, dot(x,x) 2.9e-2,
dot(x,y) 7.8e-5, np.sum(x*x) 4.4e-9. It is also BLAS-dependent - the same
sizes on this project's Windows BLAS are ~76x more accurate - so no platform
can be assumed safe.
Fixed accordingly (l2norm for sums of squares, inner for dot products, both
float64-accumulating and both returning float32 so no dtype changes):
* power_method.svd_power_method - all three norms, on projection- and
volume-sized arrays. This IS the Lipschitz constant.
* ista_algorithms - the FISTA auto-hyper power iteration. step = 1/hyper,
so an error here is an error in every FISTA step.
* pocs_algorithms - the cosine gating the TV step reduction (both ASD_POCS
and AwASD_POCS).
* hybrid_fLSQR_TV - the Golub-Kahan reorthogonalisation coefficients.
* AB/BA-GMRES - see below.
Left alone deliberately: np.sum(...) everywhere, and one np.dot(scalar,
vector) in hybrid_fLSQR_TV that is not an inner product.
AB/BA-GMRES ORTHOGONALISED WITH PYTHON'S BUILTIN sum()
`h[k,i] = sum(qk.ravel()*w[i])` iterates a CT-sized array element by element.
Measured 40x slower than np.dot and 800x less accurate. BA-GMRES's basis is
volume-space, so at 1024^3 that is ~107 s per coefficient and 210 coefficients
for niter=20 - which is most of the 4.1 h it took in the full-resolution
sweep. Its reputation as "host-bandwidth-bound" needs revisiting: the bandwidth
was never the binding constraint. The inner update also reallocated a
full-size array every step; it is in-place now.
THREE BUGS THAT DISCARDED RESULTS
decorator() ignores run_main_iter's return value and hands back self.res, so
`return <expr>` without assigning self.res silently returns the INITIAL volume
- zeros, or with the default FDK warm start an unimproved FDK image that looks
entirely plausible. Confirmed on a 64^3 phantom: ba_gmres and hybrid_lsqr
returned all-zero volumes whenever they exited early. hybrid_fLSQR_TV never
assigned self.res on ANY path, so its output was discarded outright. All exit
paths now assign. (Also removed a debug print(y[0]) from its inner loop.)
RESTART, EVERYWHERE IT WAS BROKEN
LSQR and LSMR had CGLS's exact defect - inner loop only, so the restart
`break` left run_main_iter, and the sentinel fired at i == 1. Both now have
the two-loop structure and the -1 sentinel. hybrid_LSQR, hybrid_fLSQR_TV and
AB/BA-GMRES have no restart path at all; their sentinel test could only fire
at one particular iteration and ignored a residual rise at every other one, so
they now simply stop on a genuine rise. On the phantom that is also the better
answer: ba_gmres stopping at iteration 5 scores 0.296 against 0.216 for
running to 20.
TESTS
tests/test_krylov_smoke.py (8) - each algorithm end to end on a 64^3 phantom:
finite, non-zero, and correlating at least 0.6x as well as FDK on the same
data, plus a direct check that CGLS fills more than 2 of 20 residual slots and
a source guard against builtin sum() returning. test_krylov_restart_and_norms
(8) gains inner() accuracy and a structural check that LSQR/LSMR/CGLS all
carry the two-loop fix. 20 pass on Windows/sm_86.
Not fixed: irn_tv_cgls still diverges on the phantom (corr 0.06, values 1e5),
reproducing the known 512^3 behaviour - its inner loop has no divergence check
and its lmbda=1 TV weight is untuned. That is a separate problem from these.
It diverged at 512^3 (mean 8469, max 5e7) and reproduces on a 64^3 phantom in
seconds - correlation 0.06 against FDK's 0.42, voxel values reaching 7e5. Two
causes, both fixed here.
1. L AND L^T WERE NOT AN ADJOINT PAIR (measured asymmetry 6.4e-2)
CGLS is applied to the stacked operator [A; sqrt(lmbda) W D], and it assumes
the transpose it is handed IS the adjoint. This one was wrong twice over:
* `Dxx = np.copy(img)` followed by writing only `[0:-2]` left the LAST TWO
slices holding raw image VALUES rather than differences. So D did not
annihilate a constant image - it penalised intensity there, a Tikhonov
term smuggled into part of the volume - and `D(1) != 0` is trivially
checkable, which is now a test.
* D^T was off by one at index n-2, leaving `Wx[n-2]` unreduced.
Replaced with a derived pair, D u[i] = u[i]-u[i+1] (0 at i=n-1) and its exact
transpose, applied per axis via swapaxes so the three directions share one
implementation. Asymmetry is now 0.0e0 and D(constant) is exactly zero.
2. THE INNER CGLS HAD NO DIVERGENCE GUARD
Plain CGLS has always had one; this loop ran niter_outer * niter unguarded
steps. On the phantom the residual starts rising at step 2 and compounds:
9.7e3 -> 6.9e5 over 20 steps. The guard undoes the bad step and falls out to
the next outer reweighting, which is this algorithm's natural restart.
The comparison is scoped WITHIN one outer iteration on purpose. Across a
reweighting the two residuals belong to different weighted problems, so a
first step that looks worse than the previous outer's last is not evidence of
anything - a first attempt that compared across the boundary stopped the
solver after two productive steps. A separate check stops the run when a whole
reweighting fails to improve the best residual, i.e. the outer loop stalled.
Also removed `self.res = res0` at the top of the outer loop: res0 is
reassigned inside the inner loop, so it restored the second-to-last iterate
rather than any meaningful starting point.
RESULT on the 64^3 phantom (FDK reference 0.4163):
lmbda niter/outer before after
1 5/1 0.1420 0.4264
1 10/2 0.1044 0.4197
1 20/4 0.0276 0.3902
100 10/2 0.4509 0.4697
Every configuration is now stable, monotone across outer iterations, and in a
physical value range - and at lmbda=100 it is the best method tested on this
phantom, ahead of FDK 0.4163, FISTA 0.4204, ISTA 0.4379 and ASD-POCS 0.4101.
Note lmbda=1, the default, is far too small for this scaling: ||Ax|| = 5.6e4
against ||Dx|| = 71 for the same image, so the TV block carries ~4% of the
data block's weight. That is a tuning question, not a correctness one, and is
left to the caller now that neither setting explodes.
Tests: the adjoint identity and D(constant)=0 (no GPU), and irn_tv_cgls joins
the end-to-end smoke test it previously could not pass. 22 pass.
AnderBiguri
reviewed
Sep 3, 2026
Comment on lines
+6
to
+21
|
|
||
| ``np.linalg.norm`` sums the squares in the array's OWN dtype. On the arrays | ||
| this toolbox actually handles that is not survivable: the running total | ||
| dwarfs the increments still being added to it and the answer comes back | ||
| LOW. Measured against a chunked float64 reference, on a 3072^2 detector: | ||
|
|
||
| 45 views (4.2e8 elements) 0.4 % low | ||
| 180 views (1.7e9 elements) 3 % low | ||
| 360 views (3.4e9 elements) 16 % low | ||
| 720 views (6.8e9 elements) 26 % low | ||
|
|
||
| A 1024^3 volume (1.1e9 elements) is 1.4 % low. These norms are not | ||
| only reported - they drive step sizes and stopping rules (CGLS's | ||
| ``alpha = gamma / q_norm**2``; the ASD-POCS data-fit and TV-step control), | ||
| so at production sizes the algorithms misbehave rather than merely | ||
| mis-report. |
Member
There was a problem hiding this comment.
Im not sure if all of this is good docs, it seems more like the bug descrition
Contributor
Author
There was a problem hiding this comment.
Fair point. Done in 398d3ea: both docstrings now say only what the function does and why the accumulation is float64 / the return float32. The measurements stay in the PR description.
AnderBiguri
reviewed
Sep 3, 2026
|
|
||
|
|
||
| def _norm(x, ord=2): | ||
| """Euclidean norm of a large float32 array, accumulated in float64. |
Member
There was a problem hiding this comment.
This function seems redundant, as its uniquely used for 2 norm, so it can just call l2norm instead
Contributor
Author
There was a problem hiding this comment.
Agreed. Removed in 398d3ea; all call sites use l2norm directly.
AnderBiguri
reviewed
Sep 3, 2026
Comment on lines
+86
to
+96
| The OUTER loop rebuilds r / p / gamma when orthogonality is lost; the | ||
| INNER loop iterates. This port previously had only the inner loop, so | ||
| the `break` that MATLAB uses to fall back into the outer loop and | ||
| RESTART instead left run_main_iter altogether - CGLS returned after | ||
| however few iterations preceded the first restart, reporting | ||
| "exited due to divergence". Combined with a sentinel bug in the | ||
| give-up test (`re_init_at_iteration + 1 == i` against an initial 0, | ||
| which is true at i == 1 - MATLAB compares `remember == iter` against | ||
| an initial 0 that its 1-based counter can never hit on the first | ||
| check), a single lost step at iteration 1 ended the reconstruction. | ||
| """ |
Member
There was a problem hiding this comment.
Are you sure? This code seems to change it to a more MATLAB way, but restarting existed in the code
Contributor
Author
There was a problem hiding this comment.
The re-init call is there, but it is followed directly by break (line 80 on master), and there is no outer loop for the break to fall into, so the for ends and run_main_iter returns. The rebuilt r/p/gamma are never used, and i = i - 1 has no effect on a range iterator. MATLAB's CGLS.m breaks out of for ii into while iter<niter, which rebuilds and continues; that outer loop is what the port lost. LSQR and LSMR have the same shape.
Concrete run on current master, 32³ head phantom, 30 views, cgls(proj, geo, angles, niter=5, verbose=True):
re-initilization of CGLS called at iteration:3
l2l = [2585.11, 2466.20, 2279.26, 2653.53, 0.0]
The fifth iteration is never executed and the returned volume is the reverted iteration-2 state. Same call on this branch:
re-initilization of CGLS called at iteration:3
l2l = [2585.11, 2466.20, 2279.26, 2255.83, 2178.17]
AnderBiguri
reviewed
Sep 3, 2026
Comment on lines
-72
to
-77
| if self.re_init_at_iteration + 1 == i or not self.restart: | ||
| print("CGLS exited due to divergence.") | ||
| return self.res | ||
| self.re_init_at_iteration=i | ||
| i=i-1 | ||
| self.initialize_algo() |
Contributor
Author
There was a problem hiding this comment.
See the note on line 96 for the break. The change here is the second half of it: re_init_at_iteration starts at 0, so re_init_at_iteration + 1 == i is true at i == 1, and a rise at the first iteration prints "exited due to divergence" regardless of restart. MATLAB's remember == iter starts at 0 against a 1-based counter, so it can never fire on the first check. Hence the -1.
Member
|
OK, this is a lot! for the future, I would prefer these to come independently, as its easier to review. 1- Great! this is indeed something we were trying to tackle mathematically, but this helps a lot! |
trim the l2norm/inner docstrings to what the function does and why (float64 accumulation, float32 return)
Removed unnecessary comments and improved code formatting.
Review follow-up on CERN#774: - drop the _norm wrapper in krylov_subspace_algorithms.py; every call site uses tigre.utilities.im3Dnorm.l2norm directly (all were 2-norms) - CGLS.run_main_iter docstring now describes the two-loop structure and the sentinel, not the history; the stray string literals in LSQR/LSMR become a one-line comment - im3Dnorm: finish the l2norm sentence, fix indentation in inner - tests follow the rename - krylov_subspace_algorithms.py is LF throughout again (upstream is LF)
Contributor
Author
|
Thanks for the review. Taking the points in order; the two style items are already pushed (398d3ea). 1. Glad it helps. 2. Restart. The re-init code was there, but it was followed directly by Separately the give-up test 3. Discarded results. out = algs.hybrid_flsqr_tv(proj, geo, angles, niter=5, lmbda=10)
np.abs(out).max() # 0.0 -- all zeros
out = algs.hybrid_flsqr_tv(proj, geo, angles, niter=5, lmbda=10, init="FDK")
np.array_equal(out, fdk) # True -- the FDK warm start handed back untouchedCalling 4. Thanks. 5. IRN_TV_CGLS, three changes, all inside the class:
On a 64³ phantom (FDK correlation 0.416): lmbda = 1 goes 0.142 → 0.426; lmbda = 100 goes 0.451 → 0.470. On splitting: happy to open these as five PRs (norms, restarts, results, sum→dot, IRN) if that is easier to review; say the word. |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
3 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Four commits, each self-contained and independently checkable; together they
make CGLS/LSQR/LSMR/hybrid_LSQR/AB-GMRES/BA-GMRES/hybrid_fLSQR_TV/IRN_TV_CGLS
behave at production sizes. No CUDA or
.pyxchanges — pure Python, norebuild. No new dependencies.
1. A float32
np.linalg.normcannot measure a production-sized residualnp.linalg.normaccumulates the sum of squares in the array's OWN dtype. On aprojection-sized float32 array the running sum saturates against the
increments still being added, and the answer comes back low. Measured against
a chunked float64 reference on N(0, 1e-3) data at a 3072×3072 detector:
np.linalg.normerrorA 1024³ volume (1.1e9 elements) is already ~1.4 % off. This is not a
reporting detail: every Krylov step size is a RATIO of such norms
(
alpha = gamma / q_norm**2), and the loss-of-orthogonality test comparesconsecutive residual norms — so at production sizes CGLS both steps wrongly
(volume-norm ÷ projection-norm: the two biases do not cancel) and mistakes
accumulator noise for divergence. On a real 720-view / 1024³ dataset CGLS
"exited due to divergence" at iteration 1; with this change it runs to
completion and its mean lands in family with the other iterative methods.
Two things decide severity, both worth knowing: a sum of squares is far
worse than a dot of two different vectors (nothing cancels), and the bias is
BLAS-dependent (the same sizes under a different BLAS were ~76× more
accurate), so it can look absent on one machine and be 26 % on another.
Fix:
tigre.utilities.im3Dnormgainsl2norm()andinner()that accumulatein float64 and return float32, so no caller's dtype changes (Ax/Atb reject
float64). Routed through them: every Krylov norm,
im3DNORM(...,2)(so theASD-POCS family gets it too),
svd_power_method(this IS FISTA's Lipschitzconstant) and the POCS TV-step cosine.
2. The CGLS/LSQR/LSMR restart never restarted
MATLAB's
CGLS.mbreaks out of an INNER loop into an OUTER one that rebuildsr/p/gamma. The Python port kept only the inner loop, so that
breakleftrun_main_iteraltogether — the method returned after however few iterationspreceded the first restart. Its give-up sentinel also read
re_init_at_iteration + 1 == iagainst an initial 0, which is true ati == 1. One lost step at iteration 1 ended the reconstruction. LSQR and LSMRhad the identical structure. All three now have the two-loop form and a
-1sentinel.
hybrid_LSQR, hybrid_fLSQR_TV and AB/BA-GMRES have no restart path at all — their
sentinel could only fire at one iteration and ignored a rise at every other —
so they now simply stop on a genuine rise.
3. Three algorithms discarded their result
decorator()ignores whatrun_main_iterRETURNS and hands backself.res.Any
return <expr>that does not first assignself.restherefore returns theINITIAL volume — zeros, or with an FDK warm start an unimproved FDK image that
looks entirely plausible.
ba_gmresandhybrid_lsqrreturned all-zerovolumes on every early exit;
hybrid_fLSQR_TVnever assignedself.resonany path, so its output was discarded always.
4. AB/BA-GMRES orthogonalised with Python's builtin
sum()h[k,i] = sum(qk.ravel()*w[i])iterates a CT-sized array element by element:~40× slower than
np.dot, 800× less accurate. BA_GMRES's basis isvolume-space → ~107 s per coefficient × 210 coefficients at niter 20; on a
real 1024³ problem that was most of a 4.1 h run (41.6 min after this change).
5. IRN_TV_CGLS: an operator that was not an adjoint, and no divergence guard
CGLS runs on the stacked operator
[A; √λ·W·D]and assumes the handedtranspose IS the adjoint.
Lx/Ltxhad an inner-product asymmetry of6.4e-2:
Dxx = np.copy(img)then writing only[0:-2]left the last twoslices holding raw image VALUES (so
D(constant) ≠ 0— a Tikhonov termsmuggled into part of the volume), and
Dᵀwas off by one at n-2. Plus noinner divergence guard (plain CGLS has one): the residual rose from step 2 and
compounded 9.7e3 → 6.9e5. The guard compares only WITHIN one outer iteration —
across a reweighting the residuals belong to different weighted problems. On a
64³ phantom (FDK correlation 0.4163): λ=1 → 0.142 before, 0.4264 after;
λ=100 → 0.451 before, 0.4697 after.
Tests
Python/tests/test_krylov_restart_and_norms.py(no GPU): float64 accumulationvs a chunked reference, float32 return type, restart re-entry and sentinel,
the D/Dᵀ adjoint, the ordering of
sum()vsnp.dot.Python/tests/test_krylov_smoke.py(GPU): every Krylov method runs, returns afinite non-zero volume and improves on its warm start.
Verified on Windows (RTX A6000, CUDA 13.3) against current master. Merge
conflicts with #770 are expected to be trivial (#770 changes 3 lines of
IRN_TV_CGLS scalar dtypes; this branch keeps upstream's
np.sqrt(..., dtype=np.float32)form for those lines).