Dimensionality reduction and topographic mapping of binary tensors

Abstract

In this paper, a decomposition method for binary tensors, generalized multi-linear model for principal component analysis (GMLPCA) is proposed. To the best of our knowledge at present there is no other principled systematic framework for decomposition or topographic mapping of binary tensors. In the model formulation, we constrain the natural parameters of the Bernoulli distributions for each tensor element to lie in a sub-space spanned by a reduced set of basis (principal) tensors. We evaluate and compare the proposed GMLPCA technique with existing real-valued tensor decomposition methods in two scenarios: (1) in a series of controlled experiments involving synthetic data; (2) on a real-world biological dataset of DNA sub-sequences from different functional regions, with sequences represented by binary tensors. The experiments suggest that the GMLPCA model is better suited for modelling binary tensors than its real-valued counterparts. Furthermore, we extended our GMLPCA model to the semi-supervised setting by forcing the model to search for a natural parameter subspace that represents a user-specified compromise between the modelling quality and the degree of class separation.

Appendix: Parameter estimation

To get analytical parameter updates, we use the trick of [11] and take advantage of the fact that while the model log-likelihood (7) is not convex in the parameters, it is convex in any parameter, if the others are kept fixed. This leads to an iterative estimation scheme detailed below.

The analytical updates will be derived from a lower bound on the log-likelihood (7) using [11]:

$$\log \sigma(\hat \theta) \ge -\log 2 + \frac{\hat \theta}{2} - \log \cosh \left( \frac{\theta}{2}\right) - (\hat \theta^2 - \theta^2) \,\frac{\tanh \frac{\theta}{2}}{4 \theta},$$

(31)

where $\theta$ stands for the current value of individual natural parameters $\theta_{m,{\bf i}}$ of the Bernoulli noise models ${P({\mathcal A}_{m,{\bf i}} | \theta_{m,{\bf i}})}$ and $\hat \theta$ stands for the future estimate of the parameter, given the current parameter values. Hence, from (7) we obtain⁴

$${\mathcal {L}}(\hat \Uptheta) = \sum_{m=1}^M \sum_{{\bf i} \in \Upupsilon} {\mathcal A}_{m,{\bf i}} \log \sigma(\hat \theta_{m,{\bf i}}) + (1-{\mathcal A}_{m,{\bf i}}) \log \sigma(- \hat \theta_{m,{\bf i}}) \ge \sum_{m=1}^M \sum_{{\bf i} \in \Upupsilon} {\mathcal A}_{m,{\bf i}} \left[ - \log 2 + \frac{\hat \theta_{m,{\bf i}}}{2} \right. \left.- \log \cosh \left( \frac{\theta_{m,{\bf i}}}{2}\right) - (\hat \theta_{m,{\bf i}}^2 - \theta_{m,{\bf i}}^2) \frac{\tanh \frac{\theta_{m,{\bf i}}}{2}}{4 \theta_{m,{\bf i}}} \right] + (1-{\mathcal A}_{m,{\bf i}}) \left[ -\log 2 - \frac{\hat \theta_{m,{\bf i}}}{2} \right. \left.- \log \cosh \left( \frac{\theta_{m,{\bf i}}}{2}\right) - (\hat \theta_{m,{\bf i}}^2 - \theta_{m,{\bf i}}^2) \frac{\tanh \frac{\theta_{m,{\bf i}}}{2}}{4 \theta_{m,{\bf i}}} \right]$$

(32)

$$= \,H(\hat \Uptheta, \Uptheta).$$

(33)

Denote $(\tanh \frac{\theta_{m,{\bf i}}}{2}) / \theta_{m,{\bf i}}$ by $\Uppsi_{m,{\bf i}}.$ Grouping together constant terms in (32) leads to

$$H(\hat \Uptheta, \Uptheta) = \sum_{m=1}^M \sum_{{\bf i} \in \Upupsilon} \left[ \hat \theta_{m,{\bf i}} \left({\mathcal A}_{m,{\bf i}} - \frac{1}{2}\right) - \frac{\Uppsi_{m,{\bf i}}}{4} \,\hat \theta_{m,{\bf i}}^2 \right] + \,Const. $$

(34)

Note that ${H(\hat \Uptheta, \Uptheta) = {\mathcal {L}}(\hat \Uptheta)}$ only if $\hat \Uptheta = \Uptheta.$ Therefore, by choosing $\hat \Uptheta$ that maximizes $H(\hat \Uptheta, \Uptheta)$ we guarantee ${{\mathcal {L}}(\hat \Uptheta) \ge H(\hat \Uptheta, \Uptheta) \ge H(\Uptheta, \Uptheta) = {\mathcal {L}}(\Uptheta)}$ [11].

We are now ready to constrain the Bernoulli parameters to be optimized [see (9)]:

$$\hat \theta_{m,{\bf i}} = \sum_{{\bf r} \in \rho} {\mathcal {Q}}_{m,{\bf r}} \cdot \prod_{n=1}^N u^{(n)}_{r_n,i_n} + \Updelta_{{\bf i}}.$$

(35)

We will update the model parameters so as to maximize

$${\mathcal{H}}= \sum_{m=1}^M \sum_{{\bf i} \in \Upupsilon} {\mathcal{H}}_{m,{\bf i}},$$

(36)

where

$${\mathcal{H}}_{m,{\bf i}} = \left({\mathcal A}_{m,{\bf i}} - \frac{1}{2}\right) \hat \theta_{m,{\bf i}} - \frac{\Uppsi_{m,{\bf i}}}{4} \,\hat \theta_{m,{\bf i}}^2,$$

(37)

with $\hat \theta_{m,{\bf i}}$ given by (35).

1.1 Updates for n-mode space basis

When updating the n-mode space basis {u ⁽ⁿ⁾₁ , u ⁽ⁿ⁾₂ , ..., u ⁽ⁿ⁾ _{R_n} }, the bias tensor $\Updelta$ and the expansion coefficients ${{\mathcal {Q}}_{m,{\bf r}}, }$ m = 1, 2, ...M, r  ∈ ρ, are kept fixed to their current values.

For n = 1, 2, ..., N, define

$$\Upupsilon_{-n} = \{1,2,\ldots,I_1\} \times\cdots \times \{1,2,\ldots,I_{n-1}\} \times \{ 1 \} \times \{1,2,\ldots,I_{n+1}\} \times\cdots\times \{1,2,\ldots,I_N\},$$

(38)

with obvious interpretation in the boundary cases. Given ${\bf\ {{i}}} \in \Upupsilon_{-n}$ and an n-mode index $j \in \{1,2,\ldots,I_{n}\},$ the index N-tuple (i ₁, ..., i _n-1, j, i _n+1, ..., i _N ) formed by inserting j at the nth place of i is denoted by [i, j|n].

In order to evaluate

$$\frac{\partial \,{\mathcal{H}}}{\partial \,u^{(n)}_{q,j}}, \quad q=1,2,\ldots,R_n, \,j = 1,2,\ldots,I_n,$$

we realize that $u^{(n)}_{q,j}$ is involved in expressing all $\hat \theta_{m,[{\bf i}, j|n]},$ m = 1, 2, ..., M, with ${\bf\ {{i}}} \in \Upupsilon_{-n}.$ Therefore,

$$\frac{\partial \,{\mathcal{H}}}{\partial \,u^{(n)}_{q,j}} = \sum_{m=1}^M \,\sum_{{\bf i} \in \Upupsilon_{-n}} \frac{\partial \,{\mathcal{H}}_{m,[{\bf i}, j|n] }} {\partial \,\hat \theta_{m,[{\bf i}, j|n]}} \,\frac{\partial \,\hat \theta_{m,[{\bf i}, j|n]}}{\partial \,u^{(n)}_{q,j}},$$

(39)

where

$$\frac{\partial \,{\mathcal{H}}_{m,[{\bf i}, j|n] }} {\partial \,\hat \theta_{m,[{\bf i}, j|n]}} = \left({\mathcal A}_{m,[{\bf i}, j|n]} - \frac{1}{2}\right) - \frac{\Uppsi_{m, [{\bf i}, j|n] }}{2} \,\hat \theta_{m,[{\bf i}, j|n]}$$

(40)

and from (35),

$$\frac{\partial \,\hat \theta_{m,[{\bf i}, j|n]}}{\partial \,u^{(n)}_{q,j}} = {\mathcal {B}}^{(n)}_{m,{\bf i},q} = \sum_{{\bf r} \in \rho_{-n}} {\mathcal {Q}}_{m, [{\bf r}, q|n]} \cdot \prod_{s=1, s \neq n}^N u^{(s)}_{r_s,i_s}.$$

(41)

Here, the index set $\rho_{-n}$ is defined analogously to $\Upupsilon_{-n}:$

$$\rho_{-n} = \{1,2,...,R_1\} \times\cdots\times \{1,2,\ldots,R_{n-1}\} \times \{ 1 \} \times \{1,2,\ldots,R_{n+1}\} \times\cdots\times \{1,2,\ldots,R_N\}.$$

(42)

Setting the derivative (39) to zero results in

$${\sum_{m=1}^M \,\sum_{{\bf i} \in \Upupsilon_{-n}} (2 {\mathcal A}_{m,[{\bf i}, j|n]} - 1) \,{\mathcal {B}}^{(n)}_{m,{\bf i},q} = } \sum_{m=1}^M \ \sum_{{\bf i} \in \Upupsilon_{-n}} \Uppsi_{m, [{\bf i}, j|n]} \ \hat \theta_{m,[{\bf i}, j|n]} \ {\mathcal {B}}^{(n)}_{m,{\bf i},q}.$$

(43)

Rewriting (35) as

$$\hat \theta_{m,[{\bf i},j|n]} = \sum_{t = 1}^{R_n} \ \sum_{{\bf r} \in \rho_{-n}} {\mathcal {Q}}_{m, [{\bf r},t|n] } \ u^{(n)}_{t,j} \prod_{s=1, s \neq n}^N u^{(s)}_{r_s,i_s} \\ + \ \Updelta_{[{\bf i},j|n]}$$

(44)

and applying to (43) we obtain

$$\sum_{t = 1}^{R_n} u^{(n)}_{t,j} \ {\mathcal {K}}^{(n)}_{q,t,j} = {\mathcal {S}}^{(n)}_{q,j},$$

(45)

where

$${\mathcal {S}}^{(n)}_{q,j} = \sum_{m=1}^M \sum_{{\bf i} \in \Upupsilon_{-n}} (2 {\mathcal A}_{m,[{\bf i}, j|n]} - 1 - \Uppsi_{m, [{\bf i}, j|n]} \Updelta_{[{\bf i}, j|n]}) {\mathcal {B}}^{(n)}_{m,{\bf i},q},$$

(46)

and

$${\mathcal {K}}^{(n)}_{q,t,j} = \sum_{m=1}^M \sum_{{\bf r} \in \rho_{-n}} {\mathcal {Q}}_{m, [{\bf r},t|n] } \times \sum_{{\bf i} \in \Upupsilon_{-n}} \Uppsi_{m, [{\bf i}, j|n]} \ {\mathcal {B}}^{(n)}_{m,{\bf i},q} \prod_{s=1, s \neq n}^N u^{(s)}_{r_s,i_s}.$$

(47)

For each n-mode coordinate $j \in \{ 1,2,\ldots,I_n\},$ collect the j-th coordinate values of all n-mode basis vectors into a column vector ${\bf u}^{(n)}_{:,j} = (u^{(n)}_{1,j}, u^{(n)}_{2,j}, \ldots, u^{(n)}_{R_n,j} )^T.$ Analogously, stack all the ${{\mathcal {S}}^{(n)}_{q,j}}$ values in a column vector ${{\mathcal {S}}^{(n)}_{:,j} = ({\mathcal {S}}^{(n)}_{1,j}, {\mathcal {S}}^{(n)}_{2,j}, \ldots,}$ ${\mathcal {S}}^{(n)}_{R_n,j})^T.$ Finally, we construct an $R_n \times R_n$ matrix ${\mathcal {K}}^{(n)}_{:,:,j}$ whose q-th row is $({\mathcal {K}}^{(n)}_{q,1,j}, {\mathcal {K}}^{(n)}_{q,2,j},\ldots, {\mathcal {K}}^{(n)}_{q,R_n,j}), q = 1,2,\ldots,R_n.$ The n-mode basis vectors are updated by solving I _n linear systems of size $R_n \times R_n$:

$${\mathcal {K}}^{(n)}_{:,\,:,\,j} \, {\bf u}^{(n)}_{:,\,j} = {\mathcal {S}}^{(n)}_{:,\,j},\quad j = 1,2,\ldots,I_n.$$

(48)

1.2 Updates for expansion coefficients

When updating the expansion coefficients ${{\mathcal {Q}}_{m,{\bf\ {{r}}}}, }$ the bias tensor $\Updelta$ and the basis sets $\{ {\bf\ {{u}}}^{(n)}_1, {\bf u}^{(n)}_2,\ldots, {\bf u}^{(n)}_{R_n} \}$ for all n modes n = 1, 2, ..., N are kept fixed to their current values.

For ${\bf r} \in \rho$ and ${\bf i} \in \Upupsilon$ denote $\prod_{n=1}^N u^{(n)}_{r_n,i_n}$ by $C_{{\bf r},{\bf i}}.$ For data index $\ell = 1,2,\ldots,M$ and basis index ${\bf\ {{v}}} \in \rho$ we have

$$\frac{\partial \ {\mathcal{H}}}{\partial \ {\mathcal {Q}}_{\ell,{\bf v}}} = \sum_{m=1}^M \ \sum_{{\bf i} \in \Upupsilon} \frac{\partial \ {\mathcal{H}}_{m,{\bf i}}} {\partial \ \hat \theta_{m,{\bf i}}} \ \frac{\partial \ \hat \theta_{m,{\bf i}}}{\partial \ {\mathcal {Q}}_{\ell,{\bf v}}},$$

(49)

where

$$\frac{\partial \ {\mathcal{H}}_{m,{\bf i}}} {\partial \ \hat \theta_{m,{\bf i}}} = \left({\mathcal A}_{m,{\bf i}} - \frac{1}{2}\right) - \frac{\Uppsi_{m, {\bf i} }}{2} \ \hat \theta_{m,{\bf i}}$$

(50)

and ${ \frac{\partial \ \hat \theta_{m,{\bf i}}}{\partial \ {\mathcal {Q}}_{\ell,{\bf v}}} = C_{{\bf v},{\bf i}}}$ if $m=\ell$ and ${\frac{\partial \ \hat \theta_{m,{\bf\ {{i}}}}}{\partial \ {\mathcal {Q}}_{\ell,{\bf v}}} = 0}$ otherwise.

By imposing ${\frac{\partial \ \mathcal{H}}{\partial \ {\mathcal {Q}}_{\ell,{\bf v}}} = 0, }$ we get

$${\mathcal {T}}_{{\bf v},\ell} = \sum_{{\bf r} \in \rho} {\mathcal {P}}_{{\bf v},{\bf r},\ell} \ {\mathcal {Q}}_{\ell,{\bf r}},$$

(51)

where

$${\mathcal {T}}_{{\bf v},\ell} = \sum_{{\bf i} \in \Upupsilon} (2 {\mathcal A}_{\ell,{\bf i}} - 1 - \Uppsi_{\ell,{\bf i}} \ \Updelta_{{\bf i}}) \ C_{{\bf v},{\bf i}}$$

(52)

and

$${\mathcal {P}}_{{\bf v},{\bf r},\ell} = \sum_{{\bf i} \in \Upupsilon} \Uppsi_{\ell,{\bf i}} \ C_{{\bf v},{\bf i}} \ C_{{\bf r},{\bf i}}.$$

(53)

To solve for expansion coefficients using the tools of matrix algebra, we need to vectorize tensor indices. Consider any one-to-one function $\kappa$ from $\rho$ to $\{ 1,2,\ldots,\prod_{n=1}^N R_n \}.$ For each input tensor index $\ell = 1,2,\ldots,M,$

• create a square $(\prod_{n=1}^N R_n) \times (\prod_{n=1}^N R_n)$matrix ${{\mathcal {P}}_{:,:,\ell}}$ whose $(\kappa({\bf v}), \kappa({\bf r}))$-th element is equal to ${{\mathcal {P}}_{{\bf\ {{v}}},{\bf r},\ell}, }$

• stack the values of ${{\mathcal {T}}_{{\bf v},\ell}}$ into a column vector ${{\mathcal {T}}_{:,\ell}}$ whose $\kappa({\bf\ {{v}}})$-th coordinate is ${{\mathcal {T}}_{{\bf v},\ell}, }$

• collect the expansion coefficients ${{\mathcal {Q}}_{\ell,{\bf\ {{r}}}}}$ in a column vector ${{\mathcal {Q}}_{\ell,:}}$ with $\kappa({\bf r})$-th coordinate equal to ${{\mathcal {Q}}_{\ell,{\bf r}}. }$

The expansion coefficients for the $\ell$-th input tensor ${{\mathcal A}_\ell}$ can be obtained by solving

$${\mathcal {P}}_{:,:,\ell} \ {\mathcal {Q}}_{\ell,:} = {\mathcal {T}}_{:,\ell}, \ \ \ \ell = 1,2,\ldots,M.$$

(54)

1.3 Updates for the bias tensor

As before, when updating the bias tensor $\Updelta,$ the expansion coefficients ${{\mathcal {Q}}_{m,{\bf r}}, }$ m = 1, 2, ...M, ${\bf r} \in \rho,$ and the basis sets $\{ {\bf u}^{(n)}_1, {\bf u}^{(n)}_2, \ldots, {\bf u}^{(n)}_{R_n} \}$ for all n modes n = 1, 2, ..., N are kept fixed to their current values.

Fix ${\bf j} \in \Upupsilon.$We evaluate

$$\frac{\partial \ {\mathcal{H}}}{\partial \ \Updelta_{{\bf j}}} = \sum_{m=1}^M \ \sum_{{\bf i} \in \Upupsilon} \frac{\partial \ {\mathcal{H}}_{m,{\bf i}}} {\partial \ \hat \theta_{m,{\bf i}}} \ \frac{\partial \ \hat \theta_{m,{\bf i}}}{\partial \ \Updelta_{{\bf j}}},$$

(55)

where $\frac{\partial \ \hat \theta_{m,{\bf i}}}{\partial \ \Updelta_{{\bf j}}}$ is equal to 1 if i = j and 0 otherwise.

Solving for $\frac{\partial \ \mathcal{H}}{\partial \ \Updelta_{{\bf j}}} = 0$ leads to

$$\Updelta_{{\bf j}} = \frac{\sum_{m=1}^M 2 {\mathcal A}_{m,{\bf j}} -1 - \Uppsi_{m,{\bf j}} \cdot \sum_{{\bf r} \in \rho} {\mathcal {Q}}_{m,{\bf r}} \ C_{{\bf r},{\bf j}} } {\sum_{m=1}^M \Uppsi_{m,{\bf j}}}.$$

(56)