A Gonum Tutorial (Part 1)

December 27 2023

Objective

This post is meant to provide a small introduction to gonum, one of the leading numeric libraries for the Go programming language. Rather than reviewing how a specific function works (which can be done by reading some of the tests like this in the source code), I will try to walk through how to use gonum to achieve a single task.

Example 1: Traffic Pricing

Let's consider the problem of determining prices for a set of toll roads into a city like Jersey City, New Jersey, USA. For the sake of simplicity, we will say that there are only three main roads into the city and that some helpful traffic engineers have predictions for how many cars, light trucks, and heavy trucks will use each road in a given month.

\( \begin{array}{rccc} & \text{Cars} & \text{Light Trucks} & \text{Heavy Trucks} \\ \text{Road 1} & 5.0 \times 10^4 & 1.0 \times 10^4 & 6.0 \times 10^4 \\ \text{Road 2} & 2.0 \times 10^4 & 2.5 \times 10^4 & 1.0 \times 10^3 \\ \text{Road 3} & 5.0 \times 10^4 & 5.0 \times 10^4 & 1.5 \times 10^3 \\ \text{Road 4} & 1.5 \times 10^5 & 1.5 \times 10^4 & 0.0 \\ \end{array} \)

The city wants to charge a single toll for each vehicle type that enters at each of these points (i.e., each road doesn't have the freedom to choose its own tolls). It is our goal to find those prices!

The city also has a set of target incomes that it wants to make from each road. You can think about this as accounting for the costs to maintain each unique road. For example, perhaps Road 1 runs over a bridge that needs to be maintained and thus it requires higher revenues than other roads.

The target incomes for each road are as follows:

\( \begin{array}{rc} & \text{Target Revenue} & \\ \text{Road 1} & $1,100,000 \\ \text{Road 2} & $200,000 \\ \text{Road 3} & $500,000 \\ \text{Road 4} & $350,000 \\ \end{array} \)

This problem is a classic instance of a problem that we call: The Least-Squares problem.

Translating the Example Into Math: Least-Squares Problem

The least-squares problem is defined mathematically in the following way:

\( \min_{x} \|Ax - b\|_2^2 \)

Intutitively, what this problem expresses is that we want to find the price values (i.e., the value of the vector \( x \in \mathbb{R}^n \)) that minimizes the difference between the predicted revenue and the target revenue.

There is a lot of theory that goes into how people solve this problem, but we will present the simplest way to solve it using gonum.

A Simple Implementation

Consider the following program:

        
package main

import (
	"fmt"
	"gonum.org/v1/gonum/mat"
)

func main() {
	// Define the matrix A
	A := mat.NewDense(4, 3, []float64{
		5.0e4, 1.0e4, 6.0e4,
		2.0e4, 2.5e4, 1.0e3,
		5.0e4, 5.0e4, 1.5e3,
		1.5e5, 1.5e4, 0.0,
	})

	// Define the vector b
	b := mat.NewVecDense(4, []float64{
		1.1e6,
		2.0e5,
		5.0e5,
		3.5e5,
	})

	// Solve the least-squares problem
	var x mat.VecDense
	x.SolveVec(A, b)

	// Print the solution
	fmt.Printf("x = %.4f\n", mat.Formatted(&x, mat.Prefix("    ")))

	// Print the product
	var prod mat.VecDense
	prod.MulVec(A, &x)

	fmt.Printf("A*x = %.4f\n", mat.Formatted(&prod, mat.Prefix("    ")))

        // Print the residual
	var res mat.VecDense
	res.SubVec(b, &prod)

	fmt.Printf("b - A*x = %.4f\n", mat.Formatted(&res, mat.Prefix("    ")))
}
        
    

Let's walk through this program line-by-line and we'll see some of the principles that make gonum work! Once you understand these, then you'll get an intuition for how to use other parts of the gonum module.

Creating the Matrices

First, let's construct the matrices.

        
// Define the matrix A
A := mat.NewDense(4, 3, []float64{
    5.0e4, 1.0e4, 6.0e4,
    2.0e4, 2.5e4, 1.0e3,
    5.0e4, 5.0e4, 1.5e3,
    1.5e5, 1.5e4, 0.0,
})

// Define the vector b
b := mat.NewVecDense(4, []float64{
    1.1e6,
    2.0e5,
    5.0e5,
    3.5e5,
})

        
    

Simply put, here we construct the matrix \( A \) and the vector \( b \) from the details given in the problem statement. In this case, we create them as dense matrices and vectors.

Dense matrices and vectors are a technical term. They signify matrices and vectors where we define every single element of the matrix or vector. (In contrast, in some other parts of the gonum module, you can specify sparse matrices, which you do not need to completely define.)

Now, we will solve the least-squares problem.

Solving the Least-Squares Problem

To solve the least-squares problem, we will use the SolveVec function. This function is defined in the mat package of gonum.

        
// Solve the least-squares problem
var x mat.VecDense
x.SolveVec(A, b)
        
    

The SolveVec function takes two arguments: the matrix \( A \) and the vector \( b \). It returns a vector \( x \) that is the solution to the least-squares problem.

Note that we must define the vector \( x \) as a mat.VecDense object before we call the SolveVec function.

This is because the SolveVec function will modify the vector \( x \) in-place. This is a common pattern in gonum and you will see this pattern a lot in the gonum source code.

Now, let's see how well the solution solves the problem (i.e., how far away are the expected incomes from the targets).

Evaluating the Solution

To evaluate the solution, we will use the MulVec and SubVec functions. This function is defined in the mat package of gonum.

        
// Print the solution
fmt.Printf("x = %.4f\n", mat.Formatted(&x, mat.Prefix("    ")))

// Print the solution
fmt.Printf("x = %.4f\n", mat.Formatted(&x, mat.Prefix("    ")))

// Print the product
var prod mat.VecDense
prod.MulVec(A, &x)

fmt.Printf("A*x = %.4f\n", mat.Formatted(&prod, mat.Prefix("    ")))

    // Print the residual
var res mat.VecDense
res.SubVec(b, &prod)

fmt.Printf("b - A*x = %.4f\n", mat.Formatted(&res, mat.Prefix("    ")))
        
    

Again, we'll need to define the vectors \( prod \) and \( res \) as mat.VecDense objects (the pattern returns). Also, notice that in each case the receiver of the function (i.e., prod in prod.MulVec(A, &x)) is the "result" of the function and the functions two arguments are the other "inputs" to the function.

Without going into too much detail, the code prod.MulVec(A, &x) computes the product \( A x \) and the code res.SubVec(b, &prod) computes the difference \( b - A x \).

For the sake of completeness, in this program we look at the residual which describes how far away our choice of costs gets us to the target incomes. How close to zero is the residual?

Conclusion

In this post, we have seen how to use gonum to solve a simple least-squares problem. We have seen how to construct matrices and vectors, how to solve the least-squares problem, and how to evaluate the solution.

If you have any questions, please feel free to reach out to me via email or ask a question on the gonum repository.

Acknowledgements

A big thanks to all of the contributors to gonum who enable a lot of scientific computing in Go! I've seen some of their work and thought it was so excellent that I incorporate it into MatProGo.