6  Mixed-Precision Iterative Refinement: Experiment Plan

With the corrected third task and the decision to keep the primary condition sweeps deterministic, the complete plan becomes the following.

The best place to begin remains the experiment-output infrastructure. Once the full suite is running, incomplete metadata, inconsistent filenames, or unclear seed handling would force us to regenerate results.

6.1 Overall experiment plan

The experiments can be divided into the four objectives from the project brief and a fifth group covering our FP16 robustness extension.

Group Purpose Main variables Main outputs
A. Convergence histories Examine refinement iteration by iteration precision triple, selected \(\kappa(A)\) forward error, backward error, relative correction
B. Condition-number sweeps Locate the convergence boundary \(\kappa(A)\), precision triple final/best errors, iterations, status
C. Residual-precision comparison Validate the last two rows of Table 7.1 single–double–double vs. single–double–quad forward and backward error vs. \(\kappa(A)\)
D. Direct-solve comparison Measure the benefit of refinement mixed IR vs. direct solve in working precision error vs. \(\kappa(A)\)
E. FP16 robustness extension Demonstrate range checks and residual scaling scaling on/off, \(\kappa(A)\) histories, conversion diagnostics, termination status

For the primary experiments, the default random seeds should remain fixed. This isolates the effect of the condition number and keeps the plots straightforward.

6.2 Group A: convergence histories

These experiments show how the algorithm behaves over individual refinement iterations.

Representative precision triples

A suitable selection is:

  • FP64–FP64–FP64: same-precision baseline;
  • FP32–FP64–FP64: classical mixed precision;
  • FP32–FP64–FP128: higher residual precision;
  • FP16–FP64–FP64: half-precision factorization;
  • FP16–FP64–FP128: principal FP16 configuration;
  • bfloat16–FP64–FP128: exponent-range comparison;
  • FP8–FP64–FP128: optional stress test.

Not every precision triple needs to appear in the main report. A representative subset can be selected after inspecting the results.

Condition numbers

For each factorization precision, choose condition numbers:

  • clearly inside the predicted convergence region;
  • near \(\kappa(A)u_f\approx 1\);
  • outside the predicted region.

Plots

For every selected configuration, plot against iteration:

  1. relative forward error;
  2. normwise backward error;
  3. relative correction,

\[ \frac{\lVert d_k\rVert_2}{\lVert x_k\rVert_2}. \]

The seeds remain fixed within and across comparisons. This ensures that different precision triples solve exactly the same problems.

6.3 Group B: general condition-number sweeps

These experiments study how each precision configuration behaves as conditioning deteriorates.

Use a logarithmic condition-number grid, supplemented with additional points near the expected convergence boundary \(1/u_f\).

For example:

  • FP16: concentrate around \(10^2\)–\(10^4\);
  • bfloat16: concentrate around \(10^1\)–\(10^3\);
  • FP32: concentrate around \(10^6\)–\(10^9\).

Seed policy

Only \(\kappa(A)\) should vary. Use the default fixed seeds for:

  • the matrix generator;
  • the right-hand side;
  • any generated orthogonal transformations.

Provided the generator is reset deterministically for each value of \(\kappa(A)\), this keeps the matrix orientation and right-hand side controlled while changing the prescribed spectrum.

Multiple seeds are not required. They can be added later only if the primary sweep produces unexpectedly irregular results.

A concise report formulation is:

The random seeds were fixed throughout the sweep. Thus, the matrix structure and right-hand side remained controlled while the prescribed condition number was varied. This isolates the effect of conditioning and permits direct comparisons between precision configurations.

Recorded quantities

For each run, record:

  • requested condition number;
  • measured condition number;
  • termination status;
  • final forward and backward errors;
  • best attained forward and backward errors;
  • iterations at which the best errors occurred;
  • number of completed refinement updates;
  • final relative correction.

Plots

Create:

  1. final forward error vs. \(\kappa(A)\);
  2. final backward error vs. \(\kappa(A)\);
  3. best attained error vs. \(\kappa(A)\);
  4. iteration count vs. \(\kappa(A)\);
  5. termination status vs. \(\kappa(A)\).

Show the theoretical condition

\[ \kappa(A)u_f=1 \]

as a vertical reference line.

Because the seeds are fixed, each configuration produces a single curve. No medians, ranges, or quantile bands are needed.

6.4 Group C: single–double–double versus single–double–quad

This replaces the earlier broad residual-precision study. It directly implements the experiment requested by the supervisor and validates the final two rows of Table 7.1 in the Carson–Higham paper.

Precision configurations

Compare:

  1. single–double–double:

    \[ u_f=u_{\mathrm{single}}, \qquad u=u_r=u_{\mathrm{double}}; \]

  2. single–double–quad:

    \[ u_f=u_{\mathrm{single}}, \qquad u=u_{\mathrm{double}}, \qquad u_r=u_{\mathrm{quad}}\approx u^2. \]

Only the residual precision changes. Both variants must use:

  • the same matrix;
  • the same right-hand side;
  • the same fixed seeds;
  • the same stopping criterion;
  • the same maximum iteration count;
  • the same remaining algorithm options.

Theoretical expectation

The general limiting forward-error bound is

\[ \frac{\lVert x-\widehat{x}\rVert_\infty} {\lVert x\rVert_\infty} \lesssim 4p\,u_r\,\operatorname{cond}(A,x)+u. \]

For single–double–double, \(u_r=u\), giving

\[ 4p\,u\,\operatorname{cond}(A,x)+u. \]

The conditioning-dependent term generally dominates, so the expected forward error is

\[ O\!\left(u\,\operatorname{cond}(A,x)\right). \]

For single–double–quad, \(u_r\approx u^2\), giving

\[ 4p\,u^2\,\operatorname{cond}(A,x)+u. \]

Within the convergence region of the single-precision factorization, the first term is sufficiently small that the working-precision term dominates. Therefore, the expected forward error is

\[ O(u). \]

Both variants should attain a backward error of order

\[ O(u). \]

Thus, the quad residual should improve the limiting forward error substantially, while both variants should produce approximately double-precision backward errors.

The remembered equation

\[ 4pu^2\operatorname{cond}(A,x)+u\lesssim 4pu \]

is correct, but equation (8.5) appears in the later GMRES-IR analysis. For ordinary LU-IR, the general bound above provides the more direct explanation.

Condition-number sweep

A suitable initial range is approximately

\[ 10^0\leq\kappa(A)\leq10^9, \]

with additional points between \(10^6\) and \(10^8\), near the single-precision convergence boundary.

For IEEE single precision,

\[ u_f=2^{-24}\approx5.96\times10^{-8}, \qquad u_f^{-1}\approx1.68\times10^7. \]

This is only an approximate boundary. The actual transition can occur earlier because it also depends on the stability of the computed factorization.

Required plots

The supervisor requested two log–log plots, both with condition number on the horizontal axis:

  1. forward error for both variants in one plot;
  2. backward error for both variants in one plot.

Expected behaviour:

Plot Expected result
Forward error single–double–double grows approximately as \(u\,\operatorname{cond}(A,x)\); single–double–quad remains near \(u\)
Backward error both variants remain near \(u\) while refinement converges

Useful reference lines are:

  • \(u\), representing the double-precision floor;
  • \(u\kappa(A)\), illustrating the expected single–double–double trend;
  • \(1/u_f\), marking the approximate convergence boundary.

Measurement details

Table 7.1 uses \(\kappa_\infty(A)\), whereas the forward-error bound contains \(\operatorname{cond}(A,x)\). They are related but not identical.

Each run should therefore record:

  • requested/generated condition number;
  • measured \(\kappa_2(A)\);
  • preferably measured \(\kappa_\infty(A)\);
  • \(\operatorname{cond}(A,x)\), if practical;
  • forward and backward errors;
  • termination status;
  • iteration count.

The main horizontal axis should preferably use measured \(\kappa_\infty(A)\), matching Table 7.1. Error measurements should use sufficiently high precision, ideally FP256.

This experiment can be summarized as:

A condition-number sweep compares single–double–double and single–double–quad iterative refinement. Both variants are expected to attain double-precision backward error. The quad-residual variant should also attain double-precision forward error, whereas the double-residual variant should attain forward error of order \(u\,\operatorname{cond}(A,x)\).

6.5 Group D: comparison with a direct solve

For every chosen working precision \(u\), compare:

  • a single LU solve performed entirely in \(u\);
  • mixed-precision iterative refinement whose iterates are stored in \(u\).

The principal question is whether mixed IR attains working-precision accuracy while performing the factorization in a cheaper, lower precision.

Plots

Create:

  1. relative forward error vs. \(\kappa(A)\);
  2. backward error vs. \(\kappa(A)\);
  3. optionally, iteration count for mixed IR.

Both methods must solve identical problems with fixed seeds.

6.6 Group E: FP16 robustness extension

This should form a separate report section because it goes beyond the original project requirements.

Two distinct FP16 range problems should be presented separately:

  • factorization-input overflow;
  • residual-conversion underflow.

Residual scaling mitigates the second problem. It does not repair overflow in the original factorization input.

E1. Controlled scaling comparison

Use the existing controlled problem:

  • random SPD matrix;
  • \(n=100\);
  • \(\kappa(A)=10\);
  • random-normal right-hand side;
  • FP16–FP64–FP128;
  • identical default seeds;
  • residual scaling off and on.

Only scale_residual should change.

Principal plots

Plot:

  1. forward error vs. iteration;
  2. backward error vs. iteration;
  3. relative correction vs. iteration;
  4. residual infinity norm vs. iteration.

Diagnostic plots

Also plot:

  1. number of residual components zeroed during conversion to FP16;
  2. smallest nonzero residual component vs. iteration, with the FP16 subnormal threshold shown.

This experiment should demonstrate the observed mechanism:

  • without scaling, increasingly many residual components are lost during conversion and the algorithm stagnates;
  • with scaling, the residual is normalized before conversion;
  • conversion-to-zero is avoided or greatly reduced;
  • the algorithm reaches approximately FP64 accuracy.

E2. Broader scaling validation

A simple broader validation can also keep the default seeds fixed and vary only \(\kappa(A)\). Compare scaling on and off over several condition numbers for which the FP16 factorization input remains finite.

Record:

  • termination status;
  • final and best forward errors;
  • iteration count;
  • maximum number of zeroed residual components;
  • total number of zeroing events;
  • residual-reduction factor.

Create:

  1. final forward error vs. \(\kappa(A)\);
  2. iterations vs. \(\kappa(A)\);
  3. maximum zeroed-component count vs. \(\kappa(A)\);
  4. termination status vs. \(\kappa(A)\).

Multiple seeds are optional here as well. They would strengthen a claim of general robustness, but they are not necessary for the primary report. They should only be added if time permits or if the fixed-seed results appear atypical.

E3. FP16 overflow handling

Input overflow should be demonstrated separately.

Record:

\(\kappa(A)\) Maximum \(\lvert A_{ij}\rvert\) FP16 input finite Termination status

The expected robust behaviour is termination with a status such as factorization_input_non_finite, rather than a crash or propagation of invalid values.

A table or categorical status plot is more suitable than a convergence-history plot.

6.7 Controlling the experiment space

Testing every combination of matrix family, right-hand side, precision triple, condition number, seed, dimension, and scaling mode would produce an unnecessarily large experiment suite.

A focused selection is preferable:

  • primary matrix family: random SPD;
  • primary right-hand side: random-normal;
  • primary dimension: \(n=100\);
  • fixed default seeds;
  • rotated SPD: supplementary controlled case;
  • random SVD: supplementary nonsymmetric case;
  • alternative right-hand sides: optional sensitivity study;
  • multiple seeds: optional robustness check only.

This keeps the report centred on the requested effects.

6.8 Common metadata design

The experimental utilities should be based on one shared experiment description.

Conceptually:

struct PrecisionNames {
    std::string factor;
    std::string work;
    std::string residual;
    std::string measure;
};

struct ExperimentDescription {
    ExperimentKind kind;
    MatrixFamily matrix_family;
    std::size_t dimension;
    PrecisionNames precisions;
};

The metadata writer should additionally receive:

  • TestProblemOptions;
  • MixedIROptions<T_work>;
  • the condition number or condition-number grid;
  • the variant name, such as scaled or unscaled.

Experiment identity

Record:

  • experiment kind;
  • matrix family;
  • dimension \(n\);
  • factorization precision;
  • working precision;
  • residual precision;
  • measurement precision;
  • algorithm variant.

Test-problem options

Record:

  • right-hand-side mode;
  • matrix_seed_u;
  • matrix_seed_v;
  • vector_seed;
  • rotation_theta;
  • prescribed condition number;
  • measured condition numbers.

Recording the seeds is important even though they remain fixed. It makes the deterministic problem construction reproducible.

Algorithm options

Record:

  • max_iterations;
  • effective relative-correction tolerance;
  • detect_divergence;
  • divergence_growth_factor;
  • divergence_growth_steps;
  • store_iterates;
  • scale_residual;
  • record_residual_diagnostics.

The effective tolerance should be written, not merely a sentinel value such as a nonpositive number meaning “use the default.”

Results

Record:

  • termination status;
  • convergence flag, if retained;
  • number of completed updates;
  • final relative correction;
  • final forward and backward errors;
  • best forward and backward errors;
  • iterations of the best errors.

CSV representation

I recommend repeating invariant metadata as ordinary columns in every CSV row. This is redundant but useful:

  • every file remains self-contained;
  • pandas.read_csv() works without a custom parser;
  • files can be concatenated directly;
  • the experiment remains identifiable if the file is renamed.

A JSON sidecar containing the Git revision, compiler information, timestamp, and build configuration could be added later, but it is not required for the first refactoring.

6.9 Automatic filenames

The filename should encode the experiment identity, but not every option.

A suitable pattern is:

<experiment>__<matrix>__<rhs>__uf-<factor>__u-<work>__ur-<residual>__<variant>.csv

Examples:

convergence__random-spd__random-normal__uf-fp32__u-fp64__ur-fp128.csv
residual-precision__random-spd__random-normal__uf-fp32__u-fp64__ur-comparison.csv
residual-scaling__random-spd__random-normal__uf-fp16__u-fp64__ur-fp128__comparison.csv

The filename generator and metadata writer should consume the same ExperimentDescription. This prevents the filename and CSV contents from disagreeing.

Precision names can initially be supplied explicitly. They can later be generated through traits such as precision_name<T>().

6.10 Output-directory organization

A moderate hierarchy is sufficient:

results/
├── raw/
│   ├── convergence/
│   ├── condition_sweeps/
│   ├── residual_precision/
│   ├── direct_comparison/
│   └── robustness/
│       ├── residual_scaling/
│       └── range_failures/
└── plots/
    ├── convergence/
    ├── condition_sweeps/
    ├── residual_precision/
    ├── direct_comparison/
    └── robustness/
        ├── residual_scaling/
        └── range_failures/

Plotting scripts should mirror the raw-data directory structure. This avoids accumulating many ambiguously named plots in one flat directory.