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.
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
hilb(5)
## 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
## 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
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.
cond(hilb(12))
hilb(12)*xcomp - b
# define a matrix with a bad growth factor
function bad_growth(n)
A = eye(n)
A = A - tril(ones(n,n),-1)
A[:,end] = 1
return A
end
bad_growth(5)
n = 60
A = bad_growth(n)
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
## 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
These problems never occur!
n = 60;
A = bad_growth(n) + randn(n,n)*1e-15;
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