Is there a way to write to a whole row / column of a nalgebra Matrix?

Question

I am using the DMatrix struct to allocate dynamically sized matrices, where I repeatedly overwrite each column by the normalized column vector, using the L2-norm.

// a is some DMatrix of arbitrary size
let col_0 = a.column(0);
let norm_of_col_0 = col_0.normalize();

Instead of looping over each cell within the current column:

let row = a.shape().0;
let col = a.shape().1;
for col in 0..ncols {
    let norm_of_col = a.column(col).normalize();
    for row in 0..nrows {
        *a.index_mut((row, col)) = norm_of_col()[row];
    }
}

I would like to directly overwrite the column by its normalized version. The code should semantically look something like this:

*a.index_mut((_, col)) = norm_of_col();

where (_, col) means, that I select column col and _ means the whole row.

More generally speaking, is there a way to overwrite a row or column by a new row or column of the same size and data type? Methods as insert_columns only add columns to an existing matrix.

If so, it is computationally faster to do so, or should I just write a helper method which loops over each cell to update the matrix?

Answer 1

You can do it this way with nalgebra 0.18.0:

use nalgebra::DMatrix;

fn main() {
    let mut m = DMatrix::from_vec(2, 3, (0 .. 6).map(|n| n as f64).collect());
    dbg!(&m);
    for mut col in m.column_iter_mut() {
        let normalized = col.normalize();
        col.copy_from(&normalized);
    }
    dbg!(&m);
}

I haven't measured the performance of this code compared to yours.

Note that copy_from goes over the items without checking bounds at each step, doing the checks only once before the loop instead. I haven't checked if the optimizer can perform an equivalent transformation in your code. This simple benchmark gives an edge to the solution in this answer on my machine (not sure how representative it is; usual benchmark disclaimer applies):

use criterion::{black_box, criterion_group, criterion_main, Benchmark, Criterion};
use nalgebra::DMatrix;

fn normalize_lib(m: &mut DMatrix<f64>) {
    for mut col in m.column_iter_mut() {
        let normalized = col.normalize();
        col.copy_from(&normalized);
    }
}

fn normalize_hand_rolled(a: &mut DMatrix<f64>) {
    let nrows = a.shape().0;
    let ncols = a.shape().1;
    for col in 0..ncols {
        let norm_of_col = a.column(col).normalize();
        for row in 0..nrows {
            *a.index_mut((row, col)) = norm_of_col[row];
        }
    }
}

fn benchmark(c: &mut Criterion) {
    let mut m0 = DMatrix::new_random(100, 100);
    let mut m1 = m0.clone();
    let bench = Benchmark::new("lib", move |b| b.iter(|| normalize_lib(black_box(&mut m0))))
        .with_function("hand_rolled", move |b| {
            b.iter(|| normalize_hand_rolled(black_box(&mut m1)))
        });
    c.bench("normalize", bench);
}

criterion_group!(benches, benchmark);
criterion_main!(benches);

normalize/lib           time:   [26.102 us 26.245 us 26.443 us]
normalize/hand_rolled   time:   [37.013 us 37.057 us 37.106 us]