0203 Computation on arrays ufuncs

中文版：数组计算：Ufuncs

Computation on NumPy Arrays: Universal Functions

Up until now, wehavebeen discussing someofthebasic nutsandbolts of NumPy; inthenextfew sections, wewilldiveintothe reasons that NumPy is so important in the Python data science world. Namely, it provides aneasyand flexible interface to optimized computation with arrays of data.

Computation on NumPy arrays canbeveryfast, oritcanbeveryslow. Thekeyto making itfastistouse vectorized operations, generally implemented through NumPy’s universal functions (ufuncs). This section motivates theneedforNumPy’s ufuncs, which canbeusedtomake repeated calculations on array elements much more efficient. It then introduces manyofthemost common and useful arithmetic ufuncs available intheNumPy package.

The Slowness of Loops

Python’s default implementation (known as CPython) does some operations very slowly. Thisisinpartduetothe dynamic, interpreted nature of the language: thefactthattypes are flexible, so that sequences of operations cannot be compiled down to efficient machine codeasin languages likeCand Fortran. Recently there have been various attempts to address this weakness: well-known examples are the PyPy project, a just-in-time compiled implementation of Python; the Cython project, which converts Python code to compilable C code; and the Numba project, which converts snippets of Python codetofastLLVM bytecode. Eachofthese has its strengths and weaknesses, butitissafetosaythatnoneofthethree approaches has yet surpassed the reach and popularity of the standard CPython engine.

The relative sluggishness of Python generally manifests itself in situations where many small operations are being repeated – for instance looping over arrays to operate on each element. For example, imagine wehaveanarray of values and we’dliketo compute the reciprocal of each. A straightforward approach might looklikethis:

import numpy as np
np.random.seed(0)
 
def compute_reciprocals(values):
 output = np.empty(len(values))
 foriinrange(len(values)):
 output[i] = 1.0 / values[i]
 return output
 
values = np.random.randint(1, 10, size=5)
compute_reciprocals(values)

This implementation probably feels fairly natural to someone from, say, aCorJava background. Butifwe measure the execution timeofthiscodeforalarge input, weseethatthis operation isveryslow, perhaps surprisingly so! We’ll benchmark this with IPython’s %timeit magic (discussed in Profiling and Timing Code):

big_array = np.random.randint(1, 100, size=10000)
%timeit compute_reciprocals(big_array)

It takes several seconds to compute these million operations andtostore the result! Whenevencell phones have processing speeds measured in Giga-FLOPS (i.e., billions of numerical operations per second), this seems almost absurdly slow. It turns outthatthe bottleneck hereisnotthe operations themselves, butthetype-checking and function dispatches that CPython mustdoateachcycle oftheloop. Eachtimethe reciprocal is computed, Python first examines the object’stypeanddoesa dynamic lookup of the correct function touseforthattype. Ifwewere working in compiled code instead, this type specification would be known before the code executes and the result could be computed much more efficiently.

Introducing UFuncs

Formanytypes of operations, NumPy provides a convenient interface intojustthiskindof statically typed, compiled routine. Thisisknown as a vectorized operation. Thiscanbe accomplished by simply performing an operation onthearray, which willthenbe applied to each element. This vectorized approach is designed topushtheloopintothe compiled layer that underlies NumPy, leading to much faster execution.

Compare the results of the following two:

print(compute_reciprocals(values))
print(1.0 / values)

Looking at the execution timeforourbigarray, weseethatit completes orders of magnitude faster than the Python loop:

%timeit (1.0 / big_array)

Vectorized operations in NumPy are implemented via ufuncs, whose main purpose is to quickly execute repeated operations on values in NumPy arrays. Ufuncs are extremely flexible – before wesawan operation between a scalar andanarray, butwecanalso operate between two arrays:

np.arange(5) / np.arange(1, 6)

And ufunc operations are not limited to one-dimensional arrays–theycanalsoactonmulti-dimensional arrays as well:

x = np.arange(9).reshape((3, 3))
2 ** x

Computations using vectorization through ufuncs are nearly always more efficient than their counterpart implemented using Python loops, especially as the arrays growinsize. Anytimeyouseesuchaloopina Python script, you should consider whether itcanbe replaced with a vectorized expression.

Exploring NumPy’s UFuncs

Ufuncs exist in two flavors: unary ufuncs, which operate on a single input, and binary ufuncs, which operate on two inputs. We’ll see examples ofboththese types of functions here.

Array arithmetic

NumPy’s ufuncs feel very natural to use because theymakeuseof Python’s native arithmetic operators. The standard addition, subtraction, multiplication, and division canallbeused:

x = np.arange(4)
print("x =", x)
print("x + 5 =", x + 5)
print("x - 5 =", x - 5)
print("x * 2 =", x * 2)
print("x / 2 =", x / 2)
print("x // 2 =", x // 2) # floor division

There isalsoaunary ufunc for negation, and a ** operator for exponentiation, and a % operator for modulus:

print("-x = ", -x)
print("x ** 2 = ", x ** 2)
print("x % 2 = ", x % 2)

In addition, these can be strung together however you wish, and the standard order of operations is respected:

-(0.5*x + 1) ** 2

Eachofthese arithmetic operations are simply convenient wrappers around specific functions built into NumPy; for example, the + operator is a wrapper for the add function:

np.add(x, 2)

The following table lists the arithmetic operators implemented in NumPy:

| Operator | Equivalentufunc | Description | |+ |np.add |Addition (e.g., 1 + 1 = 2) | |- |np.subtract |Subtraction (e.g., 3 - 2 = 1) | |- |np.negative |Unary negation (e.g., -2) | |* |np.multiply |Multiplication (e.g., 2 * 3 = 6) | |/ |np.divide |Division (e.g., 3 / 2 = 1.5) | |// |np.floor_divide |Floor division (e.g., 3 // 2 = 1) | |** |np.power |Exponentiation (e.g., 2 ** 3 = 8) | |% |np.mod |Modulus/remainder (e.g., 9 % 4 = 1)|

Additionally there are Boolean/bitwise operators; we will explore these in Comparisons, Masks, and Boolean Logic.

Absolute value

JustasNumPy understands Python’s built-in arithmetic operators, it also understands Python’s built-in absolute value function:

x = np.array([-2, -1, 0, 1, 2])
abs(x)

The corresponding NumPy ufunc is np.absolute, which is also available under the alias np.abs:

np.abs(x)

Trigonometric functions

NumPy provides a large number of useful ufuncs, andsomeofthemost useful forthedata scientist are the trigonometric functions. We’ll start by defining an array of angles:

theta = np.linspace(0, np.pi, 3)

Nowwecan compute some trigonometric functions on these values:

print("theta = ", theta)
print("sin(theta) = ", np.sin(theta))
print("cos(theta) = ", np.cos(theta))
print("tan(theta) = ", np.tan(theta))

The values are computed to within machine precision, which is why values that should bezerodonot always hit exactly zero. Inverse trigonometric functions are also available:

x = [-1, 0, 1]
print("x = ", x)
print("arcsin(x) = ", np.arcsin(x))
print("arccos(x) = ", np.arccos(x))
print("arctan(x) = ", np.arctan(x))

Exponents and logarithms

Another common type of operation available inaNumPy ufunc are the exponentials:

x = [1, 2, 3]
print("x =", x)
print("e^x =", np.exp(x))
print("2^x =", np.exp2(x))
print("3^x =", np.power(3, x))

The inverse of the exponentials, the logarithms, are also available. The basic np.log gives the natural logarithm; if you prefer to compute the base-2 logarithm orthebase-10 logarithm, these are available as well:

x = [1, 2, 4, 10]
print("x =", x)
print("ln(x) =", np.log(x))
print("log2(x) =", np.log2(x))
print("log10(x) =", np.log10(x))

There arealsosome specialized versions that are useful for maintaining precision withverysmall input:

x = [0, 0.001, 0.01, 0.1]
print("exp(x) - 1 =", np.expm1(x))
print("log(1 + x) =", np.log1p(x))

When x isverysmall, these functions give more precise values thaniftheraw np.log or np.exp weretobeused.

Specialized ufuncs

NumPy hasmanymore ufuncs available, including hyperbolic trig functions, bitwise arithmetic, comparison operators, conversions from radians to degrees, rounding and remainders, andmuchmore. A look through the NumPy documentation reveals alotof interesting functionality.

Another excellent source for more specialized and obscure ufuncs is the submodule scipy.special. Ifyouwantto compute some obscure mathematical function onyourdata, chances areitis implemented in scipy.special. There arefartoomany functions tolistthemall, but the following snippet shows a couple that might comeupina statistics context:

from scipy import special

# Gamma functions (generalized factorials) and related functions
x = [1, 5, 10]
print("gamma(x) =", special.gamma(x))
print("ln|gamma(x)| =", special.gammaln(x))
print("beta(x, 2) =", special.beta(x, 2))

# Error function (integral of Gaussian)
# its complement, and its inverse
x = np.array([0, 0.3, 0.7, 1.0])
print("erf(x) =", special.erf(x))
print("erfc(x) =", special.erfc(x))
print("erfinv(x) =", special.erfinv(x))

There are many, many more ufuncs available inbothNumPy and scipy.special. Because the documentation of these packages is available online, a web search along the lines of “gamma function python” will generally find the relevant information.

Advanced Ufunc Features

Many NumPy users makeuseof ufuncs without ever learning their fullsetof features. We’ll outline a few specialized features of ufuncs here.

Specifying output

For large calculations, it is sometimes useful tobeableto specify the array where the result of the calculation will be stored. Rather than creating a temporary array, thiscanbeusedtowrite computation results directly to the memory location where you’dlikethemtobe. For all ufuncs, thiscanbedoneusing the out argument of the function:

x = np.arange(5)
y = np.empty(5)
np.multiply(x, 10, out=y)
print(y)

Thiscanevenbeusedwitharray views. For example, wecanwrite the results of a computation to every other element of a specified array:

y = np.zeros(10)
np.power(2, x, out=y[::2])
print(y)

Ifwehad instead written y[::2] = 2 ** x, this would have resulted in the creation of a temporary array toholdthe results of 2 ** x, followed by a second operation copying those values into the y array. This doesn’tmakemuchofa difference forsuchasmall computation, butforverylarge arrays the memory savings from careful useofthe out argument can be significant.

Aggregates

For binary ufuncs, there are some interesting aggregates thatcanbe computed directly from the object. For example, if we’dliketo reduce an array with a particular operation, wecanusethe reduce method ofanyufunc. A reduce repeatedly applies a given operation to the elements ofanarray until only a single result remains.

For example, calling reduce on the add ufunc returns thesumofall elements inthearray:

x = np.arange(1, 6)
np.add.reduce(x)

Similarly, calling reduce on the multiply ufunc results in the product ofallarray elements:

np.multiply.reduce(x)

If we’dliketostore all the intermediate results of the computation, we can instead use accumulate:

np.add.accumulate(x)

np.multiply.accumulate(x)

Notethatforthese particular cases, there are dedicated NumPy functions to compute the results (np.sum, np.prod, np.cumsum, np.cumprod), which we’ll explore in Aggregations: Min, Max, and Everything In Between.

Outer products

Finally, any ufunc can compute the output ofallpairs of two different inputs using the outer method. This allows you, inoneline, to do things like create a multiplication table:

x = np.arange(1, 6)
np.multiply.outer(x, x)

The ufunc.at and ufunc.reduceat methods, which we’ll explore in Fancy Indexing, are very helpful as well.

Another extremely useful feature of ufuncs is the ability to operate between arrays of different sizes and shapes, asetof operations known as broadcasting. This subject is important enough thatwewill devote a whole section to it (see Computation on Arrays: Broadcasting).

Ufuncs: Learning More

More information on universal functions (including thefulllistof available functions) canbefound on the NumPy and SciPy documentation websites.

Recall thatyoucanalso access information directly from within IPython by importing the packages and using IPython’s tab-completion and help (?) functionality, as described in Help and Documentation in IPython.