RotationScaleMerge: Anderson-accelerate the per-frame scaling loops (RESULT-CHANGING)

The partials and fulls scaling loops are alternating least squares that converge linearly:
over the rc176 battery's 2527 loops a round shrinks the step by a median 0.86 (partials) and
0.93 (fulls), up to 0.98, so the 1e-3 step tolerance stopped them 0.6% / 1.4% (median) and up to
17% (fulls p90) short of their fixed point, and 313 of them were stopped as "stalled" when they
were only slow - a stalled fulls loop also skipped ShrinkToRestrained.

Each round is now accelerated in log G by type-II Anderson mixing (Walker & Ni 2011), depth 5,
each frame weighted by its observation count as in the step test; frames below the credible floor
take the plain round's scale. The history is dropped when the weighted residual grows or the set
of fitted frames changes (new AndersonMixing class, unit-tested on a synthetic slow ALS loop). The
mixed scales are gauge-pinned again before the rescale, and the stop test is unchanged: it reads
the round's residual F(x) - x, so "stalled" now means a loop that does not settle under
acceleration either.

Measured (CPU, -N 8): 4nwv 311 -> 82 iterations over 9 loops, 9v16 76 iterations over 11 loops,
the same counts as the env-gated prototype on tail-solver; no history restart fired on either.
Depth 3/5/8 give 92/91/90 iterations on 4nwv. A restart on a change of the fulls loop's
cross-validated lambda was tried and gave nothing (83 vs 82). Same space group and merging
statistics on 4nwv and 9v16.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PXc1zD1yxHqaR1DGdCTHL4

docs/ACKNOWLEDGEMENT.md

@@ -523,6 +523,13 @@ corrected once by the exact remainder, which returns the correctly rounded quoti
 P. W. Markstein, "Computation of elementary functions on the IBM RISC System/6000 processor"
 (1990), IBM J. Res. Dev. 34, 111-119.
 
+The per-frame scaling loops of the rotation merge are accelerated by Anderson acceleration in the
+type-II form analysed by Walker and Ni: H. F. Walker and P. Ni, "Anderson acceleration for fixed-point
+iterations" (2011), SIAM J. Numer. Anal. 49, 1715-1735
+[doi:10.1137/10078356X](https://doi.org/10.1137/10078356X), after D. G. Anderson, "Iterative
+procedures for nonlinear integral equations" (1965), J. ACM 12, 547-560
+[doi:10.1145/321296.321305](https://doi.org/10.1145/321296.321305).
+
 The single-precision exp the Bragg integration's profile fit shares between the CPU and an Apple GPU
 (`BraggExp.h`) is built from textbook error-free floating-point arithmetic: the argument reduction of
 W. J. Cody and W. Waite, "Software Manual for the Elementary Functions" (1980), Prentice-Hall; the

docs/CPU_DATA_ANALYSIS.md

@@ -57,6 +57,7 @@ The methods draw on, and in places reimplement, solutions from:
 - S. French & K. Wilson, "On the treatment of negative intensity observations", *Acta Cryst.* **A34** (1978), 517-525 (Bayesian amplitude estimation from intensities), and CCP4's ctruncate (C. Ballard & N. Stein), whose anisotropic Wilson prior §10.8 follows, cited through M. D. Winn et al., "Overview of the CCP4 suite and current developments", *Acta Cryst.* **D67** (2011), 235-242.
 - A. T. Brünger, "Free R value: a novel statistical quantity for assessing the accuracy of crystal structures", *Nature* **355** (1992), 472-475 (R-free cross-validation).
 - P. H. C. Eilers, "A perfect smoother", *Anal. Chem.* **75** (2003), 3631-3636, after E. T. Whittaker, "On a new method of graduation", *Proc. Edinburgh Math. Soc.* **41** (1923), 63-75 (the penalised smoother of the fulls' per-frame scale, §10.6).
+- H. F. Walker & P. Ni, "Anderson acceleration for fixed-point iterations", *SIAM J. Numer. Anal.* **49** (2011), 1715-1735, after D. G. Anderson, "Iterative procedures for nonlinear integral equations", *J. ACM* **12** (1965), 547-560 (the acceleration of the per-frame scaling loops, §10.6).
 - R. A. Fisher, "Frequency distribution of the values of the correlation coefficient in samples from an indefinitely large population", *Biometrika* **10** (1915), 507-521 (the z-transformation on which the correction surfaces' held-out half-set CC1/2 is compared).
 - M. Wojdyr, "GEMMI: A library for structural biology", *J. Open Source Softw.* **7** (2022), 4200 (model / structure-factor / map machinery used in §14).
 - J. P. Wright, "Experiences with GPU decompression for bitshuffle + LZ4 data", HDF5 User Group meeting (2021), and [github.com/jonwright/bslz4decoders](https://github.com/jonwright/bslz4decoders) (device-side decoding of bitshuffle+LZ4 images, §0).

docs/CPU_DATA_ANALYSIS_INTEGRATION.md

@@ -388,6 +388,8 @@ Before the combine a per-frame scale can also be fitted on the partials themselv
 
 How smooth that scale is over the rotation is left to the data rather than to a fixed window. Each round fits every frame on its own fulls against the current reference, giving a scale $G_f$ and its information $D_f=\sum w^2c^2$; the scale is then the curve $x=\log G$ minimising $\sum_f J_f\,(x_f-y_f)^2+\lambda\sum_f(\Delta^2 x)_f^2$, with $y_f$ the frame's own fit in log scale and $J_f$ its information carried there — a penalised (Whittaker–Eilers) smoother, solved as a five-band linear system. $\lambda$ is chosen by cross-validation: blocks of frames one rocking curve wide are left out in turn and predicted from the curve through the rest (neighbours closer than a rocking curve share their measurement, since a full sums those frames). A frame of hundreds of fulls is then followed frame by frame, a frame of two or three is carried by its neighbours, a stretch with none is bridged by a straight line, and a scale that falls a hundredfold over a few degrees — an absorbing crystal turning edge-on — is followed where a window would average across it. Once the curve settles, one free fit is shrunk toward it frame by frame by how much of each frame's deviation its neighbour shares (the lag-1 covariance), which hands back a real per-frame systematic and discards fit noise.
 
+Both per-frame scaling loops — on the partials and on the fulls — alternate between the references given the scales and the scales given the references, and converge linearly: a round typically removes only 5–15% of the remaining step, and 2% on the slowest sweeps. Each round is therefore Anderson-accelerated in $\log G$ (type II, five rounds of history, each frame weighted by its observation count): the plain round's scales are corrected by the combination of the last rounds' steps that best cancels the current one, and the history is dropped whenever the step grows or the set of fitted frames changes. This changes the path, not the fixed point; the loops stop at the same tolerance (an rms step of $10^{-3}$ in $\log G$) in about a quarter of the rounds and closer to the converged scales than the plain iteration stopped there.
+
 After scale-fulls, five **correction surfaces** are fitted on the combined fulls (rotation path, **on by default**; disable all with `--no-scaling-corrections`), each an alternating multiplicative refinement of the per-full scale against the merged reference:
 
 - **Decay.** Radiation damage weakens later frames more at higher resolution — a resolution×time (Debye–Waller) systematic the resolution-flat per-image scale cannot capture. A single global relative-$B$ rate is fitted, $\ln(I_\mathrm{ref}/I_\mathrm{obs}) = 2\,(\mathrm{d}B/\mathrm{d}n)\,(n-\bar n)\,s^2$ (frame $n$, $s^2 = 1/4d^2$), and folded into the scale. It engages only when the total relative-$B$ over the run exceeds a physical floor (2 Ų); below that the decay is negligible and "correcting" it only spreads symmetry equivalents (same $s^2$, different frames). An optional **per-batch relative-$B$** (`--relative-b[=deg]`, off unless requested; 10°-of-rotation batches by default) extends the single global rate to a smooth $B(n)$ curve — the same $s^2$-weighted decay fit solved independently over short frame batches, curvature-penalized so it cannot over-fit and cross-validated like the surfaces below — for crystals whose decay is non-linear in dose. Its cross-validation splits on **ASU-group parity**, not the frame parity the surfaces below use: a per-batch parameter owns whole frames and so cannot be scored on a held-out frame, whereas splitting the symmetry equivalents tests whether a batch's $B$ generalises to reflections it was not fitted on.

image_analysis/scale_merge/AndersonMixing.cpp

@@ -0,0 +1,47 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "AndersonMixing.h"
+
+std::vector<double> AndersonMixing::Next(const std::vector<double> &x, const std::vector<double> &fx,
+                                         const std::vector<double> &w) {
+    // Anderson acceleration, type II, after Walker & Ni (2011) SIAM J. Numer. Anal. 49, 1715-1735
+    // (Anderson (1965) J. ACM 12, 547-560): with dX and dR the differences of the states and of the
+    // residuals of consecutive rounds, gamma = argmin |r - dR gamma|_w and the next state is
+    // F(x) - (dX + dR) gamma.
+    const int n = static_cast<int>(x.size());
+    const Eigen::Map<const Eigen::VectorXd> xv(x.data(), n), fxv(fx.data(), n), wv(w.data(), n);
+    const Eigen::VectorXd r = fxv - xv;
+    std::vector<uint8_t> support(n);
+    for (int i = 0; i < n; ++i) support[i] = w[i] > 0.0 ? 1 : 0;
+    const double r_norm = (wv.array() * r.array().square()).sum();
+
+    if (x_last.size() == n && support == support_last && r_norm <= r_norm_last) {
+        dx.push_back(xv - x_last);
+        dr.push_back(r - r_last);
+        if (static_cast<int>(dx.size()) > depth) {
+            dx.pop_front();
+            dr.pop_front();
+        }
+    } else if (!dx.empty()) {
+        dx.clear();
+        dr.clear();
+        ++restarts;
+    }
+    x_last = xv;
+    r_last = r;
+    support_last = std::move(support);
+    r_norm_last = r_norm;
+
+    std::vector<double> out = fx;
+    const int m = static_cast<int>(dr.size());
+    if (m == 0) return out;
+    const Eigen::VectorXd sqrt_w = wv.cwiseMax(0.0).cwiseSqrt();
+    Eigen::MatrixXd M(n, m);
+    for (int j = 0; j < m; ++j) M.col(j) = sqrt_w.cwiseProduct(dr[j]);
+    const Eigen::VectorXd gamma = M.colPivHouseholderQr().solve(sqrt_w.cwiseProduct(r));
+    for (int j = 0; j < m; ++j)
+        for (int i = 0; i < n; ++i)
+            if (w[i] > 0.0) out[i] -= gamma[j] * (dx[j][i] + dr[j][i]);
+    return out;
+}

image_analysis/scale_merge/AndersonMixing.h

@@ -0,0 +1,40 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <cstdint>
+#include <deque>
+#include <vector>
+
+#include <Eigen/Dense>
+
+// Acceleration of a fixed-point iteration x -> F(x) that converges linearly. Each round hands in the
+// state x it started from and F(x), where one plain round took it; the residual r = F(x) - x is what
+// the plain iteration would step. The state the next round starts from is F(x) corrected by the
+// combination of the last few rounds' steps that best cancels r, so the slow directions the plain
+// iteration creeps along - a handful of them, whose contraction is close to one - are taken in a few
+// rounds. Where r = 0 the correction is 0: the fixed point is the plain iteration's.
+//
+// Each component carries a weight in the norm r is measured in; a component of weight 0 takes F(x)
+// unchanged and is not part of the fit. The history is dropped, and the round is a plain one, when
+// the residual grew - the history no longer describes the map from here - or when the set of
+// components of non-zero weight changed, since the differences of the history are then differences
+// of different things.
+class AndersonMixing {
+public:
+    explicit AndersonMixing(int depth) : depth(depth) {}
+
+    std::vector<double> Next(const std::vector<double> &x, const std::vector<double> &fx,
+                             const std::vector<double> &w);
+
+    int Restarts() const { return restarts; }
+
+private:
+    int depth;
+    int restarts = 0;
+    std::deque<Eigen::VectorXd> dx, dr;   // x and r differences of consecutive rounds, newest last
+    Eigen::VectorXd x_last, r_last;
+    std::vector<uint8_t> support_last;
+    double r_norm_last = 0.0;
+};

image_analysis/scale_merge/CMakeLists.txt

@@ -37,6 +37,8 @@ ADD_LIBRARY(JFJochScaleMerge
         WilsonOutliers.h
         ErrorModel.cpp
         ErrorModel.h
+        AndersonMixing.cpp
+        AndersonMixing.h
         SublatticeRelabel.cpp
         SublatticeRelabel.h)
 TARGET_LINK_LIBRARIES(JFJochScaleMerge Ceres::ceres Eigen3::Eigen JFJochCommon fftw3f)

image_analysis/scale_merge/RotationScaleMerge.cpp

@@ -26,6 +26,7 @@
 #include <gemmi/unitcell.hpp>
 #include <gemmi/symmetry.hpp>
 
+#include "AndersonMixing.h"
 #include "HKLKey.h"
 #include "RfreeFlags.h"
 #include "FrenchWilson.h"
@@ -60,6 +61,10 @@ namespace {
     // The per-frame scaling loop has settled when the scales moved less than this between two
     // iterations, as the rms of |log(G_new/G_old)| over the frames fitted in both.
     constexpr double SCALING_TOLERANCE = 1e-3;
+    // The rounds of history the per-frame scaling loops are accelerated with (AndersonMixing). Their slow
+    // directions are a few smooth modes over the sweep: measured on a 4nwv-class sweep, 3, 5 and 8 rounds
+    // take its ten loops in 92, 91 and 90 iterations, against 364 without.
+    constexpr int SCALING_ANDERSON_DEPTH = 5;
     // The per-frame scale is smoothed over at least this many rocking curves (see Run).
     constexpr int SMOOTH_G_MIN_ROCKING_EVENTS = 6;
 
@@ -1846,12 +1851,22 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunScalingLoop(
     // and handing any of it back re-opens it (measured: ISa 23.8 -> 22.1 on a 360 deg, 14x merge). The
     // fulls, each a whole rocking event, have their own loop (RunFullsScalingLoop).
     //
+    // The iteration converges linearly, and slowly: a round takes 5-15% off the step on most sweeps and
+    // 2% on the slowest, so a loop stopped on its step stopped several steps' worth short of where it
+    // was going. Each round's scales are therefore accelerated (AndersonMixing, in log G, each frame at
+    // its observation count): the plain round's answer corrected by the combination of the last few
+    // rounds' steps that best cancels this one. That moves no fixed point - only the path to it - and
+    // it takes the loop to its answer in a quarter of the rounds and closer: on one sweep the merged
+    // intensities ended 0.02% rms from the converged point, against 0.12% for the plain loop stopped at
+    // the same tolerance. A frame below the credible floor (below) takes the plain round's scale.
+    //
     // The loop stops when the relative scales stop moving - the rms of |log(G_new/G_old)| over the
     // frames, each weighted by its observation count and nothing else (a weight that vanishes with G
     // is what let a collapsing sweep report itself settled) - or when they stop settling: a step that
     // has not fallen below nine tenths of its smallest value for five rounds is a loop on a flat or
     // rising valley, and more rounds of it can only walk further. Both leave the frames where they are
-    // and say so.
+    // and say so. The step is the residual of the round, F(x) - x at the state the round started from,
+    // so the test reads the same whether that state was a plain round's or an accelerated one.
     ScalingLoopOutcome out;
     std::vector<uint8_t> scaled_any(n_frames, 0);
     std::vector<double> g_prev(n_frames, 0.0), g_smooth;
@@ -1859,6 +1874,8 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunScalingLoop(
     std::vector<double> ratio(n_frames, 1.0);
     double best_step = INFINITY;
     int rounds_since_best = 0;
+    AndersonMixing anderson(SCALING_ANDERSON_DEPTH);
+    std::vector<double> log_g_prev(n_frames), log_g(n_frames), weight(n_frames);
     for (int it = 0; it < scaling_iter; ++it) {
         iterate();
         out.iterations = it + 1;
@@ -1874,15 +1891,15 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunScalingLoop(
                 if (scaled_any[f] && std::isfinite(g_smooth[f]) && g_smooth[f] > 0.0)
                     g[f] = g_smooth[f];
         }
-        for (int f = 0; f < n_frames; ++f)
-            n_obs[f] = scaled_any[f] && std::isfinite(g[f]) && g[f] > 0.0 ? frame_obs_count[f] : 0;
-        const double g_ref = TypicalFrameScale(g, n_obs);
+        std::vector<uint8_t> usable(n_frames);
         for (int f = 0; f < n_frames; ++f) {
-            if (!scaled_any[f] || !(std::isfinite(g[f]) && g[f] > 0.0)) { ratio[f] = 1.0; continue; }
-            g[f] /= g_ref;
-            ratio[f] /= g[f];
+            usable[f] = scaled_any[f] && std::isfinite(g[f]) && g[f] > 0.0 ? 1 : 0;
+            n_obs[f] = usable[f] ? frame_obs_count[f] : 0;
         }
-        rescale(scaled_any, ratio);
+        const double g_ref = TypicalFrameScale(g, n_obs);
+        for (int f = 0; f < n_frames; ++f)
+            if (usable[f]) g[f] /= g_ref;
+        bool stop = false;
         if (it > 0) {
             // A frame already below the credible floor is not part of the fit's answer: the window
             // leaves it out (ComputeSmoothGWindow) and the guard after the loop drops it. Left in the
@@ -1899,28 +1916,47 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunScalingLoop(
             }
             double ss = 0.0, sw = 0.0;
             int n = 0;
-            for (int f = 0; f < n_frames; ++f)
+            for (int f = 0; f < n_frames; ++f) {
+                weight[f] = 0.0;
                 if (scaled_any[f] && g_prev[f] > 0.0 && g[f] > 0.0 && g[f] >= g_floor) {
-                    const double d = std::log(g[f] / g_prev[f]);
+                    log_g_prev[f] = std::log(g_prev[f]);
+                    log_g[f] = std::log(g[f]);
+                    const double d = log_g[f] - log_g_prev[f];
                     ss += frame_obs_count[f] * d * d;
                     sw += frame_obs_count[f];
+                    weight[f] = frame_obs_count[f];
                     ++n;
+                } else {
+                    log_g_prev[f] = log_g[f] = 0.0;
                 }
+            }
             out.step = sw > 0.0 ? std::sqrt(ss / sw) : 0.0;
             logger->Debug("Per-frame scaling ({}): iteration {}, rms |dlogG| {:.2e} over {} frames, "
-                         "gauge {:.4f}", what, it + 1, out.step, n, g_ref);
+                         "gauge {:.4f}, {} restarts", what, it + 1, out.step, n, g_ref, anderson.Restarts());
             if (out.step < SCALING_TOLERANCE) {
                 out.converged = true;
-                break;
-            }
-            if (out.step < 0.9 * best_step) {
+                stop = true;
+            } else if (out.step < 0.9 * best_step) {
                 best_step = out.step;
                 rounds_since_best = 0;
             } else if (++rounds_since_best >= 5) {
                 out.stalled = true;
-                break;
+                stop = true;
+            }
+            if (!stop) {
+                // The next round starts from the accelerated scales, gauge-pinned again.
+                const std::vector<double> next = anderson.Next(log_g_prev, log_g, weight);
+                for (int f = 0; f < n_frames; ++f)
+                    if (weight[f] > 0.0) g[f] *= std::exp(next[f] - log_g[f]);
+                const double g_ref_next = TypicalFrameScale(g, n_obs);
+                for (int f = 0; f < n_frames; ++f)
+                    if (usable[f]) g[f] /= g_ref_next;
             }
         }
+        for (int f = 0; f < n_frames; ++f)
+            ratio[f] = usable[f] ? ratio[f] / g[f] : 1.0;
+        rescale(scaled_any, ratio);
+        if (stop) break;
         g_prev = g;
     }
     frame_scaled_scratch = scaled_any;
@@ -1966,7 +2002,9 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunFullsScalingLoop(
     // zero, where DropCollapsedFullScales finds it, without one round's noise throwing it there.
     //
     // A frame below the credible floor is left out of the step, as in RunScalingLoop: it walks
-    // toward zero for as long as the loop runs, which is the answer, not a failure to settle.
+    // toward zero for as long as the loop runs, which is the answer, not a failure to settle. The
+    // rounds are accelerated as RunScalingLoop's are, in log G with each frame at its count of fulls;
+    // the frame below the floor takes the plain round's scale.
     ScalingLoopOutcome out;
     std::vector<double> y(n_frames), J(n_frames), g_new(n_frames), ratio(n_frames);
     std::vector<int64_t> n_obs(n_frames);
@@ -1977,6 +2015,8 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunFullsScalingLoop(
     }
     double best_step = INFINITY, lambda = 0.0;
     int rounds_since_best = 0;
+    AndersonMixing anderson(SCALING_ANDERSON_DEPTH);
+    std::vector<double> log_g(n_frames), log_g_new(n_frames), weight(n_frames);
     for (int it = 0; it < scaling_iter; ++it) {
         fit();
         out.iterations = it + 1;
@@ -1995,33 +2035,50 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunFullsScalingLoop(
         std::vector<double> scales;
         for (int f = 0; f < n_frames; ++f) {
             g_new[f] /= g_ref;
-            ratio[f] = g[f] / g_new[f];
             if (has_fulls[f]) scales.push_back(g_new[f]);
         }
-        rescale(has_fulls, ratio);
         const double g_floor = scales.empty() ? 0.0 : median_of(scales) * MIN_CREDIBLE_SCALE_RATIO;
         double ss = 0.0, sw = 0.0;
-        for (int f = 0; f < n_frames; ++f)
+        for (int f = 0; f < n_frames; ++f) {
+            weight[f] = 0.0;
             if (has_fulls[f] && g_new[f] >= g_floor) {
-                const double d = std::log(g_new[f] / g[f]);
+                log_g[f] = std::log(g[f]);
+                log_g_new[f] = std::log(g_new[f]);
+                const double d = log_g_new[f] - log_g[f];
                 ss += frame_obs_count[f] * d * d;
                 sw += frame_obs_count[f];
+                weight[f] = frame_obs_count[f];
+            } else {
+                log_g[f] = log_g_new[f] = 0.0;
             }
-        g = g_new;
+        }
         out.step = sw > 0.0 ? std::sqrt(ss / sw) : 0.0;
-        logger->Debug("Per-frame scaling (fulls): iteration {}, rms |dlogG| {:.2e}, lambda {:.3g}, gauge {:.4f}",
-                     it + 1, out.step, lambda, g_ref);
+        logger->Debug("Per-frame scaling (fulls): iteration {}, rms |dlogG| {:.2e}, lambda {:.3g}, gauge {:.4f}, "
+                     "{} restarts", it + 1, out.step, lambda, g_ref, anderson.Restarts());
+        bool stop = false;
         if (out.step < SCALING_TOLERANCE) {
             out.converged = true;
-            break;
-        }
-        if (out.step < 0.9 * best_step) {
+            stop = true;
+        } else if (out.step < 0.9 * best_step) {
             best_step = out.step;
             rounds_since_best = 0;
         } else if (++rounds_since_best >= 5) {
             out.stalled = true;
-            break;
+            stop = true;
         }
+        if (!stop) {
+            // The next round starts from the accelerated scales, gauge-pinned again.
+            const std::vector<double> next = anderson.Next(log_g, log_g_new, weight);
+            for (int f = 0; f < n_frames; ++f)
+                if (weight[f] > 0.0) g_new[f] *= std::exp(next[f] - log_g_new[f]);
+            const double g_ref_next = TypicalFrameScale(g_new, n_obs);
+            for (int f = 0; f < n_frames; ++f) g_new[f] /= g_ref_next;
+        }
+        for (int f = 0; f < n_frames; ++f)
+            ratio[f] = g[f] / g_new[f];
+        rescale(has_fulls, ratio);
+        g = g_new;
+        if (stop) break;
     }
     // The curve is the part of the scale the frames share over a rocking curve and more; what a
     // frame's own fulls say beyond it is handed back by how much of it is real (ShrinkToRestrained) -

tests/AndersonMixingTest.cpp

@@ -0,0 +1,115 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <cmath>
+#include <catch2/catch_all.hpp>
+
+#include "../image_analysis/scale_merge/AndersonMixing.h"
+
+namespace {
+    // A noiseless sweep of 24 frames whose reflections each span three consecutive frames, scaled by the
+    // alternating least squares of the per-frame scaling loop: the reference intensities from the
+    // scales, the scales from the references, the gauge pinned at a geometric mean of 1. A frame is
+    // tied only to its neighbours, so a smooth error of the scales over the sweep decays by a few
+    // percent a round - the slow linear convergence of the real loop.
+    struct SlowScalingLoop {
+        static constexpr int n_frames = 24, span = 3;
+        std::vector<int> first;   // reflection h is on frames first[h] .. first[h] + span - 1
+        std::vector<double> I_true, G_true;
+
+        SlowScalingLoop() {
+            for (int h = 0; h + span <= n_frames; ++h)
+                for (int rep = 0; rep < 3; ++rep) {
+                    first.push_back(h);
+                    I_true.push_back(1.0 + 0.5 * std::sin(1.7 * h + rep));
+                }
+            for (int f = 0; f < n_frames; ++f)
+                G_true.push_back(1.0 + 0.4 * std::sin(6.283 * f / n_frames) + 0.2 * std::cos(2.1 * f));
+        }
+
+        std::vector<double> Round(const std::vector<double> &g) const {
+            std::vector<double> num(n_frames, 0.0), den(n_frames, 0.0), out(n_frames);
+            for (size_t h = 0; h < first.size(); ++h) {
+                double mean = 0.0;
+                for (int k = 0; k < span; ++k) mean += G_true[first[h] + k] * I_true[h] / g[first[h] + k];
+                mean /= span;
+                for (int k = 0; k < span; ++k) {
+                    num[first[h] + k] += G_true[first[h] + k] * I_true[h] * mean;
+                    den[first[h] + k] += mean * mean;
+                }
+            }
+            double s = 0.0;
+            for (int f = 0; f < n_frames; ++f) { out[f] = num[f] / den[f]; s += std::log(out[f]); }
+            for (int f = 0; f < n_frames; ++f) out[f] /= std::exp(s / n_frames);
+            return out;
+        }
+
+        // Rounds to a step below `tolerance`, and the largest |log G - log G_true| there (gauge-pinned).
+        std::pair<int, double> Run(int depth, double tolerance) const {
+            std::vector<double> g(n_frames, 1.0), log_g(n_frames), log_g_new(n_frames), w(n_frames, 1.0);
+            AndersonMixing anderson(depth);
+            int rounds = 0;
+            for (double step = INFINITY; step >= tolerance && rounds < 100000; ++rounds) {
+                std::vector<double> g_new = Round(g);
+                double ss = 0.0;
+                for (int f = 0; f < n_frames; ++f) {
+                    log_g[f] = std::log(g[f]);
+                    log_g_new[f] = std::log(g_new[f]);
+                    ss += (log_g_new[f] - log_g[f]) * (log_g_new[f] - log_g[f]);
+                }
+                step = std::sqrt(ss / n_frames);
+                if (depth > 0) {
+                    const std::vector<double> next = anderson.Next(log_g, log_g_new, w);
+                    double s = 0.0;
+                    for (double v : next) s += v;
+                    for (int f = 0; f < n_frames; ++f) g_new[f] = std::exp(next[f] - s / n_frames);
+                }
+                g = g_new;
+            }
+            double s = 0.0, err = 0.0;
+            for (double v : G_true) s += std::log(v);
+            for (int f = 0; f < n_frames; ++f)
+                err = std::max(err, std::fabs(std::log(g[f]) - (std::log(G_true[f]) - s / n_frames)));
+            return {rounds, err};
+        }
+    };
+}
+
+TEST_CASE("AndersonMixing_SlowScalingLoop", "[RotationScale]") {
+    const SlowScalingLoop loop;
+    const auto [plain_rounds, plain_err] = loop.Run(0, 1e-6);
+    const auto [aa_rounds, aa_err] = loop.Run(5, 1e-6);
+    // The plain loop creeps: hundreds of rounds, and stopped on its step it still sits a hundred
+    // times the tolerance from the answer.
+    CHECK(plain_rounds > 300);
+    CHECK(plain_err > 5e-5);
+    // Accelerated: the same fixed point, in a tenth of the rounds, and closer to it at the same stop.
+    CHECK(aa_rounds * 10 < plain_rounds);
+    CHECK(aa_err < 0.5 * plain_err);
+}
+
+TEST_CASE("AndersonMixing_ZeroWeightTakesThePlainStep", "[RotationScale]") {
+    // x -> 0.9 x on two components, the second not part of the fit.
+    AndersonMixing anderson(5);
+    std::vector<double> x{1.0, 1.0};
+    const std::vector<double> w{1.0, 0.0};
+    for (int it = 0; it < 4; ++it) {
+        const std::vector<double> fx{0.9 * x[0], 0.9 * x[1]};
+        const std::vector<double> next = anderson.Next(x, fx, w);
+        CHECK(next[1] == fx[1]);
+        x = next;
+    }
+    // A linear map on the fitted component is solved by the first mixed step.
+    CHECK(std::fabs(x[0]) < 1e-12);
+    CHECK(anderson.Restarts() == 0);
+}
+
+TEST_CASE("AndersonMixing_RestartsWhenTheResidualGrows", "[RotationScale]") {
+    AndersonMixing anderson(5);
+    const std::vector<double> w{1.0};
+    CHECK(anderson.Next({1.0}, {0.5}, w)[0] == 0.5);    // no history: the plain step
+    anderson.Next({0.5}, {0.25}, w);                      // residual fell: history kept
+    CHECK(anderson.Restarts() == 0);
+    CHECK(anderson.Next({0.0}, {2.0}, w)[0] == 2.0);    // residual grew: history dropped, plain step
+    CHECK(anderson.Restarts() == 1);
+}

tests/CMakeLists.txt

@@ -109,6 +109,7 @@ ADD_EXECUTABLE(jfjoch_test
         DiagnosticOutputTest.cpp
         PostRefineTest.cpp
         RotationScaleWalkTest.cpp
+        AndersonMixingTest.cpp
         RugnuxLargeTest.cpp
         TestData.h
         MovingAverageTest.cpp