CHAPTER 2: MATRIX COMPUTATIONS AND LINEAR ALGEBRA
Lecture 2.4: Ill-conditioned and sparse matrices
·
Numerical
computations are always subject to rounding errors generated due to finite
machine precision, e.g. eps ~ 10-16. The rounding
error can be magnified for ill-conditioned linear systems.
A = hilb(15); % Hilbert
matrices are ill-conditioned
b1 = ones(15,1); % b1 is the vector of ones
b2 = b1 + 0.01*rand(15,1); % b2 is a small random perturbation to
b1
x1 = A\b1; y1 = x1'
x2 = A\b2; y2 = x2'
% Two solutions of linear
systems with almost equal right-hand-sides
bdif = norm(b2-b1), ydif = norm(y2-y1)
Warning: Matrix is close to singular or badly scaled.
Results may be
inaccurate. RCOND = 1.543404e-018.
y1 = 1.0e+008 *
0.0000 -0.0000
0.0007 -0.0097 0.0738
-0.3122 0.7102 -0.5709
-1.0953 3.3002 -2.9619
-0.2353 2.4214 -1.7331
0.4121
Warning: Matrix is close to singular or badly scaled.
Results may be
inaccurate. RCOND = 1.543404e-018.
y2 = 1.0e+014 *
-0.0000 0.0000
-0.0001 0.0015 -0.0157
0.0998 -0.4111 1.1384
-2.1512 2.7591 -2.3245
1.1786 -0.2709 -0.0211
0.0172
bdif = 0.0205
ydif = 4.5367e+014
·
Large
errors for ill-conditioned systems are generated by small denominators: xn
= det An/detA. If det A ~ 0, and det An
is computed with an error en, then the
absolute error is magnified: |xn – xnexact|
= en/det A ~ large. Solutions of
ill-conditioned systems become inaccurate with large absolute error.
Norms for vectors and matrices:
·
norm(x): computes norms for a vector x
n1 = norm(x) % computes a geometric length (2-norm) of a vector:
% norm(x) = norm(x,2) =
sqrt(x(1)^2 + x(2)^2)
n2 = norm(x,1) % computes
a 1-norm of a vector: norm(x,1) = |x(1)| + |x(2)|
n3 = norm(x,inf) % computes an infinity-norm of a vector:
% norm(x,inf) = max(|x(1)|,|x(2)|)
n2 = 9
n3 = 5
·
norm(A): computes norms for a matrix A
A = [ 3, 2;
1,-1 ]; % a 2-by-2 matrix
n1 = norm(A,'fro') % computes the Frobenius norm of a matrix:
% norm(A,'fro') =
sqrt(sum(diag(A'*A)))
n2 = norm(A,1) % computes the 1-norm of a matrix (the largest
column sum):
% norm(A,1) = max(sum(abs((A))))
n3 = norm(A,inf) % computes the infinity norm of a matrix (the
largest row sum):
% norm(A,inf) = max(sum(abs((X'))))
n2 = 4
n3 = 5
·
cond(A): computes condition numbers of a matrix with respect to inversion
C2 = norm(A,'fro')*norm(inv(A),'fro')
% Properties of the condition number:
% cond(A) =
norm(A)*norm(A^(-1)) in any norm
% cond(A) >= 1
C2 = 3.0000
Criteria
of ill-conditioned linear systems:
Two theoretical criteria of ill-conditioned systems
independent of the computing environment are:
·
large
values of cond(A)
·
small
values of det(A)
Two practical criteria of ill-conditioned systems in the
given computing environment are:
·
cond(A*inv(A)) ~= 1
·
det(A)*det(inv(A)) ~= 1
If
high precision is used in computational environment, solutions of linear
systems could be accurate even if the theoretical criteria predict that the
coefficient matrix A is ill conditioned. On contrary, the
practical criteria predict that the actual error in solutions of linear systems
is large in the given computational environment.
Since
inv(A) is used for solutions of the linear system, the product inv(A)*A
or A*inv(A) is used for practical measures of the error in the
given computing environment. The products inv(A)*A or A*inv(A)
must reproduce an identity matrix. The difference between the identity matrix
and inv(A)*A, if any exists, indicates a computational error in solving the linear
systems with the coefficient matrix A.
For a given ill-conditioned system, an increase in precision will reduce errors of solutions of the linear system. Double precision is used in MATLAB and it is sufficient for mildly ill-conditioned systems. Quadruple precision (not available in MATLAB) should be used for severely ill-conditioned systems.
A = hilb(n); B = inv(A); C
= A*B;
c1 = cond(A); c2 = det(A);
c3 = det(B);
d1 = c2*c3; d2 = cond(C);
fprintf('n = %2.0f c1 = %3.1e, c2 = %3.1e, d1 = %3.1e, d2
=%3.1e\n',n,c1,c2,d1,d2);
end
n = 1 c1 = 1.0e+000, c2 =
1.0e+000, d1 = 1.0e+000, d2 =1.0e+000
n = 2 c1 = 1.9e+001, c2 =
8.3e-002, d1 = 1.0e+000, d2 =1.0e+000
n = 3 c1 = 5.2e+002, c2 =
4.6e-004, d1 = 1.0e+000, d2 =1.0e+000
n = 4 c1 = 1.6e+004, c2 =
1.7e-007, d1 = 1.0e+000, d2 =1.0e+000
n = 5 c1 = 4.8e+005, c2 =
3.7e-012, d1 = 1.0e+000, d2 =1.0e+000
n = 6 c1 = 1.5e+007, c2 =
5.4e-018, d1 = 1.0e+000, d2 =1.0e+000
n = 7 c1 = 4.8e+008, c2 =
4.8e-025, d1 = 1.0e+000, d2 =1.0e+000
n = 8 c1 = 1.5e+010, c2 =
2.7e-033, d1 = 1.0e+000, d2 =1.0e+000
n = 9 c1 = 4.9e+011, c2 =
9.7e-043, d1 = 1.0e+000, d2 =1.0e+000
n = 10 c1 = 1.6e+013, c2 = 2.2e-053, d1 = 1.0e+000, d2 =1.0e+000
n = 11 c1 = 5.2e+014, c2 = 3.0e-065, d1 = 1.0e+000, d2 =1.7e+000
Warning: Matrix is close to singular or badly scaled.
Results may be
inaccurate. RCOND = 2.632091e-017.
n = 12 c1 = 1.8e+016, c2 = 2.9e-078, d1 = 9.4e-001, d2 =2.3e+002
Warning: Matrix is close to singular or badly scaled.
Results may be
inaccurate. RCOND = 2.348790e-018.
n = 13 c1 = 3.8e+018, c2 = 4.5e-092, d1 = 8.4e-001, d2 =1.4e+004
Warning: Matrix is close to singular or badly scaled.
Results may be
inaccurate. RCOND = 1.408541e-019.
n = 14 c1 = 4.1e+017, c2 = -3.2e-107, d1 = -3.3e-001, d2 =4.2e+006
Warning: Matrix is close to singular or badly scaled.
Results may be
inaccurate. RCOND = 1.543404e-018.
n = 15 c1 = 8.5e+017, c2 = -2.2e-120, d1 = 1.6e-001, d2 =2.1e+004
Remark:
In the
example above, the condition number c1 is huge and the
determinant c2 is small for n > 5. The actual
error is however invisible in the given computing environment, since d1 and
d2 do not differ from the theoretical values: d1 = d2 = 1.
The computational error becomes visible for n > 10. Solutions
of linear systems with Hilbert matrices A become inaccurate in
the given computing environment for n > 10.
Sparse
matrices:
·
Sparse
matrices have zeros in most of their entries.
·
Sparse
matrices can be represented with reduced storage in a computer
·
All
matrix multiplications with sparse matrices can be done with a reduced number
of operations
·
sparse:
converts a sparse matrix to sparse array form by squeezing out any zero
elements. MATLAB operates with sparse arrays as with original matrices
but save storage and computational time as it exploits the sparse structure of
a matrix
A =
-2*diag(ones(1,5)) + diag(ones(1,4),1) + diag(ones(1,4),-1)
S = sparse(A)
-2 1
0 0 0
1 -2
1 0 0
0 1
-2 1 0
0 0
1 -2 1
0 0
0 1 -2
S =
(1,1) -2
(2,1) 1
(1,2) 1
(2,2) -2
(3,2) 1
(2,3) 1
(3,3) -2
(4,3) 1
(3,4) 1
(4,4) -2
(5,4) 1
(4,5) 1
(5,5) -2
Q1 = b'*A*b % computations of a quadratic form with vector b
Q2 = b'*S*b % alternative
computation of the same quadratic form
% the second computation uses smaller number of floaping point
operations
Q2 = -2
% B is central-difference approximation for fourth-derivative
% u''''(x) = (u(x+2h)-4*u(x+h)+6*u(x)-4*u(x-h)+u(x-2*h))/h^4
% B = A*A ~= A.^2
% Matrix power is not the same as the pointwise power to elements
of matrix
R = S^2
% the same operation can be performed with sparse arrays to save computer time
5 -4 1 0
0
-4 6
-4 1 0
1 -4
6 -4 1
0 1
-4 6 -4
0 0
1 -4 5
R =
(1,1) 5
(2,1) -4
(3,1) 1
(1,2) -4
(2,2) 6
(3,2) -4
(4,2) 1
(1,3) 1
(2,3) -4
(3,3) 6
(4,3) -4
(5,3) 1
(2,4) 1
(3,4) -4
(4,4) 6
(5,4) -4
(3,5) 1
(4,5) -4
(5,5) 5