Linear Algebra: an Overview [1]

Concepts of Linear Algebra

Useful Libraries

import matplotlib.pyplot as plt
import numpy             as np

from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from scipy.linalg               import eig

Vector Space

x_column   = np.array([[-1], [0], [0.4], [10], [np.pi]]) # Column vector in R^5
x_row      = np.array([[-1, 0, 0.4, 10, np.pi]])         # Corresponding row vector
x_column_T = x_column.T                                  # Corresponding row vector using transpose

# Displaying the results
print('Column vector:')
print(x_column)
print('\nRow vector:')
print(x_row)
print('\nColumn vector transposed:')
print(x_column_T)

Column vector:
[[-1.        ]
[ 0.        ]
[ 0.4       ]
[10.        ]
[ 3.14159265]]

Row vector:
[[-1.          0.          0.4        10.          3.14159265]]

Column vector transposed:
[[-1.          0.          0.4        10.          3.14159265]]

Basic Vector Operations

# Define column vectors in R^5
x = np.array([-1, 0, 0.4, 10, np.pi])
y = np.array([0, 10, -0.4, -1, 3])

lamda = 0.11      # Real scalar in R
c     = x + y     # Sum of two vectors in R^5
d     = lamda * x # Scalar multiple of a vector in R^5

# Displaying the results
print('Vector x:\t', x)
print('Vector y:\t', y)
print('Scalar lamda:\t', lamda)
print('x + y:\t\t', c)
print('lamda * x:\t', d)

Vector x:	 [-1.          0.          0.4        10.          3.14159265]
Vector y:	 [ 0.  10.  -0.4 -1.   3. ]
Scalar lamda:	 0.11
x + y:		 [-1.         10.          0.          9.          6.14159265]
lamda * x:	 [-0.11        0.          0.044       1.1         0.34557519]

Important Vectors

f = np.array([0, 0, 1, 0, 0]) # Standard basis 𝑒3 in R^5
g = np.ones(5)                # Column vector of ones in R^5
h = np.zeros(5)               # Column vector of zeros in R^5

# Displaying the results
print('Standard basis:\t', f)
print('All ones:\t', g)
print('All zeros:\t', h)

Standard basis:	 [0 0 1 0 0]
All ones:	 [1. 1. 1. 1. 1.]
All zeros:	 [0. 0. 0. 0. 0.]

Matrix Space

A = np.array([[4, -1, 0], [-1, np.pi, 10]]) # Define a real-valued matrix A in R^{2x3}
B = A.T                                     # Transpose of A to get matrix B in R^{3x2}

# Displaying the matrices
print('Matrix A:')
print(A)
print('\nMatrix B (transpose of A):')
print(B)

Matrix A:
[[ 4.         -1.          0.        ]
[-1.          3.14159265 10.        ]]

Matrix B (transpose of A):
[[ 4.         -1.        ]
[-1.          3.14159265]
[ 0.         10.        ]]

Important Matrices

I = np.eye(3, 5)     # Identity matrix in R^{3x5} [notice that this matrix can be used to retrieve the standard bases]
J = np.zeros((5, 3)) # Zeros matrix in R^{5x3}
K = np.ones((6, 6))  # Ones matrix in R^{6x6}

# Displaying the matrices
print('Identity matrix:')
print(I)
print('\nZeros matrix:')
print(J)
print('\nOnes matrix:')
print(K)

Identity matrix:
[[1. 0. 0. 0. 0.]
[0. 1. 0. 0. 0.]
[0. 0. 1. 0. 0.]]

Zeros matrix:
[[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]]

Ones matrix:
[[1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1.]]

Dot Product

# Define column vectors x and y in R^3
x = np.array([-1, 0.1, np.pi])
y = np.array([0, 2, 0.3])

z_dot_product             = np.dot(x, y) # Calculate the dot product of x and y using np.dot
z_dot_product_alternative = x @ y        # Alternatively, use the @ operator for dot product

# Display the results
print('Vector x:\t\t\t', x)
print('Vector y:\t\t\t', y)
print('Dot product (using np.dot):\t', z_dot_product)
print('Dot product (using @ operator):\t', z_dot_product_alternative)

Vector x:			 [-1.          0.1         3.14159265]
Vector y:			 [0.  2.  0.3]
Dot product (using np.dot):	 1.142477796076938
Dot product (using @ operator):	 1.142477796076938

Vector Norms

Also known as $l_p$ norms. Notice how the values decrease as $p$ increases.

x = np.random.rand(5) # Generate a random vector x in R^5

# Display the results
print(f'Random vector x: {x}\n')
[print(f'l_{i + 1} norm of x:\t\t{np.linalg.norm(x, i + 1)}') for i in range(10)]
print(f'l_inf norm of x:\t{np.linalg.norm(x, np.inf)}')

Random vector x: [0.43708824 0.71746819 0.88530565 0.21570394 0.79721364]

l_1 norm of x:		3.05277964873668
l_2 norm of x:		1.4736521279127217
l_3 norm of x:		1.1848574147129909
l_4 norm of x:		1.0722496582133314
l_5 norm of x:		1.0140743389917317
l_6 norm of x:		0.979224827151816
l_7 norm of x:		0.9563783081122994
l_8 norm of x:		0.9404785071815713
l_9 norm of x:		0.9289404454462779
l_10 norm of x:		0.9203081381556143
l_inf norm of x:	0.8853056517077345

Matrix Norms

A = np.random.rand(5, 3) # Generate a random matrix A in R^{5x3}

# Calculate the spectral norm, 1-norm, inf-norm, and frobenius norm of A
b_spectral_norm  = np.linalg.norm(A, 2)
c_1_norm         = np.linalg.norm(A, 1)
d_infinity_norm  = np.linalg.norm(A, np.inf)
e_frobenius_norm = np.linalg.norm(A, 'fro')

# Display the results
print('Random matrix A:')
print(A)
print('\nSpectral norm of A:\t', b_spectral_norm)
print('1-norm of A:\t\t', c_1_norm)
print('inf-norm of A:\t\t', d_infinity_norm)
print('Frobenius norm of A:\t', e_frobenius_norm)

Random matrix A:
[[0.3647405  0.24450127 0.91209795]
[0.67139985 0.224112   0.13627266]
[0.63259974 0.57414245 0.18500187]
[0.47487561 0.81229876 0.69132714]
[0.41225556 0.42721026 0.32375644]]

Spectral norm of A:	 1.8782801278395589
1-norm of A:		 2.555871261143026
inf-norm of A:		 1.9785015058661313
Frobenius norm of A:	 2.03197240180481

Solving Linear Equations

The function np.linalg.solve solves linear equations using LAPACK [Linear Algebra PACKage] routines [2] that use $LU$ decomposition with partial pivoting and row interchanges [3], which is much faster than computing inverses.

A = np.array([[2, 1], [5, 3]]) # Coefficient matrix
b = np.array([1, 2])           # Right-hand side
x = np.linalg.solve(A, b)      # Solve for x in Ax = b

print('Solution to Ax = b:', x) # Display the results

Solution to Ax = b: [ 1. -1.]

Inverse and Pseudo-inverse

A            = np.array([[2, 1], [5, 3]]) # Create a square matrix
A_inv        = np.linalg.inv(A)           # Calculate the inverse
A_pseudo_inv = np.linalg.pinv(A)          # Calculate the pseudo-inverse

# Display the results
print('Original matrix:')
print(A)
print("\nInverse:")
print(A_inv)
print("\nPseudo-inverse:")
print(A_pseudo_inv)

Original matrix:
[[2 1]
[5 3]]

Inverse:
[[ 3. -1.]
[-5.  2.]]

Pseudo-inverse:
[[ 3. -1.]
[-5.  2.]]

Eigen/Spectral Decomposition

A               = np.random.rand(5, 5)              # Generate a random square matrix A in R^{5x5} with random entries
L, W, V         = eig(A, left = True, right = True) # Eigen decomposition Av = Lv and (w.T)A = L(w.T)

# Display the results
print('Random square matrix A:')
print(A)
print('\nEigen decomposition Av = Lv and (w.T)A = L(w.T):')
print('Eigenvalues:')
print(L)
print('\nLeft eigenvectors:')
print(W)
print('\nRight eigenvectors:')
print(V)

Random square matrix A:
[[0.8932966  0.45282265 0.36566497 0.00881553 0.44327841]
[0.09954157 0.94588855 0.94002498 0.58760723 0.52970619]
[0.96261601 0.58199552 0.79279289 0.05939403 0.2199229 ]
[0.56420436 0.82604279 0.0547052  0.97303312 0.80147121]
[0.3208708  0.52758568 0.93229637 0.50868777 0.42811888]]

Eigen decomposition Av = Lv and (w.T)A = L(w.T):
Eigenvalues:
[ 2.74036   +0.j          1.00770919+0.j         -0.02089465+0.34428413j
-0.02089465-0.34428413j  0.32685016+0.j        ]

Left eigenvectors:
[[-0.4596902 +0.j          0.64178006+0.j         -0.63170669+0.j
-0.63170669-0.j          0.55032634+0.j        ]
[-0.53394184+0.j          0.02294653+0.j         -0.14640081-0.03094337j
-0.14640081+0.03094337j -0.46617829+0.j        ]
[-0.52752528+0.j          0.46617529+0.j          0.54456807+0.26452388j
0.54456807-0.26452388j -0.48280619+0.j        ]
[-0.30276895+0.j         -0.59388954+0.j          0.1139737 +0.17332598j
0.1139737 -0.17332598j  0.07415464+0.j        ]
[-0.36556677+0.j         -0.13254291+0.j          0.0110926 -0.41094253j
0.0110926 +0.41094253j  0.49113949+0.j        ]]

Right eigenvectors:
[[-0.30888212+0.j          0.31026575+0.j         -0.32112105-0.06187647j
-0.32112105+0.06187647j  0.34953053+0.j        ]
[-0.52185555+0.j         -0.25549229+0.j          0.13859625-0.41128887j
0.13859625+0.41128887j -0.80745345+0.j        ]
[-0.37444264+0.j          0.45758844+0.j          0.0093209 +0.35691348j
0.0093209 -0.35691348j  0.10592801+0.j        ]
[-0.55098567+0.j         -0.79287212+0.j         -0.46160695+0.19740289j
-0.46160695-0.19740289j  0.36639099+0.j        ]
[-0.43412685+0.j         -0.02062692+0.j          0.57022989+0.j
0.57022989-0.j          0.28352114+0.j        ]]

Singular Value Decomposition

A               = np.random.rand(5, 5) # Generate a random square matrix A in R^{5x5} with random entries
U, S, V_T       = np.linalg.svd(A)     # Perform SVD
A_reconstructed = U @ np.diag(S) @ V_T # Reconstruct the original matrix

# Display the results
print('Random square matrix A:')
print(A)
print('\nSingular value decomposition A = US(V.T):')
print('Singular values:', S)
print('\nLeft singular vectors:')
print(U)
print('\nRight singular vectors:')
print(V_T)
print('\nReconstructed matrix:')
print(A_reconstructed)

Random square matrix A:
[[0.54507266 0.10888492 0.66770098 0.84383544 0.82635158]
[0.20988294 0.48496274 0.28872736 0.87402098 0.66010815]
[0.43101926 0.1280712  0.84323557 0.96338147 0.34364795]
[0.496449   0.73858452 0.68138546 0.16802926 0.12675549]
[0.3725026  0.17783032 0.62913034 0.90750959 0.55495563]]

Singular value decomposition A = US(V.T):
Singular values: [2.75705868 0.82044771 0.53275276 0.33354778 0.00935693]

Left singular vectors:
[[-0.51833673 -0.2350247   0.06226531 -0.78827102 -0.22548227]
[-0.42181842 -0.14829821 -0.80381593  0.32216302 -0.2239817 ]
[-0.49130256 -0.04158911  0.58065137  0.49938775 -0.41273235]
[-0.30309344  0.94394051 -0.06864311 -0.10097988  0.04692412]
[-0.46920402 -0.17325616  0.09019474  0.12350999  0.85231454]]

Right singular vectors:
[[-0.32936357 -0.22894899 -0.50194151 -0.63695291 -0.42596679]
[ 0.2765849   0.68686099  0.36489009 -0.44686102 -0.34480898]
[ 0.21590459 -0.64445595  0.58017283 -0.0381096  -0.44722353]
[-0.45248776  0.24507754 -0.00993224  0.57750659 -0.63370286]
[-0.75077325  0.02044192  0.5274569  -0.24421305  0.31316268]]

Reconstructed matrix:
[[0.54507266 0.10888492 0.66770098 0.84383544 0.82635158]
[0.20988294 0.48496274 0.28872736 0.87402098 0.66010815]
[0.43101926 0.1280712  0.84323557 0.96338147 0.34364795]
[0.496449   0.73858452 0.68138546 0.16802926 0.12675549]
[0.3725026  0.17783032 0.62913034 0.90750959 0.55495563]]

Symmetric Matrix, Trace, and Determinant

A       = np.random.rand(5, 5) # Generate a random square matrix A in R^{5x5} with random entries
B       = 0.5 * (A + A.T)      # Make a symmetric matrix B
B_trace = np.trace(B)          # Compute the trace of B
B_det   = np.linalg.det(B)     # Compute the determinant of B

# Display the results
print('Random square matrix A:')
print(A)
print('\nSymmetric matrix B:')
print(B)
print('\nTrace of B:\t\t', B_trace)
print('Determinant of B:\t', B_det)

Random square matrix A:
[[0.65668062 0.40152043 0.02255261 0.80307654 0.45529925]
[0.16087064 0.43485765 0.26740823 0.0314879  0.07877773]
[0.64471248 0.56871998 0.47943541 0.88634095 0.4488305 ]
[0.54469606 0.28753323 0.18928733 0.68736046 0.22221986]
[0.0493979  0.03626296 0.14480577 0.95558066 0.96123289]]

Symmetric matrix B:
[[0.65668062 0.28119553 0.33363255 0.6738863  0.25234858]
[0.28119553 0.43485765 0.4180641  0.15951057 0.05752034]
[0.33363255 0.4180641  0.47943541 0.53781414 0.29681813]
[0.6738863  0.15951057 0.53781414 0.68736046 0.58890026]
[0.25234858 0.05752034 0.29681813 0.58890026 0.96123289]]

Trace of B:		 3.2195670305846287
Determinant of B:	 -0.009958184122829259

Gradient Vector: Numerical Computation

f       = lambda x, y: np.exp(-x**2 - y**2)                       # Define the bivariate function f(x, y)
grad    = lambda x, y, h: np.array([(f(x + h, y) - f(x, y)) / h,
                                    (f(x, y + h) - f(x, y)) / h]) # Define the gradient function
f_value = f(2, 3)                                                 # Compute the function value at (2, 3)
g       = grad(2, 3, 0.002)                                       # Compute the numerical gradient at (2, 3)

# Display the results
print('Function Value at (2, 3):\t', f_value)
print('Numerical Gradient at (2, 3):\t', g)

Function Value at (2, 3):	 2.2603294069810542e-06
Numerical Gradient at (2, 3):	 [-9.00973323e-06 -1.34853958e-05]

Hessian Matrix: Numerical Computation

f       = lambda x, y: np.exp(-x**2 - y**2)                                                           # Define the bivariate function f(x, y)
grad    = lambda x, y, h: np.array([(f(x + h, y) - f(x, y)) / h,
                                    (f(x, y + h) - f(x, y)) / h])                                     # Define the gradient function
hess    = lambda x, y, h: np.array([[(f(x + h, y) - 2 * f(x, y) + f(x - h, y)) / h**2,
                                    (f(x + h, y + h) - f(x + h, y) - f(x, y + h) + f(x, y)) / h**2],
                                    [(f(x + h, y + h) - f(x + h, y) - f(x, y + h) + f(x, y)) / h**2,
                                    (f(x, y + h) - 2 * f(x, y) + f(x, y - h)) / h**2]])              # Define the hessian function
f_value = f(2, 3)                                                                                     # Compute the function value at (2, 3)
g       = grad(2, 3, 0.002)                                                                           # Compute the numerical gradient at (2, 3)
H       = hess(2, 3, 0.002)                                                                           # Compute the numerical hessian at (2, 3)

# Display the results
print('Function value at (2, 3):\t', f_value)
print('Numerical gradient at (2, 3):\t', g)
print('\nNumerical hessian at (2, 3):')
print(H)

Function value at (2, 3):	 2.2603294069810542e-06
Numerical gradient at (2, 3):	 [-9.00973323e-06 -1.34853958e-05]

Numerical hessian at (2, 3):
[[3.16446690e-05 5.37531470e-05]
[5.37531470e-05 7.68518599e-05]]

Visualizations of Linear Algebra

Plotting Vectors

The function np.random.rand generates an array of the specified shape, filled with random samples from a uniform distribution over $[0, 1)$ [4]. Consequently, each execution of the cell below will plot a distinct vector on the $1 \times 1$ grid.

vector = np.random.rand(2) # Create a 2D vector with random components

plt.figure()                                                          # Plot the vector
plt.quiver(0, 0, vector[0], vector[1], angles = 'xy',
        scale_units = 'xy', scale = 1, color = 'r')                # Use quiver to represent the vector as an arrow
plt.text(vector[0], vector[1], f'({vector[0]:.2f}, {vector[1]:.2f})',
        fontsize = 10, color = 'b', va = 'bottom', ha = 'left')      # Add text annotation next to the arrowhead with the vector coordinates
plt.xlim(0, 1)                                                        # Set the x-axis limits to ensure proper visualization
plt.ylim(0, 1)                                                        # Set the y-axis limits to ensure proper visualization
plt.title('Vector Plot')                                              # Set the title of the plot
plt.grid(True)                                                        # Display the grid for reference
plt.show()                                                            # Show the plot

Contour Plot

x, y         = np.meshgrid(np.arange(-2, 2.1, 0.1),
                            np.arange(-2, 2.1, 0.1)) # Create a grid of points (meshgrid)
z            = -3 * y / (x**2 + y**2 + 1)           # Define the continuous function z = f(x, y)
contour_plot = plt.contour(x, y, z, levels = 12)    # Plot contours of z = f(x, y) [change 'levels' to see different contours]plt.clabel(contour_plot, inline = True, fontsize = 8) # Label contours of z = f(x, y)

# Display the plot
plt.xlabel('x', fontsize = 15)
plt.ylabel('y', fontsize = 15)
plt.title('Contour Plot of z = f(x, y)', fontsize = 15)
plt.show()

Unit Simplex

vertices = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) # Define the vertices of the unit simplex in 3D

# Create a 3D plot
fig = plt.figure(figsize = (14, 14))
ax  = fig.add_subplot(111, projection = '3d')

ax.scatter(vertices[0], vertices[1], vertices[2],
        s = 200, label = 'vertices')                        # Plot the unit simplex vertices
ax.add_collection3d(Poly3DCollection([vertices], alpha = 0.5)) # Plot the area between the vertices

# Set labels and title
ax.set_xlabel('x', fontsize = 15)
ax.set_ylabel('y', fontsize = 15)
ax.set_zlabel('z', fontsize = 15)
ax.set_title('3D Unit Simplex', fontsize = 15)
ax.legend()

plt.show() # Show the plot

Geometry of Vector Norms

# Generating values for both the axes
x = np.linspace(-1, 1, 2001)
y = np.linspace(-1, 1, 2001)

# Generating the grid and coordinates
X, Y    = np.meshgrid(x, y)
vectors = np.vstack([X.ravel(), Y.ravel()]).T

# Calculate norms
norm_0   = np.linalg.norm(vectors, ord = 0, axis = 1)
norm_05  = np.linalg.norm(vectors, ord = 0.5, axis = 1)
norm_1   = np.linalg.norm(vectors, ord = 1, axis = 1)
norm_2   = np.linalg.norm(vectors, ord = 2, axis = 1)
norm_inf = np.linalg.norm(vectors, ord = np.inf, axis = 1)

# Rescale vectors to have unit norms [except zero norm; it is directly indexed]
vectors_0   = vectors[np.where(norm_0 == 1)]
vectors_05  = vectors[np.where(norm_05 != 0)] / norm_05[norm_05 != 0].reshape(-1, 1)
vectors_1   = vectors[np.where(norm_1 != 0)] / norm_1[norm_1 != 0].reshape(-1, 1)
vectors_2   = vectors[np.where(norm_2 != 0)] / norm_2[norm_2 != 0].reshape(-1, 1)
vectors_inf = vectors[np.where(norm_inf != 0)] / norm_inf[norm_inf != 0].reshape(-1, 1)

# Plotting
plt.figure(figsize = (10, 8))
plt.scatter(vectors_0[:, 0], vectors_0[:, 1], label = r'$||\mathbf{x}||_0 = 1$', c = 'k')
plt.scatter(vectors_05[:, 0], vectors_05[:, 1], label = r'$||\mathbf{x}||_{0.5} = 1$', c = 'red')
plt.scatter(vectors_1[:, 0], vectors_1[:, 1], label = r'$||\mathbf{x}||_1 = 1$', c = 'blue')
plt.scatter(vectors_2[:, 0], vectors_2[:, 1], label = r'$||\mathbf{x}||_2 = 1$', c = 'green')
plt.scatter(vectors_inf[:, 0], vectors_inf[:, 1], label = r'$||\mathbf{x}||_{\infty} = 1$', c = 'purple')
plt.title('Unit Vector Norms', fontsize = 15)
plt.xlabel(r'$x_1$', fontsize = 15)
plt.ylabel(r'$x_2$', fontsize = 15)
plt.legend(loc = 'upper right')
plt.grid(True)
plt.xlim([-1.5, 1.5])
plt.ylim([-1.5, 1.5])

plt.gca().set_aspect('equal', adjustable = 'box') # Set equal aspect ratio to display squares as squares
plt.show()                                        # Show the plot

Norm Balls

This cell may take longer to execute. Note that the $l_{0.5}$ ball is incomplete. Due to memory and time constraints, the arguments of np.linspace are designed to produce insufficient points. For more comprehensive visualizations, refer to resources available online [5].

# Generate vectors
x           = np.linspace(-1, 1, 101)
y           = np.linspace(-1, 1, 101)
z           = np.linspace(-1, 1, 101)
X, Y, Z     = np.meshgrid(x, y, z)
vectors     = np.vstack([X.ravel(), Y.ravel(), Z.ravel()]).T
vectors_new = vectors[(vectors != 0).all(axis = 1)]

# Calculate norms
norm_0   = np.linalg.norm(vectors, ord = 0, axis = 1)
norm_05  = np.linalg.norm(vectors_new, ord = 0.5, axis = 1)
norm_1   = np.linalg.norm(vectors_new, ord = 1, axis = 1)
norm_2   = np.linalg.norm(vectors_new, ord = 2, axis = 1)
norm_inf = np.linalg.norm(vectors_new, ord = np.inf, axis = 1)

# Rescale vectors to have unit norms [except zero norm; it is directly indexed]
vectors_0   =  vectors[np.where(norm_0 == 1)]
vectors_05  = vectors_new / norm_05.reshape(-1, 1)
vectors_1   = vectors_new / norm_1.reshape(-1, 1)
vectors_2   = vectors_new / norm_2.reshape(-1, 1)
vectors_inf = vectors_new / norm_inf.reshape(-1, 1)

fig = plt.figure(figsize = (20, 5)) # Plot the figure

# l_0 norm subplot
ax1 = fig.add_subplot(151, projection = '3d')
ax1.scatter(vectors_0[:, 0], vectors_0[:, 1], vectors_0[:, 2],
            c = norm_0[norm_0 == 1], alpha = 0.3)
ax1.set_title('$l_0$ Norm Ball')

# l_0.5 norm subplot
ax2 = fig.add_subplot(152, projection = '3d')
ax2.scatter(vectors_05[:, 0], vectors_05[:, 1], vectors_05[:, 2],
            c = norm_05, alpha = 0.3)
ax2.set_title('$l_{0.5}$ Norm Ball')

# l_1 norm subplot
ax3 = fig.add_subplot(153, projection = '3d')
ax3.scatter(vectors_1[:, 0], vectors_1[:, 1], vectors_1[:, 2],
            c = norm_1, alpha = 0.3)
ax3.set_title('$l_1$ Norm Ball')

# l_2 norm subplot
ax4 = fig.add_subplot(154, projection = '3d')
ax4.scatter(vectors_2[:, 0], vectors_2[:, 1], vectors_2[:, 2],
            c = norm_2, alpha = 0.3)
ax4.set_title('$l_2$ Norm Ball')

# l_inf norm subplot
ax5 = fig.add_subplot(155, projection = '3d')
ax5.scatter(vectors_inf[:, 0], vectors_inf[:, 1], vectors_inf[:, 2],
            c = norm_inf, alpha = 0.3)
ax5.set_title('$l_{\infty}$ Norm Ball')

plt.show() # Show the plot

References

1. Dr. Hassan Mohy Ud Din, "Lecture 0 – Linear Algebra, A Review", 2024
2. numpy.linalg.solve
3. gesv: factor and solve
4. numpy.random.rand
5. Graphing the $p$-Norm Unit Ball in 3 Dimensions