Adequacy in Solution of Simultaneous Linear Equations

Jamie Trahan, Autar Kaw, Kevin Martin 

University of South Florida

United States of America

kaw@eng.usf.edu

Introduction

The condition number allows one to quantify the accuracy in solution of [A][X]=[C], where (Subscript[[A], nxn])is an invertible square matrix, (Subscript[[X], nx1])is the solution vector, and (Subscript[[C], nx1]) is the right hand side array. Multiply the condition number by machine epsilon and compare the result to 0.5 x (10^(-m)) to find out at least how many m significant digits are at least correct in solution

0.5 x (10^(-m))  <  Cond(A)*
          machine ϵ                                  Equation (1.1)

To learn more about the relationship between condition number and the adequacy of solution, click here.

The following simulation uses three different techniques to determine the condition number of coefficient matrix (Subscript[[A], nxn]).

Section 1: Input data

Below are the input parameters to begin the simulation. This is the only section where the user can interact with the worksheet. The user can change those values that are highlighted and Mathematica will calculate the condition number of (Subscript[[A], nxn])using an exact method as well as two approximate methods.

nxn invertible square matrix, [A]

A = 
        Table[{{8, 12, 1, 7}, {1, 2, 3.001, 6}, {57, 9, 1, 4}, {1.047, 108, 
              98, 5}}];
A // MatrixForm

Number of bits used for mantissa in floating point representation

MantissaBits = 23

Order of Matrix

n = 4

Section 2: Calculating the condition number

In this section, three distinct methods are used to calculate the condition number of coefficient matrix [A]. Each has its own advantages while utilizing theorems that relate the norm of a matrix to the conditioning of the matrix. Complete details for finding the norm of a matrix and its relationship to the conditioning of a matrix are given here.

Off[General::spell1]

Method 1: Finding the exact value of the Condition Number of a matrix

The following method finds the exact condition number for a square matrix. The exact formula is given by

Cond(A) = ||A|| * ||(A^(-1))||                                Equation (2.1)

Once the condition number is calculated, it can then be used to solve for m, the number of significant digits that one can trust in solution.

Please note that, although this is the most direct method, it may not be practical in its computational time for higher order matrices because this method requires calculation of the inverse of coefficient matrix (Subscript[[A], nxn]). The problem in finding the inverse lies in solving n sets of n equations which can be computationally intensive for large coefficient matrices.

Calculating the condition number:

Calculating the inverse of [A]:

Ainverse = Inverse[A];

Calculating the norm of [A], or ||A||:

NA = Norm[A, Infinity];

Calculating the norm of [(A^(-1))], or ||(A^(-1))||:

NAInv = Norm[Ainverse, Infinity];

Calculating the condition number of [A] using Equation (2.1):

CondA = NA*NAInv

Calculating machine epsilon:

MachineEpsilon = N[2^-MantissaBits]

Now that the condition number of [A] and machine epsilon have been found, Equation (1.1) can be applied to determine the number of significant digits, m, that are at least correct in solution.

Calculating the number of significant digits that one can trust in solution:

Solve[0.5*10^(-m) == MachineEpsilon*CondA, m]

SigDigits = Floor[m /. %]

Because m can be a negative value, the following command returns the appropriate number of significant digits that can be trusted.

TrustDigits = Max[SigDigits, 0]

Method 2: Finding an approximate value of the condition number

The following numerical method finds the condition number of [A] using the inequality

\(TraditionalForm\`\(\(\ \)\(Cond(
          A) ≥ \ \(\(||\)\(Δ\ X\)\(||\)\(\ \)\(\(/\)\(\ \)\(\(||\)\(X + ΔX\)\(||\)\)\)\)\/\(\(||\)\(ΔC\)\(||\)\(\ \)\(\(/\)\(\ \)\(\(||\)\(C\)\(||\)\)\)\)\)\)\)          Equation (2.2)

However, the value of (X+ΔX) is equivalent to X'. Therefore, the inequality becomes

\(TraditionalForm\`\(\(\ \ \)\(Cond(
          A) ≥ \ \(\(||\)\(ΔX\)\(||\)\(\ \)\(\(/\)\(\ \)\(\(||\)\(X'\)\(||\)\)\)\)\/\(\(||\)\(ΔC\)\(||\)\(\ \)\(\(/\)\(\ \)\(\(||\)\(C\)\(||\)\)\)\)\)\)\)               Equation (2.3)

where ||ΔX|| / ||X'|| is the relative change in the norm of the solution vector and ||ΔC|| / ||C|| is the relative change in the norm of the right hand side vector. The ratio between these two values quantifies the conditioning of a system of equations, demonstrating the accuracy in solution. That is, for any small change made in the right hand side array, the resulting change in the solution vector will govern how accurate the system is and therefore how many significant digits one can trust in the solution of a system of simultaneous linear equations.

In this method, the condition number is calculated using Equation (2.3) by first conducting the following steps:

1) A right hand side vector [C] is chosen such that the solution vector equals 1 (i.e. [X] = [1,1,...,1]).

2) A new, unbiased right hand side vector [C'] is then generated by Mathematica. This is done by adding a random positive or negative value to each element of the old right hand side vector.

3) The new right hand side vector [C'] can then be used to calculate a new solution vector [X'].

By creating a small relative change in the right hand side array (i.e., ||ΔC||/||C|| < 1), the magnitude of the condition number will be largely influenced by the relative error in the solution vector, demonstrating how accurate the solution actually is.

NOTE: Each time the worksheet is executed, Mathematica will generate a new condition number. The user should run the worksheet several times to see if the estimate of the condition number approaches the exact value given in Method 1. The greatest of the generated values will be the most accurate approximation and will never exceed the true condition number due to the above inequality, Equation (2.3).

Parameter names:
RHS = old right hand side array [RHS]
RHS1 = new right hand side array [RHS']
X = old solution vector [X]
X1 = new solution vector [X']

Defining all vectors:

RHS = Array[0, n];
RHS1 = Array[0, n];
X = Array[0, n];
X1 = Array[0, n];

Step 1: Assigning values to [RHS] so that each element of [X] equals unity.

For[i = 1, i ≤ n, i++, 
      summ = 0; 
For[j = 1, j ≤ n, j++, 
        summ = summ + A[[i, j]]]; 
RHS[[i]] = 
        summ]; Print["RHS", "=", RHS // MatrixForm]

Solving for [X].

X = LinearSolve[A, RHS];
Print["X", "=", X // MatrixForm]

Step 2: Generating a new, unbiased right hand side vector [RHS'] by adding a random + value between -0.001*RHS[i] and +0.001*RHS[i] to create a small ΔC value.

For[i = 1, i ≤ n, i++, 
        SignVal = Random[Integer, {1, 2}]; 
RHS1[[i]] = 
          RHS[[i]] + (-1)^SignVal*0.001*RHS[[i]]*
              Random[]];
Print["RHS1", "=", RHS1 // MatrixForm]

Step 3: Calculating the new solution vector [X'].

X1 = LinearSolve[A, RHS1];
Print["X1", "=", X1 // MatrixForm]

Notice the small difference between the values of the old and new right hand side vectors. Now look at the new solution vector [X1]. Do the values deviate much from the old solution vector [X]? If not, then any small changes that are made in the right hand side array do not affect the accuracy of the solution vector, and the solution can therefore be trusted. Otherwise, the system of equations is ill-conditioned. Below, the condition number of [A] is calculated by determining the values required by Equation (2.3).

Calculating the relative change in the norm of the solution vector, ||ΔX|| / ||X'||.

NumCond = Norm[X1 - X, Infinity]/Norm[X1, Infinity]

Calculating the relative change in the norm of the right hand side array, ||ΔC|| / ||C||.

DenCond = Norm[RHS1 - RHS, Infinity]/Norm[RHS, Infinity]

Calculating the condition number using Equation (2.3).

CondAA = NumCond/DenCond

Method 3: A second technique in approximating the condition number

The following numerical technique is simpler than Method 2 as it requires fewer calculations. This method utilizes the theorem

\(TraditionalForm\`\(\(\ \)\(Cond(
          A) ≥ \ \(\(||\)\(A\)\(||\)\(\ \)\(||\)\(X\)\(||\)\)\/\(\(||\)\(C\)\(||\)\)\)\)\)                                       Equation (2.4)

The proof is as follows.

Let
[X] = [(A^(-1))] [C]

Applying norm properties, this equation becomes
\(TraditionalForm\`\(\(||\)\(X\)\(||\)\(\ \)\(\(≤\)\(\ \)\(\(||\)\(A\^\(-1\)\)\(||\)\(\ \)\(||\)\(C\)\(||\)\)\)\)\)

Multiplying ||A|| to both sides gives
||A|| ||X||  ≤ Cond(
              A)* || C || 
by the definition of condition number, and division by ||C|| results in the inequality of Equation (2.4).

In this technique, however, the right hand side values of [C] matrix are chosen to equal [+1, +1, +1] with signs generated randomly. This will result in ||C||=1, minimizing the number of calculations required to solve for Cond(A). Therefore the calculation is reduced to

 Cond(
            A) ≥  || A ||   || X ||                                   Equation (2.4)

Again, the user should reexecute the worksheet to see if the greates approximate condition number approaches the exact condition number defined by Method 1.

Randomly generating a (+1) or (-1) value for each element of the right hand side array:

For[i = 1, i ≤ n, i++, 
      signval = Random[Integer, {1, 2}]; 
      RHS[[i]] = (-1)^signval]; 
    Print["RHS", "=", RHS]

Solving for the solution vector [X].

X = LinearSolve[A, RHS];
Print["X", "=", X // MatrixForm]

Calculating the condition number.

CondAAA = Norm[A, Infinity]*Norm[X, Infinity]

Conclusion

Mathematica helped us apply our knowledge of norms and condition numbers to quantify the accuracy in a solution for a system of simultaneous linear equations.

Question 1: Choose a coefficient matrix [A] that is well-conditioned. For example,

(({{10, -7, 0}, {-3, 2.099, 6}, {5, -1, 5}}))

See how the condition number of the matrix affects the number of significant digits that one can trust in solution.

Question 2: Choose the number of bits used for the mantissa in single and double precision. See how these numbers affect the number of significant digits that one can trust in solution.

References

[1] Autar Kaw, Holistic Numerical Methods Institute, http://numericalmethods.eng.usf.edu/mws, See
System of Equations
What is Machine Epsilon?
Adequacy of Solutions