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...)