<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/calculating-a-fuzzy-kmeans-membership-matrix-with-r-and-rcpp/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://blog.revolutionanalytics.com/2017/08/kmeans-r-rcpp.html"> Revolutions</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<div>
<p><em>by <a rel="nofollow" href="https://www.linkedin.com/in/b%C5%82a%C5%BCej-moska-9a316113a/" target="_blank">Błażej Moska</a>, computer science student and data science intern</em> </p>
<p>Suppose that we have performed clustering K-means clustering in R and are satisfied with our results, but later we realize that it would also be useful to have a membership matrix. Of course it would be easier to repeat clustering using one of the fuzzy kmeans functions available in R (like <a rel="nofollow" href="https://www.rdocumentation.org/packages/cluster/versions/2.0.6/topics/fanny?tap_a=5644-dce66f&tap_s=10907-287229" target="_blank"><code>fanny</code></a>, for example), but since it is slightly different implementation the results could also be different and for some reasons we don’t want them to be changed. Knowing the equation we can construct this matrix on our own, after using the <a rel="nofollow" href="https://www.rdocumentation.org/packages/stats/versions/3.4.1/topics/kmeans?tap_a=5644-dce66f&tap_s=10907-287229" target="_blank"><code>kmeans</code></a> function. The equation is defined as follows (source: <a rel="nofollow" href="https://en.wikipedia.org/wiki/Fuzzy_clustering" target="_blank">Wikipedia</a>):</p>
<p>$$w_{ij} = \frac{1}{ \sum_ {k=1}^c ( \frac{ \| x_{i} – c_{j} \| }{ \| x_{i} – c_{k} \| }) ^{ \frac{2}{m-1} } } $$</p>
<p>\(w_{ij}\) denotes to what extent the \(i\)th object belongs to the \(j\)th cluster. So the total number of rows of this matrix equals number of observation and total number of columns equals number of variables included in clustering. \(m\) is a parameter, typically set to \(m=2\). \(w_{ij}\) values range between 0 and 1 so they are easy and convenient to compare. In this example I will use \(m = 2\) so the Euclidean distance will be calculated.</p>
<p>To make computations faster I also used the <a rel="nofollow" href="http://adv-r.had.co.nz/Rcpp.html" target="_blank">Rcpp package</a>, then I compared speed of execution of function written in R with this written in C++.</p>
<p>In implementations <code>for</code> loops were used, although it is not a commonly used method in R (see <a rel="nofollow" href="http://www.noamross.net/blog/2014/4/16/vectorization-in-r--why.html" target="_blank">this blog post</a> for more information and alternatives), but in this case I find it more convenient.</p>
<h3>Rcpp (C++ version)</h3>
<div>
<pre><span>#include &lt;Rcpp.h></span>
<span>#include &lt;math.h></span>
<span>using</span> <span>namespace</span> Rcpp;

<span>// [[Rcpp::export]]</span>
NumericMatrix <span>fuzzyClustering</span>(NumericMatrix data, NumericMatrix centers, <span>int</span> m) {  
<span>/* data is a matrix with observations(rows) and variables, </span>
<span> centers is a matrix with cluster centers coordinates, </span>
<span> m is a parameter of equation, c is a number of clusters</span>
<span>*/</span>
  
  <span>int</span> c<span>=</span>centers.rows();
  <span>int</span> rows <span>=</span> data.rows();
  <span>int</span> cols <span>=</span> data.cols(); <span>/*number of columns equals number of variables, the same as is in centers matrix*/</span>
  <span>double</span> tempDist<span>=</span><span>0</span>;        <span>/*dist and tempDist are variables storing temporary euclidean distances */</span>
  <span>double</span> dist<span>=</span><span>0</span>;      
  <span>double</span> denominator<span>=</span><span>0</span>;    <span>//denominator of “main” equation</span>
  
  NumericMatrix result(rows,c);    <span>//declaration of matrix of results</span>
  
  <span>for</span>(<span>int</span> i<span>=</span><span>0</span>;i<span></span>rows;i<span>++</span>){
    <span>for</span>(<span>int</span> j<span>=</span><span>0</span>;j<span></span>c;j<span>++</span>){
      <span>for</span>(<span>int</span> k<span>=</span><span>0</span>;k<span></span>c;k<span>++</span>){
        <span>for</span>(<span>int</span> p<span>=</span><span>0</span>;p<span></span>cols;p<span>++</span>){
          tempDist <span>=</span> tempDist<span>+</span>pow(centers(j,p)<span>-</span>data(i,p),<span>2</span>);      
         <span>//in innermost loop an euclidean distance is calculated.</span>
          dist <span>=</span> dist <span>+</span> pow(centers(k,p)<span>-</span>data(i,p),<span>2</span>);              
<span>/*tempDist is nominator inside the sum operator in the equation, dist is the denominator inside the sum operator in the equation*/</span>
        }
        tempDist <span>=</span> sqrt(tempDist);
        dist <span>=</span> sqrt(dist);
        denominator <span>=</span> denominator<span>+</span>pow((tempDist<span>/</span>dist),(<span>2</span><span>/</span>(m<span>-</span><span>1</span>)));
        tempDist <span>=</span> <span>0</span>;
        dist <span>=</span> <span>0</span>;
      }
      result(i,j) <span>=</span> <span>1</span><span>/</span>denominator;    
<span>// nominator/denominator in the main equation</span>
      denominator <span>=</span> <span>0</span>;                           
    }
  }
  <span>return</span> result;
}
</pre>
</div>
<p>We can save this in a file with .cpp extension. To compile it from R we can write:</p>
<pre>sourceCpp(&quot;path_to_cpp_file&quot;)
</pre>
<p>If everything goes right our function <code>fuzzyClustering</code> will then be available from R.</p>
<h3>R version</h3>
<div>
<pre>fuzzyClustering<span>=</span><span>function</span>(data,centers,m){

  c <span>&lt;-</span> nrow(centers)
  rows <span>&lt;-</span> nrow(data)
  cols <span>&lt;-</span> ncol(data)
  result <span>&lt;-</span> matrix(<span>0</span>,nrow<span>=</span>rows,ncol<span>=</span>c)  <span>#defining membership matrix</span>
  denominator <span>&lt;-</span> <span>0</span>
  
  <span>for</span>(i <span>in</span> <span>1</span><span>:</span>rows){
    <span>for</span>(j <span>in</span> <span>1</span><span>:</span>c){
      tempDist <span>&lt;-</span> sqrt(sum((centers[j,]<span>-</span>data[i,])<span>^</span><span>2</span>))  <span>#euclidean distance, nominator inside a sum operator </span>
      <span>for</span>(k <span>in</span> <span>1</span><span>:</span>c){
        Dist <span>&lt;-</span> sqrt(sum((centers[k,]<span>-</span>data[i,])<span>^</span><span>2</span>))    <span>#euclidean distance, denominator inside a sum operator </span>
        denominator <span>&lt;-</span> denominator <span>+</span>((tempDist<span>/</span>Dist)<span>^</span>(<span>2</span><span>/</span>(m<span>-1</span>))) <span>#denominator of an equation </span>
      }    
      result[i,j] <span>&lt;-</span> <span>1</span><span>/</span>denominator    <span>#inserting value into membership matrix</span>
      denominator <span>&lt;-</span> <span>0</span>           
    }
  }
  <span>return</span>(result);
}
</pre>
</div>
<p>Result looks as follows. Columns are cluster numbers (in this case 10 clusters were created), rows are our objects (observations). Values were rounded to the third decimal place, so the sums of rows can be slightly different than 1:</p>
<pre>          1     2     3     4     5     6     7     8     9    10  
 [1,] 0.063 0.038 0.304 0.116 0.098 0.039 0.025 0.104 0.025 0.188  
 [2,] 0.109 0.028 0.116 0.221 0.229 0.080 0.035 0.116 0.017 0.051  
 [3,] 0.067 0.037 0.348 0.173 0.104 0.066 0.031 0.095 0.018 0.062  
 [4,] 0.016 0.015 0.811 0.049 0.022 0.017 0.009 0.023 0.007 0.031  
 [5,] 0.063 0.048 0.328 0.169 0.083 0.126 0.041 0.079 0.018 0.045  
 [6,] 0.069 0.039 0.266 0.226 0.102 0.111 0.037 0.084 0.017 0.048  
 [7,] 0.045 0.039 0.569 0.083 0.060 0.046 0.025 0.071 0.015 0.046  
 [8,] 0.070 0.052 0.399 0.091 0.093 0.054 0.034 0.125 0.022 0.062  
 [9,] 0.095 0.037 0.198 0.192 0.157 0.088 0.038 0.121 0.019 0.055 
[10,] 0.072 0.024 0.132 0.375 0.148 0.059 0.025 0.081 0.015 0.067
</pre>
<h3>Performance comparison</h3>
<p>Shown below is the output of <code>Sys.time</code> for the C++ and R versions, running against a simulated matrix with 30000 observations, 3 variables and 10 clusters. <br>The hardware I used was a low-cost notebook Asus R556L with Intel Core i3-5010 2.1 GHz processor and 8 GB DDR3 1600 MHz RAM memory.</p>
<p>C++ version:</p>
<pre>  user    system  elapsed
  0.32    0.00    0.33 
</pre>
<p>R version:</p>
<pre>  user    system  elapsed
  15.75    0.02     15.94 <br>
</pre>
<p>In this example, the function written in C++ executed about 50 times faster than the equivalent function written in pure R.</p>
</div>
	<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://blog.revolutionanalytics.com/2017/08/kmeans-r-rcpp.html"> Revolutions</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 width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/xlJXCDYeMFU" height="1"></a></div>