A Tedious Problem?
Thomas J. Kennedy
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:
- Pivot
- Scale
- 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.