<p>In many applications, data naturally form an <em>n</em>-way tensor with <em>n > 2</em>, rather than a “tidy” table.
As mentioned in the beginning of <a href="http://www.alexejgossmann.com/tensor_decomposition_CP/">my last blog post</a>, 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 <em>n</em>-hyperrectangle… you get the idea…</li>
</ul><p>In this post I introduce the Tucker decomposition (<a href="https://link.springer.com/article/10.1007%2FBF02289464">Tucker (1966) “Some mathematical notes on three-mode factor analysis”</a>). The Tucker decomposition family includes methods such as</p>
<ol><li>the <em>higher-order SVD</em>, or HOSVD, which is a generalization of the matrix SVD to tensors (<a href="http://www.sandia.gov/~tgkolda/tdw2004/ldl-94-31.pdf">De Lathauwer, De Moor, and Vanderwalle (2000) “A multilinear singular value decomposition”</a>),</li>
<li>the <em>higher order orthogonal iteration</em>, or HOOI, which delivers the best approximation to a given tensor by another tensor with prescribed mode-1 rank, mode-2 rank, etc. (<a href="http://epubs.siam.org/doi/abs/10.1137/S0895479898346995?journalCode=sjmael">De Lathauwer, De Moor, and Vanderwalle (2000) “On the Best Rank-1 and Rank-(R1,R2,…,RN) Approximation of Higher-Order Tensors”</a>).</li>
</ol><p>I introduce both approaches, and in order to demonstrate the usefulness of these concepts, I present a simple data compression example using <a href="http://data.worldbank.org/data-catalog/world-development-indicators">The World Bank’s World Development Indicators dataset</a> (though I use the <a href="https://www.kaggle.com/worldbank/world-development-indicators">version available on Kaggle</a>).</p>
<p>However, before we can get started with the decompositions, we need to look at and understand the <em>k</em>-mode tensor product.</p>
<p>Throughout this post, I will also introduce the R functions from the package <a href="https://cran.r-project.org/package=rTensor"><code>rTensor</code></a>, which can be used to perform all of the presented computations.</p>
<h2>Tensor times matrix: the <em>k</em>-mode product</h2>
<p>The $k$-mode product of a tensor $X \in \mathbb{R}^{I\subscript{1} \times I\subscript{2} \times \ldots \times I\subscript{N}}$ with a matrix $A \in \mathbb{R}^{J \times I\subscript{k}}$ is written as</p>
<p>The resulting tensor $A$ is of size $I\subscript{1} \times \ldots \times I\subscript{k-1} \times J \times I\subscript{k+1} \times \ldots \times I\subscript{N}$, and contains the elements</p>
<p>It can be hard, at first, to understand what that definition really means, or to visualize it in your mind. I find that it becomes easier once you realize that the <em>k</em>-mode product amounts to multiplying each mode-<em>k</em> fiber of $X$ by the matrix $A$.</p>
<p>We can demonstrate that in R:</p>
<div><pre><code><span>library</span><span>(</span><span>rTensor</span><span>)</span><span>
</span><span>tnsr</span><span> </span><span>&lt;-</span><span> </span><span>as.tensor</span><span>(</span><span>array</span><span>(</span><span>1</span><span>:</span><span>12</span><span>,</span><span> </span><span>dim</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>2</span><span>,</span><span> </span><span>2</span><span>,</span><span> </span><span>3</span><span>)))</span><span>
</span><span>mat</span><span> </span><span>&lt;-</span><span> </span><span>matrix</span><span>(</span><span>1</span><span>:</span><span>6</span><span>,</span><span> </span><span>3</span><span>,</span><span> </span><span>2</span><span>)</span><span>
</span><span># 1-mode product performed via the function ttm in rTensor
</span><span>tnsr_times_mat</span><span> </span><span>&lt;-</span><span> </span><span>ttm</span><span>(</span><span>tnsr</span><span> </span><span>=</span><span> </span><span>tnsr</span><span>,</span><span> </span><span>mat</span><span> </span><span>=</span><span> </span><span>mat</span><span>,</span><span> </span><span>m</span><span> </span><span>=</span><span> </span><span>1</span><span>)</span><span>
</span></code></pre>
</div>
<p>Now, for example, the first slice of <code>tnsr_times_mat</code> is the same as the matrix product of <code>mat</code> with the first slice of <code>tnsr</code>:</p>
<div><pre><code><span>tnsr_times_mat</span><span>@</span><span>data</span><span>[</span><span> </span><span>,</span><span> </span><span>,</span><span> </span><span>1</span><span>]</span><span>
</span><span># [,1] [,2]
# [1,] 9 19
# [2,] 12 26
# [3,] 15 33
</span><span>mat</span><span> </span><span>%*%</span><span> </span><span>as.matrix</span><span>(</span><span>tnsr</span><span>@</span><span>data</span><span>[</span><span> </span><span>,</span><span> </span><span>,</span><span> </span><span>1</span><span>])</span><span>
</span><span># [,1] [,2]
# [1,] 9 19
# [2,] 12 26
# [3,] 15 33
</span></code></pre>
</div>
<p>You might want to play around some more with the function <code>ttm</code> in R to get a better understanding of the <em>k</em>-mode product.</p>
<p>A few important facts about the <em>k</em>-mode product:</p>
<ul><li>$X \times\subscript{m} A \times\subscript{n} B = X \times\subscript{n} B \times\subscript{m} A$ if $n \neq m$,</li>
<li>but $X \times\subscript{n} A \times\subscript{n} B = X \times\subscript{n} (BA)$ (in general $\neq X \times\subscript{n} B \times\subscript{n} A$).</li>
</ul><h2>Tucker decomposition</h2>
<p>The Tucker decomposition (<a href="https://link.springer.com/article/10.1007%2FBF02289464">Tucker (1966)</a>) decomposes a tensor into a core tensor multiplied by a matrix along each mode (i.e., transformed via a $k$-mode product for every $k = 1, 2, \ldots, N$):</p>
<p>Note that $G$ might be much smaller than the original tensor $X$ if we accept an approximation instead of an exact equality.</p>
<p>In case of three-way tensors, we can hold on to the following mental image:</p>
<p><img alt="Figure showing a Tucker decomposition of a 3-way tensor" src="http://www.alexejgossmann.com/images/Tucker_decomposition/3-modal_Tucker_decomposition.jpg"></p>
<p>It is interesting to note that the CP decomposition, that I introduced in <a href="http://www.alexejgossmann.com/tensor_decomposition_CP/">a previous blog post</a>, is a special case of the Tucker decomposition, where the core tensor $G$ is constrained to be superdiagonal.</p>
<h3>Higher-order SVD (HOSVD)</h3>
<p>So, how do you compute the Tucker decomposition?</p>
<p>Many algorithms rely on the following fundamental equivalence:</p>
<p>The above equation uses some notation that was not introduced yet:</p>
<ul><li>$\otimes$ denotes the <a href="https://en.wikipedia.org/wiki/Kronecker_product">Kronecker product</a>.</li>
<li>
<p>$X\subscript{(k)}$ is the mode-$k$ unfolding (or mode-$k$ matricization) of the tensor $X$. The mode-$k$ unfolding arranges the mode-$k$ fibers (a <em>fiber</em> is a generalization of <em>column</em> to tensors) of $X$ as columns into a matrix. The concept may be easiest to understand by looking at an example. The following R code shows a 3-way tensor and all three of its mode-$k$ unfoldings (using the <code>k_unfold</code> function from the <code>rTensor</code> package):</p>
<div><pre><code><span>tnsr</span><span> </span><span>&lt;-</span><span> </span><span>as.tensor</span><span>(</span><span>array</span><span>(</span><span>1</span><span>:</span><span>12</span><span>,</span><span> </span><span>dim</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>2</span><span>,</span><span> </span><span>3</span><span>,</span><span> </span><span>2</span><span>)))</span><span>
</span><span>tnsr</span><span>@</span><span>data</span><span>
</span><span># , , 1
#
# [,1] [,2] [,3]
# [1,] 1 3 5
# [2,] 2 4 6
#
# , , 2
#
# [,1] [,2] [,3]
# [1,] 7 9 11
# [2,] 8 10 12
</span><span>
</span><span># mode-1 unfolding:
</span><span>k_unfold</span><span>(</span><span>tnsr</span><span>,</span><span> </span><span>1</span><span>)</span><span>@</span><span>data</span><span>
</span><span># [,1] [,2] [,3] [,4] [,5] [,6]
# [1,] 1 3 5 7 9 11
# [2,] 2 4 6 8 10 12
</span><span>
</span><span># mode-2 unfolding:
</span><span>k_unfold</span><span>(</span><span>tnsr</span><span>,</span><span> </span><span>2</span><span>)</span><span>@</span><span>data</span><span>
</span><span># [,1] [,2] [,3] [,4]
# [1,] 1 2 7 8
# [2,] 3 4 9 10
# [3,] 5 6 11 12
</span><span>
</span><span># mode-3 unfolding:
</span><span>k_unfold</span><span>(</span><span>tnsr</span><span>,</span><span> </span><span>3</span><span>)</span><span>@</span><span>data</span><span>
</span><span># [,1] [,2] [,3] [,4] [,5] [,6]
# [1,] 1 2 3 4 5 6
# [2,] 7 8 9 10 11 12
</span></code></pre>
</div>
</li>
</ul><p>A straightforward approach to solve the Tucker decomposition would be to solve each mode-$k$ matricized form of the Tucker decomposition (shown in the equivalence above) for $A^{(k)}$. This approach is known as <em>higher order SVD</em>, or HOSVD. It can be regarded as a generalization of the matrix SVD, because the matrices $A^{(k)}$ are orthogonal, while the tensor $G$ is “ordered” and “all-orthogonal” (see <a href="http://www.sandia.gov/~tgkolda/tdw2004/ldl-94-31.pdf">De Lathauwer et. al. (2000)</a> for detail). The resulting algorithm is shown below.</p>
<p>In R we can perform HOSVD using the function <code>hosvd</code> from <code>rTensor</code>:</p>
<div><pre><code><span>tnsr</span><span> </span><span>&lt;-</span><span> </span><span>rand_tensor</span><span>(</span><span>modes</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>30</span><span>,</span><span> </span><span>40</span><span>,</span><span> </span><span>50</span><span>))</span><span>
</span><span>hosv_decomp</span><span> </span><span>&lt;-</span><span> </span><span>hosvd</span><span>(</span><span>tnsr</span><span>)</span><span>
</span></code></pre>
</div>
<p>Now <code>hosv_decomp$Z</code> is our matrix $G$, and <code>hosv_decomp$U</code> is a list containing all the matrices $A^{(k)}$. We can use the function <code>ttl</code>, which performs multiple <em>k</em>-mode products on multiple modes successively given a tensor and a list of matrices, to check that up to numerical error the equation
$X = G \times\subscript{1} A^{(1)} \times\subscript{2} A^{(2)} \times\subscript{3} \ldots \times\subscript{N} A^{(N)}$
is satisfied:</p>
<div><pre><code><span>HOSVD_prod</span><span> </span><span>&lt;-</span><span> </span><span>ttl</span><span>(</span><span>hosv_decomp</span><span>$</span><span>Z</span><span>,</span><span> </span><span>hosv_decomp</span><span>$</span><span>U</span><span>,</span><span> </span><span>1</span><span>:</span><span>3</span><span>)</span><span>
</span><span>error</span><span> </span><span>&lt;-</span><span> </span><span>tnsr</span><span> </span><span>-</span><span> </span><span>HOSVD_prod</span><span>
</span><span>table</span><span>(</span><span>abs</span><span>(</span><span>error</span><span>@</span><span>data</span><span>)</span><span> </span><span>&lt;</span><span> </span><span>1e-12</span><span>)</span><span>
</span><span>#
# TRUE
# 60000
</span></code></pre>
</div>
<h3>Higher order orthogonal iteration (HOOI)</h3>
<p>Note that we can also use HOSVD to compress $X$ by truncating the matrices $A^{(k)}$. The truncated HOSVD, however, is known to not give the best fit, as measured by the norm of the difference</p>
<p>The <em>higher order orthogonal iteration</em>, or HOOI, algorithm finds the optimal approximation $\widehat{X}$ (with respect to the Frobenius norm loss) by, essentially, iterating the alternating truncation and SVD until convergence. If we truncate $A^{(k)}$ to have $r\subscript{k}$ columns, then the HOOI solution can be obtained by the following algorithm.</p>
<h3>Application of HOOI to data compression</h3>
<p>The example considered below is somewhat silly, given that the tensor I’m compressing isn’t very big, and thus there isn’t much of a point in compressing it. However, I think that the example still shows off very well how the algorithm can be very useful when the data size is much bigger (or the available storage much smaller).</p>
<p>I have downloaded from Kaggle the World Development Indicators <a href="https://www.kaggle.com/worldbank/world-development-indicators">dataset</a>, originally collected and published by the The World Bank (the original dataset is available <a href="http://data.worldbank.org/data-catalog/world-development-indicators">here</a>).</p>
<p>The data can be arranged into a three-way tensor with the three modes corresponding to <em>country</em> (<a href="https://www.kaggle.com/benhamner/d/worldbank/world-development-indicators/countries-in-the-wdi-data">list of available countries</a>), <em>indicator</em> (<a href="https://www.kaggle.com/benhamner/d/worldbank/world-development-indicators/indicators-in-data">list of available indicators</a>), and <em>year</em> (1960-2014). Since I didn’t have any time to deal with NA values in any creative way, I have kept only three indicators in the dataset. And I have replaced the remaining NAs with a country-wise average value for each particular indicator. Also, I have forgotten to normalize the data :disappointed:. The preprocessing resulted in a tensor of size 247-countries-by-3-indicators-by-55-years, that looks sort of like this:</p>
<p><img alt="Figure: 3-way tensor obtained from the WDI data" src="http://www.alexejgossmann.com/images/Tucker_decomposition/3-way_tensor_small.png"></p>
<p>In particular, large stretches of the data within a given country tend to be nearly constant, or nearly piece-wise constant.</p>
<p>We use the function <code>tucker</code> from <code>rTensor</code> to obtain a Tucker decomposition via HOOI, where we set the ranks to the value 3 at each mode.</p>
<div><pre><code><span>dim</span><span>(</span><span>wdi_tnsr</span><span>)</span><span>
</span><span># [1] 247 3 55
</span><span>
</span><span>tucker_decomp</span><span> </span><span>&lt;-</span><span> </span><span>tucker</span><span>(</span><span>wdi_tnsr</span><span>,</span><span> </span><span>ranks</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>3</span><span>,</span><span> </span><span>3</span><span>,</span><span> </span><span>3</span><span>))</span><span>
</span><span>str</span><span>(</span><span>tucker_decomp</span><span>)</span><span>
</span><span># List of 7
# $ Z :Formal class 'Tensor' [package &quot;rTensor&quot;] with 3 slots
# .. ..@ num_modes: int 3
# .. ..@ modes : int [1:3] 3 3 3
# .. ..@ data : num [1:3, 1:3, 1:3] -6.60e+10 -1.13e+05 6.24e+05 -7.76e+05 -1.93e+08 ...
# $ U :List of 3
# ..$ : num [1:247, 1:3] -0.02577 -0.00065 -0.01146 -0.19637 -0.17317 ...
# ..$ : num [1:3, 1:3] -1.00 -6.97e-10 -2.08e-02 2.08e-02 -4.70e-08 ...
# ..$ : num [1:55, 1:3] -0.0762 -0.0772 -0.0785 -0.0802 -0.082 ...
# $ conv : logi TRUE
# $ est :Formal class 'Tensor' [package &quot;rTensor&quot;] with 3 slots
# .. ..@ num_modes: int 3
# .. ..@ modes : int [1:3] 247 3 55
# .. ..@ data : num [1:247, 1:3, 1:55] 9.83e+07 4.44e+06 8.81e+07 1.05e+09 8.97e+08 ...
# $ norm_percent: num 99.4
# $ fnorm_resid : num 3.9e+08
# $ all_resids : num [1:2] 3.9e+08 3.9e+08
# NULL
</span></code></pre>
</div>
<p>To see how well the tensor decomposition approximates the original tensor, we can look at the relative error</p>
<div><pre><code><span>wdi_tnsr_approx</span><span> </span><span>&lt;-</span><span> </span><span>ttl</span><span>(</span><span>tucker_decomp</span><span>$</span><span>Z</span><span>,</span><span> </span><span>tucker_decomp</span><span>$</span><span>U</span><span>,</span><span> </span><span>1</span><span>:</span><span>3</span><span>)</span><span>
</span><span>fnorm</span><span>(</span><span>wdi_tnsr</span><span> </span><span>-</span><span> </span><span>wdi_tnsr_approx</span><span>)</span><span> </span><span>/</span><span> </span><span>fnorm</span><span>(</span><span>wdi_tnsr</span><span>)</span><span>
</span><span># [1] 0.005908934
</span></code></pre>
</div>
<p>and at the percentage of the norm of the original tensor explained by the Tucker decomposition</p>
<div><pre><code><span>tucker_decomp</span><span>$</span><span>norm_percent</span><span>
</span><span># [1] 99.40911
</span></code></pre>
</div>
<p>We, observe that we indeed achieve a recovery with an accuracy of over 99%. For comparison, the original tensor contains <code>247 * 3 * 55 = 40755</code> entries, while the computed Tucker decomposition consists of only <code>127 * 3 + 3 * 3 + 55 * 3 + 3 * 3 * 3 = 582</code> numbers. That’s a reduction in size by a factor greater than 70.</p>
<p>Even though data compression does not make much sense for the size of the dataset considered here, it clearly shows potential to be very useful for purposes of data distribution and data storage, when the data size far exceeds the terabyte range.</p>