Multiscale Analysis for Higher-order Tensors

TL;DR

This paper introduces a novel multi-scale tensor decomposition method that adaptively partitions higher-order tensors into low-rank sub-tensors, improving data representation and classification accuracy for complex, multimodal datasets.

Contribution

It proposes a new multi-scale higher-order SVD approach inspired by subspace clustering, enhancing tensor approximation beyond traditional low-rank methods.

Findings

Effective dimensionality reduction demonstrated on real-life tensor data

Improved classification performance with the proposed method

Promising results in handling non-linear, multimodal datasets

Abstract

The widespread use of multisensor technology and the emergence of big datasets have created the need to develop tools to reduce, approximate, and classify large and multimodal data such as higher-order tensors. While early approaches focused on matrix and vector based methods to represent these higher-order data, more recently it has been shown that tensor decomposition methods are better equipped to capture couplings across their different modes. For these reasons, tensor decomposition methods have found applications in many different signal processing problems including dimensionality reduction, signal separation, linear regression, feature extraction, and classification. However, most of the existing tensor decomposition methods are based on the principle of finding a low-rank approximation in a linear subspace structure, where the definition of the rank may change depending on the particular decomposition. Since many datasets are not necessarily low-rank in a linear subspace, this often results in high approximation errors or low compression rates. In this paper, we introduce a new adaptive, multi-scale tensor decomposition method for higher order data inspired by hybrid linear modeling and subspace clustering techniques. In particular, we develop a multi-scale higher-order singular value decomposition (MS-HoSVD) approach where a given tensor is first permuted and then partitioned into several sub-tensors each of which can be represented as a low-rank tensor with increased representational efficiency. The proposed approach is evaluated for dimensionality reduction and classification for several different real-life tensor signals with promising results.

Figures

Tables

Table 1. TABLE I: Reconstruction error and compression rate computed for pruned tree structure obtained by applying MS-HoSVD with 2 scales to PIE data.

$\lambda$	0	0.22	0.25	0.30	0.75
Normalized error	0.0276	0.0395	0.0506	0.0530	0.0540
Compression	0.1241	0.0809	0.0377	0.0284	0.0261
Scales of subtensors	0+1+2	0+1+2	0+1+2	0+1+2	0+1

Table 2. TABLE II: Reconstruction error and compression rate computed for pruned tree structure obtained by applying MS-HoSVD with 2 scales to COIL-100 dataset.

$\lambda$	0	0.25	0.50	0.75	0.80
Normalized error	0.0857	0.0867	0.0913	0.1060	0.1207
Compression	0.0863	0.0840	0.0734	0.0526	0.0347
Scales of subtensors	0+1+2	0+1+2	0+1+2	0+1+2	0+1+2

Table 3. TABLE III: Reconstruction error and compression rate computed for pruned tree structure obtained by applying MS-HoSVD with 2 scales to Hand Gesture dataset.

$\lambda$	0	0.25	0.26	0.27	0.28
Normalized error	0.0691	0.0869	0.0913	0.0946	0.0999
Compression	0.1694	0.1056	0.0827	0.0698	0.0514
Scales of subtensors	0+1+2	0+1+2	0+1+2	0+1+2	0+1

Table 4. TABLE IV: Classification results for COIL-100 dataset over 20 trials with N f = 100 subscript 𝑁 𝑓 100 N_{f}=100 .

Training Size	Method	1-NN	Adaboost	Bayes
		mean $\pm$ std.	mean $\pm$ std.	mean $\pm$ std.
25$\%$	MS-HoSVD	93.71 $\pm$ 1.28	88.90 $\pm$ 1.24	89.88 $\pm$ 2.08
	HoSVD	93.07 $\pm$ 1.33	87.47 $\pm$ 1.53	87.36 $\pm$ 3.32
	T-Train	92.29 $\pm$ 2.23	87.26 $\pm$ 1.84	88.43 $\pm$ 1.36
50$\%$	MS-HoSVD	97.41 $\pm$ 0.69	92.54 $\pm$ 1.58	91.62 $\pm$ 1.21
	HoSVD	97.10 $\pm$ 0.83	91.29 $\pm$ 2.16	90.43 $\pm$ 1.60
	T-Train	96.99 $\pm$ 1.09	91.36 $\pm$ 1.44	91.21 $\pm$ 1.66
75$\%$	MS-HoSVD	98.38 $\pm$ 0.65	93.66 $\pm$1.35	92.60 $\pm$ 1.79
	HoSVD	98.16 $\pm$ 0.72	92.23 $\pm$1.46	92.38 $\pm$ 1.92
	T-Train	98.25 $\pm$ 0.41	93.10 $\pm$1.33	92.33 $\pm$ 1.52

Table 5. TABLE V: Classification results for hand gesture dataset over 20 trials with N f = 200 subscript 𝑁 𝑓 200 N_{f}=200 .

Training Size	Method	1-NN	Adaboost	Bayes
		mean $\pm$ std.	mean $\pm$ std.	mean $\pm$ std.
25$\%$	MS-HoSVD	75.40 $\pm$ 3.87	77.81 $\pm$ 1.87	82.89 $\pm$ 1.84
	HoSVD	75.01 $\pm$ 3.99	77.37 $\pm$ 2.51	83.11 $\pm$ 2.89
	T-Train	69.20 $\pm$ 2.63	67.53 $\pm$ 2.27	72.99 $\pm$ 4.91
50$\%$	MS-HoSVD	83.86 $\pm$ 3.12	85.46 $\pm$ 1.66	85.86 $\pm$ 2.25
	HoSVD	83.15 $\pm$ 2.90	85.14 $\pm$ 1.45	84.00 $\pm$ 1.68
	T-Train	78.97 $\pm$ 2.25	78.82 $\pm$ 2.87	81.46 $\pm$ 1.37
75$\%$	MS-HoSVD	87.47 $\pm$ 2.07	88.15 $\pm$ 2.20	86.75 $\pm$ 2.63
	HoSVD	86.93 $\pm$ 2.31	88.08 $\pm$ 2.51	85.04 $\pm$ 2.41
	T-Train	85.64 $\pm$ 2.57	82.60 $\pm$ 3.16	83.64 $\pm$ 2.31

Table 6. TABLE VI: Mean and standard deviation for reconstruction error of low-rank approximations of 𝒳 𝒳 \mathcal{X} over 20 trials.

	Reconstruction Error
0th scale rank	$(2,\mathrm{\hspace{0.33em}2},\mathrm{\hspace{0.33em}2})$	$(2,\mathrm{\hspace{0.33em}2},\mathrm{\hspace{0.33em}2})$	$(2,\mathrm{\hspace{0.33em}2},\mathrm{\hspace{0.33em}2})$
1st scale rank	$(2,\mathrm{\hspace{0.33em}2},\mathrm{\hspace{0.33em}2})$	$(4,\mathrm{\hspace{0.33em}4},\mathrm{\hspace{0.33em}4})$	$(6,\mathrm{\hspace{0.33em}6},\mathrm{\hspace{0.33em}6})$
Ground Truth	0.2502	0.0304	-
Partitioning	$\pm$ 0.0263	$\pm$ 0.0077	-
Clustering	0.3587	0.1099	0.0254
by LSA	$\pm$ 0.0583	$\pm$ 0.0391	$\pm$ 0.0152
Random	0.6095	0.3588	0.1855
Partitioning	$\pm$ 0.0398	$\pm$ 0.0298	$\pm$ 0.0251
rank	$(4,\mathrm{\hspace{0.33em}4},\mathrm{\hspace{0.33em}4})$	$(8,\mathrm{\hspace{0.33em}8},\mathrm{\hspace{0.33em}8})$	$(12,\mathrm{\hspace{0.33em}12},\mathrm{\hspace{0.33em}12})$
truncated	0.5457	0.2127	0.0733
HoSVD	$\pm$0.0449	$\pm$ 0.0195	$\pm$ 0.0129

Table 7. TABLE VII: Computed reconstruction error (mean and std.) corresponding to Theorem 1 for simulated data over 20 trials.

LHS			RHS
0th scale rank	$(2,\mathrm{\hspace{0.33em}2},\mathrm{\hspace{0.33em}2})$		rank
1st scale rank	$(2,\mathrm{\hspace{0.33em}2},\mathrm{\hspace{0.33em}2})$		$(4,\mathrm{\hspace{0.33em}4},\mathrm{\hspace{0.33em}4})$	$(6,\mathrm{\hspace{0.33em}6},\mathrm{\hspace{0.33em}6})$	$(8,\mathrm{\hspace{0.33em}8},\mathrm{\hspace{0.33em}8})$
Ground Truth	0.2618
Partitioning	$\pm$ 0.0236		2.9928	1.1835	0.4228
Clustering	0.3469		$\pm$0.5081	$\pm$0.3010	$\pm$0.0978
by LSA	$\pm$ 0.0513

Equations

$${\mathbf{Y}}_{(n)}={\mathbf{U}}^{(n)}{\mathbf{X}}_{(n)}({\mathbf{U}}^{(N)}\otimes ...{\mathbf{U}}^{(n+1)}\otimes {\mathbf{U}}^{(n-1)}...\otimes {\mathbf{U}}^{(1)}{)}^{\top },$$

$$n\text{mbox}-rank(\mathcal{X})=\left(\text{mbox}rank({\mathbf{X}}_{(1)}),\text{mbox}rank({\mathbf{X}}_{(2)}),...,\text{mbox}rank({\mathbf{X}}_{(N)})\right).$$

$$⟨\mathcal{X},\mathcal{Y}⟩=\sum _{i_1=1}^{I_1}\sum _{i_2=1}^{I_2}...\sum _{i_N=1}^{I_N}{x}_{i_1,i_2,...,i_N}{y}_{i_1,i_2,...,i_N}.$$

$$\parallel \mathcal{X}\parallel =\sqrt{⟨\mathcal{X},\mathcal{X}⟩}=\sqrt{\sum _{i_1=1}^{I_1}\sum _{i_2=1}^{I_2}...\sum _{i_N=1}^{I_N}{x}_{i_1,i_2,...,i_N}^2}.$$

$$\mathcal{X}=\mathcal{X}{\times }_n\mathbf{I}=\mathcal{X}{\times }_n((\mathbf{I}-\mathbf{P})+\mathbf{P})=\mathcal{X}{\times }_n(\mathbf{I}-\mathbf{P})+\mathcal{X}{\times }_n\mathbf{P}.$$

$$⟨\mathcal{X}{\times }_n(\mathbf{I}-\mathbf{P}),\mathcal{X}{\times }_n\mathbf{P}⟩={⟨(\mathbf{I}-\mathbf{P}){\mathbf{X}}_{(n)},\mathbf{P}{\mathbf{X}}_{(n)}⟩}_{\mathrm{F}}=\mathrm{Trace}\left({\mathbf{X}}_{(n)}^{\top }(\mathbf{I}-\mathbf{P})\mathbf{P}{\mathbf{X}}_{(n)}\right)=0.$$

$$\mathcal{X}-\mathcal{X}{\times }_n\mathbf{P}=\mathcal{X}{\times }_n(\mathbf{I}-\mathbf{P})$$

$${\Vert \mathcal{X}-\mathcal{X}{\times }_1{\mathbf{U}}^{(1)}{\times }_2{\mathbf{U}}^{(2)}...{\times }_N{\mathbf{U}}^{(N)}\Vert }^2=:{\Vert \mathcal{X}-\mathcal{X}{\text{bigtimes}}_{n=1}^N{\mathbf{U}}^{(n)}\Vert }^2=\sum _{n=1}^N{\Vert \mathcal{X}{\text{bigtimes}}_{h=1}^{n-1}{\mathbf{U}}^{(h)}{\times }_n\left(\mathbf{I}-{\mathbf{U}}^{(n)}\right)\Vert }^2.$$

$$\mathcal{X}=\mathcal{C}{\times }_1{\mathbf{U}}^{(1)}{\times }_2{\mathbf{U}}^{(2)}...{\times }_N{\mathbf{U}}^{(N)}=\mathcal{C}{\text{bigtimes}}_{n=1}^N{\mathbf{U}}^{(n)}$$

$$\mathcal{C}=\mathcal{X}{\times }_1{\left({\mathbf{U}}^{(1)}\right)}^{\top }{\times }_2{\left({\mathbf{U}}^{(2)}\right)}^{\top }...{\times }_N{\left({\mathbf{U}}^{(N)}\right)}^{\top }.$$

$$⟨{\mathcal{C}}_{i_n=\alpha },{\mathcal{C}}_{i_n=\beta }⟩=0when\alpha \ne \beta .$$

$$\parallel {\mathcal{C}}_{i_n=1}\parallel \ge \parallel {\mathcal{C}}_{i_n=2}\parallel \ge ...\ge \parallel {\mathcal{C}}_{i_n=I_n}\parallel \ge 0$$

$$\mathcal{X}=\mathcal{C}{\times }_1{\mathbf{U}}^{(1)}{\times }_2{\mathbf{U}}^{(2)}...{\times }_N{\mathbf{U}}^{(N)},$$

$${\hat{\mathcal{X}}}_0={\mathcal{C}}_0{\times }_1{\hat{\mathbf{U}}}^{(1)}{\times }_2{\hat{\mathbf{U}}}^{(2)}...{\times }_N{\hat{\mathbf{U}}}^{(N)}$$

$$r_n=\arg \underset{i}{\min }\sum _{l=1}^i{\sigma }_l^{(n)}s.t.\frac{{\sum }_{l=1}^i{\sigma }_l^{(n)}}{{\sum }_{l=1}^{I_n}{\sigma }_l^{(n)}}\ge \tau ,$$

$$\mathcal{X}={\hat{\mathcal{X}}}_0+{\mathcal{W}}_0,$$

$$\begin{array}[]{c}\mathcal{X}_{1,k}(i_{1},\;i_{2},...,\;i_{N})=\mathcal{W}_{0}(J_{1,k}^{1}(i_{1}),\;J_{1,k}^{2}(i_{2}),...,\;J_{1,k}^{N}(i_{N})),\\ \mathcal{X}_{1,k}=\mathcal{W}_{0}(J_{1,k}^{1}\times J_{1,k}^{2}\times...\times J_{1,k}^{N}),\end{array}$$

$${\hat{\mathcal{X}}}_{1,k}={\mathcal{C}}_{1,k}{\times }_1{\hat{\mathbf{U}}}_{1,k}^{(1)}{\times }_2{\hat{\mathbf{U}}}_{1,k}^{(2)}...{\times }_N{\hat{\mathbf{U}}}_{1,k}^{(N)},$$

$${\hat{\mathcal{X}}}_1(J_{1,k}^1\times J_{1,k}^2\times ...\times J_{1,k}^n)={\hat{\mathcal{X}}}_{1,k}.$$

$${\mathcal{W}}_1(J_{1,k}^1\times J_{1,k}^2\times ...\times J_{1,k}^n)={\mathcal{W}}_{1,k},$$

$$\mathcal{X}={\hat{\mathcal{X}}}_0+{\mathcal{W}}_0={\hat{\mathcal{X}}}_0+{\hat{\mathcal{X}}}_1+{\mathcal{W}}_1.$$

$$\mathcal{X}=\sum _{i=0}^j{\hat{\mathcal{X}}}_i+{\mathcal{W}}_j.$$

$$\hat{\mathcal{X}}={\mathcal{C}}_H{\times }_1{\mathbf{U}}_H^{(1)}{\times }_2{\mathbf{U}}_H^{(2)}...{\times }_N{\mathbf{U}}_H^{(N)},$$

$${\mathbf{R}}_k^{(n)}(i,j)=\{\begin{array}{ll}1,& \text{if }i=j,\text{mbox}andj\in J_{1,k}^n\\ 0,& \text{otherwise}\end{array}$$

$$\mathcal{X}{|}_k:={\mathcal{W}}_0{\text{bigtimes}}_{n=1}^N{\mathbf{R}}_k^{(n)}={\mathcal{W}}_0{\times }_1{\mathbf{R}}_k^{(1)}{\times }_2{\mathbf{R}}_k^{(2)}...{\times }_N{\mathbf{R}}_k^{(N)}.$$

$${\mathcal{W}}_0=\sum _{k=1}^K\mathcal{X}{|}_k$$

$$⟨\mathcal{X}{|}_k,\mathcal{X}{|}_j⟩=0\text{for all}j,k\in [K]\text{with}j\ne k$$

$${\mathbf{Q}}_k^{(n)}:={\hat{\mathbf{U}}}_k^{(n)}{\left({\hat{\mathbf{U}}}_k^{(n)}\right)}^{\top }.$$

$${\mathbf{Q}}_k^{(n)}:={\hat{\mathbf{U}}}_k^{(n)}{\left({\hat{\mathbf{U}}}_k^{(n)}\right)}^{\top }={\mathbf{R}}_k^{(n)}{\hat{\mathbf{U}}}_k^{(n)}{\left({\mathbf{R}}_k^{(n)}{\hat{\mathbf{U}}}_k^{(n)}\right)}^{\top }={\mathbf{R}}_k^{(n)}{\hat{\mathbf{U}}}_k^{(n)}{\left({\hat{\mathbf{U}}}_k^{(n)}\right)}^{\top }{\left({\mathbf{R}}_k^{(n)}\right)}^{\top }={\mathbf{R}}_k^{(n)}{\mathbf{Q}}_k^{(n)}{\mathbf{R}}_k^{(n)}$$

$${\mathbf{R}}_k^{(n)}{\mathbf{Q}}_k^{(n)}={\mathbf{R}}_k^{(n)}\left({\mathbf{R}}_k^{(n)}{\mathbf{Q}}_k^{(n)}{\mathbf{R}}_k^{(n)}\right)={\mathbf{R}}_k^{(n)}{\mathbf{Q}}_k^{(n)}{\mathbf{R}}_k^{(n)}={\mathbf{Q}}_k^{(n)}={\mathbf{R}}_k^{(n)}{\mathbf{Q}}_k^{(n)}{\mathbf{R}}_k^{(n)}=\left({\mathbf{R}}_k^{(n)}{\mathbf{Q}}_k^{(n)}{\mathbf{R}}_k^{(n)}\right){\mathbf{R}}_k^{(n)}={\mathbf{Q}}_k^{(n)}{\mathbf{R}}_k^{(n)}$$

Full text

I Introduction

Data in the form of multidimensional arrays, also referred to as tensors, arise in a variety of applications including chemometrics, hyperspectral imaging, high resolution videos, neuroimaging, biometrics and social network analysis [1, 2, 3]. These applications produce massive amounts of data collected in various forms with multiple aspects and high dimensionality. Tensors, which are multi-dimensional generalizations of matrices, provide a useful representation for such data. A crucial step in many applications involving higher-orders tensors is multiway reduction of the data to ensure that the reduced representation of the tensor retains certain characteristics. Early multiway data analysis approaches reformatted the tensor data as a matrix and resorted to methods developed for classical two-way analysis. However, one cannot discover hidden components within multiway data using conventional matrix decomposition methods as matrix based representations cannot capture multiway couplings focusing on standard pairwise interactions. To this end, many different types of tensor decomposition methods have been proposed in literature [4, 5, 6, 7, 8, 9, 10].

In contrast to the matrix case where data reduction is often accomplished via low-rank representations such as singular value decomposition (SVD), the notion of rank for higher order tensors is not uniquely defined. The CANDECOMP/PARAFAC (CP) and Tucker decompositions are two of the most widely used tensor decomposition methods for data reduction [11, 12]. For CP, the goal is to approximate the given tensor as a weighted sum of rank-1 tensors, where a rank-1 tensor refers to the outer product of $n$ vectors with $n$ being equal to the order of the tensor. The Tucker model allows for interactions between the factors from different modes resulting in a typically dense, but small, core tensor. This model also introduces the notion of Tucker rank or n-rank, which refers to the n-tuple of ranks corresponding to the tensor’s unfoldings along each of its modes. Therefore, low rank approximation with the Tucker model can be obtained by projections onto low-rank factor matrices. Unlike the CP decomposition, the Tucker decomposition is in general non-unique. To help obtain meaningful and unique representations by the Tucker decomposition, orthogonality, sparsity, and non-negativity constraints are often imposed on the factors yielding, e.g., the Non-Negative Tensor Factorization (NTF) and the Sparse Non-Negative Tucker Decomposition [13, 14, 15]. The Tucker decomposition with orthogonality constraints on the factors is known as the Higher-Order Singular Value Decomposition (HoSVD), or Multilinear SVD [12]. The HoSVD can be computed by simply flattening the tensor in each mode and calculating the n-mode singular vectors corresponding to that mode.

With the emergence of multidimensional big data, classical tensor representation and decomposition methods have become inadequate since the size of these tensors exceeds available working memory and the processing time is very long. In order to address the problem of large-scale tensor decomposition, several block-wise tensor decomposition methods have been proposed [6]. The basic idea is to partition a big data tensor into smaller blocks and perform tensor related operations block-wise using a suitable tensor format. Preliminary approaches relied on a hierarchical tree structure and reduced the storage of d-dimensional arrays to the storage of auxiliary three-dimensional ones such as the tensor-train decomposition (T-Train), also known as the matrix product state (MPS) decomposition, [5] and the Hierarchical Tucker Decomposition (H-Tucker) [16]. In particular, in the area of large volumetric data visualization, tensor based multiresolution hierarchical methods such as TAMRESH have attracted attention [17]. However, all of these methods are interested in fitting a low-rank model to data which lies near a linear subspace, thus being limited to learning linear structure.

Similar to the research efforts in tensor reduction, low-dimensional subspace and manifold learning methods have also been extended for higher order data clustering and classification applications. In early work in the area, Vasilescu and Terzopoulos [18] extended the eigenface concept to the tensorface by using higher order SVD and taking different modes such as expression, illumination and pose into account. Similarly, 2D-PCA for matrices has been used for feature extraction from face images without converting the images into vectors [19]. He et al. [20] extended locality preserving projections [21] to second order tensors for face recognition. Dai and Yeung [22] presented generalized tensor embedding methods such as the extensions of local discriminant embedding methods [23], neighborhood preserving embedding methods [24], and locality preserving projection methods [21] to tensors. Li et al. [25] proposed a supervised manifold learning method for vector type data which preserves local structures in each class of samples, and then extended the algorithm to tensors to provide improved performance for face and gait recognition. Similar to vector-type manifold learning algorithms, the aim of these methods is to find an optimal linear transformation for the tensor-type training data samples without vectorizing them and mapping these samples to a low dimensional subspace while preserving the neighborhood information.

In this paper, we propose a novel multi-scale analysis technique to efficiently approximate tensor type data using locally linear low-rank approximations. The proposed method consists of two major steps: 1) Constructing a tree structure by partitioning the tensor into a collection of permuted subtensors, followed by 2) Constructing multiscale dictionaries by applying HoSVD to each subtensor. The contributions of the proposed framework and the novelty in the proposed approach with respect to previously published work in [26, 27] are manifold. They include: 1) The introduction of a more flexible multi-scale tensor decomposition method which allows the user to approximate a given tensor within given memory and processing power constraints; 2) the introduction of theoretical error bounds for the proposed decomposition; 3) the introduction of adaptive pruning to achieve a better trade-off between compression rate and reconstruction error for the developed factorizations; 4) the extensive evaluation of the method for both data reduction and classification applications; and 5) a detailed comparison of the proposed method to state-of-the-art tensor decomposition methods including the HoSVD, T-Train, and H-Tucker decompositions.

The remainder of the paper is organized as follows. In Section II, basic notation and tensor operations are reviewed. The proposed multiscale tensor decomposition method along with theoretical error bounds are then introduced in Section III. Sections IV and V illustrate the results of applying the proposed framework to data reduction and classification problems, respectively.

II Background

II-A Tensor Notation and Algebra

A multidimensional array with $N$ modes $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times...\times I_{N}}$ is called a tensor, where $x_{i_{1},i_{2},..i_{N}}$ denotes the $({i_{1},i_{2},..i_{N}})^{th}$ element of the tensor $\mathcal{X}$. The vectors in $\mathbb{R}^{I_{n}}$ obtained by fixing all of the indices of such a tensor $\mathcal{X}$ except for the one that corresponds to its $n$th mode are called its mode-$n$ fibers. Let $[N]:=\{1,\dots,N\}$ for all $N\in\mathbb{N}$. Basic tensor operations are reviewed below (see, e.g., [11], [28], [29]).

Tensor addition and multiplication by a scalar: Two tensors $\mathcal{X},\mathcal{Y}\in\mathbb{R}^{I_{1}\times I_{2}\times...\times I_{N}}$ can be added using component-wise tensor addition. The resulting tensor $\mathcal{X}+\mathcal{Y}\in\mathbb{R}^{I_{1}\times I_{2}\times...\times I_{N}}$ has its entries given by $\left(\mathcal{X}+\mathcal{Y}\right)_{i_{1},i_{2},..i_{N}}=x_{i_{1},i_{2},..i_{N}}+y_{i_{1},i_{2},..i_{N}}$. Similarly, given a scalar $\alpha\in\mathbb{R}$ and a tensor $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times...\times I_{N}}$ the rescaled tensor $\alpha\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times...\times I_{N}}$ has its entries given by $\left(\alpha\mathcal{X}\right)_{i_{1},i_{2},..i_{N}}=\alpha~{}x_{i_{1},i_{2},..i_{N}}$.

Mode-$n$ products: The mode-$n$ product of a tensor $\mathcal{X}\in\mathbb{R}^{I_{1}\times...I_{n}\times...\times I_{N}}$ and a matrix ${\bf{U}}\in\mathbb{R}^{J\times I_{n}}$ is denoted as $\mathcal{Y}=\mathcal{X}\times_{n}{\bf{U}}$, $(\mathcal{Y})_{i_{1},i_{2},\ldots,i_{n-1},j,i_{n+1},\ldots,i_{N}}=\sum_{i_{n}=1}^{I_{n}}x_{i_{1},\ldots,i_{n},\ldots,i_{N}}u_{j,i_{n}}$. It is of size $I_{1}\times...\times I_{n-1}\times J\times I_{n+1}\times...\times I_{N}$. The following facts about mode-$n$ products are useful (see, e.g., [11],[29]).

Lemma 1**.**

Let $\mathcal{X},\mathcal{Y}\in\mathbb{R}^{I_{1}\times I_{2}\times...\times I_{N}}$, $\alpha,\beta\in\mathbb{R}$, and ${\bf U}^{(n)},{\bf V}^{(n)}\in\mathbb{R}^{J_{n}\times I_{n}}$ for all $n\in[N]$. The following are true:

1. ($a$)

$\left(\alpha\mathcal{X}+\beta\mathcal{Y}\right)\times_{n}{\bf U}^{(n)}=\alpha\left(\mathcal{X}\times_{n}{\bf U}^{(n)}\right)+\beta\left(\mathcal{Y}\times_{n}{\bf U}^{(n)}\right)$. 2. ($b$)

$\mathcal{X}\times_{n}\left(\alpha{\bf U}^{(n)}+\beta{\bf V}^{(n)}\right)=\alpha\left(\mathcal{X}\times_{n}{\bf U}^{(n)}\right)+\beta\left(\mathcal{X}\times_{n}{\bf V}^{(n)}\right)$. 3. ($c$)

If $n\neq m$ then $\mathcal{X}\times_{n}{\bf U}^{(n)}\times_{m}{\bf V}^{(m)}=\left(\mathcal{X}\times_{n}{\bf U}^{(n)}\right)\times_{m}{\bf V}^{(m)}=\left(\mathcal{X}\times_{m}{\bf V}^{(m)}\right)\times_{n}{\bf U}^{(n)}=\mathcal{X}\times_{m}{\bf V}^{(m)}\times_{n}{\bf U}^{(n)}$ . 4. ($d$)

If ${\bf W}\in\mathbb{C}^{P\times J_{n}}$ then $\mathcal{X}\times_{n}{\bf U}^{(n)}\times_{n}{\bf W}=\left(\mathcal{X}\times_{n}{\bf U}^{(n)}\right)\times_{n}{\bf W}=\mathcal{X}\times_{n}\left({\bf WU}^{(n)}\right)=\mathcal{X}\times_{n}{\bf WU}^{(n)}$.

Tensor matricization: The process of reordering the elements of the tensor into a matrix is known as matricization or unfolding. The mode-n matricization of a tensor $\mathcal{Y}\in\mathbb{R}^{I_{1}\times I_{2}\times...\times I_{N}}$ is denoted as ${\bf Y}_{(n)}\in\mathbb{R}^{I_{n}\times\prod_{m\neq n}I_{m}}$ and is obtained by arranging $\mathcal{Y}$’s mode-n fibers to be the columns of the resulting matrix. Unfolding the tensor $\mathcal{Y}=\mathcal{X}\times_{1}{\bf U}^{(1)}\times_{2}{\bf U}^{(2)}...\times_{N}{\bf U}^{(N)}=:{\displaystyle\mathcal{X}\bigtimes_{n=1}^{N}{\bf U}^{(n)}}$ along mode-$n$ is equivalent to

[TABLE]

where $\otimes$ is the matrix Kronecker product. In particular, (1) implies that the matricization $\left(\mathcal{X}\times_{n}{\bf U}^{(n)}\right)_{(n)}={\bf U}^{(n)}{\bf X}_{(n)}$.111Simply set ${\bf U}^{(m)}={\bf I}$ (the identity) for all $m\neq n$ in (1). This fact also easily follows directly from the definition of the mode-$n$ product.

It is worth noting that trivial inner product preserving isomorphisms exist between a tensor space $\mathbb{R}^{I_{1}\times I_{2}\times...\times I_{N}}$ and any of its matricized versions (i.e., mode-$n$ matricization can be viewed as an isomorphism between the original tensor vector space $\mathbb{R}^{I_{1}\times I_{2}\times...\times I_{N}}$ and its mode-$n$ matricized target vector space $\displaystyle\mathbb{R}^{I_{n}\times\prod_{m\neq n}I_{m}}$). In particular, the process of matricizing tensors is linear. If, for example, $\mathcal{X},\mathcal{Y}\in\mathbb{R}^{I_{1}\times I_{2}\times...\times I_{N}}$ then one can see that the mode-$n$ matricization of $\mathcal{X}+\mathcal{Y}\in\mathbb{R}^{I_{1}\times I_{2}\times...\times I_{N}}$ is $\left(\mathcal{X}+\mathcal{Y}\right)_{(n)}={\bf X}_{(n)}+{\bf Y}_{(n)}$ for all modes $n\in[N]$.

Tensor Rank: Unlike matrices, which have a unique definition of rank, there are multiple rank definitions for tensors including tensor rank and tensor n-rank. The rank of a tensor $\mathcal{X}\in\mathbb{R}^{I_{1}\times...I_{n}\times...\times I_{N}}$ is the smallest number of rank-one tensors that form $\mathcal{X}$ as their sum. The n-rank of $\mathcal{X}$ is the collection of ranks of unfoldings ${\bf X}_{(n)}$ and is denoted as:

[TABLE]

Tensor inner product: The inner product of two same sized tensors $\mathcal{X},\mathcal{Y}\in\mathbb{R}^{I_{1}\times I_{2}\times...\times I_{N}}$ is the sum of the products of their elements.

[TABLE]

It is not too difficult to see that matricization preserves Hilbert-Schmidt/Frobenius matrix inner products, i.e., $\left\langle\mathcal{X},\mathcal{Y}\right\rangle=\left\langle{\bf X}_{(n)},{\bf Y}_{(n)}\right\rangle_{\rm F}={\rm Trace}\left({\bf X}_{(n)}^{\top}{\bf Y}_{(n)}\right)$ holds for all $n\in[N]$. If $\left\langle\mathcal{X},\mathcal{Y}\right\rangle=0$, $\mathcal{X}$ and $\mathcal{Y}$ are orthogonal.

Tensor norm: Norm of a tensor $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times...\times I_{N}}$ is the square root of the sum of the squares of all its elements.

[TABLE]

The fact that matricization preserves Frobenius matrix inner products also means that it preserves Frobenius matrix norms. As a result we have that $\left\|\mathcal{X}\right\|=\left\|{\bf X}_{(n)}\right\|_{\rm F}$ holds for all $n\in[N]$. If $\mathcal{X}$ and $\mathcal{Y}$ are orthogonal and also have unit norm (i.e., have $\|\mathcal{X}\|=\|\mathcal{Y}\|=1$) we will say that they are an orthonormal pair.

II-B Some Useful Facts Concerning Mode-$n$ Products and Orthogonality

Let ${\bf I}\in\mathbb{R}^{I_{n}\times I_{n}}$ be the identity matrix. Given a (low rank) orthogonal projection matrix ${\bf P}\in\mathbb{R}^{I_{n}\times I_{n}}$ one can decompose any given tensor $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times...\times I_{N}}$ into two orthogonal tensors using Lemma 1 (b)

[TABLE]

To check that the last two summands are orthogonal one can use (1) to compute that

[TABLE]

As a result one can also verify that the Pythagorean theorem holds, i.e., that $\|\mathcal{X}\|^{2}=\|\mathcal{X}\times_{n}{\bf P}\|^{2}+\|\mathcal{X}\times_{n}({\bf I}-{\bf P})\|^{2}$.

If we now regard $\mathcal{X}\times_{n}{\bf P}$ as a low rank approximation to $\mathcal{X}$ then we can see that its approximation error

[TABLE]

is orthogonal to the low rank approximation $\mathcal{X}\times_{n}{\bf P}$, as one would expect. Furthermore, the norm of its approximation error satisfies $\|\mathcal{X}\times_{n}({\bf I}-{\bf P})\|^{2}=\|\mathcal{X}\|^{2}-\|\mathcal{X}\times_{n}{\bf P}\|^{2}$. By continuing to use similar ideas in combination with lemma 1 for all modes one can prove the following more general Pythagorean result (see, e.g., theorem 5.1 in [29]).

Lemma 2**.**

Let $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times...\times I_{N}}$ and ${\bf U}^{(n)}\in\mathbb{R}^{I_{n}\times I_{n}}$ be an orthogonal projection matrix for all $n\in[N]$. Then,

[TABLE]

II-C The Higher Order Singular Value Decomposition (HoSVD)

Any tensor $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times...\times I_{N}}$ can be decomposed as mode products of a core tensor $\mathcal{C}\in\mathbb{R}^{I_{1}\times I_{2}\times...\times I_{N}}$ with $N$ orthogonal matrices ${\bf U}^{(n)}\in\mathbb{R}^{I_{n}\times I_{n}}$ each of which is composed of the left singular vectors of ${\bf X}_{(n)}$ [12]:

[TABLE]

where $\mathcal{C}$ is computed as

[TABLE]

Let $\mathcal{C}_{i_{n}=\alpha}$ be a subtensor of $\mathcal{C}$ obtained by fixing the $n$th index to $\alpha$. This subtensor satisfies the following properties:

• •

all-orthogonality: $\mathcal{C}_{i_{n}=\alpha}$ and $\mathcal{C}_{i_{n}=\beta}$ are orthogonal for all possible values of $n$, $\alpha$ and $\beta$ subject to $\alpha\neq\beta$.

[TABLE]

• •

ordering:

[TABLE]

for $n\in\left[N\right]$.

III Multiscale Analysis of Higher-order Datasets

In this section, we present a new tensor decomposition method named Multiscale HoSVD (MS-HoSVD) for an $N$th order tensor, $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times...\times I_{N}}$. The proposed method recursively applies the following two-step approach: (i) Low-rank tensor approximation, followed by (ii) Partitioning the residual (original minus low-rank) tensor into subtensors.

A tensor $\mathcal{X}$ is first decomposed using HoSVD as follows:

[TABLE]

where the ${\bf U}^{(n)}$’s are the left singular vectors of the unfoldings ${\bf X}_{(n)}$. The low-rank approximation of $\mathcal{X}$ is obtained by

[TABLE]

where $\hat{\bf U}^{(n)}\in\mathbb{R}^{I_{n}\times r_{n}}$s are the truncated matrices obtained by keeping the first $r_{n}$ columns of ${\bf U}^{(n)}$ and $\mathcal{C}_{0}=\mathcal{X}\times_{1}\left(\hat{\bf U}^{(1)}\right)^{\top}\times_{2}\left(\hat{\bf U}^{(2)}\right)^{\top}...\times_{N}\left(\hat{\bf U}^{(N)}\right)^{\top}$. The multilinear-rank of $\hat{\mathcal{X}}_{0}$, $\left\{r_{1},...,\;r_{N}\right\}$, can either be given a priori, or an energy criterion can be used to determine the minimum number of singular values to keep along each mode as:

[TABLE]

where $\sigma^{(n)}_{l}$ is the $l$th singular value of the matrix obtained from the SVD of the unfolding ${\bf X}_{(n)}$, and $\tau$ is an energy threshold. Once $\hat{\mathcal{X}}_{0}$ is obtained, the tensor $\mathcal{X}$ can be written as

[TABLE]

where $\mathcal{W}_{0}$ is the residual tensor.

For the first-scale analysis, to better encode the details of $\mathcal{X}$, we adapted an idea similar to the one presented in [26, 30]. The $0^{\rm th}$-scale residual tensor, $\mathcal{W}_{0}$ is first decomposed into subtensors as follows. $\mathcal{W}_{0}\in\mathbb{R}^{I_{1}\times I_{2}\times...\times I_{N}}$ is unfolded across each mode yielding ${\bf W}_{0,(n)}\in\mathbb{R}^{I_{n}\times\prod_{j\neq n}I_{j}}$ whose columns are the mode-$n$ fibers of $\mathcal{W}_{0}$. For each mode, rows of ${\bf W}_{0,(n)}$ are partitioned into $c_{n}$ non-overlapping clusters using a clustering algorithm such as local subspace analysis (LSA) [31] in order to encourage the formation of new subtensors which are intrinsically lower rank, and therefore better approximated via a smaller HoSVD at the next scale. The Cartesian product of the partitioning labels coming from the $N$ modes yields $K=\prod_{i=1}^{N}c_{i}$ disjoint subtensors $\mathcal{X}_{1,k}$ where $k\in[K]$.

Let $J_{0}^{n}$ be the index set corresponding to the $n$th mode of $\mathcal{W}_{0}$ with $J_{0}^{n}=[I_{n}]$, and let $J_{1,k}^{n}$ be the index set of the subtensor $\mathcal{X}_{1,k}$ for the $n$th mode, where $J_{1,k}^{n}\subset J_{0}^{n}$ for all $k\in[K]$ and $n\in[N]$. Index sets of subtensors for the $n$th mode satisfy $\bigcup_{k=1}^{K}J_{1,k}^{n}=J_{0}^{n}$ for all $n\in[N]$. The $k$th subtensor $\mathcal{X}_{1,k}\in\mathbb{R}^{\left|J_{1,k}^{1}\right|\times\left|J_{1,k}^{2}\right|\times\dots\times\left|J_{1,k}^{N}\right|}$ is obtained by

[TABLE]

where $i_{n}\in\left[\left|J_{1,k}^{n}\right|\right]$. Low-rank approximation for each subtensor is obtained by applying HoSVD as:

[TABLE]

where $\mathcal{C}_{1,k}$ and $\hat{\bf U}^{(n)}_{1,k}\in\mathbb{R}^{|J^{n}_{1,k}|\times r^{(n)}_{1,k}}$s correspond to the core tensor and low-rank projection basis matrices of $\mathcal{X}_{1,k}$, respectively. We can then define $\hat{\mathcal{X}}_{1}$ as the $1^{\rm st}$-scale approximation of $\mathcal{X}$ formed by mapping all of the subtensors onto $\hat{\mathcal{X}}_{1,k}$ as follows:

[TABLE]

Similarly, $1^{\rm st}$ scale residual tensor is obtained by

[TABLE]

where $\mathcal{W}_{1,k}=\mathcal{X}_{1,k}-\hat{\mathcal{X}}_{1,k}$. Therefore, $\mathcal{X}$ can be rewritten as:

[TABLE]

Continuing in this fashion the $j^{th}$ scale approximation of $\mathcal{X}$ is obtained by partitioning $\mathcal{W}_{j-1,k}$s into subtensors $\mathcal{X}_{j,k}$s and fitting a low-rank model to each one of them in a similar fashion. Finally, the $j^{th}$ scale decomposition of $\mathcal{X}$ can be written as:

[TABLE]

Algorithm 1 describes the pseudo code for this approach and Fig. 1 illustrates 1-scale MS-HoSVD.

III-A Memory Cost of the First Scale Decomposition

Let $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times....\times I_{N}}$ be an $N$th order tensor. To simplify the notation, assume that the dimension of each mode is the same, i.e. $I_{1}=I_{2}=....=I_{N}=I$. Assume $\mathcal{X}$ is approximated by HoSVD as:

[TABLE]

by fixing the rank of each mode matrix as $\mbox{rank}({\bf U}^{(i)}_{H})=r_{H}$ for $i\in\left\{1,\;2,...,\;N\right\}$. Let $\mathbb{F}(\cdot)$ be a function that quantifies the memory cost, then the storage cost of ${\mathcal{X}}$ decomposed by HoSVD is $\mathbb{F}(\mathcal{C}_{H})+\sum_{i=1}^{N}(\mathbb{F}({\bf U}^{(i)}_{H}))\approx r_{H}^{N}+NIr_{H}$.

For multiscale analysis at scale 1, $\hat{\mathcal{X}}=\hat{\mathcal{X}}_{0}+\hat{\mathcal{X}}_{1}$. The cost of storing $\hat{\mathcal{X}}_{0}$ is $\mathbb{F}(\mathcal{C}_{0})+\sum_{i=1}^{N}(\mathbb{F}(\hat{\bf U}^{(i)}))\approx r_{0}^{N}+NIr_{0}$ where the rank of each mode matrix is fixed at $\mbox{rank}({\bf U}^{(i)})=r_{0}$ for $i\in\left\{1,\;2,...,\;N\right\}$. The cost of storing $\hat{\mathcal{X}}_{1}$ is the sum of the storage costs for each of the $K=\prod_{i=1}^{N}c(i)$ subtensors $\hat{\mathcal{X}}_{1,k}$. Assume $c(i)=c$ for all $i\in\left\{1,\;2,...,\;N\right\}$ yielding $c^{N}$ equally sized subtensors, and that each $\hat{\mathcal{X}}_{1,k}$ is decomposed using the HoSVD as $\hat{\mathcal{X}}_{1,k}=\mathcal{C}_{1,k}\times_{1}\hat{\bf U}^{(1)}_{1,k}\times_{2}\hat{\bf U}^{(2)}_{1,k}...\times_{N}\hat{\bf U}^{(N)}_{1,k}$. Let the rank of each mode matrix be fixed as $\mbox{rank}(\hat{\bf U}^{(i)}_{1,k})=r_{1}$ for all $i\in\left\{1,\;2,...,\;N\right\}$ and $k\in\left\{1,\;2,...,\;K\right\}$. Then, the memory cost for the first scale is $\sum_{k=1}^{K}\left(\mathbb{F}(\mathcal{C}_{1,k})+\sum_{i=1}^{N}\mathbb{F}(\hat{\bf U}^{(i)}_{1,k})\right)\approx c^{N}\left(r^{N}_{1}+\dfrac{NIr_{1}}{c}\right)$. Choosing $r_{1}\lesssim\dfrac{r_{0}}{c^{(N-1)}}$ ensures that the storage cost does not grow exponentially so that $\mathbb{F}(\hat{\mathcal{X}}_{1})<\mathbb{F}(\hat{\mathcal{X}}_{0})$ since the total cost becomes approximately equal to $r_{0}^{N}\left(1+\dfrac{1}{c^{N^{2}-2N}}\right)+2NIr_{0}$. Thus, picking $r_{0}\approx r_{H}/2$ can now provide lower storage cost for the first scale analysis than for HoSVD.

III-B Computational Complexity

The computational complexity of MS-HoSVD at the first scale is equal to the sum of computational complexity of computing HoSVD at the parent node, partitioning into subtensors and computing HoSVD for each one of the subtensors. Computational complexity of HoSVD of an N-way tensor $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times...\times I_{N}}$ where ${I_{1}=I_{2}=...=I_{N}=I}$ is $\mathcal{O}\left(NI^{(N+1)}\right)$ [32]. By assuming that the partitioning is performed using K-means (via Lloyd’s algorithm) with $c_{i}=c$ along each mode, the complexity partitioning along each mode is $\mathcal{O}\left(N{I}^{N}ci\right)$, where $i$ is the number of iterations used in Lloyd’s algorithm. Finally, the total complexity of applying the HoSVD to $c^{N}$ equally sized subtensors is $\mathcal{O}\left(c^{N}N{\left(I/c\right)}^{(N+1)}\right)$. Therefore, first scale MS-HoSVD has a total computational complexity of ${\mathcal{O}}\left(NI^{(N+1)}+N{I}^{N}ci+c^{N}N{\left(I/c\right)}^{(N+1)}\right)$. Note that this complexity is similar to that of the HoSVD whenever $ci$ is small compared to $I$. The runtime complexity of these multiscale methods can be reduced even further by computing the HoSVDs for different subtensors in parallel whenever possible, as well as by utilizing distributed and parallel SVD algorithms such as [33] when computing all the required HoSVD decompositions.

III-C A Linear Algebraic Representation of the Proposed Multiscale HoSVD Approach

Though the tree-based representation of the proposed MS-HoSVD approach used above in Algorithm 1 is useful for algorithmic development, it is somewhat less useful for theoretical error analysis. In this subsection we will develop formulas for the proposed MS-HoSVD approach which are more amenable to error analysis. In the process, we will also formulate a criterion which, when satisfied, guarantees that the proposed fist scale MS-HoSVD approach produces an accurate multiscale approximation to a given tensor.

Preliminaries: We can construct full size first-scale subtensors of the residual tensor $\mathcal{W}_{0}\in\mathbb{R}^{I_{1}\times I_{2}\times...I_{N}}$ from (12), $\mathcal{X}|_{k}\in\mathbb{R}^{I_{1}\times I_{2}\times...I_{N}}$ for all $k\in[K]$, using the index sets $J_{1,k}^{n}$ from (13) along with diagonal restriction matrices. Let ${\bf R}_{k}^{(n)}\in\left\{0,\;1\right\}^{I_{n}\times I_{n}}$ be the diagonal matrix with entries given by

[TABLE]

for all $k\in[K]$, and $n\in[N]$. We then define

[TABLE]

Thus, the $k$th subtensor $\mathcal{X}|_{k}$ will only have nonzero entries, given by $\mathcal{W}_{0}(J_{1,k}^{1}\times...\times J_{1,k}^{N})$, in the locations indexed by the sets $J_{1,k}^{n}$ from above. The properties of the index sets $J_{1,k}^{n}$ furthermore guarantee that these subtensors all have disjoint support. As a result both

[TABLE]

and

[TABLE]

will always hold.

Recall that we want to compute the HoSVD of the subtensors we form at each scale in order to create low-rank projection basis matrices along the lines of those in (14). Toward this end we compute the top $r^{(n)}_{k}\leq rank({\bf R}_{k}^{(n)})=|J^{n}_{1,k}|$ left singular vectors of the mode-n matricization of each $\mathcal{X}|_{k}$, ${\bf X|_{k}}_{(n)}\in\mathbb{R}^{I_{n}\times\prod_{m\neq n}I_{m}}$, for all $n\in[N]$. Note that ${\bf X|_{k}}_{(n)}={\bf R}_{k}^{(n)}{\bf X|_{k}}_{(n)}$ always holds for these matricizations since ${\bf R}_{k}^{(n)}$ is a projection matrix.222Here we are implicitly using (1). Thus, the top $r^{(n)}_{k}$ left singular vectors of ${\bf X|_{k}}_{(n)}$ will only have nonzero entries in locations indexed by $J^{n}_{1,k}$. Let $\hat{\bf U}^{(n)}_{k}\in\mathbb{R}^{I_{n}\times r^{(n)}_{k}}$ be the matrix whose columns are these top singular vectors. As a result of the preceding discussion we can see that $\hat{\bf U}^{(n)}_{k}={\bf R}_{k}^{(n)}\hat{\bf U}^{(n)}_{k}$ will hold for all $n\in[N]$ and $k\in[K]$. Our low-rank projection matrices ${\bf Q}^{(n)}_{k}\in\mathbb{R}^{I_{n}\times I_{n}}$ used to produce low-rank approximations of each subtensor $\mathcal{X}|_{k}$ can now be defined as

[TABLE]

As a consequence of $\hat{\bf U}^{(n)}_{k}={\bf R}_{k}^{(n)}\hat{\bf U}^{(n)}_{k}$ holding, combined with the fact that $\left({\bf R}_{k}^{(n)}\right)^{\top}={\bf R}_{k}^{(n)}$ since each ${\bf R}_{k}^{(n)}$ matrix is diagonal, we have that

[TABLE]

holds for all $n\in[N]$ and $k\in[K]$. Using (24) combined with the fact that ${\bf R}_{k}^{(n)}$ is a projection matrix, we can further see that

[TABLE]

also holds for all $n\in[N]$ and $k\in[K]$.

1-scale Analysis of MS-HoSVD: Using this linear algebraic formulation we are now able to re-express the the $1^{\rm st}$ scale approximation of $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times...I_{N}}$, $\hat{\mathcal{X}}_{1}\in\mathbb{R}^{I_{1}\times I_{2}\times...I_{N}}$, as well as the $1^{\rm st}$ scale residual tensor tensor, ${\mathcal{W}}_{1}\in\mathbb{R}^{I_{1}\times I_{2}\times...I_{N}}$, as follows (see (15) – (17)). We have that

[TABLE]

holds. Thus, we see that the residual error $\mathcal{W}_{1}$ from (17) satisfies

[TABLE]

Having derived (27), it behooves us to consider when using such a first-scale approximation of $\mathcal{X}$ is actually better than, e.g., just using a standard HoSVD-based $0^{\rm th}$-scale approximation of $\mathcal{X}$ along the lines of (12). As one might expect, this depends entirely on $(i)$ how well the $1^{\rm st}$-scale partitions (i.e., the restriction matrices utilized in (20)) are chosen, as well as on $(ii)$ how well restriction matrices of the type used in (20) interact with the projection matrices used to create the standard HoSVD-based approximation in question. Toward understanding these two conditions better, recall that $\hat{\mathcal{X}}_{0}\in\mathbb{R}^{I_{1}\times I_{2}\times...I_{N}}$ in (27) is defined as

[TABLE]

where the orthogonal projection matrices ${\bf P}^{(n)}\in\mathbb{R}^{I_{n}\times I_{n}}$ are given by ${\bf P}^{(n)}=\hat{\bf U}^{(n)}\left(\hat{\bf U}^{(n)}\right)^{\top}$ for the matrices $\hat{\bf U}^{(n)}\in\mathbb{R}^{I_{n}\times r_{n}}$ used in (10). For simplicity let the ranks of the ${\bf P}^{(n)}$ projection matrices momentarily satisfy $r_{1}=r_{2}=\dots=r_{N}=:r_{0}$ (i.e., let them all be rank $r_{0}<\max_{n}\{{\rm rank}({\bf X}_{(n)})\}$). Similarly, let all the ranks, $r^{(n)}_{k}$, of the $1^{\rm st}$ scale projection matrices ${\bf Q}^{(n)}_{k}$ in (23) be $r_{1}$ for the time being.

Motivated by, e.g., the memory cost analysis of Section III-A above, one can now ask when the multiscale approximation error, $\|\mathcal{W}_{1}\|$, resulting from (27) will be less than a standard HoSVD-based approximation error, $\|\mathcal{X}-\bar{\mathcal{X}}_{0}\|$, where

[TABLE]

and each orthogonal projection matrix $\bar{{\bf P}}^{(n)}$ is of rank $\bar{r}_{n}=r_{H}\geq 2r_{0}\geq r_{0}+c^{N-1}r_{1}$ (i.e., where each $\bar{{\bf P}}^{(n)}$ projects onto the top $r_{H}$ left singular vectors of ${\bf X}_{(n)}$). In this situation having both $\|\mathcal{W}_{1}\|<\|\mathcal{X}-\bar{\mathcal{X}}_{0}\|$ and $r_{H}\geq 2r_{0}\geq r_{0}+c^{N-1}r_{1}$ hold at the same time would imply that one could achieve smaller approximation error using MS-HoSVD than using HoSVD while simultaneously achieving better compression (recall Section III-A). In order to help facilitate such analysis we prove error bounds in Appendix -A that are implied by the choice of a good partitioning scheme for the residual tensor $\mathcal{W}_{0}$ in (20) – (22).

In particular, with respect to the question concerning how well the ${1}^{\rm st}$-scale approximation error, $\|\mathcal{W}_{1}\|$, from (27) might compare to the HoSVD-based approximation error $\|\mathcal{X}-\bar{\mathcal{X}}_{0}\|$, we can use the following notion of an effective partition of $\mathcal{W}_{0}$. The partition of $\mathcal{W}_{0}$ formed by the restriction matrices ${\bf R}_{k}^{(n)}$ in (20) – (22) will be called effective if there exists another pessimistic partitioning of $\mathcal{W}_{0}$ via (potentially different) restriction matrices $\left\{\tilde{\bf R}_{k}^{(n)}\right\}^{K}_{k=1}$ together with a bijection $f:[K]\rightarrow[K]$ such that

[TABLE]

holds for each $k\in[K]$. In (30) the $\left\{\tilde{\bf P}^{(n)}\right\}$ are the orthogonal projection matrices obtained from the HoSVD of $\mathcal{W}_{0}$ with ranks $\tilde{r}_{n}=r_{H}$ (i.e., where each $\tilde{{\bf P}}^{(n)}$ projects onto the top $\tilde{r}_{n}=r_{H}$ left singular vectors of the matricization ${\bf W}_{0,(n)}$). In Appendix -A, we show that (30) holding for $\mathcal{W}_{0}$ implies that the error $\|\mathcal{W}_{1}\|$ resulting from our $1^{\rm st}$-scale approximation in (27) is less than an upper bound of the type often used for HoSVD-based approximation errors of the form $\|\mathcal{X}-\bar{\mathcal{X}}_{0}\|$ (see, e.g., [29]). In particular, we prove the following result.

Theorem 1**.**

Suppose that (30) holds. Then, the first scale approximation error given by MS-HoSVD in (27) is bounded by

[TABLE]

where $\left\{\bar{\bf P}^{(n)}\right\}$ are low-rank projection matrices of rank $\bar{r}_{n}=\tilde{r}_{n}=r_{H}$ obtained from the truncated HoSVD of $\mathcal{X}$ as per (29).

Proof.

See Appendix -A. ∎

Theorem 1 implies that $\left\|\mathcal{W}_{1}\right\|$ may be less than $\|\mathcal{X}-\bar{\mathcal{X}}_{0}\|$ when (30) holds. It does not, however, actually prove that $\left\|\mathcal{W}_{1}\right\|\leq\|\mathcal{X}-\bar{\mathcal{X}}_{0}\|$ holds whenever (30) does. In fact, directly proving that $\left\|\mathcal{W}_{1}\right\|\leq\|\mathcal{X}-\bar{\mathcal{X}}_{0}\|$ whenever (30) holds does not appear to be easy. It also does not appear to be easy to prove the error bound in theorem 1 without an assumption along the lines of (30) which simultaneously controls both $(i)$ how well the restriction matrices utilized to partition $\mathcal{W}_{0}$ in (21) are chosen, as well as $(ii)$ how poorly (worst case) restriction matrices interact with the projection matrices used to create standard HoSVD-based approximations of $\mathcal{W}_{0}$ and/or $\mathcal{X}$. The development of simpler and/or weaker conditions than (30) which still yield meaningful error guarantees along the lines of theorem 1 is left as future work. See Appendix -A for additional details and comments, and Appendix -B below for an example illustrating Theorem 1 on an idealized tensor..

III-D Adaptive Pruning in Multiscale HoSVD for Improved Performance

In order to better capture the local structure of the tensor, it is important to look at higher scale decompositions. However, as the scale increases, the storage cost and computational complexity will increase making any gain in reconstruction error potentially not worth the additional memory cost. For this reason, it is important to carefully select the subtensors adaptively at higher scales. To help avoid the redundancy in decomposition structure we propose an adaptive pruning method across scales.

In adaptive pruning, the tree is pruned by minimizing the following cost function $\mathbb{H}=Error+\lambda\cdot Compression$ similar to the rate-distortion criterion commonly used by compression algorithms where $\lambda$ is the trade-off parameter [34]. To minimize this function we employ a greedy procedure similar to sequential forward selection [35]. First, the root node which stores $\hat{\mathcal{X}}_{0}$ is created and scale-1 subtensors $\hat{\mathcal{X}}_{1,k}$ are obtained from the 0th order residual tensor $\hat{\mathcal{W}}_{0}$ as discussed in Section III. These subtensors are stored in a list and the subtensor which decreases the cost function the most is then added to the tree structure under its parent node. Next, scale-2 subtensors belonging to the added node are created and added to the list. All of the scale-1 and scale-2 subtensors in the list are again evaluated to find the subtensor that minimizes the cost function. This procedure is repeated until the cost function $\mathbb{H}$ converges or the decrease is minimal. A pseudocode of the algorithm is given in Algorithm 2. It is important to note that this algorithm is suboptimal similar to other greedy search methods.

IV Data Reduction

In this section we demonstrate the performance of MS-HoSVD for tensor type data reduction on several real 3-mode and 4-mode datasets as compared with three other tensor decompositions: HoSVD, H-Tucker, and T-Train. The performance of tensor decomposition methods are evaluated in terms of reconstruction error and compression rate. In the tables and figures below the error rate refers to the normalized tensor approximation error $\frac{\|\mathcal{X}-\hat{\mathcal{X}}\|_{F}}{\|\mathcal{X}\|_{F}}$ and the compression rate is computed as $\frac{\#{\rm~{}total~{}bits~{}to~{}store~{}}\hat{\mathcal{X}}}{\#{\rm~{}total~{}bits~{}to~{}store~{}}\mathcal{X}}$. Moreover, we show the performance of the proposed adaptive tree pruning strategy for data reduction.

IV-A Datasets

IV-A1 PIE dataset

A 3-mode tensor $\mathcal{X}\in\mathbb{R}^{244\times 320\times 138}$ is created from PIE dataset [36]. The tensor contains 138 images from 6 different yaw angles and varying illumination conditions collected from a subject where each image is converted to gray scale. Fig. 2 illustrates the images from different frames of the PIE dataset.

IV-A2 COIL-100 dataset

The COIL-100 database contains 7200 images collected from 100 objects where the images of each object were taken at pose intervals of $5^{\circ}$. A 4-mode tensor $\mathcal{X}\in\mathbb{R}^{128\times 128\times 72\times 100}$ is created from COIL-100 dataset [37]. The constructed 4-mode tensor contains 72 images of size $128\times 128$ from 100 objects where each image is converted to gray scale. In Fig. 3, sample images of four objects taken from different angles can be seen.

IV-A3 The Cambridge Hand Gesture Dataset

The Cambridge hand gesture database consists of 900 image sequences of nine gesture classes of three primitive hand shapes and three primitive motions where each class contains 100 image sequences (5 different illuminations $\times$ 10 arbitrary motions $\times$ 2 subjects). In Fig. 4, sample image sequences collected for nine hand gestures can be seen. The created 4-mode tensor $\mathcal{X}\in\mathbb{R}^{60\times 80\times 30\times 900}$ contains 900 image sequences of size $60\times 80\times 30$ where each image is converted to gray scale.

IV-B Data Reduction Experiments

In this section, we evaluate the performance of MS-HoSVD for 1 and 2-scale decompositions compared to HoSVD, H-Tucker and T-Train decompositions. In the following experiments, tensor partitioning is performed by LSA and the cluster number along each mode is chosen as $c_{i}=2$. The rank used in HoSVD is selected adaptively using the energy criterion as per Section III’s (11). In our experiments, we performed MS-HoSVD with $\tau=0.7$ and $\tau=0.75$ and we kept $\tau$ the same for each scale. For the same compression rates as the MS-HoSVD, the reconstruction error of HoSVD, H-Tucker and T-Train models are computed.

Fig. 5 explores the interplay between compression rate and approximation error for MS-HoSVD in comparison to HoSVD, H-Tucker and T-Train for PIE, COIL-100 and hand gesture datasets. Starting from the left in Figs. 5(a), 5(b) and 5(c), the first two compression rates correspond to 1-scale MS-HoSVD with $\tau=0.7$ and $\tau=0.75$, respectively while the last two are obtained from 2-scale approximation with $\tau=0.7$ and $\tau=0.75$, respectively. As seen in Fig. 5, MS-HoSVD outperforms other approaches with respect to reducing PIE, COIL-100 and hand gesture tensors at varying compression rates. Moreover, adding the $2^{\rm nd}$ scale increases the storage requirements while decreasing the error of MS-HoSVD. Fig. 6 illustrates the influence of scale on the visual quality of the reconstructed images. As expected, introducing additional finer scales into a multiscale approximation of video data improves image detail in each frame. Moreover, the data reduction performance of T-Train is seen to be slightly better than H-Tucker in most of the experiments.

IV-C Data Reduction with Adaptive Tree Pruning

In this section, we evaluate the performance of adaptive tree pruning multiscale decompositions. In the pruning experiments, clustering is performed by LSA and the cluster number along each mode is chosen as $c_{i}=2$. The rank used in HoSVD is selected adaptively based on the energy threshold $\tau=0.7$. A pruned version of 2-scale MS-HoSVD that greedily minimizes the cost function $\mathbb{H}=Error+\lambda\cdot Compression$ for is implemented for PIE, COIL-100 and Hand Gesture datasets with varying $\lambda$ values as reported in Tables I, II and III. As $\lambda$ increases, reducing the compression rate becomes more important and the algorithm prunes the leaf nodes more. For example, a choice of $\lambda=0.75$ prunes all of the nodes corresponding to the second scale subtensors for PIE data (see Table I).

As can be seen from Tables I, II, and III, the best tradeoffs achieved between reconstruction error and compression rate occur at different $\lambda$ values for different datasets. For example, for PIE data, increasing $\lambda$ value does not provide much change in reconstruction error while increasing the compression. On the other hand, for COIL-100, $\lambda=0.75$ provides a good tradeoff between reconstruction error and compression rate. Small changes in $\lambda$ yield significant effects on pruning the subtensors of 2-scale decomposition of hand gesture data. Fig. 7 illustrates the performance of the pruning algorithm on the PIE dataset. Applying pruning with $\lambda=0.25$ increases the reconstruction error from 0.0276 to 0.0506 while reducing the compression rate by a factor of 4 (Table I). As seen in Fig. 7, the $2^{\rm nd}$-scale approximation obtained by the adaptive pruning algorithm preserves most of the facial details in the image.

Performance of the pruning algorithm reported in Tables I, II and III is also compared with HoSVD, H-Tucker and T-Train decompositions in Fig. 8. As seen in Fig. 8 (b) and (c), MS-HoSVD outperforms other approaches for compressing COIL-100 and Hand Gesture datasets at varying compression rates. However, for PIE data, the performance of MS-HoSVD and HoSVD are very close to each other while both approaches outperform H-Tucker and T-Train, as can be seen in Fig. 8 (a). In Fig. 9, sample frames of PIE data reconstructed by T-Train (top-left), H-Tucker (top-right), HoSVD (bottom-left) and pruned MS-HoSVD with 2-scales (bottom-right) are shown. It can be easily seen that the reconstructed images by H-Tucker and T-Train are more blurred than the ones obtained by HoSVD and MS-HoSVD. One can also see the facial details captured by MS-HoSVD are clearer than HoSVD although the performances of both algorithms are very similar to each other. The reason for capturing facial details better by MS-HoSVD is that the higher-scale subtensors encode facial details.

V Feature Extraction and Classification

In this section, we evaluate the features extracted from MS-HoSVD for classification of 2-mode and 3-mode tensors containing object images and hand gesture videos.

The classification accuracy of MS-HoSVD features are compared to the features extracted by HoSVD and T-Train using three different classifiers: 1-NN, Adaboost and Naive Bayes.

V-A COIL-100 Image Dataset

For computational efficiency, each image was downsampled to a gray-scale image of $32\times 32$ pixels. Number of images per object used for training data was gradually increased from 18 to 54 and selected randomly. A 3-mode tensor $\mathcal{X}^{tr}\in\mathbb{R}^{32\times 32\times I_{tr}}$ is constructed from training images where $I_{tr}\in 100\times\left\{18,\;36,\;54\right\}$ and the rest of the images are used to create the testing tensor $\mathcal{X}^{te}\in\mathbb{R}^{32\times 32\times I_{te}}$ where $I_{te}=7200-I_{tr}$.

V-B The Cambridge Hand Gesture Dataset

For computational efficiency, each image was downsampled to a gray-scale image of $30\times 40$ pixels. Number of image sequences used for training data gradually increased from 25 to 75 per gesture and selected randomly. A 4-mode tensor $\mathcal{X}^{tr}\in\mathbb{R}^{30\times 40\times 30\times I_{tr}}$ is constructed from training image sequences where $I_{tr}\in 9\times\left\{25,\;50,\;75\right\}$ and the rest of the image sequences are used to create the testing tensor $\mathcal{X}^{te}\in\mathbb{R}^{30\times 40\times 15\times I_{te}}$ where $I_{te}=900-I_{tr}$.

V-C Classification Experiments

V-C1 Training

For MS-HoSVD, the training tensor $\mathcal{X}^{tr}$ is decomposed using 1-scale MS-HoSVD as follows. Tensor partitioning is performed by LSA and the cluster number along each mode is chosen as $c={\left\{2,\;3,\;1\right\}}$ yielding 6 subtensors for COIL-100 dataset and $c={\left\{2,\;2,\;3,\;1\right\}}$ yielding 12 subtensors for hand gesture dataset. We did not partition the tensor along the last mode that corresponds to the classes to make the comparison with other methods fair. The rank used in 0th scale is selected based on the energy criterion with $\tau=0.7$, while the full rank decomposition is used for the 1st scale. The 0th scale approximation

[TABLE]

provides the $0^{\rm th}$-scale core tensor $\mathcal{C}_{0}^{tr}$, factor matrices $\hat{\bf U}^{tr,(i)}$ and residual tensor $\mathcal{W}_{0}^{tr}=\mathcal{X}^{tr}-\hat{\mathcal{X}}_{0}^{tr}$. Next, the $0^{\rm th}$-scale feature tensor $\mathcal{S}^{tr}_{0}$ for the training data is created by projecting ${\mathcal{X}}^{tr}$s onto the first $N-1$ factor matrices ${\bf U}^{tr,(i)}$ as:

[TABLE]

Subtensors of $\mathcal{W}_{0}^{tr}$ obtained by $\mathcal{X}^{tr}_{1,k}=\mathcal{W}_{0}^{tr}(J^{tr,1}_{1,k}\times J^{tr,2}_{1,k}\times...\times J^{tr,N}_{1,k})$ are used to extract $1^{st}$-order core tensors $\mathcal{C}_{1,k}$ and factor matrices ${\bf U}^{{tr},(i)}_{1,k}$ as:

[TABLE]

$1st$-order feature tensors are then created by projecting ${\mathcal{X}}^{tr}_{1,k}$s onto the first $N-1$ factor matrices ${\bf U}^{tr,(i)}_{1,k}$ as:

[TABLE]

Unfolding the feature tensors $\mathcal{S}^{tr}_{0}$ and $\mathcal{S}^{tr}_{1,k}$ along the sample mode $N$ and concatenating them to each other yields a high dimensional feature vector for each of the training samples. From these vectors, $N_{f}$ features with the highest Fisher Score [38] are selected to form the lower-dimensional feature vectors ${\bf x}^{tr}\in\mathbb{R}^{N_{f}\times 1}$ for each training sample where the number of features ($N_{f}$) is determined 100 for COIL-100 and 200 for hand gesture dataset emprically. For HoSVD and T-Train, full rank decompositions are computed and feature vectors are created by selecting $N_{f}$ features with the highest Fisher Score from the core tensors as described above. For T-Train, the procedure described in [39] is used without reducing the dimensionality.

V-C2 Testing

To create the $0^{\rm th}$-order feature tensor $\mathcal{S}^{te}_{0}$ for testing samples, first, the testing tensor $\mathcal{X}_{te}$ is projected onto $\hat{\bf U}^{tr,(i)}$ where $i\in\left[N-1\right]$ as:

[TABLE]

The $0^{\rm th}$-order residual tensor $\mathcal{W}_{0}^{te}$ of testing data is computed as $\mathcal{W}_{0}^{te}={\mathcal{X}}^{te}-{\mathcal{X}}^{te}\bigtimes_{n=1}^{n=N}\left(\hat{\bf U}^{tr,(n)}\right)^{\top}$. Then $1^{\rm st}$-order subtensors are created from $\mathcal{W}_{0}^{te}$ using the same partitioning as the 0th order training residual tensor $\mathcal{W}_{0}^{tr}$ as $\mathcal{X}^{te}_{1,k}=\mathcal{W}_{0}^{te}(J^{tr,1}_{1,k}\times J^{tr,2}_{1,k}\times...\times J^{tr,N}_{1,k})$. The $1^{\rm st}$-order feature tensors $\mathcal{S}^{te}_{1,k}$ for the testing samples are then obtained by

[TABLE]

Similar to the training step, unfolding the feature tensors $\mathcal{S}^{te}_{0}$ and $\mathcal{S}^{te}_{1,k}$ along the sample mode $N$ and concatenating them with each other yields high dimensional feature vectors for the testing samples. The features corresponding to the features selected from the training step are used to form the feature vectors for testing samples ${\bf x}^{te}\in\mathbb{R}^{N_{f}\times 1}$. A similar two-step procedure, i.e projecting the testing tensor onto training factor matrices followed by selecting $N_{f}$ features, is used to create testing feature vectors for HoSVD and T-Train. Discrimination performance of the feature vectors are evaluated using different classifiers including 1-NN, Adaboost and Naive Bayes.

Tables IV and V summarize the classification accuracy for the three methods using three different classifiers for COIL-100 and Hand gesture data sets, respectively. As it can be seen from these Tables, for both data sets and all classifiers MS-HoSVD performs the best except for a Naive Bayes Classifier trained by $25\%$ of the data to classify hand gesture dataset. As seen in Tables IV and V, the performance of HoSVD, T-Train and MS-HoSVD become close to each other as the size of the training dataset increases, as expected. The reason for the superior performance of MS-HoSVD is that MS-HoSVD captures the variations and nonlinearities across the modes such as rotation or translation better than the other methods. In both of the datasets used in this section, the images are rotated across different frames. Since these nonlinearities are encoded in the higher-scale ($1^{\rm st}$-scale) features while the average characteristics, which are the same as HoSVD, are captured by the lower scale ($0^{\rm th}$-scale) MS-HoSVD features, the classification performance of the MS-HoSVD is slightly better than HoSVD. It is also seen that T-Train features are not as good as MS-HoSVD and HoSVD features for capturing rotations and translations in the data and requires larger training set to reach the performance of MS-HoSVD and HoSVD.

VI Conclusions

In this paper, we proposed a new multi-scale tensor decomposition technique for better approximating the local nonlinearities in generic tensor data. The proposed approach constructs a tree structure by considering similarities along different fibers of the tensor and decomposes the tensor into lower dimensional subtensors hierarchically. A low-rank approximation of each subtensor is then obtained by HoSVD. We also introduced a pruning strategy to find the optimum tree structure by keeping the important nodes and eliminating redundancy in the data. The proposed approach is applied to a set of $3$-way and $4$-way tensors to evaluate its performance on both data reduction and classification applications. As it is illustrated in sections IV and V, any application involving tensor data reduction and classification would benefit from the proposed method. Some examples include hyper-spectral image compression, high-dimensional video clustering and functional connectivity network analysis in neuroscience.

Although this paper focused on the integration of a single existing tensor factorization technique (i.e., the HoSVD) into a clustering-enhanced multiscale approximation framework, we would like to emphasize that the ideas presented herein are significantly more general. In principal, for example, there is nothing impeding the development of multiscale variants of other tensor factorization approaches (e.g., PARAFAC, T-Train, H-Tucker, etc.) in essentially the same way. In this paper it is demonstrated that the use of the HoSVD as part of a multiscale approximation approach leads to improved compression and classification performance over standard HoSVD approaches. However, this paper should additionally be considered as evidence that similar improvements are also likely possible for other tensor factorization-based compression and classification schemes, as well as for other related applications.

Future work will consider automatic selection of parameters such as the number of clusters and the appropriate rank along each mode. The computational efficiency of the proposed method can also be improved through parallelization of the algorithm by, e.g., constructing the disjoint subtensors at each scale in parallel, as well as by utilizing distributed and parallel SVD algorithms such as [33] when computing their required HoSVD decompositions (see also, e.g., [40] for other related parallel implementations). Such efficient implementations will enable the computation of finer-scale decompositions for higher-order and higher-dimensional tensors.

-A Effective Partitioning, and Error Analysis

In order to facilitate error analysis for the 1-scale MS-HoSVD that is similar to the types of error analysis available for various HoSVD-based low-rank approximation strategies (see, e.g., [29]), we will engage in a more in depth discussion of condition (30) herein. Recall that the partition of $\mathcal{W}_{0}$ formed by the restriction matrices ${\bf R}_{k}^{(n)}$ in (20) – (22) is called effective if there exists another pessimistic partitioning of $\mathcal{W}_{0}$ via restriction matrices $\left\{\tilde{\bf R}_{k}^{(n)}\right\}^{K}_{k=1}$ together with a bijection $f:[K]\rightarrow[K]$ such that

[TABLE]

holds for each $k\in[K]$. In (37) the $\left\{\tilde{\bf P}^{(n)}\right\}$ are the orthogonal projection matrices obtained from the HoSVD of $\mathcal{W}_{0}$ with ranks $\tilde{r}_{n}\geq\bar{r}_{n}\geq r_{n}$ (i.e., where each $\tilde{{\bf P}}^{(n)}$ projects onto the top $\tilde{r}_{n}$ left singular vectors of the matricization ${\bf W}_{0,(n)}$). Below we will show that (37) holding for $\mathcal{W}_{0}$ implies that the error $\|\mathcal{W}_{1}\|$ resulting from our $1^{\rm st}$-scale approximation in (27) is less than an upper bound of the type given for a high-rank standard HoSVD-based approximation (29) in [29].

Considering condition (37) above, we note that experiments show that it is regularly satisfied on real datasets when $(i)$ the effective restriction matrices $\left\{{\bf R}_{k}^{(n)}\right\}^{K}_{k=1}$ in (20) – (22) are first formed by clustering the rows of each unfolding of $\mathcal{W}_{0}$ using, e.g., local subspace analysis (LSA), after which $(ii)$ pessimistic restriction matrices $\left\{\tilde{\bf R}_{k}^{(n)}\right\}^{K}_{k=1}$ are randomly generated in order to create another (random) partition of $\mathcal{W}_{0}$ into $K$ different disjoint subtensors for comparison. The bijection $f$ can then be created by, e.g., $(i)$ sorting the left-hand side errors in (37) for each $k\in[k]$, $(ii)$ sorting the right-hand side errors in (37) for each $k\in[K]$, and then $(iii)$ matching the largest left-hand and right-hand errors for comparision, the second largest left-hand and right-hand errors for comparision, etc.. When checked in this way the sorted right-hand side errors often dominate (entrywise) the sorted left-hand side errors for various reasonable ranks $\bar{r}_{n}=\tilde{r}_{n}=r_{H}=r_{0}+c^{N-1}r_{1}$ (as a function of $r_{0}$ and $r_{1}$ with, e.g., $c=2$) on every dataset considered in Section IV above, thereby verifying that (37) does indeed regularly hold.

We will now begin to prove Theorem 1 with a lemma that shows our subtensor-based approximation of $\mathcal{W}_{0}$ is accurate whenever (37) is satisfied.

Lemma 3**.**

Let $\mathcal{W}_{0}=\mathcal{X}-\hat{\mathcal{X}}_{0}\in\mathbb{R}^{I_{1}\times I_{2}\times...I_{N}}$. Suppose that $\left\{{\bf R}_{k}^{(n)}\right\}$ is a collection of effective restriction matrices that form an effective partition of $\mathcal{W}_{0}$ with respect to a pessimistic partition formed via pessimistic restriction matrices $\left\{\tilde{\bf R}_{k}^{(n)}\right\}$ as per (37) above. Similarly, let $\tilde{\bf P}^{(n)}$ be the rank $\tilde{r}_{n}\geq\bar{r}_{n}~{}\forall n$ orthogonal projection matrices from (37) obtained via the truncated HoSVD of $\mathcal{W}_{0}$ as above. Then,

[TABLE]

Proof.

We have that

[TABLE]

Applying lemmas 1 and 2 to (38) we can now see that

[TABLE]

since the ${\bf Q}^{(n)}_{k}$ matrices are orthogonal projections. Using assumption (37) we now get that

[TABLE]

where we have used the fact that the pessimistic restriction matrices $\tilde{\bf R}_{k}^{(n)}$ partition $\mathcal{W}_{0}$ in the last line. ∎

Lemma 3 indicates that the error in approximating $\mathcal{W}_{0}$ via low-rank approximations of its effective subtensors is potentially smaller than the error obtained by approximating $\mathcal{W}_{0}$ via (higher-rank) truncated HoSVDs whenever (37) holds.333That is, the upper bound on the error provided by Lemma 3 is less than or equal to the upper bound on the error for truncated HoSVDs provided by, e.g., [29] when/if (37) holds. The following theorem shows that this good error behavior extends to the entire $1^{\rm st}$-scale approximation provided by (27) whenever (37) holds.

Theorem 2** (Restatement of Theorem 1).**

Let $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}...\times I_{N}}$. Suppose that (37) holds. Then, the first-scale approximation error given by MS-HoSVD (27) is bounded by

[TABLE]

where $\left\{\bar{\bf P}^{(n)}\right\}$ are low-rank projection matrices of rank $\bar{r}_{n}\geq r_{n}$ obtained from the truncated HoSVD of $\mathcal{X}$ as per (29).

Proof.

Using (12) and (17) together with lemma 3 we can see that

[TABLE]

where $\bar{\bf Q}^{(n)}\in\mathbb{R}^{I_{n}\times I_{n}}$ is the orthogonal projection matrix of rank $\tilde{r}_{n}$ which projects onto the subspace spanned by the top $\tilde{r}_{n}$ left singular vectors of ${\bf X}_{(n)}$. Here (39) holds because the orthogonal projection matrices $\tilde{\bf P}^{(n)}$ are chosen in (37) so that $\tilde{\bf P}^{(n)}{\bf W}_{0,(n)}$ is a best possible rank $\tilde{r}_{n}$ approximation to ${\bf W}_{0,(n)}$. As a result, we have that

[TABLE]

must hold for each $n\in[N]$.

Continuing from (39) we can use the definition of $\hat{\mathcal{X}}_{0}$ in (28) to see that

[TABLE]

by lemma 1. Due to the definition of $\bar{\bf Q}^{(n)}$ together with the fact that its rank is $\tilde{r}_{n}\geq r_{n}$ we can see that $\left({\bf I}-\bar{\bf Q}^{(n)}\right){\bf P}^{(n)}={\bf 0}$. As a consequence, lemma 1 implies that $\mathcal{X}\bigtimes_{h=1}^{N}{\bf P}^{(h)}\times_{n}\left({\bf I}-\bar{\bf Q}^{(n)}\right)={\bf 0}$ for all $n\in[N]$. Continuing from (40) we now have that

[TABLE]

Again, appealing to the definition of both $\bar{\bf Q}^{(n)}$ and $\bar{\bf P}^{(n)}$ in (29), combined with the fact that $\tilde{r}_{n}\geq\bar{r}_{n}$, finally yields the desired result. ∎

We refer the reader to the strong empirical performance of MS-HoSVD in Section IV for additional evidence supporting the utility of (27) as a means of improving the compression performance of standard HoSVD-based compression techniques. In addition, we further refer the reader to Section V where it is empirically demonstrated that MS-HoSVD is also capable of selecting more informative features than HoSVD-based methods for the purposes of classification. These two facts together provide strong evidence that combining the use of clustering-enhanced multiscale approximation with existing tensor factorization techniques can lead to improved performance in multiple application domains.

-B Experiment for Error Analysis

In this experiment we evaluate the error obtained by the $1^{\rm st}$-scale MS-HoSVD analysis of a tensor along the lines of the model described in Section III. Herein we consider a three-way tensor $\mathcal{X}\in\mathbb{R}^{20\times 20\times 20}$ that is the sum of two tensors as $\mathcal{X}=\mathcal{X}_{0}+\mathcal{X}_{1}$ where $\mathcal{X}_{0}\in\mathbb{R}^{20\times 20\times 20}$ has $n$-rank $\left(2,\;2,\;2\right)$, and $\mathcal{X}_{1}\in\mathbb{R}^{20\times 20\times 20}$ is formed by concatenating $8$ subtensors $\mathcal{X}_{k}\in\mathbb{R}^{10\times 10\times 10}$ each also with $n$-rank $\left(2,\;2,\;2\right)$. Low-rank approximations for $\mathcal{X}$ and its subtensors are always obtained via the truncated HoSVD. The 1-scale MS-HoSVD is applied with the ground truth partitions ${\bf R}_{k}^{(n)}$, partitions provided by Local Subspace Analysis (LSA) clustering $\hat{\bf R}_{k}^{(n)}$, and also with randomly chosen partitions $\tilde{\bf R}_{k}^{(n)}$ of the $0^{\rm th}$-scale residual error $\mathcal{W}_{0}$ into $8$ different $10\times 10\times 10$ subtensors. For the LSA clustering, the cluster numbers are selected as $2$ along each mode also yielding $8$ subtensors.

The $0^{\rm th}$-scale $n$-rank for MS-HoSVD is selected as $\left(2,\;2,\;2\right)$, and the $1^{\rm st}$-scale ranks are varied in the experiments as shown in Table VI. The normalized reconstruction error computed for these varying $1^{\rm st}$-scale $n$-ranks can also be seen in Table VI. As seen there, using ground truth partition provides lower-rank subtensors, and using clustering as part of the $1$-scale MS-HoSVD leads to much better approximations than HoSVD does in general.

In addition, the left (LHS) and right (RHS) sides in Theorem 1 are computed where the first-scale projections ${\bf Q}^{(n)}_{k}$ each have rank $r_{1}=2$ and are obtained via both ground truth partitioning and clustering after the $0^{\rm th}$-scale ${\bf P}^{(n)}$ are obtained from the truncated HoSVD of $\mathcal{X}$ with $r_{0}=2$. For comparison the $\bar{\bf P}^{(n)}$ are also computed from the truncated HoSVD of $\mathcal{X}$ with varying ranks $\bar{r}_{n}=\kappa r_{0}$ for $\kappa\in\left\{2\;,3\;,4\right\}$. In Table VII, we report the mean value and standard deviation of both the right-hand side (RHS) and left-hand side (LHS) of the error bound in Theorem 1. As seen in Table VII, Theorem 1 holds for the $1^{\rm st}$-scale approximation of $\mathcal{X}$ via MS-HoSVD since the RHS errors based on the $\bar{\bf P}^{(n)}$ projections are larger than the LHS errors no matter whether the ${\bf Q}^{(n)}_{k}$ are obtained via ground truth partitioning or LSA clustering.