Constructing a permutation matrix

General Usage
1

Approaching the problem from a slightly different angle, one can first transform a permutation vector p into a permutation matrix, for instance:

julia> function MP(p)
           ℓ = length(p)
           vcat([ begin t=reshape(repeat([0], ℓ), 1, ℓ); t[i]=1; t; end for i in p ]...)
       end
MP (generic function with 1 method)

I’m pretty sure the implementation can be improved / this is what first worked for me.

MP assembles the corresponding (left multiplicative) permutation matrix. Therefore:

julia> MP([2, 1, 4, 3])
4×4 Matrix{Int64}:
 0  1  0  0
 1  0  0  0
 0  0  0  1
 0  0  1  0

julia> A = [1 1 1 1; 2 2 2 2; 3 3 3 3; 4 4 4 4]
4×4 Matrix{Int64}:
 1  1  1  1
 2  2  2  2
 3  3  3  3
 4  4  4  4

julia> MP([2, 1, 4, 3]) * A
4×4 Matrix{Int64}:
 2  2  2  2
 1  1  1  1
 4  4  4  4
 3  3  3  3
2

Using an explicit multiplication PA of dense matrices to permute the rows of A transforms the cost from O(n^2) to O(n^3) for an n \times n matrix.

But if you do want a permutation matrix, it’s much easier (and faster) to simply do

using LinearAlgebra
MP(p) = I(length(p))[p,:]
5 Likes
3

You’re spot-on on the algorithm complexity increase in general (for matrices without assumptions). I wonder whether there is or could there be a sort of “one-hot-sparse” matrix type (i.e., line/column permutations of identity matrices) with specialized operations that would bring down the complexity of matrix multiplications, since there would be logical ones and zeroes and trivial operations…

4

Sure. StaticPermutations.jl is relevant, however it’s in the type domain (compile-time computation), which might or might not be what you want. In Base there is already PermutedDimsArray, which I guess can be thought of as a lazy multiplication with a permutation matrix. There’s also Permutations.jl.

1 Like
5

In a sense the “sparse” representation of your permutation matrix is just the permutation vector. So you could just wrap that into a type that overloads multiplication to instead perform the permutation of rows/columns.

4 Likes
6

You could, of course, use SparseArrays.jl for this.

But representing the linear operator implicitly by the permutation vector is even more efficient. You could wrap it in a matrix-like type with LinearMaps.jl or with a custom AbstractArray type if you wanted to, though you’d need some compelling motivation to go to the trouble.

1 Like
7

The trouble is relatively minimal, at least to get the basics working:

struct PermutationMatrix <: AbstractMatrix{Bool}
    p::Vector{Int}
end

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

Base.:*(P::PermutationMatrix, A::AbstractMatrix) = A[P.p, :]
Base.:*(A::AbstractMatrix, P::PermutationMatrix) = A[:, P.p]

julia> P = PermutationMatrix([2, 1, 4, 3])
4×4 PermutationMatrix:
 0  1  0  0
 1  0  0  0
 0  0  0  1
 0  0  1  0

julia> A = Matrix(reshape(1:16, 4, 4))
4×4 Matrix{Int64}:
 1  5   9  13
 2  6  10  14
 3  7  11  15
 4  8  12  16

julia> P * A
4×4 Matrix{Int64}:
 2  6  10  14
 1  5   9  13
 4  8  12  16
 3  7  11  15

julia> A * P
4×4 Matrix{Int64}:
 5  1  13   9
 6  2  14  10
 7  3  15  11
7 Likes
8

Although the permutation should be inverse-applied (i.e., with its transpose) on one side or the other if it’s actually to be treated as a permutation matrix, which would make this a bit more complicated.

2 Likes
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
10

This package:

is neglected, but has the same ideas (matrix form, cycle counts, sparse form). Perhaps it can be alived.

3 Likes
11

From this article, it seems to me that the proposed implementation of @gdalle is already correct, isn’t it?

12

Great implementation! Concise, sweet, and readable!!!

13

Nope. It implements the transposed getindex operation, so it is equivalent to P^T. Even if you fix that so that it is the right matrix, it implements A * P as A[:, P.p], which is then equivalent to AP^T. (Even if you reinterpret what you mean by the “permutation vector” p to make his P correct, his A * P and P * A were inconsistent.) You can’t tell for your example P, since that matrix happened to be symmetric.

That’s why I gave the corrected implementation in my post.

2 Likes
14

Also, this seems like a good problem to be lazy about, by defining a PermutedMatrix in addition to PermutationMatrix.
Suppose,

P = PermutationMatrix([1,3,4,2])
A = rand(4,4)
B = rand(4,4)
M = P*A // would be a PermutedMatrix
R = B*M // == B*(P*A) == (B*P)*A
        // and performing permutation on columns
        // is easier and faster.
        // `getindex` on PermutedMatrix would go
        // through the permutation
15

Hard to say without more context. Storing the permuted matrix lazily will make many operations on it much slower.

16

My bad, I wrote it in a hurry based on C_\pi from the Wikipedia article but I did screw up one of the multiplications, thanks for fixing the issues.
Ultimately it also depends how you define the permutation matrix (row or column wise).