# Matrix Operations

Bart Massey

### Projects For Today

• Matrix multiplication

• Classical
• Strassen
• (FFT)
• Solving a system of linear equations

• LUP decomposition
• Solving the resulting system
• Matrix inversion

• Least-squares minimization

### Matrix Multiplication

• Given matrices A of size (m, n) and B of size (n, p)

• Product is C[i, j] = sum(k <- 1..n, A[i, k] * B[k, j])

• Pseudocode:

multiply A by B:
for i <- 1 .. m
for j <- 1 .. p
C[i, j] <- A[i, 1] * B[1, j]
for k <- 2 .. n
C[i, j] <- C[i, j] + A[i, k] * B[k, j]
return C

• mp space for result, mnp time

• For square matrices, these are O(n^2) and O(n^3)

### Block Multiplication

• For a square matrix, let

A = [ A[1,1] A[1,2] ]
[ A[2,1] A[2,2] ]


and similarly for B and C

• Compute the matrix product as

C[i,j] = sum(k <- 1..2, A[i, k] . B[k, j])


where "." is matrix multiplication

• Gives the same answer (easy proof), but recursive

### Strassen's Method

• "Strassen's method is not at all obvious. (This might the biggest understatement in this book)."

• Let S[1] = B[1,2] - B[2,2] and then define S[2..10] as sums or differences of blocks according to p. 80

• Let P[1] = A[1,1] . S[1] and then define P[2..7] similarly according to p. 80

• C[(1,1)..(2,2)] can be written as sums of various P according to p. 81

• Note that only the P calculations have multiply, and that there are only 7 of them.

• Traditional block multiplication would give 2x2x2 = 8 multiplies
• The P multiplies can be done by recursively applying Strassen's algorithm! This gives complexity O(n^(lg 7)) ~= O(n^2.8)

### Solving Systems Of Linear Equations

• Solving for n independent linear equations in n unknowns with unique answer

3*x1 + 5*x2 + 7*x3 = 9
4*x1 + 6*x2 + 8*x3 = 10
1*x1 + 0*x2 + 3*x3 = 4

• Watch out for underconstrained systems due to dependent equations, overconstrained systems due to unsolvability

• Matrix formulation: write

b[1] = a[1, n] * x[n] + ... + a[1,1] * x[1]
b[2] = a[2, n] * x[n] + ... + a[2,1] * x[1]
...


as

Ax = b


where b and x are vectors of length n, and A is an (n, n) matrix of coefficients

### Matrix Multiply By FFT

• Algorithm

• Take DFT of input matrices
• Then add the DFTs
• Then take inverse DFT of the result
• Most efficient known algorithm for very large n

### Gaussian Elimination: Row-Reduction To "LUP" Form

• Equations (rows) can be multiplied by a constant, added, swapped around.

• Idea: Use first coefficient in first row to get rid of first coefficient in all other rows.

a[2] -= a[1] * (a[2,1] / a[1,1])
b[2] -= b[1] * (a[2,1] / a[1,1])


Continue to eliminate successive coefficients from successive rows.

• Permutation is needed if divisor is zero or small. Preserve the permutation to keep track of which equation was which!

### Gaussian Elimination: Solving From "LUP" Form

• Note that a[n,n] * x[n] = b[p[n]] is trivially solvable for x[n]

• Substitute x[n] and subtract resulting constant from b[p[i]] to get equations like

a[n-1,n-1] * x[n-1] = ..

• Continue until all x[i] have been computed, and return the results

• Running time of decomposition, solving is O(n^2)

### Matrix Inversion

• Finding A^-1 such that A . A^-1 = I

• Can use Gaussian Elimination to find

• Set up n systems of equations, each of which solves for a column of A^-1
• Solve the systems in total time O(n^3)
• Proof of equivalence of multiplication and inversion

• Important property: two-way proof

### Least-Squares Approximation

• Practically super-important: fit a function to data points minimizing square of error in y

• See last section of the book and meditate hard for details

• Look somewhere else (e.g. Numerical Recipes, web) for implementation hints

Last modified: Monday, 12 May 2014, 4:07 PM