Skip to main content
7 of 9
Commonmark migration

How to implement the general array broadcasting method from NumPy?

A friend of mine introduced array broadcasting in the Python NumPy package which is very convenient (and also highly efficient).

The idea is perfectly shown in this picture:

scheme of the broadcasting method in NumPy

Basically, the method first checks the shape of the two arrays; if a dimension is not the same, it "broadcasts" that dimension to generate arrays of the same dimensions.

Here is an excerpt from the General Broadcasting Rules in the documentation of NumPy:

When operating on two arrays, NumPy compares their shapes element-wise. It starts with the trailing dimensions, and works its way forward. Two dimensions are compatible when

  1. they are equal, or

  2. one of them is 1

If these conditions are not met, a ValueError: frames are not aligned exception is thrown, indicating that the arrays have incompatible shapes. The size of the resulting array is the maximum size along each dimension of the input arrays.

Arrays do not need to have the same number of dimensions.

This is different from the built-in auto-threading in Mathematica. For example, Mathematica does not do this:

{1, 2} + {{1, 2}, {2, 3}, {3, 4}}

I know that there is a duplicate question. But there is no strong reason why Mathematica can't support such a technique. At least, I think that it doesn't cause any contradictions to Mathematica's existing list operation: we just need to check shape first and then "broadcast" it, which seems quite natural. And perhaps broadcasting can yield an efficiency boost because we don't need to transpose twice.

How could this technique be implemented in Mathematica?

edit

I just run a comparison for python and mma regard adding a vector to a matrix. Python's numpy is faster. The matrix is random

data=RandomReal[{0,1},{40000000,2}];

for mma

Transpose[{1.,2.}+Transpose[data];//AbsoluteTiming

Takes 1.8sec

for python

import numpy as np
import time
a=np.random.rand(40000000,2)
b=np.array([1.,2.])
start=time.time()
a+b
end=time.time()
print end-start

takes 1.08sec.

I think for mma, time is wasted in Transpose, because simply Transpose[data] takes 0.6sec

matheorem
  • 17.4k
  • 8
  • 49
  • 118