## When solving Ax=b on a computer, two things can go wrong.¶

• The first is that the $Ax=b$ problem we have may be ill-conditioned.
• The second is that we may be using a bad algorithms.

### Ill-conditioned problems¶

If Ax = b is ill-conditioned, this means that, if we make any tiny error to b, even one of order $10^{-16}$, then we'll get a totally different answer to our linear system. The only way to solve ill-conditioned linear systems accurately on a computer is to use additional precision in our arithmetic.

In [1]:
function hilb(n;T=Float64)
A = zeros(T,n,n)
for j=1:n
for i=1:n
A[i,j] = 1./(i+j)
end
end
return A
end

Out[1]:
hilb (generic function with 1 method)
In [2]:
hilb(5)

Out[2]:
5x5 Array{Float64,2}:
0.5       0.333333  0.25      0.2       0.166667
0.333333  0.25      0.2       0.166667  0.142857
0.25      0.2       0.166667  0.142857  0.125
0.2       0.166667  0.142857  0.125     0.111111
0.166667  0.142857  0.125     0.111111  0.1     
In [13]:
## Small system
n = 5
A = hilb(n)
xtrue = ones(n)
b = A*xtrue
xcomp = A\b

@printf("\n  n = %i\n\n", n);
@printf("    [ %25s %25s ]\n", "x true", "x computed");
for i=1:n
@printf("    [ %25.17f %25.17f ]\n", xtrue[i], xcomp[i]);
end

  n = 5

[                    x true                x computed ]
[       1.00000000000000000       0.99999999999925315 ]
[       1.00000000000000000       1.00000000000871170 ]
[       1.00000000000000000       0.99999999997022349 ]
[       1.00000000000000000       1.00000000003894085 ]
[       1.00000000000000000       0.99999999998275368 ]

In [7]:
## Bigger system
n = 12
A = hilb(n)
xtrue = ones(n)
b = A*xtrue
xcomp = A\b
@printf("\n  n = %i\n\n", n);
@printf("    [ %25s %25s ]\n", "x true", "x computed");
for i=1:n
@printf("    [ %25.17f %25.17f ]\n", xtrue[i], xcomp[i]);
end

  n = 12

[                    x true                x computed ]
[       1.00000000000000000       1.00000034831994822 ]
[       1.00000000000000000       0.99997498970416920 ]
[       1.00000000000000000       1.00058545274899036 ]
[       1.00000000000000000       0.99339772886215538 ]
[       1.00000000000000000       1.04233097346954851 ]
[       1.00000000000000000       0.83186032680000732 ]
[       1.00000000000000000       1.43279408301936373 ]
[       1.00000000000000000       0.26524866529664204 ]
[       1.00000000000000000       1.81683666353393236 ]
[       1.00000000000000000       0.42798869683200724 ]
[       1.00000000000000000       1.22886800007404573 ]
[       1.00000000000000000       0.96011404213052920 ]


The hilbert matrix is very badly conditioned. There is no way to accurately solve a 12-by-12 system in double precision floating point. The solution does not look wrong in floating point.

In [4]:
cond(hilb(12))

Out[4]:
6.593566388262902e16
In [8]:
hilb(12)*xcomp - b

Out[8]:
12-element Array{Float64,1}:
0.0
4.44089e-16
0.0
0.0
0.0
2.22045e-16
1.11022e-16
1.11022e-16
1.11022e-16
-1.11022e-16
0.0
1.11022e-16

In [9]:
# define a matrix with a bad growth factor
A = eye(n)
A = A - tril(ones(n,n),-1)
A[:,end] = 1
return A
end

Out[9]:
bad_growth (generic function with 1 method)
In [10]:
bad_growth(5)

Out[10]:
5x5 Array{Float64,2}:
1.0   0.0   0.0   0.0  1.0
-1.0   1.0   0.0   0.0  1.0
-1.0  -1.0   1.0   0.0  1.0
-1.0  -1.0  -1.0   1.0  1.0
-1.0  -1.0  -1.0  -1.0  1.0
In [11]:
n = 60
xtrue = (2*(rand(n,1) .< 0.5)-1)/10
b = A*xtrue
xcomp = A\b

@printf("\n  n = %i\n\n", n);
@printf("    [ %25s %25s ]\n", "x true", "x computed");
for i=1:n
@printf("    [ %25.17f %25.17f ]\n", xtrue[i], xcomp[i]);
end

  n = 60

[                    x true                x computed ]
[       0.10000000000000001       0.10000000000000001 ]
[       0.10000000000000001       0.10000000000000003 ]
[       0.10000000000000001       0.09999999999999998 ]
[      -0.10000000000000001      -0.10000000000000009 ]
[       0.10000000000000001       0.09999999999999987 ]
[       0.10000000000000001       0.09999999999999964 ]
[       0.10000000000000001       0.09999999999999964 ]
[       0.10000000000000001       0.09999999999999787 ]
[       0.10000000000000001       0.09999999999999787 ]
[       0.10000000000000001       0.09999999999999432 ]
[      -0.10000000000000001      -0.10000000000000853 ]
[       0.10000000000000001       0.09999999999996589 ]
[       0.10000000000000001       0.09999999999996589 ]
[       0.10000000000000001       0.09999999999990905 ]
[      -0.10000000000000001      -0.10000000000013642 ]
[       0.10000000000000001       0.09999999999945430 ]
[      -0.10000000000000001      -0.10000000000036380 ]
[       0.10000000000000001       0.09999999999854481 ]
[      -0.10000000000000001      -0.10000000000218279 ]
[      -0.10000000000000001      -0.10000000000582077 ]
[      -0.10000000000000001      -0.10000000000582077 ]
[       0.10000000000000001       0.09999999997671694 ]
[       0.10000000000000001       0.09999999997671694 ]
[      -0.10000000000000001      -0.10000000009313226 ]
[      -0.10000000000000001      -0.10000000009313226 ]
[       0.10000000000000001       0.09999999962747097 ]
[      -0.10000000000000001      -0.10000000055879354 ]
[      -0.10000000000000001      -0.10000000149011612 ]
[       0.10000000000000001       0.09999999776482582 ]
[      -0.10000000000000001      -0.10000000894069672 ]
[       0.10000000000000001       0.09999999403953552 ]
[      -0.10000000000000001      -0.10000002384185791 ]
[       0.10000000000000001       0.09999996423721313 ]
[      -0.10000000000000001      -0.10000014305114746 ]
[      -0.10000000000000001      -0.10000014305114746 ]
[       0.10000000000000001       0.09999942779541016 ]
[       0.10000000000000001       0.09999942779541016 ]
[      -0.10000000000000001      -0.10000228881835938 ]
[      -0.10000000000000001      -0.10000228881835938 ]
[      -0.10000000000000001      -0.10000610351562500 ]
[       0.10000000000000001       0.09999084472656250 ]
[      -0.10000000000000001      -0.10003662109375000 ]
[      -0.10000000000000001      -0.10003662109375000 ]
[       0.10000000000000001       0.09985351562500000 ]
[      -0.10000000000000001      -0.10009765625000000 ]
[       0.10000000000000001       0.09960937500000000 ]
[      -0.10000000000000001      -0.10058593750000000 ]
[      -0.10000000000000001      -0.10156250000000000 ]
[      -0.10000000000000001      -0.10156250000000000 ]
[       0.10000000000000001       0.09375000000000000 ]
[       0.10000000000000001       0.09375000000000000 ]
[       0.10000000000000001       0.06250000000000000 ]
[      -0.10000000000000001      -0.12500000000000000 ]
[       0.10000000000000001       0.00000000000000000 ]
[       0.10000000000000001       0.00000000000000000 ]
[      -0.10000000000000001      -0.50000000000000000 ]
[       0.10000000000000001       0.00000000000000000 ]
[      -0.10000000000000001       0.00000000000000000 ]
[      -0.10000000000000001       0.00000000000000000 ]
[       0.10000000000000001       0.10000000000000001 ]

In [33]:
## Using a different algorithms helps
Q,R = qr(A)
xcomp = R\(Q'*b)

@printf("\n  n = %i\n\n", n);
@printf("    [ %25s %25s ]\n", "x true", "x computed");
for i=1:n
@printf("    [ %25.17f %25.17f ]\n", xtrue[i], xcomp[i]);
end

  n = 60

[                    x true                x computed ]
[      -0.10000000000000001      -0.09999999999999772 ]
[       0.10000000000000001       0.09999999999999992 ]
[       0.10000000000000001       0.09999999999999952 ]
[      -0.10000000000000001      -0.10000000000000062 ]
[      -0.10000000000000001      -0.10000000000000013 ]
[      -0.10000000000000001      -0.09999999999999995 ]
[      -0.10000000000000001      -0.10000000000000032 ]
[      -0.10000000000000001      -0.09999999999999977 ]
[       0.10000000000000001       0.09999999999999984 ]
[      -0.10000000000000001      -0.10000000000000006 ]
[      -0.10000000000000001      -0.09999999999999973 ]
[      -0.10000000000000001      -0.09999999999999983 ]
[       0.10000000000000001       0.09999999999999992 ]
[      -0.10000000000000001      -0.10000000000000005 ]
[       0.10000000000000001       0.09999999999999987 ]
[       0.10000000000000001       0.09999999999999996 ]
[      -0.10000000000000001      -0.09999999999999973 ]
[      -0.10000000000000001      -0.09999999999999974 ]
[      -0.10000000000000001      -0.10000000000000028 ]
[      -0.10000000000000001      -0.10000000000000034 ]
[      -0.10000000000000001      -0.10000000000000009 ]
[       0.10000000000000001       0.10000000000000010 ]
[      -0.10000000000000001      -0.09999999999999966 ]
[       0.10000000000000001       0.09999999999999963 ]
[       0.10000000000000001       0.10000000000000006 ]
[      -0.10000000000000001      -0.09999999999999999 ]
[       0.10000000000000001       0.09999999999999992 ]
[      -0.10000000000000001      -0.10000000000000013 ]
[       0.10000000000000001       0.09999999999999955 ]
[       0.10000000000000001       0.10000000000000017 ]
[      -0.10000000000000001      -0.09999999999999984 ]
[       0.10000000000000001       0.09999999999999992 ]
[       0.10000000000000001       0.10000000000000009 ]
[      -0.10000000000000001      -0.09999999999999991 ]
[      -0.10000000000000001      -0.09999999999999974 ]
[      -0.10000000000000001      -0.09999999999999989 ]
[      -0.10000000000000001      -0.09999999999999996 ]
[      -0.10000000000000001      -0.09999999999999988 ]
[       0.10000000000000001       0.10000000000000005 ]
[      -0.10000000000000001      -0.10000000000000005 ]
[      -0.10000000000000001      -0.10000000000000006 ]
[      -0.10000000000000001      -0.10000000000000041 ]
[       0.10000000000000001       0.10000000000000027 ]
[       0.10000000000000001       0.09999999999999996 ]
[      -0.10000000000000001      -0.10000000000000007 ]
[       0.10000000000000001       0.09999999999999994 ]
[      -0.10000000000000001      -0.09999999999999999 ]
[       0.10000000000000001       0.10000000000000002 ]
[      -0.10000000000000001      -0.10000000000000012 ]
[       0.10000000000000001       0.10000000000000021 ]
[       0.10000000000000001       0.10000000000000005 ]
[      -0.10000000000000001      -0.09999999999999999 ]
[       0.10000000000000001       0.10000000000000006 ]
[       0.10000000000000001       0.10000000000000005 ]
[      -0.10000000000000001      -0.10000000000000005 ]
[      -0.10000000000000001      -0.10000000000000001 ]
[       0.10000000000000001       0.10000000000000006 ]
[       0.10000000000000001       0.09999999999999996 ]
[       0.10000000000000001       0.09999999999999989 ]
[       0.10000000000000001       0.10000000000000084 ]


### Why does Matlab & Julia use Gaussian Elimination then?¶

These problems never occur!

In [30]:
n = 60;
xtrue = (2*(rand(n,1) .< 0.5)-1)/10;
b = A*xtrue;
xcomp = A\b;

@printf("\n  n = %i\n\n", n);
@printf("    [ %25s %25s ]\n", "x true", "x computed");
for i=1:n
@printf("    [ %25.17f %25.17f ]\n", xtrue[i], xcomp[i]);
end

  n = 60

[                    x true                x computed ]
[      -0.10000000000000001      -0.10000000000000009 ]
[      -0.10000000000000001      -0.09999999999999978 ]
[       0.10000000000000001       0.09999999999999985 ]
[       0.10000000000000001       0.10000000000000003 ]
[       0.10000000000000001       0.09999999999999996 ]
[      -0.10000000000000001      -0.09999999999999999 ]
[       0.10000000000000001       0.10000000000000010 ]
[       0.10000000000000001       0.09999999999999984 ]
[      -0.10000000000000001      -0.10000000000000001 ]
[      -0.10000000000000001      -0.10000000000000006 ]
[       0.10000000000000001       0.10000000000000021 ]
[       0.10000000000000001       0.09999999999999995 ]
[       0.10000000000000001       0.10000000000000009 ]
[      -0.10000000000000001      -0.10000000000000009 ]
[       0.10000000000000001       0.10000000000000007 ]
[      -0.10000000000000001      -0.10000000000000009 ]
[       0.10000000000000001       0.09999999999999992 ]
[      -0.10000000000000001      -0.10000000000000007 ]
[       0.10000000000000001       0.10000000000000001 ]
[       0.10000000000000001       0.10000000000000001 ]
[      -0.10000000000000001      -0.10000000000000001 ]
[      -0.10000000000000001      -0.10000000000000009 ]
[       0.10000000000000001       0.10000000000000001 ]
[      -0.10000000000000001      -0.09999999999999998 ]
[       0.10000000000000001       0.10000000000000003 ]
[       0.10000000000000001       0.10000000000000019 ]
[       0.10000000000000001       0.09999999999999991 ]
[      -0.10000000000000001      -0.09999999999999992 ]
[       0.10000000000000001       0.10000000000000003 ]
[      -0.10000000000000001      -0.10000000000000003 ]
[      -0.10000000000000001      -0.10000000000000003 ]
[      -0.10000000000000001      -0.09999999999999995 ]
[      -0.10000000000000001      -0.10000000000000002 ]
[       0.10000000000000001       0.10000000000000005 ]
[       0.10000000000000001       0.09999999999999994 ]
[       0.10000000000000001       0.09999999999999987 ]
[      -0.10000000000000001      -0.10000000000000009 ]
[      -0.10000000000000001      -0.09999999999999999 ]
[       0.10000000000000001       0.10000000000000010 ]
[       0.10000000000000001       0.09999999999999985 ]
[       0.10000000000000001       0.09999999999999996 ]
[       0.10000000000000001       0.10000000000000009 ]
[      -0.10000000000000001      -0.09999999999999989 ]
[      -0.10000000000000001      -0.09999999999999998 ]
[      -0.10000000000000001      -0.10000000000000016 ]
[       0.10000000000000001       0.10000000000000010 ]
[       0.10000000000000001       0.10000000000000009 ]
[       0.10000000000000001       0.10000000000000007 ]
[      -0.10000000000000001      -0.10000000000000009 ]
[       0.10000000000000001       0.09999999999999985 ]
[      -0.10000000000000001      -0.09999999999999977 ]
[       0.10000000000000001       0.10000000000000014 ]
[       0.10000000000000001       0.09999999999999995 ]
[       0.10000000000000001       0.09999999999999995 ]
[       0.10000000000000001       0.09999999999999992 ]
[       0.10000000000000001       0.10000000000000009 ]
[       0.10000000000000001       0.09999999999999977 ]
[      -0.10000000000000001      -0.10000000000000003 ]
[      -0.10000000000000001      -0.09999999999999978 ]
[       0.10000000000000001       0.09999999999999996 ]