Practical 4: Practising practical Linear Algebra

This week is focussed on the key ideas from linear algebra which are the building blocks of regression, computational modelling and machine learning.

Note

Don’t be intimidated by the maths!

Even if you’re struggling with the maths try to focus on what you think the code is doing.

In the tutorial we’re going to explore the mathematical concepts from the lecture through an example.

A problem with bike parts

Background

A factory specialises in producing parts for bicycles: handlebars, brakes and wheels.

• \(x_1\) — the number of handlebars they can make
  
• \(x_2\) — the number of brakes they can make
• \(x_3\) — the number of wheels they can make

These require resources to be made: aluminium, rubber and paint.

For each handlebar they need 8 units of aluminium, 4 units of rubber and 4 units of paint. For each brake they need 1 unit of aluminium, 1 unit of rubber and 1 unit of paint. For each wheel they need 10 units of aluminium, 23 units of rubber, and 5 units of paint.

The factory is having problems with the supply of raw materials. On a given day they have:

• \(b_1=400\) units of aluminium,
  
• \(b_2=750\) units of rubber,
  
• \(b_3=200\) units of paint.

Assuming fractional quantities are allowed, what quantities of each bike part would use all three resources exactly?

This is a simplified model. Real production requires non-negative whole numbers of parts and may leave resources unused. Finding an optimal plan would also require an objective, such as maximising profit.

Linear equations

We can express the resource constraints as a system of equations. If they are producing \(x_1, x_2, x_3\) units of each part we need:

\[ 8x_1 + x_2 + 10x_3 \]

of aluminium.

\[ 4x_1 + x_2 + 23x_3 \]

of rubber.

\[ 4x_1 + x_2 + 5x_3 \]

of paint.

If all available resources are used exactly, we have:

\[ 8x_1 + x_2 + 10x_3 = b_1 = 400 \]

\[ 4x_1 + x_2 + 23x_3 = b_2 = 750 \]

\[ 4x_1 + x_2 + 5x_3 = b_3 = 200 \]

Matrix form

This can be written compactly as:

\[ \begin{pmatrix} 8 & 1 & 10 \\ 4 & 1 & 23 \\ 4 & 1 & 5 \end{pmatrix} \begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix} = \begin{pmatrix}b_1\\b_2\\b_3\end{pmatrix}. \]

That is, \[ Mx =b \]

Where \(M\) is the coefficient matrix,

\[ M = \begin{pmatrix} 8 & 1 & 10 \\ 4 & 1 & 23 \\ 4 & 1 & 5 \end{pmatrix} \]

\(x\) is the vector of the number of each part, i.e. the variables in our equation:

\[ x = \begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix} \]

and \(b\) is the resource constraint:

\[ b = \begin{pmatrix}b_1\\b_2\\b_3\end{pmatrix}. \]

Writing this in Python

Let’s start by importing the libraries we need in Python.

import numpy as np

We can write this same information in Python by creating variables.

# Coefficient matrix
# M 
M = np.array([[8, 1, 10],
              [4, 1, 23],
              [4, 1, 5]], dtype=float)

# Resource availability
# [b_1, b_2, b_3]
b = np.array([400, 750, 200], dtype=float)  

Here we create a two-dimensional NumPy array from a list of rows. M.shape is (3, 3), while b.shape is (3,): NumPy represents b as a one-dimensional array.

Question

What does the entry 23 in M represent? What do the rows and columns represent? Discuss with a neighbour before checking the answer.

Answer

Each wheel requires 23 units of rubber. The rows represent aluminium, rubber and paint, in that order. The columns represent handlebars, brakes and wheels. Each entry gives the amount of one resource required for one unit of a product.

Solving the equations

We can use in-built functions in Python to solve the system of linear equations. The library np.linalg provides us with all the matrix operations we need.

Checking the determinant

A square coefficient matrix with a non-zero determinant is invertible, so the system has a unique solution for every resource vector \(b\). If the matrix is singular (not invertible), the system can have no solution or infinitely many solutions. Let’s check the determinant of \(M\).

# Calculate the determinant of M
det_M = np.linalg.det(M)

print(det_M)

-71.99999999999994

Question

The determinant is approximately -72. Does its negative sign prevent a unique solution? What would a determinant of zero tell us?

Answer

The sign does not prevent a unique solution: the determinant is non-zero, so \(M\) is invertible. A zero determinant would mean that the matrix is singular. We would then need to examine the equations and resource vector to determine whether there are no solutions or infinitely many.

Inverting the matrix

Now let’s calculate the inverse of M.

Note

Remember: The inverse of a matrix \(M\) is denoted \(M^{-1}\) and satisfies \(MM^{-1}=I\), where \(I\) is the identity matrix.

# Invert M
M_inv = np.linalg.inv(M)

print(M_inv)

[[ 0.25       -0.06944444 -0.18055556]
 [-1.          0.          2.        ]
 [-0.          0.05555556 -0.05555556]]

The identity matrix

Let’s check that \(MM^{-1}= I\), where I is the 3x3 identity matrix:

\[ \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \]

# Multiply M by its inverse

# @ performs matrix multiplication; * multiplies element by element
MM_inv = M @ M_inv
print(MM_inv)

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

Hmmmm… that’s weird. It’s not quite the identity matrix.

Floating-point calculations can introduce small rounding errors. Here, entries around \(10^{-17}\) are effectively zero for our purposes.

Rounding the matrix we see that it is very close to the identity.

print(MM_inv.round(2))

[[ 1. -0.  0.]
 [ 0.  1. -0.]
 [ 0. -0.  1.]]

The computed product is approximately the identity matrix. Rounding changes the display; we can check numerical agreement within a tolerance using np.allclose:

print(np.allclose(MM_inv, np.eye(3)))

True

Let’s solve it!

As we saw in the lectures, we can get \(x\) as follows:

\[ x = M^{-1}b \]

We use @ again, this time to multiply the inverse matrix by the resource vector.

# Multiply the inverse matrix by the resource vector
x = M_inv @ b

print(np.round(x,2))

[11.81  0.   30.56]

We have approximately \(x_1=11.8056\), \(x_2=0\) and \(x_3=30.5556\). These fractional quantities solve the exact-use equations, but do not describe a practical whole-number production plan.

Question

Can the factory manufacture the exact solution? If it produces 11 handlebars, 0 brakes and 30 wheels instead, how much of each resource remains? Could it make any additional parts?

Answer

The exact solution includes fractions of parts, so it cannot be manufactured as whole parts. The rounded-down plan uses 388 units of aluminium, 734 of rubber and 194 of paint, leaving 12, 16 and 6 respectively. It could therefore make six additional brakes. Rounding down does not establish an optimal production plan.

x_whole = np.array([11, 0, 30])
print("Resources used:", M @ x_whole)
print("Resources remaining:", b - M @ x_whole)

Resources used: [388. 734. 194.]
Resources remaining: [12. 16.  6.]

Let’s solve it again!

The inverse helps explain the algebra. For routine calculations, use NumPy’s np.linalg.solve, which solves the system directly without explicitly constructing the inverse.

## Solve the equations using numpy's linear algebra solver
x = np.linalg.solve(M, b)

print(np.round(x, 2))

[11.81  0.   30.56]

Check the solution by substituting it back into the original equations. The residual b - M @ x should be zero apart from small numerical rounding errors.

print("Resources used:", M @ x)
print("Residual:", b - M @ x)
print(np.allclose(M @ x, b))

Resources used: [400. 750. 200.]
Residual: [0. 0. 0.]
True

Question

Suppose all three resource quantities double. Predict the new solution before calculating it. Explain your reasoning and check it in Python.

Answer

The solution doubles: approximately 23.6111 handlebars, 0 brakes and 61.1111 wheels. Since \(Mx=b\), we have \(M(2x)=2b\). This uses the same resource requirements and still allows fractional quantities.

b_double = 2 * b
x_double = np.linalg.solve(M, b_double)
print(x_double)
print(np.allclose(x_double, 2 * x))

[23.61111111  0.         61.11111111]
True

Visualising

In two and three dimensions we can visualise a system of linear equations. When there is a unique solution, two lines in two dimensions or three planes in three dimensions meet at a single point. Here each plane represents exact use of one resource, and their common intersection is our solution.

We can plot our equations in 3D space to see this.

import plotly.graph_objects as go

# Define grid for x1 (handlebars) and x2 (brakes)
x1 = np.linspace(0, 30, 50)
x2 = np.linspace(0, 30, 50)
X1, X2 = np.meshgrid(x1, x2)

# Constraint surfaces (solving for x3 = wheels)
Z_al = (b[0] - 8*X1 - 1*X2) / 10 # from aluminium
Z_rb = (b[1] - 4*X1 - 1*X2) / 23 # from rubber
Z_pn = (b[2] - 4*X1 - 1*X2) / 5 # from paint

# Interactive 3D plot
fig = go.Figure()

fig.add_surface(x=x1, y=x2, z=Z_al, opacity=0.5, colorscale="Blues", name="Aluminium")
fig.add_surface(x=x1, y=x2, z=Z_rb, opacity=0.5, colorscale="Oranges", name="Rubber")
fig.add_surface(x=x1, y=x2, z=Z_pn, opacity=0.5, colorscale="Greens", name="Paint")

# Add solution point
fig.add_trace(go.Scatter3d(x=[x[0]], y=[x[1]], z=[x[2]],
mode="markers", marker=dict(size=6, color="red"),
name="Solution"))

fig.update_layout(
scene=dict(
xaxis_title="Handlebars",
yaxis_title="Brakes",
zaxis_title="Wheels"
))

fig.show()

Note: keep M, b and x unchanged when exploring alternative resource quantities, so that the plotted planes and solution point describe the same system.

Conclusions

Here we solved a set of linear equations using algebra. Solving sets of linear equations is closely related to regression, a topic we will learn about in week 6.

Extension

If you have time you can explore implementing an ordinary least squares regression.

We have data from the factory - every day for a year they recorded how many handlebars, brakes and wheels they made and the corresponding amount of each resource used. We can load the data as:

import pandas as pd 

df = pd.read_csv('https://raw.githubusercontent.com/huanfachen/QM/refs/heads/main/sessions/L4_data/factory_data.csv')

df.head()

	x1_handlebars	x2_brakes	x3_wheels	aluminium_used	rubber_used	paint_used
0	38	12	28	601.004775	797.467413	300.493185
1	28	29	31	560.707839	862.138639	286.595862
2	14	18	18	312.346520	490.527555	180.463972
3	7	16	20	261.929106	497.427667	134.626443
4	20	18	4	201.302208	207.646869	105.584988

Choose one dependent variable (“aluminium_used”, “rubber_used” or “paint_used”), and using the values for the number of items produced try to estimate the coefficients (an intercept and three slopes).

Earlier we knew the resource requirements and solved for production quantities. Now we observe production quantities and resource use, and estimate the resource requirements: one row of \(M\), plus an intercept.

For aluminium, our regression model is

\[ \text{aluminium used} = \beta_0 + \beta_1x_1 + \beta_2x_2 + \beta_3x_3 + \varepsilon. \]

The error term \(\varepsilon\) allows observed resource use to differ from the model prediction.

Hint: You might want to use statsmodels OLS package which you can read about here.

Question

Complete the code for one resource. Interpret the three slopes and the intercept. For aluminium, compare your estimated slopes with 8, 1 and 10. Why might they differ?

import statsmodels.api as sm

# define independent variables (production quantities)
X = ??

# add constant term for intercept
X = sm.add_constant(X)

# define dependent variable (e.g. aluminium used)
y = ??

# fit OLS regression model
model = sm.OLS(y, X).fit()

# print model summary
print(model.summary())

Answer

Each slope estimates the additional resource use per extra unit of that product, holding the other production quantities fixed. The intercept estimates resource use when all three production quantities are zero; the original model assumes this is zero. Observed variation and estimation uncertainty mean the fitted coefficients need not match the original requirements exactly.

import statsmodels.api as sm

# define independent variables (production quantities)
X = df[["x1_handlebars",    "x2_brakes",    "x3_wheels"]]

# add constant term for intercept
X = sm.add_constant(X)

# define dependent variable (e.g. aluminium used)
y = df["aluminium_used"]

# fit OLS regression model
model = sm.OLS(y, X).fit()

# print model summary
print(model.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:         aluminium_used   R-squared:                       0.996
Model:                            OLS   Adj. R-squared:                  0.996
Method:                 Least Squares   F-statistic:                 3.103e+04
Date:                Mon, 26 Oct 2026   Prob (F-statistic):               0.00
Time:                        17:15:30   Log-Likelihood:                -1351.0
No. Observations:                 365   AIC:                             2710.
Df Residuals:                     361   BIC:                             2726.
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
=================================================================================
                    coef    std err          t      P>|t|      [0.025      0.975]
---------------------------------------------------------------------------------
const             2.0457      1.578      1.297      0.196      -1.057       5.148
x1_handlebars     8.0292      0.044    183.482      0.000       7.943       8.115
x2_brakes         0.9331      0.045     20.840      0.000       0.845       1.021
x3_wheels         9.9771      0.044    229.045      0.000       9.891      10.063
==============================================================================
Omnibus:                        2.119   Durbin-Watson:                   2.129
Prob(Omnibus):                  0.347   Jarque-Bera (JB):                1.793
Skew:                           0.027   Prob(JB):                        0.408
Kurtosis:                       2.661   Cond. No.                         112.
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

You’re Done!

Congratulations on completing the linear algebra practical! If you are still working on it, take your time.

Don’t worry here about understanding every detail of the maths - focus on each of the operations we looked at. Do you know what is being calculated here? Why is that a useful thing to calculate? If you’re struggling to understand this just ask!

It takes time to learn - remember practice makes perfect!