Introduction

I am writing this blog for two main reasons. First, I am using Rust and I would love to see it used more in the scientific community. However, the lack of documentation and guidance toward scientific crates is I think refraining more widespread adoption. So, let’s be the change I wanted to see! Second, by writing small documentions about scientific crates and how they interact together I am hoping to discover new crates and to familiarize myself with writing scientific code in Rust. I am usually using Python to prototype everything but it may be good to be able to do it sometime in Rust also!

In this first post I will describe the basic need we have when beginning to do scientific computation: manipulating data. Usually, when I am in Python I am using numpy which is the standard way to do it. Most Python module made for science are also using numpy which is really convenient to glue everything together. In Rust we do not have (yet!) a standard numerical crate. The reason is that the language is still young and many want to experiment with their own implementations. It is for me rather unfortunate but this is the current reality. The upside is that we have specialized crates tailored for our needs! I found two main crates that I think are of high qualities and answer different user cases: ndarray and nalgebra. ndarray is the equivalent of numpy while nalgebra aims to provide tools to manipulates matrices and vectors. I wanted to present both crates in one post but it would have resulted in a very long post. Instead my first post is about nalgebra and I will post a second one on ndarray later.

Please note that as I could obviously not describe all the different crates in this post I focused only on nalgebra and ndarray. It does not mean that I think the other crates are not worthwhile but simply I needed to make a choice and I presented what I thought was the most popular choices.

nalgebra: When you only want matrices and vectors

Having arrays with multiple dimensions is common but many times manipulating vectors and matrices is the only thing you want. In these situations the nalgebra crate is perfect for you. This crate has been created with the underlying goal to be used for general linear algebra operations and computer graphics. As such, this crate is really optimal for manipulating matrices and vector representing geometrical transformations. This crate has a really beautiful documentation and reading it is necessary to discover all the features the crate provides. However, I feel like when you are discovering the crate coming from a Python world (or other programming languages less strict than Rust) it is a bit confusing at first. How do you create and manipulate matrices? How do you use the factorizations explained in the documentaion? This post aims to give an easy start for beginners wanting to try nalgebra.

Vector/Matrix types

Creating a matrix or a vector in nalgebra is not difficult but we need to be clear on what we want to create before doing it. The most significant choice is if the matrix or vector that we want to create has a column/row size determined in compile time. The more information you give to the compiler and the more your code could be optimized. In Rust, a lot of things are done during compilation and, as such, the code we write has to be correct (which is a bothering point for lot of Rust newcomers). The remaining information needed is the type of the elements inside the matrix.

Because there are many possible type for your matrix, you have to specify manually which type you want to use. For the smallest dimension nalgebra got you covered. Indeed, there are inside the crate many alias made for you to avoid specifying the types of matrix with dimensions lesser than 6. For example, here is how to create two matrices of size 2x4 and 5x6.

extern crate nalgebra as na;
use na::{Matrix2x4, Matrix5x6};

let m1 = Matrix2x4::new(2.3, 10.9, -5.7, 0.0,
                        1.5, -1.4, 0.48, 8.9);
let m2 = Matrix5x6::new(1, 2, 3, -3, -2, -1,
                        4, 5, 6, -6, -5, -4,
                        7, 8, 9, -9, -8, -7,
                        4, 5, 6, -6, -5, -4,
                        1, 2, 3, -3, -2, -1);

And here is how to create two vectors of size 2 and 4.

extern crate nalgebra as na;
use na::{Vector2, Vector4};

let v1 = Vector2::new(1.4, -1.0);
let v2 = Vector4::new(-2, 4, 5, 0);

Note that these matrices and vectors have their dimension known at compile time. There is no alias for a matrix of 4 rows and not known at compile time column number for example.

What about bigger vector and matrices? Well you will have to define their types. And remember that a vector is just a matrix with one row or one column. Creating a new matrix type is not difficult, we need to specify the type of the data we want to store (i32, f64…), the number of rows and columns, and finally giving the buffer that will contain all the matrix components.

However there is one catch. As of today (August 4th, 2018) we cannot parametrize a type over integer values, which mean that to differentiate two matrix types we cannot use numbers (even hard-coded in compile time). This is clearly not ergonomic but people are working on it! To circumvent the problem the nalgebra crate created 127 types for you named U1 to U127.

So for example here is what we could write if we want to create a new matrix type representing matrix of integer with 4 rows and 27 columns.

extern crate nalgebra as na;
use na::{Matrix, U4, U27, MatrixArray};

type MyMatrixType = Matrix<i32, U4, U27, MatrixArray<i32, U4, U27>>;

As you can see, the buffer we used is a MatrixArray one. It is an array which means its size is determined at compile time. So what about if you only know the size of your matrix when running your program (like a matrix made from the result of another program) ?

EDIT (December 22, 2018): The author of nalgebra told me that there was an easier way to create your matrix type. Following the previous example we can omit to specify the buffer if we use the *MatrixMN* type.

extern crate nalgebra as na;
use na::{MatrixMN, U4, U27};

type MyMatrixType = MatrixMN<i32, U4, U27>;

Which is far more elegant!

If you do not know anything about the size of your vector/matrix then you can use the alias DVector / DMatrix.

use na::{DVector, DMatrix};

let v = DVector::from_row_slice(5, &[1,2,3,4,5]);
// We need the type annotation to specify the scalar type (here f32)
let m: DMatrix<f32> = DMatrix::identity(205,18);

Note that in this case we cannot simply use the new method as the compiler would have no idea of which dimension is the matrix. We instead need to initialize the matrix with different methods that requires extra parameters to specify the dimensions.

What if you know the row dimension at compile time but not the column dimension (or the opposite) ? In this case you can create your matrix type by using the Dynamic type.

extern crate nalgebra as na;
use na::{Matrix, U49, Dynamic, MatrixVec};

type MyMatrixType = Matrix<f64, U49, Dynamic, MatrixVec<f64, U49,Dynamic>>;

We created a new matrix type with 49 rows and a dynamical number of columns of 64 bits floats. Because the matrix is dynamically allocated we cannot use the MatrixArray buffer, instead we are using the MatrixVec one.

EDIT (December 22, 2018): Again, we do not need to declare the buffer type.

extern crate nalgebra as na;
use na::{MatrixMN, U49, Dynamic};

type MyMatrixType = MatrixMN<f64, U49, Dynamic>;

The different ways to initialize

There are many ways to initialize matrices and we will not cover everything here. Instead we will only see the most common one.

Zeros

The most common initialization would be to have a matrix filled of zeros. It is of course simple to do it in nalgebra but there is a small catch. Because Rust needs to know the type of the element of the matrix at compile time it will try to infer it from your code. However, if you program incrementaly you may first just initialize your matrix and later add the computation that will use it. In this case you will need to give a type annotation to tell the Rust compiler which it is dealing with.

let m1 : Matrix7x5<i32> = Matrix7x5::zeros();

// We need to precise the number of rows and columns in the dynamic case.
let m2 : DMatrix::<u8> = DMatrix::zeros(42, 51);

Identical value

When you want to initialize your matrix with the same element you can simply use the from_element method.

let m1 = Matrix3x5::from_element(1);

// We need to precise the number of rows and columns in the dynamic case.
let m2 = DMatrix::from_element(27, 162, 5.78);

Identity

A lot of time we need to have an identity matrix and it would be bothersome to write all the values by ourselves. You can instead just use the identiy method which will handle it for you. You can even use it for non square matrices but it will set only the largest square submatrix as the identity, the rest will be filled by 0.

let m1 = Matrix4::identity();
/* m1 = (1.0, 0.0, 0.0, 0.0,
         0.0, 1.0, 0.0, 0.0,
         0.0, 0.0, 1.0, 0.0,
         0.0, 0.0, 0.0, 1.0)
*/

let m2 = Matrix3x6::identity();
/* m2 = (1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
         0.0, 1.0, 0.0, 0.0, 0.0, 0.0,
         0.0, 0.0, 1.0, 0.0, 0.0, 0.0) 
*/

// We need to precise the number of rows and columns in the dynamic case.
let m3 = DMatrix::identity(4,2);
/* m3 = (1.0, 0.0,
         0.0, 1.0,
         0.0, 0.0,
         0.0, 0.0)
*/

From an iterator

In Rust iterators are really convenient and easy to use. It is possible that you will get your matrix values from one so nalgebra gives you the from_iterator method. Here is a simple example where the matrices are initialized from two iterators.

let square_iter = (1..9).map(|i| i*i);
let m1 = Matrix2x4::from_iterator(square_iter.clone());
/* m1 = (1,  9, 25, 49,
         4, 16, 36, 64)
*/

let odd_iter = (0..12).filter(|i| i%2==1);

// We need to precise the number of rows and columns in the dynamic case.
let m2 = DMatrix::from_iterator(3, 2, odd_iter);
/* m2 = (1, 7, 
         3, 9,
         5, 11)
*/

Note that the matrix is filled with column major-order.

From a function

Sometime you do not have an iterator but you have directly the function to generate the values of your matrix. In this case it could be easier to just use the from_fn method.

// We need to convert the usize type to floats
let m1 = Matrix2x3::from_fn(|r,c| r as f32 + c as f32);
/* m1 = (0, 1, 2,
        1, 2, 3)
*/
 
fn foo(r: usize, c: usize) -> i64 {
    if r > 1 { return 1 }
    else { return 0 }
}
// We need to precise the number of rows and columns in the dynamic case.
let m2 = DMatrix::from_fn(3, 5, foo);
/* m2 = (0, 0, 0, 0, 0,
         0, 0, 0, 0, 0,
         1, 1, 1, 1, 1)
*/

Accessing and modifying elements

Accessing and modifying values of matrices are fundamental operations for scientific computing. However, it could be a little bit tricky at first to do it because of the Rust compiler. Indeed, the infamous borrow checker will not let you do anything that is not safe.

It is important to be clear on if we only want to access an element to read it/copy it or if we want to modify it. If we just want to read some value of our matrix it is simple.

EDIT (December 22, 2018): The author of the crate reminded me of the tuple indexing. It means that to access the value of a matrix we can either consider the matrix as a single vector and use a single integer or see it as … a matrix and use two integers!

let m1 = Matrix3::new(0, 1, 2,
                      3, 4, 5,
		      6, 7, 8);

// Matrix seen as a vector
let a_21 = m1[1*3+0]; // a_21 = 1
let a_33 = m1[2*3+2]; // a_33 = 8

// Matrix seen as a matrix
let b_21 = m1[(0,1)]; // b_21 = 1
let b_33 = m1[(2,2)]; // b_33 = 8


let m2: DMatrix<i32> = DMatrix::identity(4,5);
/* m2 = (1, 0, 0, 0, 0,
         0, 1, 0, 0, 0,
         0, 0, 1, 0, 0,
         0, 0, 0, 1, 0)
*/

// Matrix seen as a vector
let c_44 = m2[3*4+3];       // c_44 = 1

// Matrix seen as a matrix
let d_44 = m2[(3,3)];       // d_44 = 1

let e = m2[0] + m2[3*4+3];  // e = 2

You have to remember that when considering the matrix as a single vector the indexing is done in column-major mode. On the other hand, when indexing as a matrix the indexes are in row-major mode: m[(i,j)] will be the value of the i-th row and the j-th column of m. Beware, unlike the mathematical notation both indexing method are 0-based.

But now, suppose that you want to have access to the values of a submatrix. It is a pretty common task so nalgebra gives you many functions to do it. Here we will just see how to read a rows, columns and a simple submatrix (you can have a lot of fun by defining complex submatrix). Basically what we will do is to take a slice of the matrix. The slice will also be a Matrix type (with a special buffer type) so you can manipulate it as a normal matrix except that you cannot mutates its values.

let m : Matrix4x5<f32> = Matrix4x5::identity();
/* m = (1, 0, 0, 0, 0,
        0, 1, 0, 0, 0,
        0, 0, 1, 0, 0,
        0, 0, 0, 1, 0)
*/


// Here there is no need to specify the
// resulting matrix type because a row
// and a column have known size.
let r1 = m.row(2); 
// r1 = (0, 0, 1, 0, 0)
  
let c1 = m.column(2);
/* c1 = (0,
         0,
         1,
         0)
*/

// Here we need to specify the dimensions.
// Because the sizes are known at compile
// time there are the prefix "fixed" to the methods.
let r2 = m.fixed_rows::<U3>(1);
/* r2 = (0, 1, 0, 0, 0,
         0, 0, 1, 0, 0,
         0, 0, 0, 1, 0)
*/

let c2 = m.fixed_columns::<U2>(2);
/* c2 = (0, 0,
         0, 0,
         1, 0,
         0, 1)
*/

let s1 = m.fixed_slice::<U3,U4>(1,0);
/* s = (0, 1, 0, 0, 
        0, 0, 1, 0, 
        0, 0, 0, 1)
*/

    
// Here size is supposed to be known at runtime.
// It could be initialized by the user input for example.
let size = 2;

let r3 = m.rows(0, size);
/* r3 = (1, 0, 0, 0, 0,
         0, 1, 0, 0, 0)
*/

let c3 = m.columns(1, size);
/* c3 = (0, 0,
         1, 0,
         0, 1,
         0, 0)
*/
    
// start and shape are known at runtime
let start = (1,1);
let shape = (2,3);
let s2 = m.slice(start, shape);
/* s2 = (1, 0, 0,
         0, 1, 0)
*/

Note that as always, because the slices we are taking are in fact matrices we also need to know if their sizes are determined at compile time or not.

So what about modifying the values of our matrix ? You simply need to tell the compiler that you want your matrix to be mutable. However, remember that a slice cannot be mutated as it is just a view of a part of another matrix.

let mut m1 = Matrix3::identity();
m1[0] = 5;
/* m1 = (5, 0, 0
         0, 1, 0,
         0, 0, 1)
*/

let mut m2 = DMatrix::from_row_slice(3,4,&[1.5, -4.7, 0.2, -5.9,
                                           5.8, 7.81, 0.8, 17.2,
                                           1.4, -2.1, 6.6, -4.1])
for i in 0..12 {m2[i]=i as f32}
/* m2 = (0.0, 3.0, 6.0, 9.0,
         1.0, 4.0, 7.0, 10.0,
         2.0, 5.0, 8.0, 11.0)
*/

Last but not least, what if you want to change the size of your matrix ? Again, nalgebra is making a difference between compile time and runtime. If you know at compile time how many rows or columns you want to add or remove then your matrix will stay statically-sized (except if it is already dynamically sized of course). Conversely, if you do not know at compile time how you want to resize your matrix then the output will be a dynamically sized matrix. Because there are many different methods we will just show a subset of them. You can read all in the documentation.

let m = Matrix4x5::<i32>::identity();
/* m = (1, 0, 0, 0, 0,
        0, 1, 0, 0, 0,
        0, 0, 1, 0, 0,
        0, 0, 0, 1, 0)
*/

let m = m.remove_row(2);
/* m = (1, 0, 0, 0, 0,
        0, 1, 0, 0, 0,
        0, 0, 0, 1, 0)
*/
    
let m = m.insert_column(3,-15);
/* m = (1, 0, 0, -15, 0, 0,
        0, 1, 0, -15, 0, 0,
        0, 0, 0, -15, 1, 0)
*/

// Here n is initialized at runtime
let n = 2;
let m = m.remove_columns(1,n);
/* m = (1, -15, 0, 0,
        0, -15, 0, 0,
        0, -15, 1, 0)
*/

let m = m.insert_rows(2, n, 188);
/* m = (  1, -15,   0,   0,
          0, -15,   0,   0,
        188, 188, 188, 188,
        188, 188, 188, 188,
          0, -15,   1,   0)
*/

Note that during the resizing you need to use the let keyword to mean that you redefine the type of your variable. Indeed, because the size of your matrix is changing its type is also changing.

Matrix operations

Having to specify the types of the matrices we are using can be a pain but it begins to be useful when we are manipulating them. It is simple, if you do a mistake by manipulating two matrices with non-suitable dimensions it will not compile! Unlike for example numpy in Python where the matrix are all dynamically sized and where you discover your mistake by running your program (which can be time expensive) in Rust you will know it immediately! Of course, there is no magic. If your matrices are dynamically sized you will also have to debug at runtime. Therefore, the more information you give to the compiler and the easier it will protect you.

Addition and multiplication

Now that we have seen all of the basics we can begin to use our matrices for real! Additioning and multiplicating matrices are two operations really easy in nalgebra, just consider them as simple scalar values and it will work out of the box!

let a = Matrix2::new(1, 2,
                     3, 4);
let b = Matrix2::new(5, 0,
                     0, 0);
let c = a+b;
/* c = (6, 2,
        3, 4)
*/
    
let d = a*b;
/* d = ( 5, 0,
        15, 0)
*/

// We can use slices for operations
let e = a*b.column(0);
/* e = ( 5,
        15)
*/

Besides, it is worth insisting on the protection against matrix shapes mismatch. The Rust compiler will prevent you to do not compatible operation at compile time but only if you give enough information to it! It is always advantageous to you use methods requiring compile time information because first it will protect you again stupid mistakes and second it will be faster (the compiler could use the stack instead of the heap)!

Linear system

Most of the interest of using nalgebra is to do…algebra! And a very common operation in algebra is to solve a linear system. However, nalgebra gives you many ways to solve it. Indeed, if your linear system is in the form Ax=b then computing A^-1 may not be the fastest solution! It is not in the scope of this article to explain about all the different algebra tricks but the core idea is to exploit matrices’ structures to reduce the computation time.

If you just want the inverse of a function, you can call the try_invert_to function. The function takes two arguments, an input matrix and an output matrix. By performing a LU decomposition on the input matrix the function is trying to find an inverse and store the result in the output matrix. Be careful, if no inverse is found your output matrix may contain invalid data (so check the result).

let mut m = Matrix2::new(1.0, 2.0,
                         3.0, 4.0);
let mut result = Matrix2::zeros();
    
na::linalg::try_invert_to(m, &mut result);
/* result = ( -2,    1,
             1.5, -0.5)
*/

Now what can we do if we want to solve the classical linear system Ax=b ? Depending on the property of your matrix you can use different decompositions to speed up the resolution. I will let you take care of the choice but here I will just show some examples. The following code solves a linear system in three different ways. Once by inversing the matrix A and multiplying it by b (very unefficient). Once by computing the LU decomposition and asking nalgebra to solve it. And finally once by computing the Cholesky decomposition and using nalgebra to solve it (faster than the LU decomposition !).

let mut A = Matrix2::new(1.0, -1.0,
                         -1.0, 4.0);

let mut b = Vector2::new(1.0,5.0);
    
let mut A_inv = Matrix2::zeros();
na::linalg::try_invert_to(A, &mut A_inv);
let A_lu = A.lu();
let A_cho = A.cholesky();

println!("{}",A_inv*b);
println!("{}",A_lu.solve(&b).unwrap());
println!("{}",A_cho.unwrap().solve(&b));

/* result = (3,
             2)
*/

Please note that the in the case of the LU decomposition the result of the solve function is an Option. It means we have to check if a solution has been found. In the case of the Cholesky decomposition it is the decomposition itself that is an Option (because the matrix may not be symmetric definite positive) but the solve function outputs a Matrix.

Other

I cannot write about all the possibilities that nalgebra offers (that is the documentation duty) because this post is already too long. Just for you to know, nalgebra allows you to perform many factorization of your matrices (QR, SVD…), finding the determinant and the eigenvalues, and to do a lot of computer graphic operations (rotation, projection…). As always, you can check it in the documention.