Floating points: Make Gaussian elimination go wrong You know that some decimal numbers can't be expressed as IEEE 754 floating points. You also know that arithmetic with floating points can give results that seem to be wrong:


System.out.print(4.0-3.1); (source)
0.2 + 0.04 = 0.24000000000000002 (source)


While one wrong number with an error that is this small might not be significant, more of those error might be.

Task

Find a system of linear equations in form of a matrix A \in R^{n \times n} and b \in R^n, such that the result vector is as wrong as possible. Reference is the script below, executed with Python 2.7.

Points

#!/usr/bin/env python
# -*- coding: utf-8 -*-

def pprint(A):
    n = len(A)
    for i in range(0, n):
        line = ""
        for j in range(0, n+1):
            line += str(A[i][j]) + "\t"
            if j == n-1:
                line +=  "| "
        print(line)
    print("")

def gauss(A):
    n = len(A)

    for i in range(0,n):
        # Search for maximum in this column
        maxEl = abs(A[i][i])
        maxRow = i
        for k in range(i+1,n):
            if abs(A[k][i]) &gt; maxEl:
                maxEl = A[k][i]
                maxRow = k

        # Swap maximum row with current row (column by column)
        for k in range(i,n+1):
            tmp = A[maxRow][k]
            A[maxRow][k] = A[i][k]
            A[i][k] = tmp

        # Make all rows below this one 0 in current column
        for k in range(i+1,n):
            c = -A[k][i]/A[i][i]
            for j in range(i,n+1):
                if i==j:
                    A[k][j] = 0
                else:
                    A[k][j] += c * A[i][j]

    # Solve equation Ax=b for an upper triangular matrix A
    x=[0 for i in range(n)]
    for i in range(n-1,-1,-1):
        x[i] = A[i][n]/A[i][i]
        for k in range(i-1,-1,-1):
            A[k][n] -= A[k][i] * x[i]
    return x;

if __name__ == "__main__":
    from fractions import Fraction
    datatype = float # Fraction 
    n = input()

    A = [[0 for j in range(n+1)] for i in range(n)]

    # Read input data
    for i in range(0,n):
        line = map(datatype, raw_input().split(" "))
        for j, el in enumerate(line):
            A[i][j] = el
    raw_input()

    line = raw_input().split(" ")
    lastLine = map(datatype, line)
    for i in range(0,n):
        A[i][n] = lastLine[i]

    # Print input
    pprint(A)

    # Calculate solution
    x = gauss(A)

    # Print result
    line = "Result:\t"
    for i in range(0,n):
       line += str(x[i]) + "\t"
    print(line)


Let x be the correct answer and x' the answer that the script above gives.

Your points are (||x-x'||^2)/n, where || . || is the euclidean distance.

Your solution gets 0 points if it gives nan as result.

Example

bad.in:

2
1 4
2 3.1

1 1


The first line is n, the next n lines describe the matrix A and the last line describes b.

$ python gauss.py &lt; bad.in 
1.0 4.0 | 1.0   
2.0 3.1 | 1.0   

Result: 0.183673469388  0.204081632653  


The answer should be:

Result: 9/49    10/49


So the score would be:

((0.183673469388-9/49)^2+(0.204081632653-10/49)^2)/2
= 3.186172428154... × 10^-26


The last step was calculated with Wolfram|Alpha.
Taking this as a question about numerical inaccuracy in calculations rather than in conversion from numbers which can't be exactly represented in 64 bits to 64-bit representations, we're interested in nearly singular matrices.

Let's take the simplest possible case:

[a b] [x] = [p]
[c d] [y]   [q]


has solution

[x] = [(dp - bq) / (ad - bc)]
[y]   [(aq - cp) / (ad - bc)]


so that if ad - bc is very small we find that errors in dp - bq and aq - cp are magnified.

We also find that scaling p and q scales the absolute error linearly, which points out a massive flaw in the scoring system of this question.

How can we get a nearly singular matrix? In order of ill-conditionedness:


a=1+e, b=c=d=1 gives ad - bc = e. This allows us to get a determinant of 1 ulp (which for IEEE 64-bit is on the order of 10^-16).
a=1+e, d=1-e, b=c=1 gives ad - bc = -e^2. This allows us to get a determinant on the order of 10^-32.

2
1.000000000000001 1
1 0.999999999999999

1 1


gives output:

1.0 1.0 | 1.0   
1.0 1.0 | 1.0   

Result: -9.0    10.0


rather than x = 100000000000000, y = -100000000000000

Similarly,

2
1.000000000000001 1
1 0.999999999999999

1E307 1E307


gives

1.0 1.0 | 1e+307    
1.0 1.0 | 1e+307    

Result: -9.1120238836e+307  1.01120238836e+308  


rather than x = 1E321, y = -1E321. This seems to roundly beat Johannes' solution (it gives a score on the order of 10^642 vs 10^588), and moreover it would continue to generate about the same score even if the input numbers were tweaked to be exactly representable in IEEE 64-bit, whereas his would then score 0.
a=d=e, bc=0 gives ad - bc = e^2 but this time we can choose e to be a very small number rather than 1 ULP. In principle this allows us to get a determinant on the order of 10^-600 (i.e. underflowing IEEE 64-bit). However, so far I haven't managed to get a decent score without an infinity in the answer, and I gather from your chat with Johannes that even though the question doesn't rule out infinite scores, in practice you do.