Scientific computing: a Rust adventure [Part 0 - Vectors]

February 2019 · 22 minute read

It wasn’t always so clear, but the Rust programming language is fundamentally about empowerment: no matter what kind of code you are writing now, Rust empowers you to reach farther, to program with confidence in a wider variety of domains than you did before.
[Foreword - The Rust Programming Book]

A personal foreword

My daily work revolves around building Machine Learning applications, while a lot of my evenings have been spent experimenting with Rust, getting more and more fascinated and in love with the language.

It couldn’t be helped: I started to have a look at what the Rust ecosystem had to offer for Machine Learning, Big Data and scientific computing at large. I quickly found out that there is a lot to be done and a lot of potential (see here or here). It got me really fired up 🔥

Looking back at the Python world, one of the many reasons behind the explosion of the scientific ecosystem was the availability of a small pocket of solid foundational libraries - the SciPy stack: NumPy, SciPy, Pandas.
NumPy and SciPy, in particular, provided at the same time a reliable toolset to get work done and a concrete promise of interoperability: most projects are built, directly or indirectly, on NumPy’s arrays or SciPy’s sparse arrays (e.g. Scikit-learn, Pandas, PyMC3) - an informal “interface” you can rely on, in line with Python’s duck typing philosophy.

There is no reason to think that things are going to be any different for Rust: building a handful of solid core crates will go a long way in promoting the growth of a healthy ecosystem.

That’s how I came in contact with ndarray (GitHub), a fantastic Rust crate written by bluss and maintained together with jturner314.
You can think of ndarray as Rust’s equivalent of NumPy: it provides a data structure to represent and manipulate n-dimensional arrays.

I started contributing with a bunch of PRs that eventually evolved into a stand-alone crate, ndarray-stats, now close to its second release.
At the same time I have been looking at ndarray’s GitHub issues as well as reading discussions from the #science-and-ai channel on Discord: yes, we have to build new capabilities and new crates, but at the same time we really need a good introduction to the existing tools to get more hands on board and increase the reach of the good crates we already have.
I guess we all agree that Rust wouldn’t be where it currently is without The Rust Programming Book and the awesome work done on documentation.

This series of blog posts aims to do, on a smaller scale, something similar for ndarray: I want to provide a gentle, problem-based introduction, speaking to people with a background in scientific computing who have just got into Rust (let’s say they have skimmed most of the book once) and want to have a look at what is possible, and how.
I will quickly recall some core Rust concepts, but this is not going to be by any means a complete introduction to Rust as a language - the book is and remain the best reference to get started with Rust.

This foreword was meant to be short and concise, but whatevs, context is important. Now, let’s get started for real.

Scientific Computing: A Rust adventure (TOC)

Vectors

We won’t actually be touching n-dimensional arrays in this first post ¯\_(ツ)_/¯
We will instead spend some time to get familiar with their one-dimensional counterparts: Vec<T>, vectors.
Vectors are an easier starting point and we can use them to briefly recall a lot of concepts that will become extremely important when we get to the multidimensional case.

All code examples will have a link to the Rust Playground, so that you can compile and run all code snippets even if you don’t have Rust installed on your machine (or you are reading the article from your smartphone!).

Vec<T> is part of Rust’s standard library: vectors are used to store several values of the same type, T, in a single data structure. All values are placed next to each other in memory and can be accessed using indexes (Rust Playground):

1// A vector of unsigned integers, created using the vec! macro
2let v: Vec<usize> = vec![1, 2, 3];
3// Its first element (indexing is 0-based, like Python)
4// `first` has type `usize` (unsigned integer)
5// No need to type annotate it, it's inferred by the compiler
6// given that v is a vector with usize elements
7let first = v[0];
8assert_eq!(first, 1);
9

What manipulation can we do using a vector?
First of all, we might be interested in getting a subset of its elements, a slice.

This can be achieved with a syntax similar to indexing (RP):

1let v: Vec<usize> = vec![1, 2, 3, 4, 5, 6];
2// Let's take the first three elements
3// 0..3 is a range - indexes from 0 (included) to 3 (not included)
4let first_half = &v[0..3];
5assert_eq!(first_half.len(), 3);
6

You might have noticed the ampersand & in front v[0..3]. Slices are nothing more than references to a subset of vector elements.

References play an important role in Rust programming due to the concept of ownership: each value has a variable that’s called its owner and there can only be one owner at a time.
Ownership is a key language feature and it allows the Rust compiler to take care of allocation and de-allocation of memory without requiring user intervention or garbage collection (as in Python, Go or Julia). The ownership system has far-reaching consequences, including memory safety and fearless concurrency - we’ll learn more as we get deeper into the language.

In our case, v owns its elements. In Rust assigning a value to a variable using = moves its ownership1 to the variable (it’s called move semantics).
It is not allowed to move values out of a structure or a collection, such as a vector. You can either move the ownership of the whole vector over to a new variable (e.g. let w = v;) or you can take a reference, using the & symbol: the vector retains ownership, but you get access to some of its elements as long as the original owner is in scope.

If we had wanted to copy over the whole vector, we could have used the clone method (RP):

 1// A new vector
 2let v = vec![0, 1, 2, 3];
 3// A reference to v - v's elements are not copied
 4let reference = &v;
 5// v's elements are copied over
 6// We now have two copies of v's elements around:
 7// v is the owner of the first set of elements,
 8// cloned owns the second set of elements.
 9// There is no ownership relationship of any kind
10// between v's and cloned's elements.
11let cloned = v.clone();
12// Cloning a reference gives us a new reference
13// It does not create a copy of v
14let cloned_reference = reference.clone();
15

Let’s have a look at what we can do using a vector.

A naive scalar product

As we said before, a vector is just a data structure: it is not a 1-dimensional array in the mathematical sense. The best analogy is probably a list in Python.
It is thus educational to try to implement a very naive version of scalar product:
$$
v\cdot w=\sum_{i=0}^{n-1} v_i w_i
$$ where $v$ and $w$ are vectors with $n$ elements.

Let’s start with vectors of signed integers. A first attempt to write down the function signature could roughly look like this:

1fn scalar_product(v: Vec<i32>, w: Vec<i32>) -> i32
2

We take two vectors as input, both containing signed integers (i32), and we return an integer as output. Let’s flesh it out (RP).

 1fn scalar_product(v: Vec<i32>, w: Vec<i32>) -> i32
 2{
 3    // The number of elements we are dealing with
 4    let length = v.len();
 5    // We expect both vectors to have the same number of elements
 6    assert_eq!(length, w.len());
 7
 8    // We initialise a new variable to 0
 9    // We will use it to accumulate the scalar product
10    let sum = 0;
11    // We iterate over the range of indexes
12    for index in 0..length {
13        // We add the product of elements corresponding to the same index
14        sum = sum + v[index] * w[index];
15    }
16    // We return the scalar product
17    // The final expression in a function is used as return value in Rust
18    sum
19}
20

Using an assertion is not an idiomatic way to handle errors, but we can live with it for now.
All in all, this implementation looks legit. If you try to compile it, though, you will get an error.

110 |let sum = 0;
2   |    ---
3   |    |
4   |    first assignment to `sum`
5   |    help: make this binding mutable: `mut sum`
6...
713 |    sum = sum + v[index] * w[index];
8   |    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
9   |    cannot assign twice to immutable variable

This is an important point: variables in Rust are immutable by default. You have to explicitly mark them with the mut keyword to make them mutable.
As the compiler suggests (and where would we be, without its suggestions?), we can quickly fix it using let mut sum = 0 instead of let sum = 0.

Done, now it compiles2.
Let’s try it out! (RP)

1let v = vec![1, 1, 1];
2let w = vec![1, 0, 1];
3let prod = scalar_product(v, w);
4assert_eq!(prod, 2); // Spoiler: it works
5

Cool, everything is working out smoothly. Let’s tweak the problem statement a little bit: we fix a vector, v, and we want to find out which one among two other vectors, let’s say w and z, gives us the highest number when we take its scalar product against v.
It sounds straight-forward: (RP)

 1let v = vec![0, 1, 2];
 2let w = vec![1, 2, 3];
 3let z = vec![-1, 0, 6];
 4
 5let v_dot_w = scalar_product(v, w);
 6let v_dot_z = scalar_product(v, z);
 7
 8if v_dot_w > v_dot_z {
 9    println!("The solution is w!")
10} else {
11    if v_dot_w < v_dot_z {
12        println!("The solution is z!")
13    } else {
14        println!("There is not difference")
15    }
16}
17

We launch it and… 🚨 the compiler is not happy! 🚨

1error[E0382]: use of moved value: `v`
2  --> src/main.rs:25:34
3   |
424 |     let v_dot_w = scalar_product(v, w);
5   |                                  - value moved here
625 |     let v_dot_z = scalar_product(v, z);
7   |                                  ^ value used here after move

To understand what is happening we need to take a step back and look again at our scalar_product’s function signature:

1fn scalar_product(v: Vec<i32>, w: Vec<i32>) -> i32
2

We are taking both arguments, v and w, by value - this triggers the move semantics we mentioned a couple of paragraphs above (with the = operator): we are giving scalar_product ownership over its arguments, v and w. This means that after calling scalar_product the first time, we can no longer use v and w in our main function.
When we try to access v a second time, to compute its product with z, we are de facto trying to access something that is no longer there!
And the compiler, patient as it always is, does not fail to remind us of our misstep.
How do we fix it?

We could change

1let v_dot_w = scalar_product(v, w)
2

into

1let v_dot_w = scalar_product(v.clone(), w)
2

This would create a copy of v on the fly using clone and we would be giving scalar_product ownership over that copy, instead of v itself. v can be thus consumed in the next function call to compute v_dot_z.
It’s easy, and it works, but it’s wasteful: scalar_product does not need to own its arguments, it’s enough for scalar_product to have read-only access to their elements to do its job!

We can change scalar_product’s signature to avoid triggering Rust’s move semantics.
As we said before, we use references when we need to access elements without taking ownership. We can rephrase scalar_product as follows:

1fn scalar_product(v: &Vec<i32>, w: &Vec<i32>) -> i32
2

Nothing has to be changed into the function body, but we do need to tweak our main program to pass in references instead of values: (RP)

1[...]
2// A couple of ampersands to rule them all!
3let v_dot_w = scalar_product(&v, &w);
4let v_dot_z = scalar_product(&v, &z);
5[...]
6

Painless, it works and we don’t have to copy data around.
It’s not perfect yet. What happens if we try to call scalar_product on a vector slice?
Something like this perhaps: (RP)

1let v = vec![0, 1, 2, 3, 4, 5];
2let head = &v[0..3];
3let tail = &v[3..6];
4let head_vs_tail = scalar_product(head, tail);
5

The compiler is not pleased

1error[E0308]: mismatched types
2  --> src/main.rs:23:39
3   |
423 |     let head_vs_tail = scalar_product(head, tail);
5   |                                       ^^^^
6   |      expected struct `std::vec::Vec`, found slice
7   |
8   = note: expected type `&std::vec::Vec<i32>`
9              found type `&[{integer}]`

How do we read this? We saw that Vec<i32> is a vector of integers, but I didn’t tell you what type is used to represent slices: it’s &[i32]. In general, from Vec<T> we can take slices of type &[T].
Armed with this knowledge we can interpret what the compiler is complaining about: scalar_product expects a vector of integers, not a slice!
Once again, we can solve this conondrum by changing scalar_product’s function signature: (RP)

1fn scalar_product(v: &[i32], w: &[i32]) -> i32
2

No change is required to the function body.

But how would we deal with vectors now? Do we need to keep two separate versions of scalar_product to cope with both vectors and slices?

Absolutely not.
We can use the as_slice method: given a vector, it returns a slice pointing to all its elements - the only difference being that the slice does not have ownership of the elements it points to. (RP)

1let v = vec![1, 1, 1];
2let w = vec![1, 0, 1];
3let prod = scalar_product(v.as_slice(), w.as_slice());
4assert_eq!(prod, 2);
5

Things are actually simpler thanks to some syntactic sugar built into the Rust compiler: a reference to a Vec<i32> will be automatically interpreted as a slice if required. (RP)

1let v = vec![1, 1, 1];
2let w = vec![1, 0, 1];
3// All good!
4let prod = scalar_product(&v, &w);
5assert_eq!(prod, 2);
6

What are the performance of this uber-naive scalar product implementation?
Let’s run a quick and dirty benchmark3 computing the scalar product of two random vector of 100000 integers, magnitude between -1000 and +1000.

Language Time (us) Notes
Python 1720 Same function, using Python’s lists
Python 12.2 NumPy, using np.array - (v * w).sum()
Python 6.1 NumPy, using np.array - v.dot(w)
Rust 1.8 Super naive

It’s roughly ~3 times faster than NumPy’s dot method.
Is it starting to get interesting? We have just scratched the surface.

Bubble sort

Let’s take a look at another simple algorithm we can implement to play around a little more with vectors: sorting.
Truth be told, the standard library includes a sort method for vectors, but it’s useful to do it on our own for educational purposes.

Speed is not our main concern so we will roll with a very simple bubble sort.
Let’s put down some rules: let’s say that, given a vector of integers, we want to sort it in increasing order without allocating a copy of its elements.

We can use move semantics for the argument of our function!
Let’s sketch it out: (RP)

 1fn bubble_sort(v: Vec<i32>) -> Vec<i32> {
 2    let length = v.len();
 3
 4    for _ in 0..length {
 5        for j in 0..(length - 1) {
 6            // If two consecutive elements are not properly ordered
 7            if v[j] > v[j+1]  {
 8                // We swap their values
 9                let temp = v[j+1];
10                v[j+1] = v[j];
11                v[j] = temp;
12            }
13        }
14    }
15    // We return the sorted array
16    v
17}
18

We try to compile it and… 💥! Error!

1error[E0596]: cannot borrow `v` as mutable, as it is not declared as mutable
2  --> src/main.rs:10:17
3   |
41  | fn bubble_sort(v: Vec<i32>) -> Vec<i32> {
5   |                - help: consider changing this to be mutable: `mut v`
6...
710 |                 v[j+1] = v[j];
8   |                 ^ cannot borrow as mutable

Remember when I said that variables are immutable by default? That holds as well for function arguments. If we want to mutate them, then we need to add the mut keyword in front of them. The signature thus becomes (RP)

1fn bubble_sort(mut v: Vec<i32>) -> Vec<i32> {
2

and the compiler is happy.
We can check out that our implementation works with a very simple example: (RP)

1let v = vec![6, 5, 4, 3];
2assert_eq!(bubble_sort(v), vec![3, 4, 5, 6]); // It works!
3

Let’s take it one step further: we only want to sort the first half of a vector, leaving the rest untouched. A first approach would be to get the first half as a slice and pass it over to bubble_sort: (RP)

1let v = vec![6, 5, 4, 3];
2let midpoint_index = v.len() / 2;
3// We must take a reference
4// We can't move elements out of their owner, v
5let first_half = &v[0..midpoint_index];
6assert_eq!(bubble_sort(first_half), vec![5, 6, 4, 3]);
7

Do I even have to say it? Rust compiler said no 😞

 1error[E0308]: mismatched types
 2  --> src/main.rs:22:28
 3   |
 422 |     assert_eq!(bubble_sort(first_half), vec![5, 6, 4, 3]);
 5   |                            ^^^^^^^^^^
 6   |                            |
 7   |     expected struct `std::vec::Vec`, found reference
 8   |     help: try using a conversion method: `first_half.to_vec()`
 9   |
10   = note: expected type `std::vec::Vec<i32>`
11              found type `&[{integer}]`

Its objection is reasonable: our bubble_sort function takes a Vec<i32> as input, not a slice.
We could use the to_vec method, as the compiler suggests: this would clone the elements first_half is pointing to, creating a new vector and passing its ownership to bubble_sort.

But we said that we do not want to allocate new memory ⛔
I said what I said

Let’s try to change the function signature of bubble_sort to work with slices: (RP)

1fn bubble_sort(v: &[i32]) -> &[i32]
2

If we try to compile it, this is what we get:

 1error[E0594]: cannot assign to `v[..]` which is behind a `&` reference
 2  --> src/main.rs:10:17
 3   |
 41  | fn bubble_sort(v: &[i32]) -> &[i32] {
 5   |                   ------ help:
 6   | consider changing this to be a mutable reference: `&mut [i32]`
 7...
 810 |                 v[j+1] = v[j];
 9   |                 ^^^^^^^^^^^^^
10   | `v` is a `&` reference, so the data it refers to cannot be written

References, just like variables, are immutable by default: they cannot mutate the values they point to!

This brings us to last piece in the puzzle of ownership: mutable references.
We can take a mutable reference to a value using the &mut prefix: (RP)

 1// The owner of the elements, v, must be marked as mutable
 2// to allow for a mutable reference
 3let mut v = vec![0, 1, 2];
 4// We take a mutable reference using &mut
 5let mutable_reference = &mut v;
 6// We change one of the entries using the mutable reference...
 7mutable_reference[0] = -1;
 8// ...and the original vector is modified!
 9assert_eq!(v, vec![-1, 1, 2]);
10

Mutable references are powerful, but with great power comes great responsability 🕸️
At any point in time, there can only be one active mutable reference pointing to a certain value: (RP)

1let mut v = vec![0, 1, 2];
2let first_mutable_reference = &mut v; // All good
3let second_mutable_reference = &mut v; // Compiler error
4second_mutable_reference[1] = 5;
5first_mutable_reference[0] = 7;
6

Why is that?
Well, it’s quite undesirable to work with second_mutable_reference, assuming that certain invariants and properties hold, just to find out that somewhere else first_mutable_reference was used to modify the data you were working with behind your back - that would be a recipe for disaster, wouldn’t it? (data races, anyone?)

You don’t have this problem when working with immutable reference, using &. (RP)

1let v = vec![0, 1, 2];
2let first_immutable_reference = &v; // All good
3let second_immutable_reference = &v; // All good too!
4println!("My first immutable reference: {:?}", first_immutable_reference);
5println!("My second immutable reference: {:?}", second_immutable_reference);
6

How do mutable and immutable references interact with each other? Let’s look at an example: (RP)

1// v is marked as mutable, to allow a mutable reference
2let mut v = vec![0, 1, 2];
3let immutable_reference = &v;
4let mutable_reference = &mut v;
5// We modify v using mutable_reference
6mutable_reference[0] = 1;
7// We print v using an immutable reference
8println!("My vector: {:?}", immutable_reference);
9

The compiler bails:

 1error[E0502]: cannot borrow `v` as mutable because
 2it is also borrowed as immutable
 3 --> src/main.rs:5:29
 4  |
 54 |     let immutable_reference = &v;
 6  |                               --
 7  |                immutable borrow occurs here
 8  |
 95 |     let mutable_reference = &mut v;
10  |                             ^^^^^^
11  |                   mutable borrow occurs here
12...
139 |     println!("My vector: {:?}", immutable_reference);
14  |                                 -------------------
15  |                           immutable borrow later used here

The reasoning goes along the same lines of what we discussed when talking about multiple mutable references: if you have an immutable reference, you don’t want the values you are referring (or borrowing, using Rust’s terminology) to be changed under your nose by a different part of the code using a mutable reference.
At any point in time, you can have either4:
- exactly one active (=in scope) mutable reference;
- an infinite number of immutable references.

Once again, the Ownership chapter of the book is the best place to go to get a deeper understanding of the ownership system.

Let’s grab the power of mutable references to write down the (final?) version of our bubble_sort function: (RP)

 1fn bubble_sort(v: &mut [i32]) {
 2    let length = v.len();
 3
 4    for _ in 0..length {
 5        for j in 0..(length - 1) {
 6            // If two consecutive elements are not properly ordered
 7            if v[j] > v[j+1]  {
 8                // We swap their values
 9                let temp = v[j+1];
10                v[j+1] = v[j];
11                v[j] = temp;
12            }
13        }
14    }
15}
16

Given that we are modifying the input in place there is not need to return the slice reference at all.
It compiles! And our little example, with some tweaks to adapt to the new signature of bubble_sort, (RP)

 1// v is marked as mut, to allow for mutable references
 2let mut v = vec![6, 5, 4, 3];
 3let midpoint_index = v.len() / 2;
 4let first_half = &mut v[0..midpoint_index];
 5// bubble_sort modifies first_half in place,
 6// no need for the function to return an output
 7bubble_sort(first_half);
 8// But first_half is a mutable reference to a subset of v
 9// We can thus check that v itself has been modified
10// as expected!
11assert_eq!(v, vec![5, 6, 4, 3]);
12

runs without any issues.

Sorting is a good example of the two language problems in Python: we could use a NumPy array to hold our data, but unless we rely on built-in methods (such as np.sort), we cannot implement a fast sorting algorithm working in the Python world (because for cycles are veeeery slow!). We need to do it in C and then call the C function from Python itself (which is what NumPy does under the hood).
In Rust, instead, we can express the bubble sort algorithm using a syntax that is very similar to a high-level language without having to compromise on performance. The best of both worlds?

Wrapping up

Hopefully, this brief tour of vectors has provided a simple but clear introduction to Rust and its peculiarities.

If you had already been exposed to Rust (by skimming the book or by experimenting with it a little bit) you have probably mumbled “I know this already!” at multiple points in the article, but ownership takes a while to get used to and I believe it’s beneficial to have a closer look at what it means when working with vectors.

We will draw a lot of analogies in the next episode between Vec and ArrayBase, ndarray’s data structure for n-dimensional arrays. It will allow us to introduce the complexity required to deal with multiple dimensions without losing sight of what we are trying to accomplish and how all the pieces are supposed to work together.

Notes, suggestions and feedbacks to improve this introduction are more than welcome. You can reach me either on LinkedIn or on Twitter (@algo_luca - we want more beautiful people to try Rust out for scientific computing, don’t we?).

Scientific Computing: A Rust adventure (TOC)

Appendix

A quick table to sum up what we said about ownership and vectors:

Type Name Notes
Vec<i32> Vector It owns its elements and it can be modified (using the mut keyword).
&Vec<i32> Immutable vector reference It points to a Vec<i32>. It does not own its elements and it cannot modify the vector it points to.
&[i32] Immutable vector slice It points to a subset of a Vec<i32>. An immutable vector slice does not own its elements and it cannot modify the elements it points to.
&mut Vec<i32> Mutable vector reference It points to a Vec<i32>. A mutable vector reference does not own its elements but it can modify the elements of the vector it points to. There can only be at most one active mutable vector reference around at any point in time and the original vector variable has to be mut to allow a mutable reference.
&mut [i32] Mutable vector slice It points to a subset of a Vec<i32>. A mutable vector slice does not own its elements but it can modify the elements it points to. There can only be at most one active mutable slice around at any point in time and the original vector variable has to be mut to allow mutable slices.

Now that we have a complete picture, it makes sense to point out that you seldomly have a reason to take &Vec<i32> as type argument for a function: it doesn’t allow you to accept slices, while using &[i32] you get both vectors and slices, as we showed in some of the code snippets above (using as_slice or just & + syntactic sugar).


  1. Unless the value is of a type that implements the Copy trait: instead of being moved, its value is copied over. This is the case, for instance, for all primitive numeric types (signed/unsigned integers, floats, etc.). We will discuss traits (and generics) in the article - bear with me for the time being.
    [return]
  2. A quick observation: apart from the function signature, we didn’t have to annotate with a type any of the variables we declared (e.g. length, sum or index). Even though Rust is statically typed, the compiler is usually able to infer the right type for all your variables starting from the signatures of the functions you are using. You will seldom need to type-annotate individual variables in your code.
    [return]
  3. The code used to measure Rust’s execution time can be found here - using version 0.2 of criterion on a AMD Ryzen 7 1700. Python’s benchmarks have been computed using %timeit.
    [return]
  4. This is not 100% true, given non-lexical lifetimes, an improvement to the compiler that has been recently released. The compiler will not complain if it can determine that the actual usage scopes of your mutable and immutable references are not going to overlap. It makes your life easier, without compromising Rust’s memory guarantees.
    [return]