Constructing a permutation matrix

General Usage
9

Yes, the getindex should be

Base.getindex(P::PermutationMatrix, i::Integer, j::Integer) = P.p[i] == j

and the right-multiplication should be:

function Base.:*(A::AbstractMatrix, P::PermutationMatrix)
    AP = similar(A)
    AP[:, P.p] = A   # computes AP = (PᵀAᵀ)ᵀ = (P⁻¹Aᵀ)ᵀ
    return AP
end

(The test permutation P = PermutationMatrix([2, 1, 4, 3]) happens to be symmetric, which hides these bugs. Try e.g. PermutationMatrix([2, 4, 3, 1]) instead.)

If you wanted to be fancy, you could define all sorts of specialized routines, most obviously for inv(P) and P \ x (ldiv!), but also for things like eigenvectors and eigenvalues (which can be computed from the cycles of the permutation, and moreover you can then multiply by the eigenvectors using FFTs, giving you fast algorithms for things like matrix exponentials and logarithms of permutations etcetera), and more.

It’s a fair amount of work to get something truly full-featured, but it might be fun for someone. Still I’m not sure of the practical motivation.

3 Likes