A Tedious Problem?

Thomas J. Kennedy

Contents:

Thus far we have worked through small problems (i.e., is problems with 5 or fewer points). What if we were to try using more points?

x f(x)
0 0
1 1
2 4
3 9
4 16
5 25
6 36
7 49
8 64
9 81
10 100

This is another set of points generated using $f(x) = x^2$. Even with this function… 11 points requires quite a few multiplications. However, we are Computer Scientists!

1 Drafting Pseudocode

Let us draft some matrix multiplication pseudo code!

# Compute XT * X
for i in 0...2
    for j in 0...length(X)
        for k in 0...length(XT)
          XTX[i][j] += XT[i][k] * X[k][j]


# Compute XT * Y
for i in 0...2
    for j in 0...length(X)
        for k in 0...length(XT)
          XTY[i][j] += XT[i][k] * Y[k][j]

Of course, we have violated the Do Not Repeat Yourself (D.R.Y) rule. We wrote the same matrix multiplication code twice. Let us fix that by defining a multiplication function.

def multiply(lhs, rhs):

    # let result be a lhs.rows by rhs.columns matrix

    for i in 0...result.rows
        for j in 0...result.cols
            for k in 0...lhs.cols
              result[i][j] += lhs[i][k] * rhs[k][j]

    return result

That definition allows us to write:

 
XTX <- multiply(XT, X)
XTY <- multiply(XT, Y)

2 Write the Code

Take the matrix multiplication pseudocode and implement it in a language of your choosing. I strongly encourage you to use C++ or Python 3.6+. Once you have implemented the multiplication algorithm test it on our current set of points:

x f(x)
0 0
1 1
2 4
3 9
4 16
5 25
6 36
7 49
8 64
9 81
10 100

You will need to implement the logic to construct $X$, $X^T$, and $Y$. If your implementation is correct, you will end up with:

$$ X = \left[\begin{array}{rr} 1 & 0\\ 1 & 1\\ 1 & 2\\ 1 & 3\\ 1 & 4\\ 1 & 5\\ 1 & 6\\ 1 & 7\\ 1 & 8\\ 1 & 9\\ 1 & 10\\ \end{array}\right] $$

$$ Y = \left[\begin{array}{r} 0\\ 1\\ 4\\ 9\\ 16\\ 25\\ 36\\ 49\\ 64\\ 81\\ 100\\ \end{array}\right] $$

$$ X^T = \left[\begin{array}{rrrrrrrrrrr} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\ 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10\\ \end{array}\right] $$

$$ X^TX = \left[\begin{array}{rr} 11 & 55\\ 55 & 385\\ \end{array}\right] $$

$$ X^TY \left[\begin{array}{r} 385\\ 3025\\ \end{array}\right] $$

3 Solve the Problem By Hand

We need to start with the augmented matrix $[X^TX|X^TY]$

$$ \left[\begin{array}{rr|r} 11 & 55 & 385\\ 55 & 385 & 3025\\ \end{array}\right] $$

We need systematic process to solve the system. For this example, we will use, in order:

  1. Pivot
  2. Scale
  3. Elimination

Let us start by swapping $r_0$ and $r_1$. We want the largest value in (row 0, column 0). This represents the pivot step.

$$ \left[\begin{array}{rr|r} 55 & 385 & 3025\\ 11 & 55 & 385\\ \end{array}\right] $$

Next, we need to scale $r_0$ by $\frac{1}{55}$.

$$ \left[\begin{array}{rr|r} 1 & 7 & 55\\ 11 & 55 & 385\\ \end{array}\right] $$

The third step, elimination, is where $11r_0$ is subtracted from $r_1$.

$$ \left[\begin{array}{rr|r} 1 & 7 & 55\\ 0 & -22 & -220\\ \end{array}\right] $$

That completes the first iteration (i.e., iteration 0). We will repeat the same steps for the second row, $r_1$.

$$ \left[\begin{array}{rr|r} 1 & 7 & 55\\ 0 & -22 & -220\\ \end{array}\right] $$

Since $r_1$ is the last row, the pivot step is a no-op (i.e., nothing changes). We then need to scale $r_1$ by $-\frac{1}{22}$.

$$ \left[\begin{array}{rr|r} 1 & 7 & 55\\ 0 & 1 & 10\\ \end{array}\right] $$

Since $r_1$ is the last row, the elimination step is also a no-op (i.e., nothing changes).

$$ \left[\begin{array}{rr|r} 1 & 7 & 55\\ 0 & 1 & 10\\ \end{array}\right] $$

Our final step is to backsolve using $r_0 = r_0 - 7r_1$.

$$ \left[\begin{array}{rr|r} 1 & 0 & -15\\ 0 & 1 & 10\\ \end{array}\right] $$

3.1 The Result

That leaves us with the coefficients:

$$c_0 = -15$$ $$c_1 = 10$$

and

$$\hat{\varphi} = -15 + 10x^1$$

as our approximation function.

4 I Automated Those Steps!

I did not actually solve this problem by hand. I implemented a Gaussian Elimination based matrix solver. This solver is instrumented and prints the augmented matrix after each step using LaTeX.

You will implement a similar solver in a future exercise.