Tensor Completion from Regular Sub-Nyquist Samples

TL;DR

This paper demonstrates that regular sampling can effectively reconstruct high-dimensional tensor signals, such as images and videos, with sample complexity depending on tensor rank, and applies this to accelerate fMRI imaging.

Contribution

It introduces a novel framework for tensor signal reconstruction from regular samples, emphasizing tensor rank over bandwidth, and applies it to practical fMRI acceleration.

Findings

Reconstruction from regular samples is feasible for tensors.

Sample complexity depends on tensor rank, not bandwidth.

The method significantly accelerates fMRI sampling without loss of accuracy.

Abstract

Signal sampling and reconstruction is a fundamental engineering task at the heart of signal processing. The celebrated Shannon-Nyquist theorem guarantees perfect signal reconstruction from uniform samples, obtained at a rate twice the maximum frequency present in the signal. Unfortunately a large number of signals of interest are far from being band-limited. This motivated research on reconstruction from sub-Nyquist samples, which mainly hinges on the use of random / incoherent sampling procedures. However, uniform or regular sampling is more appealing in practice and from the system design point of view, as it is far simpler to implement, and often necessary due to system constraints. In this work, we study regular sampling and reconstruction of three- or higher-dimensional signals (tensors). We show that reconstructing a tensor signal from regular samples is feasible. Under the proposed framework, the sample complexity is determined by the tensor rank---rather than the signal bandwidth. This result offers new perspectives for designing practical regular sampling patterns and systems for signals that are naturally tensors, e.g., images and video. For a concrete application, we show that functional magnetic resonance imaging (fMRI) acceleration is a tensor sampling problem, and design practical sampling schemes and an algorithmic framework to handle it. Numerical results show that our tensor sampling strategy accelerates the fMRI sampling process significantly without sacrificing reconstruction accuracy.

Tables

Table 1. TABLE I : Reconstruction performance of the competing algorithms.

Algorithm	RETSINA	k-t Focuss	k-t SLR	IDFT
NRE	0.124	0.339	1.41	0.8156
${\text{NRE}}_2$	0.081	0.286	0.073	0.7376
runtime	12min	25.6min (48sec/coil)	480min (15min/coil)	14sec

Table 2. TABLE II : NMSE performance of MS-RETSINA .

$s\backslash r$	2	3	4	5	6	7	8	9	10
$2$	0.16	0.17	0.18	0.19	0.20	0.19	0.19	0.19	0.21
$3$	0.18	0.19	0.20	0.19	0.20	0.22	0.21	0.20	0.22
$4$	0.20	0.20	0.20	0.22	0.21	0.21	0.23	0.22	0.23

Equations

$$\underline{\mathbf{X}}=\sum _{f=1}^F{\mathbf{a}}_f\circ {\mathbf{b}}_f\circ {\mathbf{c}}_f.$$

$${\mathbf{X}}^{(3)}:=[\text{vec}(\underline{\mathbf{X}}(:,:,1)),\text{vec}(\underline{\mathbf{X}}(:,:,2)),\dots ,\text{vec}(\underline{\mathbf{X}}(:,:,K))],$$

$${\mathbf{X}}^{(3)}=(\mathbf{B}\odot \mathbf{A}){\mathbf{C}}^T\in {\mathbb{C}}^{IJ\times K},$$

$$\widetilde{\underline{\mathbf{X}}}=\underline{\mathbf{X}}{\times }_1{\mathbf{P}}_1{\times }_2{\mathbf{P}}_2{\times }_3{\mathbf{P}}_3$$

$$\underline{\tilde{\mathbf{X}}}=[[{\mathbf{P}}_1\mathbf{A},{\mathbf{P}}_2\mathbf{B},{\mathbf{P}}_3\mathbf{C}]],$$

$$\mathbf{y}=\mathsf{Sample}(\underline{\mathbf{X}}),$$

$$\underline{\mathbf{X}}({\mathcal{S}}_r,{\mathcal{S}}_c,{\mathcal{S}}_f)=[[\mathbf{A}({\mathcal{S}}_r,:),\mathbf{B}({\mathcal{S}}_c,:),\mathbf{C}({\mathcal{S}}_f,:)]].$$

$${\underline{\mathbf{Y}}}_1=\underline{\mathbf{X}}\left({\mathcal{S}}_h,:,:\right)=\underline{\mathbf{X}}{\times }_1{\mathbf{P}}_1^{(1)}$$

$${\underline{\mathbf{Y}}}_2=\underline{\mathbf{X}}\left(:,:,{\mathcal{S}}_f\right)=\underline{\mathbf{X}}{\times }_3{\mathbf{P}}_3^{(2)}$$

$${\underline{\mathbf{Y}}}_1$$

$${\underline{\mathbf{Y}}}_2$$

$$|{\mathcal{S}}_r^{(d)}|,|{\mathcal{S}}_c^{(d)}|\ge 2$$

$$\bigcup _{d=1}^D{\mathcal{S}}_r^{(d)}=\{1,\dots ,I\},\bigcup _{d=1}^D{\mathcal{S}}_c^{(d)}=\{1,\dots ,J\}$$

$$\bigcap _{d=1}^D\bigcup _{d^{\prime }\ne d}{\mathcal{S}}_c^{(d)}\cap {\mathcal{S}}_c^{(d^{\prime })}\bigcup {\mathcal{S}}_r^{(d)}\cap {\mathcal{S}}_r^{(d^{\prime })}\ne \varnothing ,$$

$${\underline{\mathbf{Y}}}_d=$$

$$=$$

$$|{\mathcal{S}}_r^{(d)}|,|{\mathcal{S}}_c^{(d)}|,|{\mathcal{S}}_f^{(d)}|\ge 2$$

$$\bigcup _{d=1}^D{\mathcal{S}}_r^{(d)}=\{1,\dots ,I\},\bigcup _{d=1}^D{\mathcal{S}}_c^{(d)}=\{1,\dots ,J\},$$

$$\bigcup _{d=1}^D{\mathcal{S}}_f^{(d)}=\{1,\dots ,K\}$$

$$\bigcup _{m\ne m^{\prime }}\bigcap _{d=1}^D\bigcup _{d^{\prime }\ne d}\left\{{\mathcal{S}}_m^{(d)}\cap {\mathcal{S}}_m^{(d^{\prime })}\bigcap {\mathcal{S}}_{m^{\prime }}^{(d)}\cap {\mathcal{S}}_{m^{\prime }}^{(d^{\prime })}\right\}\ne \varnothing$$

$$\bigcup _m\left\{|{\mathcal{S}}_m^{(d)}\cap {\mathcal{S}}_m^{(d^{\prime })}|\ge 2\right\}\ne \varnothing ,$$

$${\underline{\mathbf{Y}}}_d=$$

$$=$$

$${\mathbf{Y}}_1^{(3)}=(\mathbf{B}\odot {\mathbf{P}}_1^{(1)}\mathbf{A}){\mathbf{C}}^T$$

$${\mathbf{A}}_d={\mathbf{P}}_1^{(d)}\mathbf{A}{\mathbf{\Pi }}^{(d)}{\mathbf{\Lambda }}_1^{(d)}=\mathbf{A}({\mathcal{S}}_r^{(d)},:){\mathbf{\Pi }}^{(d)}{\mathbf{\Lambda }}_1^{(d)},$$

$${\mathbf{B}}_d={\mathbf{P}}_2^{(d)}\mathbf{B}{\mathbf{\Pi }}^{(d)}{\mathbf{\Lambda }}_2^{(d)}=\mathbf{B}({\mathcal{S}}_c^{(d)},:){\mathbf{\Pi }}^{(d)}{\mathbf{\Lambda }}_2^{(d)},$$

$${\mathbf{C}}_d={\mathbf{P}}_3^{(d)}\mathbf{C}{\mathbf{\Pi }}^{(d)}{\mathbf{\Lambda }}_3^{(d)}=\mathbf{C}({\mathcal{S}}_f^{(d)},:){\mathbf{\Pi }}^{(d)}{\mathbf{\Lambda }}_3^{(d)},$$

$${\mathbf{C}}_d^{(d-d^{\prime })}=\mathbf{C}({\mathcal{S}}_f^{(d-d^{\prime })},:){\mathbf{\Pi }}^{(d)}{\mathbf{\Lambda }}_3^{(d)}$$

$${\mathbf{C}}_{d^{\prime }}^{(d-d^{\prime })}=\mathbf{C}({\mathcal{S}}_f^{(d-d^{\prime })},:){\mathbf{\Pi }}^{(d^{\prime })}{\mathbf{\Lambda }}_3^{(d^{\prime })}$$

$${\bar{\mathbf{C}}}_d^{(d-d^{\prime })}={\mathbf{C}}_d^{(d-d^{\prime })}{\mathbf{G}}_d^{-1},{\bar{\mathbf{C}}}_{d^{\prime }}^{(d-d^{\prime })}={\mathbf{C}}_{d^{\prime }}^{(d-d^{\prime })}{\mathbf{G}}_{d^{\prime }}^{-1},$$

$${\bar{\mathbf{C}}}_d^{(d-d^{\prime })}={\bar{\mathbf{C}}}_{d^{\prime }}^{(d-d^{\prime })}{\mathbf{\Pi }}^{(d^{\prime }{)}^{-1}}{\mathbf{\Pi }}^{(d)}={\bar{\mathbf{C}}}_{d^{\prime }}^{(d-d^{\prime })}\bar{\mathbf{\Pi }}$$

$${\mathbf{A}}_d^{(d-d^{\prime })}={\mathbf{A}}_{d^{\prime }}^{(d-d^{\prime })}{\mathbf{\Lambda }}_1^{(d^{\prime }{)}^{-1}}{\mathbf{\Lambda }}_1^{(d)}$$

$${\mathbf{B}}_d^{(d-d^{\prime })}={\mathbf{B}}_{d^{\prime }}^{(d-d^{\prime })}{\mathbf{\Lambda }}_3^{(d^{\prime }{)}^{-1}}{\mathbf{\Lambda }}_3^{(d)}$$

$${\mathbf{C}}_d^{(d-d^{\prime })}={\mathbf{C}}_{d^{\prime }}^{(d-d^{\prime })}{\mathbf{\Lambda }}_3^{(d^{\prime }{)}^{-1}}{\mathbf{\Lambda }}_3^{(d)}$$

$${\mathbf{\Lambda }}_1^{(d)}{\mathbf{\Lambda }}_2^{(d)}{\mathbf{\Lambda }}_3^{(d)}=\mathbf{I},{\mathbf{\Lambda }}_1^{(d^{\prime })}{\mathbf{\Lambda }}_2^{(d^{\prime })}{\mathbf{\Lambda }}_3^{(d^{\prime })}=\mathbf{I}$$

$$\mathbf{A}({\mathcal{S}}_r^{(d)},:)\leftarrow {\mathbf{A}}_d,\mathbf{B}({\mathcal{S}}_c^{(d)},:)\leftarrow {\mathbf{B}}_d,\mathbf{C}({\mathcal{S}}_f^{(d)},:)\leftarrow {\mathbf{C}}_d,\forall d.$$

$$\underset{\mathbf{A},\mathbf{B},\mathbf{C}}{\mathrm{minimize}}$$

$$\text{NRE}=\frac{{\sum }_{k=1}^K\Vert \widehat{\underline{\mathbf{X}}}(:,:,k)-\underline{\mathbf{X}}(:,:,k){\Vert }_F}{{\sum }_{k=1}^K\Vert \underline{\mathbf{X}}(:,:,k){\Vert }_F}$$

Full text

I Introduction

Sampling and reconstruction of signals is a fundamental problem in signal processing. In the first half of the 20th century, Whittaker, Nyquist, Kotelnikov and Shannon [1, 2, 3, 4] laid the foundation of sampling theory. It guarantees perfect reconstruction of a signal from uniformly spaced samples, if sampling is performed at a rate at least twice the maximum frequency present in the signal. The Shannon-Nyquist theorem applies to both continuous and discrete signals. It capitalizes on the band-limitedness property and is the first, and one of the very few results, that allow perfect reconstruction of a signal under a uniform, or more generally, regular sampling process. The challenge is that applying Shannon-Nyquist sampling to wideband signals requires very high sampling rates—which can be quite costly in practice.

In the early 2000’s compressive sensing (CS) [5, 6, 7] emerged as an alternative, enabling reconstruction from a set of measurements, sampled or compressed below the Nyquist rate. CS works under two basic premises: the signal of interest must have a sparse representation in a known transform domain; and the sampling pattern should be ‘incoherent’. Under these assumptions, tractable algorithms are shown to recover the signal of interest. Compared to the Shannon-Nyquist sampling theorem, CS leverages signal sparsity, rather than bandlimitedness. This result is significant, since some wideband signals of practical interest, are sparse in certain domains [8, 9]. On the downside, CS entails higher reconstrunction complexity than sinc function interpolation, and relies on incoherent/random sampling – thus losing the simplicity of regular/uniform sampling. A few exceptions exist, e.g., [10, 11], but the results are probabilistic and/or usually quite restrictive in practice.

Following the ideas of CS, low-rank matrix completion (LRMC) techniques were proposed for reconstructing matrix signals from a set of samples [12, 13]. This line of research utilizes the rank of the matrix as complexity measure for sampling and has attracted significant attention, since it is related to a number of important applications such as recommender systems [14]. However, similar to CS, LRMC is based on incoherent sampling. Furthermore the reconstruction guarantees in both CS and LRMC are probabilistic, contrary to the Shannon-Nyquist theorem which deterministically guarantees signal reconstruction.

Our work is motivated by the following question. Is there a sub-Nyquist sampling mechanism that works under regular sampling for certain signals of interest? This research question is very intriguing: regular sampling is efficient, friendly to implementation and often mandatory, and sub-Nyquist sampling is desired since numerous real-world signals are far from being bandlimitted.

In this work, we offer an affirmative to the above research question for a large variety of multidimensional signals. The proposed approach guarantees recoverability under realistic conditions, both generic and deterministic. Specifically, we focus on the problem of sampling tensor signals, i.e., signals whose entries are indexed by three or more coordinates [15, 16]. Tensor signals naturally arise in a large number of areas such as machine learning and data analytics [17], signal processing and communications [18], image processing and remote sensing [19, 20, 21], medical imaging [22], genomics [23], chemometrics [24], just to name a few. Hence, considering sampling and reconstruction of tensor signals is of broad interest. The problem is challenging, since various tensor signals are neither bandlimited, sparse, nor low-rank matrices (via ‘unfolding’)—and thus existing sampling techniques are not always applicable.

The reconstruction of sampled tensor signals, known in the literature as tensor completion, has been studied in machine learning and computer vision [25, 26]. The majority of existing works [26, 27, 28, 29, 30] focus on the algorithmic aspect of tensor completion. The few that provide recovery guarantees [31, 32] are based on random sampling schemes and/or LRMC ideas, which are not tailored to the tensor specifics. The work that is closest to ours is [33], which offers reconstruction conditions when the tensor ‘fibers’ are sampled. However, the conditions are somehow restrictive, since the rank is constrained to be lower than the fiber dimension, and a variety of other interesting types of regular tensor sampling have not been considered.

Contributions: In this work, we study the task of sampling and reconstruction of signals that are tensors—or tensor sampling in short. We propose a tensor sampling framework that is flexible and easy to implement. Generic as well as deterministic theoretical conditions are derived, under which reconstruction is guaranteed. Similar to matrix completion, the sample complexity for tensor signal reconstruction is mainly affected by the canonical polyadic rank and the tensor size—instead of signal bandwidth or sparsity. Unlike CS and LRMC, the proposed approach does not require incoherent sampling. Therefore, regular, equispaced and highly structured sampling strategies can be adopted—which has a much broader spectrum of applications in practice.

Our second major contribution lies in designing accelerated acquisition schemes for functional magnetic resonance imaging (fMRI) utilizing the proposed tensor sampling principles. Note that traditional fMRI acquisition is considered an “agonizingly slow” scanning process, which strongly motivates exploring appropriate sampling techniques for acceleration. However, due to hardware limitations, random or incoherent sampling strategies are considered impractical for this task [34]. Nevertheless, the proposed tensor sampling framework fits this task very well as fMRI signals are naturally tensors. Extensive simulations using synthetically generated data show that the proposed tensor sampling schemes are promising. More importantly, experiments using real fMRI data demonstrate remarkable acceleration compared to traditional fMRI scanning approaches, without sacrificing reconstruction accuracy.

An early version of part of this work appears in conference form in Proc. ICASSP 2019 [35]. In this journal version, we consider additional sampling schemes, include deterministic reconstruction conditions along with thorough model analysis and proofs. We also design multi-slice fMRI acceleration schemes and conduct detailed experiments.

II Tensor Algebra Preliminaries

In this work we heavily use tensor algebra. To facilitate the upcoming discussion we briefly present some preliminary tensor algebra concepts. The reader is referred to [15, 16] for further details.

A third-order tensor $\underline{\bm{X}}\in\mathbb{C}^{I\times J\times K}$ is a three-way array indexed by $i,j,k$ with elements $\underline{\bm{X}}(i,j,k)$. It consists of three modes: columns $\underline{\bm{X}}(i,:,k)$, rows $\underline{\bm{X}}(:,j,k)$, fibers $\underline{\bm{X}}(i,j,:)$; and three types of slabs: horizontal $\underline{\bm{X}}(i,:,:)$, vertical $\underline{\bm{X}}(:,j,:)$ and frontal $\underline{\bm{X}}(:,:,k)$ – see Fig. 1, 2, respectively. A rank-one tensor $\underline{\bm{Z}}\in\mathbb{C}^{I\times J\times K}$ is the outer product of three vectors, $\bm{a}\in\mathbb{C}^{I},~{}\bm{b}\in\mathbb{C}^{J},~{}\bm{c}\in\mathbb{C}^{K}$, denoted as $\underline{\bm{Z}}=\bm{a}\circ\bm{b}\circ\bm{c}$, where $\circ$ is the outer product operator. Any tensor can be realized as a sum of three way outer products (rank one tensors), i.e.

[TABLE]

The above expression is known as the polyadic decomposition (PD) of a third-order tensor. If $F$ denotes the minimum number of outer products needed to synthesize $\underline{\bm{X}}$, then $F$ is called tensor rank or CP rank and the decomposition is known as canonical polyadic decomposition (CPD) or parallel factor analysis (PARAFAC) [36]. The CPD elementwise representation can be written as $\underline{\bm{X}}(i,j,k)=\sum_{f=1}^{F}{\bm{A}}(i,f){\bm{B}}(j,f){\bm{C}}(k,f),$ where $\bm{A}=[\bm{a}_{1},\dots,\bm{a}_{F}]\in\mathbb{C}^{I\times F},~{}\bm{B}=[\bm{b}_{1},\dots,\bm{b}_{F}]\in\mathbb{C}^{J\times F},~{}\bm{C}=[\bm{c}_{1},\dots,\bm{c}_{F}]\in\mathbb{C}^{K\times F}$ are called the low rank factors of the tensor. A third-order tensor can be fully characterized by its latent factors, thus we adopt the notation \underline{\bm{X}}=\left\llbracket{\bm{A}},{\bm{B}},{\bm{C}}\right\rrbracket to represent the tensor.

A striking property of tensors is that the CPD is essentially unique under mild conditions. A generic result on the uniqueness of the CPD is as follows.

Theorem 1

[37]** Let \underline{\bm{X}}=\left\llbracket{\bm{A}},{\bm{B}},{\bm{C}}\right\rrbracket with $\bm{A}:I\times F$, $\bm{B}:J\times F$, and $\bm{C}:K\times F$. Assume that ${\bm{A}}$, ${\bm{B}}$ and ${\bm{C}}$ are drawn from some joint absolutely continuous distribution. Also assume $I\geq J\geq K$ without loss of generality. If $F\leq 2^{\lfloor\log_{2}J\rfloor+\lfloor\log_{2}K\rfloor-2}$, then the decomposition of $\underline{\bm{X}}$ in terms of $\bm{A},\bm{B}$, and $\bm{C}$ is essentially unique, almost surely.

Here, essential uniqueness means that if $\tilde{\bm{A}},\tilde{\bm{B}},\tilde{\bm{C}}$ also satisfy $\underline{\bm{X}}=\llbracket\tilde{\bm{A}},\tilde{\bm{B}},\tilde{\bm{C}}\rrbracket$, then ${\bm{A}}=\tilde{\bm{A}}{\bm{\Pi}}{\bm{\Lambda}}_{1}$, ${\bm{B}}=\tilde{\bm{B}}{\bm{\Pi}}{\bm{\Lambda}}_{2}$, and ${\bm{C}}=\tilde{\bm{C}}{\bm{\Pi}}{\bm{\Lambda}}_{3}$, where ${\bm{\Pi}}$ is a permutation matrix and ${\bm{\Lambda}}_{i}$ is a full rank diagonal matrix such that ${\bm{\Lambda}}_{1}{\bm{\Lambda}}_{2}{\bm{\Lambda}}_{3}={\bm{I}}$.

As far as deterministic identifiability is concerned, we have:

Theorem 2

[38]** Let \underline{\bm{X}}=\left\llbracket{\bm{A}},{\bm{B}},{\bm{C}}\right\rrbracket with $\bm{A}:I\times F$, $\bm{B}:J\times F$, and $\bm{C}:K\times F$. The decomposition \underline{\bm{X}}=\left\llbracket{\bm{A}},{\bm{B}},{\bm{C}}\right\rrbracket is essentially unique with CP rank $F$ if $k_{\bm{A}}+k_{\bm{B}}+k_{\bm{C}}\geq 2F+2$.

Here $k_{\bm{A}}$ denotes the Kruskal rank of a matrix, i.e., the largest integer $k_{\bm{A}}$ such that any $k_{\bm{A}}$ columns of $\bm{A}$ are linearly independent.

A tensor can be represented in a matrix form using the matricization operation. There are three common ways to matricize (or unfold) a third-order tensor, by stacking columns, rows, or fibers of the tensor to form a matrix. For example the following operation stacks the fibers of tensor $\underline{\bm{X}}$ as rows of matrix $\bm{X}^{(3)}$:

[TABLE]

where ‘${\rm vec}(\cdot)$’ is the vectorization operator. One can see that:

[TABLE]

where $\odot$ denotes the Khatri-Rao (column-wise Kronecker) product. The superscript $(3)$ denotes that the unfolding is performed on the third mode of the tensor, i.e. fibers are stacked together.

Finally we will need the mode product operation in tensor analytics. The mode product operator multiplies a matrix to a tensor in a single mode. A third order tensor has three modes (rows, columns, fibers), thus three different mode products are defined. A joint mode-1, mode-2, and mode-3 product of a third-order tensor is represented by the following notation:

[TABLE]

where “$\times_{1}$” denotes the operation that multiplies each column of $\underline{\bm{X}}$ with ${\bm{P}}_{1}$, “$\times_{2}$” denotes multiplying each row of $\underline{\bm{X}}$ with ${\bm{P}}_{2}$, and “$\times_{3}$” denotes multiplying each fiber of $\underline{\bm{X}}$ with ${\bm{P}}_{3}$. The mode product is reflected in the polyadic decomposition of the tensor, i.e., the outcome of (3) results in a tensor $\underline{\tilde{\bm{X}}}$ with polyadic decomposition:

[TABLE]

The above decomposition is essentially unique under some conditions—this point will turn out to be crucial in the upcoming discussion.

III Tensor Sampling Mechanisms

The core of this work discusses the sampling and reconstruction of third-order tensors. The main claim is fundamental: roughly speaking, any third-order tensor that does not have very high rank can always be recovered from a sufficient number of regular samples. The sampling is not constrained to follow a randomized or incoherent process. On the contrary, we focus on regular and highly structured schemes. Various regular sampling strategies are considered. They involve sampling whole slabs in different modes (slab sampling), certain fibers in a single or multiple modes (fiber sampling) and entries in a systematic manner (entry sampling). Exposition and development use third-order tensors, but all the techniques can be naturally extended to higher-order tensors in a conceptually straightforward way. Similar to the case of matrices, even if a tensor is high-rank in the strict mathematical sense, it can often be approximated using low rank, in which case it can be approximately recovered using the proposed sampling and reconstruction schemes, as we will see.

III-A General Strategy and Insight

Let us consider the following general form of tensor sampling:

[TABLE]

where ${\sf Sample}(\cdot):\mathbb{C}^{I\times J\times K}\rightarrow\mathbb{C}^{L}$ is a down-sampling operator with $L\ll IJK$. Our goal is to study under what conditions and sampling strategies, recovering $\underline{\bm{X}}$ from $\bm{y}$ is possible. This is an inverse problem like in CS [5, 6, 7] and LRMC [12, 13]. However, unlike in [5, 6, 7, 12, 13, 27], we do not consider random/incoherent down-sampling operators but highly structured ones—which model a plethora of engineering applications, are easier for practical system implementation and computationally more efficient.

Our work rests upon two basic ideas. The first utilizes the uniqueness property of the CPD. Recall that every tensor admits a CPD, and the CPD is essentially unique if the CP rank is not very large (in many cases: not full-rank). The second exploits the relation between a sampled sub-tensor and the original tensor $\underline{\bm{X}}=\llbracket{\bm{A}},{\bm{B}},{\bm{C}}\rrbracket$. If we sample ${\cal S}_{r}\subseteq\{1,\ldots,I\}$ rows, ${\cal S}_{c}\subseteq\{1,\ldots,J\}$ columns and ${\cal S}_{f}\subseteq\{1,\ldots,K\}$ fibers and form a sub-tensor $\underline{\bm{X}}({\cal S}_{r},{\cal S}_{c},{\cal S}_{f})$, then:

[TABLE]

One key observation is that the above sub-tensor can be decomposed to a sum of rank-one terms of number equal to the rank of the original tensor. Furthermore, the latent factors share certain rows with the original latent factors. Intuitively, if ${\rm rank}(\underline{\bm{X}})$ is not huge, there is a good chance that the sub-tensor admits a unique CPD, and part of the information of $\bm{A}$, $\bm{B}$, and $\bm{C}$ can be extracted from the sub-tensor. Hence, by judiciously sampling and constructing sub-tensors, it seems viable to recover the entire $\bm{A}$, $\bm{B}$, and $\bm{C}$, and thus reconstruct $\underline{\bm{X}}$. This is the main idea.

Despite this conceptual simplicity, however, fleshing out this task is nontrivial. First, when factoring the sub-tensors, there are always permutation and scaling ambiguities—even if every sub-tensor admits unique CPD, reconstruction of the whole tensor is not guaranteed. Thus the sampling mechanisms need to be carefully designed to address this issue. Second, balancing the sampling ratio with the ability to identify the original tensor is a key consideration and needs attentive thinking and design. In the remaining section, we propose a series of sampling mechanisms that take into consideration both design challenges. The considered sampling schemes are practical and motivated by real engineering applications, particularly in the field of medical imaging.

III-B Slab sampling

First, we study the task of reconstructing a third-order tensor from slab samples, taken from two different modes. Recovering tensor signals from sampled slabs finds applications in fMRI acceleration [39, 40, 41] and image/video inpainting [42, 43]. However, there is no unified characterization for recoverability under regular sampling patterns, to our best knowledge.

Let $\underline{\bm{X}}\in\mathbb{C}^{I\times J\times K}$ be the original full tensor, which is not fully accessible or is subject to sampling. Instead we sample/observe a subset of slabs in one mode, e.g., horizontal slabs, $\mathcal{S}_{r}\subseteq\{1,\ldots,I\},$ and a subset of slabs in a different mode, e.g. frontal slabs, $\mathcal{S}_{f}\subseteq\{1,\ldots,K\}$. If $|\mathcal{S}_{r}|=I_{1}\geq 2$ and $|\mathcal{S}_{f}|=K_{2}\geq 2$, two separate sampled tensors are formed, i.e., $\underline{\bm{Y}}_{1}\in\mathbb{C}^{I_{1}\times J\times K}$ and $\underline{\bm{Y}}_{2}\in\mathbb{C}^{I\times J\times K_{2}}$, which represent the subset of observable horizontal and frontal slabs of $\underline{\bm{X}}$ respectively. Apparently, $\underline{\bm{Y}}_{1}$ can be written as the mode $1$ multiplication of tensor $\underline{\bm{X}}$ with selection matrix $\bm{P}_{1}^{(1)}\in\mathbb{R}^{I_{1}\times I}$, i.e.

[TABLE]

and $\underline{\bm{Y}}_{2}$ as a mode $3$ multiplication with matrix $\bm{P}_{3}^{(2)}\in\mathbb{R}^{K_{2}\times K}$, i.e.

[TABLE]

The sampling matrices $\bm{P}_{1}^{(1)},~{}\bm{P}_{3}^{(2)}$ perform slab selection in a single mode of $\underline{\bm{X}}$, thus $I_{1}<I,~{}K_{2}<K$ (they are ‘fat’) and also have full row rank. A schematic illustration of the tensor slab sampling model is given in Fig. 3. Note that, $\bm{P}_{1}^{(1)},~{}\bm{P}_{3}^{(2)}$ are not constrained to be randomly drawn in our framework. On the contrary, the sampling process is allowed to be regular or highly structured, see Fig. 3. Assuming \underline{\bm{X}}=\left\llbracket{\bm{A}},{\bm{B}},{\bm{C}}\right\rrbracket, following (4), (5), the sub-tensors $\underline{\bm{Y}}_{1},~{}\underline{\bm{Y}}_{2}$ can be expressed in a PD form:

[TABLE]

Using (6) identifiability of $\underline{\bm{X}}$ from $({\bm{Y}}_{1},{\bm{Y}}_{2})$ can be established:

Theorem 3

Let $\underline{\bm{X}}\in\mathbb{C}^{I\times J\times K}$ be the original tensor signal to recover, with CPD \underline{\bm{X}}=\left\llbracket{\bm{A}},{\bm{B}},{\bm{C}}\right\rrbracket of rank $F$. Assume that ${\bm{A}}$, ${\bm{B}}$ and ${\bm{C}}$ are drawn from a joint absolutely continuous distribution over $\mathbb{C}^{(I+J+K)F}$, and that ${\bm{A}}^{\star},{\bm{B}}^{\star},{\bm{C}}^{\star}$ satisfy the equations in (6). Then, $\underline{\hat{\bm{X}}}(i,j,k)=\sum_{f=1}^{F}{\bm{A}}^{\star}(i,f){\bm{B}}^{\star}(j,f){\bm{C}}^{\star}(k,f)$ recovers the ground-truth $\underline{\bm{X}}$ almost surely if $\min\left\{{I_{1}J},{JK},{I_{1}K},4JK_{2}\right\}\geq 4\lceil F\rceil$ or $\min\left\{{IJ},{JK_{2}},{IK_{2}},4I_{1}J\right\}\geq 4\lceil F\rceil$; where $\lceil x\rceil$ denotes the closest power of $2$ which is greater than $x$.

The proof is presented in Appendix A. The intuition is that if $\underline{\bm{Y}}_{1}$ or $\underline{\bm{Y}}_{2}$ admit a unique CPD, under Theorem 1, the factors $\bm{B},~{}\bm{C}$ or $\bm{A},~{}\bm{B}$ respectively can be identified. Then $\bm{A}$ or $\bm{C}$ are recovered from the other tensor, where $\bm{A}$ or $\bm{C}$ have been left uncompressed. Note that in slab sampling only one sub-tensor is required to admit a unique CPD. The reason is that identifying the latent factors of one sub-tensor, directly estimates two original latent factors. Then, the remaining factor can be obtained via solving a linear system of equations. Furthermore, permutation and scaling ambiguities are automatically resolved, since $\underline{\bm{Y}}_{1},~{}\underline{\bm{Y}}_{2}$ sample common rows of $\underline{\bm{X}}$. Overall reconstruction is performed as $\underline{\hat{\bm{X}}}(i,j,k)=\sum_{f=1}^{F}{\bm{A}}^{\star}(i,f){\bm{B}}^{\star}(j,f){\bm{C}}^{\star}(k,f)$. Deterministic conditions can also be derived, and this discussion is postponed to the following section.

The previous analysis can be easily extended to the case where slab sampling is performed in all 3 modes of the tensor.

III-C Fiber sampling

Next, we consider the reconstruction of tensor $\underline{\bm{X}}$ from a subset of fibers, sampled along a single mode of the tensor. Fiber sampling is also of interest to a number of applications in chemometrics [29, 44] nuclear magnetic resonance spectroscopy [45] and fMRI acceleration (see Sec. V). To make the discussion concrete, consider the scenario illustrated in Fig. 4, where $D=3$ fiber patterns appear in the tensor sampling scheme.

A pattern will be defined as a subset of rows $\mathcal{S}_{r}^{(d)}\subseteq\{1,\ldots,I\}$ and columns $\mathcal{S}_{c}^{(d)}\subseteq\{1,\ldots,J\}$ for which every point $\underline{\bm{X}}(i,j,k),~{}i\in\mathcal{S}_{r},~{}j\in\mathcal{S}_{c},~{}k\in\left\{1,\dots,K\right\}$ belongs to the pattern. In the illustrated scenario, each pattern (blue, d=1; red, d=2; green, d=3) samples fibers defined by the following subset of rows and columns: $\mathcal{S}_{r}^{(1)}=\{1,4,7,10\},~{}\mathcal{S}_{c}^{(1)}=\{1,4,7\},~{}\mathcal{S}_{r}^{(2)}=\{2,5,8,11\},\mathcal{S}_{c}^{(2)}=\{2,5,7,8\},~{}\mathcal{S}_{r}^{(3)}=\{3,6,9,12\},\mathcal{S}_{c}^{(3)}=\{3,6,8\}$. Rearranging the order of the columns results in the model shown in Fig. 5.

In the general case, the proposed fiber sampling framework entails each pattern forming a third-order tensor, i.e., $|\mathcal{S}_{r}^{(d)}|,~{}|\mathcal{S}_{c}^{(d)}|\geq 2$ and that samples are taken from every row and column of the tensor. The latter is a necessary condition for every factorization based completion approach, since completely unobserved slabs are impossible to recover. Furthermore, each pattern is required to sample from a common row or column with at least one more, thus creating an overlapping chain between patterns. The reason is that for pairwise mutually exclusive patterns, there exists a non-trivial scaling ambiguity, which cannot be determined. Formally the necessary sampling rules, for the proposed fiber sampling framework, are expressed as:

[TABLE]

where $d,d^{\prime}\in\left\{1,\dots D\right\}$. The rules in 7 handle a plethora of sampling schemes. Specifically, each pattern is allowed to be equispaced, regular, random etc. This shows that reconstruction from regular samples is indeed doable. The sampling in Fig. 4, 5, for example, is regular and each pattern consists of equispaced rows and deterministically spaced columns.

Following similar analysis as in slab sampling, let $\underline{\bm{Y}}_{d}\in\mathbb{C}^{I_{d}\times J_{d}\times K}$ be the sampled subtensor, formed by pattern $d$. Also let ${\bm{P}}_{1}^{(d)}\in\mathbb{R}^{I_{d}\times I},~{}{\bm{P}}_{2}^{(d)}\in\mathbb{R}^{J_{d}\times J}$ be the row and column selection matrices determining the $d$ pattern. Then $\underline{\bm{Y}}_{d}$ is written as follows.

[TABLE]

Using the equation in (III-C) we can establish generic identifiability of fiber sampling as:

Theorem 4

Let $\underline{\bm{X}}\in\mathbb{C}^{I\times J\times K}$ be the original tensor signal, fiber sampled according to (7), with CPD \underline{\bm{X}}=\left\llbracket{\bm{A}},{\bm{B}},{\bm{C}}\right\rrbracket of rank $F$. Assume that ${\bm{A}}$, ${\bm{B}}$ and ${\bm{C}}$ are drawn from a joint absolutely continuous distribution over $\mathbb{C}^{(I+J+K)F}$, and that ${\bm{A}}^{\star},{\bm{B}}^{\star},{\bm{C}}^{\star}$ satisfy the equations in (III-C). Then, $\underline{\hat{\bm{X}}}(i,j,k)=\sum_{f=1}^{F}{\bm{A}}^{\star}(i,f){\bm{B}}^{\star}(j,f){\bm{C}}^{\star}(k,f)$ recovers the ground-truth $\underline{\bm{X}}$ almost surely if $\min_{d}\left\{{I_{d}J_{d}},{J_{d}K},{I_{d}K}\right\}\geq 4\lceil F\rceil$; where $\lceil x\rceil$ denotes the closest power of $2$ which is greater than $x$.

The proof is relegated to Appendix B. In contrast to the previous case of slab sampling, where identifiability of one sampled tensor $\underline{\bm{Y}}_{i}$ is enough, fiber sampling requires all $\underline{\bm{Y}}_{i}$’s to admit a unique CPD model—otherwise certain rows of $\bm{A},~{}\bm{B}$ would be impossible to identify. The claim is simple and intuitive: The number of samples required to recover a fiber sampled tensor are proportional to the rank of the tensor.

Remark 1

Theorem 4 studies general tensors where factor $\bm{C}$ is not required to have full column rank, and thus $K<F$ can be easily handled. Fiber sampling and recovery of tensors with $\bm{C}$ having full column rank, is extensively studied in [33]. Compared to our work, the sampling strategy therein has to follow rules (7b), (7c), whereas (7a) can be relaxed. On the other hand, the results of this paper are tailored to cases where the sampling process exhibits some regularity and $K<F$ is allowed. Note that $\bm{C}$ being full column rank, which is mandatory in [33], is a quite restrictive condition and prohibitive for several applications, e.g., fMRI acceleration as we will see next.

III-D Entry sampling

So far we have discussed slab, fiber sampling of third-order tensors and provided conditions under which reconstruction is guaranteed. In this subsection, we move a step further and study the more general problem of tensor reconstruction from a subset of entries, sampled in a regular fashion along the tensor. Entry sampling is another important sampling mechanism, which along with fiber sampling will prove very useful in accelerating the fMRI scan acquisition (see Sec. V).

We are interested in cases where the sampling process can be viewed as a series of patterns. A pattern is defined, similarly to fiber sampling, as a subset of rows $\mathcal{S}_{r}^{(d)}\subseteq\{1,\ldots,I\},$ columns $\mathcal{S}_{c}^{(d)}\subseteq\{1,\ldots,J\}$ and fibers $\mathcal{S}_{f}^{(d)}\subseteq\{1,\ldots,K\}$, for which every point $\underline{\bm{X}}(i,j,k),~{}i\in\mathcal{S}_{r},j\in\mathcal{S}_{c},k\in\mathcal{S}_{f}$ belongs to the pattern. For example consider the scenario illustrated in Fig. 6.

The number of patterns is $D=3$ and $\mathcal{S}_{r}^{(1)}=\mathcal{S}_{c}^{(1)}=\{1,3,5,7\},~{}\mathcal{S}_{f}^{(1)}=\{1,4\},~{}\mathcal{S}_{r}^{(2)}=\mathcal{S}_{c}^{(2)}=\{2,4,6,8\},\mathcal{S}_{f}^{(2)}=\{2,5\},~{}\mathcal{S}_{r}^{(3)}=\mathcal{S}_{c}^{(3)}=\{3,4,5,6\},\mathcal{S}_{f}^{(3)}=\{3,6\}$. In general, the proposed framework requires samples to be taken from all rows, columns and fibers of the tensor, thus $D\geq 2$ and each pattern should include $2$ rows, $2$ columns and $2$ fibers at minimum. Furthermore, the patterns need to overlap like a domino, i.e. any pattern should sample from some common rows, columns or fibers with another. Furthermore, the overlap between two patterns is required to involve 2 elements in one mode and 1 element in a different mode. For example, a pair of patterns should sample 2 common rows and 1 common column. The latter is a necessary condition resulting from the inherent permutation and scaling ambiguity of the CPD. Formally the rules of entry sampling are:

[TABLE]

where $d,d^{\prime}\in\left\{1,\dots D\right\},~{}m,m^{\prime}\in\left\{c,r,f\right\}$. Following similar analysis as in single mode fiber sampling, let $\underline{\bm{Y}}_{d}\in\mathbb{C}^{I_{d}\times J_{d}\times K_{d}}$ be the sampled subtensor representation of pattern $d$. Also let ${\bm{P}}_{1}^{(d)}\in\mathbb{R}^{I_{d}\times I},~{}{\bm{P}}_{2}^{(d)}\in\mathbb{R}^{J_{d}\times J},~{}{\bm{P}}_{3}^{(d)}\in\mathbb{R}^{J_{d}\times J}$ be the row, column and fiber selection matrices determining pattern $d$. Then $\underline{\bm{Y}}_{d}$ is written as:

[TABLE]

The model in (III-D) is identifiable, under generic conditions presented in the following theorem.

Theorem 5

Let $\underline{\bm{X}}\in\mathbb{C}^{I\times J\times K}$ be the original tensor signal, sampled according to (9), with CPD \underline{\bm{X}}=\left\llbracket{\bm{A}},{\bm{B}},{\bm{C}}\right\rrbracket of rank $F$. Assume that ${\bm{A}}$, ${\bm{B}}$ and ${\bm{C}}$ are drawn from a joint absolutely continuous distribution over $\mathbb{C}^{(I+J+K)F}$, and that ${\bm{A}}^{\star},{\bm{B}}^{\star},{\bm{C}}^{\star}$ satisfy the equations in (III-D). Then, $\underline{\hat{\bm{X}}}(i,j,k)=\sum_{f=1}^{F}{\bm{A}}^{\star}(i,f){\bm{B}}^{\star}(j,f){\bm{C}}^{\star}(k,f)$ recovers the ground-truth $\underline{\bm{X}}$ almost surely if $\min_{d}\left\{{I_{d}J_{d}},{J_{d}K_{d}},{I_{d}K_{d}}\right\}\geq 4\lceil F\rceil$; where $\lceil x\rceil$ denotes the closest power of $2$ which is greater than $x$.

The proof is presented in Appendix B. Similar to fiber sampling, reconstruction of a tensor from entries, sampled as described in (9), is guaranteed, if all the sub-sampled tensors formed by the emerging patterns admit a unique CPD.

IV Deterministic Identifiability and insights

The sampling mechanisms, discussed so far, can be realized as separate, yet coupled, sub-sampled versions of the original third-order tensor $\underline{\bm{X}}$. Recovery of $\underline{\bm{X}}$, under various sampling mechanisms, was established by applying generic identifiability results on the CPD of the sub-tensors. However, reconstruction of the original tensor is also guaranteed under purely deterministic conditions. In the case of slab sampling we have the following theorem.

Theorem 6

Let $\underline{\bm{X}}\in\mathbb{C}^{I\times J\times K}$ be the original tensor signal to recover, with CPD \underline{\bm{X}}=\left\llbracket{\bm{A}},{\bm{B}},{\bm{C}}\right\rrbracket of rank $F$. Assume that ${\bm{A}}^{\star},{\bm{B}}^{\star},{\bm{C}}^{\star}$ satisfy the equations in (6). Then, $\underline{\hat{\bm{X}}}(i,j,k)=\sum_{f=1}^{F}{\bm{A}}^{\star}(i,f){\bm{B}}^{\star}(j,f){\bm{C}}^{\star}(k,f)$ recovers the ground-truth $\underline{\bm{X}}$ if $2F+2\leq k_{{\bm{P}_{1}^{(1)}\bm{A}^{\star}}}+k_{\bm{B}^{\star}}+k_{\bm{C}^{\star}}$ and $\bm{B}^{\star}\odot{\bm{P}}_{3}^{(2)}\bm{C}^{\star}$ has full column rank, or if $2F+2\leq k_{{\bm{A}^{\star}}}+k_{{\bm{B}^{\star}}}+k_{\bm{P}_{3}^{(2)}\bm{C}^{\star}}$ and $\bm{B}^{\star}\odot{\bm{P}}_{1}^{1}\bm{A}^{\star}$ has full column rank.

When fiber or entry sampling is employed, we have:

Theorem 7

Let $\underline{\bm{X}}\in\mathbb{C}^{I\times J\times K}$ be the original tensor signal, fiber or entry sampled according to (7) or (9) respectively. Also let \underline{\bm{X}}=\left\llbracket{\bm{A}},{\bm{B}},{\bm{C}}\right\rrbracket denote the rank-$F$ CPD of $\underline{\bm{X}}$. Assume that ${\bm{A}}^{\star},{\bm{B}}^{\star},{\bm{C}}^{\star}$ have no repeated entries and satisfy the equations in  (III-C), (III-D), according to the sampling mechanism. Then, $\underline{\hat{\bm{X}}}(i,j,k)=\sum_{f=1}^{F}{\bm{A}}^{\star}(i,f){\bm{B}}^{\star}(j,f){\bm{C}}^{\star}(k,f)$ recovers the ground-truth $\underline{\bm{X}}$ if $~{}2F+2\leq\min_{d}\left\{k_{{\bm{P}_{1}^{(d)}\bm{A}^{\star}}}+k_{{\bm{P}_{2}^{(d)}\bm{B}^{\star}}}+k_{{\bm{P}_{3}^{(d)}\bm{C}^{\star}}}\right\}$, where ${\bm{P}}_{3}^{(d)}=\bm{I}$ for fiber sampling.

Proof of both theorems is presented in Appendix C.

The implication of Theorems 3- 7 is significant and intuitive. Recovery of $\underline{\bm{X}}$ is based on two basic principles: identifiability of the factors of the sub-sampled tensors and ability to reconcile for the permutation and scaling ambiguities. The first is a property of the both the signal of interest and the sampling mechanism. In particular the rank of the tensor signal, along with Kruskal or generic conditions on the factors determine the number of samples required for full tensor recovery. Note that the number and structure of samples varies according to the applied sampling mechanism. Hence, there is a clear correlation between the rank of the tensor and the number of samples needed—higher ranks require higher number of samples. The second principle is a necessary property of the sampling mechanism. Although slab sampling automatically handles permutation and scaling ambiguities, fiber and entry sampling schemes have to be carefully designed to satisfy (7) or (9) and eliminate permutation and scaling mismatches.

In a nutshell, one can learn the rows of factors $\bm{A},~{}\bm{B},~{}\bm{C}$ from the sub-tensors, up to column permutation and scaling and resolve the mismatches using common information between the sub-tensors. Then reconstruction of the original tensor is attained as \underline{\bm{X}}=\left\llbracket{\bm{A}},{\bm{B}},{\bm{C}}\right\rrbracket.

To build some further intuition on the theoretical conditions, consider the following example. Let $\underline{\bm{X}}\in\mathbb{C}^{512\times 512\times 513}$ be the tensor with CP rank $F=1000$ which is subject to sampling. First we sample the $I_{1}$ equispaced horizontal slabs and $K_{2}$ equispaced frontal slabs. Following Theorem 3, reconstruction of $\underline{\bm{X}}$ is guaranteed if we sample at least $I_{1}=8$ horizontal slabs and $K_{2}=2$ frontal slabs and vise versa. This results in sampling ratio $r=\frac{\#\text{observed entries}}{\#\text{total entries}}=0.019$. Next we sample fibers of the tensor in a regular fashion, similarly to Fig. 4. According to theorem 4, $33,280$ fibers are sufficient to recover the original tensor, which gives sampling ratio $r=0.13$. Finally, entries are sampled in a regular fashion, as shown in Fig. 6. The total number of entries required to reconstruct the original tensor is $2,097,664$, according to theorem 5 , which results in sampling ratio $r=0.016$. Note that for smaller rank, e.g, $F=250$ the total number of samples can be significantly reduced, giving sampling ratio $r_{slab}=0.008$, $r_{fiber}=0.063$, $r_{entry}=0.004$.

V Application to parallel fMRI acceleration

Interestingly, the previously described sampling mechanisms find application in accelerating fMRI scan acquisition. fMRI is used to measure brain activity associated with changes in blood oxygen levels. MRI acquisitions typically use a set of coils (sensors), that in parallel collect a series of frames focusing on different parts of the brain. In fMRI, the three-dimensional (3D) volume covering the whole brain is typically acquired using multiple two-dimensional (2D) slices. These are discrete-space signals, sampled along a particular trajectory in the k-space, which is a 2-D frequency domain ($k_{x},~{}k_{y}$), for each brain slice. Therefore an fMRI scan, can be represented by a five-way array with coil, $k_{x},~{}k_{y}$, slice and time (frame) dimensions.

Acquiring high spatial resolution fMRI is challenging due to time restrictions. On the one hand, the $k$-space sampling has to follow the Shannon-Nyquist theorem to avoid artifacts, when inverse Fourier transform is used for reconstruction. On the other hand sampling at a Nyquist rate leads to prolonged scan acquisition time for each frame, which is prohibitive for high temporal resolution, required in fMRI and neuroscience research. The objective is therefore two-fold: accelerate the scanning process and capture fast brain activity changes. Since the scan acquisition time is proportional to the number of $k$-space samples, ongoing efforts focus on sampling part of the $k$-space ($k_{y}$ frequencies) of each slice and/or measuring the $k$-space of only a subset of slices. The majority of work is mainly proposed for MRI scans. Classic methods use learning and calibration type techniques [46, 47, 48], while others employ the CS framework [47, 49, 50] or LRMC [50, 51, 52, 53] to perform the reconstruction.

While MRI offers significant freedom in designing the $k$-space sampling trajectories for each frame, fMRI acquisition is more restrictive. Specifically, fMRI is performed using a special fast imaging acquisition, called echo planar imaging, that is practically only used with equispaced sub-sampling patterns due to restrictions associated with magnetic field inhomogenities and Eddy currents [34]. In simple words, the $k_{y}$ frequencies sampled for each frame have to be equispaced and all coils need to measure the same frequencies. For example the sampling scheme illustrated in Fig. 7 is typical in fMRI and performs 3-fold acceleration.

In general, a fully sampled scan is acquired first, which is beneficial for calibration purposes. Then $n$-fold acceleration is achieved by sampling $1/n$ of equispaced $k_{y}$ frequencies. The frequencies to sample for each frame can be the same for the whole procedure or can be circularly shifted as in Fig. 7. For the benefit of our method we propose to circularly shift between $n$ equispaced set of frequencies in order to capture the temporal behavior of the brain accurately.

Sampling the $k_{y}$ dimension is one way to accelerate the fMRI scanning process. Another idea that is being used is to observe the $k$-space of only a subset of slices at each time slot. In this work we propose to combine these two ideas to further reduce scanning time. Specifically at each time instance sub-sampled $k$-space measurements are acquired for only a subset of slices, instead of the complete set. To design a sampling mechanism that fits our tensor models and fMRI constraints, we first need to acquire a fully sampled scan for every slice. Then for each frame $1/r$ of equispaced $k_{y}$ frequencies is observed for $1/s$ of the brain slices, so that at the $(rs)$-th frame we have measured every frequency for every slice. This results in $(rs)$-fold acceleration. Fig. 8 illustrates an fMRI acceleration technique, where $r,s=2$.

Note that the two aforementioned sampling procedures can be tricky for classic techniques. On the one hand, calibration-based techniques such as GRAPPA [46] are linear and suffer from noise amplification at high acceleration rates. On the other hand, CS and LRMC schemes have difficulties in operating with regular samples, since their success rests upon incoherent sampling.

On the contrary, the proposed tensor sampling and reconstruction framework is exactly designed to handle these highly structured and constrained sampling schemes used in fMRI acquisitions. In particular, the single-slice fMRI acceleration task, as illustrated in Fig. 7, can be cast as a tensor fiber sampling mechanism, analyzed in subsection III-C. As mentioned earlier the raw fMRI scan is originally a five-way array and thus each slice is a four-way array. Although the previous analysis could be easily extended to tensors of order higher than three, we choose to work with third-order ones. Specifically the $k$-space is processed in a single dimension (mode) by concatenating $k_{x}$ and $k_{y}$. The reason is that the relation between $k_{x}$ and $k_{y}$ is often hard to be captured by a multilinear tensor model. As a result one fMRI slice is modeled as a third-order tensor $\underline{\bm{X}}\in\mathbb{C}^{I\times J\times K}$, where $I=m_{x}m_{y}$ with $m_{x},~{}m_{y}$ representing the number of frequencies in $k_{x},~{}k_{y}$ space respectively, $J$ represents the total number of frames (time slots) and $K$ the number of coils. Following the analysis of subsection III-C we have the following result.

Proposition 1

Let $\underline{\bm{X}}\in\mathbb{C}^{I\times J\times K}$ be the a single-slice fMRI tensor with rank $F$, modeled as previously explained. Under the assumptions of Theorem 4, $n$-fold acceleration can be achieved if $n\leq\min\left\{\sqrt{\frac{IJ}{4\lceil F\rceil}},\frac{JK}{4\lceil F\rceil},\frac{IK}{4\lceil F\rceil}\right\}$.

Similarly, the proposed multi-slice fMRI acceleration scheme, which performs joint $k$-space and slice sampling, is cast as an entry tensor sampling procedure, introduced in III-D. To do so, the $k$-space is considered as a single mode, as before, and we also concatenate coils and slices in one dimension. The resulting third-order tensor $\underline{\bm{X}}\in\mathbb{C}^{I\times J\times K}$ has the $k$-space in the first mode, i.e, $I=m_{x}m_{y}$, $J$ represents the total number of frames and the third mode includes the concatenation of coils and slices, i.e., $K=m_{s}m_{c}$ with $m_{s},~{}m_{c}$ being the number of slices and coils respectively. Following the analysis of subsection III-D we have:

Proposition 2

Let $\underline{\bm{X}}\in\mathbb{C}^{I\times J\times K}$ be the a multi-slice fMRI tensor with rank $F$, modeled as previously explained. Under the assumptions of Theorem 5, ${(rs)}$-fold acceleration can be achieved if $rs\leq\frac{1}{4\lceil F\rceil}\min\left\{IK,\frac{JK}{s},\frac{IJ}{r}\right\}$.

We should mention that tensor approaches have also been proposed in medical imaging [54, 55, 56, 57, 58]. The work in [54], for example, uses a tensor model to approach the MRI sampling and reconstruction problem in an on-line fashion. However, [54] works under different sampling schemes, which are not regular and appropriate for fMRI, and reconstruction guarantees are not discussed. Moreover, a tensor model is also used in [55], in the context of MRI denoising, which is different from the fMRI or MRI acceleration problem. Finally the works in [56, 57, 58] adopt a Tucker model [15], and also require auxiliary acquisition data to estimate the bases in non-spatial modes prior to reconstruction, which are not available in most fMRI acquisitions.

VI General Algorithmic framework for Tensor Sampling

Previously we studied sampling and reconstruction of third-order tensors under different, identifiable sampling mechanisms. In the current section, the algorithmic component of our approach is discussed. In general, reconstruction of a tensor from a subset of entries falls under the framework of tensor completion. A plethora of algorithms have been proposed, e.g. [29, 44]. The idea is to use the CPD factors, computed from the incomplete tensor, and reconstruct the original one. Popular methods approach the problem as a system of non-linear equations and handle it using descent direction approaches, such as gradient descent, alternating optimization, or the Gauss-Newton method.

Existing tensor completion works could be employed to approach the recovery task of a regularly sampled tensor. However, the unique characteristics of our models would be ignored. To put it in context, the special structure of regular sampling allows tensor completion by computing the factors of complete tensors. Note that CPD computation of a complete tensor is a considerably easier task than that of an incomplete one. Several polynomial time algebraic algorithms [59, 60] have been shown to retrieve the original factors,under certain conditions, or effectively initialize optimization approaches with significant success.

We propose a three step approach to tackle the completion task, which follows the insights of Theorems 3-7. The first step solves the CPD of the sub-sampled tensors independently, the second reconciles for permutation/scaling ambiguities and gets an initial estimate of the factors, and the third solves the coupled CPD problem. Detailed analysis follows.

VI-A Step 1: Computing the CPD of sub-tensors

First, the CPD of the sub-sampled tensors $\underline{\bm{Y}}_{i}$ is computed. This step is guided from the requirements of each sampling mechanism. To be more precise, CPD of $\underline{\bm{Y}}_{i}$ is computed if the reconstruction conditions require $\underline{\bm{Y}}_{i}$ to admit an identifiable CPD. The slab sampling model, for instance, requires only $\underline{\bm{Y}}_{1}$ or $\underline{\bm{Y}}_{2}$ to admit unique CPD. Therefore the CPD of only one sub-tensor is needed. On the other hand when fiber or entry sampling is considered, the CPD computation of every sub-tensor is performed, following the previous identifiability analysis.

VI-B Step 2: Initializing the factors

After computing the CPD of the sub-tensors, step 2 computes an initial estimate of the $\bm{A},\bm{B},\bm{C}$ factors, after resolving possible permutation and scaling mismatches. We distinguish between 2 different cases:

Case 1, slab sampling: As mentioned earlier slab sampling automatically reconciles for permutation and scaling ambiguities. Furthermore, two of the factors have been already computed from step 1 (e.g., $\bm{A}$,\bm{B}$$\leftarrowCPD($\underline{\bm{Y}}_{2}$)). What remains to be obtained is the third factor (e.g., $\bm{C}$), which is revealed by the other sub-tensor ($\underline{\bm{Y}}_{1}$), via solving a linear system of equations:

[TABLE]

Case 2, fiber and entry sampling: Contrary to slab sampling, the permutation and scaling ambiguity is an important issue when fiber or entry sampling is applied. To be more precise, let \underline{\bm{Y}}_{d}=\left\llbracket{\bm{A}}_{d},{\bm{B}}_{d},{\bm{C}}_{d}\right\rrbracket,~{}d\in\{1,\dots D\}, be the sub-tensors formed after fiber or entry sampling. Then:

[TABLE]

where ${\bm{\Pi}}^{(d)}\neq{\bm{\Pi}}^{({d^{\prime}})}$ are permutation matrices and ${\bm{\Lambda}}_{i}^{(d)}\neq{\bm{\Lambda}}_{i}^{({d^{\prime}})}$ are full rank diagonal matrices such that ${\bm{\Lambda}}_{1}^{(d)}{\bm{\Lambda}}_{2}^{(d)}{\bm{\Lambda}}_{3}^{(d)}={\bm{I}}$, $d,d^{\prime}\in\{1,\dots,D\}$. Clearly, in order to synthesize $\bm{A},\bm{B},\bm{C}$ from $\bm{A}_{d},\bm{B}_{d},\bm{C}_{d}$ and reconstruct $\underline{\bm{X}}$, the permutation and scaling mismatch should be resolved, i.e., ${\bm{\Pi}}^{(d)}={\bm{\Pi}}^{({d^{\prime}})}$, ${\bm{\Lambda}}_{i}^{(d)}={\bm{\Lambda}}_{i}^{({d^{\prime}})}$ for every $d,~{}d^{\prime}$.

To overcome this issue, the common information between sub-tensors is utilized. In simple words, (7) or (9) require $\bm{C}_{d}-\bm{C}_{d^{\prime}}$ (or $\bm{A}_{d}-\bm{A}_{d^{\prime}}$, or $\bm{B}_{d}-\bm{B}_{d^{\prime}}$) to share some common rows, i.e., $|\mathcal{S}_{f}^{(d)}\cap\mathcal{S}_{f}^{(d^{\prime})}:=\mathcal{S}_{f}^{(d-d^{\prime})}|\geq 2$. Now let:

[TABLE]

and normalize $\bm{C}_{d}^{(d-d^{\prime})},~{}\bm{C}_{d^{\prime}}^{(d-d^{\prime})}$, such that the share a common row, e.g.:

[TABLE]

where $\bm{G}_{d}=\text{diag}\left({\bm{C}}_{d}^{(d-d^{\prime})}(1,:)\right)$, $\bm{G}_{d^{\prime}}=\text{diag}\left({\bm{C}}_{{d^{\prime}}}^{(d-d^{\prime})}(1,:)\right)$ and $\text{diag}(\bm{x})$ is the diagonal matrix of row vector $\bm{x}$. Then:

[TABLE]

and we can solve for $\bar{\bm{\Pi}}$ using the Hungarian algorithm [61]. This procedure resolves the permutation mismatch between the factors, i.e., ${\bm{\Pi}}^{(d)}={\bm{\Pi}}^{({d^{\prime}})}$. Note that in case of fiber sampling $\mathcal{S}_{f}^{(d)}=\mathcal{S}_{f}^{(d^{\prime})}=\mathcal{S}_{f}^{(d-d^{\prime})}=\{1\dots K\}$.

To reconcile for the scaling ambiguity we require extra information coming from the factors that were not involved in permutation match, i.e., $\bm{A}_{d}-\bm{A}_{d^{\prime}}$ or $\bm{B}-\bm{B}_{d^{\prime}}$ in our example. The necessary rules (7c), (9c) enforce that there is at least one common row between $\bm{A}_{d}-\bm{A}_{d^{\prime}}$ or $\bm{B}-\bm{B}_{d^{\prime}}$. Then the scaling mismatch can be solved by the following set of equations:

[TABLE]

${\bm{A}}_{d}^{(d-d^{\prime})},~{}{\bm{A}}_{d^{\prime}}^{(d-d^{\prime})}$ and ${\bm{B}}_{d}^{(d-d^{\prime})},~{}{\bm{B}}_{d^{\prime}}^{(d-d^{\prime})}$ represent the common rows between $\bm{A}_{d}-\bm{A}_{d^{\prime}}$ and $\bm{B}-\bm{B}_{d^{\prime}}$ respectively. Next, an initial estimate of the factors is extracted by reading out the appropriate rows from the sub-tensor factors to synthesize $\bm{A},~{}\bm{B},~{}\bm{C}$ i.e.,

[TABLE]

VI-C Step 3: Coupled CPD

Finally, $\bm{A},~{}\bm{B},~{}\bm{C}$ are jointly computed as a classic tensor factorization problem with missing entries, which is equivalent to the following coupled CPD estimator:

[TABLE]

In slab sampling, ${\bm{P}}_{2}^{(1)},{\bm{P}}_{3}^{(1)},{\bm{P}}_{1}^{(2)},{\bm{P}}_{2}^{(2)}$ are identity matrices and in fiber sampling ${\bm{P}}_{3}^{(d)}$ is always the identity. There are several ways to handle the above non-convex problem. We choose to employ the tensorlab [62] toolbox, which uses a Gauss Newton approach to solve this nonlinear least squares (NLS) problem. After obtaining the estimates of $\bm{A},\bm{B}$ and $\bm{C}$, $\underline{\bm{X}}$ can be reconstructed by:

$\underline{\hat{\bm{X}}}(i,j,k)=\sum_{f=1}^{F}\hat{\bm{A}}(i,f)\hat{\bm{B}}(j,f)\hat{\bm{C}}(k,f)$.

VI-D REgular Tensor Sampling and INterpolation Algorithm (RETSINA)

As mentioned earlier the accelerated fMRI acquisition can be cast as a tensor sampling and reconstruction task. Therefore it falls under the class of problems that the previously described framework can handle. However we choose to follow a different initialization approach tailored to the specific application. We design two algorithms, one for single-slice fMRI and one for multi-slice fMRI.

Single slice acceleration: The REgular Tensor Sampling and INterpolation Algorithm (RETSINA) is presented in Algorithm 1. We follow a 3 step procedure. In step 1 (initialization), for $n$-fold acceleration we sum every $n$ vertical slabs (frames), where the missing $k$-space measurements are considered zeros, and obtain a tensor $\underline{\bm{X}}_{n}\in\mathbb{C}^{I\times J/n\times K}$ without missing entries. Then we compute the CPD of $\underline{\bm{X}}_{n}$ and get a rough estimate of $\bm{A},~{}\bm{C}$ factors. In step 2 (refinement), we compute the CPD of $\underline{\bm{Y}}_{1}$ (initialized by step 1) and the CPD of $\{\underline{\bm{Y}}_{i}\}_{i\neq 1}$ with known $\bm{C}$. Finally, in step 3, we compute the final factors by solving (15) with tensorlab’s Gauss-Newton algorithm. Compared to the previously presented general framework, RETSINA empirically yields enhanced reconstruction accuracy and reduces the operational time.

Multi slice acceleration: The Multi-Slice RETSINA (MS-RETSINA) is presented in Algorithm 2. Compared to RETSINA the initialization step has been modified to this specific case and the refinement step is skipped, since it has been observed that it does not improve the overall performance.

VII Simulations

In this section, we showcase the effectiveness of the proposed tensor sampling framework using numerical experiments. The experiments involve synthetically generated data as well as fMRI scans in the $k$-space. All simulations are performed in Matlab on a Linux server with 3.6GHz cores and 32GB RAM, except part C which is performed on a Linux server with 2GHz cores and 128GB RAM.

VII-A Synthetic Experiments

First synthetically generated experiments are conducted to examine the validity of our claims and the performance of the proposed framework. In particular, a tensor $\underline{\bm{X}}\in\mathbb{R}^{500\times 500\times 500}$ is generated as \underline{\bm{X}}=\left\llbracket{\bm{A}},{\bm{B}},{\bm{C}}\right\rrbracket. The elements of the factor matrices are drawn from an independent identically distributed (i.i.d.) zero mean, unit variance Gaussian distribution. We regularly sample $\underline{\bm{X}}$ according to the previously presented sampling mechanisms, i.e,. slab sampling, fiber sampling, and entry sampling, as shown in Figs. 3- 6. Specifically, for slab sampling we sample equispaced frontal and horizontal slabs. Regarding fiber sampling, each tensor $\underline{\bm{Y}}_{i}$ is a set of fibers, defined by equispaced rows and columns of $\underline{\bm{X}}$. Note that one vertical slab is fully observed to reconcile for permutation and scaling ambiguities. Equivalently entry sampling is designed to observe different sets of equispaced entries plus a fully sampled vertical slab.

For the experiments, we define the sampling ratio, ${r=\frac{\#\text{of sampled entries}}{IJK}}$, and vary it from $0.75$ to $0.001$. We also alter the tensor rank $F$ from $5$ to $1000$. To evaluate the performance of tensor reconstruction, we measure the normalized reconstruction error, i.e.

[TABLE]

Fig. 9 presents the results for the three sampling schemes.

As expected the reconstruction accuracy is deteriorating as the rank $F$ increases or the sampling ratio $r$ decreases. For reasonably small ranks and high number of samples the reconstruction is perfect.

VII-B Accelerated parallel fMRI

Next, the tensor sampling and reconstruction framework is tested in a real and important problem, that of parallel fMRI acceleration. First, we test the performance of the proposed RETSINA with fMRI scans, fully sampled in the $k$-space, obtained from the Center for Magnetic Resonance Research (CMRR) at the University of Minnesota.

The single slice raw scan is originally a fourth-order tensor of size $104\times 104\times 32\times 490$ and we unfold it as a third-order tensor $\underline{\bm{X}}\in\mathbb{C}^{10816\times 32\times 490}$. We apply 3-fold acceleration by observing $1/3$ of the $k_{y}$ frequencies, as shown in Fig. 7. We choose $F=100$ and run step 1, 2, 3 for 50, 2 and 5 iterations respectively. The baseline algorithms used for comparison are k-t Focuss [50], which is a CS type algorithm, k-t SLR [52] which combines ideas from both LRMC and CS, and the zero padding inverse discrete Fourier transform (IDFT). The performance of IDFT is an indicator on how difficult the reconstruction is. Note that k-t SLR directly reconstructs the fMRI signal in the absolute $(x-y)$-time-coil space. To be more precise it reconstructs signal $\underline{\bm{W}}=|\mathcal{Q}(\underline{\bm{X}})|$, where $\mathcal{Q}$ denotes the inverse Fourier transform from the $k_{x}-k_{y}$ to the $x-y$ space and $|\cdot|$ is the absolute value. Thus we also measure the NRE of signal $\underline{\bm{W}}$, denoted as $\text{NRE}_{2}$, for fair comparisons. For both k-t Focuss and k-t SLR the publicly available code was used. Note that k-t Focuss and k-t SLR are single coil algorithms in their original implementation, thus we treated each coil separately. k-t SLR requires parameter tuning and so we used a validation step to tune effectively.

The results are presented in Table I, which includes the NRE in $k$-space and absolute $x-y$ space as well as runtime. The proposed RETSINA achieves highest reconstruction quality in the $k$-space and works comparably well (but markedly faster) with k-t SLR in reconstructing the signal magnitude in the $x-y$ space. This is expected, since RETSINA reconstructs both the magnitude and phase (which is very important in images) in the $k$-space, whereas k-t SLR reconstructs the magnitude in the $x-y$ space. In terms of runtime RETSINA works faster than k-t SLR and k-t Focuss, but slower than IDFT. However IDFT exhibits very poor reconstruction performance. It is worth noting, that k-t Focuss and k-t SLR qualify for parallel implementation, which could potentially speed up operation time. This would be held, however, at the cost of computational resources.

Fig. 10 shows the reconstructed fMRI scans at different time frames produced by RETSINA along with the fully sampled data. The quality of the reconstruction is significantly high, rendering the proposed RETSINA a good alternative for fMRI acceleration.

Finally, in Fig. 11 we illustrate the reconstruction performance at a single frame for the competing algorithms. IDFT gives an illustration of the downsampled image, RETSINA works the best and k-t SLR work comparably well, although being slightly off in contrast.

VII-C Accelerated multi-slice parallel fMRI

Finally, the proposed framework is tested in the task of accelarated multi-slice fMRI acquisition. Recall that acceleration is performed at 2 levels, since at each time slot we measure the sampled $k$-space of only a subset of slices. The multi-slice fMRI raw scan is a fifth-order tensor of size $104\times 104\times 32\times 490\times 8$, where the number of slices is $8$. We unfold it as a third-order tensor $\underline{\bm{X}}\in\mathbb{C}^{10816\times 490\times 256}$ and observe $1/r$ of the $k_{y}$ frequencies and $1/s$ of slices, which leads to $rs$-fold acceleration, as illustrated in Fig. 8. The tensor rank used in MS-RETSINA is $F=100$. Table II shows the performance of the proposed MS-RETSINA in terms of NRE for various values of $r$ and $s$. An illustrative example of the reconstruction performance is presented in Fig. 12.

VIII Conclusion

In this paper we studied the sampling and reconstruction of tensors under various schemes. Compared to CS, LRMC, as well as other tensor works, we provide concrete conditions, deterministic and generic, under which exact reconstruction of the tensor was proven to be attainable from regular or even equispaced samples. Furthermore, we cast the fMRI acceleration task as regular tensor sampling process and provided an efficient algorithmic framework to approach the problem. Simulations with synthetic data as well as fMRI scans in the $k$-space show the validity and effectiveness of our approach.

Appendix A Proof of theorem 3

First, we adjust lemma 1 from [19] for complex numbers and selection matrices, which is essential for the proofs.

Lemma 1

[19]** Let $\tilde{\bm{Z}}=\bm{Q}\bm{Z}$, where the elements of $\bm{Z}$ are drawn from an absolutely continuous joint distribution with respect to the Lebesgue measure in $\mathbb{C}^{IF}$ and $\bm{\bm{Q}}\in\mathbb{R}^{I^{\prime}\times I}$ is a row selection matrix with full row rank. Then the joint distribution of the elements in $\tilde{\bm{Z}}$ is absolutely continuous with respect to the Lebesgue measure in $\mathbb{C}^{I^{\prime}F}$

It follows that $\bm{P}_{1}^{(1)}\bm{A},\bm{P}_{3}^{(2)}{\bm{C}}$ are drawn from non-singular absolutely continuous distributions. Then Theorem 1 determines the conditions under which the factors of $\underline{\bm{Y}}_{1}$ or $\underline{\bm{Y}}_{2}$ can be identified. Let’s consider the case where $\underline{\bm{Y}}_{1}$ is identifiable. Under the conditions of Theorem 3, \underline{\bm{Y}}_{1}=\left\llbracket{\bm{A}}_{1},{\bm{B}}_{1},{\bm{C}}_{1}\right\rrbracket is essentially unique and from (4) holds that:

[TABLE]

where ${\bm{\Pi}}$ is a permutation matrix and ${\bm{\Lambda}}_{i}$ is a full rank diagonal matrix such that ${\bm{\Lambda}}_{1}{\bm{\Lambda}}_{2}{\bm{\Lambda}}_{3}={\bm{I}}$. Recall that $\underline{\bm{Y}}_{2}$ admits a (possibly non-unique) PD \underline{\bm{Y}}_{2}=\left\llbracket\bm{A},{\bm{P}}_{2}\bm{B},{\bm{C}}\right\rrbracket. Matricizing $\underline{\bm{Y}}_{2}$ and plugging in $\bm{B}_{1}$ and $\bm{C}_{1}$ leads to:

[TABLE]

Since $JK_{2}\geq\lceil F\rceil$, then $\bm{C}\odot\bm{(}P_{2}\bm{B})$ has full column rank almost surely [63], and $\bm{A}_{2}=\bm{A}\bm{\Lambda}_{1}\bm{\Pi}$ can be identified from (17). Therefore \hat{\underline{\bm{Y}}}=\left\llbracket{\bm{A}}_{2},{\bm{B}}_{1},{\bm{C}}_{1}\right\rrbracket reconstructs signal $\underline{\bm{Y}}$.

The proof shares insights with that of Theorem 2 [19]. The basic difference lies in the fact that $\bm{P}_{1}^{(1)},~{}\bm{P}_{3}^{(2)}$ are full row rank selection matrices instead of full rank dense matrices.

Appendix B Proof of theorems 4,5

To begin, we use Lemma 1 and observe that $\bm{P}_{1}^{(d)}\bm{A},$ $\bm{P}_{2}^{(d)}{\bm{B}},\bm{P}_{3}^{(d)}{\bm{C}}$ are drawn from non-singular absolutely continuous distributions. Note that $\bm{P}_{1}^{(d)},\bm{P}_{2}^{(d)},\bm{P}_{3}^{(d)}$ have full row rank by construction. Then we use Theorem 1 to claim identifiability of the factors of each sub-tensor $\underline{\bm{Y}}_{d}$. Under the conditions of Theorems 4, 5, $\bm{P}_{1}^{(d)}\bm{A},\bm{P}_{2}^{(d)}{\bm{B}},\bm{P}_{3}^{(d)}{\bm{C}}$ can be identified, which corresponds to identifying all the rows of $\bm{A},{\bm{B}},{\bm{C}}$, up to column permutation and scaling. The caveat is that the rows of the factors are subject to column permutation and scaling mismatch, since they are obtained by the CPD of independent sub-sampled tensors. For example, let \underline{\bm{Y}}_{d}=\left\llbracket{\bm{A}}_{d},{\bm{B}}_{d},{\bm{C}}_{d}\right\rrbracket and \underline{\bm{Y}}_{d^{\prime}}=\left\llbracket{\bm{A}}_{d^{\prime}},{\bm{B}}_{d^{\prime}},{\bm{C}}_{d^{\prime}}\right\rrbracket. Then, from equations (11), it becomes clear that in order to obtain $\bm{A},\bm{B},\bm{C}$ from $\bm{A}_{d},\bm{B}_{d},\bm{C}_{d},~{}d=1,\dots,D$ and complete $\underline{\bm{X}}$, the permutation and scaling mismatch should be resolved, i.e., ${\bm{\Pi}}^{(d)}={\bm{\Pi}}^{({d^{\prime}})}$, ${\bm{\Lambda}}_{i}^{(d)}={\bm{\Lambda}}_{i}^{({d^{\prime}})}$ for every $d,~{}d^{\prime}$. To do so the following lemma is being used:

Lemma 2

Assume the entries of $\bm{C}\in\mathbb{C}^{K\times F}$ are jointly drawn from an absolutely continuous distribution over $\mathbb{C}^{KF}$. Then $\bm{C}(i,f)\neq\bm{C}(i^{\prime},f^{\prime})$ almost surely.

Proof:

The proof is similar to the proof of Corollary 1 in [63] and uses the fact that $\bm{C}(i,f)-\bm{C}(i^{\prime},f^{\prime})$ is a non-trivial analytic function of the entries of $\bm{C}$ and thus $\bm{C}(i,f)-\bm{C}(i^{\prime},f^{\prime})\neq 0$ almost surely. ∎

Finally, to resolve the mismatch, we utilize the rules in (7), (9). Particularly, when the original tensor is fiber sampled $\bm{C}_{d}=\bm{C}_{d^{\prime}}={\bm{C}},~{}\forall d,~{}d^{\prime}$ up to column permutation and scaling. Then, column permutation can be fixed to be the same for all $\underline{\bm{Y}}_{d}$s, since the entries of $\bm{C}$ are not equal almost surely. In order to reconcile for scaling mismatch, (7c), guarantees that there exist at least one row of $\bm{A}$ or $\bm{B}$ that is identified (up to permutation and scaling) from 2 different sub-sampled tensors $\underline{\bm{Y}}_{d}$. This is sufficient to resolve the CPD scaling mismatch between every $\underline{\bm{Y}}_{d}$-$\underline{\bm{Y}}_{d^{\prime}}$ couple, due to Lemma 1. The entry sampling mechanism, differs to the fiber sampling one, in the fact that $\bm{C}$ can only be partially identified from each sub-sampled version $\underline{\bm{Y}}_{d}$. Following same principles as before, permutation and scaling mismatch on the CPD of different $\underline{\bm{Y}}_{d}$s is resolved by (9c), (9d) along with Lemma 1.

Appendix C Proof of Theorems 6,7

The proof is similar to that of Theorem 3, 4, 5. The main difference lies in the fact that Theorem 2 is now employed, to establish identifiability on the CPD of each sub-sampled tensor and therefore recoverability of the original tensor. Furthermore, permutation and scaling alignment is performed using the rows of the latent factors which are common among the sub-tensors. This is accomplished, since factors with repeated entries are not allowed.