a logo for the course

Matrix Computations

David Gleich

Purdue University

Fall 2017

Course number CS-51500

Tuesday and Thursday, 10:30-11:45am

Location Forney B124


Homework 4

Homework 4

Please answer the following questions in complete sentences in a typed, clearly prepared manuscript and submit the solution by the due date on Blackboard (early morning on Monday, September 25th, 2017)

Updates

Problem 0: Homework checklist

Problem 1: Backsolve

For this problem, you must use Julia.

  1. Implement backsolve and forward solve as functions in Julia. Show and document your code.

  2. Construct a random upper-triangular linear system via:

    A = triu(randn(n,n));
    b = randn(n);
    

    Compare the performance of your backsolve to Julia's backslash method to solve a linear system.

  3. Write a routine to solve a linear system of equations using the LU factorization in Julia both with and without pivoting.
    (Use lufact(A, Val{false}) for no pivoting in LU in Julia.)

  4. Use your backsolve code, along with Julia's lu factorization in order to implement your own linear solver. Present a paragraph or two (and a figure or two) comparing it's speed and accuracy to Julia's own solver and focus on the effect of pivoting.

Problem 2: Poisson's equation

In this problem, we'll meet one of the most common matrices studied in numerical linear algebra: the -Laplacian. We arrive at this matrix by discretizing a partial differential equation. Poisson's equation is: where is a continuous function defined over the unit-plane (i.e. ), is a continuous function defined over the same region, and is the Laplacian operator:

Given a function , we want to find a function that satifies this equation. There are many approaches to solve this problem theoeretically and numerically. We'll take a numerical approach here.

Suppose we discretize the function at regular points , and where so that we have: For this discretization, note that What we've done here is use the approximation: for both partial terms.
We need this equation to hold at each point . But note that there are some issues with this equation at the boundary values (where x=0 or 1, or where y=0 or 1).
For this problem, we'll make it very simple and set:

Now, we'll do what we always do! Turn this into some type of matrix equation!

Let be an matrix that we'll index from zero instead of one: where . At this point, we are nearly done. What we are going to do is turn Poisson's equation into a linear system. This will be somewhat like how we turned image resampling into a matrix vector equation in the first homework.

There was a typo here before with row-vs-column indexing. Either is okay for this problem. Just be consistent and mention what you are using.

In order to write as a vector, we'll keep the convention from last time: Let be the vector of elements here. Note that our approximation to , just involved a linear combination of the elements of . This means we have a linear system: where the rows of and correspond to equations of the form:

  1. Let . Write down the linear equation for including all the boundary conditions. Note that you can encode the boundary conditions by adding a row of where: .

  2. Write a Julia or Matlab/Python code to construct the matrix and vector when and . Here's some pseudo-code to help out:

    function laplacian(n::Integer, f::Function)
        N = (n+1)^2
        A = zeros(N,N)
        fvec = zeros(N)
        # 2017-09-21 added transpose to mirror the row-indexing in the problem
        # but see note above, you can do either, just be consistent.
        G = reshape(1:N, n+1, n+1)' # index map, like we saw before; 
        h = 1.0/(n)
        for i=0:n
            for j=0:n
                row = G[i+1,j+1]
                if i==0 || j == 0 || i == n || j == n
                    # we are on a boudnary
                    fvec[row] = 0.0
                    # fill in A[row,:]
                else
                    fvec[row] = f(i*h, j*h)*h^2
                    # fill in A[row,:]
                end
            end
        end
        return A, fvec
    end
    A, f = laplacian(10, (x,y) -> 1.0)
    
  3. Solve for using Julia's or Matlab's backslash solver, and show the result using the mesh function (Matlab) or surface function (Plots.jl in Julia).

Problem 3: The Schur Complement

Suppose we wish to solve and further suppose that you KNOW some of the values of .

Let us permute and partition to be a block system: where is what you know.

  1. Show how to solve for given . Under what conditions is this possible?

  2. Now, suppose that you have a very kind, but very confused dog that happened to eat your flash memory stick holding the piece of that you knew. However, you had saved your computed on your Purdue account, and so you have a backup. (This means you can assume that computing from is possible for this problem if you determined it wasn't always possible above.) Can you get back?

  3. Combine these two parts to derive a single linear system to compute without computing . The system you'll derive is called the Schur complement.

Problem 4: Ranking with Linear Systems

Watch the video on sports linear systems available on Blackboard. (this was recorded in 2014). No questions about it, but this material is now fair-game for the Midterm. Your answer to this homework problem should be a question about the material presented in the lecture.