<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/machine-learning-explained-vectorization-and-matrix-operations/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://enhancedatascience.com/2018/05/07/machine-learning-explained-vectorization-matrix-operations/"> Enhance Data Science</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Today in <a rel="nofollow" href="http://enhancedatascience.com/tag/machine-learning-explained/" target="_blank">Machine Learning Explained</a>, we will tackle a central (yet under-looked) aspect of Machine Learning: <strong>vectorization</strong>. Let’s say you want to compute the sum of the values of an array. The naive way to do so is to loop over the elements and to sequentially sum them. This naive way is slow and tends to get even slower with large amounts of data and large data structures.<br>
With vectorization these operations can be seen as matrix operations which are often more efficient than standard loops. Vectorized versions of algorithm are several orders of magnitudes faster and are easier to understand from a mathematical perspective.</p>
<h2> A basic exemple of vectorization </h2>
<h3>Preliminary exemple – Python </h3>
<p> <br>
Let’s compare the naive way and the vectorized way of computing the sum of the elements of an array. To do so, we will create a large (100,000 elements) Numpy array and compute the sum of its element 1,000 times with each algorithm. The overall computation time will then be compared.</p>
<pre>
import numpy as np
import time
W=np.random.normal(0,1,100000)
n_rep=1000
</pre>
<p>The naive way to compute the sum iterates over all the elements of the array and stores the sum:</p>
<pre>
start_time = time.time()
for i in range(n_rep):
    loop_res=0
    for elt in W:
        loop_res+=elt
time_loop = time.time() - start_time
</pre>
<p>If <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-61e78bd9c83a78e47b60a97ff56ab794_l3.png?resize=21%2C13" alt="mathbf{W}" width="21" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="13"> is our vector of interest, the sum of its elements can be expressed as the dot product <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-bdb7db25bae87d00c16dc718d125ac7a_l3.png?resize=109%2C20" alt="mathbf{1}^T mathbf{W} = sum W_i" width="109" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="20">:</p>
<pre>
start_time = time.time()
for i in range(n_rep):
    one_dot=np.ones(W.shape)
    vect_res=one_dot.T.dot(W)
time_vect = time.time() - start_time
</pre>
<p>Finally, we can check that both methods yield the same results and compare their runtime. The vectorized version run approximately 100 to 200 times faster than the naive loop.</p>
<pre>
print(np.abs(vect_res-loop_res)&lt;10e-10)
print(time_loop/time_vect)
</pre>
<p>Note: The same results can be obtained with <em>np.sum</em>.The numpy version has a very similar runtime to our vectorized version. Numpy being very optimized, this show that our vectorized sum is reasonably fast.</p>
<pre>
start_time = time.time()
for i in range(n_rep):
    vect_res=np.sum(W)
time_vect_np = time.time() - start_time
</pre>
<h3>Preliminary exemple – R </h3>
<p>The previous experiments can be replicated in R:</p>
<pre>
##Creation of the vector
W=matrix(rnorm(100000))
n_rep=10000
#Naive way:
library(tictoc)
tic('Naive computation')
for (rep in 1:n_rep)
{
  res_loop=0
  for (w_i in W)
    res_loop=w_i+res_loop
}
toc()

tic('Vectorized computation')
# vectorized way
for (rep in 1:n_rep)
{
  ones=rep(1,length(W))
  res_vect= crossprod(ones,W)
}
toc()

tic('built-in computation')
# built-in way
for (rep in 1:n_rep)
{
  res_built_in= sum(W)
}
toc()
</pre>
<p>In R, the vectorized version is only an order of magnitude faster than the naive way. The built-in way achieves the best performances and is an order of magnitude faster than our vectorized way.</p>
<h3>Preliminary exemple – Results </h3>
<p><span>   </span><span>   </span><img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-40f2631dae68da4963fc2dc6b9da6463_l3.png?resize=320%2C68" width="320" alt="[begin{tabular}{|l|c|c|r|}   hline   Lang &amp; Naive &amp; Vectorized &amp; Built in \   hline   Python &amp; 115 s &amp; 0.55 s &amp; 0.28 s \   R &amp; 25 s &amp; 3.8 s &amp; 0.9 s \   hline end{tabular}]" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="68"></p>
<p><span>   </span><span>   </span><img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-593f085906083f8cd848e735d395f3f6_l3.png?resize=578%2C18" width="450" alt="[text{Computation time of the sum of an array of 100,000 elements (10,000 rep)}]" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="18"></p>
<p>Vectorization divides the computation times by several order of magnitudes and the difference with loops increase with the size of the data. Hence, if you want to deal with large amount of data, rewriting the algorithm as matrix operations may lead to important performances gains.</p>
<h2> Why vectorization is (often) faster </h2>
<ul>
<li> R and Python are <strong>interpreted language</strong>, this means that your instructions are analyzed and interpreted at each execution. Since they are not statically typed, the loops have to assess the type of the operands at each iteration which leads to a computational overhead.</li>
<li> R and Python linear algebra relies on <strong>optimized back-end for matrix operations and linear algebra</strong>. These back-end are written in C/C++ and can process loops efficiently. Furthermore, while loops are sequential, these back-end can run operations in parallel which improves the computation speed on modern CPU.</li>
</ul>
<p><strong>Note 1</strong>: Though vectorization is often faster, it requires to allocate the memory of the array. If your amount of RAM is limited or if the amount of data is large, loops may be required.<br>
<strong>Note 2</strong>: When you deal with large arrays or computationally intensive algebra ( like inversion, diagonalization, eigenvalues computations, ….) computations on GPU are even order of magnitudes faster than on CPU. To write efficient GPU code, the code needs to be composed of matrix operations. Hence, having vectorized code maked it easier to translate CPU code to GPU (or tensor-based frameworks).</p>
<h2> A small zoo of matrix operations </h2>
<p>The goal of this part is to show some basic matrix operations/vectorization and to end on a more complex example to show the thought process which underlies vectorization.</p>
<h3> Column-wise matrix sum</h3>
<p>The column wise sum (and mean) can be expressed as a matrix product. Let <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-1a62e250fd3bdd137e35e78ef18c41f6_l3.png?resize=19%2C13" alt="mathbf{M}" width="19" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="13"> be our matrix of interest. Using the matrix multiplication formula, we have: <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-b0854e8d457807277a1e27a5a9fc17bf_l3.png?resize=242%2C21" alt="[mathbf{1}^Tmathbf{M}]_j= sum_i mathbf{1}_i mathbf{M}_{i,j} =sum_i mathbf{M}_{i,j}" width="242" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="21">. Hence, the column-wise sum of <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-1a62e250fd3bdd137e35e78ef18c41f6_l3.png?resize=19%2C13" alt="mathbf{M}" width="19" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="13"> is <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-28265f3a34cf221b1d8bf887029bd018_l3.png?resize=39%2C16" alt="mathbf{1}^Tmathbf{M}" width="39" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="16">.</p>
<p><strong>Python code:</strong></p>
<pre>
def colWiseSum(W):
    ones=np.ones((W.shape[0],1))
    return ones.T.dot(W)
</pre>
<p><strong>R code:</strong></p>
<pre>
colWiseSum=function(W)
{
  ones=rep(1,nrow(W))
  t(W)%*%ones
}
</pre>
<h3> Row-wise matrix sum</h3>
<p>Similarly, the row-wise sum is <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-155c785d01b3f5d919d8b45539065d07_l3.png?resize=28%2C13" alt="mathbf{M}mathbf{1}" width="28" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="13">.</p>
<p><strong>Python code:</strong></p>
<pre>
def rowWiseSum(W):
    ones=np.ones((W.shape[1],1))
    return W.dot(ones)
</pre>
<p><strong>R code:</strong></p>
<pre>
rowWiseSum=function(W)
{
  ones=rep(1,ncol(W))
  W%*%ones
}
</pre>
<h3> Matrix sum</h3>
<p>The sum of all the elements of a matrix is the sum of the sum of its rows. Using previous expression, the sum of all the terms of <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-1a62e250fd3bdd137e35e78ef18c41f6_l3.png?resize=19%2C13" alt="mathbf{M}" width="19" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="13"> is <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-1755436353fa2d214147f98c45a53371_l3.png?resize=48%2C16" alt="mathbf{1}^Tmathbf{M}mathbf{1}" width="48" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="16">.</p>
<p><strong>Python code:</strong></p>
<pre>
def matSum(W):
    rhs_ones=np.ones((W.shape[1],1))
    lhs_ones=np.ones((W.shape[0],1))
    return lhs_ones.T.dot(W).dot(rhs_ones)
</pre>
<p><strong>R code:</strong></p>
<pre>
matSum=function(W)
{
  rhs_ones=rep(1,ncol(W))
  lhs_ones=rep(1,nrow(W))
  t(lhs_ones) %*% W%*% rhs_ones
}
</pre>
<h3>Similarity matrix (Gram matrix)</h3>
<p>Let’s say we have a set of words and for each of this words we want to find the most similar words from a dictionary. We assume that the words have been projected in space of dimension <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-5faad0904f612a3fa5b27faafb8dc903_l3.png?resize=10%2C12" alt="p" width="10" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12"> (using word2vect). Let <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-4c568e5c3b26b080d475ca653af29edd_l3.png?resize=15%2C14" alt="mathbf{X}" width="15" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="14"> (our set of words) and <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-9e00900a1ce3c771d11a9667e4a2b002_l3.png?resize=16%2C14" alt="mathbf{Y}" width="16" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="14"> (our dictionary) be two matrices resp. in <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-9cc4c6dcd3969da91fdc3f429ec62688_l3.png?resize=39%2C12" alt="mathbb{R}^{n_1,p}" width="39" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12"> and <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-1580fff38b29b362dbb6aceaf3848b04_l3.png?resize=39%2C12" alt="mathbb{R}^{n_2,p}" width="39" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12">. To compute the similarity of all the observations of <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-4c568e5c3b26b080d475ca653af29edd_l3.png?resize=15%2C14" alt="mathbf{X}" width="15" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="14"> and <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-9e00900a1ce3c771d11a9667e4a2b002_l3.png?resize=16%2C14" alt="mathbf{Y}" width="16" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="14"> we simply need to compute <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-27b8ead2d17deb9fee1dae7b7f6fa0c6_l3.png?resize=41%2C16" alt="mathbf{X} mathbf{Y}^T" width="41" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="16">.</p>
<p><strong>Python code:</strong></p>
<pre>
def gramMatrix(X,Y):
    return X.dot(Y.T)
</pre>
<p><strong>R code:</strong></p>
<pre>
gramMatrix=function(X,Y)
{
  X %*% t(Y)
}
</pre>
<h3>L2 distance</h3>
<p>We want to compute the pair-wise distance between two sets of vector. Let <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-4c568e5c3b26b080d475ca653af29edd_l3.png?resize=15%2C14" alt="mathbf{X}" width="15" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="14"> and <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-9e00900a1ce3c771d11a9667e4a2b002_l3.png?resize=16%2C14" alt="mathbf{Y}" width="16" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="14"> be two matrix in <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-9cc4c6dcd3969da91fdc3f429ec62688_l3.png?resize=39%2C12" alt="mathbb{R}^{n_1,p}" width="39" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12"> and <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-1580fff38b29b362dbb6aceaf3848b04_l3.png?resize=39%2C12" alt="mathbb{R}^{n_2,p}" width="39" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12">. For each vector of <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-996ff7036e644e89f8ac379fa58d0cf7_l3.png?resize=16%2C12" alt="X" width="16" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12"> we need to compute the distance with all the vectors of <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-42ae22abcaa05c2d6c2fdc3746446019_l3.png?resize=14%2C12" alt="Y" width="14" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12">.Hence, the output matrix should be of size <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-ae233d43ca66c51be13421bea93227cc_l3.png?resize=56%2C18" alt="(n_1,n_2)" width="56" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="18">.</p>
<p>If <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-7e5fbfa0bbbd9f3051cd156a0f1b5e31_l3.png?resize=10%2C8" alt="x" width="10" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="8"> and <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-38461fc041e953482219abf5d4cce1cb_l3.png?resize=9%2C12" alt="y" width="9" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12"> are two vectors, their <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-8076e0263b06f61c109bb1ccc449875c_l3.png?resize=19%2C15" alt="mathbf{L}^2" width="19" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="15"> distance is: </p>
<p><span>   </span><span>   </span><img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-e86be021d246f942e94fb980545ab82e_l3.png?resize=494%2C25" width="450" alt="[sum (x_i - y_i)^2 = sum x_i^2 + sum y_i^2 - 2 sum x_i y_i = x^T x + y^T y - 2 x^T y]" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="25"></p>
<p>. To compute all pairwise distance, some work is required on the last equality. All the matrices should be of size <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-22f8f43e9b8a414418aa63539dc8896a_l3.png?resize=87%2C18" alt="(n_1 times n_2,p)" width="87" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="18">, then the output vector of distance will be of size <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-ff787d07760647fbbeb083a3f05d0a83_l3.png?resize=87%2C18" alt="(n_1 times n_2,1)" width="87" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="18">, which can be reshaped into a vector of size <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-ae233d43ca66c51be13421bea93227cc_l3.png?resize=56%2C18" alt="(n_1,n_2)" width="56" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="18">.</p>
<p>The first two terms <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-983f9238d6c727156929a43dae0ace79_l3.png?resize=31%2C15" alt="x^T x" width="31" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="15"> and <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-1ad66e3357236df8577b111d796a8040_l3.png?resize=29%2C19" alt="y^T y" width="29" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="19"> need to be computed for each <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-7e5fbfa0bbbd9f3051cd156a0f1b5e31_l3.png?resize=10%2C8" alt="x" width="10" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="8"> and <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-38461fc041e953482219abf5d4cce1cb_l3.png?resize=9%2C12" alt="y" width="9" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12">. <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-af5d4f8f969e54eca37a884b2865575d_l3.png?resize=52%2C15" alt="mathbf{X} odot mathbf{X}" width="52" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="15"> is the element-wise multiplication of X with itself (its elements are <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-2ecb66724b0b3467aedcbb2d04726265_l3.png?resize=25%2C23" alt="x_{i,j}^2" width="25" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="23">). Hence, the i-th element of <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-47545f93883b0ca31f436d24e1a7d810_l3.png?resize=123%2C20" alt="x_{sq}=(mathbf{X} odot mathbf{X}) mathbf{1}" width="123" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="20"> is the squared sum of the coordinate of the i-th observation, <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-0f6051f7b740a2981d8197ec1c4c0fee_l3.png?resize=54%2C23" alt="sum_j x_{i,j}^2" width="54" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="23">. </p>
<p>However, this is a vector and its shape is <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-bb3a05ccaf386861bc993dc008ac4a4f_l3.png?resize=47%2C18" alt="(n_1,1)" width="47" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="18">. By replicating each of its elements <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-bf1e42c248eee22a0911c24c95fe28f0_l3.png?resize=18%2C11" alt="n_2" width="18" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="11"> times, we will get a matrix of size <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-ff787d07760647fbbeb083a3f05d0a83_l3.png?resize=87%2C18" alt="(n_1 times n_2,1)" width="87" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="18">. The replication can be done easily, if we consider the right matrix <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-edb823aac216b9ebb7bf4c7c5eb24909_l3.png?resize=30%2C18" alt="mathbf{Delta}_{n_2}^{n_1}" width="30" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="18">.</p>
<p>Let <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-1f763864e6567ca3d811ed4b71c7af7a_l3.png?resize=23%2C16" alt="mathbf{1}_{n_2}" width="23" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="16"> be a vector of one of size <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-bf1e42c248eee22a0911c24c95fe28f0_l3.png?resize=18%2C11" alt="n_2" width="18" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="11">. Let <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-edb823aac216b9ebb7bf4c7c5eb24909_l3.png?resize=30%2C18" alt="mathbf{Delta}_{n_2}^{n_1}" width="30" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="18"> be the matrix of size <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-2db3c4a0dd2aedcd5ad751b95922eefe_l3.png?resize=96%2C18" alt="(n_1times n_2, n_1)" width="96" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="18"> with <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-cebc0a013985f2695aeb53ded9e7afb1_l3.png?resize=17%2C12" alt="n_1" width="17" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12"> repetitions of <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-1f763864e6567ca3d811ed4b71c7af7a_l3.png?resize=23%2C16" alt="mathbf{1}_{n_2}" width="23" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="16"> on the “diagonal”:</p>
<p><span>   </span><span>   </span><img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-555e11c3115179d51f2353aa0838cc67_l3.png?resize=201%2C91" width="201" alt="[mathbf{Delta}_{n_2}^{n_1}=begin{matrix}                 mathbf{1}_{n_2} &amp;  0  &amp; ldots &amp; 0 \                 0  &amp; mathbf{1}_{n_2} &amp; ldots &amp; 0 \                 vdots &amp; vdots &amp; ddots &amp; vdots \                  0  &amp;   0       &amp;ldots &amp; mathbf{1}_{n_2} end{matrix}]" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="91"></p>
<p>Then, our final vector is <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-216eb2a24b936ffa1eb7ddb2cd64d523_l3.png?resize=54%2C18" alt="Delta_{n_2}^{n_1}x_{sq}" width="54" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="18"> (The same expression holds for <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-9e00900a1ce3c771d11a9667e4a2b002_l3.png?resize=16%2C14" alt="mathbf{Y}" width="16" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="14">). We denote <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-21f36758b04341c7980aa18b13ced720_l3.png?resize=12%2C14" alt="Phi" width="12" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="14"> a reshape operator (used to transform the previous vector in matrix). With previous part on similarity matrix, we get the following expression of the pairwise distance:</p>
<p><span>   </span><span>   </span><img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-cdf500fb0f5e1324f057f83770bbffd7_l3.png?resize=409%2C24" width="409" alt="[Phi(mathbf{Delta}_{n_2}^{n_1} ((mathbf{X} odot mathbf{X}) mathbf{1}_{p}))+Phi(mathbf{Delta}_{n_1}^{n_2} ((mathbf{Y} odot mathbf{Y}) mathbf{1}_{p})) -2 mathbf{X} mathbf{Y}^T]" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="24"></p>
<p>The previous expression can seem complex, but this will help us a lot to code the pairwise distance. We only have to do the translation from maths to Numpy or R.</p>
<p><strong>Python code:</strong></p>
<pre>
def L2dist(X,Y):
    n_1=X.shape[0]
    n_2=Y.shape[0]
    p=X.shape[1]
    ones=np.ones((p,1))
    x_sq=(X**2).dot(ones)[:,0]
    y_sq=(Y**2).dot(ones)[:,0]
    delta_n1_n2=np.repeat(np.eye(n_1),n_2,axis=0)
    delta_n2_n1=np.repeat(np.eye(n_2),n_1,axis=0)
    return np.reshape(delta_n1_n2.dot(x_sq),(n_1,n_2))+np.reshape(delta_n2_n1.dot(y_sq),(n_2,n_1)).T-2*gramMatrix(X,Y)
</pre>
<p><strong>R Code:</strong></p>
<pre>
L2dist=function(X,Y)
{
  n_1=dim(X)[1]
  n_2=dim(Y)[1]
  p=dim(X)[2]
  ones=rep(1,p)
  x_sq=X**2 %*% ones
  x_sq=t(matrix(diag(n_1) %x% rep(1, n_2) %*% x_sq, n_2,n_1))
  
  y_sq=Y**2 %*% ones
  y_sq=matrix(diag(n_2) %x% rep(1, n_1) %*% y_sq,n_1,n_2)
  
  x_sq+y_sq-2*gramMatrix(X,Y)
}
</pre>
<h3>L2 distance (improved)</h3>
<p>Actually the previous <em>L2dist</em> is not completely optimized <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-26a10d09b56bd54c851cafdfbcfed8f0_l3.png?resize=129%2C20" alt="mathbf{Delta}_{n_2}^{n_1} ((mathbf{X} odot mathbf{X}) mathbf{1}_{p})" width="129" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="20"> requires a lot of memory since <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-5d602d301016105a5ae8868fdf4b2565_l3.png?resize=29%2C18" alt="Delta}_{n_2}^{n_1}" width="29" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="18"> has <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-bf72cf53942d0108ebafb4c2b9b82657_l3.png?resize=58%2C21" alt="n_2 times n_1^2" width="58" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="21"> cells and is mostly empty. Using Numpy built-in function, we can circumvent this multiplication by directly repeating the vector (which reduces the memory footprints by a factor <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-510cf31ae7126a36d3b51aa03ef62e33_l3.png?resize=94%2C18" alt="max(n_1,n2)" width="94" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="18">):</p>
<p><strong>Python code:</strong></p>
<pre>
def L2dist_improved(X,Y):
    n_1=X.shape[0]
    n_2=Y.shape[0]
    p=X.shape[1]
    ones=np.ones((p,1))
    x_sq=(X**2).dot(ones)[:,0]
    y_sq=(Y**2).dot(ones)[:,0]
##Replace multiplication by a simple repeat
    X_rpt=np.repeat(x_sq,n_2).reshape((n_1,n_2))
    Y_rpt=np.repeat(y_sq,n_1).reshape((n_2,n_1)).T
    return X_rpt+Y_rpt-2*gramMatrix(X,Y)
</pre>
<p><strong>R code:</strong></p>
<pre>
L2dist_improved=function(X,Y)
{
  n_1=dim(X)[1]
  n_2=dim(Y)[1]
  p=dim(X)[2]
  ones=rep(1,p)
  x_sq=X**2 %*% ones
  x_sq=t(matrix(rep(x_sq,each=n_2),n_2,n_1))
  
  y_sq=Y**2 %*% ones
  y_sq=matrix(rep(y_sq,each=n_1),n_1,n_2)
  
  x_sq+y_sq-2*gramMatrix(X,Y)
}
</pre>
<h3>L2 distance – benchmark</h3>
<p>To show the interest of our previous work, let’s compare the computation speed of the vectorized L2 distance, the naive implementation and the scikit-learn implementation. The experiments are run on different size of dataset with 100 repetitions.</p>
<figure class="wp-caption"><img data-orig-file="https://i0.wp.com/enhancedatascience.com/wp-content/uploads/2018/05/computationTimeL2.png?fit=398%2C278" data-orig-size="398,278" data-medium-file="https://i0.wp.com/enhancedatascience.com/wp-content/uploads/2018/05/computationTimeL2.png?fit=300%2C210" data-attachment-id="2906" data-image-meta="{" data-image-description="" src="https://i0.wp.com/enhancedatascience.com/wp-content/uploads/2018/05/computationTimeL2.png?resize=398%2C278" data-permalink="http://enhancedatascience.com/2018/05/07/machine-learning-explained-vectorization-matrix-operations/computationtimel2/" data-large-file="https://i0.wp.com/enhancedatascience.com/wp-content/uploads/2018/05/computationTimeL2.png?fit=398%2C278" alt="compute time L2" data-image-title="computationTimeL2" sizes="(max-width: 398px) 100vw, 398px" width="398" data-comments-opened="1" class="wp-image-2906" data-recalc-dims="1" height="278"><figcaption class="wp-caption-text">Computation time of the L2 distance</figcaption></figure>
<p>
The vectorized implementation is 2 to 3 orders of magnitude faster than the naive implementation and on par with the scikit implementation.</p>
<p>The post <a rel="nofollow" href="http://enhancedatascience.com/2018/05/07/machine-learning-explained-vectorization-matrix-operations/" target="_blank">Machine Learning Explained: Vectorization and matrix operations</a> appeared first on <a rel="nofollow" href="http://enhancedatascience.com/" target="_blank">Enhance Data Science</a>.</p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://enhancedatascience.com/2018/05/07/machine-learning-explained-vectorization-matrix-operations/"> Enhance Data Science</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/Mq_tN5S30Co" height="1"></a></div>