What is the mathematical explanation of adding 2d Array and 1d Array?

Question

I can't seem to replicate this numpy arithmetic. I'm using Julia, but want to know the mathematical explanation for this code. It seems to break what I know about Linear Algebra.

X = np.arange(-5, 5, 0.2).reshape(-1, 1)

X.shape  ## (50, 1)

test = np.sum(X**2, 1).reshape(-1, 1) + np.sum(X**2, 1)

test.shape  ## (50, 50)

In Julia, I would write

X = reshape(collect(range(-5, stop=5, length=N)), :, 1);

size(X)  ## (50, 1)

test = sum(X.^2, dims=2) + vec(sum(X.^2, dims=2));

size(test) ## (50, 1)

I'm trying to think how a 50x50 matrix would be the result of adding two vectors? I know numpy uses a lot of broadcasting under the hood, but it doesn't seem clear to me what this is doing.

What is the mathematical notation or Julia equivalent for what numpy is doing here?

Answer 1

Your are doing a lot of stuff that really obscures your point, which, I believe, concerns how to add arrays of different shapes.

Python:

In [21]: x = np.random.rand(5, 1)                                               

In [22]: x.shape                                                                
Out[22]: (5, 1)

In [23]: y = np.random.rand(1, 4)                                               

In [24]: y.shape                                                                
Out[24]: (1, 4)

In [25]: (x + y).shape                                                          
Out[25]: (5, 4)

Julia:

julia> x = rand(5);

julia> y = rand(1, 4);

julia> x + y
ERROR: DimensionMismatch("dimensions must match")

julia> x .+ y
5×4 Array{Float64,2}:
 1.95779  1.31897   1.23345   1.32423 
 1.78126  1.14244   1.05692   1.14771 
 1.08306  0.444243  0.35872   0.449509
 1.69756  1.05874   0.97322   1.06401 
 1.18661  0.547789  0.462265  0.553054

julia> size(x .+ y)
(5, 4)

As you can tell, Python broadcasts arrays by default, while Julia requires that you specifically ask for it, by using the dot operator, ..

It is exactly because it does not make sense to add two arrays of different shapes, that Julia does not broadcast by default. Similarly, with multiplication, * and .* differ:

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

julia> B = [4 5; 6 7]
2×2 Array{Int64,2}:
 4  5
 6  7

julia> A * B
2×2 Array{Int64,2}:
 16  19
 36  43

julia> A .* B
2×2 Array{Int64,2}:
  4  10
 18  28

The ordinary * is matrix multiplication, while the latter is elementwise array multiplication.

Another example:

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

julia> b = [7, 8]
2-element Array{Int64,1}:
 7
 8

julia> A * b
ERROR: DimensionMismatch("matrix A has dimensions (2,3), vector B has length 2")

julia> A .* b
2×3 Array{Int64,2}:
  7  14  21
 32  40  48