Assignment into array of arrays by broadcasting

General Usage
,
1

I want to simulate a process with 3 discrete outcomes. Let’s say they have probabilities .2, .3 and .5 and I want to run nSim = 10 simulations, each of which is a series of nDraw = 5 draws from the discrete distribution:

using Distributions
d = Multinomial(1,[.2;.3;.5])
M = rand(d, 5, 10)

I get an nDraw x nSim matrix, each element of which is a 3-vector with 2 zeros and 1 one in the slot to show which of the 3 outcomes occurred for the particular draw. Here’s an example:

5×10 Matrix{Vector{Int64}}:
 [1, 0, 0]  [0, 0, 1]  [0, 0, 1]  [0, 0, 1]  [0, 1, 0]  [1, 0, 0]  [0, 0, 1]  [0, 1, 0]  [0, 0, 1]  [0, 0, 1]
 [0, 0, 1]  [0, 1, 0]  [0, 0, 1]  [0, 0, 1]  [0, 0, 1]  [0, 0, 1]  [0, 1, 0]  [0, 1, 0]  [0, 0, 1]  [0, 0, 1]
 [0, 1, 0]  [0, 1, 0]  [0, 1, 0]  [0, 0, 1]  [0, 1, 0]  [0, 1, 0]  [0, 0, 1]  [1, 0, 0]  [0, 0, 1]  [0, 1, 0]
 [1, 0, 0]  [0, 0, 1]  [0, 1, 0]  [0, 0, 1]  [0, 1, 0]  [0, 0, 1]  [0, 0, 1]  [0, 0, 1]  [1, 0, 0]  [1, 0, 0]
 [0, 1, 0]  [1, 0, 0]  [0, 0, 1]  [0, 1, 0]  [0, 1, 0]  [0, 1, 0]  [0, 1, 0]  [0, 0, 1]  [0, 0, 1]  [0, 1, 0]

Anyway, all of this is just to say that, because I’m using Distributions.Multinomial and want to create all the sims at once, I’m stuck with a matrix of vectors…

THE QUESTION (finally): Why can’t I say something like a[1,:] .= [0, 1, 0] to assign that 3 vector to all the elements in the first row of M? When I do so I get:

ERROR: DimensionMismatch: array could not be broadcast to match destination
Stacktrace:
 [1] check_broadcast_shape
   @ ./broadcast.jl:540 [inlined]
 [2] check_broadcast_axes
   @ ./broadcast.jl:543 [inlined]
 [3] instantiate
   @ ./broadcast.jl:284 [inlined]
 [4] materialize!
   @ ./broadcast.jl:871 [inlined]
 [5] materialize!(dest::SubArray{Vector{Int64}, 1, Matrix{Vector{Int64}}, Tuple{Int64, Base.Slice{Base.OneTo{Int64}}}, true}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(identity), Tuple{Vector{Int64}}})
   @ Base.Broadcast ./broadcast.jl:868
 [6] top-level scope
   @ REPL[159]:1

This is weird to me. I know I can do the assignment by looping but I like broadcasting and also want to understand what I’m doing wrong here…

EDIT: The moral of the story is “Think before you type” :<
I should have written the assignment as a[1,:] .= Ref([0, 1, 0]).

In honor of the world cup: Screwed by the Refs once again…

2

You could also try to use an SVector from StaticArrays for the three elements in the multinomial distribution, because I think then the draws might be SVectors as well, which would be beneficial in terms of memory layout and allocations.

2 Likes
3

Good idea, thanks.

6

try this:

M[1,:] .= [[0, 1, 0]]
7

yup, that works also. very interesting… what’s the logic here? normally for a matrix M of elements we’d say M[row,:] = new element rather than M[row,:] .= [new element]. why is this different syntax demanded by the language? in hindsight, the same question applies to my belated solution of M[row,:] .= Ref(new element).

8

Nope, that errors:

julia> M = rand(3, 3)
3×3 Matrix{Float64}:
 0.450693  0.230774   0.282827
 0.46242   0.992411   0.114117
 0.697477  0.0835279  0.832331

julia> M[3, :] = 1.0
ERROR: ArgumentError: indexed assignment with a single value to possibly many locations is not supported; perhaps use broadcasting `.=` instead?
Stacktrace:
...

You need to use broadcasting assignment .= to fill a slice of an array with a single element:

julia> M[3, :] .= 1.0
3-element view(::Matrix{Float64}, 3, :) with eltype Float64:
 1.0
 1.0
 1.0

However, when the elements of your array are themselves broadcastable collections, such as vectors in your case, this won’t do what you want—it tries to broadcast the two collections to match each other’s shape, as you observed in your original post.

julia> M = [rand(3) for _ in 1:2, _ in 1:2]
2×2 Matrix{Vector{Float64}}:
 [0.0224723, 0.512887, 0.379622]  [0.078725, 0.355573, 0.403171]
 [0.366589, 0.296196, 0.33523]    [0.314711, 0.156093, 0.214503]

julia> M[2, :] .= [1.0, 0.0, -1.0]
ERROR: DimensionMismatch: array could not be broadcast to match destination
Stacktrace:
...

Wrapping a Ref or a vector [] around the element protects the element itself from participating in broadcasting because the broadcast will apply to the outer wrapper instead, thus accomplishing what you wanted.

julia> M[2, :] .= Ref([1.0, 0.0, -1.0])
2-element view(::Matrix{Vector{Float64}}, 2, :) with eltype Vector{Float64}:
 [1.0, 0.0, -1.0]
 [1.0, 0.0, -1.0]

If you find this confusing, another alternative is to use the fill! function, though in that case you have to remember to use @view, otherwise you’ll be filling a throwaway copy of the row.

julia> fill!(@view(M[2, :]), [1.0, 0.0, -1.0])
2-element view(::Matrix{Vector{Float64}}, 2, :) with eltype Vector{Float64}:
 [1.0, 0.0, -1.0]
 [1.0, 0.0, -1.0]
1 Like
9

Just one comment/caveat for future reference, in case someone is looking for replacing the entries of the vectors instead of the vectors themselves:

In the above version, all replaced entries of M are referring to the same vector (which I think is what the original question was after)

julia> M = [rand(3) for _ in 1:2, _ in 1:2]
2×2 Matrix{Vector{Float64}}:
 [0.288746, 0.731025, 0.580022]  [0.503281, 0.276495, 0.138057]
 [0.117934, 0.800984, 0.197072]  [0.332391, 0.803246, 0.0430873]

julia> M[2, :] .= Ref([1.0, 0.0, -1.0]);

julia> M
2×2 Matrix{Vector{Float64}}:
 [0.288746, 0.731025, 0.580022]  [0.503281, 0.276495, 0.138057]
 [1.0, 0.0, -1.0]                [1.0, 0.0, -1.0]

julia> M[2,1][1] = 100.0;

julia> M
2×2 Matrix{Vector{Float64}}:
 [0.288746, 0.731025, 0.580022]  [0.503281, 0.276495, 0.138057]
 [100.0, 0.0, -1.0]              [100.0, 0.0, -1.0]

If the elements of the vectors should be replaced without allocating a new copy of the vector every time, one (the only?) possibility would be to propagate the broadcast one level lower:

julia> M = [rand(3) for _ in 1:2, _ in 1:2]
2×2 Matrix{Vector{Float64}}:
 [0.582125, 0.759812, 0.985753]  [0.883305, 0.0040317, 0.0707623]
 [0.430306, 0.963027, 0.484389]  [0.0516542, 0.957012, 0.649871]

julia> for (target, source) in zip(M[2, :], Iterators.repeated([1.0, 0.0, -1.0]))
           target .= source
       end

julia> M
2×2 Matrix{Vector{Float64}}:
 [0.582125, 0.759812, 0.985753]  [0.883305, 0.0040317, 0.0707623]
 [1.0, 0.0, -1.0]                [1.0, 0.0, -1.0]

julia> M[2,1][1] = 100.0;

julia> M
2×2 Matrix{Vector{Float64}}:
 [0.582125, 0.759812, 0.985753]  [0.883305, 0.0040317, 0.0707623]
 [100.0, 0.0, -1.0]              [1.0, 0.0, -1.0]
10

Thanks, Daniel. The fact that I wrote M[row,:] = new element rather than M[row,:] .= new element was just a typo on my part :slight_smile: Sorry for the misdirection…

But I still don’t really understand why M[1,:] .= 3 should work if M is a matrix of scalars but M[1,:] .= [3,0] doesn’t work if M is a matrix of 2-vectors. Why must it be M[1,:] .= [[3,0]] ?

More generally, what information is conveyed to the compiler by the syntax
.= [new element] rather than the more obvious .= new element? I understand and can use the recipe now, but I’d like to understand why it’s required - it’ll help me understand Julia better.

11

I think this is an unfortunate artifact from the fact that in Julia numbers themselves are iterable and count as collections containing themselves.

julia> for i in 1; println(i); end
1

julia> 1[1]
1

julia> 1[2]
ERROR: BoundsError
Stacktrace:
 [1] getindex(x::Int64, i::Int64)
   @ Base ./number.jl:98
 [2] top-level scope
   @ REPL[6]:1
12

But wouldn’t the fact that numbers are iterable suggest that M[1,:] .= 3 shouldn’t work and M[1,:].=[3] would be required even in the case of a matrix of scalars? The fact that numbers are iterable is interesting but seems to make the required syntax even less motivated to me :<

13

This is what I was trying to explain, sorry if the punchline drowned in the rather lengthy post. It’s because in the assignment M[1,:] .= [3, 0], both sides of the .= are broadcastable containers, so both are expanded in an attempt to match up their dimensions. Think about what would happen if you had a function f(x, y) and tried to apply it as f.(M[1, :], [3, 0])—once you understand what happens there, you should be able to see what happens here. (Try it out with a function that can take both vectors and scalars, such as f(x, y) = sum(x) + sum(y), and vary the length of the second argument, i.e., try f.(M[1, :], [3]), f.(M[1, :], [3, 0, -3]), and so on.)

When you wrap [3, 0] in another container, such as [[3, 0]] or Ref([3, 0]) or ([3, 0],) it’s the outer container that participates in broadcasting, and since this only has a singleton dimension (or no dimension at all in the case of Ref), it expands to match any other container, matching up its only element with every element in the other container. This is exactly what you need when you want a broadcastable object to be broadcast as a scalar, which is what you’re looking for here.

Numbers are treated as scalars under broadcasting even if they are iterable and indexable. From the documentation:

By default, only some argument types are considered scalars, including (but not limited to) Numbers, Strings, Symbols, Types, Functions and some common singletons like missing and nothing. All other arguments are iterated over or indexed into elementwise.

https://docs.julialang.org/en/v1/manual/arrays/#Broadcasting

14

Based on this description, I would have thought that following code is more appropriate:

julia> d = Multinomial(5, [0.2,0.3,0.5]); # each simulation has 5 draws

julia> rand(d, 10) # 10 simulations
3×10 Matrix{Int64}:
 1  1  1  0  0  1  2  0  2  0
 3  3  2  4  2  1  1  1  1  2
 1  1  2  1  3  3  2  4  2  3

Each column represents the outcome of 5 draws, and there are 10 simulations. The first has 1 of the first outcome, 3 of the second, and 1 of the third, etc.

I don’t know how this works with the rest of your setup, but it avoids the awkward ‘matrix-of-vector’ data, which also seems quite wasteful in terms of representation.

1 Like
15

@DNF: Your way is much cleaner, thanks.

@danielwe: Got it, finally. Thanks for elaborating.

Thanks to the others, too. I’m good now…