Passing Matrix slice to function

nalgebra
1

Hi,

Firstly i’m new to rust and come from the python world. There is a good chance my question is unrelated to nalgebra but rather a miss-understanding of the language.

here’s my problem:

I want to iterate over the rows of a matrix and pass the returned slice to a function, however I couldn’t figure out what type annotation to use in the function signature. I ended up just doing what the compiler told me to do… however my understanding is lacking and I’m sure there is a better way.

This is what I came up with:

extern crate nalgebra as na;


use na::{Matrix, U1, U2, U100, ArrayStorage, SliceStorage};
use rand::Rng;

fn real_fn(X: Matrix<f64, U1, U2, SliceStorage<f64, U1, U2, U1, U100>>)
           -> na::Matrix<f64, na::U1, na::U1, na::ArrayStorage<f64, na::U1, na::U1>> {
    2.0 * X.column(0) - 3.4 * X.column(1)
}


fn main() {
    type Matrix2x100 = Matrix<f64, U100, U2, ArrayStorage<f64, U100, U2, >>;

    let dist = rand::distributions::Normal::new(20.0, 1.5);
    let X = Matrix2x100::from_iterator(rand::thread_rng().sample_iter(&dist
    ));


    for i in 0..100 {
        real_fn(X.row(i));
    }
}
2

Hi!

Here are various cleaner version of what you are trying to achieve here. The code bellow is commented but let me know if you need more explanations:


use na::{U1, U2, U100, MatrixSlice1x2, MatrixMN, RowVector2, Dim};
use rand::Rng;

// The `.row` method will give you a matrix slice with one row and two columns.
// Such a slice has type `MatrixSlice1x2` (which is an alias to what you did already).
// The two dimension arguments `U1` and `U100` are the row-stride and column-stride of the slice. The
// column stride must be `U100` here because your original matrix has 100 rows.
fn real_fn(X: MatrixSlice1x2<f64, U1, U100>) -> f64 {
    // Note that you can access the vector component using indexing directly.
    2.0 * X[0] - 3.4 * X[1]
}

// A variant where we pass an owned vector (i.e. without reference to the original matrix) instead of a slice.
// You can achieve this by calling `X.row(i).into()` which will do the conversion.
fn real_fn_owned(X: RowVector2<f64>) -> f64 {
    // Note that you because this is an owned vector, you can access its components directly
    // instead of using indexing.
    2.0 * X.x - 3.4 * X.y
}


// This is similar to the `real_fn` method, except that it is generic wrt
// the column-stride of the original matrix. This will allow you to reuse
// this function with an original matrix of arbtrary dimension.
fn real_fn_generic<CStride: Dim>(X: MatrixSlice1x2<f64, U1, CStride>) -> f64 {
    // Note that you can access the vector component using indexing directly.
    2.0 * X[0] - 3.4 * X[1]
}


fn main() {
    // The type alias should be called `Matrix100x2` instead of
    // `Matrix2x100` because the usual mathematical convention is
    // to put the number of rows first.
    //
    // We use the `MatrixMN` alias here. It will fill the storage type for you.
    type Matrix100x2 = MatrixMN<f64, U100, U2>;

    let dist = rand::distributions::Normal::new(20.0, 1.5);
    let X = Matrix100x2::from_iterator(rand::thread_rng().sample_iter(&dist));


    for i in 0..100 {
        real_fn(X.row(i));
        real_fn_owned(X.row(i).into());
        real_fn_generic(X.row(i));
    }
}
3 Likes
3

Thank you so much, that is very helpful. Much appreciated

4

Hi, sorry for necrobumping this post, but I’m a little confused about RStride and CStride of SliceStorages, why in MatrixSlice1x2<f64, U1, U100> the CStride is U100, I’m not understanding what it means in depth, what do you actually mean with stride, why the column stride is equal to the number of rows and why the row stride is equal to U1?

5

@lans98 Hi!

nalgebra represents matrices as arrays of elements arranged in column-major orders. For example the matrix:

let matrix = Matrix3x4::new(
   1.0, 4.0, 7.0, 10.0,
   2.0, 5.0, 8.0, 11.0,
   3.0, 6.0, 9.0, 12.0
)

is represented as an array equal to:

[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0]

As you can see all the columns are stashed one after the other in this array. Note that:

  • In this array, each element of the same column are contiguous. In other words, the component matrix[(i + 1, j)] is right after the component matrix[(i, j)] in the array. This index increment between two consecutive elements of the same matrix columns is called the row-stride (or RStride in nalgebra). It is equal to 1 here.
  • In this array, each element of the same row are not contiguous. The component matrix[(i, j + 1)] is 3 elements after matrix[(i, j)] in the array. This index increment between two consecutive elements of the same matrix row is called the column-stride (or CStride in nalgebra). It is equal to 3 here.

So for owned column-major matrices, the row-stride is 1 and column stride is equal to the number of rows (which is 3 in this example). For matrix slices, those strides may be different if one want to take a sub-matrix. For example a 2x2 slice with the same component array as the previous example but with row-stride equal to 2 and column-stride equal to 9 could be the matrix equal to:

    1.0, 10.0,
    3.0, 12.0,

Let me known if this is not clear enough!

3 Likes
6

Got it! It would be nice to add this explanation on the web page reference, while I was searching on the Internet for this concept I found it in the numpy documentation, it’s basically what you explained here, now I totally got it ! Thank you