Swap cols/rows of a matrix

General Usage
1

Hi,

I need to swap columns or rows of a matrix. I have two ideas:
(1). Create a new matrix and copy cols/rows to the new matrix. We may use @view.
(2). Multiply by a permutation matrix. Of course, we also need to create that permutation matrix (but it’s sparse)

To my understanding, multiplication takes O(n^3) given a square matrix, but not sure how multiplying with a sparse matrix will make a difference or not.
I’m wondering which way is more preferable in terms of performance?

Your help is appreciated in advance!

2

I think it depends also on what you have to do with yr matrix later… With @view you have no allocations but you also have a structure where the data is no longer contigue, so if you have to use it “often” is penalising…

4

I am not sure if I understand what you want, but for swapping cols, for example, I would do this, which does not allocate anything:

julia> function swapcol!(x,i,j)
         for k in axes(x,1) # edited according to next answer
           idata = x[k,i]
           x[k,i] = x[k,j]
           x[k,j] = idata
         end
       end
swapcol! (generic function with 1 method)

julia> x = [ 1 2 3 ; 4 5 6 ]
2×3 Array{Int64,2}:
 1  2  3
 4  5  6

julia> swapcol!(x,1,2)

julia> x
2×3 Array{Int64,2}:
 2  1  3
 5  4  6

julia> @btime swapcol!($x,1,2)
  4.634 ns (0 allocations: 0 bytes)
5

agreed, thanks!

6

Interesting. Why don’t we need to allocate memory for this local temp variable "idate’'?

7

You don’t need the temporary variable:

function _swapcol!(x,i,j)
    for k in axes(x, 1)  # <- give dimension as input to axes function
        x[k, i], x[k, j] = x[k, j], x[k, i]
    end
end

If you are reshuffling the columns (or rows) you can supply a vector of indices:

julia> x = collect(reshape(1:12, 3, 4))
3×4 Array{Int64,2}:
 1  4  7  10
 2  5  8  11
 3  6  9  12

julia> x[:, [2,3,1,4]]
3×4 Array{Int64,2}:
 4  7  1  10
 5  8  2  11
 6  9  3  12
8

idata is just a scalar value, so doesn’t need any memory allocation, unlike an array. But as I said, you don’t really need to create that variable explicitly (though, of course, x, y = y, x will do that ‘under the hood’.)

1 Like
9

This is a good question which I do not completely understand. Because the scope of that scalar is only local to the loop, I guess it is only created in a local memory during computation that does not count as an allocation. I do not know the details of that.

As DNF pointed out, in his version in which this scalar is not created explicitly, it is created anyway. Effectively the @code_typed of the two versions are identical and, if I understand what is there (not sure, though), this part shows that %16 creates a copy of element i, which is then used in the third line of instruction %17.

│    %16 = Base.arrayref(false, x, %14, i)::Int64
│    %17 = Base.arrayref(false, x, %14, j)::Int64
│          Base.arrayset(false, x, %17, %14, i)::Array{Int64,2}
│          Base.arrayset(false, x, %16, %14, j)::Array{Int64,2}

ps: I think that the counter of the loop must be created in a similar way and does not count as an allocation either.

1 Like
10

Sorry for reviving this old thread, but is there something builtin for this now?

11

There’s nothing magic about built-in functions. The suggestions provided above above are just a couple of lines and are basically as good as anything.

There is the internal and undocumented LinearAlgebra.rcswap! that swaps a pair of rows and columns. But I would discourage using that because it is internal (it could change or be removed at any time). Also, it doesn’t perform the originally-requested task of swapping only rows or only columns.

12

There is also Base.swapcols!, but again this is an internal function that may disappear in the future. But since it is only a 10-line function that doesn’t depend on any other Base internals it is perfectly reasonable to copy-and-paste it into your code as needed:

# swap columns i and j of a, in-place
function swapcols!(a::AbstractMatrix, i, j)
    i == j && return
    cols = axes(a,2)
    @boundscheck i in cols || throw(BoundsError(a, (:,i)))
    @boundscheck j in cols || throw(BoundsError(a, (:,j)))
    for k in axes(a,1)
        @inbounds a[k,i],a[k,j] = a[k,j],a[k,i]
    end
end

It’s very similar to the suggestions above, except it adds an @inbounds annotation and consequently requires bounds checks before the loop. You could also do e.g.

if isdefined(Base, :swapcols!)
    using Base: swapcols! # assumes the API hasn't changed
else
    function swapcols!(...)
        ...
    end
end

Inspecting Base code can be a useful exercise, even if you end up copy-pasting it for safety, since it usually has been reviewed by multiple people and tries extra hard to be robust and reasonably efficient.

5 Likes
13 
swap!(m,c1,c2)= begin m[:,c1],m[:,c2]= m[:,c2],m[:,c1];nothing end
14

This allocates two temporary vectors, though.

1 Like
15

Just for fun:

swap2cols(M,i,j) = M[:, replace(axes(M,2), i=>j, j=>i)]
2 Likes
16

Same idea, but for swapping in place (not tested):

swap2cols!(M,i,j) = Base.permutecols!!(M, replace(axes(M,2), i=>j, j=>i))
17

Thanks, this is what I do generally, it’s just one of those cases where it feels a bit clunky to have some unrelated utils.jl file in my package to have this sort of simple-but-unrelated-to-the-actual package functionality.

Maybe a bit like the discussion on whether mean should be available when people could just write sum(x)/length(x) - I get these are all language design decisions which can be reasoned about endlessly and I don’t intend to do that here, just wanted to check in.