5 Implementing the Experiments & Robustness Improvements

We want to plot the forward error and the backward error as function of the refinement step \(i\).

forward error:

\[ \text{forward error}_i = \frac{\norm{x_i - x_{\text{true}}}}{\norm{x_{\text{true}}}} \]

Backward error:

\[ \eta_{i} = \frac{\norm{b - Ax_i}}{\norm{A}\norm{x_i} + \norm{b}} \]

in mixed_ir.hpp the result struct currently only stores rel_correction vector, that records

\[ \frac{\norm{d_i}}{\norm{x_i}} \]

We have to modify mixed_ir to store the iterates \(x_i\), so that we can compute forward and backward errors. We add the following member variable to the MixedIRResult struct:

std::vector<hdnum::Vector<T_work>> iterates;

We append the solutions to the iterates vector before the loop and after each loop iteration with the line:

if (options.store_iterates) {
    result.iterates.push_back(result.x);
}

5.1 Temporary Matrix Generator

We start with a simple controlled family, a rotated SPD matrix:

\[ A = Q\Lambda Q^T \]

where \(\Lambda\) has eigenvalues \(1/\sqrt{\kappa} \dots \sqrt{\kappa}\). Then

\[ \kappa_2(A) = \kappa \]

This family of matrices give a known true (exact) solution:

\[ x_{\text{true}} = (1, \dots, 1)^T \]

So it can be used temporarily to test the pipeline.

The implementation is in

code/include/test_matrices.hpp

5.2 Matrix Generator Supervisor Version

The temporary generator constructs

\[ A = Q_{\text{block}}\Lambda Q_{\text{block}}^T \]

where \(Q_{\text{block}}\) consists only of independent \(2 \times 2\) rotations. Consequently \(A\) is block diagonal:

\[ A = \begin{pmatrix} A_1 & & & 0 \\ & A_2 & & \\ & & \ddots & \\ 0 & & & A_m \end{pmatrix}, \qquad A_i \in \mathbb{R}^{2 \times 2}. \]

Even though the entire block-diagonal matrix had condition number \(\kappa\), the individual \(2 \times 2\) blocks is only approximately:

\[ \kappa^{1/(n-1)} \]

Supervisors randspd constructs

\[ A = Q\Lambda Q^T \]

where \(Q\) is a dense random orthogonal matrix. Consequently, \(A\) is generally dense, and every variable is coupled to every other variable, producing a genuinelly ill-conditioned system.

Nonsymmetric Generator

Iterative refinement with LU is not restricted to SPD matrices, so randsvd is also important. The wrapper

template<class T>
LinearSystem<T>
make_random_svd_problem(std::size_t n,
                        double kappa,
                        unsigned int seed_u = 42,
                        unsigned int seed_v = 137)
{
    hdnum::DenseMatrix<T> A(n, n);
    hdnum::Vector<T> x_true(n);
    hdnum::Vector<T> b(n);

    hdnum::randsvd(
        A,
        scalar_cast<T>(kappa),
        seed_u,
        seed_v
    );

    for (std::size_t i = 0; i < n; ++i) {
        x_true[i] = T(1);
    }

    A.mv(b, x_true);

    return LinearSystem<T>{
        A,
        b,
        x_true,
        kappa
    };
}

produces

\[ A = U\Sigma V^T \]

with singular values \(1, \dots 1/\kappa\). It is nonsymmetric and a more general test of the iterative LU algorithm.

5.3 Configurable test-problem Construction

The test-matrix infrastructure contained in ‘test_matrices.hpp’ supports three different ways of constructing the right-hand side \(b\) and solution of

\[ Ax = b \]

The matrix generator and the right-hand-side strategy are separated:

A matrix generator constructs \(A\)
A shared complete_problem function constructs \(b\).
The reference solution is computed in a separate, higher precision.

Three Right-hand-side Modes

There are three availabe strategies represented by:

enum class RightHandSideMode {
    ones_solution,
    random_sign_solution,
    random_normal_rhs
};

ones_solution:

\[ \begin{aligned} x &= (1, \dots, 1)^T \\ b &= \text{fl}(Ax) \end{aligned} \]

random_sign_solution:

Construct a reproducible random vector with

\[ x_i \in \{-1, 1\} \]

and form

\[ b = \text{fl}(Ax) \]

random_normal_rhs:

Here we generate the right-hand-side directly (randomly):

\[ b_i \sim \mathcal{N}(0, 1) \]

In all these cases the reference (true) solution is obtained by solving

\[ Ax_{\text{ref}} = b \]

in a high precision (called measure or reference precision).

Test-problem Options

The choices and random sees are collected in:

struct TestProblemOptions {
    RightHandSideMode rhs_mode =
        RightHandSideMode::ones_solution;

    unsigned int matrix_seed_u = 42;
    unsigned int matrix_seed_v = 137;
    unsigned int vector_seed = 2718;

    double rotation_theta = 0.3;
};

rhs_mode: selects how \(b\) is constructed
matrix_seed_u: seed for the first random orthogonal matrix (used by randspd)
matrix_seed_v: seed for the second random orthogonal matrix (used by randsvd)
vector_seed: seed for the random-sign solution or random-normal right-hand side.
rotation_theta: rotation angle used by the structured rotated-SPD generator.

The default mode is ones_solution.

Linear-system Representation

The system stores the problem and reference solution (potentially in different precisions):

template<class T_data, class T_reference = T_data>
struct LinearSystem {
    hdnum::DenseMatrix<T_data> A;
    hdnum::Vector<T_data> b;
    hdnum::Vector<T_reference> x_true;
    double kappa;
};

T_data: the precision in which \(A\) and \(b\) are stored, normally T_work.
T_reference: precision used for the reference solution, normally T_measure.

Shared `complete_problem` function

We have three matrix generators,

make_rotated_spd_problem
make_random_spd_problem
make_random_svd_problem

each created only \(A\) and delegates the remaining work to create_problem, the signature:

template<class T_data, class T_reference>
LinearSystem<T_data, T_reference>
complete_problem(
    hdnum::DenseMatrix<T_data> A,
    double kappa,
    const TestProblemOptions& options);

The logical structure is:

if (options.rhs_mode ==
    RightHandSideMode::ones_solution) {

    // Construct x = (1,...,1).
    // Compute b = A*x in T_data.
}
else if (options.rhs_mode ==
         RightHandSideMode::random_sign_solution) {

    // Construct reproducible x_i in {-1,1}.
    // Compute b = A*x in T_data.
}
else {
    // Construct reproducible b_i ~ N(0,1)
    // directly in T_data.
}

After \(b\) is constructed the reference solution is computed:

auto x_reference =
    high_precision_solve<T_reference>(A, b);

The reference solution is recomputed

Even in the first two modes, even if we explicitly construct \(x\), we recompute

\[ x_{\text{true}} = \text{solve}_{\text{T}_{\text{ref}}}(A, b) \]

because the matrix-vector product

\[ b = \text{fl}(Ax_{\text{constructed}}) \]

is founded. Therefore, the exact solution of the stored pair \((A, b)\) can differ slightly from \(x_{\text{constructed}}\).

High-precision reference solve

The helper

template<class T_reference, class T_data>
hdnum::Vector<T_reference>
high_precision_solve(
    const hdnum::DenseMatrix<T_data>& A,
    const hdnum::Vector<T_data>& b);

contained in reference_solve.hpp computes a reference solution using hdnums LU-solver functions in T_reference precision.

For example, T_reference = hdnum::FP256 provides a sufficiently accurate reference solution.

Matrix-generator interfaces

As mentioned there are three generators, all with compatible interfaces:

template<class T_data, class T_reference = T_data>
LinearSystem<T_data, T_reference>
make_rotated_spd_problem(
    std::size_t n,
    double kappa,
    const TestProblemOptions& options = {});

template<class T_data, class T_reference = T_data>
LinearSystem<T_data, T_reference>
make_random_spd_problem(
    std::size_t n,
    double kappa,
    const TestProblemOptions& options = {});

template<class T_data, class T_reference = T_data>
LinearSystem<T_data, T_reference>
make_random_svd_problem(
    std::size_t n,
    double kappa,
    const TestProblemOptions& options = {});

Each generator

constructs its particular matrix \(A\)
calls complete_problem
returns LinearSystem<T_data, T_reference> struct.

Which enables plugging in different matrix families in an experiment.

Usage Examples

Define the precisions:

using T_work = hdnum::FP64;
using T_measure = hdnum::FP256;

all-ones \(x\):

mpir::TestProblemOptions problem_options;

problem_options.rhs_mode =
    mpir::RightHandSideMode::ones_solution;

auto problem =
    mpir::make_random_spd_problem<T_work, T_measure>(
        n,
        kappa,
        problem_options
    );

random-sign ones \(x\):

mpir::TestProblemOptions problem_options;

problem_options.rhs_mode =
    mpir::RightHandSideMode::random_sign_solution;

problem_options.vector_seed = 2718;

auto problem =
    mpir::make_random_spd_problem<T_work, T_measure>(
        n,
        kappa,
        problem_options
    );

random-normal \(b\):

mpir::TestProblemOptions problem_options;

problem_options.rhs_mode =
    mpir::RightHandSideMode::random_normal_rhs;

problem_options.vector_seed = 1;

auto problem =
    mpir::make_random_svd_problem<T_work, T_measure>(
        n,
        kappa,
        problem_options
    );

only the generator call needs to change, the resulting problems are used identically:

auto result =
    mpir::mixed_ir<T_factor, T_work, T_residual>(
        problem.A,
        problem.b,
        ir_options
    );

5.4 Error Metrics Helper

We implement functions to compute

\(\normsub{\cdot}{\infty}\) for matrices and vectors,
relative forward error with the infinity norm
relative backward error with the inifinity norm

contained in the file

code/include/error_metrics.hpp

5.5 First Convergence Test

We test our pipeline with a first test that uses the precision tripple:

\(u_f\) : FP16
\(u\): FP64
\(u_r\): FP128

For matrices of size \(100 \times 100\) with \(\kappa\) values: \(1, 10, 100, 1000, 10000\).

The test writes for each matrix the convergence history to csv file. Specifically it writes the values:

kappa,
iteration,
relative forward error
relative backward error
relative correction

in a single lines for each iteration step, for each matrix. This gives us the convergence histories.

It is implemented in the file

code/experiments/exp_convergence.cc

Following was the result of the experiment:

kappa,iteration,forward_error_inf,backward_error_inf,rel_correction
1,0,0,0,0
1,1,0,0,0
10,0,0.000976562,0.000441793,0
10,1,9.53674e-07,4.31651e-07,0.00047591
10,2,5.96046e-08,5.1214e-09,3.1641e-07
10,3,5.96046e-08,5.1214e-09,1.68587e-08
10,4,5.96046e-08,5.1214e-09,1.68587e-08
10,5,5.96046e-08,5.1214e-09,1.68587e-08
10,6,5.96046e-08,5.1214e-09,1.68587e-08
10,7,5.96046e-08,5.1214e-09,1.68587e-08
10,8,5.96046e-08,5.1214e-09,1.68587e-08
10,9,5.96046e-08,5.1214e-09,1.68587e-08
10,10,5.96046e-08,5.1214e-09,1.68587e-08
10,11,5.96046e-08,5.1214e-09,1.68587e-08
10,12,5.96046e-08,5.1214e-09,1.68587e-08
10,13,5.96046e-08,5.1214e-09,1.68587e-08
10,14,5.96046e-08,5.1214e-09,1.68587e-08
10,15,5.96046e-08,5.1214e-09,1.68587e-08
10,16,5.96046e-08,5.1214e-09,1.68587e-08
10,17,5.96046e-08,5.1214e-09,1.68587e-08
10,18,5.96046e-08,5.1214e-09,1.68587e-08
10,19,5.96046e-08,5.1214e-09,1.68587e-08
10,20,5.96046e-08,5.1214e-09,1.68587e-08
100,0,0.000976562,0.000481752,0
100,1,9.53674e-07,4.70694e-07,0.000444832
100,2,1.78814e-07,1.42242e-09,2.52178e-07
100,3,1.78814e-07,1.42242e-09,0
1000,0,0.000976562,0.000460731,0
1000,1,9.53674e-07,3.07722e-07,0.000591975
1000,2,8.9407e-07,6.88472e-10,4.02364e-07
1000,3,8.9407e-07,5.65293e-10,1.3328e-08
1000,4,8.9407e-07,5.65293e-10,1.19209e-08
1000,5,8.9407e-07,5.65293e-10,1.19209e-08
1000,6,8.9407e-07,5.65293e-10,1.19209e-08
1000,7,8.9407e-07,5.65293e-10,1.19209e-08
1000,8,8.9407e-07,5.65293e-10,1.19209e-08
1000,9,8.9407e-07,5.65293e-10,1.19209e-08
1000,10,8.9407e-07,5.65293e-10,1.19209e-08
1000,11,8.9407e-07,5.65293e-10,1.19209e-08
1000,12,8.9407e-07,5.65293e-10,1.19209e-08
1000,13,8.9407e-07,5.65293e-10,1.19209e-08
1000,14,8.9407e-07,5.65293e-10,1.19209e-08
1000,15,8.9407e-07,5.65293e-10,1.19209e-08
1000,16,8.9407e-07,5.65293e-10,1.19209e-08
1000,17,8.9407e-07,5.65293e-10,1.19209e-08
1000,18,8.9407e-07,5.65293e-10,1.19209e-08
1000,19,8.9407e-07,5.65293e-10,1.19209e-08
1000,20,8.9407e-07,5.65293e-10,1.19209e-08
10000,0,0.000976562,0.000375237,0
10000,1,3.33786e-06,4.73515e-08,0.000488379
10000,2,9.53674e-07,3.8223e-10,5.28939e-07
10000,3,9.53674e-07,1.29324e-10,1.03238e-08
10000,4,9.53674e-07,1.29324e-10,0

5.6 Robustness Improvements to `mixed_ir`

We encountered two failures at \(\kappa(A)=10^9\). Both caused a floating-point exception, but for different reasons.

1. FP16–FP64–FP128: overflow during input conversion

Problem

We converted \(A\) and \(b\) from FP64 to the FP16 factorization precision. FP16 has a maximum finite value of approximately

\[ 65504. \]

Some entries exceeded this limit. The converted data therefore contained non-finite values. LU factorization and the initial solve propagated these values into \(x_0\).

The program crashed when compute_residual attempted to convert the non-finite iterate to the GMP-backed FP128 type. This type cannot represent NaN or infinity.

Detection

The GDB backtrace located the crash in the GMP conversion performed by compute_residual.

We traced the invalid values backward:

\[ A_f,b_f \longrightarrow LU_f \longrightarrow x_{0,f} \longrightarrow x_0 \longrightarrow \text{GMP conversion}. \]

Inspection showed that the first non-finite values appeared during conversion to FP16.

The bfloat16 experiment supported this conclusion. Bfloat16 has a wider exponent range and did not overflow at the same point.

Solution

We chose to detect the overflow and terminate safely. We did not implement scaling or equilibration.

The algorithm now terminates with factorization_input_non_finite if \(A_f\) or \(b_f\) contains a non-finite value.

Scaling could prevent this specific overflow, but it would not change \(\kappa(A)\). It would therefore not make FP16 iterative refinement reliable at \(\kappa(A)=10^9\).

Source-Code Changes

Termination statuses

We introduced statuses that distinguish the main outcomes:

enum class MixedIRStatus {
    converged,
    max_iterations,
    factorization_input_non_finite,
    non_finite,
    diverged
};

These represent:

successful convergence;
reaching the iteration limit;
non-finite factorization input;
a non-finite value produced later;
detected divergence.

Finite-value checks

We added all_finite helpers for vectors and matrices.

IEEE-like types and CPFloat values are checked with std::isfinite. GMP-backed hdnum::FP<m> values require separate handling because they cannot represent NaN or infinity.

The main rule is:

Check IEEE-like values before converting them to a GMP-backed type.

We kept all checks inside mixed_ir. We did not modify HDNUM, lowprec_cpfloat.hh, compute_residual, or the triangular solvers.

Guard placement

We added finite-value checks at the following points:

Location	Purpose
After converting \(A,b\) to \(A_f,b_f\)	Detect factorization-input overflow
After the initial solve	Detect an invalid initial solution
After converting \(x_{0,f}\) to working precision	Validate the initial iterate
Before `compute_residual`	Protect the GMP conversion
After computing \(r_r\)	Support generic residual types
After converting \(r_r\) to \(r_f\)	Detect narrowing overflow
After solving for \(d_f\)	Detect an invalid correction
After converting \(d_f\) to \(d_w\)	Validate the working-precision correction
After converting the relative correction to `double`	Detect norm or conversion overflow
After forming \(x_{k+1}\)	Detect overflow during the update

We now form the new iterate as a candidate:

hdnum::Vector<T_work> x_next(result.x);

for (std::size_t i = 0; i < n; ++i) {
    x_next[i] += d_w[i];
}

if (!all_finite(x_next)) {
    result.status = MixedIRStatus::non_finite;
    return result;
}

result.x = x_next;

This preserves the last valid iterate if the update fails. We increment result.iterations only after accepting the update.

Divergence options

We added three options:

bool detect_divergence = true;

T_work divergence_growth_factor = T_work(10);

std::size_t divergence_growth_steps = 3;

They control:

whether divergence detection is enabled;
the minimum excessive-growth factor;
the required number of consecutive growth steps.

Valid user-supplied values should satisfy:

divergence_growth_factor > T_work(1)
divergence_growth_steps > 0

Correction-norm history

We added diagnostic storage to MixedIRResult:

std::vector<double> correction_norms;

It stores each finite value of \(\lVert d_k\rVert_2\).

The history also includes a correction that triggers divergence but is not applied. We store the norms even when divergence detection is disabled.

Detector state and comparison

Before the refinement loop, we introduced:

T_work previous_correction_norm = T_work(0);
bool have_previous_correction_norm = false;
std::size_t consecutive_growth_steps = 0;

After computing and validating \(d_w\), we calculate:

const T_work correction_norm_w = hdnum::norm(d_w);

We then test for excessive growth:

const bool excessive_growth =
    previous_correction_norm > T_work(0)
    &&
    correction_norm_w
        > options.divergence_growth_factor
            * previous_correction_norm;

Excessive growth increments the counter. Any other result resets it. Once the counter reaches the configured limit, the algorithm returns

result.status = MixedIRStatus::diverged;

before applying the correction.

Result

With divergence detection enabled, the FP32–FP64–FP128 experiment terminates early with diverged.

With divergence detection disabled, the finite-value guards still prevent a crash. In our test, the algorithm completed 19 updates before detecting a non-finite value:

kappa = 1e+09
converged = 0
iterations = 19
final_rel_correction = 37.2908

The two safeguards have separate roles:

The divergence detector stops persistent growth early.
The finite-value checks prevent unsafe operations and crashes.

FP16 residual underflow and stagnation problem

We tested three-precision iterative refinement with

\[ T_{\mathrm{factor}}=\mathrm{FP16},\qquad T_{\mathrm{work}}=\mathrm{FP64},\qquad T_{\mathrm{residual}}=\mathrm{FP128}. \]

The LU factorization and correction solve used FP16. Iterates were stored and updated in FP64. Residuals were computed in FP128.

For \(n=100\), the convergence sweep produced unexpected results:

For \(\kappa=10\), the method did not converge after 20 iterations. The forward error stagnated near \(5.96\times10^{-8}\), and the relative correction stagnated near \(1.69\times10^{-8}\).
For \(\kappa=100\), the method appeared to converge after only three iterations because the computed correction became exactly zero.
Similar nonmonotone behavior occurred for larger condition numbers.

This was inconsistent with the usual refinement condition. Since

\[ u_{\mathrm{FP16}}=2^{-11}\approx4.88\times10^{-4}, \]

we have, for \(\kappa=10\),

\[ \kappa u_{\mathrm{FP16}} \approx4.88\times10^{-3}\ll1. \]

The \(\kappa=10\) problem should therefore have been favorable for iterative refinement.

The repeated values

\[ 5.96046\times10^{-8} \quad\text{and}\quad 1.19209\times10^{-7} \]

were an important clue. These are respectively the smallest positive FP16 subnormal number and twice that number.

Empirical diagnosis

Comparison with bfloat16

We replaced FP16 with bfloat16. Bfloat16 has fewer significand bits than FP16 but a much larger exponent range.

The relevant results were:

\(\kappa\)	Converged	Iterations	Final relative correction
1	Yes	2	\(5.07389\times10^{-17}\)
10	Yes	14	\(4.71855\times10^{-17}\)
100	No	20	\(6.68688\times10^{-15}\)
1000	No	20	\(2.24847\)
10000	No	20	\(12.5079\)

In particular, bfloat16 converged for \(\kappa=10\), while FP16 did not. Because bfloat16 has lower precision but greater range, this supported the hypothesis that the FP16 failure was caused by its restricted exponent range, not by insufficient significand precision alone.

Conceptual hypothesis

At iteration \(k\), the residual is computed accurately in FP128:

\[ r_k=b-Ax_k. \]

The original implementation then converted it directly to FP16:

\[ r_k^{(128)} \longrightarrow r_k^{(16)}. \]

For FP16,

\[ x_{\min,\mathrm{normal}} = 2^{-14} \approx6.1035\times10^{-5}, \]

and

\[ x_{\min,\mathrm{subnormal}} = 2^{-24} \approx5.9605\times10^{-8}. \]

Under round-to-nearest, values below approximately

\[ 2^{-25}\approx2.9802\times10^{-8} \]

round to zero.

We therefore hypothesized that, as refinement reduced the residual:

Some FP128 residual components became zero when converted to FP16.
Surviving components were rounded to a small number of FP16 subnormal values.
The FP16 correction solve received an increasingly inaccurate right-hand side.
The correction eventually stagnated or became exactly zero.

An exactly zero correction could also cause false convergence because the relative-correction stopping criterion would evaluate to zero.

Diagnostic implementation

We added per-iteration residual-conversion diagnostics to mixed_ir:

\(\lVert r_k\rVert_\infty\), computed in residual precision;
the smallest nonzero component of \(r_k\);
the number of nonzero FP128 residual components;
the number of nonzero components that became exactly zero during FP128-to-FP16 conversion.

The important distinction was:

nonzero_components described the original FP128 residual;
zeroed_by_conversion described information lost in the FP16 copy.

Thus,

nonzero_components = 100
zeroed_by_conversion = 20

meant that all 100 components were nonzero in FP128, but only 80 remained nonzero after conversion to FP16.

Diagnostic experiment

We tested one fixed random symmetric positive-definite system with

\[ n=100,qquad \kappa(A)=10, \]

using FP16–FP64–FP128 refinement.

The initial diagnostic results were:

Iteration	\(\lVert r_k\rVert_\infty\)	Smallest nonzero component	Zeroed by FP16 conversion
0	\(1.3474\times10^{-2}\)	\(3.3279\times10^{-5}\)	0
1	\(5.4754\times10^{-5}\)	\(3.4203\times10^{-8}\)	0
2	\(4.8336\times10^{-7}\)	\(3.0833\times10^{-9}\)	9
3	\(4.1084\times10^{-7}\)	\(2.0158\times10^{-9}\)	20

The first two refinement steps reduced the residual substantially:

\[ \frac{\lVert r_0\rVert_\infty} {\lVert r_1\rVert_\infty} \approx246, \qquad \frac{\lVert r_1\rVert_\infty} {\lVert r_2\rVert_\infty} \approx113. \]

At iteration 2, residual components first became zero during conversion. At essentially the same point, residual reduction stopped.

For \(k\ge2\),

\[ 2.97\times10^{-7} \le \lVert r_k\rVert_\infty \le 5.00\times10^{-7}, \]

with mean approximately

\[ 4.17\times10^{-7}. \]

Between 9 and 20 components were lost per iteration, with an average of approximately 14.5.

Even before components became zero, conversion was already inaccurate. At iteration 1, the smallest component was

\[ 3.4203\times10^{-8}. \]

It was above the rounding-to-zero threshold but below the smallest FP16 subnormal, so it was rounded to approximately

\[ 5.9605\times10^{-8}. \]

Therefore, zeroed_by_conversion = 0 did not imply an accurate conversion.

After iteration 2, the entire residual was only several FP16 subnormal increments away from zero. The residual was consequently represented by a very small set of possible FP16 values. This produced severe quantization even for components that remained nonzero.

The experiment confirmed the hypothesis:

The FP128 residual remained fully nonzero, but conversion to unscaled FP16 discarded or severely quantized its components. This information loss began at the same point as residual stagnation.

Implemented solution

Conceptual solution

We introduced infinity-norm residual scaling:

\[ \theta_k=\lVert r_k\rVert_\infty, \qquad \widehat r_k=\frac{r_k}{\theta_k}. \]

Instead of solving directly with the small residual, we solve

\[ A\widehat d_k=\widehat r_k \]

in factor precision and then recover the actual correction:

\[ d_k=\theta_k\widehat d_k. \]

In exact arithmetic, this is equivalent to

\[ Ad_k=r_k. \]

The normalization ensures

\[ \lVert\widehat r_k\rVert_\infty=1. \]

The residual conversion therefore uses the available FP16 range effectively rather than placing the entire vector near zero.

Scaling does not increase FP16 precision. It removes the irrelevant absolute magnitude of the residual before conversion.

Implementation in `mixed_ir`

We added an optional setting:

bool scale_residual = false;

The default preserves the original algorithm.

When scaling is enabled, the implementation:

Computes \(\lVert r_k\rVert_\infty\) in T_residual.
Preserves the original residual r_r for diagnostics.

Creates a separate correction right-hand side:

hdnum::Vector<T_residual> correction_rhs_r(r_r);

Normalizes it when the residual is nonzero:

correction_rhs_r[i] =
    correction_rhs_r[i] / residual_scale_r;

Converts the normalized vector, rather than the original residual, to factor precision:
```
convert(r_f, correction_rhs_r);
```
Performs the existing FP16 LU solve without changing the factors or solve routine.
Converts the normalized correction to working precision.
Rescales it in FP64:
```
d_w[i] = residual_scale_w * d_w[i];
```
Performs the original update:

\[ x_{k+1}=x_k+d_k. \]

The scale is applied in working precision, not factor precision. Thus, the small value \(\theta_k\) never needs to be represented in FP16.

The zero-residual case skips normalization, preventing division by zero. We also retained the existing finite-value checks after rescaling.

The diagnostic semantics were updated accordingly:

residual magnitudes still describe the original \(r_k\);
zeroed_by_conversion describes conversion of the actual correction right-hand side:
- \(r_k\) in unscaled mode;
- \(\widehat r_k\) in scaled mode.

Validation experiment

Experimental design

We generated the problem once and ran:

run("unscaled", false);
run("scaled", true);

Both runs used exactly the same:

matrix;
right-hand side;
random seed;
dimension \(n=100\);
condition number \(\kappa(A)=10\);
precision triple FP16–FP64–FP128;
stopping criterion and iteration limit.

The identical initial residual,

\[ \lVert r_0\rVert_\infty = 1.34744462924774981\times10^{-2}, \]

confirmed that scaling did not alter the matrix factorization or initial solution.

The unscaled results also reproduced the previous baseline. Therefore, adding the optional scaling path did not change the original behavior when scaling was disabled.

Unscaled result

The unscaled method:

first lost components at iteration 2;
lost between 9 and 20 components per later iteration;
stagnated near \(4.17\times10^{-7}\);
failed to converge after 20 updates;
ended with

\[ \text{final relative correction} = 1.1685\times10^{-7}. \]

Scaled result

The scaled run produced:

Iteration	\(\lVert r_k\rVert_\infty\)	Reduction factor
0	\(1.3474\times10^{-2}\)	–
1	\(5.3874\times10^{-5}\)	\(250.1\)
2	\(2.8349\times10^{-7}\)	\(190.0\)
3	\(2.1471\times10^{-9}\)	\(132.0\)
4	\(1.1636\times10^{-11}\)	\(184.5\)
5	\(7.5515\times10^{-14}\)	\(154.1\)
6	\(4.3479\times10^{-16}\)	\(173.7\)

No normalized residual component became zero during FP16 conversion.

The residual decreased by

\[ \frac{\lVert r_0\rVert_\infty} {\lVert r_6\rVert_\infty} \approx3.10\times10^{13}, \]

or approximately 13.5 decimal orders.

The method converged after seven updates with

\[ \text{final relative correction} = 9.7687\times10^{-17}. \]

This was below FP64 unit roundoff:

\[ u_{\mathrm{FP64}} = 2^{-53} \approx1.1102\times10^{-16}. \]

Interpretation

The observed reduction factors, between approximately 132 and 250, correspond to contraction factors between roughly

\[ 4.0\times10^{-3} \quad\text{and}\quad 7.6\times10^{-3}. \]

This agrees well with the theoretical scale

\[ \kappa(A)u_{\mathrm{FP16}} \approx4.88\times10^{-3}. \]

Once the residual-range problem was removed, the algorithm behaved as expected from mixed-precision iterative-refinement theory.

Conclusion

The investigation established the following sequence:

FP16 refinement failed unexpectedly for a well-conditioned problem.
Bfloat16 succeeded despite having lower significand precision, pointing to an exponent-range problem.
Direct diagnostics showed that nonzero FP128 residual components were zeroed or severely quantized during FP16 conversion.
Conversion loss began at the same point as residual stagnation.
Infinity-norm scaling removed the irrelevant absolute residual magnitude before conversion.
A controlled scaled-versus-unscaled experiment eliminated all observed conversion-to-zero events.
The scaled method restored regular convergence and reached FP64 accuracy in seven iterations.

This is strong validation for the tested system. Broader validation should repeat the comparison over multiple dimensions, condition numbers, and random seeds. Residual scaling addresses residual underflow only; it does not address separate FP16 overflow caused by matrix entries exceeding the FP16 range.

5.7 Plotting Convergence Experiments

The script

code/scripts/plot_convergence.py

generates convergence-history and summary plots from one mixed-precision iterative-refinement CSV file.

It is designed to be reusable for different:

matrix generators;
right-hand-side or solution-generation modes;
precision combinations.

The input CSV is selected using a command-line argument, so a separate plotting script is not needed for each experiment.

Expected CSV format

The CSV file must contain the columns:

kappa
iteration
forward_error_inf
backward_error_inf
rel_correction
converged
total_iterations
final_rel_correction

Each row describes one iterate for one condition number.

The usual relationship is:

iteration 0  → initial solution x0
iteration 1  → solution x1 after the first correction
iteration 2  → solution x2 after the second correction
...

Generated plots

For every input CSV, the script generates six plots.

Convergence histories

Relative forward error against refinement iteration.
Normwise backward error against refinement iteration.
Relative correction norm against refinement iteration.

Each condition number is displayed as a separate curve.

Summary plots

Final forward and backward errors against condition number.
Best forward and backward errors attained at any iteration.
Total refinement iterations against condition number.

The best-attained-error plot is useful when an iteration initially improves and later deteriorates. For example, a difficult () problem may attain its best forward error around iteration 7 but have a worse result at iteration 20.

The script also prints a numerical summary containing:

condition number;
convergence flag;
total iterations;
final forward error;
final backward error;
final relative correction;
best forward and backward errors;
iterations at which the best errors occurred.

Logarithmic plots and zero values

The error histories use logarithmic vertical axes. Since zero cannot be displayed on a logarithmic scale, the script replaces nonpositive plotting values with NaN.

This affects only plotting. The original CSV data is not modified.

Consequently, a curve or point with an exactly zero error may disappear from a logarithmic plot.

Command-line interface

The general invocation is:

python code/scripts/plot_convergence.py CSV_FILE \
    --label "Description of experiment"

If the terminal is already inside code/, use:

python scripts/plot_convergence.py CSV_FILE \
    --label "Description of experiment"

The script first checks whether CSV_FILE is a direct path. If it is not found there, it looks under:

code/results/raw/

The plots are written to:

code/results/plots/

Default behavior

The CSV argument is optional. If no filename is provided, the script uses:

convergence_fp16_fp64_fp128.csv

Therefore:

python code/scripts/plot_convergence.py

is equivalent to:

python code/scripts/plot_convergence.py \
    convergence_fp16_fp64_fp128.csv

Output tags

The output tag is normally derived automatically from the CSV filename.

For example, from:

convergence_random_spd_fp16_fp64_fp128.csv

the script removes the initial convergence_ and obtains:

random_spd_fp16_fp64_fp128

It then generates filenames such as:

convergence_forward_error_random_spd_fp16_fp64_fp128.png
convergence_backward_error_random_spd_fp16_fp64_fp128.png
convergence_rel_correction_random_spd_fp16_fp64_fp128.png
summary_final_errors_random_spd_fp16_fp64_fp128.png
summary_best_errors_random_spd_fp16_fp64_fp128.png
summary_iterations_random_spd_fp16_fp64_fp128.png

A custom output tag can be supplied with:

--tag custom_name

For example:

python code/scripts/plot_convergence.py \
    convergence_random_spd_fp16_fp64_fp128.csv \
    --tag random_spd_test

Descriptive plot labels

The optional --label argument is added to plot titles:

--label "Random dense SPD"

This produces titles such as:

Random dense SPD: relative forward error, infinity norm vs. refinement iteration

The label does not affect the output filenames.

Matrix-family examples

Block-rotated SPD

python code/scripts/plot_convergence.py \
    convergence_fp16_fp64_fp128.csv \
    --label "Block-rotated SPD"

Random dense SPD

python code/scripts/plot_convergence.py \
    convergence_random_spd_fp16_fp64_fp128.csv \
    --label "Random dense SPD"

Random SVD

python code/scripts/plot_convergence.py \
    convergence_random_svd_fp16_fp64_fp128.csv \
    --label "Random SVD"

Encoding the solution or RHS mode

The test-problem infrastructure supports three modes:

RightHandSideMode::ones_solution
RightHandSideMode::random_sign_solution
RightHandSideMode::random_normal_rhs

These correspond to:

construct (x=(1,,1)^T) and form (b=Ax);
construct random (x_i) and form (b=Ax);
generate (b_iN(0,1)) and compute a high-precision reference solution.

The third mode uses a normal distribution, not a uniform distribution.

It is sensible to encode this choice in the CSV filename. A consistent naming scheme is:

convergence_<matrix-family>_<rhs-mode>_<factor>_<work>_<residual>.csv

For example:

convergence_random_spd_x_ones_fp16_fp64_fp128.csv
convergence_random_spd_x_random_signs_fp16_fp64_fp128.csv
convergence_random_spd_b_random_normal_fp16_fp64_fp128.csv

The plotting script automatically retains the RHS mode as part of the output tag. No code change is needed.

Examples for the three modes

All-ones solution

python code/scripts/plot_convergence.py \
    convergence_random_spd_x_ones_fp16_fp64_fp128.csv \
    --label "Random SPD, ones solution"

This produces files such as:

convergence_forward_error_random_spd_x_ones_fp16_fp64_fp128.png

Random-sign solution

python code/scripts/plot_convergence.py \
    convergence_random_spd_x_random_signs_fp16_fp64_fp128.csv \
    --label "Random SPD, random-sign solution"

This produces files such as:

convergence_forward_error_random_spd_x_random_signs_fp16_fp64_fp128.png

Random-normal right-hand side

python code/scripts/plot_convergence.py \
    convergence_random_spd_b_random_normal_fp16_fp64_fp128.csv \
    --label "Random SPD, normally distributed right-hand side"

This produces files such as:

convergence_forward_error_random_spd_b_random_normal_fp16_fp64_fp128.png

What the script detects automatically

The script automatically derives an output tag from the filename. It does not semantically interpret components such as:

random_spd
x_ones
x_random_signs
b_random_normal
fp16_fp64_fp128

Therefore:

output filenames are generated correctly automatically;
the same plots are generated for every mode;
a human-readable title should currently be supplied through --label.

Automatic semantic interpretation could be added later, but it is not necessary for the current workflow.

Recording metadata inside the CSV

Encoding the experiment configuration in the filename is useful, but the CSV should eventually also record metadata such as:

matrix_family
rhs_mode
matrix_seed_u
matrix_seed_v
vector_seed

For example:

matrix_family,rhs_mode,matrix_seed_u,matrix_seed_v,vector_seed,kappa,...
random_spd,x_random_signs,42,137,2718,10,...

This makes the experiment self-describing even if the file is renamed.

Later comparison plots

plot_convergence.py plots one CSV dataset at a time. It should remain focused on that task.

A later script, for example

code/scripts/plot_compare_generators.py

code/scripts/plot_compare_rhs_modes.py

can place multiple experiments on the same figures.

Possible comparisons include:

block-rotated SPD versus random SPD versus random SVD;
all-ones versus random-sign solutions;
prescribed (x) versus random-normal (b);
different precision triples;
different random seeds.

Thus, the intended organization is:

plot_convergence.py
    → detailed plots for one experiment

plot_compare_generators.py
    → compare matrix families

plot_compare_rhs_modes.py
    → compare solution/RHS-generation strategies

5.1 Temporary Matrix Generator

5.2 Matrix Generator Supervisor Version

Nonsymmetric Generator

5.3 Configurable test-problem Construction

Three Right-hand-side Modes

Test-problem Options

Linear-system Representation

Shared complete_problem function

High-precision reference solve

Matrix-generator interfaces

Usage Examples

5.4 Error Metrics Helper

5.5 First Convergence Test

5.6 Robustness Improvements to mixed_ir

1. FP16–FP64–FP128: overflow during input conversion

Problem

Detection

Solution

2. FP32–FP64–FP128: divergence of iterative refinement

Problem

Detection

Solution

Source-Code Changes

Termination statuses

Finite-value checks

Guard placement

Divergence options

Correction-norm history

Detector state and comparison

Result

FP16 residual underflow and stagnation problem

Empirical diagnosis

Comparison with bfloat16

Conceptual hypothesis

Diagnostic implementation

Diagnostic experiment

Implemented solution

Conceptual solution

Implementation in mixed_ir

Validation experiment

Experimental design

Unscaled result

Scaled result

Interpretation

Conclusion

5.7 Plotting Convergence Experiments

Expected CSV format

Generated plots

Convergence histories

Summary plots

Logarithmic plots and zero values

Command-line interface

Default behavior

Output tags

Descriptive plot labels

Matrix-family examples

Block-rotated SPD

Random dense SPD

Random SVD

Encoding the solution or RHS mode

Examples for the three modes

All-ones solution

Random-sign solution

Random-normal right-hand side

What the script detects automatically

Recording metadata inside the CSV

Later comparison plots

Shared `complete_problem` function

5.6 Robustness Improvements to `mixed_ir`

Implementation in `mixed_ir`