# 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.

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]) > 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)]

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

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 < 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. • Hooray! Now there's a geeky comic we can link to whenever this question comes up! – mob Commented Jun 8, 2013 at 2:13 ## 2 Answers 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: 1. 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). 2. 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. 3. 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. • I've wanted this challenge to be a numerical one, so you've answered what I wanted to (but failed to) ask. Thanks! Commented Jun 5, 2013 at 13:22 • +1 Good work. I just used a flaw in the question, but you actually solved it. Commented Jun 5, 2013 at 13:30 Ok, this feels kind of cheating. # Tcl, score 1234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234569876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876543209876544 gets stdin a puts$a
for {set i 0} {$i<$a} {incr i} {
set r {}
for {set j 0} {$j<$a} {incr j} {
lappend r [expr {$i==$j?1:0}]
}
puts $r } puts "" puts [join [lrepeat$a [string repeat 8 308]]]

• Commented May 29, 2013 at 19:52
• Ok, seems to be correct: pastebin.com/JgKrhLe1 I didn't think it was so easy to get such a high number of points. Congratulations. I will wait until tomorrow, but I think you've won this challenge. Commented May 29, 2013 at 20:28
• Actually it is 1.2345*10^588 Commented May 29, 2013 at 20:42
• It is kind of cheating. You're not exploiting any instability of Gaussian elimination at all, but merely the fact that Python will parse arbitrary bigints and convert them to doubles. Commented Jun 5, 2013 at 10:52
• @PeterTaylor Yeah, it feels that way. This question asks for the biggest difference, not a relative difference. So using high numbers almost win the challenge. Relative error would be better for this challenge. Commented Jun 5, 2013 at 13:02