8. Estimation¶
8.1. Objective and weights¶
Source: rietx.optimize.least_squares
minimised over the residual rows of (1.5) [Rietveld, 1969]. Weights come from the data file’s esd column whenever it is present; Poisson \(\sigma = \sqrt{\max(y, 1)}\) is only the fallback for bare counts.
The Jacobian is assembled column by column, taking exact work wherever it exists: linear background columns, analytic peak-chain columns (everything flowing through per-peak position, width, mixing and intensity), analytic site-DOF columns for coordinates and anisotropic ADPs over the frozen operator subsets, analytic axial columns for the FCJ apertures, and plain forward differences as the fallback. Under a differentiable backend the same residual is traced and columns come from forward-mode autodiff, whose cost scales with the parameter count [Nocedal and Wright, 2006]; every backend is held to per-column agreement with the analytic Jacobian.
8.2. Agreement statistics¶
Defined per Toby [Toby, 2006]:
Source: rietx.optimize.statistics
with \(\mathrm{GoF} = \sqrt{\chi^2_{\mathrm{red}}} = R_{wp}/R_{\exp}\), plus the background-subtracted \(R_{wp}\) variant Toby recommends when the background carries much of the raw intensity. The Durbin-Watson statistic on weighted residuals [Hill and Flack, 1987] flags serial correlation (\(d \approx 2\) ⇒ uncorrelated).
Across the shipped corrections, \(\Delta R_{wp}\) recurs as a poor judge of a physical improvement (Intensity corrections, Reading a paper against its own numbers). These indices measure agreement, and correctness is a separate question.
8.3. Structure agreement indices¶
Those indices compare profiles. The two that compare the structure are unweighted sums over reflections [McCusker et al., 1999] (their eqs 13 and 14):
Source: rietx.optimize.statistics.structure_r_factors
with \(m\) the reflection multiplicity, so \(|F| = \sqrt{I/m}\) and \(R_F\) is the index a single-crystal \(R\) is comparable with. Powder data measure neither sum directly. \(I_o\) is the observed profile partitioned in proportion to \(I_c\), the same partition Intensities performs for Le Bail extraction, evaluated once on the converged structural model. Two consequences follow from the definition itself.
Both indices are biased towards the model being tested, because a wrong model receives the intensity it predicted. The paper introduces them for monitoring a structure’s improvement, and not for judging one in isolation. Both are also unweighted, so a reflection the weighted fit barely constrains counts as much as one that dominates it. No weighted variant is computed here, and a trace phase’s \(R_B\) does not compare with the major phase’s.
8.4. The effective observation count¶
\(N\) in (8.2) counts profile steps, and the steps across one peak are repeated measurements of one number. Only the integrated intensities of individual reflections are unique observations [McCusker et al., 1999], and overlap reduces even those: two reflections at one \(2\theta\) are one observation, and two that partly overlap lie between one and two, because the profile shape still says how to split them. Altomare et al. [Altomare et al., 1995] make that a count,
Source: rietx.optimize.statistics.effective_observations
where \(G\) is the symmetric profile of Peak profiles and \(\chi_k\) is the part of reflection \(k\)’s own interval (\(\pm\)2.0 FWHM) on which some overlapping reflection stands higher. An isolated line contributes 1 and the weaker of an exactly coincident pair contributes 0, so the pair is one observation. The guideline built on it asks for at least 3.0 and preferably 5.0 effective observations per structural parameter. Those are the atomic ones, since peak positions rather than intensities pay for the cell, profile and background terms.
The estimate is not a theorem, and both papers say so: the approach “may not have a rigorous basis”. Its own \(\alpha\) is a case in point. The paper’s check at \(\alpha = 4\) lands 6.5 % lower on average, so the value tabulated at \(\alpha = 2\), reported here for comparability with it, runs a little generous.
8.5. Esds and the Bérar-Lelann inflation¶
Source: rietx.optimize.least_squares.covariance_estimates
where the second factor is the Bérar-Lelann serial-correlation inflation [Bérar and Lelann, 1991]. Consecutive same-sign weighted residuals are summed coherently, \(\chi'^2 = \sum_{\mathrm{runs}} (\sum_{i\in\mathrm{run}} \delta_i)^2 \ge \chi^2\), because serially correlated neighbours do not carry independent information. The estimator is conservative: even white residuals land at an expected factor ≈1.51, so read it as an upper bound on the serial-correlation esd damage. Andreev’s serial-correlations figure of merit [Andreev, 1994] removes that bias by carrying the correlation into the minimised quantity itself.
Reported esds carry the inflation. The correlation matrix does not. It is the true Pearson matrix, so a genuinely degenerate pair reports \(|\rho| \approx 1\) and the 0.98 high-correlation guard means what it says. Values are quoted with two-significant-figure su’s per the IUCr convention [Schwarzenbach et al., 1989].
8.6. Esds of derived quantities¶
A bond length, an angle or a weight fraction is a function of the refined parameters rather than one of them, and its esd is the quadratic form
Source: rietx.model.geometry
over the whole covariance of (8.5): “the whole correlation matrix, not just the diagonal elements, should be included in the calculation” [McCusker et al., 1999]. Dropping the off-diagonal terms looks conservative and is not. Measured across the 88 interatomic distances of an 11-BM NAC structural refinement, the diagonal-only number runs from 0.86 to 1.41 times the full one, so it is as often too small as too large. Both are reported, so the difference is visible. It exists only where the coordinates refine, since a quantity depending on one free parameter has no off-diagonal term to drop.
A derived esd is absent rather than zero whenever that form cannot be evaluated honestly. Four cases qualify:
no covariance at all (an evaluate-only pass, a replayed history node);
no free parameter the quantity depends on;
a \(g^\top \mathrm{Cov}\, g\) that reaches zero by cancellation, detected as a variance below 1e-12 of its own form’s absolute terms. A symmetry-fixed 90° angle has exactly zero variance while its partials against \(x\), \(y\) and \(z\) do not vanish, so the quadratic form lands on roundoff;
an angle within 2.6e-03° of 0° or 180°, where the angle is a stationary point of the coordinates and the linearisation the propagation rests on does not hold at all.
The quantity itself is exact in all four cases. Only its uncertainty is withheld.
The partials are the restraint derivative chain of (9.11) evaluated at \(\sigma = w = 1\). A geometry row is a restraint row, so the two cannot drift apart, and the weight scale of (9.12) is kept out of them for the same reason.
8.7. A derived esd describes one basin¶
Every esd of (8.5) and (8.6) is built from the curvature of χ² at the converged point, so it describes the basin the fit stopped in. A trace phase’s weight fraction is where that stops being enough. The phase’s scale trades against its width: broadened far enough, its peaks become a hump the background shares, and the scale can then grow at almost no cost in χ². That ridge can hold several basins at one χ², each with ordinary curvature and a tight esd.
The width profile finds them. It pins one width parameter \(p\) of the phase on a grid, refits every other parameter at each value, and admits a point when
Source: rietx.strategy.fraction_profile.profile_fraction
with \(\chi^2\) the data’s own weighted sum, \(\chi^2_{\min}\) the lowest found (the fit’s included), and the reduced χ² and Bérar-Lelann factor those of the fit. Each pinned refit is the constrained hypothesis of Hamilton’s test with a single constraint [Hamilton, 1965], whose statistic is asymptotically \(\chi^2_1\), and \(c\) = 3.84 is its 95 % point. The two factors put the cut on the scale of the esds themselves, so on a single quadratic basin the admissible fractions reproduce \(W \pm 1.96\,\sigma_W\). Every admissible point’s fraction lies inside the 95 % profile interval for \(W\), so the range the profile reports is an inner bound, and a finer grid can only widen it. The profile is reported as undetermined when an admissible fraction lies more than 2.0 times \(1.96\,\sigma_W\) from the fit’s.
The pinned parameter is the width rather than the scale. At a pinned scale the refit can reach the hump basin only by broadening the phase below the noise, and there the stage holds the phase’s structure (a phase the data cannot see is held rather than bounded), so the scan would stop short of the basin it is looking for. At a pinned width the rest of the problem is close to linear in the scale, and each refit has one answer.
8.8. Staged strategy and series¶
Parameter groups are freed cumulatively in the IUCr-guideline order [McCusker et al., 1999]: scale and background first, then peak positions, then profile widths. The discrete model state is regenerated between stages and frozen within them. A series of related patterns (an in-situ ramp, a parametric sweep) is chained by warm starts, and the result is a parameter trajectory, path-dependent by construction. The chain can be run in both directions, and parameters the two runs disagree on are flagged. That is the one check separating a measured trajectory from an ordering artefact. True parametric refinement across patterns [Stinton and Evans, 2007] is out of scope.
8.9. Solvers¶
The default driver is scipy’s Trust Region Reflective. The bounded Levenberg-Marquardt alternative implements Coelho’s adaptive Marquardt constant [Coelho, 2018] over the bound-constrained conjugate-gradient solve of the normal equations [Coelho, 2005] (conjugate gradients per [Hestenes and Stiefel, 1952; Polak, 1971]), with the system diagonally pre-conditioned to \(A_{ii} = 1\). That pre-conditioning makes λ dimensionless and lets the published constants transfer. Conventions: \(A = J^\top J\), \(b = -J^\top r\), and the paper’s objective \(S = r^\top r\) is \(\chi^2\).
The driver earns its place on constraint vocabulary rather than speed. Box
bounds are enforced inside the linear solve, and linear inequalities on
functionals of θ are available as rows \(T\theta \ge 0\), the shape of the
Stephens positivity cone (Microstructure), which no per-parameter box
can express. An answer that pressed the cone says so. Per-stage truncation
counts are recorded, and a CONSTRAINT_ACTIVE diagnostic fires for the
answer-producing stage, which is the one signal that a declared constraint was
active. Speed was measured at 0.74–1.04× against TRF, as expected: the
normal-equation solve is a minority of the runtime, so solver work is
Amdahl-bounded at ≈1.25× here. Reading a paper against its own numbers works through two places where
the Coelho papers disagree with their own text, and how measurement settled
them.
8.10. Convergence¶
The Rietveld literature judges convergence by the final cycle’s parameter shifts against their own esds, rather than by a cost decrement [McCusker et al., 1999] (their §7):
Source: rietx.optimize.least_squares.run_least_squares
with the paper’s band \(\epsilon_{\mathrm{conv}} =\) 0.1 quoted, never tuned. Both sides are measured in external parameter units. \(\Delta\theta\) is decoded exactly through the transform chain of Parameterisation and constraints, and the esd is the chain-ruled physical one of (8.5), because at finite step size an internal-space ratio is a different number.
The measured value is Statistics.max_shift_over_esd, computed by the solver
from the final accepted step and copied onto the answer-producing stage’s
statistics. Nothing else derives it, and it gates nothing. A converged Trust
Region Reflective solve at ftol \(10^{-9}\) satisfies the criterion a fortiori
(measured \(\sim 3\times 10^{-4}\) on the synthetic LaB₆ round trip), so the
information is on the other branch: a stage stopped on its iteration budget
reports how far it was still moving in esd units, and the same stage starved to
one iteration measures \(\approx 14\). The value is withheld wherever it cannot
be measured: with no accepted step, no esds, an evaluate-only replay, or the
joint multi-pattern residual.
8.11. The fp64 floor¶
The residual used for cost and statistics, and the parameter solve and covariance, are always fp64 on host. A GPU backend may compute Jacobian columns in fp32. Three properties set that asymmetry.
The residual cancels. \(\sqrt{w}(y_{\mathrm{obs}} - y_{\mathrm{calc}})\) subtracts numbers of order 10⁵ counts to leave order 10², and fp32’s ~7 digits put an absolute error of order 10 counts into everything that reads it, which is ~10 % of the quantity just formed.
The solve squares the conditioning [Higham, 2002]:
Source: rietx.backend.linalg64
so a routine Rietveld \(\operatorname{cond}(J) \sim 10^4\) leaves \(\operatorname{cond}(J^\top J) \sim 10^8\), which fp32 cannot invert at all. The bounded LM forms \(J^\top J\) explicitly, and is the direct illustration.
Columns are relative-accuracy tolerant. A column enters through a descent direction and a curvature estimate, and the trust region re-measures every step against a fresh fp64 cost. Measured on real hardware, an Apple-GPU refinement with every column in fp32 lands 3.5×10⁻⁸ Å from the numpy fp64 cell.
8.12. From fit to report¶
The FitReport reads the converged state in three layers. Layer 0 is model-free:
cumulative-χ² breakpoints localise where misfit lives [David, 2004], and
unindexed peaks are flagged against the tick positions of every emission line
(an impurity must clear 8.0σ). Layer 1 attributes per-region
misfit to profile shape derivatives, under four gates (resolvability on the
scale-normalised Gram, a validity radius, local-χ² significance, global
maturity), so a collinear pair is declared non-separable rather than resolved
into a confident wrong singleton. Layer 2 turns attributions into actions. The
actions are advisory and carry no statistical test of their own. Whether a
freed parameter pays for itself is ΔBIC [Schwarz, 1978], predicted before
the fit by Refinement.suggest and measured after it by
report.compare_freed (next section). The report’s thresholds are versioned
(currently 1.10).
8.13. Which parameter to free next¶
Refinement.suggest() ranks every held-but-refinable parameter by the χ²
reduction one Gauss-Newton solve would obtain from freeing it, at the cost
of a single Jacobian evaluation and no solve. With \(F\) the currently-free
columns, \(P_F\) the orthogonal projector onto their span, \(r\) the weighted
residual and \(J_j\) a held parameter’s column,
Source: rietx.optimize.statistics.one_parameter_gains
which is Rao’s score statistic [Rao, 1948] applied to the linearised model, computed through the Frisch-Waugh-Lovell projection identity [Frisch and Waugh, 1933; Lovell, 1963]. It is exactly the drop in \(\sum w\Delta^2\) that a least-squares solve of \([F \mid J_j]\) achieves over \(F\) alone. It is invariant under any rescaling of the column, so no per-parameter step heuristics are needed. At a converged minimum \(J^\top r \approx 0\) makes every gain vanish. GSAS-II answers the same recipe problem [Toby, 2024] by ±δ finite differences with per-type δ heuristics and a sign-consistency test, because its analytic derivatives are locked inside Hessian assembly. Exact columns at the current state make all three workarounds unnecessary.
Under the null hypothesis a gain is distributed as
\(\chi^2_1 \cdot \chi^2_{\mathrm{red}}\), so a candidate is quotable only above
a noise floor of 9.0 · max(χ²_red, 1). That is the 3σ point
of \(\chi^2_1\), with the same floor-at-one convention as the covariance scale.
Two gates keep the ranking honest, the Layer-1 discipline one call over. The
free block absorbs some candidates’ columns, and those are reported
non-separable and carry no score; the same projection caps the \(1/(1-R^2)\)
inflation of near-collinear gains. Candidates whose projected columns are
pairwise indistinguishable come back as one unresolved group carrying a joint
gain, which is a tie. As with indexing there is no .best, and
best_or_none() answers None whenever the evidence does not choose one
parameter.
A gain that clears the floor is then read as a model-selection answer, ΔBIC [Schwarz, 1978], and that needs a premise the gain does not: Schwarz’s \(N\) counts independent observations. A powder residual is serially correlated, and at raw channel counts the reward term \(N\ln(\chi^2_r/\chi^2_f)\) outvotes the \(\ln N\) penalty for almost any improvement. The fit already measures the correlation, as the Bérar-Lelann factor \(f\) every esd is inflated by [Bérar and Lelann, 1991], and the penalty is charged at
Source: rietx.optimize.statistics.effective_sample_size
the count at which, for one added parameter, the reward term is that parameter’s \(t^2\) at its inflated esd, so ΔBIC is positive exactly when \(t^2 > \ln N_{\mathrm{eff}}\). It is a heuristic and not a theorem, and it is conservative the way the esds are: white residuals still give \(f \approx 1.51\). Hamilton’s test [Hamilton, 1965] takes the same \(N_{\mathrm{eff}}\), so the two verdicts are read off one count. Measured on four ~49 500-channel synchrotron fits of one occupancy, raw \(N\) gave ΔBIC +36 to +211 to a parameter each fit’s own esd put within 0.76-1.89σ of zero, and \(N_{\mathrm{eff}}\) turned all four negative.
After the fit the same count prices the measured gain. report.compare_freed
takes the fit without a parameter and the fit with it, and returns
(8.11) beside each freed parameter’s \(t\). The two agree for one
parameter, so their disagreement is information: a multi-parameter block whose
joint ΔBIC is positive while no member’s \(t\) clears 2, or an \(f\) that moved
between the fits.