Devlog 0x11 - Gaussian Elimination in Zig

Today, I was forced to use a Jupyter Notebook to do some row reduction problems. Juptyer Notebooks are not my style, so I wrote a terminal based Gaussian Elimination program in Zig.

Before we begin, if you'd like to use zrref(zig reduced row echelon form), you can checkout the github page: zrref github

Why?

My Linear Algebra course has it's third worksheet of the semester due tomorrow, and in one of the problems it provides a .jpynb file. After looking up the file extension and seeing that it's a Jupyter Notebook file, I decided to open it up in vim to see what I was working with.

I open the file, and to my horror this is what I find:

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "# Problem 2 part (a): Row reducing a matrix\n",
        "\n",
        "Run the code below and use the output to answer the homework questions."
      ],
      "metadata": {
        "id": "5CeJK40jjeuO"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Row Reduction in Python\n",
        "\n",
        "# Importing necessary libraries -- you need to do this once for each notebook\n",
        "import numpy as np\n",
        "from sympy import Matrix\n",
        "\n",
        "# Step 1: Use numpy to generate random numbers.\n",
        "np.random.seed(0)  # For reproducibility -- the \"random numbers\" should be the same for everyone\n",
        "A = np.random.randint(-5, 6, size=(3, 4))  # a 3x4 matrix full of random integers between -5 and 5. (Including -5, but not including 6; the upper bound is exluded. Weird, I know.)\n",
        "\n",
        "print(\"Original Matrix A:\")\n",
        "print(A)\n",
        "\n",
        "# Step 2: Use \"Matrix\" and sympy commands to reduce the matrix\n",
        "A_sympy = Matrix(A)   # Create a matrix our of the array A\n",
        "A_echelon = A_sympy.echelon_form()       # perform row reduction to echelon form\n",
        "A_ref, pivot_columns = A_sympy.rref()   # perform row reduction to reduced echelon form\n",
        "\n",
        "# Print results\n",
        "print(\"\\nEchelon Form (EF) of A:\")\n",
        "print(np.array(A_echelon))\n",
        "\n",
        "print(\"\\nReduced Echelon Form (REF) of A:\")\n",
        "print(np.array(A_ref))\n",
        "\n",
        "# Isn't this nice?  Just be careful with the labeling.  Column \"zero\" is really the first column; everything is shifted back 1.\n",
        "print(\"\\nPivot Columns:\")\n",
        "print(pivot_columns)\n"
      ],
      "metadata": {
        "id": "LE7fGwtK0gpH"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Printing the results for better readability (This code supports problem 2(a) in HW3)\n",
        "\n",
        "Okay, it works!  But, the formatting of the printing might make the numbers misaligned (e.g. the EF of the matrix is not great).  Here's how we can change the printing so that the numbers in the matrix are aligned.  You may need to adjust the \"8\" to a larger number if the numbers in the matrix get big."
      ],
      ...

This gives me a few options, but none of them are particularly good.

  1. Use Jupyter Notebook
  2. Install a Jupyter Notebook Neovim Plugin
  3. Make my own Gaussian Elimination program

What is Row Reduction?

Row reduction (also known as Gaussian elimination) is a systematic method for solving systems of linear equations by transforming them into what's called "row echelon form." The process gets its name from Carl Friedrich Gauss and involves three basic operations:

  1. Swapping rows
  2. Multiplying a row by a non-zero scalar
  3. Adding a multiple of one row to another row

For example, consider this system of equations:

2x + y - z = 8
-3x - y + 2z = -11
-2x + y + 2z = -3

Which we can write as an augmented matrix:

[ 2  1  -1 | 8  ]
[-3 -1   2 |-11 ]
[-2  1   2 |-3  ]

And the row reduction, or Gaussian elimination, process transforms this into:

[1 0 0 | 2]  → x = 2
[0 1 0 | 3]  → y = 3
[0 0 1 |-1]  → z = -1

Gaussian Elimination preserves the solution set of the original matrix while making solutions easy to see. This is one of the key fundamentals of Linear Algebra and because of this fact(and I didn't want to use a Jupyter Notebook) I created a row reducing program in Zig.

I created this implementation following the SymPy implementation available on the SymPy Github.

Row Reduction in Zig

One choice I made while creating this program was to use an epsilon value for determining if a solution was a "good one". Essentially, floating point operations are inaccurate and don't work well for comparisons, so we just see if the absolute value of the difference between the divisor(leading term in the row) and 1 is less than some epsilon value. If this case holds true then it means that the row has successfully been reduced, and the constant factor needed for reduction has been found.

Here are my functions for swapping and reducing rows:

pub fn swapRows(self: *Self, row1: usize, row2: usize) void {
          var i: usize = 0;
          while (i < self.cols) : (i += 1) {
              const temp = self.at(row1, i).*;
              self.at(row1, i).* = self.at(row2, i).*;
              self.at(row2, i).* = temp;
          }
      }

pub fn rowReduce(self: *Self) void {
    print("\n{s}Starting row reduction...{s}\n", .{ YELLOW, RESET });
    
    var lead: u32 = 0;
    var row: usize = 0;

    while (row < self.rows) : (row += 1) {
        if (lead >= self.cols) break;

        // Find row with non-zero entry in lead column
        var i = row;
        while (@abs(self.at(i, lead).*) < 1e-10) {
            i += 1;
            if (i == self.rows) {
                i = row;
                lead += 1;
                if (lead == self.cols) return;
            }
        }

        // Swap rows if necessary
        if (i != row) {
            print("{s}Swapping row {d} with row {d}{s}\n", .{ CYAN, row + 1, i + 1, RESET });
            self.swapRows(i, row);
            self.printMatrix();
        }

        // Scale the row
        const divisor = self.at(row, lead).*;
        if (@abs(divisor - 1) > 1e-10) {
            print("{s}Scaling row {d} by {d:.3}{s}\n", .{ MAGENTA, row + 1, 1.0 / divisor, RESET });
            var col: usize = 0;
            while (col < self.cols) : (col += 1) {
                self.at(row, col).* /= divisor;
            }
            self.printMatrix();
        }

        // Eliminate in other rows
        var r: usize = 0;
        while (r < self.rows) : (r += 1) {
            if (r != row) {
                const factor = self.at(r, lead).*;
                if (@abs(factor) > 1e-10) {
                    print("{s}Subtracting {d:.3} times row {d} from row {d}{s}\n", 
                        .{ GREEN, factor, row + 1, r + 1, RESET });
                    var col: usize = 0;
                    while (col < self.cols) : (col += 1) {
                        self.at(r, col).* -= factor * self.at(row, col).*;
                    }
                    self.printMatrix();
                }
            }
        }

        lead += 1;
    }
    print("{s}Row reduction complete!{s}\n", .{ BLUE, RESET });
}

And that's about it! With these two functions we're able to reduce any matrix into row echelon form, and potentially even reduced row echelon form. Thanks for reading and I'll see you next time :)