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

x = [ 4 5 ]; % a row 2-vector

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

n1 =    6.4031

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

n1 =    3.8730

n2 =     4

n3 =     5

·         cond(A): computes condition numbers of a matrix with respect to inversion

C1 = cond(A,'fro')

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

C1 =      3.0000

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.

for n = 1 : 15

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)

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

b = ones(5,1);

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

Q1 =    -2

Q2 =    -2

B = A^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

B =

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