Diffusion vs listabilité
NumPy la diffusion vous permet d'effectuer, de manière efficace, des opérations élément par élément sur des tableaux, tant que les dimensions de ces tableaux sont considérées comme "compatibles" dans un certain sens.
Mathematica dispose également d'un tel mécanisme. Quelques Mathematica les fonctions sont Listable et vous permettent également d'effectuer des opérations élément par élément sur des listes imbriquées avec des dimensions "compatibles" dans un certain sens. Les fonctions listables intégrées sont optimisées pour les tableaux compressés et, de la même manière que NumPy , vous donnera une efficacité de "niveau C".
En plus de ça Mathematica vous permet de Compile fonctionne avec Listable RuntimeAttributes ce qui vous donne un contrôle supplémentaire sur la "compatibilité" des tableaux. Listable les fonctions compilées peuvent également être facilement parallélisées.
Il existe deux différences importantes entre la façon dont NumPy et Mathematica La listabilité (compilée ou non) de détermine si les tableaux sont "compatibles" :
- ordre dans lequel les dimensions sont comparées,
- que se passe-t-il lorsque certaines dimensions sont égales à 1 ?
Dimensions de début et de fin
Diffusion
NumPy commence par les dimensions de fin, Mathematica - avec en tête.So NumPy peut par ex. ajouter des tableaux de dimensions {8,5,7,4} et {7,4} prêt à l'emploi :
import numpy as np
(np.zeros((8,5,7,4))+np.ones((7,4))).shape
# (8, 5, 7, 4)
Dans Mathematica cela conduirait à une erreur :
Array[0 &, {8, 5, 7, 4}] + Array[1 &, {7, 4}];
(* Thread::tdlen: Objects of unequal length in ... cannot be combined. *)
Pour utiliser la listabilité, nous pouvons transposer l'un des tableaux pour mettre les dimensions "compatibles" à l'avant et après l'ajout, transposer à l'arrière :
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} *)
Lisabilité
En revanche Mathematica peut, prêt à l'emploi, ajouter des tableaux de dimensions {4,7,5,8} et {4,7} :
Array[0 &, {4, 7, 5, 8}] + Array[1 &, {4, 7}] // Dimensions
(* {4, 7, 5, 8} *)
ce qui conduirait à une erreur dans 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)
De même pour utiliser la diffusion nous pourrions transposer nos tableaux :
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)
Je ne sais pas si c'est la "bonne" façon de le faire dans NumPy .Pour autant que je sache, contrairement à Mathematica , NumPy ne copie pas un tableau lors de la transposition, il renvoie un view d'un tableau, c'est-à-dire un objet avec des informations sur la façon dont les données de base array devrait être accessible.Je pense donc que ces transpositions sont beaucoup moins chères que dans Mathematica .
Je doute qu'il soit possible de répliquer NumPy l'efficacité de , sur les tableaux qui sont "incompatibles avec la listabilité", en utilisant uniquement Mathemaica de niveau supérieur code.
Comme indiqué dans le commentaire, par @LLlAMnYP, la décision de conception de partir des dimensions principales fait, dans Mathematica , plus de sens, puisque la listabilité s'applique non seulement aux tableaux complets, mais aussi aux listes imbriquées arbitraires.
Listabilité compilée
Étant donné que les fonctions compilées n'acceptent que des tableaux complets avec un rang spécifié, la compilation vous permet de "diviser" les rangs de tableaux complets en deux parties. Dernières dimensions données par les rangs dans la liste des arguments de Compile sera géré à l'intérieur du corps de votre fonction compilée, et les dimensions principales restantes seront gérées par Listable attribut de la fonction compilée.
Pour les tests, compilons une simple fonction listable acceptant deux tableaux de réels de rang 2 :
cPlus22 = Compile[{{x, _Real, 2}, {y, _Real, 2}}, x + y, RuntimeAttributes -> {Listable}]
Maintenant, les deux dernières dimensions doivent être égales car elles sont gérées par Plus à l'intérieur du corps de la fonction compilée. Les dimensions restantes seront gérées par des règles de listabilité ordinaires en commençant par les premières :
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} *)
Traitement des dimensions égales à 1
Diffusion
Lors de la comparaison de dimensions consécutives NumPy la diffusion de les traite comme "compatibles" s'ils sont égaux, ou si l'un d'eux vaut 1.Mathematica La listabilité de traite les dimensions comme "compatibles" uniquement si elles sont égales.
Dans NumPy nous pouvons faire
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)
ce qui donne un produit extérieur généralisé.
Externe
Mathematica a une fonction intégrée pour faire ce genre de tâches :Outer (comme indiqué dans le commentaire de @Sjoerd), qui est "efficace au niveau C" lorsqu'il est donné Plus , Times et List fonctions et tableaux compressés.Mais Outer a ses propres règles de "compatibilité" de dimension, pour répliquer NumPy Selon les conventions de diffusion de , toutes les dimensions égales par paire doivent être déplacées vers la fin, et les dimensions égales à un, qui sont censées être diffusées, doivent être supprimées. Cela nécessite en général d'accéder à Part s de tableaux et de transpositions (qui dans Mathematica impose la copie).
(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} *)
Listabilité compilée
Utilisation de différents rangs dans la liste d'arguments de Compile se traduit par une sorte de produit extérieur à. Les dimensions de fin "excessives" d'un tableau de rang supérieur ne doivent pas nécessairement être compatibles avec les dimensions d'un tableau de rang inférieur, car elles se retrouveront ajoutées aux dimensions et des résultats.
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} *)
Pour émuler la diffusion dans ce cas, les dimensions sont égales à 1 doivent être supprimés, les dimensions à diffuser d'un tableau doivent être déplacées au début, et de l'autre - à la fin. La fonction compilée doit avoir un argument avec un rang égal au nombre de dimensions compatibles, car cet argument, tableau avec des dimensions à diffuser au début, doit être passé. L'autre argument doit avoir un rang égal au rang du tableau avec des dimensions à diffuser à la fin.
(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} *)
Étant donné que les dimensions compatibles ne doivent pas être gérées dans le corps de la fonction compilée, mais peuvent être gérées par Listable attribut, il existe différents classements possibles. Chaque dimension compatible peut être déplacée du milieu des dimensions du premier tableau vers le début, et le rang des deux arguments de la fonction compilée peut être diminué de un pour chacune de ces dimensions.
(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} *)
Diffusion générale dans Mathematica
Ci-dessous, je présente trois approches de diffusion dans Mathematica , avec une généralité et une efficacité différentes.
-
Procédure de niveau supérieur code.
C'est simple, complètement général (fonctionne pour un nombre arbitraire de listes et une fonction arbitraire), mais c'est lent.
-
LibraryLink statique fonction.
C'est très rapide, fonctionne actuellement pour l'ajout d'un nombre arbitraire de tableaux réels avec des dimensions arbitraires.
-
LibraryLink JIT fonction compilée.
C'est le plus rapide, à partir des solutions présentées, et assez général (fonctionne pour une fonction compilable arbitraire et un nombre arbitraire de tableaux empaquetables arbitraires avec des dimensions arbitraires), mais il est compilé séparément pour chaque fonction et chaque "type" d'arguments.
1. Procédure de niveau supérieur
Cette implémentation utilise les dimensions des tableaux d'entrée pour construire le Table approprié expression qui crée le tableau résultant en un seul appel en extrayant les éléments appropriés des tableaux d'entrée.
Une fonction d'assistance qui construit le 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
]
Exemple d'expression de table (avec head remplacé par Hold ) pour trois tableaux de dimensions :{4, 1, 5} , {7, 4, 3, 1} et {1, 5} ressemble à ceci :
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}] *)
Et maintenant la fonction finale :
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 = [email protected][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
]
]
Cela fonctionne pour n'importe quelle fonction et n'importe quelle liste de tableaux complets non nécessaires :
broadcasted[f, {a, {b, c}}, {{1}, {2}}]
(* {{f[a, 1], f[{b, c}, 1]}, {f[a, 2], f[{b, c}, 2]}} *)
Pour les tableaux complets donne les mêmes résultats que 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} *)
Si les dimensions ne sont pas diffusables le message est imprimé et la fonction reste non évaluée :
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 statique
Voici une fonction LibraryLink qui gère un nombre arbitraire de tableaux de réels avec des dimensions arbitraires.
/* 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;
}
Enregistrer le code ci-dessus dans broadcasting.c fichier dans le même répertoire que le bloc-notes actuel ou collez-le sous forme de chaîne au lieu de {"broadcasting.c"} , comme premier argument de CreateLibrary dans le code ci-dessous. Passe, en "CompileOptions" , les drapeaux d'optimisation appropriés pour votre compilateur, ceux ci-dessous sont pour GCC .
Needs["CCompilerDriver`"]
SetDirectory[NotebookDirectory[]];
broadcastingLib =
CreateLibrary[
{"broadcasting.c"}, "broadcasting",
(* "CompileOptions" -> "-Wall -march=native -O3" *)
];
LibraryFunction::plusBroadcastedDims =
"Given arrays could not be broadcasted together.";
Une fonction d'assistance qui charge la fonction de bibliothèque appropriée pour un nombre donné d'arguments de tableau.
ClearAll[loadPlusBroadcastedReal]
loadPlusBroadcastedReal[argc_] := loadPlusBroadcastedReal[argc] =
Quiet[
LibraryFunctionLoad[
broadcastingLib,
"plusBroadcastedReal",
ConstantArray[{Real, _, "Constant"}, argc],
{Real, _}
],
LibraryFunction::overload
]
Maintenant, la fonction finale qui accepte un nombre arbitraire de tableaux avec des dimensions arbitraires, charge la fonction de bibliothèque nécessaire et l'utilise.
ClearAll[plusBroadcastedReal]
plusBroadcastedReal[arrays__] :=
loadPlusBroadcastedReal[[email protected]{arrays}][arrays]
Cela fonctionne comme prévu :
plusBroadcastedReal[{1., 2.}, {{3., 4.}, {5., 6.}, {7., 8.}}]
(* {{4., 6.}, {6., 8.}, {8., 10.}} *)
Si des tableaux donnés ont des dimensions incompatibles, alors une erreur est générée :
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] *)
Le message complet a dépassé la taille maximale autorisée, il continue donc dans la deuxième réponse.
Il sera extrêmement difficile de remettre en question les performances de NumPy et, par conséquent, l'effort de mise en œuvre n'en vaudra probablement pas la peine. La raison en est que la méthode de transposition multiple, même si elle a une certaine surcharge, est déjà un très bon moyen d'accomplir ce type d'opération dans Mathematica :
mat = RandomReal[1., {40000000, 2}];
vec = {1., 2.};
Transpose[vec + Transpose[mat]]; // AbsoluteTiming (* -> 1.812500 seconds *)
Votre code Python prend 1,484375 secondes sur mon ordinateur, donc Mathematica perd environ 25 %, et non les 70 % que vous affichez. D'après mon expérience, on rencontre rarement des situations dans lesquelles une différence de performance mineure comme celle-ci changerait quoi que ce soit de manière significative.
Prenons ce qui suit comme 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"
];
Ce sont les fonctions compilées les plus rapides que j'ai pu écrire. Mais leurs performances sont loin d'être proches :
mat = RandomReal[1., {40000000, 2}];
vec = {1., 2.};
cf1[vec, mat]; // AbsoluteTiming (* -> 4.546875 seconds *)
cf2[vec, mat]; // AbsoluteTiming (* -> 4.421875 seconds *)
Ils ont également une consommation de mémoire intermédiaire très importante.
Si nous ne pouvons pas progresser avec le code compilé, que pouvons-nous faire ? La prochaine étape serait probablement d'écrire son propre code LibraryLink en C++ (en utilisant par exemple Eigen) ou en Fortran (en utilisant un fournisseur BLAS, tel que MKL). Bien sûr, ces bibliothèques sont destinées aux applications d'algèbre linéaire, et offrent donc des fonctions limitées ou inexistantes pour manipuler des tableaux de plus grande dimension. Pourtant, on peut le faire efficacement, sinon nécessairement directement, en utilisant les fonctions matricielles et vectorielles soigneusement réglées et hautement performantes comme primitives.
Cependant, Mathematica utilise également les routines BLAS et est lié au MKL. Certaines des fonctions sont exposées dans le LinearAlgebra`BLAS` contexte (et plus dans LinearAlgebra`LAPACK` , pour des opérations d'algèbre linéaire de plus haut niveau plutôt qu'une simple arithmétique matrice-vecteur). Il s'agit alors de choisir une opération adaptée parmi celles qui sont disponibles.
GER semble utile :
$$ \mathrm{GER} :\alpha, \vec{x}, \vec{y}, \mathbf{A} :\mathbf{A} \leftarrow \alpha \vec{x} {\vec{y}} ^\mathrm{T} + \mathbf{A} $$
Comme vous pouvez le voir, il s'agit d'une opération plus générale que la somme par colonne recherchée, elle pourrait donc être adaptée à d'autres fins au-delà de celle-ci sans pénalité de performance supplémentaire. Mais notez qu'il écrase son entrée, de sorte que pour un test équitable, nous devrions d'abord faire une copie. Nous pouvons l'utiliser comme suit :
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 *)
Ainsi, nous avons égalé (ou même légèrement battu) NumPy. Mais, c'est loin d'être une opération à usage général. L'intention de cette réponse est de montrer qu'il est extrêmement difficile de rivaliser avec les performances de NumPy en utilisant uniquement Mathematica de haut niveau. , simplement parce que NumPy est conçu pour effectuer ces opérations efficacement, alors que Mathematica n'a pas le même design et nous, en tant qu'utilisateurs, ne sommes pas libres de le reconcevoir. Au lieu de cela, nous devons utiliser les outils mis à disposition, dont la plupart ne se rapprochent pas non plus de l'implémentation C pure de NumPy car ils ne sont pas conçus pour cette tâche particulière. Je soupçonne que, dans la plupart des cas, il ne sera tout simplement pas possible d'obtenir des performances comparables sans utiliser des approches de niveau inférieur en C/C++ ou Fortran.
Mathematica ne le fait pas parce que c'est ambigu. Notez que Mathematica se contente parfaitement de faire de la "diffusion", comme vous l'appelez, si le deuxième tableau est transposé :
In[1]:= {1, 2} + {{1, 2, 3}, {2, 3, 4}}
Out[1]= {{2, 3, 4}, {4, 5, 6}}
Ceci, en fait, vous donne un moyen d'obtenir le résultat que vous souhaitez :
In[2]:= Transpose[{1, 2} + [email protected]{{1, 2, 3}, {2, 3, 4}}]
Out[2]= {{2, 4}, {3, 5}, {4, 6}}
Quant à savoir pourquoi l'un fonctionne et l'autre pas, eh bien, qu'est-ce que Mathematica est censé faire si vous ajoutez un vecteur de longueur $ 2 $ à une matrice $ 2 \times 2 $ ? Devrait
$$ [un\ ; b] + \begin{bmatrice} x &y \\ z &w \end{bmatrice} $$
retour
$$ \begin{bmatrice} a + x &a + y \\ b + z &b + w \end{bmatrice} $$
ou
$$ \begin{bmatrice} a + x &b + y \\ a + z &b + w \end{bmatrice} $$
Dans Mathematica, nous pouvons compter sur le fait qu'il renvoie le premier :
In[3]:= {a, b} + {{x, y}, {z, w}}
Out[3]= {{a + x, a + y}, {b + z, b + w}}
Cela signifie que nous n'avons pas besoin d'une règle difficile à retenir pour les cas particuliers, et si vous voulez l'autre comportement, vous devez le demander explicitement, d'une manière ou d'une autre.
