Hello World of Programming with Linear Algebra

Imagine this simplified code for inventory modeling. (When I say "simplified" i mean really simplified. We're using floats for prices (bad), we store the data in global state, the architecture is far from even a simple web application. But the model is familiar enough to a typical software developer.)

(def products {:banana {:price 1.3 :id :banana}
               :mango {:price 2.0 :id :mango}
               :pineapple {:price 1.9 :id :pineapple}
               :pears {:price 1.8 :id :pears}})

This (imaginary) application exists to track sales. A customer puts the desired products into a cart, we calculate the total price, and, later, perform the delivery. Each cart only stores the products' identifiers and quantities (in unspecified units; yes, it's super-simplified).

(def cart1 {:banana 10
            :pineapple 7
            :pears 3})

(def cart2 {:pineapple 3
            :mango 9})

Having defined product and cart data, we write a function that, given the products "database" and a cart, calculates the total price of the products in the cart. The cart-price function reduces all [product quantity] pairs in the cart, by retrieving the appropriate product map in the product-db, and taking the value associated with its :price key. It multiplies that price with the quantity, and accumulates it in total.

(defn cart-price [product-db cart]
  (reduce (fn [total [product quantity]]
            (+ total (* (:price (product-db product)) quantity)))
          0
          cart))

Let's call this function with the available carts, and see it in action.

(cart-price products cart1)

=> 31.699999999999996

(cart-price products cart2)

=> 23.7

We hopefully have more than one order. Our code can easily process sequences of carts, and compute the total revenue.

(reduce + (map (partial cart-price products) [cart1 cart2]))

=> 55.39999999999999

It's all good; but what does it have to do with linear algebra?

In the previous implementation, we entangled the specifics of data storage and the algorithm that computes the total price. In this simple model, it's not much of a problem, but if the data model is more complex, and the algorithm not as simple as the straightforward map/reduce, this quickly leads to (at least) two problems:

• code becomes too complicated
• program performance degrades quickly

Let's first tackle the code complexity by extracting the computation logic from the domain into the abstract mathematical notion of vectors and operations on these vectors. In this particular example, vectors help us encapsulate a bunch of numbers as one atomic unit.

(def product-prices [1.3 2.0 1.9 1.8])
(def cart-vec-1 [10 0 7 3])
(def cart-vec-2 [0 9 3 0])

We recognize that the logic we've already developed for computing the total price matches a simple and well known mathematical operation, known as the dot product, a scalar product of two vectors.

(defn dot-product-vec [xs ys]
  (reduce + (map * xs ys)))

Given two vectors, [1 2 3] and [4 5 6], the dot product computes one number, a scalar, that represent a scalar product of these two vectors. Right now, we don't even care about theoretical details of the dot product; we recognize that it technically computes the same thing that we need in our domain, and it seems useful. There are other ways to multiply vectors, which return non-scalar structures.

(dot-product-vec [1 2 3] [4 5 6])

=> 32

We can see that, when applied to the vectors holding product prices and quantities, it returns the correct results that we've already seen.

(dot-product-vec product-prices cart-vec-1)

=> 31.699999999999996

(dot-product-vec product-prices cart-vec-2)

=> 23.7

Getting the total price requires another map/reduce, but we will quickly see that this, too, can be generalized.

(reduce + (map (partial dot-product-vec product-prices)
               [cart-vec-1 cart-vec-2]))

=> 55.39999999999999

When we abstract away the specifics of the domain, we end up with a number of general operations that can be reused over and over, and combined into more complex, but still general, operations. Countless such operations have been studied and theoretically developed by various branches of mathematics and related applied disciplines for a long time. What's more, many have been implemented and optimized for popular hardware and software ecosystems, so our main task is to learn how to apply that vast resource to the specific domain problems.

Linear algebra is particularly well supported in implementations. Whenever we need to process arrays of numbers, it is likely that at least some part of this processing, if not all of it, can be described through vector, matrix, or tensor operations.

Let's introduce a bit more complication. Say that we want to support different discounts for each product, in the form of multipliers. That gets us the price reductions, that we should subtract from the price. An alternative way, shown in the following snippets is to subtract the discount coefficients from 1.0 to get the direct multiplier that gets us to the reduced price.

(def discounts (dv [0.07 0 0.33 0.25]))
(def ones (entry! (dv 4) 1))

We can subtract two vectors by the axpy function. axpy stands for "scalar a times x plus y".

(axpy -1 discounts ones)

=>
#RealBlockVector[double, n:4, offset: 0, stride:1]
[   0.93    1.00    0.67    0.75 ]

The mul function multiplies its vector, matrix, or tensor arguments element-wise, entry by entry.

(mul (axpy -1 discounts ones) product-prices)

=>
#RealBlockVector[double, n:4, offset: 0, stride:1]
[   1.21    2.00    1.27    1.35 ]

The following code seamlessly incorporates this new part of the algorithm into the implementation that we already have.

(asum (mv carts (mul (axpy -1 discounts ones) product-prices)))

=> 46.87

Why stop here? Suppose that we want to simulate the effects of multiple discount combinations on the total price. As earlier, we put all these hypothetical discount vectors into a matrix, in this case, three samples that we'd like to investigate.

(def discount-mat (dge 4 3 [0.07 0 0.33 0.25
                            0.05 0.30 0 0.1
                            0 0 0.20 0.40]))

=>
#RealGEMatrix[double, mxn:4x3, layout:column, offset:0]
▥       ↓       ↓       ↓       ┓
→       0.07    0.05    0.00
→       0.00    0.30    0.00
→       0.33    0.00    0.20
→       0.25    0.10    0.40
┗                               ┛

We have to subtract these numbers from 1.0. Instead of populating the matrix with 1.0, we will demonstrate the outer product operation, implemented by the function rk. Given two vectors, it produces a matrix that holds all combinations of the product of the entries of the vectors.

(rk ones (subvector ones 0 3))

=>
#RealGEMatrix[double, mxn:4x3, layout:column, offset:0]
▥       ↓       ↓       ↓       ┓
→       1.00    1.00    1.00
→       1.00    1.00    1.00
→       1.00    1.00    1.00
→       1.00    1.00    1.00
┗                               ┛

We can also utilize rk to "lift" product prices vector to a matrix whose shape matches the shape of the discount combinations matrix.

(def ones-3 (subvector ones 0 3))
(def discounted-prices (mul (axpy -1 discount-mat (rk ones ones-3))
                            (rk product-prices ones-3)))

The result is a matrix of hypothetical discount prices that we'd like to simulate.

=>
#RealGEMatrix[double, mxn:4x3, layout:column, offset:0]
▥       ↓       ↓       ↓       ┓
→       1.21    1.23    1.30
→       2.00    1.40    2.00
→       1.27    1.90    1.52
→       1.35    1.62    1.08
┗                               ┛

Now, the most interesting part: how do we calculate the totals from this matrix and the matrix of carts we've produced earlier. (Not) surprisingly, just a single operation, matrix multiplication, completes this task!

(mm carts discounted-prices)

=>
#RealGEMatrix[double, mxn:2x3, layout:column, offset:0]
▥       ↓       ↓       ↓       ┓
→      25.05   30.51   26.88
→      21.82   18.30   22.56
┗                               ┛

Now we only need to sum the columns up to get the three final totals. We won't do this column-by-column. Instead, we'll use the "mv with ones" approach we've already encountered. Note that we need to transpose the matrix to match the desired structure.

(trans (mm carts discounted-prices))

=>
#RealGEMatrix[double, mxn:3x2, layout:row, offset:0]
   ▤       ↓       ↓       ┓
   →      25.05   21.82
   →      30.51   18.30
   →      26.88   22.56
   ┗                       ┛

And, the final answer is…

(mv (trans (mm carts discounted-prices)) (subvector ones 0 2))

=>
#RealBlockVector[double, n:3, offset: 0, stride:1]
[  46.87   48.81   49.44 ]

Given three (or three million) possible discount combinations, we get a vector of the total revenue amounts. Of course, being a toy example, this code doesn't take into account that lower prices would (likely) induce more sales; let's not carry a Hello World example too far.

So, the first major benefit of using a library, such as Neanderthal (), based on many decades of numerical computing research and development is that we have access to a large treasure trove of useful, well thought, general functions for developing useful, general or customized, number processing algorithms.

Another major improvement is performance. Although toy examples may be implemented in any way you'd like, and they'd still work reasonably well, real-world data processing almost always involves either many data points, or many computation steps, or, often – both.

I'm not talking about a couple dozen percentages, but improvements of many orders of magnitude. But that's a story that has to be looked at in more depth. I've written two books where I go into much, much, more depth and width on this. You can also check some earlier articles on this blog; there's many examples that demonstrate how fast this approach is!