Question about generic programming

nalgebra
1

Hello,

I’ve implemented matrix centering and scaling (as in the R function scale) as follow:

    let column_means = mat.column_mean();

    let centered_mat = mat.map_with_location(|row, col, n| {
        n - column_means.get(row).unwrap()
    });

    let variance = mat.map(|n| n * n).sum() / mat.len() as f64;
    let stddev = variance.sqrt();
    let scaled_mat = centered_mat.map(|n| n / stddev);

Now I would like to implement both center(matrix) -> matrix and a scale(matrix) -> matrix generic over any matrix.

So far I was not able to find which type my functions should be generic over, any help on how to define those would be really appreciated!

Thanks!

2

Hi!

Here is what this function can look like:

use nalgebra::{Matrix, MatrixMN, RealField, Dim, DefaultAllocator};
use nalgebra::allocator::Allocator;
use nalgebra::storage::Storage;


fn center_scale<N, R, C, S>(mat: &Matrix<N, R, C, S>) -> MatrixMN<N, R, C>
where N: RealField,
      R: Dim,
      C: Dim,
      S: Storage<N, R, C>,
      DefaultAllocator: Allocator<N, R, C> + Allocator<N, R> {
    let column_means = mat.column_mean();

    let centered_mat = mat.map_with_location(|row, col, n| {
        n - column_means[row]
    });

    let variance = mat.map(|n| n * n).sum() / nalgebra::convert(mat.len() as f64);
    let stddev = variance.sqrt();
    centered_mat.map(|n| n / stddev)
}
  • The type parameter N means that the matrix is generic wrt. the scalar type. Here it can be any float type (f32 or f64).
  • The type parameter R means that the matrix is generic wrt. the number of rows.
  • The type parameter C means that the matrix is generic wrt. the number of columns.
  • The type parameter S means that the matrix is generic wrt. its internal storage type. This allows this function to be used on both matrixces that own their data as well as matrix slices.
  • The return type is MatrixMN<N, R, C> does not include a parameter for its storage type. This is because its storage type must be deduced automatically by the DefaultAllocator.
  • The DefaultAllocator: Allocator<N, R, C> is here to deduce the storage type of the result matrix.
  • The DefaultAllocator: Allocator<N, R> is required by mat.column_mean().
  • I replaced column_means.get(row).unwrap() by direct indexing column_means[row] because you unwrap it anyways.
  • I replaced mat.len() as f64 by nalgebra::convert(mat.len() as f64) which will convert usize -> f64 -> N.

For even more genericity, you can do this:

use nalgebra::{Matrix, MatrixMN, SimdComplexField, Dim, DefaultAllocator};
use nalgebra::allocator::Allocator;
use nalgebra::storage::Storage;


fn center_scale<N, R, C, S>(mat: &Matrix<N, R, C, S>) -> MatrixMN<N, R, C>
where N: SimdComplexField,
      R: Dim,
      C: Dim,
      S: Storage<N, R, C>,
      DefaultAllocator: Allocator<N, R, C> + Allocator<N, R> {
    let column_means = mat.column_mean();

    let centered_mat = mat.map_with_location(|row, _, n| {
        n - column_means[row]
    });

    let zero: N::SimdRealField = nalgebra::zero();
    let variance: N::SimdRealField = mat.fold(zero, |acc, n| acc + n.simd_modulus_squared()) / nalgebra::convert(mat.len() as f64);
    let stddev = variance.simd_sqrt();
    centered_mat.unscale(stddev)
}

This will allow your function to work on matrices containing complex numbers too, as well as SIMD types (useful if you want to center four matrices simultaneously for example). Note that I replaced your mat.map(...).sum() by a mat.fold. You can do this on the first piece of code I mentioned too and it will avoid the allocation made by the map.

2 Likes
3

Thank you for the quick and outstanding quality response!

This is very instructive and exactly what I was looking for, maybe this kind of tutorial could be documented so it can be easier to find.

What I failed to grasp was that I could define the return type without specifying the allocator explicitly.