4.3. Opérations#
Il est possible d’effectuer des opérations sur les tableaux numpy soit termes à termes, soit des opérations vectorielles ou matricielles.
En général, ces opérations sont à privilégier car beaucoup plus efficace (elles sont programmées dans un langage bas niveau, ce qui améliore les performances).
import numpy as np
x_4 = np.random.rand(4)
x_5 = np.random.rand(5)
A = np.random.rand(5, 4)
Aint = np.random.randint(10, size=(5, 4))
xint = np.random.randint(10, size=(4,))
4.3.1. Opérations termes à termes#
Toutes les opérations classiques sur les nombres sont possibles en utilisant les mêmes opérateurs (+, -, *, /, **, <, <=, >, >=, //, %…)
print(A)
print(A + A)
print(A * A)
[[0.01504516 0.81878973 0.33503881 0.0072234 ]
[0.21056661 0.29536325 0.66424324 0.01018915]
[0.90465145 0.95643119 0.21984221 0.33036782]
[0.7626912 0.40803054 0.68687705 0.50660092]
[0.80056929 0.52541155 0.98316403 0.90037446]]
[[0.03009033 1.63757945 0.67007763 0.01444679]
[0.42113322 0.5907265 1.32848647 0.0203783 ]
[1.8093029 1.91286238 0.43968442 0.66073563]
[1.5253824 0.81606108 1.37375409 1.01320184]
[1.60113859 1.05082311 1.96632805 1.80074892]]
[[2.26356951e-04 6.70416617e-01 1.12251006e-01 5.21774637e-05]
[4.43382975e-02 8.72394485e-02 4.41219076e-01 1.03818801e-04]
[8.18394250e-01 9.14760620e-01 4.83305978e-02 1.09142894e-01]
[5.81697867e-01 1.66488920e-01 4.71800076e-01 2.56644490e-01]
[6.40911196e-01 2.76057300e-01 9.66611503e-01 8.10674171e-01]]
print(A)
print(2.1*A)
print(A**2)
print(A[0, 0]**2, (A**2)[0, 0])
print(A+1)
[[0.01504516 0.81878973 0.33503881 0.0072234 ]
[0.21056661 0.29536325 0.66424324 0.01018915]
[0.90465145 0.95643119 0.21984221 0.33036782]
[0.7626912 0.40803054 0.68687705 0.50660092]
[0.80056929 0.52541155 0.98316403 0.90037446]]
[[0.03159484 1.71945843 0.70358151 0.01516913]
[0.44218988 0.62026282 1.39491079 0.02139722]
[1.89976805 2.0085055 0.46166864 0.69377242]
[1.60165152 0.85686413 1.4424418 1.06386193]
[1.68119552 1.10336426 2.06464446 1.89078637]]
[[2.26356951e-04 6.70416617e-01 1.12251006e-01 5.21774637e-05]
[4.43382975e-02 8.72394485e-02 4.41219076e-01 1.03818801e-04]
[8.18394250e-01 9.14760620e-01 4.83305978e-02 1.09142894e-01]
[5.81697867e-01 1.66488920e-01 4.71800076e-01 2.56644490e-01]
[6.40911196e-01 2.76057300e-01 9.66611503e-01 8.10674171e-01]]
0.00022635695149511118 0.00022635695149511118
[[1.01504516 1.81878973 1.33503881 1.0072234 ]
[1.21056661 1.29536325 1.66424324 1.01018915]
[1.90465145 1.95643119 1.21984221 1.33036782]
[1.7626912 1.40803054 1.68687705 1.50660092]
[1.80056929 1.52541155 1.98316403 1.90037446]]
print(A)
print(A+A - 2*A)
print(x_4)
print(A+x_4) # Attention, dans ce cas, `x_4` est augmenté pour avoir 2 dimensions
print(A+x_5) # ne fonctionne pas : problème de dimension
[[0.01504516 0.81878973 0.33503881 0.0072234 ]
[0.21056661 0.29536325 0.66424324 0.01018915]
[0.90465145 0.95643119 0.21984221 0.33036782]
[0.7626912 0.40803054 0.68687705 0.50660092]
[0.80056929 0.52541155 0.98316403 0.90037446]]
[[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]]
[0.09540411 0.59695434 0.92759781 0.97849711]
[[0.11044927 1.41574406 1.26263662 0.98572051]
[0.30597072 0.89231758 1.59184105 0.98868626]
[1.00005556 1.55338552 1.14744002 1.30886493]
[0.85809531 1.00498487 1.61447486 1.48509803]
[0.89597341 1.12236589 1.91076184 1.87887157]]
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[4], line 5
3 print(x_4)
4 print(A+x_4) # Attention, dans ce cas, `x_4` est augmenté pour avoir 2 dimensions
----> 5 print(A+x_5) # ne fonctionne pas : problème de dimension
ValueError: operands could not be broadcast together with shapes (5,4) (5,)
Il est possible de répéter un tableau selon un certain axe. Les commandes repeat et tile font sensiblement la même chose mais en ne répétant pas de la même manière. Par exemple
####################### repeat #######################
print(x_4)
print(x_4.repeat(2)) # répétition de chaque élément deux fois
x_41 = x_4[:, None] # augmentation d'une dimension (vecteur->matrice(4,1))
print(x_41)
print(x_41.repeat(3, axis=1)) # duplication de la colonne trois fois
x_14 = x_4[None, :] # augmentation d'une dimension (vecteur->matrice(1,4))
print(x_14)
print(x_14.repeat(3, axis=0)) # duplication de la ligne trois fois
[0.09540411 0.59695434 0.92759781 0.97849711]
[0.09540411 0.09540411 0.59695434 0.59695434 0.92759781 0.92759781
0.97849711 0.97849711]
[[0.09540411]
[0.59695434]
[0.92759781]
[0.97849711]]
[[0.09540411 0.09540411 0.09540411]
[0.59695434 0.59695434 0.59695434]
[0.92759781 0.92759781 0.92759781]
[0.97849711 0.97849711 0.97849711]]
[[0.09540411 0.59695434 0.92759781 0.97849711]]
[[0.09540411 0.59695434 0.92759781 0.97849711]
[0.09540411 0.59695434 0.92759781 0.97849711]
[0.09540411 0.59695434 0.92759781 0.97849711]]
####################### tile #######################
print(x_4)
print(np.tile(x_4, 2)) # répétition de chaque élément deux fois
print(np.tile(x_4, (3, 1))) # duplication de la ligne trois fois
x_41 = x_4[:, None] # augmentation d'une dimension (vecteur->matrice(4,1))
print(x_41)
print(np.tile(x_41, (1, 3))) # duplication de la colonne quatre fois
[0.09540411 0.59695434 0.92759781 0.97849711]
[0.09540411 0.59695434 0.92759781 0.97849711 0.09540411 0.59695434
0.92759781 0.97849711]
[[0.09540411 0.59695434 0.92759781 0.97849711]
[0.09540411 0.59695434 0.92759781 0.97849711]
[0.09540411 0.59695434 0.92759781 0.97849711]]
[[0.09540411]
[0.59695434]
[0.92759781]
[0.97849711]]
[[0.09540411 0.09540411 0.09540411]
[0.59695434 0.59695434 0.59695434]
[0.92759781 0.92759781 0.92759781]
[0.97849711 0.97849711 0.97849711]]
print(A)
print(A < 0.5)
B = A * (A < 0.5) - A * (A > 0.5)
print(B)
[[0.01504516 0.81878973 0.33503881 0.0072234 ]
[0.21056661 0.29536325 0.66424324 0.01018915]
[0.90465145 0.95643119 0.21984221 0.33036782]
[0.7626912 0.40803054 0.68687705 0.50660092]
[0.80056929 0.52541155 0.98316403 0.90037446]]
[[ True False True True]
[ True True False True]
[False False True True]
[False True False False]
[False False False False]]
[[ 0.01504516 -0.81878973 0.33503881 0.0072234 ]
[ 0.21056661 0.29536325 -0.66424324 0.01018915]
[-0.90465145 -0.95643119 0.21984221 0.33036782]
[-0.7626912 0.40803054 -0.68687705 -0.50660092]
[-0.80056929 -0.52541155 -0.98316403 -0.90037446]]
print(Aint)
print((2*Aint+1) % (Aint+1)) # cela donne bien le même résultat que Aint !
[[7 9 7 3]
[8 5 2 7]
[8 8 4 4]
[9 5 9 9]
[4 1 5 4]]
[[7 9 7 3]
[8 5 2 7]
[8 8 4 4]
[9 5 9 9]
[4 1 5 4]]
Test de vitesse
Un petit test de vitesse qui n’est pas une mesure précise mais plutôt qui donne un ordre de grandeur.
Nota Bene
Il est presque toujours plus intéressant de faire les calculs directements sur les tableaux numpy plutôt que de le faire avec des boucles !
%timeit B = 2*A
649 ns ± 6.32 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)
%%timeit
B = np.zeros(A.shape)
for i in range(A.shape[0]):
for j in range(B.shape[1]):
B[i, j] = 2*A[i, j]
4.14 µs ± 33.5 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
4.3.2. Opérations vectorielles#
Il est assez facile d’identifier les tableaux numpy à des objets algébriques comme les vecteurs ou les matrices. Attention toutefois aux vecteurs qui sont des éléments d’espaces vectoriels de dimension 1. Il est également possible d’identifier ces vecteurs à des matrices lignes ou colonnes.
Petit rappel d’algèbre linéaire
Si l’on considère un espace vectoriel \(E\) muni d’une base \((e_1,\ldots, e_n)\) et un espace vectoriel \(F\) muni d’une base \((f_1,\ldots, f_n)\) (donc \(E\) est de dimension \(n\) et \(F\) de dimension \(m\)), une application linéaire \(\phi:E\to F\) est caractérisée par l’image de la base \((e_i)\) qui s’écrit dans la base des \((f_j)\) :
\[\phi(e_i) = \sum_{j=1}^m A_{i, j} f_j, \qquad 1\leq i\leq n.\]
On appelle matrice de \(\phi\) le tableau des nombres \(A_{i,j}\).
Un vecteur \(x\) de \(E\) (resp. \(y\) de \(F\)) s’écrit dans la base des \((e_i)\) (resp. \((f_j)\))
\[
x = \sum_{i=1}^n x_i e_i, \qquad y = \sum_{j=1}^m y_j f_j.
\]
On assimile donc le vecteur \(x\) (resp. \(y\)) au tableau monodimensionnel \((x_1,\ldots, x_n)\) (resp. \((y_1,\ldots, y_m)\)).
Cependant, étant donné un vecteur \(x\) de \(E\), il est possible de considérer la forme linéaire \(\phi_x : E\to \mathbb{R}\) définie par
\[
\phi_x(a) = \sum_{i=1}^n x_i a_i, \qquad a\in E.
\]
Le vecteur \(x\) est alors assimilé à une matrice ligne avec \(x_{1, i}=x_i\), \(1\leq i\leq n\). Il est également possible de plonger \(\mathbb{R}\) dans \(E\) en utilisant l’application linéaire \(\psi_x : \mathbb{R}\to E\) définie par
\[
\psi_x(a) = (ax_1, ax_2,\ldots ax_n), \qquad a\in\mathbb{R}.
\]
Le vecteur \(x\) est alors assimilé à une matrice colonne avec \(x_{i, 1} = x_i\), \(1\leq i\leq n\).
Produit matrice vecteur
Soit \(A=(A_{i, j})_{1\leq i,j \leq N}\) et \(x=(x_i)_{1\leq i\leq N}\), le produit \(Ax\) est le vecteur défini par
\[ (Ax)_i = \sum_{j=1}^N A_{i, j} x_j, \qquad 1\leq i\leq N.\]
Attention
En python, des tableaux de forme (n, ), (n, 1) et (1, n) peuvent être utilisés pour représenter des vecteurs mathématiques (et l’on peut facilement les transformer pour passer d’une version à l’autre). Mais il fait faire attention à ce que l’on fait car les opérations dessus ne sont pas les mêmes…
Sommes et produits
La somme de deux tableaux (lorsque les tailles sont compatibles) a déjà été vue : c’est la même opération que termes à termes. On utilise le +.
Pour le produit, il faut utiliser la fonction np.matmul() ou de manière équivalente le @. Tant que l’on n’utilise pas de tableaux avec plus de 3 dimensions, on peut aussi utiliser de manière équivalente la fonction np.dot.
Pour les vecteurs, on a aussi np.inner pour le produit scalaire, np.outer pour le produit extérieur :
\[
x\cdot y = \sum_{i=1}^N x_i y_i,
\qquad
(x\otimes y)_{i, j} = x_i y_j, \quad 1\leq i,j\leq N.
\]
print(Aint)
print(xint)
print(Aint @ xint) # produit matrice-vecteur
print(np.dot(Aint, xint)) # produit matrice-vecteur
print(xint @ xint) # produit scalaire de deux vecteurs
print(np.dot(xint, xint)) # produit scalaire de deux vecteurs
print(np.inner(xint, xint)) # produit scalaire de deux vecteurs
print(np.outer(xint, xint)) # produit tensoriel de deux vecteurs
[[7 9 7 3]
[8 5 2 7]
[8 8 4 4]
[9 5 9 9]
[4 1 5 4]]
[4 4 9 7]
[148 119 128 200 93]
[148 119 128 200 93]
162
162
162
[[16 16 36 28]
[16 16 36 28]
[36 36 81 63]
[28 28 63 49]]
Transposition
On calcule la transposée d’une matrice à l’aide de la fonction np.transpose. Lorsque c’est possible (lorsque python estime que c’est possible), une vue est renvoyée sans faire de copie…
print(A)
print(A.T)
print(A.shape)
print(A.transpose().shape) # transposition d'une matrice
print(x_4, x_4.T)
x_4bis = x_4[None, :]
print(x_4bis)
print(x_4bis.T)
print(x_4bis.shape)
print(x_4bis.T.shape)
print(x_4.shape)
print(x_4.transpose().shape) # ne fait rien sur les tableaux avec une seule dimension
[[0.01504516 0.81878973 0.33503881 0.0072234 ]
[0.21056661 0.29536325 0.66424324 0.01018915]
[0.90465145 0.95643119 0.21984221 0.33036782]
[0.7626912 0.40803054 0.68687705 0.50660092]
[0.80056929 0.52541155 0.98316403 0.90037446]]
[[0.01504516 0.21056661 0.90465145 0.7626912 0.80056929]
[0.81878973 0.29536325 0.95643119 0.40803054 0.52541155]
[0.33503881 0.66424324 0.21984221 0.68687705 0.98316403]
[0.0072234 0.01018915 0.33036782 0.50660092 0.90037446]]
(5, 4)
(4, 5)
[0.09540411 0.59695434 0.92759781 0.97849711] [0.09540411 0.59695434 0.92759781 0.97849711]
[[0.09540411 0.59695434 0.92759781 0.97849711]]
[[0.09540411]
[0.59695434]
[0.92759781]
[0.97849711]]
(1, 4)
(4, 1)
(4,)
(4,)
Le module linalg
Le module linalg est un sous-module de numpy qui apporte un grand nombre de méthodes numériques adaptées à l’algèbre linéaire. Citons par exemple :
np.linalg.detpour calculer le déterminant d’une matrice ;np.linalg.solvepour résoudre un système linéaire ;np.linalg.invpour calculer l’inverse d’une matrice…
4.3.3. Sommation d’Einstein#
Le module numpy permet d’utiliser la convention d’Einstein (version implicite ou explicite) sur les tableaux. Cette convention consiste à sommer les indices répétés. Voici quelques exemples mathématiques : nous noterons \(x\), \(y\) des vecteurs et \(A\), \(B\) des matrices (les objets à plus de dimensions se traitent pareil)
\[\begin{split}
\begin{aligned}
x_i y_i &\leftrightarrow \sum_{i=0}^{n} x_i y_i, &&\text{produit scalaire}\\
A_{i,j} x_j &\leftrightarrow \sum_{j=0}^{n} A_{i,j} x_j, &&\text{produit matrice-vecteur}\\
A_{i,j} y_i &\leftrightarrow \sum_{i=0}^{n} x_i A_{i,j}, &&\text{produit vecteur-matrice}\\
A_{i,j} B_{j,k} &\leftrightarrow \sum_{j=0}^{n} A_{i,j}B_{j,k}, &&\text{produit matrice-matrice}\\
A_{i,j} A_{j,i} &\leftrightarrow \sum_{i,j=0}^{n} A_{i,j}A_{j,i}, &&\text{trace de $AA^T$}
\end{aligned}
\end{split}\]
La commande np.einsum() permet de faire ce genre d’opérations d’une manière très agréable et lisible.
help(np.einsum)
Help on _ArrayFunctionDispatcher in module numpy:
einsum(*operands, out=None, optimize=False, **kwargs)
einsum(subscripts, *operands, out=None, dtype=None, order='K',
casting='safe', optimize=False)
Evaluates the Einstein summation convention on the operands.
Using the Einstein summation convention, many common multi-dimensional,
linear algebraic array operations can be represented in a simple fashion.
In *implicit* mode `einsum` computes these values.
In *explicit* mode, `einsum` provides further flexibility to compute
other array operations that might not be considered classical Einstein
summation operations, by disabling, or forcing summation over specified
subscript labels.
See the notes and examples for clarification.
Parameters
----------
subscripts : str
Specifies the subscripts for summation as comma separated list of
subscript labels. An implicit (classical Einstein summation)
calculation is performed unless the explicit indicator '->' is
included as well as subscript labels of the precise output form.
operands : list of array_like
These are the arrays for the operation.
out : ndarray, optional
If provided, the calculation is done into this array.
dtype : {data-type, None}, optional
If provided, forces the calculation to use the data type specified.
Note that you may have to also give a more liberal `casting`
parameter to allow the conversions. Default is None.
order : {'C', 'F', 'A', 'K'}, optional
Controls the memory layout of the output. 'C' means it should
be C contiguous. 'F' means it should be Fortran contiguous,
'A' means it should be 'F' if the inputs are all 'F', 'C' otherwise.
'K' means it should be as close to the layout as the inputs as
is possible, including arbitrarily permuted axes.
Default is 'K'.
casting : {'no', 'equiv', 'safe', 'same_kind', 'unsafe'}, optional
Controls what kind of data casting may occur. Setting this to
'unsafe' is not recommended, as it can adversely affect accumulations.
* 'no' means the data types should not be cast at all.
* 'equiv' means only byte-order changes are allowed.
* 'safe' means only casts which can preserve values are allowed.
* 'same_kind' means only safe casts or casts within a kind,
like float64 to float32, are allowed.
* 'unsafe' means any data conversions may be done.
Default is 'safe'.
optimize : {False, True, 'greedy', 'optimal'}, optional
Controls if intermediate optimization should occur. No optimization
will occur if False and True will default to the 'greedy' algorithm.
Also accepts an explicit contraction list from the ``np.einsum_path``
function. See ``np.einsum_path`` for more details. Defaults to False.
Returns
-------
output : ndarray
The calculation based on the Einstein summation convention.
See Also
--------
einsum_path, dot, inner, outer, tensordot, linalg.multi_dot
einops :
similar verbose interface is provided by
`einops <https://github.com/arogozhnikov/einops>`_ package to cover
additional operations: transpose, reshape/flatten, repeat/tile,
squeeze/unsqueeze and reductions.
opt_einsum :
`opt_einsum <https://optimized-einsum.readthedocs.io/en/stable/>`_
optimizes contraction order for einsum-like expressions
in backend-agnostic manner.
Notes
-----
.. versionadded:: 1.6.0
The Einstein summation convention can be used to compute
many multi-dimensional, linear algebraic array operations. `einsum`
provides a succinct way of representing these.
A non-exhaustive list of these operations,
which can be computed by `einsum`, is shown below along with examples:
* Trace of an array, :py:func:`numpy.trace`.
* Return a diagonal, :py:func:`numpy.diag`.
* Array axis summations, :py:func:`numpy.sum`.
* Transpositions and permutations, :py:func:`numpy.transpose`.
* Matrix multiplication and dot product, :py:func:`numpy.matmul` :py:func:`numpy.dot`.
* Vector inner and outer products, :py:func:`numpy.inner` :py:func:`numpy.outer`.
* Broadcasting, element-wise and scalar multiplication, :py:func:`numpy.multiply`.
* Tensor contractions, :py:func:`numpy.tensordot`.
* Chained array operations, in efficient calculation order, :py:func:`numpy.einsum_path`.
The subscripts string is a comma-separated list of subscript labels,
where each label refers to a dimension of the corresponding operand.
Whenever a label is repeated it is summed, so ``np.einsum('i,i', a, b)``
is equivalent to :py:func:`np.inner(a,b) <numpy.inner>`. If a label
appears only once, it is not summed, so ``np.einsum('i', a)`` produces a
view of ``a`` with no changes. A further example ``np.einsum('ij,jk', a, b)``
describes traditional matrix multiplication and is equivalent to
:py:func:`np.matmul(a,b) <numpy.matmul>`. Repeated subscript labels in one
operand take the diagonal. For example, ``np.einsum('ii', a)`` is equivalent
to :py:func:`np.trace(a) <numpy.trace>`.
In *implicit mode*, the chosen subscripts are important
since the axes of the output are reordered alphabetically. This
means that ``np.einsum('ij', a)`` doesn't affect a 2D array, while
``np.einsum('ji', a)`` takes its transpose. Additionally,
``np.einsum('ij,jk', a, b)`` returns a matrix multiplication, while,
``np.einsum('ij,jh', a, b)`` returns the transpose of the
multiplication since subscript 'h' precedes subscript 'i'.
In *explicit mode* the output can be directly controlled by
specifying output subscript labels. This requires the
identifier '->' as well as the list of output subscript labels.
This feature increases the flexibility of the function since
summing can be disabled or forced when required. The call
``np.einsum('i->', a)`` is like :py:func:`np.sum(a, axis=-1) <numpy.sum>`,
and ``np.einsum('ii->i', a)`` is like :py:func:`np.diag(a) <numpy.diag>`.
The difference is that `einsum` does not allow broadcasting by default.
Additionally ``np.einsum('ij,jh->ih', a, b)`` directly specifies the
order of the output subscript labels and therefore returns matrix
multiplication, unlike the example above in implicit mode.
To enable and control broadcasting, use an ellipsis. Default
NumPy-style broadcasting is done by adding an ellipsis
to the left of each term, like ``np.einsum('...ii->...i', a)``.
To take the trace along the first and last axes,
you can do ``np.einsum('i...i', a)``, or to do a matrix-matrix
product with the left-most indices instead of rightmost, one can do
``np.einsum('ij...,jk...->ik...', a, b)``.
When there is only one operand, no axes are summed, and no output
parameter is provided, a view into the operand is returned instead
of a new array. Thus, taking the diagonal as ``np.einsum('ii->i', a)``
produces a view (changed in version 1.10.0).
`einsum` also provides an alternative way to provide the subscripts
and operands as ``einsum(op0, sublist0, op1, sublist1, ..., [sublistout])``.
If the output shape is not provided in this format `einsum` will be
calculated in implicit mode, otherwise it will be performed explicitly.
The examples below have corresponding `einsum` calls with the two
parameter methods.
.. versionadded:: 1.10.0
Views returned from einsum are now writeable whenever the input array
is writeable. For example, ``np.einsum('ijk...->kji...', a)`` will now
have the same effect as :py:func:`np.swapaxes(a, 0, 2) <numpy.swapaxes>`
and ``np.einsum('ii->i', a)`` will return a writeable view of the diagonal
of a 2D array.
.. versionadded:: 1.12.0
Added the ``optimize`` argument which will optimize the contraction order
of an einsum expression. For a contraction with three or more operands this
can greatly increase the computational efficiency at the cost of a larger
memory footprint during computation.
Typically a 'greedy' algorithm is applied which empirical tests have shown
returns the optimal path in the majority of cases. In some cases 'optimal'
will return the superlative path through a more expensive, exhaustive search.
For iterative calculations it may be advisable to calculate the optimal path
once and reuse that path by supplying it as an argument. An example is given
below.
See :py:func:`numpy.einsum_path` for more details.
Examples
--------
>>> a = np.arange(25).reshape(5,5)
>>> b = np.arange(5)
>>> c = np.arange(6).reshape(2,3)
Trace of a matrix:
>>> np.einsum('ii', a)
60
>>> np.einsum(a, [0,0])
60
>>> np.trace(a)
60
Extract the diagonal (requires explicit form):
>>> np.einsum('ii->i', a)
array([ 0, 6, 12, 18, 24])
>>> np.einsum(a, [0,0], [0])
array([ 0, 6, 12, 18, 24])
>>> np.diag(a)
array([ 0, 6, 12, 18, 24])
Sum over an axis (requires explicit form):
>>> np.einsum('ij->i', a)
array([ 10, 35, 60, 85, 110])
>>> np.einsum(a, [0,1], [0])
array([ 10, 35, 60, 85, 110])
>>> np.sum(a, axis=1)
array([ 10, 35, 60, 85, 110])
For higher dimensional arrays summing a single axis can be done with ellipsis:
>>> np.einsum('...j->...', a)
array([ 10, 35, 60, 85, 110])
>>> np.einsum(a, [Ellipsis,1], [Ellipsis])
array([ 10, 35, 60, 85, 110])
Compute a matrix transpose, or reorder any number of axes:
>>> np.einsum('ji', c)
array([[0, 3],
[1, 4],
[2, 5]])
>>> np.einsum('ij->ji', c)
array([[0, 3],
[1, 4],
[2, 5]])
>>> np.einsum(c, [1,0])
array([[0, 3],
[1, 4],
[2, 5]])
>>> np.transpose(c)
array([[0, 3],
[1, 4],
[2, 5]])
Vector inner products:
>>> np.einsum('i,i', b, b)
30
>>> np.einsum(b, [0], b, [0])
30
>>> np.inner(b,b)
30
Matrix vector multiplication:
>>> np.einsum('ij,j', a, b)
array([ 30, 80, 130, 180, 230])
>>> np.einsum(a, [0,1], b, [1])
array([ 30, 80, 130, 180, 230])
>>> np.dot(a, b)
array([ 30, 80, 130, 180, 230])
>>> np.einsum('...j,j', a, b)
array([ 30, 80, 130, 180, 230])
Broadcasting and scalar multiplication:
>>> np.einsum('..., ...', 3, c)
array([[ 0, 3, 6],
[ 9, 12, 15]])
>>> np.einsum(',ij', 3, c)
array([[ 0, 3, 6],
[ 9, 12, 15]])
>>> np.einsum(3, [Ellipsis], c, [Ellipsis])
array([[ 0, 3, 6],
[ 9, 12, 15]])
>>> np.multiply(3, c)
array([[ 0, 3, 6],
[ 9, 12, 15]])
Vector outer product:
>>> np.einsum('i,j', np.arange(2)+1, b)
array([[0, 1, 2, 3, 4],
[0, 2, 4, 6, 8]])
>>> np.einsum(np.arange(2)+1, [0], b, [1])
array([[0, 1, 2, 3, 4],
[0, 2, 4, 6, 8]])
>>> np.outer(np.arange(2)+1, b)
array([[0, 1, 2, 3, 4],
[0, 2, 4, 6, 8]])
Tensor contraction:
>>> a = np.arange(60.).reshape(3,4,5)
>>> b = np.arange(24.).reshape(4,3,2)
>>> np.einsum('ijk,jil->kl', a, b)
array([[4400., 4730.],
[4532., 4874.],
[4664., 5018.],
[4796., 5162.],
[4928., 5306.]])
>>> np.einsum(a, [0,1,2], b, [1,0,3], [2,3])
array([[4400., 4730.],
[4532., 4874.],
[4664., 5018.],
[4796., 5162.],
[4928., 5306.]])
>>> np.tensordot(a,b, axes=([1,0],[0,1]))
array([[4400., 4730.],
[4532., 4874.],
[4664., 5018.],
[4796., 5162.],
[4928., 5306.]])
Writeable returned arrays (since version 1.10.0):
>>> a = np.zeros((3, 3))
>>> np.einsum('ii->i', a)[:] = 1
>>> a
array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])
Example of ellipsis use:
>>> a = np.arange(6).reshape((3,2))
>>> b = np.arange(12).reshape((4,3))
>>> np.einsum('ki,jk->ij', a, b)
array([[10, 28, 46, 64],
[13, 40, 67, 94]])
>>> np.einsum('ki,...k->i...', a, b)
array([[10, 28, 46, 64],
[13, 40, 67, 94]])
>>> np.einsum('k...,jk', a, b)
array([[10, 28, 46, 64],
[13, 40, 67, 94]])
Chained array operations. For more complicated contractions, speed ups
might be achieved by repeatedly computing a 'greedy' path or pre-computing the
'optimal' path and repeatedly applying it, using an
`einsum_path` insertion (since version 1.12.0). Performance improvements can be
particularly significant with larger arrays:
>>> a = np.ones(64).reshape(2,4,8)
Basic `einsum`: ~1520ms (benchmarked on 3.1GHz Intel i5.)
>>> for iteration in range(500):
... _ = np.einsum('ijk,ilm,njm,nlk,abc->',a,a,a,a,a)
Sub-optimal `einsum` (due to repeated path calculation time): ~330ms
>>> for iteration in range(500):
... _ = np.einsum('ijk,ilm,njm,nlk,abc->',a,a,a,a,a, optimize='optimal')
Greedy `einsum` (faster optimal path approximation): ~160ms
>>> for iteration in range(500):
... _ = np.einsum('ijk,ilm,njm,nlk,abc->',a,a,a,a,a, optimize='greedy')
Optimal `einsum` (best usage pattern in some use cases): ~110ms
>>> path = np.einsum_path('ijk,ilm,njm,nlk,abc->',a,a,a,a,a, optimize='optimal')[0]
>>> for iteration in range(500):
... _ = np.einsum('ijk,ilm,njm,nlk,abc->',a,a,a,a,a, optimize=path)
n = 5
x = np.random.randint(10, size=(n,))
A = np.random.randint(10, size=(n, n))
print(x)
print(A)
[1 5 8 1 7]
[[1 0 4 7 3]
[3 2 2 2 7]
[6 9 3 2 1]
[4 9 2 2 8]
[3 7 9 0 8]]
# somme des coefficients
print(np.einsum('i -> ', x))
print(np.sum(x))
# produit scalaire
print(np.einsum('i, i', x, x))
print(np.inner(x, x))
# produit tensoriel
print(np.einsum('i, j', x, x))
print(np.outer(x, x))
22
22
140
140
[[ 1 5 8 1 7]
[ 5 25 40 5 35]
[ 8 40 64 8 56]
[ 1 5 8 1 7]
[ 7 35 56 7 49]]
[[ 1 5 8 1 7]
[ 5 25 40 5 35]
[ 8 40 64 8 56]
[ 1 5 8 1 7]
[ 7 35 56 7 49]]
# produit matrice-vecteur à droite
print(np.einsum('ij, j', A, x))
print(A@x)
# produit matrice-vecteur à gauche
print(np.einsum('ji, j', A, x))
print(x@A)
# produit matrice-matrice
print(np.einsum('ij, jk', A, A))
print(A@A)
print(np.einsum('ij, kj', A, A))
print(A@A.T)
print(np.einsum('ij, ij', A, A))
print(np.trace(A@A.T))
print(np.linalg.norm(A, 'fro')**2)
[ 61 80 84 123 166]
[ 61 80 84 123 166]
[ 89 140 103 35 110]
[ 89 140 103 35 110]
[[ 62 120 57 29 87]
[ 50 89 89 33 97]
[ 62 70 64 70 108]
[ 75 110 116 54 157]
[102 151 125 53 131]]
[[ 62 120 57 29 87]
[ 50 89 89 33 97]
[ 62 70 64 70 108]
[ 75 110 116 54 157]
[102 151 125 53 131]]
[[ 75 46 35 50 63]
[ 46 70 53 94 97]
[ 35 53 131 123 116]
[ 50 94 123 169 157]
[ 63 97 116 157 203]]
[[ 75 46 35 50 63]
[ 46 70 53 94 97]
[ 35 53 131 123 116]
[ 50 94 123 169 157]
[ 63 97 116 157 203]]
648
648
648.0
# trace de la matrice
print(np.einsum('ii', A))
# diagonale de la matrice
diagA = np.einsum('ii->i', A)
print(diagA) # vue de la diagonale
diagA[:] = 0 # on modifie la diagonale
print(A)
16
[1 2 3 2 8]
[[0 0 4 7 3]
[3 0 2 2 7]
[6 9 0 2 1]
[4 9 2 0 8]
[3 7 9 0 0]]
# transposition
print(np.einsum('ij->ji', A))
print(A.T)
[[0 3 6 4 3]
[0 0 9 9 7]
[4 2 0 2 9]
[7 2 2 0 0]
[3 7 1 8 0]]
[[0 3 6 4 3]
[0 0 9 9 7]
[4 2 0 2 9]
[7 2 2 0 0]
[3 7 1 8 0]]