4.5. Errors in Linear Systems#

4.5.1. Learning Objectives#

After studying this notebook, completing the activities, and asking questions in class, you should be able to:

  • Interpret the rank of a matrix. What does it mean about the corresponding linear system?

  • Propagate uncertainty through the solution of a linear system using condition number.

import numpy as np

4.5.2. Errors in Linear Systems#

We will now learn about condition number and singular value decomposition as tools to analyze solutions to linear systems. These are helpful for diagnosing if a linear system is rank deficient, i.e., the model is not fully specified (not enough equations or information given to solve).

4.5.3. Rank of the Matrix#

The rank of a matrix is the number of linear independent equations encoded in it. Let’s explore a few examples.

4.5.4. Example 1#

Let’s take a look at the same \(\mathbf{A}\) matrix we used as an example in this class and the previous class. We know this matrix corresponds to a system of three independent linear equations:

\[3 x_1 + 2 x_2 + x_3 = 6\]
\[-x_1 + 4 x_2 + 5x_3 = 8\]
\[2x_1 -8 x_2 + 10 x_3 = 4\]
\[\begin{split}\underbrace{\begin{pmatrix} 3 & 2 & 1\\ -1 & 4 & 5\\ 2& -8 & 10\end{pmatrix}}_{\mathbf{A}} ~ \underbrace{\begin{pmatrix} x_1\\x_2\\x_3\end{pmatrix}}_{\mathbf{x}} = \underbrace{\begin{pmatrix}6\\8\\4\end{pmatrix}}_{\mathbf{b}} \end{split}\]

Note: whether or not each row (equation) is independent does not depend on the right hand side. In other words, we are only focusing on \(\mathbf{A}\).

A1 = np.array([(3.0,2,1),(-1,4,5),(2,-8,10)])
print("A = \n",A1)

A = 
 [[ 3.  2.  1.]
 [-1.  4.  5.]
 [ 2. -8. 10.]]

print("rank(A) =",np.linalg.matrix_rank(A1))

rank(A) = 3

This makes sense. There are three linearly independent equations, so the rank of the matrix is 3.

4.5.5. Example 2#

Now let’s consider a system where equation 2 and 3 are the same:

\[1 x_1 + 1 x_2 + 1 x_3 = 0\]
\[2 x_1 + 1 x_2 + 0 x_3 = 0\]
\[2 x_1 + 1 x_2 + 0 x_3 = 0\]

Clearly, there are only two unique linear equations. Let’s verify the rank is two.

A2 = np.array([[1, 1, 1], [2, 1, 0], [2, 1, 0]])
print("A2 = \n",A2)

print("\nrank(A2) =",np.linalg.matrix_rank(A2))

A2 = 
 [[1 1 1]
 [2 1 0]
 [2 1 0]]

rank(A2) = 2

4.5.6. Example 3#

Let’s consider a more subtle example.

\[1 x_1 + 0 x_2 + 0 x_3 + 1 x_4 = 1\]
\[0 x_1 + 1 x_2 + 0 x_3 + 1 x_4 = 2\]
\[0 x_1 + 0 x_2 + 1 x_3 - 1 x_4 = 3\]
\[2 x_1 + 2 x_2 + 2 x_3 + 2 x_4 = 4\]

Home Activity

Convert the linear system above into a matrix. Store the left hand side in the numpy array (or matrix) A3. Then calculate the rank of A3 and store in r3. Finally, write a sentence as to why the rank makes sense.

# Add your solution here

### HIDDEN TESTS
A3_ans = np.array([[1, 0, 0, 1],[0, 1, 0, 1], [0, 0, 1, -1], [2, 2, 2, 2]])

assert np.shape(A3) == (4,4), "A3 should be a 4x4 matrix"

for i in range(4):
    for j in range(4):
        assert np.abs(A3_ans[i,j] - A3[i,j]) < 1E-6, "Check the ("+str(i)+","+str(j)+") element in A3"

### HIDDEN TESTS

Please see your linear algebra textbook for a more formal definition of rank. This wikipedia article is also quite good: https://en.wikipedia.org/wiki/Rank_(linear_algebra)

4.5.7. Singular Value Decomposition#

Singular Value Decomposition (SVD) is an incredibly important tool in numeric analysis and engineering. SVD is at the core of image compression, singular processing, control and many other important applications. There is a chapter on SVD in your linear algebra textbook. I encourage you to read it.

For our purposes, we will think of SVD as another type of matrix factorization. We can represent a \(m \times n\) matrix as

\[\mathbf{A} = \mathbf{U} \cdot \mathbf{\Sigma} \cdot \mathbf{V}^T\]

where the columns of \(\mathbf{U}\) and \(\mathbf{V}\) are the left and right singular vectors:

\[\mathbf{U} = \begin{bmatrix} \vec{u}_1 & \vec{u}_2 & ... & \vec{u}_m \end{bmatrix}\]
\[\mathbf{V} = \begin{bmatrix} \vec{v}_1 & \vec{v}_2 & ... & \vec{v}_n \end{bmatrix}\]

Moreover, \(\mathbf{\Sigma}\) is a \(m \times n\) matrix of zeros with the singular values on the diagonal. For \(m > n\), we have:

\[\begin{split}\mathbf{\Sigma} = \begin{bmatrix} \sigma_1 & 0 & \dots & 0 \\ 0 & \sigma_2 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \sigma_n \\ 0 & 0 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & 0 \end{bmatrix}\end{split}\]

Thus, we can rewrite the SVD factorization as follows (assuming \(m \geq n\)):

\[\mathbf{A} = \sum_{i=1}^{n} \vec{u}_i \cdot \sigma_i \cdot \vec{v}_i^T\]

Class Activity

What are the dimensions of A, U, V, and \(\mathbf{\Sigma}\)? What changes if \(n < m\)?

Now let’s look at SVD in Python:

A = np.array([(3.0,2,1),(-1,4,5),(2,-8,10),(-2,-8,10)])
print("A =\n",A,"\n")

U,S,V = np.linalg.svd(A)
print("U =\n",U,"\n")
print("S =\n",S,"\n")
print("V =\n",V,"\n")

A =
 [[ 3.  2.  1.]
 [-1.  4.  5.]
 [ 2. -8. 10.]
 [-2. -8. 10.]] 

U =
 [[ 0.02399156 -0.33975594 -0.698885    0.62892771]
 [-0.08707759 -0.93759246  0.25472139 -0.2201247 ]
 [-0.70322927  0.04758521 -0.50005219 -0.50314217]
 [-0.70520244  0.056762    0.4434237   0.55031174]] 

S =
 [18.17773759  6.60309487  4.23898503] 

V =
 [[ 0.00896694  0.60332736 -0.79744322]
 [-0.01514879 -0.7973018  -0.60339071]
 [-0.99984504  0.01749087  0.00199033]]

Notice that Python only returns the first \(n\) elements of the diagonal \(\mathbf{\Sigma}\). Because \(m > n\), we know the remaining \(m-n = 1\) elements are zero.

Let’s check the rank of this matrix.

np.linalg.matrix_rank(A)

3

The matrix rank is 3 and \(n=3\), thus we say this matrix is full rank. Notice SVD returned 3 non-zero singular values.

Now let’s try SVD on a rank deficient matrix.

Anew = np.array([[1, 1, 0],[1, -1, 2],[2, 0, 2]])
print("Anew = \n",Anew,"\n")
print("and has rank",np.linalg.matrix_rank(Anew))

Anew = 
 [[ 1  1  0]
 [ 1 -1  2]
 [ 2  0  2]] 

and has rank 2

U,S,V = np.linalg.svd(Anew)
print("U =\n",U,"\n")
print("S =\n",S,"\n")
print("V =\n",V,"\n")

U =
 [[-0.13550992  0.8051731  -0.57735027]
 [-0.6295454  -0.51994159 -0.57735027]
 [-0.76505532  0.28523152  0.57735027]] 

S =
 [3.64575131e+00 1.64575131e+00 6.69376379e-17] 

V =
 [[-0.6295454   0.13550992 -0.76505532]
 [ 0.51994159  0.8051731  -0.28523152]
 [ 0.57735027 -0.57735027 -0.57735027]]

Notice the almost zero third singular value. This matrix is numerically rank deficient. The number of non-zero singular values is the rank of a matrix.

4.5.8. Condition Number#

The condition number is the ratio of the largest to smallest singular values of a matrix:

\[\kappa = \frac{\sigma_{\mathrm{max}}}{\sigma_{\mathrm{min}}}\]

In some sense, the condition number \(\kappa\) measures how close a matrix is to singular. The condition number also bounds the relative error in solving a linear system:

\[\frac{||\vec{\Delta x}||}{||\vec{x}||} \leq \kappa \frac{||\vec{\Delta b}||}{||\vec{b}||}\]

where \(\vec{\Delta b}\) is a perturbation (error) in the vector \(\vec{b}\) and \(\vec{\Delta x}\) is the corresponding or induced error in \(\vec{x}\). The operator \(|| \cdot ||\) is the \(\ell_2\) norm (i.e., length of a vector) you learned about in physics:

\[|| \vec{a} || = \sqrt{\sum_{i=1}^{N} a_i^2}\]

where \(a\) is an arbitrary \(N \times 1\) vector.

4.5.9. Example#

Let’s revisit the motivating example from last class.

\[3 x_1 + 2 x_2 + x_3 = 6\]
\[-x_1 + 4 x_2 + 5x_3 = 8\]
\[2x_1 -8 x_2 + 10 x_3 = 4\]
# Declare matrix
A33 = np.matrix([(3.0,2,1),(-1,4,5),(2,-8,10)])
b3 = np.matrix([6,8,4])

print("A = \n",A33)
print("\nb =",b3)

A = 
 [[ 3.  2.  1.]
 [-1.  4.  5.]
 [ 2. -8. 10.]]

b = [[6 8 4]]

## Calculate condition number
kappa = np.linalg.cond(A33)
print("condition number =",kappa)

condition number = 3.99784371803463

## Calculate the norm of a vector, such as b3
print("l2 norm of b =",np.linalg.norm(b3))

l2 norm of b = 10.770329614269007

Class Activity

Let’s say in this linear system \(\vec{b}\) comes from an experiment measurement with \(\pm 0.1\) absolute error for each element, thus \(\vec{\Delta b} = [0.1, 0.1, 0.1]^T\). What is the corresponding error when we solve the linear system for \(\vec{x}\). In other words, what is \(||\vec{\Delta x}||\)?

# Add your solution here

4.5.10. Take Home Message#

Rank and condition number tell us sensitivity of solutions to a linear system:

  • Rank tells us how many linear independent equations are encoded in a matrix

  • Singular value decomposition does many things… we’ll use to tell if a matrix is near-singular, i.e., numerically rank deficient

  • Condition number tells us how \(\vec{x}\), the solution to a linear system, changes with uncertainty in \(\vec{b}\)