Ch 3.4. Inverse matrices and linear systems

• Prev: square-matrices-determinant-trace
• Next: vector-spaces-subspaces-rank

1. System of linear equations

Where there is a system of linear equations (i.e. no terms like $x^2,xy,\sqrt{x},\sin (x)$ etc), the equations can be represented as matrix-vector multiplication.

$$\begin{array}{c}\underset{\text{Unknown variables}}{\underbrace{\begin{array}{c}\\ \\ xyz\end{array}}}\underset{\text{Equations}}{\underbrace{\begin{array}{c}2x+5y+3z=-3\\ 4x+0y+8z=0\\ 1x+3y+0z=2\end{array}}}\to \stackrel{\text{A}}{\overbrace{\underset{\text{Coefficients}}{\underbrace{\left[\begin{array}{ccc}2& 5& 3\\ 4& 0& 8\\ 1& 3& 0\end{array}\right]}}}}\stackrel{\text{x}}{\overbrace{\underset{\text{Variables}}{\underbrace{\left[\begin{array}{c}x\\ y\\ z\end{array}\right]}}}}=\stackrel{\text{b}}{\overbrace{\underset{\text{Constants}}{\underbrace{\left[\begin{array}{c}-3\\ 0\\ 2\end{array}\right]}}}}\end{array}$$

• To solve $\text{Ax}=\text{b}$, means to find the vector $\text{x}$. To do this, we must find the inverse of $\text{A}$. (This is not always possible).

Geometric intuition

There are two ways to interpret the above system of equations, and find the solution, $\text{x}$:

1. Vector interpretation - Find the values of $x,y,z$ that correctly scale the columns of matrix $\text{A}$ to produce the vector $\text{b}$
   1. Matrix $\text{A}$ linearly transforms the ${\mathbb{R}}^{n=3}$ vector space in such a way that the vector $\text{x}$ gets transformed into $\text{b}$.
   2. i.e. find the linear combination of the columns of $\text{A}$ that equals $\text{b}$:

$$\begin{array}{r}x\left[\begin{array}{c}2\\ 4\\ 1\end{array}\right]+y\left[\begin{array}{c}5\\ 0\\ 3\end{array}\right]+z\left[\begin{array}{c}3\\ 8\\ 0\end{array}\right]=\left[\begin{array}{c}-3\\ 0\\ 2\end{array}\right]\end{array}$$

2. Planar interpretation - Find the intersection of three planes. Each equation represents a plane in ${\mathbb{R}}^3$.
   1. The solution to the system of equations is the set of points that lie on all three planes.
   2. If the solution is unique, then the planes intersect at a single point.
   3. Otherwise, the planes are parallel (no solution) or coincide as a line or plane (infinitely many solutions).

$$\begin{array}{r}\begin{array}{c}2x+5y+3z=-3\\ 4x+0y+8z=0\\ 1x+3y+0z=2\end{array}\end{array}$$

2. Matrix inverse

2.1. Is it invertible?

• All non-square matrices are not invertible (they have no determinant). Some square matrices are also not invertible:
• Singular matrices are those which have no inverses (like how 0 has no inverse)
  • $\text{det}(\text{A})=0$: No inverse, because the matrix collapses the vector space (dimensionality is lost)
• Non-singular matrices are those which have an inverse
  • $\text{det}(\text{A})\ne 0$: Inverse exists, because the matrix preserves the vector space (dimensionality is maintained)
  • The inverse of a matrix is unique: An invertible matrix only has one inverse.

2.1.1. Caveat:

• Solutions may exist even when the matrix is singular (i.e. $\text{det}(\text{A})=0$).
  • In this case, the system of equations is underdetermined (i.e. there are more unknowns than equations).
  • Geometrically, despite a singular matrix collapsing the vector space, the following solutions may still exist:
    • If a 2D vector space is collapsed into a 1D line, the solution space is a 1D line
    • or if a 3D vector space is collapsed into a 2D plane (or 1D line), the solution space is a 2D plane (or a 1D line)

2.2. Finding the inverse

• Only square, non-singular matrices (${\mathbb{R}}^{n\times n}$) are invertible.
• $\text{A}$ is invertible if there exists a matrix ${\text{A}}^{-1}$ such that $\text{A}{\text{A}}^{-1}={\text{A}}^{-1}\text{A}={\text{I}}_n$.
  • akin to how how $3\times \frac{1}{3}=1$
• Geometric intuition: If you apply $\text{A}$ and then apply its inverse ${\text{A}}^{-1}$, you get back to the original vector. In other words, you did “nothing” to the vector, hence the identity matrix ${\text{I}}_n$.
• For a $2\times 2$ matrix, the inverse is:

$$\begin{array}{r}{\text{A}}^{-1}={\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]}^{-1}=\frac{1}{|\text{A}|}\times \text{adj(}\text{A}\text{)}=\frac{1}{ad-bc}\left[\begin{array}{cc}d& -b\\ -c& a\end{array}\right]\end{array}$$

• Complicated: Where $\text{adj(}\text{A}\text{)}$ is the adjugate matrix of a matrix $\text{A}$ (i.e. the transpose of the cofactor matrix of $\text{A}$). See link for computation details: Adjugate matrix

2.3. Using the inverse to solve a system of linear equations

Using the earlier example:

$$\begin{array}{r}\stackrel{\text{A}}{\overbrace{\underset{\text{Coefficients}}{\underbrace{\left[\begin{array}{ccc}2& 5& 3\\ 4& 0& 8\\ 1& 3& 0\end{array}\right]}}}}\stackrel{\text{x}}{\overbrace{\underset{\text{Variables}}{\underbrace{\left[\begin{array}{c}x\\ y\\ z\end{array}\right]}}}}=\stackrel{\text{b}}{\overbrace{\underset{\text{Constants}}{\underbrace{\left[\begin{array}{c}-3\\ 0\\ 2\end{array}\right]}}}}\end{array}$$

We can solve for $\text{x}$ by multiplying both sides by the inverse of $\text{A}$:

$$\begin{array}{rl}\text{A}\text{x}& =\text{b}\\ {\text{A}}^{-1}\text{A}\text{x}& ={\text{A}}^{-1}\text{b}\\ \text{x}& ={\text{A}}^{-1}\text{b}\end{array}$$

Finding the inverse of $\text{A}$:

$$\begin{array}{rl}{\text{A}}^{-1}& =\frac{1}{|\text{A}|}\times \text{adj(}\text{A}\text{)}=\text{np.linalg.inv(A)}\end{array}$$

2.3.1. Notes on computational complexity and numerical stability

• Computational complexity and memory-intensive: Inverting a matrix is computationally expensive. Inverting large matrices may cause memory issues. Use the inverse of a matrix only when necessary.
• Numerical stability: Inverting a matrix can lead to numerical instability. Use other methods (e.g. LU decomposition, iterative methods) for solving systems of linear equations.

import numpy as np
from numpy.linalg import det, inv
 
A = np.array([[2, 5, 3], [4, 0, 8], [1, 3, 0]])
b = np.array([[-3], [0], [2]])
# b = np.array([-3, 0, 2]) # This is also valid because python is smart and can infer the required array shape
print("Problem statement: Ax = b. Find x, given:\nA =\n", A, "\nb =\n", b)
 
# Calculate the determinant of A
print(f"\ndet(A) = {det(A):.1f}")  # A has non-zero determinant so is invertible
 
# Find the inverse of A (round to 2 dp)
print("\nInv(A): A^-1 =\n", np.round(inv(A), 3))
 
# Demonstrate that A^-1 * A = I, and also A * A^-1 = I (round to 2 dp)
print("\nDemonstrate that A^-1 is in fact the inverse of A:")
print("A^-1 * A =\n", np.round(inv(A) @ A, 10))
print("\n(And also):")
print("A * A^-1 =\n", np.round(A @ inv(A), 10))
 
# Find x in Ax = b
x = inv(A) @ b
print("\nRecall, problem statement: Ax = b. Therefore:")
print("x = A^-1 * b =\n", np.round(x, 3))
 
# Check the solution
print("\nCheck the solution: Ax =\n", A @ x)
 

Problem statement: Ax = b. Find x, given:
A =
 [[2 5 3]
 [4 0 8]
 [1 3 0]] 
b =
 [[-3]
 [ 0]
 [ 2]]
 
det(A) = 28.0
 
Inv(A): A^-1 =
 [[-0.857  0.321  1.429]
 [ 0.286 -0.107 -0.143]
 [ 0.429 -0.036 -0.714]]
 
Demonstrate that A^-1 is in fact the inverse of A:
A^-1 * A =
 [[ 1. -0.  0.]
 [ 0.  1.  0.]
 [ 0.  0.  1.]]
 
(And also):
A * A^-1 =
 [[ 1.  0. -0.]
 [ 0.  1.  0.]
 [-0.  0.  1.]]
 
Recall, problem statement: Ax = b. Therefore:
x = A^-1 * b =
 [[ 5.429]
 [-1.143]
 [-2.714]]
 
Check the solution: Ax =
 [[-3.]
 [ 0.]
 [ 2.]]

# Demonstrate what happens when square matrix A is singular (non-invertible, i.e. determinant is zero)
P = np.array([[0, 1, 0], [0, 0, 0], [1, 0, 1]])
print("det(p):", det(P))
print("Inv P:\n", inv(P))  # <-- LinAlgError thrown because P is Singular (not invertible)
 

det(p): 0.0

---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
Cell In[2], line 4
      2 P = np.array([[0, 1, 0], [0, 0, 0], [1, 0, 1]])
      3 print("det(p):", det(P))
----> 4 print("Inv P:\n", inv(P))  # <-- LinAlgError thrown because P is Singular (not invertible)
 
File ~/Code/notes/.venv/lib/python3.13/site-packages/numpy/linalg/_linalg.py:609, in inv(a)
    606 signature = 'D->D' if isComplexType(t) else 'd->d'
    607 with errstate(call=_raise_linalgerror_singular, invalid='call',
    608               over='ignore', divide='ignore', under='ignore'):
--> 609     ainv = _umath_linalg.inv(a, signature=signature)
    610 return wrap(ainv.astype(result_t, copy=False))
 
File ~/Code/notes/.venv/lib/python3.13/site-packages/numpy/linalg/_linalg.py:104, in _raise_linalgerror_singular(err, flag)
    103 def _raise_linalgerror_singular(err, flag):
--> 104     raise LinAlgError("Singular matrix")
 
LinAlgError: Singular matrix

2.4. Ill-Conditioned Matrices

• An ill-conditioned matrix is one which is close to being singular
  • Its determinant will be close to 0 (problematic in the same way dividing by a tiny number is)
  • Computation errors (overflow, underflow, round-off errors) may occur
• Condition number: Higher number means the matrix is more ill-conditioned (i.e. closer to being singular)

# Compute the condition numbers, determinants, and matrix inverses for the following 3 matrices A_1, A_2, A_3
 
from numpy.linalg import cond
 
A_1 = np.array([[1, 1], [2, 3]])
A_2 = np.array([[4.1, 2.8], [9.7, 6.6]])
A_3 = np.array([[4.1, 2.8], [9.6760, 6.6080]])
 
print("\nMatrix A_1:")
print("Condition number(A_1):", cond(A_1))
print("det(A_1):", det(A_1))
print("Inv A_1:\n", inv(A_1))
 
print("\nMatrix A_2:")
print("Condition number(A_2):", cond(A_2))
print("det(A_2):", det(A_2))
print("Inv A_2:\n", inv(A_2))
 
print("\nMatrix A_3 -- This matrix is close to singular or badly scaled. It is ill-conditioned:")
print("Condition number(A_3):", cond(A_3))  # Note that the condition number is very high
print("det(A_3):", det(A_3))  # Note that the determinant is very close to zero
print("Inv A_3:\n", inv(A_3))  # Note that the inverse is very large
 

Matrix A_1:
Condition number(A_1): 14.933034373659225
det(A_1): 1.0
Inv A_1:
 [[ 3. -1.]
 [-2.  1.]]
 
Matrix A_2:
Condition number(A_2): 1622.9993838565622
det(A_2): -0.0999999999999999
Inv A_2:
 [[-66.  28.]
 [ 97. -41.]]
 
Matrix A_3 -- This matrix is close to singular or badly scaled. It is ill-conditioned:
Condition number(A_3): 1.1149221731402912e+16
det(A_3): -4.461167435465537e-15
Inv A_3:
 [[-1.53781451e+15  6.51616316e+14]
 [ 2.25179981e+15 -9.54152463e+14]]