Crafting a Numerical Instability Dataset

1. The Strategy: Creating Instability by Shifting Data

From our first module, you'll recall that the naive single-pass algorithm is vulnerable to catastrophic cancellation when calculating variance. This happens when the two main terms in its formula are very large and nearly equal:

$s^2=\frac{\stackrel{\text{Large number 1}}{\overbrace{\sum _{i=1}^n{x}_i^2}}-\stackrel{\text{Large number 2}}{\overbrace{\frac{({\sum }_{i=1}^n{x}_i{)}^2}{n}}}}{n-1}$

The problem is most severe when the data's standard deviation is much smaller than its mean. So, how do we create a dataset with this property on demand?

The strategy is surprisingly simple:

1. Start with a collection of numbers that have a small mean and a small variance.
2. Add a very large constant value (a "shift") to every number in the collection.

Let's consider the effect of this shift. If we take a dataset $X=\{x_1,x_2,...,x_n\}$ and create a new dataset $Y=\{x_1+C,x_2+C,...,x_n+C\}$ by adding a large constant $C$, the statistical properties change in a predictable way:

• The Mean: The new mean becomes ${\mu }_Y={\mu }_X+C$. It increases by the value of the shift.
• The Variance: The variance remains unchanged, ${\sigma }_Y^2={\sigma }_X^2$, because variance measures the spread of the data, which is unaffected by shifting all points together.

By choosing a large $C$, we create a dataset with a very large mean but the original small variance—exactly the conditions required to trigger catastrophic cancellation.

2. A Proven Recipe for Instability

This data-shifting technique is a standard method for testing the numerical stability of statistical algorithms. Researchers use it to show how different implementations behave under stress.

The paper "A Closer Look at Variance Implementations In Modern Database Systems" describes this exact method for their experiments. Let's see how they did it.

As the paper describes, the recipe is:

1. Generate base data: Create a set of data points from a Uniform(0,1) distribution. Data from this distribution has a theoretical mean of 0.5 and a theoretical variance of $1/12\approx 0.0833$, which is quite small.
2. Apply a shift: Add a large number, like $10^9$ or $10^{12}$, to each data point.

The resulting dataset has a mean very close to the shift value (e.g., $10^9$) but a variance that should still be very close to $1/12$. This is the "poison pill" for the naive algorithm.

Another paper, "Scalable and Numerically Stable Descriptive Statistics," uses a similar technique, generating data in ranges like [1000000.0, 1000000.5). This is equivalent to starting with data from Uniform(0, 0.5) and adding a shift of $10^6$, confirming that this is a standard approach.

3. Implementing the Dataset Generator

Now, let's turn this recipe into a reusable Python function. Since you have a background in data science, you'll be familiar with NumPy, which is perfect for this task.

Our function will take the number of data points (n) and the desired shift as inputs and return the unstable dataset.

import numpy as np

def create_unstable_dataset(n, shift, seed=42):
    """
    Generates a dataset designed to cause numerical instability in the
    naive single-pass variance algorithm.

    The dataset consists of numbers from a Uniform(0,1) distribution
    with a large constant 'shift' added to each.

    Args:
        n (int): The number of data points in the dataset.
        shift (float): A large constant value to add to each data point.
        seed (int): A seed for the random number generator for reproducibility.

    Returns:
        numpy.ndarray: The generated dataset with mean ~shift and variance ~1/12.
    """
    # Set the seed for the random number generator to ensure our results
    # are the same every time we run the code.
    np.random.seed(seed)
    
    # 1. Generate n data points from a Uniform(0,1) distribution.
    base_data = np.random.rand(n)
    
    # 2. Add the large shift to each data point.
    shifted_data = base_data + shift
    
    return shifted_data

# Example of how to use it:
# This creates 10,000 data points with a mean around 1 billion
# and a variance that should be around 1/12.
unstable_data = create_unstable_dataset(n=10000, shift=1e9)
print(f"First 5 values: {unstable_data[:5]}")
print(f"Dataset mean: {unstable_data.mean()}")
print(f"Dataset variance (using stable numpy.var): {unstable_data.var(ddof=1)}")

Notice in the example output that the mean is enormous, while the variance is tiny. This is the numerical trap we've set.

4. A Deeper Look at the Failure Mechanism

Why does this shifted data cause the naive algorithm to fail so spectacularly? The paper we looked at earlier provides an excellent, concrete illustration.

As the paper explains, when the shift is large, the two numbers being subtracted, $\sum x_i^2$ and $\frac{(\sum x_i{)}^2}{n}$, become dominated by the shift value. They grow to be massive and extremely close to each other.

A standard double-precision floating-point number has about 15-17 decimal digits of precision. If the two numbers you are subtracting agree to 16 digits, but the actual difference is in the 17th digit and beyond, that information is lost. The result of the subtraction might be zero or just numerical "noise". This is exactly what Table 3 in the paper shows: as the shift exponent increases, more of the significant bits of the two terms become identical, until the subtraction wipes out all the meaningful information about the variance.

This is why the two-pass and Welford's algorithms are superior. They don't subtract two massive, nearly-identical numbers. Instead, they work with the differences from the mean, $(x_i-\mu )$, which are small numbers, preserving precision.

Conclusion

In this lesson, we have successfully created the primary tool for our final experiment. You now understand not only how to construct a dataset that breaks the naive variance algorithm but also why it works at the level of floating-point arithmetic.

Key Takeaways:

• You can create a numerically unstable dataset for variance calculation by adding a large constant (a shift) to a low-variance set of numbers.
• This technique creates a dataset with a large mean but a small variance, the ideal condition to trigger catastrophic cancellation.
• We've implemented a function, create_unstable_dataset, that puts this theory into practice, giving us a way to generate these datasets on demand.

Preview of the Next Lesson:

The stage is now set for the final act. In our next and final lesson, we will bring everything together. We will use the create_unstable_dataset function to generate a challenging dataset, run all three of our algorithms (Welford, naive, and two-pass) on it, and compare their accuracy. This will provide a clear, empirical demonstration of the power and importance of numerically stable algorithms like Welford's.