I would like to concatenate anonymous calls in Julia

user_11305741 0
Viewed 104

So, I'm learning more about Julia and I would like to do the following:

I have a 3 row by 2 columns matrix, which is fixed, A = rand(2,3)

julia> A = rand(2,3)
2×3 Matrix{Float64}:
 0.705942  0.553562  0.731246
 0.205833  0.106978  0.131893

Then, I would like to have a anonymous function, which does the following:

a = ones(1,3);
a[2] = rand();

Finally, I would like to broadcast broadcast(+, ones(1,3) => a[2]=rand(), A)

So I have the middle column of A, i.e., A[:,2], added by two different random numbers, and in the rest of the columns, we add ones.

EDIT: If I add a, as it is:

julia> a = ones(1,3)
1×3 Matrix{Float64}:
 1.0  1.0  1.0

julia> a[2] = rand()
0.664824196431979

julia> a
1×3 Matrix{Float64}:
 1.0  0.664824  1.0

I would like that this a were dynamic, and a function.

So that:

broadcast(+, a, A)

Would give:

julia> broadcast(+, a, A)
2×3 Matrix{Float64}:
 1.70594  0.553562 + rand()  (correct)          1.73125
 1.20583  0.106970 + rand()  (different rand()) 1.13189

Instead of:

julia> broadcast(+, a, A)
2×3 Matrix{Float64}:
 1.70594  1.21839  (0.553562 + -> 0.664824)   1.73125
 1.20583  0.771802 (0.106978 + -> 0.664824)   1.13189

So, I thought of this pseudo-code:

broadcast(+, a=ones(1,3) => a[2]=rand(), A)

Formalizing:

broadcast(+, <anonymous-fucntion>, A)

Second EDIT:

Rules/Constrains:

  • Rule 1: the call must be data-transparent. That is, A must not change state, just like when we call f.(A).
  • Rule 2: not creating an auxiliary variable (a must not exist). The only vector that must exist, before and after, the call is A.
  • Rule 3: f.(A) must be anonymous; that is, you can't use define f as function f(A) ... end
julia dataframe anonymous-function
4 Answers

Firstly I'd like to say that the approach taken in the other answers is the most performant one. It seems like you want the entire matrix at the end, in that case for the best performance it is generally good to get data (like randomness) in big batches and to not "hide" data from the compiler (especially type information). A lot of interesting things can be achieved with higher level abstractions but since you say performance is important, let's establish a baseline:

function approach1(A)
    a = ones(2,3)
    @. a[:, 2] = rand()
    broadcast(+, a, A)
end

julia> A = rand(2,3)                         
 2×3 Matrix{Float64}:                         
   0.199619   0.273481  0.99254                
   0.0927839  0.179071  0.188591                                                                  julia> @btime approach1($A)                    
   65.420 ns (2 allocations: 256 bytes)         
 2×3 Matrix{Float64}:                           
  1.19962  0.968391  1.99254                    
  1.09278  1.14451   1.18859

With that out of the way let's try some other solutions.

If a single row with lazy elements doesn't count as an auxiliary variable this seems like a good starting point:

function approach2(A)                           
    a = Matrix{Function}(undef, 1, 3)
    fill!(a, ()->1.0)
    a[2] = rand
    broadcast((a,b)->a() + b, a, A)                 
end

We get a row a = [()->1.0 rand ()->1.0] and evaluate each function when the broadcast gets that element.

julia> @btime approach2($A)
  1.264 μs (24 allocations: 960 bytes)

The performance is 20 times worse, why? We've hidden type information from the compiler, it can't tell that a() is a Float64 by just asserting this (changing the last row to broadcast((a,b)->a()::Float64 + b, a, A) increase the performance almost tenfold:

julia> @btime approach2($A)                   
    164.108 ns (14 allocations: 432 bytes)

If this is acceptable we can make it cleaner: introduce a LazyNumber type that keeps track of the return type, and has promote rules/operators so we can get back to broadcast(+, ...). However, we are still 2-3 times slower, can we do better?

An approach that could allow us to squeeze out some more would be to represent the whole array lazily. Something like a Fill type, a LazySetItem that applies on top of a matrix. Once again actually creating the array will be cheaper unless you can avoid getting parts of the array

user_152812240

With the caveat that I don't know how much you really learn by setting artificial rules like this, some tidier ways are:

julia> A = [ 0.705942  0.553562  0.731246
             0.205833  0.106978  0.131893 ];  # as given

julia> r = 0.664824196431979;  # the one random number

julia> (A' .+ (1, r, 1))'  # no extra vector
2×3 adjoint(::Matrix{Float64}) with eltype Float64:
 1.70594  1.21839   1.73125
 1.20583  0.771802  1.13189

julia> mapslices(row -> row .+ (1, r, 1), A; dims=2)  # one line, but slow
2×3 Matrix{Float64}:
 1.70594  1.21839   1.73125
 1.20583  0.771802  1.13189

julia> B = A .+ 1; @views B[:, 2] .+= (-1 + r); B  # fast, no extra allocations
2×3 Matrix{Float64}:
 1.70594  1.21839   1.73125
 1.20583  0.771802  1.13189

I can't tell from your question whether you want one random number or two different ones. If you want two, then you can do this:

julia> using Random

julia> Random.seed!(1); mapslices(row -> row .+ (1, rand(), 1), A; dims=2)
2×3 Matrix{Float64}:
 1.70594  0.675436  1.73125
 1.20583  0.771383  1.13189

julia> Random.seed!(1); B = A .+ 1; @views B[:, 2] .+= (-1 .+ rand.()); B 
2×3 Matrix{Float64}:
 1.70594  0.675436  1.73125
 1.20583  0.771383  1.13189

Note that (-1 .+ rand.()) isn't making a new array on the right, it's fused by .+= into one loop over a column of B. Note also that B[:,2] .= stuff just writes into B, but B[:, 2] .+= stuff means B[:, 2] .= B[:, 2] .+ stuff and so, without @views, the slice B[:, 2] on the right would allocate a copy.

user_137067350

I agree that it is not very clear what you are trying to achieve, and even if what you want to learn is how to achieve something or how theoretically something works. If all you want is just to add a random vector to a matrix column (and 1 elsewhere), it is as simple as... add a random vector to the desired matrix column and 1 elsewhere:

julia> A = rand(2,3)
2×3 Matrix{Float64}:
 0.94194   0.691855  0.583107
 0.198166  0.740017  0.914133

julia>  A[:,[1,3]] .+= 1
2×2 view(::Matrix{Float64}, :, [1, 3]) with eltype Float64:
 1.94194  1.58311
 1.19817  1.91413

julia> A[:,2] += rand(size(A,1))
2-element Vector{Float64}:
 1.0306116987831297
 0.8757712661515558

julia> A
2×3 Matrix{Float64}:
 1.94194  1.03061   1.58311
 1.19817  0.875771  1.91413
user_15868600

Why not just have a be the same size as A, and then you don't even need broadcasting or any weird tricks:

julia> A = rand(2,3)
2×3 Matrix{Float64}:
 0.564824  0.765611  0.841353
 0.566965  0.322331  0.109889

julia> a = ones(2,3);

julia> a[:, 2] .= [rand() for _ in 1:size(a, 1)] #in every row, make the second column's value a different rand() result

julia> a
2×3 Matrix{Float64}:
 1.0  0.519228   1.0
 1.0  0.0804104  1.0

julia> A + a
2×3 Matrix{Float64}:
 1.56482  1.28484   1.84135
 1.56696  0.402741  1.10989
sundar_81270
Related