# 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: 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 Mathematica regarding adding a vector to a matrix. Python's Numpy is faster. The matrix is random:

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


For Mathematica:

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


Takes 1.8 sec

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.08 sec.

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

• Is there some reason not to just use Map? That has the advantage of avoiding edge-case issues (what to do if multiple dimensions have equal lengths, etc.) Nov 11 '15 at 16:56
• From the description of the broadcasting process in NumPy (emphasis mine): "Broadcasting provides a means of vectorizing array operations so that looping occurs in C instead of Python". Mathematica already has a means of doing that, i.e. Map and its siblings, as @Daniel mentioned. I would also be somewhat wary of "implicit mapping" operations: I'd rather have a system that pointed out a mismatch in array dimensions and let me deal with it in whatever way is correct for my application. Nov 11 '15 at 17:20
• Outer[Plus, list1,list2] will handle case 3 just fine. Nov 12 '15 at 6:20
• I've changed the question a bit and voted to reopen it. numpys broadcasting have well defined general rules which are not met by current answers to this or linked questions. In my opinion it would be useful, for people coming from python+numpy, to know how to exactly duplicate broadcasting mechanism. Nov 12 '15 at 13:02
• @jkuczm Granted, Python has nice rules, as I've learned today, but Mathematica has its own specifics. It only really makes sense to start not with the trailing dimensions, but with the leading ones in MMA, as the contents of a given array may easily be lists or other nested structures. This is especially so when doing operations on symbolic constructs which may later on turn out to have more complicated nested structures as their -1 level elements. One needs to make sure first, that arrays worked on are proper "rectangular" matrices, not some ragged constructs. Nov 12 '15 at 23:00

NumPy broadcasting lets you perform, in efficient way, element-wise operations on arrays, as long as dimensions of those arrays are considered "compatible" in some sense.

Mathematica also has such mechanism. Some Mathematica functions are Listable and also allow you to perform element-wise operations on nested lists with dimensions "compatible" in some sense. Built-in listable functions are optimized for packed arrays and, similarly to NumPy's broadcasting, will give you "C-level" efficiency.

In addition to that Mathematica allows you to Compile functions with Listable RuntimeAttributes which gives you some additional control over "compatibility" of arrays. Listable compiled functions can also be easily parallelized.

There are two important differences between how NumPy's broadcasting and Mathematica's listability (compiled and not) determine if arrays are "compatible":

1. order in which dimensions are compared,
2. what happens when certain dimensions are equal 1.

NumPy starts with trailing dimensions, Mathematica - with leading. So NumPy can e.g. add arrays with dimensions {8,5,7,4} and {7,4} out of the box:

import numpy as np
(np.zeros((8,5,7,4))+np.ones((7,4))).shape
# (8, 5, 7, 4)


In Mathematica this would lead to an error:

Array[0 &, {8, 5, 7, 4}] + Array[1 &, {7, 4}];
(* Thread::tdlen: Objects of unequal length in ... cannot be combined. *)


To use listability we can transpose one of arrays to put "compatible" dimensions to the front and after addition transpose back:

Transpose[
Transpose[Array[0 &, {8, 5, 7, 4}], {3, 4, 1, 2}] +
Array[1 &, {7, 4}], {3, 4, 1, 2}
] // Dimensions
(* {8, 5, 7, 4} *)


### Listability

In contrast Mathematica can, out of the box, add arrays with dimensions {4,7,5,8} and {4,7}:

Array[0 &, {4, 7, 5, 8}] + Array[1 &, {4, 7}] // Dimensions
(* {4, 7, 5, 8} *)


which would lead to an error in NumPy

import numpy as np
(np.zeros((4,7,5,8))+np.ones((4,7)))
# Traceback (most recent call last):
#   File "<stdin>", line 1, in <module>
# ValueError: operands could not be broadcast together with shapes (4,7,5,8) (4,7)


Similarly to use broadcasting we could transpose our arrays:

import numpy as np
(np.zeros((4,7,5,8)).transpose(2,3,0,1)+np.ones((4,7))).transpose(2,3,0,1).shape
# (4, 7, 5, 8)


I don't know if this is the "correct" way to do it in NumPy. As far as I know, in contrast to Mathematica, NumPy is not copying an array on transposition, it returns a view of an array i.e. an object with information on how data from base array should be accessed. So I think that those transpositions are much cheaper than in Mathematica.

I doubt that it's possible to replicate NumPy's efficiency, on arrays which are "listability incompatible", using only top-level Mathemaica code.

As noted in comment, by @LLlAMnYP, design decision to start from leading dimensions makes, in Mathematica, more sense, since listability applies not only to full arrays, but to arbitrary nested lists.

### Compiled Listability

Since compiled functions accept only full arrays with specified rank, Compilation allows you to "split" ranks of full arrays into two parts. Last dimensions given by ranks in arguments list of Compile will be handled inside body of your compiled function, and remaining leading dimensions will be handled by Listable attribute of compiled function.

For tests let's compile simple listable function accepting two rank 2 arrays of reals:

cPlus22 = Compile[{{x, _Real, 2}, {y, _Real, 2}}, x + y, RuntimeAttributes -> {Listable}]


Now last two dimensions need to be equal since they are handled by Plus inside body of compiled function. Remaining dimensions will be handled by ordinary listability rules starting with leading ones:

cPlus22[Array[0 &, {4, 7, 5, 8}], Array[1 &, {5, 8}]] // Dimensions
(* {4, 7, 5, 8} *)
cPlus22[Array[0 &, {4, 7, 5, 8}], Array[1 &, {4, 5, 8}]] // Dimensions
(* {4, 7, 5, 8} *)
cPlus22[Array[0 &, {4, 7, 5, 8}], Array[1 &, {4, 7, 5, 8}]] // Dimensions
(* {4, 7, 5, 8} *)
cPlus22[Array[0 &, {4, 7, 5, 8}], Array[1 &, {4, 7, 3, 5, 8}]] // Dimensions
(* {4, 7, 3, 5, 8} *)


## Treating Dimensions equal to 1

When comparing consecutive dimensions NumPy's broadcasting treats them as "compatible" if they are equal, or one of them is 1. Mathematica's listability treats dimensions as "compatible" only if they are equal.

In NumPy we can do

import numpy as np
(np.zeros((1,8,1,3,7,1))+np.ones((2,1,5,3,1,4))).shape
# (2, 8, 5, 3, 7, 4)


which gives a generalized outer product.

### Outer

Mathematica has a built-in to do this kind of tasks: Outer (as noted in comment by @Sjoerd), which is "C-level efficient" when given Plus, Times and List functions and packed arrays. But Outer has its own rules for dimension "compatibility", to replicate NumPy's broadcasting conventions, all pairwise equal dimensions should be moved to the end, and dimensions equal one, that are supposed to be broadcasted, should be removed. This in general requires accessing Parts of arrays and transpositions (which in Mathematica enforces copying).

(a = Transpose[Array[0 &, {1, 8, 1, 3, 7, 1}][[1, All, 1, All, All, 1]], {1, 3, 2}]) // Dimensions
(* {8, 7, 3} *)
(b = Transpose[Array[1 &, {2, 1, 5, 3, 1, 4}][[All, 1, All, All, 1]], {1, 2, 4, 3}]) // Dimensions
(* {2, 5, 4, 3} *)
Transpose[Outer[Plus, a, b, 2, 3], {2, 5, 1, 3, 6, 4}] // Dimensions
(* {2, 8, 5, 3, 7, 4} *)


### Compiled Listability

Using different ranks in argument list of Compile results in a kind of outer product to. "Excessive" trailing dimensions of higher rank array don't have to be compatible with any dimensions of lower rank array since they will end up appended at the and of dimensions of result.

cPlus02 = Compile[{x, {y, _Real, 2}}, x + y, RuntimeAttributes -> {Listable}];
cPlus02[Array[0 &, {4, 7, 5, 8}], Array[1 &, {3, 9}]] // Dimensions
(* {4, 7, 5, 8, 3, 9} *)
cPlus02[Array[0 &, {4, 7, 5, 8}], Array[1 &, {4, 3, 9}]] // Dimensions
(* {4, 7, 5, 8, 3, 9} *)
cPlus02[Array[0 &, {4, 7, 5, 8}], Array[1 &, {4, 7, 3, 9}]] // Dimensions
(* {4, 7, 5, 8, 3, 9} *)
cPlus02[Array[0 &, {4, 7, 5, 8}], Array[1 &, {4, 7, 5, 3, 9}]] // Dimensions
(* {4, 7, 5, 8, 3, 9} *)
cPlus02[Array[0 &, {4, 7, 5, 8}], Array[1 &, {4, 7, 5, 8, 3, 9}]] // Dimensions
(* {4, 7, 5, 8, 3, 9} *)
cPlus02[Array[0 &, {4, 7, 5, 8}], Array[1 &, {4, 7, 5, 8, 2, 3, 9}]] // Dimensions
(* {4, 7, 5, 8, 2, 3, 9} *)


To emulate broadcasting in this case dimensions equal 1 should be removed, dimensions to be broadcasted from one array should be moved to beginning, and from other - to the end. Compiled function should have an argument with rank equal to number of compatible dimensions, as this argument, array with dimensions to be broadcasted at beginning, should be passed. Other argument should have rank equal to rank of array with dimensions to be broadcasted at end.

(a = Transpose[Array[0 &, {1, 8, 1, 3, 7, 1}][[1, All, 1, All, All, 1]], {1, 3, 2}]) // Dimensions
(* {8, 7, 3} *)
(b = Transpose[Array[1 &, {2, 1, 5, 3, 1, 4}][[All, 1, All, All, 1]], {2, 3, 1, 4}]) // Dimensions
(* {3, 2, 5, 4} *)
cPlus14 = Compile[{{x, _Real, 1}, {y, _Real, 4}}, x + y, RuntimeAttributes -> {Listable}];
Transpose[cPlus14[a, b], {2, 5, 4, 1, 3, 6}] // Dimensions
(* {2, 8, 5, 3, 7, 4} *)


Since compatible dimensions don't have to be handled inside body of compiled function, but can be handled by Listable attribute, there are different orderings possible. Each compatible dimension can be moved from the middle of dimensions of first array to the beginning, and rank of both arguments of compiled function can be decreased by one for each such dimension.

(a = Transpose[Array[0 &, {1, 8, 1, 3, 7, 1}][[1, All, 1, All, All, 1]], {2, 1, 3}]) // Dimensions
(* {3, 8, 7} *)
(b = Transpose[Array[1 &, {2, 1, 5, 3, 1, 4}][[All, 1, All, All, 1]], {2, 3, 1, 4}]) // Dimensions
(* {3, 2, 5, 4} *)
cPlus03 = Compile[{x, {y, _Real, 3}}, x + y, RuntimeAttributes -> {Listable}];
Transpose[cPlus03[a, b], {4, 2, 5, 1, 3, 6}] // Dimensions
(* {2, 8, 5, 3, 7, 4} *)


Below I present three approaches to broadcasting in Mathematica, with different generality and efficiency.

1. Top-level Procedural code.

It's straightforward, completely general (works for arbitrary number of lists and arbitrary function), but it's slow.

It's very fast, currently works for addition of arbitrary number of real arrays with arbitrary dimensions.

It's fastest, from presented solutions, and quite general (works for arbitrary compilable function and arbitrary number of arbitrary packable arrays with arbitrary dimensions), but it's compiled separately for each function and each "type" of arguments.

## 1. Top-level Procedural

This implementation uses dimensions of input arrays to construct proper Table expression that creates resulting array in one call by extracting proper elements from input arrays.

A helper function that constructs the Table expression:

ClearAll[broadcastingTable]
broadcastingTable[h_, f_, arrays_, dims_, maxDims_] :=
Module[{inactive, tableVars = Table[Unique["i"], Length[maxDims]]},
Prepend[
inactive[h] @@ Transpose[{tableVars, maxDims}],
inactive[Part][#1, Sequence @@ #2] &,
{
arrays,
If[#1 === 1, 1, #2] &,
] & /@ dims
}
]
] /. inactive[x_] :> x
]


Example table expression (with head replaced by Hold) for three arrays with dimensions: {4, 1, 5}, {7, 4, 3, 1} and {1, 5} looks like this:

broadcastingTable[Hold, Plus,
{arr1, arr2, arr3},
{{4, 1, 5}, {7, 4, 3, 1}, {1, 5}},
{7, 4, 3, 5}
]
(* Hold[arr1[[i4, 1, i6]] + arr2[[i3, i4, i5, 1]] + arr3[[1, i6]], {i3, 7}, {i4, 4}, {i5, 3}, {i6, 5}] *)


And now the final function:

ClearAll[broadcasted]
broadcasted::incompDims = "Objects with dimentions 1 can't be broadcasted.";
Module[{listOfLists, dims, dimColumns},
listOfLists = {lists};
dims = Dimensions /@ listOfLists;
broadcastingTable[Table, f, listOfLists, dims, Max /@ dimColumns] /;
If[MemberQ[dimColumns, dimCol_ /; ! SameQ @@ DeleteCases[dimCol, 1]],
False
(* else *),
True
]
]


It works for any function and any lists not necessary full arrays:

broadcasted[f, {a, {b, c}}, {{1}, {2}}]
(* {{f[a, 1], f[{b, c}, 1]}, {f[a, 2], f[{b, c}, 2]}} *)


For full arrays gives same results as NumPy:

broadcasted[Plus, Array[a, {2}], Array[b, {10, 2}]] // Dimensions
(* {10, 2} *)

broadcasted[Plus, Array[a, {3, 4, 1, 5, 1}], Array[b, {3, 1, 2, 1, 3}]] // Dimensions
(* {3, 4, 2, 5, 3} *)

broadcasted[Plus, Array[a, {10, 1, 5, 3}], Array[b, {2, 1, 3}], Array[# &, {5, 1}]] // Dimensions
(* {10, 2, 5, 3} *)


If dimensions are not broadcastable message is printed and function remains unevaluated:

broadcasted[Plus, Array[a, {3}], Array[b, {4, 2}]]
(* During evaluation of In[]:= broadcasted::incompDims: Objects with dimentions {{3},{4,2}} can't be broadcasted. *)
{a, a, a},
{{b[1, 1], b[1, 2]}, {b[2, 1], b[2, 2]}, {b[3, 1], b[3, 2]}, {b[4, 1], b[4, 2]}}
] *)


Here is a LibraryLink function that handles arbitrary number of arrays of reals with arbitrary dimensions.

/* broadcasting.c */
#include "WolframLibrary.h"

DLLEXPORT mint WolframLibrary_getVersion() {
return WolframLibraryVersion;
}
DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData) {
return LIBRARY_NO_ERROR;
}
DLLEXPORT void WolframLibrary_uninitialize(WolframLibraryData libData) {}

WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res
) {
switch (Argc) {
case 0:
/* At least one argument is needed. */
return LIBRARY_FUNCTION_ERROR;
case 1:
/* If one argument is given just return it. */
MArgument_setMTensor(Res, MArgument_getMTensor(Args));
return LIBRARY_NO_ERROR;
}

mint i, j;

/* ranks[i] is rank of i-th argument tensor. */
mint ranks[Argc];

/* dims[i][j] is j-th dimension of i-th argument tensor. */
const mint *(dims[Argc]);

/* data[i][j] is j-th element of i-th argument tensor. */
double *(data[Argc]);

/* Rank of result tensor. */
mint resultRank = 1;

for (i = 0; i < Argc; i++) {
MTensor tmpT = MArgument_getMTensor(Args[i]);

if (libData->MTensor_getType(tmpT) != MType_Real) {
return LIBRARY_TYPE_ERROR;
}

ranks[i] = libData->MTensor_getRank(tmpT);
dims[i] = libData->MTensor_getDimensions(tmpT);
data[i] = libData->MTensor_getRealData(tmpT);

if (resultRank < ranks[i]) {
resultRank = ranks[i];
}
}

/*
* Array of dimensions of argument tensors, with rows,
* for tensors with ranks lower than rank of result,
* filled with 1s from the beginning.
*/
mint extendedDims[Argc][resultRank];

/*
* Array of strides of argument tensors, with rows,
* for tensors with ranks lower than rank of result,
* filled with product of all tensor dimensions from the beginning.
*/
mint strides[Argc][resultRank];

/* Array of indices enumerating element of argument tensors. */
mint indices[Argc];

for (i = 0; i < Argc; i++) {
mint rankDiff = resultRank - ranks[i];

extendedDims[i][resultRank - 1] = dims[i][ranks[i] - 1];
strides[i][resultRank - 1] = extendedDims[i][resultRank - 1];
for (j = resultRank - 2; j >= rankDiff; j--) {
extendedDims[i][j] = dims[i][j - rankDiff];
strides[i][j] = strides[i][j + 1] * extendedDims[i][j];
}
for (j = rankDiff - 1; j >= 0; j--) {
extendedDims[i][j] = 1;
strides[i][j] = strides[i][rankDiff];
}

indices[i] = 0;
}

/* Dimensions of result tensor. */
mint resultDims[resultRank];

/*
* jumps[i][j] is jump of index of i-th argument tensor when index in j-th
* dimension of result tensor is incremented.
*/
mint jumps[Argc][resultRank];

/* Total number of elements in result tensor. */
mint resultElementsNumber = 1;

/* Array of indices enumerating elements of result tensor one index per dimension. */
mint resultIndices[resultRank];

for (i = resultRank - 1; i >= 0; i--) {
resultDims[i] = 1;
for (j= 0; j < Argc; j++) {
if (extendedDims[j][i] == 1) {
/*
* i-th dimension of j-th argument tensor is 1,
* so it should be broadcasted.
*/
jumps[j][i] = 1 - strides[j][i];
} else if (resultDims[i] == 1 || resultDims[i] == extendedDims[j][i]) {
/*
* i-th dimension of j-th argument tensor is not 1,
* but it's equal to all non-1 i-th dimensions of previous argument tensors,
* so i-th dimension of j-th argument tensor should be i-th dimension
* of result and it shouldn't be broadcasted.
*/
resultDims[i] = extendedDims[j][i];
jumps[j][i] = 1;
} else {
/*
* i-th dimension of j-th argument tensor is not 1,
* i-th dimension  of at least one of previous argument tensors was not 1
* and those dimensions are not equal, so tensors are not broadcastable.
*/
return LIBRARY_DIMENSION_ERROR;
}
}

resultElementsNumber *= resultDims[i];
resultIndices[i] = 0;
}

/* Returned tensor. */
MTensor resultT;
libData->MTensor_new(MType_Real, resultRank, resultDims, &resultT);

/* Actual data of returned tensor. */
double *result;
result = libData->MTensor_getRealData(resultT);

/*
* We use single loop over all elements of result array.
* resultIndices array is updated inside loop and contains indices
* corresponding to current result element as if it was accessed using one
* index per dimension, i.e. result[i] is like
* result[resultIndices][resultIndices]...[resultIndices[resultRank-1]]
* for multidimensional array.
*/
for (i = 0; i < resultElementsNumber; i++) {
mint k = resultRank - 1;
resultIndices[k]++;
while (resultIndices[k] >= resultDims[k] && k >= 1) {
resultIndices[k] = 0;
k--;
resultIndices[k]++;
}
/*
* If result would be accessed using one index per dimension,
* then current value of k would correspond to dimension which
* index was incremented in this iteration.
*/

/* At this point we know that we have at least two argument tensors. */
result[i] = data[indices] + data[indices];
indices += jumps[k];
indices += jumps[k];
for (j = 2; j < Argc; j++) {
result[i] += data[j][indices[j]];
indices[j] += jumps[j][k];
}
}

MArgument_setMTensor(Res, resultT);

return LIBRARY_NO_ERROR;
}


Save above code in broadcasting.c file in same directory as current notebook, or paste it as a string, instead of {"broadcasting.c"}, as first argument of CreateLibrary in code below. Pass, in "CompileOptions", appropriate optimization flags for your compiler, the ones below are for GCC.

Needs["CCompilerDriver"]
SetDirectory[NotebookDirectory[]];
CreateLibrary[
(* "CompileOptions" -> "-Wall -march=native -O3" *)
];
"Given arrays could not be broadcasted together.";


A helper function that loads appropriate library function for given number of array arguments.

ClearAll[loadPlusBroadcastedReal]
Quiet[
ConstantArray[{Real, _, "Constant"}, argc],
{Real, _}
],
]


Now final function that accepts arbitrary number of arrays with arbitrary dimensions, loads necessary library function, and uses it.

ClearAll[plusBroadcastedReal]


It works as expected:

plusBroadcastedReal[{1., 2.}, {{3., 4.}, {5., 6.}, {7., 8.}}]
(* {{4., 6.}, {6., 8.}, {8., 10.}} *)


If given arrays have incompatible dimensions, then an error is generated:

plusBroadcastedReal[RandomReal[{0, 1}, {4}], RandomReal[{0, 1}, {2, 3}]]
(* During evaluation of In[]:= LibraryFunction::plusBroadcastedDims: Given arrays could not be broadcasted together. >> *)
(* During evaluation of In[]:= LibraryFunction::dimerr: An error caused by inconsistent dimensions or exceeding array bounds was encountered evaluating the function plusBroadcastedReal. >> *)
(* LibraryFunctionError["LIBRARY_DIMENSION_ERROR", 3] *)


Full post exceeded maximum allowed size, so it's continued in second answer.

• This is an extremely nice answer. Two points: most NumPy functions have the parameter axis, which avoids the need to transpose in many cases; and also, I am not sure that Transpose really does copy in Mathematica, or at least not in all cases. Nov 16 '15 at 0:57
• Takes a great while to digest your answer, still not completely comprehended. But come to say a big thank you, great answer ! Especially the LibraryLink part. How long it takes you to write such a function? It is so fast. Nov 16 '15 at 16:51
• @OleksandrR. I just found that python has another crazy package "numexpr". See here, code.google.com/p/numexpr Easy to use, just import numexpr as ne, then ne.evaluate will accelerate just in time. On my computer, 40000000x2 reduced from 0.84 sec to 0.47 sec, while libfunc takes 0.74 sec. I must say, I am shocked by python's capability now. Nov 16 '15 at 16:58
• @OleksandrR. I don't know internal details of Transpose, but in all my experiments MaxMemoryUsed@Transpose[array] was equal to ByteCount[array] for packed arrays, and equal to size of new array of pointers for symbolic array. I'm not sure how conclusive is that about creating a copy. Nov 16 '15 at 22:58
• @matheorem "How long it takes you to write such a function?" I'm not sure how much it took precisely, preparing this answer was spread across whole weekend. I think, in total, it took few hours, with significant part being writing and debugging this library function. Nov 16 '15 at 23:01

Challenging NumPy's performance will be extremely difficult, and thus the effort of implementing this is not likely to be worthwhile. The reason is that the multiple-transpose method, even though it has some overhead, is already a very good way to accomplish this type of operation in Mathematica:

mat = RandomReal[1., {40000000, 2}];
vec = {1., 2.};
Transpose[vec + Transpose[mat]]; // AbsoluteTiming (* -> 1.812500 seconds *)


Your Python code takes 1.484375 seconds on my computer, so Mathematica loses out by about 25%, and not the 70% you show. In my experience, one rarely encounters situations in which a minor performance difference such as this would change anything significantly.

Let's take the following as an alternative:

cf1 = Compile[{{vec, _Real, 1}, {mat, _Real, 2}},
Table[i + vec, {i, mat}],
CompilationTarget -> "C"
];
cf2 = Compile[{{vec, _Real, 1}, {mat, _Real, 2}},
Block[{res = mat}, Do[res[[i]] += vec, {i, Length[res]}]; res],
CompilationTarget -> "C"
];


These are the fastest compiled functions that I could write. But their performance doesn't even come close:

mat = RandomReal[1., {40000000, 2}];
vec = {1., 2.};
cf1[vec, mat]; // AbsoluteTiming (* -> 4.546875 seconds *)
cf2[vec, mat]; // AbsoluteTiming (* -> 4.421875 seconds *)


They also have very significant intermediate memory consumption.

If we can't make progress with compiled code, what can we do? The next step would probably be to write one's own LibraryLink code in C++ (using e.g. Eigen), or Fortran (using a vendor BLAS, such as MKL). Of course, these libraries are intended for linear algebra applications, and so offer limited or no functions for manipulating higher-dimensional arrays. Still, one can do so efficiently, if not necessarily straightforwardly, using the carefully tuned and high-performing matrix and vector functions as primitives.

However, Mathematica also uses BLAS routines and is linked to the MKL. Some of the functions are exposed in the LinearAlgebraBLAS context (and more in LinearAlgebraLAPACK, for higher-level linear algebra operations rather than simple matrix-vector arithmetic). It is then a matter of choosing a suitable operation from the ones that are available.

GER looks useful:

$$\mathrm{GER}: \alpha, \vec{x}, \vec{y}, \mathbf{A} : \mathbf{A} \leftarrow \alpha \vec{x} {\vec{y}}^\mathrm{T} + \mathbf{A}$$

As you can see, it is a more general operation than the column-wise sum sought, so it could be adapted to other purposes beyond this one with no additional performance penalty. But note that it overwrites its input, so that for a fair test, we should first make a copy. We can use it as follows:

A = RandomReal[1., {40000000, 2}];
alpha = 1.;
x = ConstantArray[1., Length[A]];
y = {1., 2.};
Block[{A = A}, LinearAlgebraBLASGER[alpha, x, y, A]; A]; // AbsoluteTiming
(* -> 1.390625 seconds *)


Thus, we have matched (or even slightly beaten) NumPy. But, this is far from a general-purpose operation. The intention of this answer is to show that rivalling NumPy's performance is extremely difficult using only top-level Mathematica, simply because NumPy is designed to do these operations efficiently, whereas Mathematica doesn't have the same design and we as users aren't at liberty to redesign it. Instead, we have to use the tools that are made available, most of which similarly do not come close to NumPy's pure C implementation because they are not designed for this particular task. I suspect that, in most cases, to achieve comparable performance simply won't be possible without using lower-level approaches in C/C++ or Fortran.

• Of course, the design of Mathematica can't now be changed to support array broadcasting because it would seriously break existing code. NumPy supported this operation from the beginning and made a lot of design choices with that in mind. You can ask your friend who you say designed the array broadcasting functionality in NumPy what kind of trade-offs were involved. Nov 15 '15 at 18:22
• Thank you so much for this answer! And especially for introducing BLAS and LAPACK in mma, never heard this before, I should learn this. But on the other hand, I still have hope in mma to improve this. I saw this video "wolfram.com/broadcast/video.php?c=101&v=170" several days ago, Rob Knapp pointed out many times that packed array is identical to C in every aspect. And we have seen this is really the case, if not, transpose approach would not be that close to numpy. I hope mma could bring more numeric specific built-in to enhance its power in numeric calc, so we can done all in mma. Nov 16 '15 at 1:10
• I mean seamlessly switch between symbolic and numeric calculation without that worrying about the efficiency in numeric calculation Nov 16 '15 at 1:12
• It is strange. BLAS is much slower than double transpose on my computer. For 40000000x2, transpose takes 1.8sec, BLAS taks 4.8sec. What is wrong? Nov 16 '15 at 16:46
• @matheorem I tried it on a different computer and the transpose method takes 0.76s while BLAS takes 0.37s. Again BLAS is much faster. What is your computer? Nov 16 '15 at 17:42

Mathematica doesn't do that because it's ambiguous. Note that Mathematica is perfectly happy to do "broadcasting", as you call it, if the second array is transposed:

In:= {1, 2} + {{1, 2, 3}, {2, 3, 4}}
Out= {{2, 3, 4}, {4, 5, 6}}


This, in fact, gives you one way to get the result that you want:

In:= Transpose[{1, 2} + Transpose@{{1, 2, 3}, {2, 3, 4}}]
Out= {{2, 4}, {3, 5}, {4, 6}}


As to why one works, and the other doesn't, well, what is Mathematica supposed to do if you're adding a length $2$ vector to a $2 \times 2$ matrix? Should

$$[a\; b] + \begin{bmatrix} x & y \\ z & w \end{bmatrix}$$

return

$$\begin{bmatrix} a + x & a + y \\ b + z & b + w \end{bmatrix}$$

or

$$\begin{bmatrix} a + x & b + y \\ a + z & b + w \end{bmatrix}$$

In Mathematica, we can rely on the fact that it returns the former:

In:= {a, b} + {{x, y}, {z, w}}
Out= {{a + x, a + y}, {b + z, b + w}}


This means we don't need a hard-to-remember rule for special cases, and if you want the other behavior, you need to ask for it explicitly, in one way or another.

• Thank you. I already know the transpose trick. But you indeed provides an example contradicts with array broadcasting in {a, b} + {{x, y}, {z, w}}. It seems that Plus use MapThread scheme, and if mma also supports array broadcasting, then it will not know how to do {a, b} + {{x, y}, {z, w}}. Nov 12 '15 at 1:18
• And I have to complain about Tuples it only supports generating {{x, y}, {z, w},{h,k},...}, I should call this row first. While MapThead scheme used by Plus and many other built-in is column first. So we have to use Map or transpose. But Map is much slower, and also for large list Transpose also takes time which can't be neglected. Nov 12 '15 at 1:24

This is second part of answer 99553 that exceeded maximum allowed size.

In this approach we'll create functions generating optimized SymbolicC expression for given function and given "broadcasting types" of argument tensors. Here specific "broadcasting type" of tensor means: data type (Integer, Real, or Complex), rank of tensor, and positions of dimensions that are supposed to be broadcasted.

We start with patterns matching data- and tensor types and function converting Mathematica types to corresponding C types.

ClearAll[typePatt, tensorTypePatt, typeToCType]
typePatt = Integer | Real | Complex;
tensorTypePatt = {typePatt, _Integer, {___Integer} | PatternSequence[]};

typeToCType[Integer] = "mint";
typeToCType[Real] = "mreal";
typeToCType[Complex] = "mcomplex";


Simple helper function we'll use to normalize tensors "broadcasting types".

ClearAll[thirdOrEmptyList]
thirdOrEmptyList[_, _, third_: {}] := third


A function generating names of consecutive "indexed" variables:

ClearAll[indexedVars]
indexedVars = Table[#1 <> ToString[i], {i, 0, #2 - 1}] &;


Few helper functions generating common SymbolicC expressions:

Needs["SymbolicC"]

ClearAll[unequalC]
unequalC[first_, rest__] :=
COperator[Or, COperator[Unequal, {first, #}] & /@ {rest}]

ClearAll[declareAssignArrayC]
declareAssignArrayC[type_, var_, values : {__}] := {
CDeclare[type, CArray[var, Length[values]]],
MapIndexed[CAssign[CArray["resultDims", First[#2] - 1], #1] &, values]
}

ClearAll[declareAssignTensorDimsC]
declareAssignTensorDimsC[dimVar_, tensVar_] :=
CDeclare[
{"const", CPointerType["mint"]},
CAssign[
dimVar,
CCall[CPointerMember["libData", "MTensor_getDimensions"], {tensVar}]
]
]

ClearAll[declareAssignTensorDataC]
declareAssignTensorDataC[dataType:typePatt, dataVar_, tensVar_] :=
CDeclare[
CPointerType[typeToCType[dataType]],
CAssign[
dataVar,
CCall[CPointerMember["libData", "MTensor_get" <> ToString[dataType] <> "Data"], {tensVar}]
]
]

ClearAll[declareNewTensorC]
declareNewTensorC[tensorVar_, dataType:typePatt, rank_, dims_] := {
CDeclare["MTensor", tensorVar],
CCall[
CPointerMember["libData", "MTensor_new"],
{
"MType_" <> ToString[dataType],
rank,
dims,
}
]
}

ClearAll[multidimensionalTensorIndexC]
multidimensionalTensorIndexC[{}] := 0
multidimensionalTensorIndexC[{{index_, _}}] := index
multidimensionalTensorIndexC[indexDimPairs : {{_, _} ..}] :=
COperator[Plus,
MapIndexed[
COperator[Times,
Prepend[indexDimPairs[[First[#2] + 1 ;;, 2]], #1]
] &,
indexDimPairs[[All, 1]]
]
]

ClearAll[doC]
doC[body_, {iMax_}] := doC[body, {Unique["autoGeneratedIndex"], iMax}]
doC[body_, {i_, iMin_: 0, iMax_}] := {
CDeclare["mint", i],
CFor[CAssign[i, iMin], COperator[Less, {i, iMax}], COperator[PreIncrement, i],
body
]
}
doC[body_, iterators : Repeated[_List, {2, Infinity}]] :=
Fold[doC, body, Reverse@{iterators}]
doC[body_] := body

ClearAll[cMessageReturnError]
cMessageReturnError[messageTag_String, errorCode_] := {
CCall[CPointerMember["libData", "Message"], {"\"" <> messageTag <> "\""}],
CReturn[errorCode]
}


Now "core" function generating CFunction expression for our broadcasting function. It accepts name, that function will have in C code, a SymbolicC expression that will be used in inner loop to assign value to element of result tensor, list with "broadcasting types" of argument tensors, data type of result tensor, and, optionally, a boolean indicating whether dimensions of argument tensors should be tested.

ClearAll[broadcastedCFunction]
name_String,
resultElementAssignment_,
argType : {tensorTypePatt..},
retType : typePatt,
testDimensions : True | False : True
] :=
Module[
{
}
,
(* Integer being umber of argument tensors. *)
argc = Length[argType];
(* List of integers being ranks of argument tensors. *)
argRanks = argType[[All, 2]];
(* Integer being rank of result tensor. *)
resultRank = Max[argRanks];
(* List of lists of integers.
of dimensions that will be broadcasted,
in list of dimensions of i-th argument tensor. *)
(* List of lists of integers.
of dimensions that will not be broadcasted,
in list of dimensions of i-th argument tensor,
padded on the left to match rank of result. *)
Delete[Range[#1], Transpose[{#2}]] + resultRank - #1 &,
];
(* List of integers being positions,
of dimensions that will not be broadcasted in at least one argument tensor,
in list of dimensions of result tensor. *)
(* List of strings with names of variables holding data of argument tensors. *)
argDataVars = indexedVars["data", argc];
(* List of strings with names of variables holding arrays of dimentions of argument tensors. *)
argDimsVars = indexedVars["dims", argc];
(* List of lists of CArray expressions.
argNonBroadDimVarsElements[[i]] is list of symbolic C expressions,
representing non-braoadcasted dimensions, of all argument tensors,
corresponding to i-th dimension of result. *)
With[{argDimsVarRankPairs = Transpose[{argDimsVars, argRanks}]},
Function[i,
CArray[#1, i + #2 - resultRank - 1] & @@@
Extract[
argDimsVarRankPairs,
]
];
(* List of lists.
loopVarResultDimElementPairs[[i]] is a pair containing:
string with name of loop variable enumerating elements in i-th dimension of result tensor
and CArray expression representing i-th dimension of result tensor. *)
loopVarResultDimElementPairs =
MapIndexed[
{#1, CArray["resultDims", First[#2] - 1]} &,
indexedVars["i", resultRank]
];

CFunction[
{"DLLEXPORT", "int"},
name,
{
{"WolframLibraryData", "libData"},
{"mint", "Argc"},
{CPointerType["MArgument"], "Args"},
{"MArgument", "Res"}
}
,
{
CStatement@{CDeclare["MTensor", "argT"]},
(* Assign dimensions and data of argument tensors to variables. *)
CStatement@{
CAssign["argT", CCall["MArgument_getMTensor", {CArray["Args", #1]}]],
declareAssignTensorDimsC[#2, "argT"],
declareAssignTensorDataC[#4, #3, "argT"]
} &,
{Range[0, argc - 1], argDimsVars, argDataVars, argType[[All, 1]]}
],
If[testDimensions,
CStatement@{
(* Test whether all dimensions, which are expected to be 1, are really 1 *)
CIf[
unequalC[
1,
Function[pos, CArray[#1, pos - 1]] /@ #2 &,
]
],
]
(* else *),
{}
],
(* Test whether all dimensions, not declared as 1, are compatible. *)
CIf[
COperator[Or,
]
(* else *),
{}
]
]
}
(* else *),
CStatement[]
],
CStatement@declareAssignArrayC["mint", "resultDims",
ReplacePart[ConstantArray[1, resultRank],
]
],
CStatement@{
declareNewTensorC["resultT", retType, resultRank, "resultDims"],
declareAssignTensorDataC[retType, "result", "resultT"]
},
CStatement@doC[
resultElementAssignment[
CArray[#1, multidimensionalTensorIndexC[loopVarResultDimElementPairs[[#2]]]] &,
]
],
]
],
CCall["MArgument_setMTensor", {"Res", "resultT"}],
CReturn["LIBRARY_NO_ERROR"]
}
]
]
"One of dimensions intended for broadcasting is not equal to one.";
"Given arrays could not be broadcasted together.";


Function generating code of function that is supposed to be applied to elements of arrays. It returns list containing three elements: CFunction expression with definition of "element function", Function with SymbolicC expression that assigns output of element function to element of result tensor, and data type of result.

Needs["CCodeGenerator"]
ClearAll[elementFunctionSymbolicCGenerate]
elementFunctionSymbolicCGenerate[f_, argType : {tensorTypePatt ..}, opts : OptionsPattern[Compile]] :=
Module[{inactive, compileArgs, cf, retType},
compileArgs = {Unique["arg"], Blank[#1]} & @@@ argType;
cf =
inactive[Compile][compileArgs,
inactive[f] @@ compileArgs[[All, 1]],
CompilationTarget -> "MVM",
opts
] /. inactive[x_] :> x;
{
SymbolicCGenerate[cf, "elementFunction", "CodeTarget" -> "NestedFunction"] /.
CCall[setMacroName : "MArgument_setInteger" | "MArgument_setReal" | "MArgument_setComplex", {"Res", val_}] :> (
retType = Symbol[StringDrop[setMacroName, 13]];
CAssign[CDereference["Res"], val]
),
CCall["elementFunction", {"libData", Sequence @@ #2, CAddress[#1]}] &,
retType
}
]


Function generating SymbolicC code, of full library, handling particular function and argument types.

ClearAll[$libVerInitUninit]$libVerInitUninit = {
CInclude["math.h"],
CInclude["WolframRTL.h"],
CFunction[{"DLLEXPORT", "mint"}, "WolframLibrary_getVersion", {}, CReturn["WolframLibraryVersion"]],
CFunction[{"DLLEXPORT", "int"}, "WolframLibrary_initialize", {{"WolframLibraryData", "libData"}}, CReturn["LIBRARY_NO_ERROR"]],
CFunction[{"DLLEXPORT", "void"}, "WolframLibrary_uninitialize", {{"WolframLibraryData", "libData"}}, {}]
};

elFuncDef_,
name_String,
resultElementAssignment_,
argType : {tensorTypePatt ..},
retType : typePatt,
testDimensions : True | False : True
] :=
CProgram[{
$libVerInitUninit, elFuncDef, broadcastedCFunction[name, resultElementAssignment, argType, retType, testDimensions] }]  Function that compiles generated library code and loads a library function. ClearAll[compileBroadcasted] Options[compileBroadcasted] = Join[ {"TestDimensions" -> True, "LibraryName" -> Automatic}, FilterRules[Options[Compile], Except[CompilationTarget]], Options[CreateLibrary], Options[LibraryFunctionLoad] ]; compileBroadcasted[f_, argType : {tensorTypePatt ..}, opts : OptionsPattern[]] := compileBroadcasted[f, argType, opts] = Module[{name, elFuncDef, resultElementAssignment, retType}, name = Replace[OptionValue["LibraryName"], Automatic :> "broadcast_" <> IntegerString[ Hash[{f, argType, OptionValue@Options[compileBroadcasted][[All, 1]]}, "MD5"], 16, 32 ] ]; {elFuncDef, resultElementAssignment, retType} = elementFunctionSymbolicCGenerate[f, argType, FilterRules[{opts, Options[compileBroadcasted]}, Options[Compile]]]; LibraryFunctionLoad[ CreateLibrary[ ToCCodeString@broadcastingLibSymbolicCGenerate[ elFuncDef, name, resultElementAssignment, argType, retType, OptionValue["TestDimensions"] ], name, Complement[ FilterRules[{opts, Options[compileBroadcasted]}, Options[CreateLibrary]], Options[CreateLibrary] ] ], name, {#1, #2, "Constant"} & @@@ argType, {retType, Max[argType[[All, 2]]]}, FilterRules[{opts, Options[compileBroadcasted]}, Options[LibraryFunctionLoad]] ] ]  Set, in "CompileOptions", appropriate optimization flags for your compiler, the ones below are for GCC. SetOptions[compileBroadcasted, "RuntimeOptions" -> "Speed", "CompileOptions" -> "-Wall -march=native -O3" ];  Function returning "broadcasting type" of given array, i.e. a list containing: data type, depth, and positions of dimensions equal to one. ClearAll[arrayBroadcastingType] arrayBroadcastingType[array_] := With[{dims = Dimensions[array]}, { Head[Extract[array, dims]], Length[dims], Join @@ Position[dims, 1, {1}, Heads -> False] } ]  Now the final function, accepting a function and arbitrary number of packable arrays. It inspects arrays "broadcasting types", compiles a library function appropriate for given function and array types (if it's not already compiled), and uses this library function. ClearAll[broadcastedJIT] broadcastedJIT[f_, arrays__List?(ArrayQ[#, _, NumberQ] &)] := compileBroadcasted[f, arrayBroadcastingType /@ {arrays}][arrays]  Basic test that it works as expected (first evaluation may take a moment since appropriate function needs to be compiled): broadcastedJIT[Plus, {1., 2.}, {{3., 4.}, {5., 6.}, {7., 8.}}] (* {{4., 6.}, {6., 8.}, {8., 10.}} *) broadcastedJIT[Plus, {1., 2., 3.}, {{3., 4.}, {5., 6.}, {7., 8.}}] (* LibraryFunction::broadcastIncompDims: Given arrays could not be broadcasted together. >> *) (* LibraryFunction::dimerr: An error caused by inconsistent dimensions or exceeding array bounds was encountered evaluating the function broadcast_da399e884b724e81dec7fa89cbacae6d. >> *) (* LibraryFunctionError["LIBRARY_DIMENSION_ERROR", 3] *)  ## Tests Few tests that broadcasted, plusBroadcastedReal and broadcastedJIT give the same results: Replace[ SameQ[ (* When library function encounters incompatible dimensions it returns LibraryFunctionError, top level function remains unevaluated. In both cases replace result with$Failed. *)
Replace[broadcasted[Plus, ##], _broadcasted -> $Failed], Replace[plusBroadcastedReal[##], _LibraryFunctionError ->$Failed],
Replace[broadcastedJIT[Plus, ##], _LibraryFunctionError -> $Failed] ], False :> (Print[Dimensions /@ {##}]; False) ] & @@@ Map[ RandomReal[1, #] &, Join[ (* All possible pairs of array dimension lists containing from one to three elements being 1, 2 or 3. *) Tuples[Join @@ (Tuples[Range, #] & /@ Range), 2], (* Manually picked pairs and triples of array dimensions. *) { {{5, 1, 1}, {2, 8, 5, 1, 3}}, {{5, 1, 3}, {2, 8, 5, 1, 3}}, {{1, 5}, {4, 5}}, {{5, 1, 3}, {2, 8, 5, 1, 1}}, {{3, 4, 1, 5, 1}, {3, 1, 2, 1, 3}}, {{1, 3, 1, 5, 1, 7, 8, 9}, {2, 1, 4, 1, 6, 1, 8, 9}}, {{10, 1, 5, 3}, {2, 1, 3}, {5, 1}} } ], {2} ]; And @@ % (* Error messages from all functions when arrays have incompatible dimensions. *) (* True *)  ## Timings Basic timings for example data given by OP: $HistoryLength = 0;
offset = {1., 2.} // DeveloperToPackedArray;
data = RandomReal[{0, 1}, {4 10^7, 2}];
Transpose[offset + Transpose[data]] // MaxMemoryUsed // AbsoluteTiming
Outer[Plus, {offset} // DeveloperToPackedArray, data, 1] // First // MaxMemoryUsed // AbsoluteTiming
plusBroadcastedReal[offset, data] // MaxMemoryUsed // AbsoluteTiming
broadcastedJIT[Plus, offset, data] // MaxMemoryUsed // AbsoluteTiming
broadcasted[Plus, offset, data] // MaxMemoryUsed // AbsoluteTiming
(* {1.05674, 1280000488} *)
(* {1.03718, 1280000392} *)
(* {0.508851, 640000216} *)
(* {0.371809, 640000216} *)
(* {14.7058, 640003408} *)


Two approaches using built-ins: Transpose and Outer, have similar efficiency. My LibraryLink function is two times faster and uses half of memory used by built-ins, so it seems it can outperform NumPy. JIT compiled library function is almost three times faster than built-ins (note that compilation time is not included since appropriate compilation was performed in our "basic test"). Top-level procedural code is much slower, which is not surprising.

In main loop of JIT compiled library function, iterations are independent of each other, so I believe it can be easily parallelized, but I don't have hardware on which effects of parallelization could be tested. I hope that after parallelization it could compete with numpy + numexpr packages.

• Wow, you are still working on it ! This is a really a fantastic answer. Dec 17 '15 at 2:39
• But I encountered problem running it. It shows that "CreateLibrary::instl: The compiler installation directive "CompilerInstallation" -> None does not indicate a usable installation of Generic C Compiler. >>". What is wrong? I can compile with CompilationTarget->"C" properly, so I don't understand why it can't find the installation Dec 17 '15 at 2:40
• @matheorem What does Options[CreateLibrary, "CompilerInstallation"] return? Dec 17 '15 at 8:08
• @matheorem It seems CreateLibrary is handling its options in some non-standard way and default options from Options[CreateLibrary] can lead to errors with some compilers. I've changed definition of compileBroadcasted function. Does new definition work in your setup? Dec 17 '15 at 10:42
• Now, it works. Thank you so much for your work. I think it will take me considerable time to fully understand your implementation : ) Dec 17 '15 at 12:08

A top-level alternative to using Transpose twice on a big matrix is to create an appropriate matrix from the vector, and then to add the two matrices. This seems to be faster by about 30+%:

m = RandomReal[1, {4 10^7, 2}];
v = DeveloperToPackedArray[{1.,2.}];

r1 = Transpose[v + Transpose[m]]; //RepeatedTiming
r2 = ConstantArray[v, Length[m]] + m; //RepeatedTiming

r1 === r2


{1.4, Null}

{0.90, Null}

True