Counting: Easy as 1, 2, 3…

As an illustration, consider a 1-dimensional vector ofTrue andFalse for which you want to count the number of “False to True” transitions in the sequence:

>>>np.random.seed(444)>>>x=np.random.choice([False,True],size=100000)>>>xarray([ True, False,  True, ...,  True, False,  True])

With a Pythonfor loop, one way to do this would be to evaluate, in pairs, thetruth value of each element in the sequence along with the element that comes right after it:

>>>defcount_transitions(x)->int:...count=0...fori,jinzip(x[:-1],x[1:]):...ifjandnoti:...count+=1...returncount...>>>count_transitions(x)24984

In vectorized form, there’s no explicitfor loop or direct reference to the individual elements:

>>>np.count_nonzero(x[:-1]<x[1:])24984

How do these two equivalent functions compare in terms of performance? In this particular case, the vectorized NumPy call wins out by a factor of about 70 times:

>>>fromtimeitimporttimeit>>>setup='from __main__ import count_transitions, x; import numpy as np'>>>num=1000>>>t1=timeit('count_transitions(x)',setup=setup,number=num)>>>t2=timeit('np.count_nonzero(x[:-1] < x[1:])',setup=setup,number=num)>>>print('Speed difference:{:0.1f}x'.format(t1/t2))Speed difference: 71.0x

Buy Low, Sell High

Here’s another example to whet your appetite. Consider the following classictechnical interview problem:

Given a stock’s price history as a sequence, and assuming that you are only allowed to make one purchase and one sale, what is the maximum profit that can be obtained? For example, givenprices = (20, 18, 14, 17, 20, 21, 15), the max profit would be 7, from buying at 14 and selling at 21.

(To all of you finance people: no, short-selling is not allowed.)

There is a solution with n-squaredtime complexity that consists of taking every combination of two prices where the second price “comes after” the first and determining the maximum difference.

However, there is also an O(n) solution that consists of iterating through the sequence just once and finding thedifference between each price and a running minimum. It goes something like this:

>>>defprofit(prices):...max_px=0...min_px=prices[0]...forpxinprices[1:]:...min_px=min(min_px,px)...max_px=max(px-min_px,max_px)...returnmax_px>>>prices=(20,18,14,17,20,21,15)>>>profit(prices)7

Can this be done in NumPy? You bet. But first, let’s build a quasi-realistic example:

# Create mostly NaN array with a few 'turning points' (local min/max).>>>prices=np.full(100,fill_value=np.nan)>>>prices[[0,25,60,-1]]=[80.,30.,75.,50.]# Linearly interpolate the missing values and add some noise.>>>x=np.arange(len(prices))>>>is_valid=~np.isnan(prices)>>>prices=np.interp(x=x,xp=x[is_valid],fp=prices[is_valid])>>>prices+=np.random.randn(len(prices))*2

Here’s what this looks like withmatplotlib. The adage is to buy low (green) and sell high (red):

>>>importmatplotlib.pyplotasplt# Warning! This isn't a fully correct solution, but it works for now.# If the absolute min came after the absolute max, you'd have trouble.>>>mn=np.argmin(prices)>>>mx=mn+np.argmax(prices[mn:])>>>kwargs={'markersize':12,'linestyle':''}>>>fig,ax=plt.subplots()>>>ax.plot(prices)>>>ax.set_title('Price History')>>>ax.set_xlabel('Time')>>>ax.set_ylabel('Price')>>>ax.plot(mn,prices[mn],color='green',**kwargs)>>>ax.plot(mx,prices[mx],color='red',**kwargs)

What does the NumPy implementation look like? While there is nonp.cummin() “directly,” NumPy’suniversal functions (ufuncs) all have anaccumulate() method that does what its name implies:

>>>cummin=np.minimum.accumulate

Extending the logic from the pure-Python example, you can findthe difference between each price and a running minimum (element-wise), and then take the max of this sequence:

>>>defprofit_with_numpy(prices):..."""Price minus cumulative minimum price, element-wise."""...prices=np.asarray(prices)...returnnp.max(prices-cummin(prices))>>>profit_with_numpy(prices)44.2487532293278>>>np.allclose(profit_with_numpy(prices),profit(prices))True

How do these two operations, which have the sametheoretical time complexity, compare in actual runtime? First, let’s take a longer sequence. (This doesn’t necessarily need to be a time series of stock prices at this point.)

>>>seq=np.random.randint(0,100,size=100000)>>>seqarray([ 3, 23,  8, 67, 52, 12, 54, 72, 41, 10, ..., 46,  8, 90, 95, 93,       28, 24, 88, 24, 49])

Now, for a somewhat unfair comparison:

>>>setup=('from __main__ import profit_with_numpy, profit, seq;'...' import numpy as np')>>>num=250>>>pytime=timeit('profit(seq)',setup=setup,number=num)>>>nptime=timeit('profit_with_numpy(seq)',setup=setup,number=num)>>>print('Speed difference:{:0.1f}x'.format(pytime/nptime))Speed difference: 76.0x

Above, treatingprofit_with_numpy() as pseudocode (without considering NumPy’s underlying mechanics), there are actually three passes through a sequence:

• cummin(prices) has O(n) time complexity
• prices - cummin(prices) is O(n)
• max(...) is O(n)

This reduces to O(n), because O(3n) reduces to just O(n)–then “dominates” asn approaches infinity.

Therefore, these two functions have equivalentworst-case time complexity. (Although, as a side note, the NumPy function comes with significantly more space complexity.) But that is probably the least important takeaway here. One lesson is that, while theoretical time complexity is an important consideration, runtime mechanics can also play a big role. Not only can NumPy delegate to C, but with some element-wise operations andlinear algebra, it can also take advantage of computing within multiple threads. But there are a lot of factors at play here, including the underlying library used (BLAS/LAPACK/Atlas), and those details are for a whole ‘nother article entirely.

Clustering Algorithms

Machine learning is one domain that can frequently take advantage of vectorization and broadcasting. Let’s say that you have the vertices of a triangle (each row is anx, y coordinate):

>>>tri=np.array([[1,1],...[3,1],...[2,3]])

Thecentroid of this “cluster” is an(x, y) coordinate that is the arithmetic mean of each column:

>>>centroid=tri.mean(axis=0)>>>centroidarray([2.    , 1.6667])

It’s helpful to visualize this:

>>>trishape=plt.Polygon(tri,edgecolor='r',alpha=0.2,lw=5)>>>_,ax=plt.subplots(figsize=(4,4))>>>ax.add_patch(trishape)>>>ax.set_ylim([.5,3.5])>>>ax.set_xlim([.5,3.5])>>>ax.scatter(*centroid,color='g',marker='D',s=70)>>>ax.scatter(*tri.T,color='b',s=70)

Manyclustering algorithms make use of Euclidean distances of a collection of points, either to the origin or relative to their centroids.

In Cartesian coordinates, the Euclidean distance between pointsp andq is:

[source: Wikipedia]

So for the set of coordinates intri from above, the Euclidean distance of each point from the origin (0, 0) would be:

>>>np.sum(tri**2,axis=1)**0.5# Or: np.sqrt(np.sum(np.square(tri), 1))array([1.4142, 3.1623, 3.6056])

You may recognize that we are really just finding Euclidean norms:

>>>np.linalg.norm(tri,axis=1)array([1.4142, 3.1623, 3.6056])

Instead of referencing the origin, you could also find the norm of each point relative to the triangle’s centroid:

>>>np.linalg.norm(tri-centroid,axis=1)array([1.2019, 1.2019, 1.3333])

Finally, let’s take this one step further: let’s say that you have a 2d arrayX and a 2d array of multiple(x, y) “proposed” centroids. Algorithms such asK-Means clustering work by randomly assigning initial “proposed” centroids, then reassigning each data point to its closest centroid. From there, new centroids are computed, with the algorithm converging on a solution once the re-generated labels (an encoding of the centroids) are unchanged between iterations. A part of this iterative process requires computing the Euclidean distance ofeach point from each centroid:

>>>X=np.repeat([[5,5],[10,10]],[5,5],axis=0)>>>X=X+np.random.randn(*X.shape)# 2 distinct "blobs">>>centroids=np.array([[5,5],[10,10]])>>>Xarray([[ 3.3955,  3.682 ],       [ 5.9224,  5.785 ],       [ 5.9087,  4.5986],       [ 6.5796,  3.8713],       [ 3.8488,  6.7029],       [10.1698,  9.2887],       [10.1789,  9.8801],       [ 7.8885,  8.7014],       [ 8.6206,  8.2016],       [ 8.851 , 10.0091]])>>>centroidsarray([[ 5,  5],       [10, 10]])

In other words, we want to answer the question,to which centroid does each point withinX belong? We need to do some reshaping to enable broadcasting here, in order to calculate the Euclidean distance betweeneach point inX and each point incentroids:

>>>centroids[:,None]array([[[ 5,  5]],       [[10, 10]]])>>>centroids[:,None].shape(2, 1, 2)

This enables us to cleanly subtract one array from another using acombinatoric product of their rows:

>>>np.linalg.norm(X-centroids[:,None],axis=2).round(2)array([[2.08, 1.21, 0.99, 1.94, 2.06, 6.72, 7.12, 4.7 , 4.83, 6.32],       [9.14, 5.86, 6.78, 7.02, 6.98, 0.73, 0.22, 2.48, 2.27, 1.15]])

In other words, the NumPy shape ofX - centroids[:, None] is(2, 10, 2), essentially representing two stacked arrays that are each the size ofX. Next, we want thelabel (index number) of each closest centroid, finding the minimum distance on the 0th axis from the array above:

>>>np.argmin(np.linalg.norm(X-centroids[:,None],axis=2),axis=0)array([0, 0, 0, 0, 0, 1, 1, 1, 1, 1])

You can put all this together in functional form:

>>>defget_labels(X,centroids)->np.ndarray:...returnnp.argmin(np.linalg.norm(X-centroids[:,None],axis=2),...axis=0)>>>labels=get_labels(X,centroids)

Let’s inspect this visually, plotting both the two clusters and their assigned labels with a color-mapping:

>>>c1,c2=['#bc13fe','#be0119']# https://xkcd.com/color/rgb/>>>llim,ulim=np.trunc([X.min()*0.9,X.max()*1.1])>>>_,ax=plt.subplots(figsize=(5,5))>>>ax.scatter(*X.T,c=np.where(labels,c2,c1),alpha=0.4,s=80)>>>ax.scatter(*centroids.T,c=[c1,c2],marker='s',s=95,...edgecolor='yellow')>>>ax.set_ylim([llim,ulim])>>>ax.set_xlim([llim,ulim])>>>ax.set_title('One K-Means Iteration: Predicted Classes')

Amortization Tables

Vectorization has applications in finance as well.

Given an annualized interest rate, payment frequency (times per year), initial loan balance, and loan term, you can create an amortization table with monthly loan balances and payments, in a vectorized fashion. Let’s set some scalar constants first:

>>>freq=12# 12 months per year>>>rate=.0675# 6.75% annualized>>>nper=30# 30 years>>>pv=200000# Loan face value>>>rate/=freq# Monthly basis>>>nper*=freq# 360 months

NumPy comes preloaded with a handful offinancial functions that, unlike theirExcel cousins, are capable of producing vector outputs.

The debtor (or lessee) pays a constant monthly amount that is composed of a principal and interest component. As the outstanding loan balance declines, the interest portion of the total payment declines with it.

>>>periods=np.arange(1,nper+1,dtype=int)>>>principal=np.ppmt(rate,periods,nper,pv)>>>interest=np.ipmt(rate,periods,nper,pv)>>>pmt=principal+interest# Or: pmt = np.pmt(rate, nper, pv)

Next, you’ll need to calculate a monthly balance, both before and after that month’s payment, which can be defined as thefuture value of the original balance minus the future value of an annuity (a stream of payments), using a discount factord:

Functionally, this looks like:

>>>defbalance(pv,rate,nper,pmt)->np.ndarray:...d=(1+rate)**nper# Discount factor...returnpv*d-pmt*(d-1)/rate

Finally, you can drop this into a tabular format with aPandas DataFrame. Be careful with signs here.PMT is an outflow from the perspective of the debtor.

>>>importpandasaspd>>>cols=['beg_bal','prin','interest','end_bal']>>>data=[balance(pv,rate,periods-1,-pmt),...principal,...interest,...balance(pv,rate,periods,-pmt)]>>>table=pd.DataFrame(data,columns=periods,index=cols).T>>>table.index.name='month'>>>withpd.option_context('display.max_rows',6):...# Note: Using floats for $$ in production-level code = bad...print(table.round(2))...         beg_bal     prin  interest    end_balmonth1      200000.00  -172.20  -1125.00  199827.802      199827.80  -173.16  -1124.03  199654.643      199654.64  -174.14  -1123.06  199480.50...          ...      ...       ...        ...358      3848.22 -1275.55    -21.65    2572.67359      2572.67 -1282.72    -14.47    1289.94360      1289.94 -1289.94     -7.26      -0.00

At the end of year 30, the loan is paid off:

>>>final_month=periods[-1]>>>np.allclose(table.loc[final_month,'end_bal'],0)True

Image Feature Extraction

In one final example, we’ll work with an October 1941image of the USS Lexington (CV-2), the wreck of which was discovered off the coast of Australia in March 2018. First, we can map the image into a NumPy array of its pixel values:

>>>fromskimageimportio>>>url=('https://www.history.navy.mil/bin/imageDownload?image=/'...'content/dam/nhhc/our-collections/photography/images/'...'80-G-410000/80-G-416362&rendition=cq5dam.thumbnail.319.319.png')>>>img=io.imread(url,as_grey=True)>>>fig,ax=plt.subplots()>>>ax.imshow(img,cmap='gray')>>>ax.grid(False)

For simplicity’s sake, the image is loaded in grayscale, resulting in a 2d array of 64-bit floats rather than a 3-dimensionalMxNx4 RGBA array, with lower values denoting darker spots:

>>>img.shape(254, 319)>>>img.min(),img.max()(0.027450980392156862, 1.0)>>>img[0,:10]# First ten cells of the first rowarray([0.8078, 0.7961, 0.7804, 0.7882, 0.7961, 0.8078, 0.8039, 0.7922,       0.7961, 0.7961])>>>img[-1,-10:]# Last ten cells of the last rowarray([0.0784, 0.0784, 0.0706, 0.0706, 0.0745, 0.0706, 0.0745, 0.0784,       0.0784, 0.0824])

One technique commonly employed as an intermediary step in image analysis ispatch extraction. As the name implies, this consists of extracting smaller overlapping sub-arrays from a larger array and can be used in cases where it is advantageous to “denoise” or blur an image.

This concept extends to other fields, too. For example, you’d be doing something similar by taking “rolling” windows of a time series with multiple features (variables). It’s even useful for buildingConway’s Game of Life. (Although,convolution with a3x3 kernel is a more direct approach.)

Here, we will find themean of each overlapping 10x10 patch withinimg. Taking a miniature example, the first3x3 patch array in the top-left corner ofimg would be:

>>>img[:3,:3]array([[0.8078, 0.7961, 0.7804],       [0.8039, 0.8157, 0.8078],       [0.7882, 0.8   , 0.7961]])>>>img[:3,:3].mean()0.7995642701525054

The pure-Python approach to creating sliding patches would involve anestedfor loop. You’d need to consider that the starting index of the right-most patches will be at indexn - 3 + 1, wheren is the width of the array. In other words, if you were extracting 3x3 patches from a 10x10 array calledarr, the last patch taken would be fromarr[7:10, 7:10]. Also keep in mind that Python’srange() does not include itsstop parameter:

>>>size=10>>>m,n=img.shape>>>mm,nn=m-size+1,n-size+1>>>>>>patch_means=np.empty((mm,nn))>>>foriinrange(mm):...forjinrange(nn):...patch_means[i,j]=img[i:i+size,j:j+size].mean()>>>fig,ax=plt.subplots()>>>ax.imshow(patch_means,cmap='gray')>>>ax.grid(False)

With this loop, you’re performing a lot of Python calls.

An alternative that will be scalable to larger RGB or RGBA images is NumPy’sstride_tricks.

An instructive first step is to visualize, given the patch size and image shape, what a higher-dimensional array of patches would look like. We have a 2d arrayimg with shape(254, 319)and a(10, 10) 2d patch. This means our output shape (before taking the mean of each “inner”10x10 array) would be:

>>>shape=(img.shape[0]-size+1,img.shape[1]-size+1,size,size)>>>shape(245, 310, 10, 10)

You also need to specify thestrides of the new array. An array’s strides is a tuple of bytes to jump in each dimension when moving along the array. Each pixel inimg is a 64-bit (8-byte) float, meaning the total image size is254 x 319 x 8 = 648,208 bytes.

>>>img.dtypedtype('float64')>>>img.nbytes648208

Internally,img is kept in memory as one contiguous block of 648,208 bytes.strides is hence a sort of “metadata”-like attribute that tells us how many bytes we need to jump ahead to move to the next positionalong each axis. We move in blocks of 8 bytes along the rows but need to traverse8 x 319 = 2,552 bytes to move “down” from one row to another.

>>>img.strides(2552, 8)

In our case, the strides of the resulting patches will just repeat the strides ofimg twice:

>>>strides=2*img.strides>>>strides(2552, 8, 2552, 8)

Now, let’s put these pieces together with NumPy’sstride_tricks:

>>>fromnumpy.libimportstride_tricks>>>patches=stride_tricks.as_strided(img,shape=shape,strides=strides)>>>patches.shape(245, 310, 10, 10)

Here’s the first10x10 patch:

>>>patches[0,0].round(2)array([[0.81, 0.8 , 0.78, 0.79, 0.8 , 0.81, 0.8 , 0.79, 0.8 , 0.8 ],       [0.8 , 0.82, 0.81, 0.79, 0.79, 0.79, 0.78, 0.81, 0.81, 0.8 ],       [0.79, 0.8 , 0.8 , 0.79, 0.8 , 0.8 , 0.82, 0.83, 0.79, 0.81],       [0.8 , 0.79, 0.81, 0.81, 0.8 , 0.8 , 0.78, 0.76, 0.8 , 0.79],       [0.78, 0.8 , 0.8 , 0.78, 0.8 , 0.79, 0.78, 0.78, 0.79, 0.79],       [0.8 , 0.8 , 0.78, 0.78, 0.78, 0.8 , 0.8 , 0.8 , 0.81, 0.79],       [0.78, 0.77, 0.78, 0.76, 0.77, 0.8 , 0.8 , 0.77, 0.8 , 0.8 ],       [0.79, 0.76, 0.77, 0.78, 0.77, 0.77, 0.79, 0.78, 0.77, 0.76],       [0.78, 0.75, 0.76, 0.76, 0.73, 0.75, 0.78, 0.76, 0.77, 0.77],       [0.78, 0.79, 0.78, 0.78, 0.78, 0.78, 0.77, 0.76, 0.77, 0.77]])

The last step is tricky. To get a vectorized mean of each inner10x10 array, we need to think carefully about the dimensionality of what we have now. The result should collapse the last two dimensions so that we’re left with a single245x310 array.

One (suboptimal) way would be toreshapepatches first, flattening the inner 2d arrays to length-100 vectors, and then computing the mean on the final axis:

>>>veclen=size**2>>>patches.reshape(*patches.shape[:2],veclen).mean(axis=-1).shape(245, 310)

However, you can also specifyaxis as a tuple, computing a mean over the last two axes, which should be more efficient than reshaping:

>>>patches.mean(axis=(-1,-2)).shape(245, 310)

Let’s make sure this checks out by comparing equality to our looped version. It does:

>>>strided_means=patches.mean(axis=(-1,-2))>>>np.allclose(patch_means,strided_means)True

If the concept of strides has you drooling, don’t worry: Scikit-Learn has alreadyembedded this entire process nicely within itsfeature_extraction module.