How to implement the general array broadcasting method from NumPy?

Question

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

Answer 1

Broadcasting vs Listability

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.

Leading vs Trailing Dimensions

Broadcasting

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

Broadcasting

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} *)

---

General Broadcasting in Mathematica

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.

2. LibraryLink static function.

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

3. LibraryLink JIT compiled function.

   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[f] @@ MapThread[
                inactive[Part][#1, Sequence @@ #2] &,
                {
                    arrays,
                    MapThread[
                        If[#1 === 1, 1, #2] &,
                        {#, PadLeft[tableVars, Length[#]]}
                    ] & /@ 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.";
broadcasted[f_, lists__] :=
    Module[{listOfLists, dims, dimColumns},
        listOfLists = {lists};
        dims = Dimensions /@ listOfLists;
        dimColumns = Transpose@PadLeft[dims, Automatic, 1];
        broadcastingTable[Table, f, listOfLists, dims, Max /@ dimColumns] /;
            If[MemberQ[dimColumns, dimCol_ /; ! SameQ @@ DeleteCases[dimCol, 1]],
                Message[broadcasted::incompDims, dims];
                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. *)
(* broadcasted[Plus,
       {a[1], a[2], a[3]},
       {{b[1, 1], b[1, 2]}, {b[2, 1], b[2, 2]}, {b[3, 1], b[3, 2]}, {b[4, 1], b[4, 2]}}
   ] *)

2. LibraryLink static

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) {}

DLLEXPORT int plusBroadcastedReal(
        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[0]));
            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.
                 */
                libData->Message("plusBroadcastedDims");
                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[0]][resultIndices[1]]...[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[0][indices[0]] + data[1][indices[1]];
        indices[0] += jumps[0][k];
        indices[1] += jumps[1][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[]];
broadcastingLib =
    CreateLibrary[
        {"broadcasting.c"}, "broadcasting",
        (* "CompileOptions" -> "-Wall -march=native -O3" *)
    ];
LibraryFunction::plusBroadcastedDims =
    "Given arrays could not be broadcasted together.";

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

ClearAll[loadPlusBroadcastedReal]
loadPlusBroadcastedReal[argc_] := loadPlusBroadcastedReal[argc] =
    Quiet[
        LibraryFunctionLoad[
            broadcastingLib,
            "plusBroadcastedReal", 
            ConstantArray[{Real, _, "Constant"}, argc],
            {Real, _}
        ],
        LibraryFunction::overload
    ]

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

ClearAll[plusBroadcastedReal]
plusBroadcastedReal[arrays__] :=
    loadPlusBroadcastedReal[Length@{arrays}][arrays]

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.

Answer 2

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 LinearAlgebra`BLAS` context (and more in LinearAlgebra`LAPACK`, 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}, LinearAlgebra`BLAS`GER[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.

Answer 3

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]:= {1, 2} + {{1, 2, 3}, {2, 3, 4}}
Out[1]= {{2, 3, 4}, {4, 5, 6}}

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

In[2]:= Transpose[{1, 2} + Transpose@{{1, 2, 3}, {2, 3, 4}}]
Out[2]= {{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[3]:= {a, b} + {{x, y}, {z, w}}
Out[3]= {{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.

Answer 4

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

3. LibraryLink JIT

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, 
            CAddress[tensorVar]
        }
    ]
}

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]
broadcastedCFunction[
    name_String,
    resultElementAssignment_,
    argType : {tensorTypePatt..},
    retType : typePatt,
    testDimensions : True | False : True
] :=
   Module[
       {
           argc, argRanks, resultRank, argBroadDims, argNonBroadPadDims,
           resultNonBroadDims, argDataVars, argDimsVars, argNonBroadDimsElements, loopVarResultDimElementPairs
       }
       ,
       (* 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.
          argBroadDims[[i]] is list of positions,
          of dimensions that will be broadcasted,
          in list of dimensions of i-th argument tensor. *)
       argBroadDims = thirdOrEmptyList @@@ argType;
       (* List of lists of integers.
          argNonBroadExtDims[[i]] is list of positions,
          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. *)
       argNonBroadPadDims =
           MapThread[
               Delete[Range[#1], Transpose[{#2}]] + resultRank - #1 &,
               {argRanks, argBroadDims}
           ];
       (* 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. *)
       resultNonBroadDims = Union @@ argNonBroadPadDims;
       (* 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. *)
       argNonBroadDimsElements =
           With[{argDimsVarRankPairs = Transpose[{argDimsVars, argRanks}]},
               Function[i,
                   CArray[#1, i + #2 - resultRank - 1] & @@@
                       Extract[
                           argDimsVarRankPairs,
                           Position[argNonBroadPadDims, dims_ /; MemberQ[dims, i], {1}, Heads -> False]
                       ]
               ] /@ resultNonBroadDims
           ];
       (* 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. *)
               MapThread[
                   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 *)
                       If[MemberQ[argBroadDims, Except[{}]],
                           CIf[
                               unequalC[
                                   1,
                                   Sequence @@ Join @@ MapThread[
                                       Function[pos, CArray[#1, pos - 1]] /@ #2 &,
                                       {argDimsVars, argBroadDims}
                                   ]
                               ],
                               cMessageReturnError["broadcastNonOneDim", "LIBRARY_DIMENSION_ERROR"]
                           ]
                       (* else *),
                           {}
                       ],
                       (* Test whether all dimensions, not declared as 1, are compatible. *)
                       With[{argNonBroadDimsElementsAtLeastTwo = Select[argNonBroadDimsElements, Length[#] > 1 &]},
                           If[argNonBroadDimsElementsAtLeastTwo =!= {},
                               CIf[
                                   COperator[Or, 
                                       unequalC @@@ argNonBroadDimsElementsAtLeastTwo],
                                   cMessageReturnError["broadcastIncompDims", "LIBRARY_DIMENSION_ERROR"]
                               ]
                           (* else *),
                               {}
                           ]
                       ]
                   }
               (* else *),
                   CStatement[]
               ],
               CStatement@declareAssignArrayC["mint", "resultDims",
                   ReplacePart[ConstantArray[1, resultRank], 
                       Thread[resultNonBroadDims -> argNonBroadDimsElements[[All, 1]]]
                   ]
               ],
               CStatement@{
                   declareNewTensorC["resultT", retType, resultRank, "resultDims"],
                   declareAssignTensorDataC[retType, "result", "resultT"]
               },
               With[{loopVarDimVarPairsNonBroad = loopVarResultDimElementPairs[[resultNonBroadDims]]},
                   CStatement@doC[
                       resultElementAssignment[
                           CArray["result", multidimensionalTensorIndexC[loopVarDimVarPairsNonBroad]],
                           MapThread[
                               CArray[#1, multidimensionalTensorIndexC[loopVarResultDimElementPairs[[#2]]]] &,
                               {argDataVars, argNonBroadPadDims}
                           ]
                       ],
                       Sequence @@ loopVarDimVarPairsNonBroad
                   ]
               ],
               CCall["MArgument_setMTensor", {"Res", "resultT"}],
               CReturn["LIBRARY_NO_ERROR"]
            }
        ]
    ]
LibraryFunction::broadcastNonOneDim = 
    "One of dimensions intended for broadcasting is not equal to one.";
LibraryFunction::broadcastIncompDims = 
    "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"}}, {}]
};

ClearAll[broadcastingLibSymbolicCGenerate]
broadcastingLibSymbolicCGenerate[
    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[3], #] & /@ Range[3]), 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.} // Developer`ToPackedArray;
data = RandomReal[{0, 1}, {4 10^7, 2}];
Transpose[offset + Transpose[data]] // MaxMemoryUsed // AbsoluteTiming
Outer[Plus, {offset} // Developer`ToPackedArray, 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.

Answer 5

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 = Developer`ToPackedArray[{1.,2.}];

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

r1 === r2

{1.4, Null}

{0.90, Null}

True

Answer 6

In v13.1 (now available on Wolfram Cloud) Threaded is introduced. Sadly as a function marked as EXPERIMENTAL, the efficiency is still a problem:

data = RandomReal[{0, 1}, {4 10^5, 2}];

tst1 = Transpose[{1., 2.} + Transpose@data]; // RepeatedTiming

tst2 = Threaded@{1., 2.} + data;// RepeatedTiming

tst3 = Threaded@Developer`ToPackedArray@{1., 2.} + data;// RepeatedTiming

tst1 == tst2 == tst3