<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/understanding-the-candecompparafac-tensor-decomposition-aka-cp-with-r-code/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.alexejgossmann.com/tensor_decomposition_CP/"> Alexej's blog</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>A tensor is essentially a multi-dimensional array:</p>
<ul>
<li>a tensor of order one is a vector, which simply is a column of numbers,</li>
<li>a tensor of order two is a matrix, which is basically numbers arranged in a rectangle,</li>
<li>a tensor of order three looks like numbers arranged in rectangular box (or a cube, if all modes have the same dimension),</li>
<li>an <em>n</em>th order (or <em>n</em>-way) tensor looks like numbers arranged in… an $n$-hyperrectangle or something like that… you get the idea…</li>
</ul>
<p>In many applications, data naturally form an <em>n</em>-way tensor with <em>n > 2</em>, rather than a “tidy” table.</p>
<p>So, what’s a tensor decomposition?</p>
<p>Well, there are several types of tensor decomposition, but in this blog post I will introduce only the CANDECOMP/PARAFAC decomposition. Following <a target="_blank" rel="nofollow" href="http://epubs.siam.org/doi/abs/10.1137/07070111X?journalCode=siread&">Kolda &amp; Bader (2009)</a> I will refer to it as <em>CP decomposition</em>. But before spelling it out in mathematical terms, let’s start with a simple toy example using the R language.</p>
<h2>Rank-1 approximation to a 3-way tensor (toy example)</h2>
<p>Assume that we observe spatio-temporal<sup><a target="_blank" rel="nofollow" href="http://www.alexejgossmann.com/tensor_decomposition_CP/#fn:1">1</a></sup> data that, apart from random noise, correspond to the following three classes:</p>
<p><img alt="Three classes for rank-1 tensor approximation toy example" src="https://i0.wp.com/www.alexejgossmann.com/images/CP_tensor_decomposition/cases_for_rank-1_approx.png?w=456" data-recalc-dims="1"></p>
<p>As you can see from the above figures, the essential difference between the three cases is that they behave differently with respect to time. The spatial component is just a Gaussian curve, while the temporal component is piecewise constant with a sudden jump at time 50, which differs in magnitude and in direction between the three classes.</p>
<p>Now assume we have collected samples that correspond to the three classes above (with some added noise). <em>But we don’t know which sample falls in what class, how many classes there are, and how they differ</em> :confused:. In an unsupervised fashion we want to learn the different classes and their differences with respect to time and space.</p>
<p>We can arrange the samples in a 3-way tensor, <em>sample-by-space-by-time</em>. For simplicity, however, assume that the samples are already grouped according to their class within the tensor (but the algorithm doesn’t know that!), so that the resulting tensor looks like this:</p>
<p><img alt="3-way tensor" src="https://i1.wp.com/www.alexejgossmann.com/images/CP_tensor_decomposition/3-modal_tensor.png?w=456" data-recalc-dims="1"></p>
<p>In R the above tensor (let’s call it <code>X</code>) can be generated with the following lines of code:</p>
<div>
<pre><code><span>space_index</span><span> </span><span>&lt;-</span><span> </span><span>seq</span><span>(</span><span>-1</span><span>,</span><span> </span><span>1</span><span>,</span><span> </span><span>l</span><span> </span><span>=</span><span> </span><span>100</span><span>)</span><span>
</span><span>bell_curve</span><span> </span><span>&lt;-</span><span> </span><span>dnorm</span><span>(</span><span>space_index</span><span>,</span><span> </span><span>mean</span><span> </span><span>=</span><span> </span><span>0</span><span>,</span><span> </span><span>sd</span><span> </span><span>=</span><span> </span><span>0.5</span><span>)</span><span>
</span><span>case1</span><span> </span><span>&lt;-</span><span> </span><span>matrix</span><span>(</span><span>rep</span><span>(</span><span>bell_curve</span><span>,</span><span> </span><span>10</span><span>),</span><span> </span><span>100</span><span>,</span><span> </span><span>100</span><span>)</span><span>
</span><span>case2</span><span> </span><span>&lt;-</span><span> </span><span>matrix</span><span>(</span><span>rep</span><span>(</span><span>bell_curve</span><span>,</span><span> </span><span>10</span><span>),</span><span> </span><span>100</span><span>,</span><span> </span><span>100</span><span>)</span><span>
</span><span>case3</span><span> </span><span>&lt;-</span><span> </span><span>matrix</span><span>(</span><span>rep</span><span>(</span><span>bell_curve</span><span>,</span><span> </span><span>10</span><span>),</span><span> </span><span>100</span><span>,</span><span> </span><span>100</span><span>)</span><span>
</span><span>case2</span><span>[</span><span> </span><span>,</span><span> </span><span>51</span><span>:</span><span>100</span><span>]</span><span> </span><span>&lt;-</span><span> </span><span>case2</span><span>[</span><span> </span><span>,</span><span> </span><span>51</span><span>:</span><span>100</span><span>]</span><span> </span><span>+</span><span> </span><span>0.1</span><span>
</span><span>case3</span><span>[</span><span> </span><span>,</span><span> </span><span>51</span><span>:</span><span>100</span><span>]</span><span> </span><span>&lt;-</span><span> </span><span>case3</span><span>[</span><span> </span><span>,</span><span> </span><span>51</span><span>:</span><span>100</span><span>]</span><span> </span><span>-</span><span> </span><span>0.1</span><span>
</span><span>X</span><span> </span><span>&lt;-</span><span> </span><span>array</span><span>(</span><span>NA</span><span>,</span><span> </span><span>dim</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>90</span><span>,</span><span> </span><span>100</span><span>,</span><span> </span><span>100</span><span>))</span><span>
</span><span>for</span><span>(</span><span>i</span><span> </span><span>in</span><span> </span><span>1</span><span>:</span><span>30</span><span>)</span><span> </span><span>{</span><span>
</span><span>X</span><span>[</span><span>i</span><span>,</span><span> </span><span>,</span><span> </span><span>]</span><span> </span><span>&lt;-</span><span> </span><span>case1</span><span> </span><span>+</span><span> </span><span>matrix</span><span>(</span><span>rnorm</span><span>(</span><span>10000</span><span>,</span><span> </span><span>sd</span><span> </span><span>=</span><span> </span><span>0.1</span><span>),</span><span> </span><span>100</span><span>,</span><span> </span><span>100</span><span>)</span><span>
</span><span>X</span><span>[</span><span>i</span><span>+30</span><span>,</span><span> </span><span>,</span><span> </span><span>]</span><span> </span><span>&lt;-</span><span> </span><span>case2</span><span> </span><span>+</span><span> </span><span>matrix</span><span>(</span><span>rnorm</span><span>(</span><span>10000</span><span>,</span><span> </span><span>sd</span><span> </span><span>=</span><span> </span><span>0.1</span><span>),</span><span> </span><span>100</span><span>,</span><span> </span><span>100</span><span>)</span><span>
</span><span>X</span><span>[</span><span>i</span><span>+60</span><span>,</span><span> </span><span>,</span><span> </span><span>]</span><span> </span><span>&lt;-</span><span> </span><span>case3</span><span> </span><span>+</span><span> </span><span>matrix</span><span>(</span><span>rnorm</span><span>(</span><span>10000</span><span>,</span><span> </span><span>sd</span><span> </span><span>=</span><span> </span><span>0.1</span><span>),</span><span> </span><span>100</span><span>,</span><span> </span><span>100</span><span>)</span><span>
</span><span>}</span><span>
</span></code></pre>
</div>
<p>Using the excellent R package <a target="_blank" rel="nofollow" href="https://cran.r-project.org/package=rTensor"><code>rTensor</code></a> we obtain the CP decomposition with one component per mode of the tensor:</p>
<div>
<pre><code><span>library</span><span>(</span><span>rTensor</span><span>)</span><span>
</span><span>cp_decomp</span><span> </span><span>&lt;-</span><span> </span><span>cp</span><span>(</span><span>as.tensor</span><span>(</span><span>X</span><span>),</span><span> </span><span>num_components</span><span> </span><span>=</span><span> </span><span>1</span><span>)</span><span>
</span><span>str</span><span>(</span><span>cp_decomp</span><span>$</span><span>U</span><span>)</span><span>
</span><span># List of 3
# $ : num [1:90, 1] 0.0111 0.0111 0.0111 0.0111 0.0112 ...
# $ : num [1:100, 1] -0.00233 -0.00251 -0.00271 -0.00292 -0.00314 ...
# $ : num [1:100, 1] -0.00996 -0.00994 -0.00996 -0.00993 -0.00997 ...
# NULL
</span></code></pre>
</div>
<p>Visualizing the three components, we get the following figures:</p>
<p><img alt="components of a rank-1 approximation to a 3-way tensor" src="https://i1.wp.com/www.alexejgossmann.com/images/CP_tensor_decomposition/rank-1_approx.png?w=456" data-recalc-dims="1"></p>
<p>We can clearly see that</p>
<ol>
<li>The sample-specific component has correctly separated the samples into three groups (in fact, <code>1:30</code> correspond to case 1, <code>31:60</code> correspond to case 2, <code>61:90</code> correspond to case 3).</li>
<li>The spatial component is bell-shaped, just as the input data with respect to the spatial dimension.</li>
<li>The temporal component has correctly picked up a change at time 50, which is where the three classes differ.</li>
</ol>
<p>So, what did just happen? :astonished:</p>
<h2>CP decomposition (quick summary of the math behind it)</h2>
<p>The CP decomposition factorizes a tensor into a sum of outer products of vectors. For example, for a 3-way tensor $X$, the CP decomposition can be written as</p>
<p></p>
<p>where $R>0$ and $u\subscript{r}$, $v\subscript{r}$, $w\subscript{r}$ are vectors of appropriate dimensions, and where the notation “$\circ$” denotes the outer product for tensors, i.e.,</p>
<p></p>
<p>In case of the previous toy example we have that $R = 1$. Thus, the corresponding CP decomposition has the following form:</p>
<p><img alt="Rank-1 CP decomposition cartoon" src="https://i1.wp.com/www.alexejgossmann.com/images/CP_tensor_decomposition/rank-1_decomposition_cartoon.png?w=456" data-recalc-dims="1"></p>
<p>I hope that explains, why the components <code>u</code>, <code>v</code>, and <code>w</code> in the toy example look the way they do! :bowtie:</p>
<p>Now, how do you solve for the components $u\subscript{r}$, $v\subscript{r}$, $w\subscript{r}$ ($r = 1, 2, \ldots, R$)? You need to solve the following optimization problem:</p>
<p></p>
<p>where $\lVert \cdot \rVert$ is the Frobenius norm.</p>
<p>The simplest way to do it is via an alternating least squares approach, where we would regard certain components as fixed while solving for others, and then iterate while alternating the components regarded as fixed… For much more rigour and detail see <a target="_blank" rel="nofollow" href="http://epubs.siam.org/doi/abs/10.1137/07070111X?journalCode=siread&">Kolda &amp; Bader (2009) <em>Tensor Decompositions and Applications</em></a>.</p>
<h2>A higher-rank approximation via CP (toy example cont.)</h2>
<p>Since in the previous toy example, there are no differentiating features between the three classes, apart from a jump in the temporal component, it makes perfect sense to set $R = 1$ in CP. In order to try out CP with more than one component per mode, I generated data with a more complex structure, and with further differentiation between the three groups with respect to their temporal and spatial makeup:</p>
<ol>
<li>The three classes still have a bell-shaped component in the spatial mode, but now each class has a different mean.</li>
<li>In the temporal mode the data is shaped like a sine wave, with different scaling per class.</li>
<li>As before, there is a sudden jump in the temporal mode at time 50.</li>
</ol>
<p>The three classes in the resulting dataset have the following means.</p>
<p><img alt="Three classes for rank-3 tensor approximation toy example" src="https://i1.wp.com/www.alexejgossmann.com/images/CP_tensor_decomposition/cases_for_rank-3_approx.png?w=456" data-recalc-dims="1"></p>
<p>As before, we generate a tensor <code>X</code> of dimensions 90 × 100 × 100, with 30 samples per class obscured with random noise.</p>
<p></p>
<p>We use a CP decomposition in order to obtain a rank-3 approximation to that tensor:</p>
<div>
<pre><code><span>cp_decomp</span><span> </span><span>&lt;-</span><span> </span><span>cp</span><span>(</span><span>as.tensor</span><span>(</span><span>X</span><span>),</span><span> </span><span>num_components</span><span> </span><span>=</span><span> </span><span>3</span><span>,</span><span> </span><span>max_iter</span><span> </span><span>=</span><span> </span><span>100</span><span>)</span><span>
</span></code></pre>
</div>
<p>Here, we increase <code>max_iter</code> to 100, in order to ensure convergence, as can be checked with the <code>conv</code> attribute:</p>
<div>
<pre><code><span>cp_decomp</span><span>$</span><span>conv</span><span>
</span><span># [1] TRUE
</span></code></pre>
</div>
<p>Since we set <code>num_components = 3</code>, the solution now has three components per mode, organized in a three-column matrix for each mode:</p>
<div>
<pre><code><span>str</span><span>(</span><span>cp_decomp</span><span>$</span><span>U</span><span>)</span><span>
</span><span># List of 3
# $ : num [1:90, 1:3] 0.00131 0.00137 0.00141 0.0014 0.00135 ...
# $ : num [1:100, 1:3] 0.000926 0.001345 0.001799 0.002228 0.002755 ...
# $ : num [1:100, 1:3] 0.00969 0.0097 0.00974 0.0097 0.00971 ...
# NULL
</span></code></pre>
</div>
<p>And we can even check the percentage of the Frobenius norm of $X$ explained by the rank-3 approximation $\widehat{X}$:</p>
<div>
<pre><code><span>cp_decomp</span><span>$</span><span>norm_percent</span><span>
</span><span># [1] 83.13865
</span></code></pre>
</div>
<p>83% isn’t too bad!</p>
<p>:grinning: Let’s look at a visualization of the obtained components!</p>
<p><img alt="components of a rank-3 approximation to a 3-way tensor" src="https://i2.wp.com/www.alexejgossmann.com/images/CP_tensor_decomposition/rank-3_approx.png?w=456" data-recalc-dims="1"></p>
<p>Indeed, we observe that,</p>
<ul>
<li>the sample-specific components $u\subscript{r}$ clearly distinguish between the three groups of samples (<code>1:30</code>, <code>31:60</code>, and <code>61:90</code>),</li>
<li>the spatial components $v\subscript{r}$ clearly picks up the Gaussian shapes with the three different means (at -0.5, 0, and 0.5),</li>
<li>the temporal components $w\subscript{r}$ clearly show the sudden jump at time 50, as well as a sine wave.</li>
</ul>
<p>That’s it for now. In the next couple of weeks I am planning to write a couple blog posts on other types of tensor decompositions and tensor regression methods, as I am learning about them.</p>
<hr>
<div>
<ol>
<li>
<p>The two data modes can correspond to many types of measurements, other than space and time. Here, I use space and time for example purposes only because those are very familiar concepts. I am neither suggesting that specifically spatio-temporal data should be analyzed in this way, nor that tensor decomposition is generally a good approach for spatio-temporal data (I actually have no idea). <a target="_blank" rel="nofollow" href="http://www.alexejgossmann.com/tensor_decomposition_CP/#fnref:1">↩</a></p>
</li>
</ol>
</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://www.alexejgossmann.com/tensor_decomposition_CP/"> Alexej's blog</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 title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/qFOHP6qbIn4"></a></div>