These notes were prepared by a student attending the optional lecture on matrices. As such, they may contain a small error here or there (standard disclaimer). They are provided here as a service to those who may have missed the lecture.

First, here is the code that we wrote: matrices.py

Second, here is the header for that file, which explains the overall goal and approach:

# Partly written in class on Thu 6-Oct with many additions in an
# optional lecture on Sat 8-Oct, and then one more addition (changing
# floats to fractions) by motivated students after the lecture.  :-)
# As usual, with code written in real-time in class, this may have
# some stylistic shortcomings (such as lacking comments) and may
# even have a bug lurking here or there.

# This optional case study solves each of the following problems.  Each
# solution uses the solutions to the preceding problems.

# 1) Multiply two matrices
# 2) Invert a matrix
# 3) Solve a system of linear equations
# 4) Fit an exact polynomial to a list of points
# 5) Find the coefficients of the polynomial for a power sum

# This last problem needs a little explaining.  Consider this power sum:
# 1 + 2 + ... + n.  We know this equals n(n+1)/2.  We can represent this
# quadratic function by its coefficients alone, as [1/2, 1/2, 0] (ignoring
# integer division in this explanation).

# Now, the sum of the squares is: 1**2 + 2**2 + ... + n**2.  This is a cubic
# equation.  It happens to be [1/3, 1/2, 1/6, 0].  That is:
# 1**2 + ... + n**2 = (1/3)n**3 + (1/2)n**2 + (1/6)n

# In fact, for positive integer k, 1**k + ... + n**k is a (k+1)-degree
# polynomial.  The problem we are solving is to write a function which takes
# this number k, and returns the coefficients of that polynomial.

# For some test cases, you can see the expected polynomials for k=1 to k=10
# here: http://mathworld.wolfram.com/PowerSum.html

# Also, our approach for matrix multiplication needs some slight explaining.
# We basically use Gauss-Jordan Elimination to get reduced row echelon form
# while applying the same elementary row operations to an identity matrix -- in
# so doing, the identity matrix is transformed into the inverse of the original
# matrix.  Amazing!
# You can see this in action, interactively entering your own NxN matrix and
# watch the inverse matrix being derived step-by-step, here:
# http://www.math.odu.edu/~bogacki/cgi-bin/lat.cgi?c=inv

# Finally, since our solutions are all fractions, we opted to the fractions
# module in Python.  This hardly changed our original floating-point solution,
# so the code should be readable even if you are unfamiliar with the
# fractions module.

Ultimately we’re going to be able to solve for the coefficients of a power sum, but let’s first go over Matrix Multiplication:

· “K*c^th row”: if K is in column c in the 1^st matrix, it will grab the c^th row of the 2^nd matrix when multiplying

· “into r^th row”: whichever row of the 1^st matrix K appears in will be the row K will be added onto in the answer

m = multiplyRowOfSquareMatrix(m, diagonalRow, k) #make sure to assign this to a variable!

inverse = multiplyRowOfSquareMatrix(inverse, diagonalRow, k) #we could swap rows here instead

assert(m[diagonalRow][col] != 0) #make sure it’s a singular, non-invertible matrix

for targetRow in xrange(n): #targetRow is the variable marching through other rows, NOT sourceRow

Now that we know how to use matrix multiplication to invert a matrix, we can use inverse matrices to solve systems of linear equations!

#Keep in mind that variables must be in the same order across all lines and equations must be linear

A = [ [3, 2, -2], #A consists of the constant coefficients of the left hand side of the system of equations

b = [ [ 10 ], #And b is just the 2d matrix of the right hand side of the above system of equations!

Next, we realize that given a set of points, we can find a polynomial that fits them.

Where the variables x & y are now the constants in the equations and the coefficients are now the variables.

If we have three points, there must be a unique parabola that goes through them.

(To prove this we’d have to show the unique coefficients, and since the coefficients are unique, so is the parabola.)

Also remember to keep in mind that we need three points to describe a polynomial of size 2!

Notice here that the 2^nd matrix holds all the coefficients! This is great because we know from earlier that we can solve for the values in that matrix by multiplying both sides of the equation by the inverse of the 1^st matrix. Thus, we’d be solving for the coefficients.

If we follow this pattern, we know that we can generate as many points as needed. We can then pass these points to our curve fitting function, which solves for the coefficients by solving a system of equations.

# We argued by handwaving-ish-calculusy-stuff that f(n) is a polynomial of degree (k+1)

for x in xrange(1,k+3): #we removed 0 as the first element to avoid a row swap, so this becomes k+3 instead of k+2

So we originally set out to find the power sum coefficients and now have the ability to take a set of points and fit them to a curve. We can also take a polynomial with n equations and n unknowns and use matrix multiplication and matrix inverses to solve the system of linear equations. Pretty cool!