Legendre decomposition for tensors^*

Legendre decomposition always finds the best approximation of a given tensor $\newcommand{\ten}[1]{\mathcal{#1}} \ten{P}$. Our strategy is to choose an index set $B \subseteq [I_1] \times [I_2] \times \dots \times [I_N]$ as a decomposition basis, where we assume $(1, 1, \dots, 1) \not\in B$ for a technical reason, and approximate the normalized tensor $\newcommand{\ten}[1]{\mathcal{#1}} \ten{P}$ by a multiplicative combination of parameters associated with B.

First we introduce the relation '$\leqslant$' between index vectors $u = (\,j_1, \dots, j_N)$ and $v = (i_1, \dots, i_N)$ as $u \leqslant v$ if $j_1 \leqslant i_1$, $j_2 \leqslant i_2$, ..., $j_N \leqslant i_N$. It is clear that this relation gives a partial order (Gierz et al 2003); that is, $\leqslant$ satisfies the following three properties for all $u, v, w$: (1) $v \leqslant v$ (reflexivity), (2) $u \leqslant v$, $\newcommand{\R}{\mathbb{R}} v \leqslant u \Rightarrow u = v$ (antisymmetry), and (3) $u \leqslant v$, $\newcommand{\R}{\mathbb{R}} v \leqslant w \Rightarrow u \leqslant w$ (transitivity). Each tensor is treated as a discrete probability mass function with the sample space $\Omega \subseteq [I_1] \times \dots \times [I_N]$. While it is natural to set $\Omega = [I_1] \times \dots \times [I_N]$, our formulation allows us to use any subset $\Omega \subseteq [I_1] \times \dots \times [I_N]$. Hence, for example, we can remove unnecessary indices such as missing or zero entries of an input tensor $\newcommand{\ten}[1]{\mathcal{#1}} \ten{P}$ from $\Omega$. We define $\Omega^+ = \Omega \setminus \{(1, 1, \dots, 1)\}$.

We define a tensor $\newcommand{\R}{\mathbb{R}} \newcommand{\ten}[1]{\mathcal{#1}} \ten{Q} \in \R^{I_1 \times I_2 \times \dots \times I_N}_{\geqslant 0}$ as fully decomposable with $B \subseteq \Omega^+$ if each entry $q_{v} \in \Omega$ is represented in the form of

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\exp}{\mathop{{\rm exp}}} \newcommand{\re}{{\rm Re}} \displaystyle \label{eq:theta} q_v = \frac{1}{\exp(\psi(\theta))} \prod_{u \in {\downarrow}v} \exp(\theta_u),\quad {\downarrow}v = \{u \in B | u \leqslant v\}, \nonumber \end{align} \tag{ 1 }$$

using $|B|$ parameters $(\theta_v)_{v \in B}$ with $\newcommand{\R}{\mathbb{R}} \theta_{v} \in \R$ and the normalizer $\newcommand{\R}{\mathbb{R}} \psi(\theta) \in \R$, which is always uniquely determined from the parameters $(\theta_v)_{v \in B}$ as

$$\begin{align*} \newcommand{\e}{{\rm e}} \renewcommand{\exp}{\mathop{{\rm exp}}} \newcommand{\re}{{\rm Re}} \renewcommand{\log}{\mathop{{\rm log}}} \displaystyle \psi(\theta) = \log \sum_{v \in \Omega} \prod_{u \in {\downarrow}v} \exp(\theta_{u}). \nonumber \end{align*}$$

This normalization does not have any effect on the decomposition performance, but rather it is needed to formulate our decomposition as an information geometric projection, as shown in the next subsection. There are two extreme cases for a choice of a basis B: if $\newcommand{\e}{{\rm e}} B = \emptyset$, a fully decomposable $\newcommand{\ten}[1]{\mathcal{#1}} \ten{Q}$ is always uniform; that is, $q_{v} = 1 / |\Omega|$ for all $v \in \Omega$. In contrast, if $B = \Omega^+$, any input $\newcommand{\ten}[1]{\mathcal{#1}} \ten{P}$ itself becomes decomposable.

We now define Legendre decomposition as follows: given an input tensor $\newcommand{\R}{\mathbb{R}} \newcommand{\ten}[1]{\mathcal{#1}} \ten{P} \in \R^{I_1 \times I_2 \times \dots \times I_N}_{\geqslant 0}$, the sample space $\Omega \subseteq [I_1] \times [I_2] \times \dots \times [I_N]$, and a parameter basis $B \subseteq \Omega^+$, Legendre decomposition finds the fully decomposable tensor $\newcommand{\R}{\mathbb{R}} \newcommand{\ten}[1]{\mathcal{#1}} \ten{Q} \in \R^{I_1 \times I_2 \times \dots \times I_N}_{\geqslant 0}$ with a B that minimizes the KL divergence $\newcommand{\ten}[1]{\mathcal{#1}} \newcommand{\KL}{D_{{\rm KL}}} \newcommand{\re}{{\rm Re}} \renewcommand{\log}{\mathop{{\rm log}}} \KL(\ten{P}, \ten{Q}) = \sum_{v \in \Omega} p_v \log (\,p_v / q_v)$ (figure 1(a)). In the next subsection, we introduce an additional parameter $\newcommand{\e}{{\rm e}} (\eta_v)_{v \in B}$ and show that this decomposition is always possible via the dual parameters $\newcommand{\e}{{\rm e}} (\,(\theta_v)_{v \in B}, (\eta_v)_{v \in B}\,)$ with information geometric analysis. Since $\theta$ and $\newcommand{\e}{{\rm e}} \eta$ are connected via Legendre transformation, we call our method Legendre decomposition.

Zoom In Zoom Out Reset image size

Legendre decomposition for second-order tensors (that is, matrices) can be viewed as a low rank approximation not of an input matrix $\newcommand{\ten}[1]{\mathcal{#1}} \ten{P}$ but of its entry-wise logarithm $\newcommand{\ten}[1]{\mathcal{#1}} \newcommand{\re}{{\rm Re}} \renewcommand{\log}{\mathop{{\rm log}}} \log \ten{P}$. This is why the rank of $\newcommand{\ten}[1]{\mathcal{#1}} \newcommand{\re}{{\rm Re}} \renewcommand{\log}{\mathop{{\rm log}}} \log\ten{Q}$ with the fully decomposable matrix $\newcommand{\ten}[1]{\mathcal{#1}} \ten{Q}$ coincides with the parameter matrix $\newcommand{\R}{\mathbb{R}} \newcommand{\ten}[1]{\mathcal{#1}} \ten{T} \in \R^{I_1 \times I_2}$ such that $t_{v} = \theta_{v}$ if $v \in B$, $\newcommand{\re}{{\rm Re}} \renewcommand{\exp}{\mathop{{\rm exp}}} \newcommand{\e}{{\rm e}} t_{(1, 1)} = 1 / \exp(\psi(\theta))$, and $t_{v} = 0$ otherwise. Thus we fill zeros into entries in $\Omega^+ \setminus B$. Then we have $\newcommand{\re}{{\rm Re}} \renewcommand{\log}{\mathop{{\rm log}}} \log q_{v} = \sum_{u \in {\downarrow}v} t_{u}$, meaning that the rank of $\newcommand{\ten}[1]{\mathcal{#1}} \newcommand{\re}{{\rm Re}} \renewcommand{\log}{\mathop{{\rm log}}} \log\ten{Q}$ coincides with the rank of $\newcommand{\ten}[1]{\mathcal{#1}} \ten{T}$. Therefore if we use a decomposition basis B that includes only l rows (or columns), $\newcommand{\rank}{\mathop{{\rm rank}}} \newcommand{\ten}[1]{\mathcal{#1}} \newcommand{\re}{{\rm Re}} \renewcommand{\log}{\mathop{{\rm log}}} \rank(\log\ten{Q}) \leqslant l$ always holds.

We solve the Legendre decomposition by formulating it as a convex optimization problem. Let us assume that $B = \Omega^+ = \Omega \setminus \{(1, 1, \dots, 1)\}$, which means that any tensor is fully decomposable. Our definition in equation (1) can be re-written as

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\log}{\mathop{{\rm log}}} \newcommand{\re}{{\rm Re}} \displaystyle \label{eq:loglinear} \log q_v = \sum_{u \in \Omega^+} \zeta(u, v) \theta_u - \psi(\theta) = \sum_{u \in \Omega} \zeta(u, v) \theta_u,\quad \zeta(u, v) = \left\{\begin{array}{@{}ll@{}} 1 & {\rm if}~ u \leqslant v, \nonumber \\ 0 & {\rm otherwise} \end{array} \right. \nonumber \end{align} \tag{ 2 }$$

with $- \psi(\theta) = \theta_{(1, 1, \dots, 1)}$, and the sample space $\Omega$ is a poset (partially ordered set) with respect to the partial order '$\leqslant$' with the least element $\bot = (1, 1, \dots, 1)$. Therefore our model belongs to the log-linear model on posets introduced by Sugiyama et al (2016, 2017), which is an extension of the information geometric hierarchical log-linear model (Amari 2001, Nakahara and Amari 2002). Each entry $q_v$ and the parameters $(\theta_v)_{v \in \Omega^+}$ in equation (2) directly correspond to those in equation (8) in Sugiyama et al (2017). According to theorem 2 in Sugiyama et al (2017), if we introduce $\newcommand{\e}{{\rm e}} (\eta_v)_{v \in \Omega^+}$ such that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:eta} \eta_v = \sum_{u \in {\uparrow}v} q_u = \sum_{u \in \Omega} \zeta(v, u) q_u,\quad {\uparrow}v = \{u \in \Omega | u \geqslant v\}, \nonumber \end{align} \tag{ 3 }$$

for each $v \in \Omega^+$ (see figure 1(b)), the pair $\newcommand{\e}{{\rm e}} (\,(\theta_v)_{v \in \Omega^+}, (\eta_v)_{v \in \Omega^+}\,)$ is always a dual coordinate system of the set of normalized tensors $\newcommand{\Sc}{\mathcal{S}} \newcommand{\ten}[1]{\mathcal{#1}} \boldsymbol{\Sc} = \{\ten{P} \mid 0 < p_v < 1~{\rm and}~\sum_{v \in \Omega} p_v = 1\}$ with respect to the sample space $\Omega$, as they are connected via Legendre transformation. Hence $\newcommand{\Sc}{\mathcal{S}} \boldsymbol{\Sc}$ becomes a dually flat manifold (Amari 2009).

Here we formulate Legendre decomposition as a projection of a tensor onto a submanifold. Suppose that $B \subseteq \Omega^+$ and let $\newcommand{\Pf}{\boldsymbol{\mathcal{S}}} \boldsymbol{\Pf}_{B}$ be the submanifold such that

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\ten}[1]{\mathcal{#1}} \newcommand{\Pf}{\boldsymbol{\mathcal{S}}} \displaystyle \boldsymbol{\Pf}_{B} = \{\ten{Q} \in \boldsymbol{\Pf} | \theta_v = 0 ~{\rm for~all}~ v \in \Omega^+ \setminus B\}, \nonumber \end{align*}$$

which is the set of fully decomposable tensors with B and is an e-flat submanifold as it has constraints on the $\theta$ coordinate (Amari (2016), chapter 2.4). Furthermore, we introduce another submanifold $\newcommand{\Pf}{\boldsymbol{\mathcal{S}}} \newcommand{\ten}[1]{\mathcal{#1}} \boldsymbol{\Pf}_{\ten{P}}$ for a tensor $\newcommand{\Pf}{\boldsymbol{\mathcal{S}}} \newcommand{\ten}[1]{\mathcal{#1}} \ten{P} \in \boldsymbol{\Pf}$ and $A \subseteq \Omega^+$ such that

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\ten}[1]{\mathcal{#1}} \newcommand{\Pf}{\boldsymbol{\mathcal{S}}} \displaystyle \boldsymbol{\Pf}_{\ten{P}} = \{\ten{Q} \in \boldsymbol{\Pf} | \, \eta_v = \hat{\eta}_v ~ {\rm for~all}~ v \in A\}, \nonumber \end{align*}$$

where $\newcommand{\e}{{\rm e}} \hat{\eta}_v$ is given by equation (3) by replacing q_u with p _u, which is an m-flat submanifold as it has constraints on the $\newcommand{\e}{{\rm e}} \eta$ coordinate.

The dually flat structure of $\newcommand{\Pf}{\boldsymbol{\mathcal{S}}} \boldsymbol{\Pf}$ with the dual coordinate systems $\newcommand{\e}{{\rm e}} (\,(\theta_v)_{v \in \Omega^+}, (\eta_v)_{v \in \Omega^+}\,)$ leads to the following strong property: if A  =  B, that is, $(\Omega^+ \setminus B) \cup A = \Omega^+$ and $\newcommand{\e}{{\rm e}} (\Omega^+ \setminus B) \cap A = \emptyset$, the intersection $\newcommand{\Pf}{\boldsymbol{\mathcal{S}}} \newcommand{\ten}[1]{\mathcal{#1}} \boldsymbol{\Pf}_{B} \cap \boldsymbol{\Pf}_{\ten{P}}$ is always a singleton; that is, the tensor $\newcommand{\ten}[1]{\mathcal{#1}} \ten{Q}$ such that $\newcommand{\Pf}{\boldsymbol{\mathcal{S}}} \newcommand{\ten}[1]{\mathcal{#1}} \ten{Q} \in \boldsymbol{\Pf}_{B} \cap \boldsymbol{\Pf}_{\ten{P}}$ always uniquely exists, and $\newcommand{\ten}[1]{\mathcal{#1}} \ten{Q}$ is the minimizer of the KL divergence from $\newcommand{\ten}[1]{\mathcal{#1}} \ten{P}$ (Amari (2009), theorem 3):

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\KL}{D_{{\rm KL}}} \newcommand{\argmin}{\mathop{{\rm argmin}}} \newcommand{\ten}[1]{\mathcal{#1}} \newcommand{\Pf}{\boldsymbol{\mathcal{S}}} \displaystyle \label{eq:optimization} \ten{Q} = \argmin_{\ten{R} \in \boldsymbol{\Pf}_{B}} \KL(\ten{P}, \ten{R}). \nonumber \end{align} \tag{ 4 }$$

The transition from $\newcommand{\ten}[1]{\mathcal{#1}} \ten{P}$ to $\newcommand{\ten}[1]{\mathcal{#1}} \ten{Q}$ is called the m-projection of $\newcommand{\ten}[1]{\mathcal{#1}} \ten{P}$ onto $\newcommand{\Pf}{\boldsymbol{\mathcal{S}}} \boldsymbol{\Pf}_{B}$, and Legendre decomposition coincides with m-projection (figure 2). In contrast, if some fully decomposable tensor $\newcommand{\Pf}{\boldsymbol{\mathcal{S}}} \newcommand{\ten}[1]{\mathcal{#1}} \ten{R} \in \boldsymbol{\Pf}_{B}$ is given, finding the intersection $\newcommand{\Pf}{\boldsymbol{\mathcal{S}}} \newcommand{\ten}[1]{\mathcal{#1}} \ten{Q} \in \boldsymbol{\Pf}_{B} \cap \boldsymbol{\Pf}_{\ten{P}}$ is called the e-projection of $\newcommand{\ten}[1]{\mathcal{#1}} \ten{R}$ onto $\newcommand{\Pf}{\boldsymbol{\mathcal{S}}} \newcommand{\ten}[1]{\mathcal{#1}} \boldsymbol{\Pf}_{\ten{P}}$. In practice, we use e-projection because the number of parameters to be optimized is $|B|$ in e-projection while it is $|\Omega \setminus B|$ in m-projection, and $|B| \leqslant |\Omega \setminus B|$ usually holds.

Zoom In Zoom Out Reset image size

The e-projection is always convex optimization as the e-flat submanifold $\newcommand{\Pf}{\boldsymbol{\mathcal{S}}} \boldsymbol{\Pf}_B$ is convex with respect to $(\theta_v)_{v \in \Omega^+}$. More precisely,

$$\begin{align*} \newcommand{\e}{{\rm e}} \renewcommand{\log}{\mathop{{\rm log}}} \newcommand{\re}{{\rm Re}} \newcommand{\KL}{D_{{\rm KL}}} \newcommand{\ten}[1]{\mathcal{#1}} \displaystyle \KL(\ten{P}, \ten{Q}) = \sum_{v \in \Omega} p_v \log \frac{p_v}{q_v} = \sum_{v \in \Omega} p_v \log p_v - \sum_{v \in \Omega} p_v \log q_v = - \sum_{v \in \Omega} p_v \log q_v - H(\ten{P}), \nonumber \end{align*}$$

where $\newcommand{\ten}[1]{\mathcal{#1}} H(\ten{P})$ is the entropy of $\newcommand{\ten}[1]{\mathcal{#1}} \ten{P}$ and independent of $(\theta_v)_{v \in \Omega^+}$. Since we have

$$\begin{align*} \newcommand{\e}{{\rm e}} \renewcommand{\exp}{\mathop{{\rm exp}}} \newcommand{\re}{{\rm Re}} \renewcommand{\log}{\mathop{{\rm log}}} \displaystyle - \sum_{v \in \Omega} p_v \log q_v = \sum_{v \in \Omega} p_v \left(\sum_{u \in B} \zeta(u, v) (- \theta_u) + \psi(\theta)\right),\quad \psi(\theta) = \log\sum_{v \in \Omega}\exp\left(\sum_{u \in B} \zeta(u, v)\theta(u)\right), \nonumber \end{align*}$$

$\psi(\theta)$ is convex and $\newcommand{\ten}[1]{\mathcal{#1}} \newcommand{\KL}{D_{{\rm KL}}} \KL(\ten{P}, \ten{Q})$ is also convex with respect to $(\theta_v)_{v \in \Omega^+}$.

Here we present our two gradient-based optimization algorithms to solve the KL divergence minimization problem in equation (4). Since the KL divergence $\newcommand{\ten}[1]{\mathcal{#1}} \newcommand{\KL}{D_{{\rm KL}}} \KL(\ten{P}, \ten{Q})$ is convex with respect to each $\theta_v$, the standard gradient descent shown in algorithm 1 can always find the global optimum, where $\newcommand{\re}{{\rm Re}} \renewcommand{\epsilon}{\varepsilon} \newcommand{\e}{{\rm e}} \epsilon > 0$ is a learning rate. Starting with some initial parameter set $(\theta_v)_{v \in B}$, the algorithm iteratively updates the set until convergence. The gradient of $\theta_w$ for each $w \in B$ is obtained as

$$\begin{align*} \newcommand{\e}{{\rm e}} \renewcommand{\log}{\mathop{{\rm log}}} \newcommand{\re}{{\rm Re}} \newcommand{\KL}{D_{{\rm KL}}} \newcommand{\ten}[1]{\mathcal{#1}} \displaystyle \frac{\partial}{\partial \theta_w} \KL(\ten{P}, \ten{Q}) &= - \frac{\partial}{\partial \theta_w} \sum_{v \in \Omega} p_v \log q_v = - \frac{\partial}{\partial \theta_w} \sum_{v \in \Omega} p_v \sum_{u \in B} \zeta(u, v) \theta_u + \frac{\partial}{\partial \theta_w} \sum_{v \in \Omega} p_v \psi(\theta)\nonumber \\ &= - \sum_{v \in \Omega} p_v \zeta(w, v) + \frac{\partial \psi(\theta)}{\partial \theta_w} = \eta_w - \hat{\eta}_w, \nonumber \end{align*}$$

where the last equation uses the fact that $\newcommand{\e}{{\rm e}} \partial \psi(\theta) / \partial \theta_w = \eta_w$ in theorem 2 in Sugiyama et al (2017). This equation also shows that the KL divergence $\newcommand{\ten}[1]{\mathcal{#1}} \newcommand{\KL}{D_{{\rm KL}}} \KL(\ten{P}, \ten{Q})$ is minimized if and only if $\newcommand{\e}{{\rm e}} \eta_v = \hat{\eta}_v$ for all $v \in B$. The time complexity of each iteration is $O(|\Omega| |B|)$, as that of computing $\newcommand{\ten}[1]{\mathcal{#1}} \ten{Q}$ from $(\theta_v)_{v \in B}$ (line 5 in algorithm 1) is $O(|\Omega| |B|)$ and computing $\newcommand{\e}{{\rm e}} (\eta_v)_{v \in B}$ from $\newcommand{\ten}[1]{\mathcal{#1}} \ten{Q}$ (line 6 in algorithm 1) is $O(|\Omega|)$. Thus the total complexity is $O(h |\Omega| |B|^2)$ with the number of iterations h until convergence.

Algorithm 1. Legendre decomposition by gradient descent.

Although gradient descent is an efficient approach, in Legendre decomposition, we need to repeat 'decoding' from $(\theta_v)_{v \in B}$ and 'encoding' to $\newcommand{\e}{{\rm e}} (\eta_v)_{v \in B}$ in each iteration, which may lead to a loss of efficiency if the number of iterations is large. To reduce the number of iterations to gain efficiency, we propose to use a natural gradient (Amari 1998), which is the second-order optimization method, shown in algorithm 2. Again, since the KL divergence $\newcommand{\ten}[1]{\mathcal{#1}} \newcommand{\KL}{D_{{\rm KL}}} \KL(\ten{P}, \ten{Q})$ is convex with respect to $(\theta_v)_{v \in B}$, a natural gradient can always find the global optimum. More precisely, our natural gradient algorithm is an instance of the Bregman algorithm applied to a convex region, which is well known to always converge to the global solution (Censor and Lent 1981). Let $B = \{v_1, v_2, \dots, v_{|B|}\}$, $\boldsymbol{\theta} = (\theta_{v_1}, \dots, \theta_{v_{|B|}}){}^T$, and $\newcommand{\e}{{\rm e}} \boldsymbol{\eta} = (\eta_{v_1}, \dots, \eta_{v_{|B|}}){}^T$. In each update of the current $\boldsymbol{\theta}$ to $\boldsymbol{\theta}_{{\rm next}}$, the natural gradient method uses the relationship,

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \Delta \boldsymbol{\eta} = - {\bf G} \Delta \boldsymbol{\theta},\quad \Delta \boldsymbol{\eta} = \boldsymbol{\eta} - \hat{\boldsymbol{\eta}} ~{\rm and}~ \Delta \boldsymbol{\theta} = \boldsymbol{\theta}_{{\rm next}} - \boldsymbol{\theta}, \nonumber \end{align*}$$

which leads to the update formula

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{\theta}_{{\rm next}} = \boldsymbol{\theta} - {\bf G}^{-1} \Delta \boldsymbol{\eta}, \nonumber \end{align*}$$

where $\newcommand{\R}{\mathbb{R}} {\bf G} = (g_{uv}) \in \R^{|B| \times |B|}$ is the Fisher information matrix such that

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\log}{\mathop{{\rm log}}} \newcommand{\re}{{\rm Re}} \newcommand{\Exp}{{\bf E}} \displaystyle \label{eq:Fisher_theta} g_{uv}(\theta) = \frac{\partial \eta_u}{\partial \theta_v} = \Exp\left[\frac{\partial\log p_w}{\partial\theta_u}\frac{\partial\log p_w}{\partial\theta_v}\right] = \sum_{w \in \Omega}\zeta(u, w)\zeta(v, w)\,p_w - \eta_u\eta_v \nonumber \end{align} \tag{ 5 }$$

as given in theorem 3 in Sugiyama et al (2017). Note that natural gradient coincides with Newton's method in our case as the Fisher information matrix ${\bf G}$ corresponds to the (negative) Hessian matrix:

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\KL}{D_{{\rm KL}}} \newcommand{\ten}[1]{\mathcal{#1}} \displaystyle \frac{\partial^2}{\partial \theta_{u} \partial \theta_{v}} \KL(\ten{P}, \ten{Q}) = - \frac{\partial\eta_{u}}{\partial \theta_{v}} = - g_{uv}. \nonumber \end{align*}$$

Algorithm 2. Legendre decomposition by natural gradient.

The time complexity of each iteration is $O(|\Omega| |B| + |B|^3)$, where $O(|\Omega| |B|)$ is needed to compute $\newcommand{\ten}[1]{\mathcal{#1}} \ten{Q}$ from $\theta$ and O(|B|³) to compute the inverse of ${\bf G}$, resulting in the total complexity $O(h' |\Omega| |B| + h'|B|^3)$ with the number of iterations $h'$ until convergence.

We demonstrate interesting relationships between Legendre decomposition and statistical models, including the exponential family and the Boltzmann (Gibbs) distributions, and show that our decomposition method can be viewed as a generalization of Boltzmann machine learning (Ackley et al 1985). Although the connection between tensor decomposition and graphical models has been analyzed by Yılmaz et al (2011), Yılmaz and Cemgil (2012), and Chen et al (2018), our analysis adds a new insight as we focus on not the graphical model itself but the sample space of distributions generated by the model.

2.4.1. Exponential family.

We show that the set of normalized tensors $\newcommand{\R}{\mathbb{R}} \newcommand{\Pf}{\boldsymbol{\mathcal{S}}} \newcommand{\ten}[1]{\mathcal{#1}} \boldsymbol{\Pf} = \{\ten{P} \in \R_{> 0}^{I_1 \times I_2 \times \dots \times I_N} \mid \sum\nolimits_{v \in \Omega} p_{v} = 1\}$ is included in the exponential family. The exponential family is defined as

$$\begin{align*} \newcommand{\e}{{\rm e}} \renewcommand{\exp}{\mathop{{\rm exp}}} \newcommand{\re}{{\rm Re}} \displaystyle p(x, \boldsymbol{\theta}) = \exp\left(\sum \theta_i k_i(x) + r(x) - C(\boldsymbol{\theta})\right) \nonumber \end{align*}$$

for natural parameters $\boldsymbol{\theta}$. Since our model in equation (1) can be written as

$$\begin{align*} \newcommand{\e}{{\rm e}} \renewcommand{\exp}{\mathop{{\rm exp}}} \newcommand{\re}{{\rm Re}} \displaystyle p_{v} = \frac{1}{\exp(\psi(\theta))}\prod_{u \in {\downarrow}v} \exp(\theta_u) = \exp\left(\sum_{u \in \Omega^+}\theta_u\zeta(u, v) - \psi(\theta)\right) \nonumber \end{align*}$$

with $\theta_u = 0$ for $u \in \Omega^+ \setminus B$, it is clearly in the exponential family, where $\zeta$ and $\psi(\theta)$ correspond to k and $C(\boldsymbol{\theta})$, respectively, and $r(x) = 0$. Thus, the $(\theta_v)_{v \in B}$ used in Legendre decomposition are interpreted as natural parameters of the exponential family. Moreover, we can obtain $\newcommand{\e}{{\rm e}} (\eta_v)_{v \in B}$ by taking the expectation of $\zeta(u, v)$:

$$\begin{align*} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle {\bf E}\bigl[\zeta(u, v)\bigr] = \sum_{v \in \Omega}\zeta(u, v)\,p_{v} = \eta_u. \nonumber \end{align*}$$

Thus Legendre decomposition of $\newcommand{\ten}[1]{\mathcal{#1}} \ten{P}$ is understood to find a fully decomposable $\newcommand{\ten}[1]{\mathcal{#1}} \ten{Q}$ that has the same expectation with $\newcommand{\ten}[1]{\mathcal{#1}} \ten{P}$ with respect to a basis B.

2.4.2. Boltzmann machines.

A Boltzmann machine is represented as an undirected graph $G = (V, E)$ with a vertex set $V$ and an edge set $E \subseteq V \times V$, where we assume that $V = [N] = \{1, 2, \dots, N\}$ without loss of generality. This $V$ is the set of indices of N binary variables. A Boltzmann machine G defines a probability distribution P, where each probability of an N-dimensional binary vector $\boldsymbol{x} \in \{0, 1\}^{N}$ is given as

$$\begin{align*} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \renewcommand{\exp}{\mathop{{\rm exp}}} \newcommand{\re}{{\rm Re}} \displaystyle p(\boldsymbol{x}) = \frac{1}{Z(\theta)} \prod_{{i \in V}}\exp\big(\theta_i x_i\big) \prod_{{\{i, j\} \in E}} \exp\big(\theta_{ij} x_i x_j\big), \nonumber \end{align*}$$

where $\theta_i$ is a bias, $\theta_{ij}$ is a weight, and $Z(\theta)$ is a partition function.

To translate a Boltzmann machine into our formulation, let $\Omega = \{1, 2\}^N$ and suppose that

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle B(V) &= \{(i_{1}^a, \dots, i_{N}^a) \in \Omega | a \in V\},\qquad & i_{l}^a = \left\{\begin{array}{@{}ll@{}} 2 & {\rm if}~ l = a, \nonumber \\ 1 & {\rm otherwise}, \end{array} \right.\nonumber \\ B(E) &= \{(i_{1}^{ab}, \dots, i_{N}^{ab}) \in \Omega | \{a, b\} \in E\},\qquad & i_{l}^{ab} = \left\{\begin{array}{@{}ll@{}} 2 & {\rm if}~ l \in \{a, b\}, \nonumber \\ 1 & {\rm otherwise}. \end{array} \right. \nonumber \end{align*}$$

Then it is clear that the set of probability distributions, or Gibbs distributions, that can be represented by the Boltzmann machine G is exactly the same as $\newcommand{\Pf}{\boldsymbol{\mathcal{S}}} \boldsymbol{\Pf}_{B}$ with $B = B(V) \cup B(E)$ and $\newcommand{\re}{{\rm Re}} \renewcommand{\exp}{\mathop{{\rm exp}}} \newcommand{\e}{{\rm e}} \exp(\psi(\theta)) = Z(\theta)$; that is, the set of fully decomposable Nth-order tensors defined by equation (1) with the basis $B(V) \cup B(E)$ (figure 3). Moreover, let a given Nth-order tensor $\newcommand{\R}{\mathbb{R}} \newcommand{\ten}[1]{\mathcal{#1}} \ten{P} \in \R_{\geqslant 0}^{2 \times 2 \times \dots \times 2}$ be an empirical distribution estimated from data, where each $p_{v}$ is the probability for a binary vector $v - (1, \dots, 1) \in \{0, 1\}^N$. The tensor $\newcommand{\ten}[1]{\mathcal{#1}} \ten{Q}$ obtained by Legendre decomposition with $B = B(V) \cup B(E)$ coincides with the distribution learned by the Boltzmann machine $G = (V, E)$. The condition $\newcommand{\e}{{\rm e}} \eta_v = \hat{\eta}_v$ in the optimization of the Legendre decomposition corresponds to the well-known learning equation of Boltzmann machines, where $\newcommand{\e}{{\rm e}} \hat{\eta}$ and $\newcommand{\e}{{\rm e}} \eta$ correspond to the expectation of the data distribution and that of the model distribution, respectively.

Zoom In Zoom Out Reset image size

Therefore our Legendre decomposition is a generalization of Boltzmann machine learning in the following three aspects:

• 1.  
  The domain is not limited to binary but can be ordinal; that is, {0,1}^N is extended to $[I_1] \times [I_2] \times \dots \times [I_N]$ for any $\newcommand{\N}{\mathbb{N}} I_1, I_2, \dots, I_N \in \N$.
• 2.  
  The basis B with which parameters $\theta$ are associated is not limited to $B(V) \cup B(E)$ but can be any subset of $[I_1] \times \dots \times [I_N]$, meaning that higher-order interactions (Sejnowski 1986) can be included.
• 3.  
  The sample space of probability distributions is not limited to {0,1}^N but can be any subset of $[I_1] \times [I_2] \times \dots \times [I_N]$, which enables us to perform efficient computation by removing unnecessary entries such as missing values.

Hidden variables are often used in Boltzmann machines to increase the representation power, such as in restricted Boltzmann machines (Smolensky 1986, Hinton 2002) and deep Boltzmann machines (Salakhutdinov and Hinton 2009, 2012). In Legendre decomposition, including a hidden variable corresponds to including an additional dimension. Hence if we include H hidden variables, the fully decomposable tensor $\newcommand{\ten}[1]{\mathcal{#1}} \ten{Q}$ has the order of N  +  H. This is an interesting extension to our method and an ongoing research topic, but it is not a focus of this paper since our main aim is to find a lower dimensional representation of a given tensor $\newcommand{\ten}[1]{\mathcal{#1}} \ten{P}$.

In the learning process of Boltzmann machines, approximation techniques of the partition function $Z(\theta)$ are usually required, such as annealed importance sampling (Salakhutdinov and Murray 2008) or variational techniques (Salakhutdinov 2008). This requirement is because the exact computation of the partition function requires the summation over all probabilities of the sample space $\Omega$, which is always fixed to $2^V$ with the set $V$ of variables in the learning process of Boltzmann machines, and which is not tractable. Our method does not require such techniques as $\Omega$ is a subset of indices of an input tensor and the partition function can always be directly computed.