Broadcast type mismatch

General Usage
, ,
1

In the following code

import Base: ^

struct cccbbb
    v::Matrix{Float64}
end

^(cb::cccbbb, p::Int) = cb.v^p
Base.Broadcast.broadcastable(cb::cccbbb) = Ref(cb)

n = 5
r = LinRange(0.0, 10.0, 5)
cb = cccbbb(rand(n,n))

tmp = cb.^2
@show typeof(cb.^2), typeof(cb.^2 ./r), typeof(tmp./r)

the result is

(typeof(cb .^ 2), typeof(cb .^ 2 ./ r), typeof(tmp ./ r)) = (Matrix{Float64}, Vector{Matrix{Float64}}, Matrix{Float64})

The type of (cb .^ 2 ./ r) is not as expected. Is this a bug?

2

Yeah, that’s a little confusing, but I think the result is correct. Here’s a slightly simplified version of your example:

import Base: ^

struct A
    m::Matrix{Float64}
end

^(a::A, p::Int) = a.m ^ p

Base.broadcastable(a::A) = Ref(a)

m = [1 2
     3 4]

a = A(m)

julia> m ^ 2
2×2 Array{Int64,2}:
  7  10
 15  22

julia> m .^ 2
2×2 Array{Int64,2}:
 1   4
 9  16

julia> m .^ 2 ./ [1, 2]
2×2 Array{Float64,2}:
 1.0  4.0
 4.5  8.0

julia> a .^ 2 ./ [1, 2]
2-element Array{Array{Float64,2},1}:
 [7.0 10.0; 15.0 22.0]
 [3.5 5.0; 7.5 11.0]

The Base.broadcastable(a::A) = Ref(a) line tells Julia to treat an object of type A as a scalar for broadcasting. So

a .^ 2 ./ [1, 2]

is effectively equivalent to

[(a^2)/1, (a^2)/2]

And they way you’ve defined ^ on your new type means that it is doing a matrix power (m * m), not an element by element power.

Can you elaborate on what output you expected and/or what output you are trying to achieve?

3 Likes
3

As in your example,

julia> b = a.^2
2×2 Matrix{Float64}:
  7.0  10.0
 15.0  22.0

and

julia> b ./ [ 1, 2]
2×2 Matrix{Float64}:
 7.0  10.0
 7.5  11.0

The result is the Matrix{Float64} type as expected. But the

a .^ 2 ./ [1, 2]

type is

2-element Array{Array{Float64,2},1}

which is not expected.

Can the Base.broadcastable(a::A) = Ref(a) just affect the operator ^?

4

What we’re trying to tell you is that b ./ [1, 2] and a .^ 2 ./ [1, 2] are two distinct computations. An important feature of the dot syntax is that vectorized functions are fused into a single loop. In the second case, you’re thinking Julia would first broadcast exponentiation to obtain the matrix, but as @CameronBieganek noted in his reply, the operations are fused into [(a^2)/1, (a^2/2)]. This is the expected behavior, since element-wise exponentiation is different from exponentiation after all (take, matrices, for example). The behavior here is more confusing than usual because you’re wrapping a matrix, but also tell Julia to treat your struct as a scalar.

If you want the same result with the first case, try a ^ 2 ./ [1,2] instead.

3 Likes