To model the above tensor, previous works (Kolda and Bader,, 2009) assume that

• •

  There are a set of rank-one tensor components $\mathcal{E}=\{\mathcal{E}_{1},\dots,$ $\mathcal{E}_{R}\}$, which are learnable;

• •

  The tensor data sample/slice $\mathcal{T}^{(n)}$ admits a low-rank structure and can be represented as a weighted sum of these tensor components $\mathcal{E}$, where $\mathbf{x}^{(n)}$ denotes its coefficient vector;

• •

  On top of the low-rank structure, each data sample $\mathcal{T}^{(n)}$ also contains additional element-wise i.i.d. Gaussian noise due to real-world distortion (e.g., physical noise in signal measurements).

In the context of downstream classifications, we further assume that each sample $\mathcal{T}^{(n)}$ is semantically associated to one of the latent classes $p\in\{1,\dots,P\}$, and we let $\mathcal{D}_{p}$ be the sample distribution of class-$p$. Thus, the data sample can be formulated as, $\forall n$,

where

Here $\mathbf{a}_{r},\mathbf{b}_{r},\mathbf{c}_{r}$ are the learnable parameters, which produces the rank-one component $\mathcal{E}_{r}$, and they are the column vectors of $\{\mathbf{A}\in\mathbb{R}^{I\times R},\mathbf{B}\in\mathbb{R}^{J\times R},\mathbf{C}\in\mathbb{R}^{K\times R}\}$, referred as "bases".

Given the above tensor model, standard CPD only captures the i.i.d. Gaussian noise by minimizing the negative log-likelihood (NLL), which results in the following standard fitness/reconstruction loss,

Here, the Kruskal product $\llbracket\cdot\rrbracket$ outputs a fourth-order reconstructed tensor from four input factor matrices. For consistency, if the first input is a vector, the output is considered as a third-order tensor.

This paper considers two design principles for good bases. The first principle is fitness, which requires a low-rank tensor reconstruction with the bases. Second, data samples associated with the same latent class should be decomposed into similar coefficient vectors, with the bases, while the vectors should be dissimilar if the samples are from different latent classes. This principle is called alignment, which is important for classification but not considered in the standard tensor decomposition. In this paper, we assess the quality of the learned bases by the performance of downstream classification, where the coefficient vectors (obtained using the bases via decomposition) are the feature inputs (into the downstream classifier).

To put it succinctly, the paper tackles an unsupervised learning problem while using downstream supervised classification for evaluation. The procedures are briefly outlined:

• •

  First, we learn the bases $\{\mathbf{A},\mathbf{B},\mathbf{C}\}$ from a large set of unlabeled data. The loss function is developed in consideration of the fitness and alignment (defined in the next section) principles.

• •

  Then, we construct the following feature extractor given $\{\mathbf{A},\mathbf{B},\mathbf{C}\}$. The feature vector of a new sample is obtained by the closed-form solution of the least squares problem ($\alpha>0$ is a hyperparameter),

  $\displaystyle\mathbf{f}(\mathcal{T}^{(new)};\mathbf{A},\mathbf{B},\mathbf{C})=\operatorname*{arg\,min}_{\mathbf{x}\in\mathbb{R}^{1\times R}}\left(\left\|\mathcal{T}^{(new)}-\llbracket\mathbf{x},\mathbf{A},\mathbf{B},\mathbf{C}\rrbracket\right\|_{F}^{2}+\alpha\|\mathbf{x}\|_{2}^{2}\right).$

  (2)

  Note that, when $\mathbf{f}(\cdot)$ is applied to a batch of samples, it outputs a coefficient matrix.

• •

  Next, we evaluate the feature extractor with a set of labeled data. Given $\mathbf{f}(\cdot)$, we first apply it on the labeled data to extract their features and then train an additional logistic regression model (as the downstream classifier) on top of the extracted features, so that the result of classifications will implicitly reflect how good the bases are.

In this paper, we consider using the law of total probability to construct the negative sampler in a statistical way. Formally, assume $r_{p}$ is the label rate of latent class-$p$ (thus, we have $\sum_{p}r_{p}=1$), we apply the law of total probability and the following holds,

Here, the usage of $\mathbb{E}[\cdot]$ means that the expectation is taken over four interdependent random variables, i.e., $p,q,\mathcal{X}_{p},\mathcal{Y}_{q}$, while $\mathbb{E}_{p,q}[\cdot]$ means that $p,q$ is fixed and thus it is only taken over two random variables, i.e., $\mathcal{X}_{p},\mathcal{Y}_{q}$. The result shows that the negative sampler can be equivalently replaced by a weighted combination of the random sampler and positive sampler. Here we do not have access to $r_{p},~{}\forall p$ with unlabeled data, this issue is dealt with later.

Consequently, we can reformulate our self-supervised loss as,

From Eqn. (5) to Eqn. (6), we use the results in Eqn. (4)

The above form still requires label rate information, i.e., $r_{p},~{}\forall p$, we therefore consider using the following approximation to the above loss $\mathcal{L}_{ss}$,

Here, $\gamma\geq 0$ is a hyperparameter, while $\mathcal{L}_{ss}$ can be bounded as (derivations in appendix E),

The equivalence is established when $C_{1}=C_{2}$, i.e., the class labels are balanced. To simplify the derivation, we ignore $\lambda$ in the following and let $\gamma$ be a new hyperparameter. Also, the constants $C_{1}$ and $C_{2}$ can be absorbed into a weight hyperparameter $\beta$, given in the next section. This bound implies that, an easy-to-compute $\beta\mathcal{L}_{ss}^{\Theta}(\gamma)$ is often a good approximation of $\mathcal{L}_{ss}$ for some $\beta$. The next section specifies how to compute $\beta\mathcal{L}_{ss}^{\Theta}(\gamma)$ unsupervisedly as our empirical self-supervised loss.

We obtain an empirical estimator of $\mathcal{L}^{\Theta}_{ss}$ with Monte Carlo method. Suppose $\mathcal{T}$ and $\tilde{\mathcal{T}}$ are the input tensor and the augmented tensor respectively, and $\mathbf{X}=\mathbf{f}(\mathcal{T}),\tilde{\mathbf{X}}=\mathbf{f}(\tilde{\mathcal{T}})\in\mathbb{R}^{N\times R}$ are the coefficient/feature matrices. We use the row vectors of $\mathbf{X},\tilde{\mathbf{X}}$ to estimate Eqn. (7).

The first term $\mathbb{E}\left[\mbox{sim}\left(\mathbf{f}\left(\mathcal{X}_{p}\right),\mathbf{f}\left({\mathcal{Y}_{q}}\right)\right)\right]$ essentially requires a random sampler, which is approximated by the average cosine similarity of all possible pairs of non-corresponding row vectors of $\mathbf{X},\tilde{\mathbf{X}}$, while the second term $\mathbb{E}\left[\mbox{sim}\left(\mathbf{f}\left(\mathcal{X}_{p}\right),\mathbf{f}\left({\mathcal{Y}_{q}}\right)\right)\mid p=q\right]$ requires a positive sampler, which is estimated by the average cosine similarity of all pairs of corresponding row vectors,

where $\mathbf{D}\left(\mathbf{X}\right)=\mathrm{diag}\left(\frac{1}{\|\mathbf{x}^{(1)}\|_{2}},\cdots,\frac{1}{\|\mathbf{x}^{(N)}\|_{2}}\right)$ is the row-wise scaling matrix and

Note that, the form in Eqn. (9) is significantly different from tensor discriminant analysis (Jia et al.,, 2014; Tao et al.,, 2007), which integrates the actual label information as a regularizer and is also different from graph regularized tensor decomposition (Cai et al.,, 2010), which incorporates side information, such as geometrical positions (Maki et al.,, 2018). Compared to standard noise contrastive estimation (NCE) (Gutmann and Hyvärinen,, 2010; Chen et al.,, 2020) in the area of contrastive SSL, our SSL form in Eqn. (9) considers a subtraction form instead of the softmax formulation, making it amenable to quadratic optimization, as we will show in Sec. 3.3.

According to Eqn. (8), the self supervised loss ${\mathcal{L}}_{ss}$ is bounded by ${\mathcal{L}}^{\Theta}_{ss}(\gamma)$, while the constants (i.e., $C_{1},C_{2}$) can be absorbed into a weight hyperparameter $\beta$. We let the empirical self-supervised loss, $\tilde{\mathcal{L}}_{ss}{=}\beta\tilde{\mathcal{L}}^{\Theta}_{ss}(\gamma)$. Our objective follows both the fitness (i.e., CPD reconstruction loss) and alignment (i.e., self-supervised loss) principles, while also considering Tikhonov regularization (Golub and Von Matt,, 1997) to constrain the scale of all parameters,

where

The objective has (i) three hyperparameters, $\gamma,\alpha,\beta>0$; (ii) five factor matrices, $\mathbf{A},\mathbf{B},\mathbf{C}$, $\mathbf{X}$, $\tilde{\mathbf{X}}$.

We use the subproblem of $\mathbf{X}$ as an example. Given $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$, $\tilde{\mathbf{X}}$ fixed, we want to minimize Eqn. (10) by finding the optimal solution, denotd as $\mathbf{X}^{*}$,

• •

  First, we are interested to find that the matrix-formed problem in Eqn. (12) can be decomposed into row-wise subproblems and to obtain the solution of Eqn. (12), it is suffice to solve each subproblem independently. Let us consider the subproblem of the $k$-th row, which is

  $\displaystyle\operatorname*{arg\,min}_{\mathbf{x}}\left(\left\|{\mathcal{T}}^{(k)}-\llbracket{\mathbf{x}},\mathbf{A},\mathbf{B},\mathbf{C}\rrbracket\right\|_{F}^{2}+\alpha\|{\mathbf{x}}\|_{F}^{2}+\beta\mbox{Tr}\left(\frac{\mathbf{x}^{\top}}{\|\mathbf{x}\|_{2}}\mathbf{g}^{(k)}\mathbf{D}(\tilde{\mathbf{X}})\tilde{\mathbf{X}}\right)\right),$

  (13)

  where ${\mathcal{T}}^{(k)}$ is the $k$-th slice of ${\mathcal{T}}$, and $\mathbf{g}^{(k)}$ is the $k$-th row of $\mathbf{G}(\gamma)$.

• •

  Here, Eqn. (13) is still non-convex. We let the derivative of Eqn. (13) objective to be zero and arrange the terms, which yields,

  $\displaystyle{\mathbf{x}}$

  $\displaystyle=\mathbf{v}_{1}\mathbf{V}_{3}-\frac{\beta\mathbf{v}_{2}}{2\|\mathbf{x}\|_{2}}\left(\mathbf{I}-\frac{\mathbf{x}^{\top}\mathbf{x}}{\|\mathbf{x}\|^{2}_{2}}\right)\mathbf{V}_{3},$

  (14)

  where

  $\displaystyle\mathbf{v}_{1}=\mathbf{T}^{(k)}_{1}(\mathbf{A}\odot\mathbf{B}\odot\mathbf{C}),~{}~{}~{}~{}\mathbf{v}_{2}=\mathbf{g}^{(k)}\mathbf{D}(\tilde{\mathbf{X}})\tilde{\mathbf{X}},~{}~{}~{}~{}\mathbf{V}_{3}=\left(\mathbf{A}^{\top}\mathbf{A}*\mathbf{B}^{\top}\mathbf{B}*\mathbf{C}^{\top}\mathbf{C}+\alpha\mathbf{I}\right)^{-1}.$

  Here $\mathbf{x},\mathbf{v}_{1},\mathbf{v}_{2}$ are row vectors and $\mathbf{V}_{3}$ is a matrix. $\mathbf{T}^{(k)}_{1}$ is the $1$-mode unfolding of $\mathcal{T}^{(k)}$, $\odot$ is the Khatri-Rao product and $*$ is the Hadamard product (i.e., element-wise product).

• •

  We consider the following iterative rule and the fixed point is a solution for Eqn. (14), which is a stationary point of Eqn. (13),

  $\displaystyle{\mathbf{x}}_{\mathrm{impr}}$

  $\displaystyle=\mathbf{v}_{1}\mathbf{V}_{3}-\frac{\beta\mathbf{v}_{2}}{2\|\mathbf{x}_{\mathrm{init}}\|_{2}}\left(\mathbf{I}-\frac{\mathbf{x}^{\top}_{\mathrm{init}}\mathbf{x}_{\mathrm{init}}}{\|\mathbf{x}_{\mathrm{init}}\|^{2}_{2}}\right)\mathbf{V}_{3}.$

  (15)

  We use an initial guess ${\mathbf{x}}_{\mathrm{init}}$ (obtained by solving Eqn. (13) with while $\beta=0$, which is a least square problem) to start and then we repeat Eqn. (15) for each row (i.e., each $k$) and let the improved guess be the initial guess, ${\mathbf{x}}_{\mathrm{init}}\leftarrow{\mathbf{x}}_{\mathrm{impr}}$, to iteratively improve the result.

Theorem 1 (proof in appendix A) ensures that Eqn. (15) converges linearly if $\beta$ is chosen to be sufficiently small. In appendix B, we verify the liner convergence and also empirically show that one round of Eqn. (15) is sufficient in our experiment, where $\beta=2$. The non-convex subproblem of $\tilde{\mathbf{X}}$ can be solved in the same way.

To minimize Eqn. (11), we alternatively update $\mathbf{A},\mathbf{B},\mathbf{C}$, $\mathbf{X}$, and $\tilde{\mathbf{X}}$, where each subproblem involves only solving least squares problems with closed-form solutions. With large-scale tensors (as in the experiments), we optimizes the factors in mini-batches. Between mini-batches, the basis factors $\mathbf{A,B,C}$ are shared and updated incrementally. We summarize the algortihms in appendix C. The computation head of the algorithm is matricized tensor times Khatri-Rao product (MTTKRP). The complexity of our optimization algorithm is asymptotically the same as applying CP-ALS, which costs $O(NIJKR)$ to sweep over the whole tensor once.

We use four real-world datasets: (i) Sleep-EDF (Kemp et al.,, 2000), which contains EOG, EMG and EEG Polysomnography recordings; (ii) human activity recognition (HAR) (Anguita et al.,, 2013) with smartphone accelerometer and gyroscope data; (iii) Physikalisch Technische Bundesanstalt large scale cardiology database (PTB-XL) (Alday et al.,, 2020) with 12-lead ECG signals; (iv) Massachusetts General Hospital (MGH) (Biswal et al.,, 2018) datasets with multi-channel EEG waves. All datasets are split into three disjoint sets (i.e., unlabeled, training and test) by subjects, while training and test sets have labels. Basic statistics are shown in Table 1. All models (baselines and our ATD) use the same augmentation techniques: (a) jittering, (b) bandpass filtering, and (c) 3D position rotation. We provide an ablation study on the augmentation methods in appendix D.5.

We include the following comparison models from different perspectives:

• •

  Tensor based models: ${\bf\texttt{ATD}}_{ss-}$ is our variant, which removes the self-supervised loss from the objective in Eqn. (10); Stochastic alternating least squares (SALS) applies on the the CPD objective with Tikhonov regularizer, which works on large tensors; Graph regularized SALS (GR-SALS) augments the objective of SALS with a graph regularizer (Maki et al.,, 2018; Cai et al.,, 2010), define as $\mbox{Tr}\left(\mathbf{X}^{\top}\mathbf{G}\mathbf{X}\right)$.

• •

  Self-supervised models: SimCLR-$r$ (Chen et al.,, 2020) and BYOL-$r$ (Grill et al.,, 2020) are two popular SSL models with their own objective functions, where $r$ indicates the size of the output representation.

• •

  Auto-encoder models: AE-$r$ denotes a CNN based autoencoder with mean square error (MSE) reconstruction loss, and $\mbox{AE}_{ss}$-$r$ denotes the same autoencoder model with standard NCE loss in the bottleneck layer, where $r$ is the representation size.

All models use the unlabeled set to train a feature encoder and use training and test sets to evaluate. Note that, deep neural network models use the same CNN backbone. In appendix D.8, we have also compared with two recent supervised tensor learning models, which shows the usefulness of our ATD and the large unlabeled set, especially in low-label rate scenarios.

We evaluate model performance mainly based on classification accuracy, where we train an additional logistic classifier (He et al.,, 2020) on top of the feature encoder. Also, for different models, we compare their number of learnable parameters. The experiments are implemented by Python 3.8.5, Torch 1.8.0+cu111 on a Linux workstation with 256 GB memory, 32 core CPUs (3.70 GHz, 128 MB cache), two RTX 3090 GPUs (24 GB memory each). All training is performed on the GPU. For tensor based models, we use $R=32$ and implement the pipeline in CUDA manually, instead of using torch-autograd.