3.4. Chapter 4: Geometric Aspects of Linear Algebra#

Before Class:

Read Chapter 4 in Savov (2020) and take notes
Watch the following videos and take notes:
Compile a list of questions to bring to class

During and After Class:

Take notes (on paper or a tablet computer)
Complete this notebook, submit you answer via Gradescope

import sys
if "google.colab" in sys.modules:
    # If running on colab, install ipympl
    !pip install -q ipympl

    # Enable 3rd party output/widgets
    from google.colab import output
    output.enable_custom_widget_manager()

# % means a "magic command" that is special
# for Jupyter notebooks.
# This specific command enables interactive plots
%matplotlib ipympl

# Here are our standard imports
import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as la

3.4.1. Visualize Linear Systems in 3D#

The code below visualizes the rows on a N x 3 matrix in three dimensions.

from mpl_toolkits.mplot3d import Axes3D

def visualize_linear_system(A, b, xlim=(-2,2), ylim=(-2,2), return_ax = False):
    """ Visualize rows of linear system as planes in 3D
    
    Arguments:
        A: N x 3 matrix
        B: N x 1 vector
        xlim: tuple of bounds for x dimension
        ylim: tuple of bounds for y dimension
    """

    assert A.shape[1] == 3, "Matrix A must have 3 columns"
    assert A.shape[0] == len(b), "Matrix A and vector b must have the same number of rows"
    
    N = len(b)

    for i in range(N):
        assert np.abs(A[i,2]) > 1E-8, "This function does not support zeros in A[:,2]"

    # Create empty plot
    fig = plt.figure()
    ax = Axes3D(fig, auto_add_to_figure=False)
    fig.add_axes(ax)

    # Create a grid of points for plotting
    X, Y = np.meshgrid(np.linspace(xlim[0],xlim[1],5), np.linspace(ylim[0],ylim[1],5))

    # Loop over equations
    for i in range(N):

        # Define function to calculate z
        func_z = lambda x,y : (b[i] - A[i,0]*x - A[i,1]*y)/A[i,2]
        
        # Calculate values of z over meshgrid
        Z = func_z(X, Y)


        s = "{0:.2f}x + {1:.2f}y + {2:.2f}z = {3:.2f}".format(A[i,0], A[i,1], A[i,2], b[i])
        surf = ax.plot_surface(X, Y, Z, alpha=0.5, label=s)

        # Needed for a legend
        # https://stackoverflow.com/questions/27449109/adding-legend-to-a-surface-plot
        surf._edgecolors2d = surf._edgecolor3d
        surf._facecolors2d = surf._facecolor3d

    # If the linear system is square and full rank
    if A.shape[0] == 3 and np.linalg.matrix_rank(A) == 3:
        x = np.linalg.solve(A,b)
        ax.scatter(x[0], x[1], x[2], marker='x', label="solution", color="k")


    ax.set_xlabel('x',fontsize=18)
    ax.set_ylabel('y',fontsize=18)
    ax.set_zlabel('z',fontsize=18)

    if return_ax:
        return ax
    else:
        plt.legend()
        plt.show()

Let’s visualize the equations in the rows of the linear system:

\[\begin{split} \begin{bmatrix} 1 & 1 & 1 \\ 2 & -1 & 1 \end{bmatrix} {\bf x} = \begin{bmatrix} 0 \\ 2 \end{bmatrix} \end{split}\]

A = np.array([[1,1,1],[2,-1,1]])
b = np.array([0, 2])

visualize_linear_system(A, b)

Geometrically describe the solution of this linear system in a few words.

Your Answer:

Let’s visualize the equations in the rows of the linear system:

\[\begin{split} \begin{bmatrix} 1 & 2 & 1 \\ -1 & 0 & -1 \\ 0 & 3 & 1 \end{bmatrix} {\bf x} = \begin{bmatrix} 0 \\ 2 \\ -1 \end{bmatrix} \end{split}\]

# Add your solution here

3.4.2. Projections#

3.4.2.1. Vector Projected onto a Vector#

Let’s start by projecting vector ${\bf u}$ onto the line defined by vector ${\bf v}$.

\[ \Pi_{\bf{v}} {\bf u}= \frac{ {\bf u} \cdot {\bf v} }{|| {\bf v} ||^2 } {\bf v} \]

Now let’s perform an example calculation in Python and visualize the result in 3D. We’ll reuse some plotting code from Chapter 2: Intro to Linear Algebra.

# define vectors
u = np.array([1,1,1])
v = np.array([2,3,1])

# calculate the l2-norm of v
norm_v = np.linalg.norm(v, ord=2)

proj_u_onto_v = np.dot(u, v)/ norm_v**2 * v

# Create 3D figure and axes
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Set labels and limits
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_xlim([-1,2])
ax.set_ylim([-1,3])
ax.set_zlim([0,2])

# Plot vectors
def plot_vector(vector, color, label, p = [0,0,0], linestyle="-"):
    ''' Plot a 3D vector
        Arguments:
            vector: np array
            color: 'r', 'b', 'g', etc.
            label: string for legend
            p: point to start vector
        
    '''
    assert len(vector) == 3, "Must be a 3D vector"
    assert len(p) == 3, "p must have length 3"

    ax = plt.gca()
    ax.quiver(p[0], p[1], p[2], vector[0], vector[1], vector[2], color=color, arrow_length_ratio=0.1, label=label, linestyle=linestyle)

plot_vector(u, 'r', r'${\bf u}$')
plot_vector(v, 'b', r'${\bf v}$')
plot_vector(proj_u_onto_v, 'purple', r'$\Pi_{\bf v} {\bf u}$')

plt.legend(loc='best',ncol=3)
plt.show()

print("v =",v)
print("u =",u)
print("u projected onto v =",proj_u_onto_v)

Verify this calculation with pencil and paper. (It is excellent quiz practice.)

3.4.2.2. Vector Projected onto a Plane#

Projecting vector $\bf{u}$ onto a plane is a little more involved, requiring three steps:

Calculate the vector ${\bf n}$ orthagonal to the plane.
Project $\bf{u}$ onto ${\bf n}$.
Subtract this projection from $\bf{u}$.

Let’s do it in Python and visualize.

## Step 1: Orthagonal vector for plane
# Consider the plane defined by
A = np.array([[1, 1, 2]])
b = np.array([-1])

# The vector orthagonal to this plane is simply
n = A.copy().flatten()

# And we'll draw our vectors around the point
p = np.array([0, 0, b[0]/A[0,2]])

## Step 2: Projection

# calculate the l2-norm of v
norm_n = np.linalg.norm(n, ord=2)

# calculate the project
proj_u_onto_n = np.dot(u, n)/ norm_n**2 * n

## Step 3: Subtract
proj_u_onto_P = u - proj_u_onto_n

## Step 4: Plot

# Plane
ax = visualize_linear_system(A, b, return_ax=True)

# u
plot_vector(u, 'r', r'${\bf u}$', p=p)

# Orthagonal vector
plot_vector(n, 'b', r'${\bf n}$', p=p)

# Projection (step 2)
plot_vector(proj_u_onto_n, 'purple', r'$\Pi_{\bf n} {\bf u}$', p=p)

# Projection (step 3)
plot_vector(proj_u_onto_P, 'g', r'$\Pi_{P} {\bf u}$', p=p)

plt.legend()
plt.show()

print("u =",u)
print("n =",n)
print("u projected onto n =",proj_u_onto_n)
print("u projected onto P =",proj_u_onto_P)

Verify these calculations with pencil and paper (again, good quiz practice).

3.4.3. Visualizing Spaces of a Matrix#

We’ll consider the matrix

\[\begin{split} {\bf A} = \begin{bmatrix} 4 & -4 & 0 \\ 1 & 1 & -2 \\ 2 & -6 & 4 \end{bmatrix} \end{split}\]

Example 1 on pg. 248-250 of Savov (2020) walks through how to calculate the column space, row space, and null space of this matrix. Please read these few pages. In this notebook, we’ll perform these calculations numerically and then visualize in 3D.

# Define matrix
A = np.array([[4, -4, 0], [1, 1, -2], [2, -6, 4]])

# In order to plot the planes of the linear system
# encoded in A we also need to define b.
# Our choice is arbitrary and just informs plotting.
b = np.zeros(3)

3.4.3.1. Basis for the Column Space#

A basis for the range of ${\bf A}$ is also a basis for the column space of ${\bf A}$.

basis_column_space = la.orth(A)
print(basis_column_space)

3.4.3.2. Basis for the Row Space#

The basis of the row space can be calculated from the basis of the range of ${\bf A}^T$.

basis_row_space = la.orth(A.T)
print(basis_row_space)

3.4.3.3. Basis for the Null Space#

basis_null_space = la.null_space(A)
print(basis_null_space)

3.4.3.4. Analysis#

The row space basis, column space basis, and null space basis calculated above do not match the book exactly. Why does this make sense? Explain in a few sentences.

Your Answer:

Is the basis of the row space orthagonal to the basis of the null space?

print("first dot product =",np.dot(basis_row_space[:,0], basis_null_space)[0])
print("second dot product =",np.dot(basis_row_space[:,1], basis_null_space)[0])

Your Answer:

3.4.3.5. Visualize#

### Let's customize our function from earlier

xlim = [-1, 1]
ylim = [-1, 1]

colors = ["red","blue","green"]

assert A.shape[1] == 3, "Matrix A must have 3 columns"
assert A.shape[0] == len(b), "Matrix A and vector b must have the same number of rows"

N = len(b)

# Create empty plot
fig = plt.figure()
ax = Axes3D(fig, auto_add_to_figure=False)
fig.add_axes(ax)

### Plot rows (planes)

# Create a grid of points for plotting
X, Y = np.meshgrid(np.linspace(xlim[0],xlim[1],5), np.linspace(ylim[0],ylim[1],5))

# Loop over equations
for i in range(N):

    s = "row "+str(i)

    # Okay to invert for z
    if np.abs(A[i,2]) > 1E-8:

        # Define function to calculate z
        func_z = lambda x,y : (b[i] - A[i,0]*x - A[i,1]*y)/A[i,2]
        
        # Calculate values of z over meshgrid
        Z = func_z(X, Y)

        # s = "{0:.2f}x + {1:.2f}y + {2:.2f}z = {3:.2f}".format(A[i,0], A[i,1], A[i,2], b[i])
        
        surf = ax.plot_surface(X, Y, Z, alpha=0.5, label=s, color=colors[i])

    # Instead invert for x
    else:
        func_x = lambda z,y : (b[i] - A[i,2]*z - A[i,1]*y)/A[i,0]

        # Calculate values of x over meshgrid
        Z_ = X
        X_ = func_x(Z_, Y)

        surf = ax.plot_surface(X_, Y, Z_, alpha=0.5, label=s, color=colors[i])

    # Needed for a legend
    # https://stackoverflow.com/questions/27449109/adding-legend-to-a-surface-plot
    surf._edgecolors2d = surf._edgecolor3d
    surf._facecolors2d = surf._facecolor3d

    # Grab orthagonal vector for each plane
    n = A[i,:].copy().flatten()
    n = n / la.norm(n, ord=2)

    # This works out for this example
    p = [0, 0, 0]
    
    plot_vector(n, colors[i], None, p=p)

### Plot basis of column space
for i in range(basis_column_space.shape[1]):
    plot_vector(basis_column_space[:,i], 'black', "column space basis", p=p, linestyle='--')

### Plot basis of row space
for i in range(basis_row_space.shape[1]):
    plot_vector(basis_row_space[:,i], 'black', "row space basis", p=p, linestyle=':')

### Plot basis of null space
plot_vector(basis_null_space, 'grey', "null space basis", linestyle='-.')

### Finish formating

ax.set_xlabel('x',fontsize=18)
ax.set_ylabel('y',fontsize=18)
ax.set_zlabel('z',fontsize=18)

plt.legend()
plt.show()

Observations:

Basis of row space is orthagonal to the basis of the null space.
Because ${\bf b} = {\bf 0}$ for this example, the basis of the null space points along the line at the intersection of the planes. Thus the null space described the infinite number of solutions to this linear system.
Each plane corresponds to a row in the linear system. The rows are further visualized as the vectors orthagonal to each plane.
We can see these three vectors (each row) fall on a plane, thus these rows are not linearly independent.
Moreover, the basis of the row space also falls on the same plane as the three colored vectors.

3.4.4. Balance Chemical Formula#

Have you ever wondered how to systematically balance a chemical formula? For example,

$a$ CH$_4$ + $b$ O$_2$ $\rightarrow$ $c$ CO$_2$ + $d$ H$_2$O

We know that the number of C, H, and O atoms needs to be conserved (balanced), which gives the following system of linear equations:

\[a = c, \quad 4a = 2d, \quad 2b = 2c + d\]

Let’s write this in matrix notation: $$ \underbrace{\begin{bmatrix} 1 & 0 & -1 & 0 \\ 4 & 0 & 0 & -2 \\ 0 & 2 & -2 & -1 \end{bmatrix}}_{\bf A} \underbrace{\begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix}}_{\bf x} = {\bf 0} $$

There is an infinite number of choices for $\bf x$ because we have three equations (from three atom balances) but there are four coefficients. What should we do?

Let’s calculate the basis of the null space of $\bf A$. Please review pg. 364 and 365 for how to do this with pencil and paper. We’ll use Python here:

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

basis_null_space_chemical = la.null_space(A)
print("\nnull space =\n", basis_null_space_chemical)

We are free to choose any value of $\bf x$ that lines within the null space. Because the null space is dimension one, we’ll choose an multiple of the vector above. For example,

\[\begin{split}\bf{x} = \begin{bmatrix} 1 \\ 2 \\ 1 \\ 2 \end{bmatrix}\end{split}\]

which corresponds to:

1 CH$_4$ + 2 O$_2$ $\rightarrow$ 1 CO$_2$ + 2 H$_2$O

Now it’s your turn to try. Balance the following chemical formula,

$a$ Fe(OH)$_3$ + $b$ HCl $\rightarrow$ $c$ FeCl$_3$ + $d$ H$_2$ O

using the following steps:

On paper, write the atom balances to obtain a system of linear equations.
On paper, write this system in matrix form.
Using Python, calculate a basis for the null space. Hint: If you find the null space is empty, check the coefficients in your linear system. You made a mistake.
On paper, verify the null space calculation. This is good practice.
Interpret the basis of the null space to balance the chemical formula. By chemistry convention, avoid fractions for the unknown coefficients.

# Add your solution here

Typseset your answer here by replacing $a$, …, $d$ with whole numbers.

$a$ Fe(OH)$_3$ + $b$ HCl $\rightarrow$ $c$ FeCl$_3$ + $d$ H$_2$ O

3.4.5. More Practice with Linear Systems#

This question is from Prof. David Leighton.

In your laboratory, you measure UV absorbance ( units: [-]/(g/liter) ) to infer protein concentrations in a solution. For dilute solutions, the protein absorbance is proportional to the concentration, and the total absorbance is just the sum of that resulting from each species independently. Suppose you have calibrated the absorbance at three different wavelengths (200 nm, 300 nm, 400 nm) for three samples of pure protein species (BSA, BHb, Insulin):

Species	200 nm	300 nm	400 nm
BSA	0.1	0.2	0.1
BHb	0.1	0.3	0.4
Insulin	0.1	0.1	0.2

All entries in this table have units [-]/(g/liter). We are using [-] to represent dimensionless units for absorbance.

For a new sample, you measure the absorbance at 200 nm, 300 nm, and 400 nm to be 0.08 [-], 0.16 [-], and 0.20 [-], respectively. What is the concentration of each protein in your new sample?

3.4.5.1. Linear Model#

Develop a linear model on paper in canonical form (i.e., $\mathbf{A} \vec{x} = \vec{b}$) to calculate the concentration of each protein. Express your final answer in matrix form and turn in via Gradescope (pencil and paper submission).

Hint: What is unknown and what is known? What should be $\vec{x}$ and what should be $\vec{b}$?

3.4.5.2. Solve#

Solve the model in Python. Store your answers for protein concentration in the numpy array x.

# Set up the system of equations in matrix form

proteins = ['BSA','BHb','Insulin']

# Add your solution here

# Print solutions
print('The concentration of each protein in the sample is as follows:')
for i in range(len(x)):
    print("Concentration of {0} = {1:0.3f} g/liter.".format(proteins[i],x[i]))

3.4.6. Submission#

Please submit two files for this assignment:

Submit your written answers (e.g., pencil and paper calculation) to Gradescope (one assignment)
Download this notebook as a .ipynb file and submit to to Gradescope (second assignment)

Data and Computing for Chemical Engineers

Chapter 4: Geometric Aspects of Linear Algebra

Contents