Is there a way to swap columns in O(1) in Julia?

I picked up Julia to do some numerical analysis stuff and was trying to implement a full pivot LU decomposition (as in, trying to get an LU decomposition that is as stable as possible). I thought that the best way of doing so was finding the maximum value for each column and then resorting the columns in descending order of their maximum values.

Is there a way of avoiding swapping every element of two columns and instead doing something like changing two references/pointers?

Upvotes: 4

Views: 2358

Answers (4)

Lakshya
Lakshya

Reputation: 11

Why not use something as Simple as this? Is there a computational load that this incurs, which I'm missing?

julia> y = randn(4,4)
4×4 Matrix{Float64}:
  0.382866   1.95684     0.679879  -1.0213
 -0.644027   0.0286492   0.907008  -0.0116672
  1.13972   -1.66938     0.731684   0.39922
 -0.508178  -1.37069    -0.770984   0.058949

julia> temp1 = y[:,1]
4-element Vector{Float64}:
  0.3828664714656943
 -0.6440268000802201
  1.1397177369324478
 -0.5081782475166619

julia> temp2 = y[:,2]
4-element Vector{Float64}:
  1.9568436266456812
  0.0286492205511295
 -1.6693793024656984
 -1.3706938971653277

julia> y[:,1]=temp2
4-element Vector{Float64}:
  1.9568436266456812
  0.0286492205511295
 -1.6693793024656984
 -1.3706938971653277

julia> y[:,2]=temp1
4-element Vector{Float64}:
  0.3828664714656943
 -0.6440268000802201
  1.1397177369324478
 -0.5081782475166619````

Upvotes: 1

carstenbauer
carstenbauer

Reputation: 10147

Just for "completeness", in case you actually want to swap columns in-place,

function swapcols!(X::AbstractMatrix, i::Integer, j::Integer)
    @inbounds for k = 1:size(X,1)
        X[k,i], X[k,j] = X[k,j], X[k,i]
    end
end

is simple and fast.

In fact, in an individual benchmark for small matrices this is even faster than the view approach mentioned in the other answers (views aren't always free):

julia> A = rand(1:10,4,4);

julia> @btime view($A, :, $([3,2,1,4]));
  31.919 ns (3 allocations: 112 bytes)

julia> @btime swapcols!($A, 1,3);
  8.107 ns (0 allocations: 0 bytes)

Upvotes: 5

tholy
tholy

Reputation: 12179

Following up on @longemen3000's answer, you can use views to swap columns. For example:

julia> A = reshape(1:12, 3, 4)
3×4 reshape(::UnitRange{Int64}, 3, 4) with eltype Int64:
 1  4  7  10
 2  5  8  11
 3  6  9  12

julia> V = view(A, :, [3,2,4,1])
3×4 view(reshape(::UnitRange{Int64}, 3, 4), :, [3, 2, 4, 1]) with eltype Int64:
 7  4  10  1
 8  5  11  2
 9  6  12  3

That said, whether this is a good strategy depends on access patterns. If you'll use elements of V once or a few times, this view strategy is a good one. In contrast, if you access elements of V many times, you may be better off making a copy or moving values in-place, since that's a price you pay once whereas here you pay an indirection cost every time you access a value.

Upvotes: 4

longemen3000
longemen3000

Reputation: 1313

in julia there is the @view macro, that allows you to create an array that is just a reference to another array, for example:


A = [1 2;3 4]
Aview = @view A[:,1] #view of the first column
Aview[1,1] = 10

julia> A
2×2 Array{Int64,2}:
 10  2
 3  4

with that said, when working with concrete number types (Float64,Int64,etc), julia uses contiguous blocks of memory with the direct representation of the number type. that is, a julia array of numbers is not an array of pointers were each element of an array is a pointer to a value. if the values of an array can be represented by a concrete binary representation (an array of structs, for example) then an array of pointers is used.

I'm not a computer science expert, but i observed that is better to have your data tightly packed that using a lot of pointers when doing number crunching.

Another different case is Sparse Arrays. the basic julia representation of an sparse array is an array of indices and an array of values. here you can simply swap the indices instead of copying the values

Upvotes: 1

Related Questions