The module numpy.linalg¶

author: Parin Chaipunya affil: KMUTT

We have already compute REF/RREF manually, and from there we could get the rank, nullity, solve linear systems, etc. These tasks and many more (except computing explicity the REF/RREF numerically) are actually built-in with the module numpy.linalg.

In this notebook, we will go through some of these functions.

Rank¶

To compute the rank of a matrix (i.e. a numpy array), we use the function matrix_rank from the module numpy.linalg.

In [1]:
import numpy as np
In [2]:
# Define a matrix A.
A = np.array([
    [1, 3, -1, 0],
    [0, -1, 3, -2],
    [1, 0, 1, -1],
    [1, 1, 1, -2],
    [1, 0, 4, -4]
])
rank_A = np.linalg.matrix_rank(A)
print(f"A = \n{A}\nand rank(A) = {rank_A}.")

A = 
[[ 1  3 -1  0]
 [ 0 -1  3 -2]
 [ 1  0  1 -1]
 [ 1  1  1 -2]
 [ 1  0  4 -4]]
and rank(A) = 4.

Determinant¶

The function det from numpy.linalg can be used to compute the determinant of a matrix.

In [3]:
A = np.array([
    [1, 0, 0, -1],
    [1, 1, 0, 4],
    [1, 0, -2, 0],
    [1, 1, 1, 0]
])
print(f"A = \n{A}\nand det(A) = {np.linalg.det(A)}.")

A = 
[[ 1  0  0 -1]
 [ 1  1  0  4]
 [ 1  0 -2  0]
 [ 1  1  1  0]]
and det(A) = 7.000000000000001.

Inverse¶

An inverse matrix can also be computed. To do so, use inv.

In [4]:
A = np.array([
    [1, 0, 0, -1],
    [1, 2, 0, 5],
    [1, 1, -2, 0],
    [1, 1, 1, 0]
])
A_inv = np.linalg.inv(A)
print(f"The inverse of A is A^-1 = \n{A_inv}")

The inverse of A is A^-1 = 
[[ 1.25        0.25       -0.16666667 -0.33333333]
 [-1.25       -0.25        0.5         1.        ]
 [ 0.          0.         -0.33333333  0.33333333]
 [ 0.25        0.25       -0.16666667 -0.33333333]]

Let us double check that the inverse $A^{-1}$ computed above is correct by multiplyting it to $A$. This should result in an identity matrix.

In [5]:
# Compute A*inv(A)
print(A@A_inv)

[[1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 1.00000000e+00 1.11022302e-16 2.22044605e-16]
 [0.00000000e+00 0.00000000e+00 1.00000000e+00 1.11022302e-16]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00]]
In [6]:
# Compute inv(A)*A
print(A_inv@A)

[[ 1.00000000e+00  5.55111512e-17  0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  1.00000000e+00  0.00000000e+00  0.00000000e+00]
 [-5.55111512e-17 -5.55111512e-17  1.00000000e+00  0.00000000e+00]
 [ 5.55111512e-17  5.55111512e-17  0.00000000e+00  1.00000000e+00]]

Eigenvalues and eigenvectors¶

To get eigenvalues of a matrix $A$ and the associated eigenvectors, we use eig. This function eig will return both the eigenvalues and eigenvectors.

The eig function returns two outputs. The first is the list of all eigenvalues. The second is the matrix whose $i^{\text{th}}$ column is a normalized eigenvector associated to the $i^{\text{th}}$ eigenvalue from the first output.

In [7]:
A = np.array([[1, 1, 4], [0, 2, 2], [0, 0, 2]])
eigval, eigvec = np.linalg.eig(A)
In [8]:
# Print eigenvalues λ_1, λ_2, λ_3
print(eigval)

[1. 2. 2.]
In [9]:
# Print a matrix whose columns are eigenvectors 
# corresponding to the eigenvalues above.
print(eigvec)

[[ 1.00000000e+00  7.07106781e-01 -7.07106781e-01]
 [ 0.00000000e+00  7.07106781e-01 -7.07106781e-01]
 [ 0.00000000e+00  0.00000000e+00  1.57009246e-16]]

Solving a linear system¶

Consider a linear system $Ax = b$. This linear system can be solved with the solve function.

The solve function works only when $A$ is a square matrix and the system has a unique solution. Here are different results based on the number of solutions of $Ax = b$:

  1. $Ax = b$ has a unique solution.
    • solve will properly return the solution.
  2. $Ax = b$ has no solution.
    • solve fails.
  3. $Ax = b$ has multiple solutions.
    • solve fails.

Unique solution of $Ax = b$¶

The uniqueness of a solution is guaranteed if $\mathrm{rank}(A) = n$, where $n$ is the number of variables.

In [10]:
A = np.array([[1, 2, 1], [0, 1, 0], [0, 1, 2]])
b = np.array([[1, 1, 0]]).T
Ab = np.hstack((A,b))

print(f"rank(A) = {np.linalg.matrix_rank(A)}, \
    rank([A b]) = {np.linalg.matrix_rank(Ab)}")

rank(A) = 3,     rank([A b]) = 3
In [11]:
x = np.linalg.solve(A,b)
print(f"The unique solution of Ax = b is \nx = \n{x}.")

The unique solution of Ax = b is 
x = 
[[-0.5]
 [ 1. ]
 [-0.5]].
In [12]:
# Double check that Ax = b.
A@x == b
Out[12]:
array([[ True],
       [ True],
       [ True]])

No solution of $Ax = b$¶

This case is characterized by $\mathrm{rank}(A) < \mathrm{rank}([A\; b])$.

In [13]:
A = np.array([[1, 1, 1], [1, 1, 0], [1, 1, 0], [0, 1, 0]])
b = np.array([[1, 0, 1, 0]]).T
Ab = np.hstack((A,b))
In [14]:
np.linalg.matrix_rank(A) == np.linalg.matrix_rank(Ab)
Out[14]:
np.False_
In [15]:
# Try solving Ax = b. This should produce an error.
x = np.linalg.solve(A,b)

---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
Cell In[15], line 2
      1 # Try solving Ax = b. This should produce an error.
----> 2 x = np.linalg.solve(A,b)

File /usr/lib/python3.13/site-packages/numpy/linalg/_linalg.py:457, in solve(a, b)
    384 """
    385 Solve a linear matrix equation, or system of linear scalar equations.
    386 
   (...)    454 
    455 """
    456 a, _ = _makearray(a)
--> 457 _assert_stacked_square(a)
    458 b, wrap = _makearray(b)
    459 t, result_t = _commonType(a, b)

File /usr/lib/python3.13/site-packages/numpy/linalg/_linalg.py:264, in _assert_stacked_square(*arrays)
    261     raise LinAlgError('%d-dimensional array given. Array must be '
    262             'at least two-dimensional' % a.ndim)
    263 if m != n:
--> 264     raise LinAlgError('Last 2 dimensions of the array must be square')

LinAlgError: Last 2 dimensions of the array must be square

Multiple solutions of $Ax = b$¶

This case is characterized by $\mathrm{rank}(A) = \mathrm{rank}([A\; b]) < n$.

In [16]:
A = np.array([[1, 1, 1], [1, 0, 0]])
b = np.array([[0, 2]]).T
Ab = np.hstack((A,b))
In [17]:
(np.linalg.matrix_rank(A) == np.linalg.matrix_rank(Ab)) \
and np.linalg.matrix_rank(A) < np.shape(b)[0]
Out[17]:
np.False_
In [18]:
# Try solving Ax = b. This should produce an error.
x = np.linalg.solve(A,b)

---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
Cell In[18], line 2
      1 # Try solving Ax = b. This should produce an error.
----> 2 x = np.linalg.solve(A,b)

File /usr/lib/python3.13/site-packages/numpy/linalg/_linalg.py:457, in solve(a, b)
    384 """
    385 Solve a linear matrix equation, or system of linear scalar equations.
    386 
   (...)    454 
    455 """
    456 a, _ = _makearray(a)
--> 457 _assert_stacked_square(a)
    458 b, wrap = _makearray(b)
    459 t, result_t = _commonType(a, b)

File /usr/lib/python3.13/site-packages/numpy/linalg/_linalg.py:264, in _assert_stacked_square(*arrays)
    261     raise LinAlgError('%d-dimensional array given. Array must be '
    262             'at least two-dimensional' % a.ndim)
    263 if m != n:
--> 264     raise LinAlgError('Last 2 dimensions of the array must be square')

LinAlgError: Last 2 dimensions of the array must be square