Numerically Unstable

Numerical instability refers to algorithms that are prone to large errors in their output due to accumulated rounding errors and finite precision arithmetic. A small change in input can lead to a large change in output.

Some common unstable numerical algorithms:

Naive summation of partial sums
Gaussian elimination without pivoting
Gram-Schmidt process for orthonormalization

Java example - Naive summation:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
// Unstable - sum in ascending order
double sum(double[] arr) {
  double total = 0;
  for (int i = 0; i < arr.length; i++) {
    total += arr[i]; 
  }
  return total;
}

// Stable - sum in descending order
double sum(double[] arr) {
  double total = 0;
  for (int i = arr.length-1; i >= 0; i--) {
    total += arr[i];
  }
  return total; 
}

C++ example - Gaussian elimination:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
// Unstable - no pivoting
void gaussianElimination(vector<vector<double>> A) {
  int n = A.size();

  for (int i = 0; i < n; i++) {
    // Elimination step 
  }

}

// Stable - with partial pivoting
void gaussianElimination(vector<vector<double>>& A) {
  int n = A.size();
  
  for (int i = 0; i < n; i++) {
    // Pivot step
  }

  // Elimination
}

Python example - Polynomial evaluation:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
# Unstable 
def evaluate(coeffs, x):
  result = coeffs[0]
  for i in range(1, len(coeffs)):
    result += coeffs[i] * x**i
  return result

# Stable - Horner's method
def evaluate(coeffs, x):
  result = coeffs[-1]
  for i in range(len(coeffs)-2, -1, -1):
    result = coeffs[i] + x * result
  return result

Numerical instability leads to incorrect results. Awareness and pivoting, normalization etc. help.

An algorithm is considered numerically unstable if small errors or changes in the input or intermediate calculations can lead to large errors in the output. Such algorithms can produce inaccurate or unreliable results, especially when dealing with floating-point arithmetic or iterative methods. Numerical instability can manifest in various ways, such as loss of precision, overflow, or even wildly incorrect results.

For example, let’s consider summing up an array of floating-point numbers. One might think that the order of addition doesn’t matter, but it can in the context of numerical instability.

Example Code

Java

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
public class NumericallyUnstableExample {
    public static void main(String[] args) {
        double[] numbers = {1e10, -1e10, 1, 1};
        double sum = 0;
        for(double num : numbers) {
            sum += num;
        }
        System.out.println("Sum: " + sum);  // Ideally should be 2
    }
}

C++

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
#include <iostream>
#include <vector>
using namespace std;

int main() {
    vector<double> numbers = {1e10, -1e10, 1, 1};
    double sum = 0;
    for(auto num : numbers) {
        sum += num;
    }
    cout << "Sum: " << sum << endl;  // Ideally should be 2
    return 0;
}

Python

1
2
3
numbers = [1e10, -1e10, 1, 1]
sum_result = sum(numbers)
print(f"Sum: {sum_result}")  # Ideally should be 2

In these examples, we are trying to sum up the numbers [1e10, -1e10, 1, 1]. Mathematically, the sum should be 2. However, because of the large magnitudes involved and floating-point arithmetic, the sum might not be exactly 2. This example demonstrates how even simple calculations can become numerically unstable.