Sort with indices

Some linear algebra functions require the output to be sorted. For example, svd returns (u, s, v), where u and v sometimes need to be sorted by s. nalgebra currently doesn't offer this feature so I coded it myself:

let svd = matrix.svd(true, true);

// Add an index to all values and sort the tuples, so we know how it was sorted
let mut s: Vec<(f64, usize)> = svd.singular_values.into_iter().enumerate().map(
    |(idx, &v)| (v, idx)
).collect();
s.sort_unstable_by(|a, b| b.partial_cmp(a).unwrap());
let order: Vec<_> = s.iter().map(|t| t.1).collect();

let s: Array1<_> = s.into_iter().map(|t| t.0).collect();
// also get u and v using order

I feel this is allocating too much:

  1. the zipped array s
  2. s to order
  3. s to singular values again.

In a "better" world, we would only be allocating order. Do you think there's a better way?

1 Like

Use sort instead of sort_unstable so that you won't need to add the indices yourself.

There's unzip to make two vectors in one go.

Sorry, I don't understand your answer. Isn't sort the same as sort_unstable, except for equal elements? I don't see how sort can return the indices. Can you please be clearer?

If you want maximum efficiency, then I think the most straightforward approach will be to fork Rust quicksort implementation and change its signature to:

fn quicksort2<T, T2, F>(v: &mut [T], v2: &mut [T2] mut is_less: F)
    where F: FnMut(&T, &T) -> bool
{
    assert_eq!(v.len(), v2.len());
    // ..
}

And for all v swaps in the code add the same swaps for v2.

So, partly in response to this, I started typing up an obscenely large post on my favorite tool for solving problems like these (tl;dr permutations are really cool, they have incredibly useful properties, and I think they should be in every programmer's problem-solving toolkit). Eventually it reached the point where I started drafting it into a series of blog posts.

Unfortunately, as it grew, it occurred to me that I cannot possibly fit the in the bit that actually responds to your primary concern in my first post. So I'll sever the umbilical cord and just cut to the chase here:


  • You just performed an SVD, a horrendously expensive linear algebra operation.
  • You're worried about a fixed number of allocations. Like, half a dozen.
  • Unless your matrices are something ridiculously tiny (like 3x4), stop worrying!!!
    • And if your matrices ARE ridiculously tiny, profile before worrying!!!

The only way to completely avoid allocations would be to do @newpavlov's idea on sort_unstable (because sort allocates) and write an implementation that swaps the rows of v and the columns of u alongside the elements of s. But repeatedly swapping entire rows and columns sounds like an extremely bad idea and I can't imagine why anyone would want to do this.

Actually, forget sort_unstable. If your matrices are truly tiny enough for this to matter, write a bubble sort! It'll probably be faster.


For comparison, here's the number of allocations I would incur using my own "beloved" toolkit: (you know, the one I wasted so much time not-finishing-a-blog-post about)

// I'm just going to pretend that Vec<_>-backed nalgebra matrices have
//
//     fn rows(&self) -> usize;
//     fn cols(&self) -> usize;
//     fn into_column_major_vec(self) -> Vec<T>;
//     fn from_column_major_vec(usize, usize, Vec<T>) -> Self;
//
// because I still can't find anything by searching their docs.

// NOTE: NOT TESTED

// helper to construct a perm that permutes the first perm.len() elements
// by perm, and leaves the rest of the elements alone
let extend_perm = |mut perm, n| -> Perm {
    assert!(perm.len() <= n); // can't truncate a perm!
    perm.append_mut(&Perm::eye(n - perm.len()))
};

// independently permute the rows and columns of an nalgebra matrix
let permute_nalgebra = |row_perm, col_perm, m| {
    // complete perm is an outer product
    let flat_perm = &col_perm.with_inner(&row_perm);
    
    let (rows, cols) = (m.rows(), m.cols());
    let data = m.into_column_major_vec();
    let data = data.permuted_by(&full_perm);
    MatrixMN::from_column_major_vec(rows, cols, data)
};

// permutation of the singular values
let sv_perm = Perm::argsort(&s);
let s = s.permuted_by(&sv_perm); // sorts s

// rearrange u and v to match
let col_perm = extend_perm(sv_perm.clone(), u.cols());
let row_perm = Perm::eye(u.rows());
let u = permute_nalgebra(row_perm, col_perm, u);

let row_perm = extend_perm(sv_perm.clone(), v.rows());
let col_perm = Perm::eye(v.cols());
let v = permute_nalgebra(row_perm, col_perm, v);

Let's rattle it off:

  • One unavoidable allocation of a Vec<usize> in argsort.
  • Three unavoidable allocations for the output of permute operations. (one for u, one for s, one for v)
  • One pointless and entirely avoidable allocation of an identity perm of length abs(n - m) in one of the calls to extend_perm
    • Avoidable by adding an append_eye(&mut self, n: usize) method to Perm
  • Two avoidable allocations of identity perms of length n and m (u's row_perm and v's col_perm)
    • Avoidable by adding special methods for doing outer products with identity perms
  • Two technically avoidable allocations in the construction of outer products
    • I'd have to hand-write functions for permuting just the rows or colums.

Even the four "unavoidable" allocations listed above can be theoretically "armortized" by keeping around extra vectors, double-buffer style. So look how wasteful that seems; nine avoidable allocations, and that's just how I'd write it.

But in all my time using it, I've never seen Perm anywhere on the top page of a perf report. I've got much bigger fish to fry.


(wait, nine allocations and they're all avoidable? So I guess I contradicted myself on the whole "The only way to completely avoid allocations" bit...)

3 Likes

Thank you all for your answers. I thought that I was missing something obvious but it's now clear that there's no easy way out of this. It's not the first time I want to do this, in Rust and other langages, it would be nice if there was an acceptable solution, like the one proposed by newpavlov.

I really like newpavlov idea but I can't really send a PR to nalgebra with a copy pasted quicksort algorithm :slight_smile: Also, as ExpHP wrote, SVD is super expensive so worrying about 3-4 copies is probably useless.

Well, you could make a crate with copy-pasted code and use it in your PR. :wink: BTW you could even generalize function with the following signature or something similar:

fn quicksort2<T, F>(v: &mut [T], mut is_less: F, swap_callback: F2)
    where F: FnMut(&T, &T) -> bool, F2: FnMut(i: usize, j: usize)
{ .. }