Use PermutationMatrix instead of indices

[akleeman]

The sparse Gaussian process implementations used to store a vector of indices which represented a permutation matrix. The reason it did this was that the PermutationMatrix isn't directly serializable, so you end up needing to copy the indices when you deserialize one of these models. Storing the matrices this way led to a number of places with homegrown permutation matrix application, you'd see things like:

  for (Eigen::Index i = 0; i < permutation_indices.size(); ++i) {
    lhs.row(i) = rhs.row(permutation_indices.coeff(i));
  }

These loops are not easy to interpret. In this case it's computing lhs = rhs * P.transpose(), but the transpose is not evident from the code. It's also really easy to introduce a bug in these situations (if the rhs isn't the same shape as how lhs was allocated you can encounter some really tricky memory handling errors).

This PR adds a serialization method for a raw PermutationMatrix<> and changes that to be the type stored in sparse GP fits. While perhaps slightly less efficient in the few situations where we now need to copy indices, that should only really matter when the size of the indices is huge, in which case these copy operations will probably be dwarfed by large linear algebra operations.

[peddie]

This is a great cleanup, thanks.