Python-Workshop-3-NumPy.md

@@ -230,6 +230,9 @@ Problem](https://www.cs.odu.edu/~tkennedy/cs417/s21/Public/approximationExample2
 The code snippets in this section are extracted from
 [least_squares.py](https://git-community.cs.odu.edu/tkennedy/python-workshop/-/blob/master/NumPy/least_squares.py).
 
+
+## Getting Started
+
 The first line is the NumPy import. By convention the `numpy` module is given
 the alias `np`.
 
@@ -295,21 +298,20 @@ After setting everything up... it is time for a couple matrix multiplications.
 The augmented `XTX|XTY` matrix is what we are going to solve. *Note that the
 `@` operator is often used instead of `np.matmul`.*
 
----
-
-```python
-def _backsolve(matrix_XTX, matrix_XTY):
 
-    num_rows, _ = matrix_XTX.shape
+## The Solver
 
-    for i in reversed(range(1, num_rows)):
-        for j in reversed(range(0, i)):
-            s = matrix_XTX[j, i]
+A Gaussian Elimination Matrix solver can viewed as [4 operations](https://www.cs.odu.edu/~tkennedy/cs417/s21/Assts/matrixSolverExercise/):
 
-            matrix_XTX[j, i] -= (s * matrix_XTX[i, i])
-            matrix_XTY[j] -= (s * matrix_XTY[i])
+  1. Pivot
+  2. Scale
+  3. Eliminate
+  4. Backsolve
 
+For our implementation... these operatrionrs will be wrapped in a
+`solve_matrix` function...
 
+```python
 def solve_matrix(matrix_XTX, matrix_XTY):
     """
     Solve a matrix and return the resulting solution vector
@@ -319,6 +321,23 @@ def solve_matrix(matrix_XTX, matrix_XTY):
     num_rows, num_columns = matrix_XTX.shape
 
     for i in range(0, num_rows):
+```
+
+Note the use of the `.shape` tuple.
+
+
+### Pivot
+
+**Pseudocode**
+
+```
+    idx <- find_largest_row_by_col(A, col_index, num_rows)
+    swap_row(A, i, idx)
+```
+
+**Implementation**
+
+```python
         # Find column with largest entry
         largest_idx = i
         current_col = i
@@ -327,16 +346,54 @@ def solve_matrix(matrix_XTX, matrix_XTY):
             if matrix_XTX[largest_idx, i] < matrix_XTX[j, current_col]:
                 largest_idx = j
 
-        # Swap
         if largest_idx != current_col:
             matrix_XTX[[i, largest_idx], :] = matrix_XTX[[largest_idx, i], :]
             matrix_XTY[[i, largest_idx]] = matrix_XTY[[largest_idx, i]]
+```
+
+
+### Swap
+
+**Pseudocode**
+
+```
+def scale_row(A, row_idx, num_cols, s):
+
+    for every j in 0 to num_cols:
+        A[row_idx][j] = A[row_idx][j] / s
+```
+
+**Implementation***
 
         # Scale
         scaling_factor = matrix_XTX[i, i]
         matrix_XTX[i, :] /= scaling_factor
         matrix_XTY[i] /= scaling_factor
+```
+
+### Eliminate
+
+**Pseudocode**
+
+```
+def eliminate(A, src_row_idx, num_cols, num_rows):
+
+    start_col <- src_row_idx
 
+    for every i in (src_row_idx + 1) to num_rows:
+
+        s <- A[i][start_col]
+
+        for every j in start_col to num_cols:
+            A[i][j] = A[i][j] - s * A[src_row_idx][j]
+
+        A[i][start_col] = 0
+
+```
+
+**Implementation**
+
+```python
         # Eliminate
         for row_i in range(i + 1, num_rows):
             s = matrix_XTX[row_i][i]
@@ -350,41 +407,42 @@ def solve_matrix(matrix_XTX, matrix_XTY):
     _backsolve(matrix_XTX, matrix_XTY)
 
     return matrix_XTY
+```
 
+### Backsolve
 
-def main():
+**Pseudocode**
 
-    # Set up input data points, X, Y, and XT
-    points = [(0., 0.), (1., 1.), (2., 4.)]
+```
+def back_solve(matrix):
 
-    matrix_X = np.array([[1., 0., 0.],
-                         [1., 1., 1.],
-                         [1., 2., 4.]])
+    augColIdx <- matrix.cols()  # the augmented column
+    lastRow = matrix.rows() - 1
 
-    matrix_Y = np.array([0,
-                         1,
-                         4])
+    for i in lastRow down to 1:
+        for j <- (i - 1) down to 0:
+            s <- matrix[j][i]
 
-    matrix_XT = matrix_X.transpose()
+            matrix[j][i] -= (s * matrix[i][i])
+            matrix[j][augCol] -= (s * matrix[i][augColIdx])
+```
 
-    # Compute XTX and XTY
-    matrix_XTX = np.matmul(matrix_XT, matrix_X)
-    matrix_XTY = np.matmul(matrix_XT, matrix_Y)
+**Implementation**
 
-    print_matrices(matrix_XTX, matrix_XTY)
+```python
+def _backsolve(matrix_XTX, matrix_XTY):
 
-    print()
-    print("{:-^40}".format("Solution"))
-    solution = solve_matrix(matrix_XTX, matrix_XTY)
-    print(solution)
+    num_rows, _ = matrix_XTX.shape
 
+    for i in reversed(range(1, num_rows)):
+        for j in reversed(range(0, i)):
+            s = matrix_XTX[j, i]
 
-if __name__ == "__main__":
-    main()
+            matrix_XTX[j, i] -= (s * matrix_XTX[i, i])
+            matrix_XTY[j] -= (s * matrix_XTY[i])
 ```
 
-
-## Statistics
+# Statistics