Determination of latent dimensionality in international trade flow

Throughout this work vectors are denoted with lowercase bold letters, x, matrices are denoted with uppercase letters X, and tensors are denoted with uppercase script letters, $\mathcal{X}$. Mode-1 multiplication between an order three tensor and a matrix is defined $(\mathcal{X}\times_1 {Y})_{i,\,j,\,k} = \sum_{l = 1} {Y}_{i,\,l} \mathcal{X}_{l,\,j,\,k}$, and similarly for modes 2 and 3. The Frobenius norm of a matrix or tensor is the square root of the sum of the squares of the elements, or $|| {X}||_F = \sqrt{\sum_{i,\,j} {X}_{i,\,j}}$, $|| \mathcal{X}||_F = \sqrt{\sum_{i,\,j,\,k} \mathcal{X}_{i,\,j,\,k}}$. A single subscript on a matrix, or tensor indicates the slice of the object along the last index, so X_i indicates the i^th column of the matrix, and $\mathcal{X}_i$ indicates the i^th matrix along the third mode. Additionally, a superscript in parentheses is used to enumerate items in a set or ensemble, $\{ {X^{\,(1)}},{X^{\,(2)}},,{X^{\,(p)}}\}$.

Here we describe some relevant works that are working with a similar decomposition model. The first proposed model was the single domain Decomposition into Directional Components (DEDICOM) model [28]. DEDICOM is capable of analyzing asymmetric data in marketing research and it has both matrix and tensor versions without constraints on its factors:

$$\begin{align*} &\text{two-way DEDICOM} & {X} = {A} {R} {A}^\top & \text{where} {A} \in \mathbb{R}^{n \times r}, {R} \in \mathbb{R}^{r \times r} \\ &\text{three-way DEDICOM} & \mathcal{X}_k = {A} \mathcal{D}_k {R} \mathcal{D}_k {A}^\top &\text{for }k = 1,\dots,T \end{align*}$$

In the three-way DEDICOM model, $\mathcal{X}_k$ is the k^th frontal slice of the tensor $\mathcal{X}$, and $\mathcal{D}_k \in \mathbb{R}^{r \times r}$ is a diagonal matrix. These models describe a single domain in the sense that they require the row space to be the same as the column space. Since there is no constraint on the factors, these decompositions can be estimated similarly as SVD.

Later, an alternating algorithm, Alternating Simultaneous Approximation Least Square and Newton (ASALSAN) was proposed to perform a three-way DEDICOM [29] that can be applied to large and sparse data. Moreover, the non-negative version of three-way DEDICOM was also introduced using a multiplicative update algorithm. The model was then applied to analyze email network data and export/import data. However, the latent dimension was selected without being justified, and the meaning of the extracted factors was not analyzed in details.

A relaxed version of three-way DEDICOM is the RESCAL model [5]. The model decomposes a three-way tensor $\mathcal{X}$ as follows:

$$\begin{align*} \mathcal{X}_k = A \mathcal{R}_k A^\top \text{for }k = 1,\dots,m \end{align*}$$

In [5], the factors A and $\mathcal{R}_k$ minimized the l₂-regularized minimization problem, and can be estimated using a variation of ASALSAN:

$$\begin{align*} \min_{A, \mathcal{R}_k} \dfrac{1}{2}\left(\sum_k || \mathcal{X}_k - A \mathcal{R}_k A^\top ||^2_F\right) + \dfrac{1}{2}\lambda \left( ||A||^2_F + \sum_k || \mathcal{R}_k ||^2_F \right) \end{align*}$$

The model was then applied to different datasets, and the authors suggested to find latent dimension using cross-validation. Instead, they chose a rather large dimension (k = 20), and then used k-means clustering method to cluster the matrix A into 6 groups.

Following the initial work, reference [30] introduced updating schemes for different variants of the non-negative RESCAL mode, including least-squares with l₂ regularization, KL-divergence with l₁ regularization, and other.

The RESCAL model is a tensor decomposition that takes advantage of the inherent structure of relational properties such as those expected in a dynamic network. More specifically, RESCAL decomposes an order three tensor by finding a common low dimensional latent space for the first two modes such that

$$\begin{equation*} \mathcal{X} = \mathcal{R} \times_1 A \times_2 A \end{equation*}$$

where A is an n × r matrix containing the features, $\mathcal{R}$ is an r × r × T tensor capturing the mixing relations between the features, and ×_i is the mode-i product defined above. In practice, a decomposition is always approximated by solving the optimization problem

$$\begin{align} \mathrm{argmin}_{{A}, \mathcal{R}} || \mathcal{X} - \mathcal{R} \times_1 {A} \times_2 {A}||_F^2 \end{align} \tag{ 1 }$$

When applied to a dynamic networks, RESCAL is interpretable in the way that each column of A represents a group of objects or nodes, and $\mathcal{R}_t$ represents the relations among these groups at time t. An equivalent problem statement demonstrates that RESCAL simultaneously decomposes each slice of an order three tensor with a rank r factorization,

$$\begin{align} \mathrm{argmin}_{{A}, {R}_t} \sum_t || \mathcal{X}_t - {A} \mathcal{R}_t {A}^{\top} ||_F^2\;. \end{align} \tag{ 2 }$$

With this interpretation RESCAL attempts to find a common set of features that can be simultaneously used to represent both the row space and the column space of the matrices $\mathcal{X}_t$. To remove a scaling ambiguity, RESCAL also can constrain the columns of A to be unit norm, $|| A_i || = 1$ for 1 ≤ i ≤ r.

In reality, we can easily find experimental data that are non-negative, such as signal/imaging data, trading amount (in currency) between countries, DNA microarray, etc quite often for the non-negative data, and especially with a correct latent dimension, a non-negative decomposition provides more physically meaningful and easier-to-interpret features [14]. Nonnegative RESCAL provides the same advantage by decomposing the data with non-negative features, and non-negative mixing matrices with the optimization

$$\begin{equation} \begin{aligned} \mathrm{argmin}_{{A}, \mathcal{R}_t} & \quad \sum_t || \mathcal{X}_t - {A} \mathcal{R}_t {A}^{\top} ||_F^2\\ \mathrm{subject\ to} & \quad \begin{cases} \sum_j A_{ij} = 1, \mathrm{for\ } 1 \leq j \leq r \\ A \geq 0\\ \mathcal{R} \geq 0 \end{cases}.\end{aligned} \end{equation} \tag{ 3 }$$

Zoom In Zoom Out Reset image size

To solve the non-negative constraint minimization problem in equation (3), we use the multiplicative update scheme similar to the one used for DEDICOM model [29], and for non-negative RESCAL [30]. Starting from non-negative random initialization of A and R, the following update steps are performed until the relative error converges to a predefined tolerance:

$$\begin{array}{l} {({{\cal R}_t})_{ij}} = {({{\cal R}_t})_{ij}}\frac{{{{[{A^ \top }{{\cal X}_t}A]}_{ij}}}}{{{{[{A^ \top }A{{\cal R}_t}{A^ \top }A]}_{ij}} + \varepsilon }}\;\;{\rm{fort = 1}},...,{\rm{T}}\\ {{\rm{A}}_{{\rm{ij}}}}{\rm{ = }}{{\rm{A}}_{{\rm{ij}}}}\frac{{{{\left[ {\sum\limits_{{\rm{t = 1}}}^{\rm{T}} {{{\cal X}_{\rm{t}}}} {\rm{A}}{\cal R}_{\rm{t}}^ \top {\rm{ + }}{\cal X}_{\rm{t}}^ \top {\rm{A}}{{\cal R}_{\rm{t}}}} \right]}_{{\rm{ij}}}}}}{{{{\left[ {{\rm{A}}\left( {{{\cal R}_{\rm{t}}}{{\rm{A}}^ \top }{\rm{A}}{\cal R}_{\rm{t}}^ \top {\rm{ + }}{\cal R}_{\rm{t}}^ \top {{\rm{A}}^ \top }{\rm{A}}{{\cal R}_{\rm{t}}}} \right)} \right]}_{{\rm{ij}}}}{\rm{ + }}\varepsilon }}\;. \end{array} \tag{ 4 }$$

After the multiplicate update scheme converges, A and $\mathcal{R}$ are appropriately scaled such that columns of A have a sum equal to one. Notice that the decomposition can be scaled without affecting the reconstruction error.

To select an appropriate latent dimension, we adapt the clustering procedure that has found success in NMF [23]. For each explored latent dimension, k, our procedure applies three steps: (i) bootstrapping by resampling the data, (ii) decomposing the bootstrapped data, and then (iii) analyzing the cluster stability of the solutions. Figure 3 shows a diagram of the model selection scheme. To resample the data, we construct an ensemble of tensors $\{ {{\cal X}^{(1)}},{{\cal X}^{(2)}},....,{{\cal X}^{(P)}}\}$ sampled from $\mathcal{X}^{(p)}_{i,\,j,\,k} \sim {U}(1-\epsilon, 1+\epsilon) * \mathcal{X}_{i,\,j,\,k}$ for 1 ≤ p ≤ P where U(a, b) is a uniform distribution between a and b. Each resampling introduces some variability in the data which mitigates the possibility of overfitting.

Zoom In Zoom Out Reset image size

We use a custom clustering algorithm designed to exploit the stability of the NMF's solutions corresponding to the resampled data. The custom property of the clustering is that each cluster should contain precisely one feature vector from each A^(p). Our algorithm is based on k-means but iterates over each A^(p) to assign each vector to an appropriate centroid. The assignment is done by a greedy algorithm applied to the cosine similarity between the columns of A^(p) and the current centroids.

To analyze the quality of the clustering, we use silhouette statistics [27]. The silhouette of a single point has a value between −1 and 1 relating how close it is to other points in the same cluster and how far is this point to the closest of the other clusters. A high silhouette score indicates that the clusters are compact and well separated, while a low silhouette score indicates that the clusters are not well separated. We use both the mean silhouette score, as well as the minimum silhouette score of a group as cluster quality metrics.

The selection of the correct latent dimension is accomplished by considering both the silhouette statistics, which measures the stability of the solutions, and the relative error. We expect that with the correct latent dimension the solutions will cluster well with a small relative error. With a too small latent space, the relative error will be too large, and with too large latent space there will be latent features representing the noise in the system that are not stable and hence will not cluster well, resulting in a low silhouette score. We consider the correct latent dimension to be the largest dimension that generates low relative errors and high silhouette scores.

Here we test the performance of our protocol for model selection in determining the ground truth latent dimension of synthetic datasets with different levels of noise. First, a synthetic data tensor $\mathcal{X}$ with dimensions n × n × T and the latent dimension k is generated as follows:

• Elements of matrix $A_{n \times k}$ are randomly sampled from U[0, 10). The matrix A is only selected if rank(A) = k. Otherwise, it is regenerated.
• Elements of matrix A is sparsified by setting elements smaller than Thresh_A to zeros.
• The tensor $\mathcal{R}$ is also generated from U[0, 10).
• Tensor $\mathcal{R}$ is also sparsified by setting elements smaller than Thresh_R to zeros.
• $\mathcal{X} = \mathcal{R} \times_1 A \times_2 A + {\epsilon}$ where elements of are sampled from U[0, NoiseFactor).
• The noise level of the data tensor can be adjusted by the value of NoiseFactor.
• The noise level is computed as $\dfrac{|| \mathcal{X} - A \mathcal{R}A^\top||_F}{||A \mathcal{R}A^\top||_F}$.

The latent dimension determination procedure described in 2.3 is then applied to the simulated data. More specifically, the sample perturbation is 0.03 and the number of iterations P = 50.

In these scenarios, for a range of NoiseFactor ∈ {1, 10, 50, 100, 150, 200}, a tensor of dimension 10 × 10 × 100 is simulated with the ground truth dimension k_true = 4. The plots of relative reconstruction error, minimum and average silhouette scores are shown in figure 4. As the noise factor increases, the noise level also increases and gradually distorts the L-shape of the reconstruction error curves. However, the silhouette score curves are consistent up to k = 4, and after that, the silhouette score curves start diverging, demonstrating that the clusters are no longer compact and well separated.

Zoom In Zoom Out Reset image size

International trade has been shown to generate mutual benefit between countries by allowing them to focus on specialization and exchange their produced goods and services. Unsurprisingly, it has increasingly contributed to the Gross Domestic Product (GDP). In fact, according to the World Bank, in 2017, the world GDP is about $80 trillion. It has been shown that the economic growth of a country is strongly associated with its role in the world trading network. Therefore, understanding the trade pattern is an essential step before we can discuss the effects of international trade, or even suggest policy changes.

Here we show that by applying the non-negative RESCAL model, with our latent dimension determination method, we can decompose the international trade network into different groups of countries whose exports and imports are tightly linked together and their interactions over time are captured in the resulting interacting tensor. More interestingly, the groups' activities are also meaningful in the sense that they are consistent with stylized economic facts.

3.2.1. International trade flows data

3.2.2. Determine the latent dimension & optimal factors

To determine the latent dimension for the model, we resampled as described in section 2.3 with the uniform distribution U(0.9, 1.1) so that each value in our ensemble has the measured value ±10% error. We resampled from this distribution 50 times to construct our ensemble and calculated the silhouette statistic for a range of dimensions $k \in \{3,\dots,8\}$. Figure 5 shows the relative reconstruction errors, the average and minimum silhouette statistics leading us to conclude that that the true latent dimension is 5.

Zoom In Zoom Out Reset image size

To determine the optimal factors for the analysis, we first ran 100 iterations from random initialization on the original data with the selected dimension k = 5, and select the decomposition, which provided the lowest reconstruction error. For both purposes, the stopping criterion is the relative convergence rate = 1 × 10⁻⁸.

3.2.3. Interpretation of the decomposition

Here we show the interpretation and analysis of the optimal decomposition. First, since the columns of A are normalized to have a sum equal to one, the entry A_ij can be interpreted as how much the country i contributes into group j. Figure 6 shows that the latent factors, which here are the country contribution profile in each group, approximately correspond to five economic regions: Asia and Pacific (without China), Europe, NAFTA (Canada, Mexico, and the U.S.), U.S., and China. These geo-economic regions are essential participants in the international trade, verifying that the model selection procedure extracts meaningful latent factors.

Zoom In Zoom Out Reset image size

Next, we analyze the interacting tensor $\mathcal{R}$ to determine its economic meanings. We first look at the aggregate export level for each economic group over time by summing across the rows of each $\mathcal{R}_t$ without the diagonal elements. By doing this, we exclude the contribution of the groups themselves to their aggregate export activities. We then check how well these approximated export activities agree with four global and local economic recession periods, which are identified by the National Bureau of Economic Research (NBER).

As shown in figure 7, the approximated activities match well with the international trading trends in all considered periods. For example, during the Great Recession (12/2007−06/2009), in which international trade was dramatically decreasing, the activities of all groups dropped substantially. Secondly, during early 2000s Recession, which was partially caused by the dotcom bubble and September 11, the activities represent the fact that the trading trend of all groups, except Asia, were dropping. Thirdly, during the Asia Financial Crisis (1997−199), only the activities of the group Asia is going down, representing the fact that this Recession mostly affected Asian countries. Lastly, during the European Debt Crisis, which peaked in 2010−2012, while the economy of other regions was recovering from the Great Recession, only European countries struggled with their high government debt. This styled fact is replicated in the approximated activity of the Europe-group. Overall, the interacting tensor captured the connection between the export and economic health of different economic regions.

Zoom In Zoom Out Reset image size

Finally, we analyze the interaction tensor $R \in \mathbb{R}_{5 \times 5 \times 420}$ by summing the tensor R over time, which gives us the matrix $S \in \mathbb{R}_{5 \times 5}$ describing how strongly these groups interact. Rs is then normalized by its maximum value for better visualization.

$$\begin{equation*} \bar{S} = S/\mathrm{max}(S)\ \mathrm{where}\ S_{\,i,\,j} = \sum_{t = 1}^{420} R_{i,\,j,\,t} \end{equation*}$$

We makes three observations from figure 8. First, the European group has the strongest interaction with itself, which makes sense since this is a group of advanced economies, and they have formed the European Union since 1995. Second, the strong two-way connection between NAFTA and the United States and between China and Asia, are also matched with statistics. Third, by comparing the row and column for China and United States, we can see that China is a trade surplus (the amount of exports is greater than the number of imports), and the United States is a trade deficit. Overall, the interacting tensor did extract meaningful information from the data in the sense that they agree with international economics empirical facts, indicating that our model selection procedure can capture meaningful activities from the data.

Zoom In Zoom Out Reset image size