# Sparse Nonnegative Tucker Decomposition and Completion under Noisy
Observations

Xiongjun Zhang111School of Mathematics and Statistics and Hubei Key Laboratory
of Mathematical Sciences, Central China Normal University, Wuhan 430079, China
(e-mail: xjzhang@mail.ccnu.edu.cn). The research of this author was supported
in part by the National Natural Science Foundation of China under Grant No.
12171189 and the Knowledge Innovation Project of Wuhan under Grant No.
2022010801020279.  and Michael K. Ng222Department of Mathematics, The
University of Hong Kong, Pokfulam, Hong Kong (e-mail: mng@maths.hku.hk). The
research of this author was supported in part by the Hong Kong Research Grant
Council GRF 12300218, 12300519, 17201020, 17300021, C1013-21GF, C7004-21GF and
Joint NSFC-RGC N-HKU76921.

${\bf Abstract}$

Tensor decomposition is a powerful tool for extracting physically meaningful
latent factors from multi-dimensional nonnegative data, and has been an
increasing interest in a variety of fields such as image processing, machine
learning, and computer vision. In this paper, we propose a sparse nonnegative
Tucker decomposition and completion method for the recovery of underlying
nonnegative data under noisy observations. Here the underlying nonnegative
data tensor is decomposed into a core tensor and several factor matrices with
all entries being nonnegative and the factor matrices being sparse. The loss
function is derived by the maximum likelihood estimation of the noisy
observations, and the $\ell_{0}$ norm is employed to enhance the sparsity of
the factor matrices. We establish the error bound of the estimator of the
proposed model under generic noise scenarios, which is then specified to the
observations with additive Gaussian noise, additive Laplace noise, and Poisson
observations, respectively. Our theoretical results are better than those by
existing tensor-based or matrix-based methods. Moreover, the minimax lower
bounds are shown to be matched with the derived upper bounds up to logarithmic
factors. Numerical examples on both synthetic and real-world data sets
demonstrate the superiority of the proposed method for nonnegative tensor data
completion.

Key Words: Sparse nonnegative Tucker decomposition, maximum likelihood
estimation, noisy observations, error bound

Mathematics Subject Classification 2020: 15A69, 90C26

##  1 Introduction

Tensors, also called multi-dimensional data, are high-order generalizations of
vectors and matrices, and have a variety of applications including signal
processing, machine learning, computer vision, and so on [27, 9]. The data
required from real-world applications are usually represented as high-order
tensors, for example, color images are third-order tensors, color videos are
fourth-order tensors. In order to explore the intrinsic structure of tensor
data, tensor decomposition is a powerful and effective approach to represent
the data, and has been widely used in several domains [8, 44]. Besides, tensor
decomposition can identify the inherent structure of tensor data and extract
meaningful and explanatory information by a lower-dimensional set of latent
factors, which is due to the fact that many high-dimensional data reside in a
low-dimensional subspace [10].

The observation tensor data and latent factors are naturally nonnegative in
several real-world applications such as images, video volumes, and text data.
Nonnegative tensor decomposition can be employed to extract the intrinsic
structure of a tensor by the decomposition formulation [32]. For example, many
high-dimensional data, such as nonnegative hyperspectral images or video
images, are factorized to find meaningful latent nonnegative components [40].
Besides dimensionality reduction, nonnegative tensor decomposition can be
modeled and interpreted better by means of nonnegative and sparse components,
which can achieve a unique additive parts-based representation.

In particular, when the order of a tensor is second, nonnegative tensor
decomposition reduces to nonnegative matrix factorization (NMF), which has
been attracted much attention in the past decades, see [29, 15, 39, 17] and
references therein. The NMF, which requires the factors of the low-rank
decomposition to be componentwise nonnegative, is capable of giving physically
meaningful and more interpretable results, and has the ability of leaning the
local parts of objects [29]. There are plenty of real-world applications for
NMF [15]. For example, the NMF method can extract parts of faces, such as
eyes, noses, and lips, in a series of facial images [42], identify topics in a
set of documents [2], and extract materials and their proportions in
hyperspectral images unmixing [34]. Moreover, sparse NMF was introduced
further because it enhances the ability of NMF to learn a parts-based
representation and produces more easily interpretable factors. For example, in
facial feature extraction, sparsity leads to more localized features and
identifiable solutions, while fewer features are used to reconstruct each
input image [16]. Besides, Hoyer [23] proposed a sparse NMF method with sparse
constraints to improve the decomposition for the observations with additive
Gaussian noise, and designed a projected gradient descent algorithm to solve
the resulting model. In theory, Soni et al. [46] proposed a noisy matrix
completion method under a class of general noise models, and established the
error bound of the estimator of their proposed model, which can reduce to
sparse NMF and completion under nonnegative constraints. Moreover, Sambasivan
et al. [43] showed the error bound in [46] achieves minimax error rates up to
multiplicative constants and logarithmic factors. However, for multi-
dimensional data, the matrix based method may destroy the structure of a
tensor via unfolding a tensor into a matrix.

For nonnegative tensor decomposition, the key issue is the decomposition
formulation of a tensor. There are some popular decompositions for tensors
such as CANDECOMP/PARAFAC (CP) decomposition [19], Tucker decomposition [48],
tensor train decomposition [38], tensor decomposition via tensor-tensor
product [25]. Due to the nonnegativity of a tensor, the nonnegative CP
decomposition provides an interpretable, low tensor rank representation of the
data and has been used in a variety of applications related to sparse image
representation and image processing. For example, Zhang et al. [52] employed a
nonnegative CP decomposition for hyperspectral unmixing with the goal to
identify materials and material abundance in the data and to compress the
data. Moreover, Hong et al. [20] proposed a generalized CP low-rank tensor
decomposition that allows other loss functions besides squared error, and
presented a gradient descent type algorithm to solve the resulting model,
which can handle missing data by a similar approach in [1]. Besides, Jain et
al. [24] proposed a CP decomposition and completion method with one sparse
factor under a general class of noise models for third-order tensors, and
established the error bound of the estimator of their proposed model, which
was further derived to the error bound for the observations with additive
Gaussian noise. The method in [24] can reduce to nonnegative CP decomposition
with one sparse factor under nonnegative constraint of factor matrices.
However, the previous methods need to give the CP rank of a tensor in advance,
which is NP-hard in general [18].

Based on the algebra framework of tensor-tensor product in [25], Soltani et
al. [45] proposed a nonnegative tensor decomposition method with one sparse
factor for third-order tensors under additive Gaussian noise, which was
constructed by a tensor dictionary prior from training data for tomographic
image reconstruction. Furthermore, Newman et al. [37] proposed a sparse
nonnegative tensor patch-based dictionary approach based on tensor-tensor
product for image compression and deblurring, where the observations were
corrupted by additive Gaussian noise. Recently, Zhang [54] proposed a sparse
nonnegative tensor factorization and completion based on the algebra framework
of tensor-tensor product for third-order tensors, where one factor tensor was
sparse and the observations were corrupted by a general class of noise models.
Then the error bound of the estimator in [54] was established, and the minimax
lower bound was derived, which matched with the established upper bounds up to
a logarithmic factor of the sizes of the underlying tensor. However, for
$n$-dimensional $m$-th order tensors, there are $n^{m-2}$ $n$-by-$n$ matrix
singular value decompositions to be computed [25] and its computational cost
may be quite high in the tubular singular value decomposition approach. For
nonnegative TT decomposition, Lee et al. [30] proposed a hierarchical
alternating least squares algorithm for feature extraction and clustering.
Although the performance of nonnegative TT decomposition was better than that
of standard nonnegative Tucker decomposition, the TT decomposition typically
lacks interpretability and the efficient implementation depends on the tensor
already being in TT format [37].

Tucker decomposition is to decompose a given tensor into the product of a core
tensor with smaller dimensions and a series of factor matrices. And the best
low-rank approximation of Tucker decomposition of a tensor was discussed and
studied in [13, 12]. Nonnegative Tucker decomposition (NTD) provides an
additional multiway structured representation of tensor data with all entries
of the core tensor and factor matrices being nonnegative in Tucker
decomposition [40], which has a wide range of applications including image
denoising and hyperspectral image restoration [32, 3, 56]. Moreover, Mørup et
al. [36] proposed two multiplicative update algorithms for sparse NTD with the
observations including additive Gaussian noise and Poisson observations,
respectively, which yields a parts-based representation and is a more
interpretable decomposition. Liu et al. [33] presented a fast and flexible
algorithm for sparse NTD with a special core tensor based on columnwise
coordinate descent, where the observations were corrupted by additive Gaussian
noise and the factor matrices were nonnegative and sparse. Besides, Xu [51]
designed an alternating proximal gradient algorithm for sparse NTD and
completion with global convergence guarantee, where the observations were
corrupted by additive Gaussian noise. However, there is no theoretical result
about the error bounds of the models for sparse NTD with missing values in the
previous work. Also a general class of noise models is not discussed and
studied for sparse NTD and completion in the existing literature.

In this paper, we propose a sparse NTD and completion approach for the
observations with a general class of noise models. The loss function is
derived by the maximum likelihood estimation of observations, and the
nonnegativity for the entries of the core tensor and the factor matrices in
Tucker decomposition are imposed. The $\ell_{0}$ norm is used to characterise
the sparsity of the factor matrices. Moreover, the error bound of the
estimator of the proposed model is established, and then error bounds of the
estimators for the observations with additive Gaussian noise, additive Laplace
noise, and Poisson observations are derived, respectively. Besides, the
minimax lower bound of general class of noisy observations is established,
which matches to the upper bound up to a logarithmic factor. Then an
alternating direction method of multipliers (ADMM) is designed to solve the
resulting model. Numerical experiments on synthetic data and image data
demonstrate the effectiveness of the proposed model compared with the matrix
based method [46] and the sparse nonnegative tensor factorization and
completion via tensor-tensor product [54]. We summarize the existing
literature and our contribution for sparse nonnegative tensor
decomposition/factorization in Table 1.

Table 1: Comparisons of sparse NTD and completion with existing methods based on different decompositions. Methods | Rank | Order | Noise type | Error bounds
---|---|---|---|---
Matrix based method [46] |  matrix rank $r$ | second-order | general | $O(\frac{rn_{1}+\|A^{*}\|_{0}}{m}\log(n_{1}\vee n_{2}))$
Sparse NCPD [24] |  CP rank $r$ | third-order | general | $O(\frac{(n_{1}+n_{2})r+\|C\|_{0}}{m}\log(\max\\{n_{1},n_{2},n_{3},r\\}))$
NTTD [30] | tensor train rank |  $d$th-order ($d\geq 3$) | Gaussian | N/A
Sparse NTF [54] |  tubal rank $r$ | third-order | general | $O(\frac{rn_{1}n_{3}+\|\mathcal{B}\|_{0}}{m}\log(\max\\{n_{1},n_{2}\\}))$
Sparse NTD (this paper) |  Tucker rank $(r_{1},\ldots,r_{d})$ |  $d$th-order ($d\geq 3$) | general | $O(\frac{r_{1}\cdots r_{d}+\sum_{i=1}^{d}\|A_{i}\|_{0}}{m}\log(\max\\{n_{1},\ldots,n_{d}\\}))$

The remaining parts of this paper are organized as follows. In Section 2, some
preliminaries about tensor Tucker decomposition and related information
theoretical quality are provided. In Section 3, a sparse NTD and completion
model is proposed under a general class of noisy observations. An upper error
bound of the estimator of the proposed model is established in Section 4.
Moreover, the error bounds of the estimators are further derived for the
observations with special noise models. Section 5 is devoted to the minimax
lower bound of the estimator. In Section 6, an ADMM based algorithm is
developed to solve the resulting model. In Section 7, numerical experiments
are presented to demonstrate the effectiveness of the proposed model.
Discussions and list of future work are given in Section 8. All proofs of
theorems are provided in the appendices.

##  2 Preliminaries

Let $\mathbb{R}^{n}$ and $\mathbb{R}_{+}^{n_{1}\times n_{2}\times\cdots\times
n_{d}}$ denote the $n$-dimensional Euclidean space with real numbers and the
set of $n_{1}\times n_{2}\times\cdots\times n_{d}$ tensors with nonnegative
real entries, respectively. Scalars, vectors, matrices, and tensors are
represented by lowercase letters, lowercase boldface letters, uppercase
letters, and capital Euler script letters, respectively, e.g.,
$x,\mathbf{x},X$, $\mathcal{X}$. $\log$ refers to natural logarithm throughout
this paper. $\lfloor x\rfloor$ and $\lceil x\rceil$ are the integer part of
$x$ and smallest integer that is larger or equal to $x$, respectively. We
denote $n_{m}=\max\\{n_{1},n_{2},\ldots,n_{d}\\}$, $a\vee b=\max\\{a,b\\}$,
and $a\wedge b=\min\\{a,b\\}$, respectively. For any positive integer $n$, we
let $[n]=\\{1,2,\ldots,n\\}$. For an arbitrary set $\mathfrak{C}$,
$|\mathfrak{C}|$ denotes the number of entries in $\mathfrak{C}$.
$\delta_{\mathfrak{C}}(\cdot)$ denotes the indicator function over
$\mathfrak{C}$, i.e., $\delta_{\mathfrak{C}}(x)=0$ if $x\in\mathfrak{C}$, and
$\delta_{\mathfrak{C}}(x)=\infty$ if $x\neq\mathfrak{C}$.

Let $\mathcal{X}_{i_{1}i_{2}\cdots i_{d}}$ denote the
$(i_{1},i_{2},\ldots,i_{d})$th entry of a tensor
$\mathcal{X}\in\mathbb{R}^{n_{1}\times n_{2}\times\cdots\times n_{d}}$. The
number of ways or dimensions of a tensor is called the order [27]. The
mode-$i$ unfolding of a tensor $\mathcal{X}\in\mathbb{R}^{n_{1}\times
n_{2}\times\cdots\times n_{d}}$, denoted by $\mathcal{X}_{(i)}$, arranges the
mode-$n$ fibers to be the columns of the resulting matrix, where the
$(i_{1},i_{2},\ldots,i_{d})$th entry of $\mathcal{X}$ maps to matrix entry
$(i_{n},j)$ with

$j=1+\sum_{k=1\atop k\neq n}^{d}(i_{k}-1)J_{k}\ \textup{and}\
J_{k}=\prod_{l=1\atop l\neq n}^{k-1}n_{l}.$

Here a fiber is defined by fixing every index but one. The vectorization of a
tensor $\mathcal{X}$ is a vector, which is obtained by stacking all mode-$1$
fibers of $\mathcal{X}$ and denoted by vec($\mathcal{X}$).

The Tucker decomposition of a $d$th-order tensor
$\mathcal{X}\in\mathbb{R}^{n_{1}\times n_{2}\times\cdots\times n_{d}}$ is
given by [48]

$\mathcal{X}=\mathcal{C}\times_{1}A_{1}\times_{2}\cdots\times_{d}A_{d},$

where $\mathcal{C}\in\mathbb{R}^{r_{1}\times r_{2}\times\cdots\times r_{d}}$
is called the core tensor, $A_{i}\in\mathbb{R}^{n_{i}\times r_{i}}$ are the
factor matrices, $i\in[d]$, and $(r_{1},r_{2},\ldots,r_{d})$ is called the
Tucker rank of $\mathcal{X}$. Here the $n$-mode product of a tensor
$\mathcal{C}$ by a matrix $A\in\mathbb{R}^{n\times r_{i}}$, denoted by
$\mathcal{C}\times_{i}A$, is a $r_{1}\times\cdots\times r_{i-1}\times n\times
r_{i+1}\times\cdots\times r_{d}$-tensor whose entries are given by [12,
Definition 3]

$(\mathcal{C}\times_{i}A)_{j_{1}j_{2}\cdots j_{i-1}k_{i}j_{i+1}\cdots
j_{d}}=\sum_{j_{i}=1}^{r_{i}}\mathcal{C}_{j_{1}j_{2}\cdots
j_{i-1}j_{i}j_{i+1}\cdots j_{d}}A_{k_{i}j_{i}}.$

It is easy to verify that if
$\mathcal{X}=\mathcal{C}\times_{1}A_{1}\times_{2}\cdots\times_{d}A_{d}$, then

$\textup{vec}(\mathcal{X})=(\otimes_{i=d}^{1}A_{i})\textup{vec}(\mathcal{C}),$
(1)

where $\otimes_{i=d}^{1}A_{i}:=A_{d}\otimes A_{d-1}\otimes\cdots\otimes A_{1}$
and $A\otimes B$ denotes the Kronecker product of $A$ and $B$. Moreover, the
mode-$n$ unfolding of $\mathcal{X}$ is given by [27]

$\mathcal{X}_{(n)}=A_{n}\mathcal{C}_{(n)}(\otimes_{i=d,i\neq
n}^{1}A_{i})^{T}.$ (2)

For an arbitrary vector $\mathbf{x}\in\mathbb{R}^{n}$, $\|\mathbf{x}\|$ and
$\|\mathbf{x}\|_{1}$ denote the Euclidean norm and the $\ell_{1}$ norm of
$\mathbf{x}$, respectively, where $\|\mathbf{x}\|_{1}=\sum_{i=1}^{n}|x_{i}|$
and $x_{i}$ is the $i$th component of $\mathbf{x}$, $i\in[n]$. $\|X\|_{0}$
denotes the number of nonzero entries of a matrix $X$. The tensor Frobenius
norm of $\mathcal{X}\in\mathbb{R}^{n_{1}\times n_{2}\times\cdots\times n_{d}}$
is defined as
$\|\mathcal{X}\|_{F}=\sqrt{\langle\mathcal{X},\mathcal{X}\rangle}$, where the
inner product of two same-sized tensors is defined as
$\langle\mathcal{X},\mathcal{Y}\rangle=\sum_{i_{1},i_{2},\ldots,i_{d}}\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}\mathcal{Y}_{i_{1}i_{2}\cdots i_{d}}$. The tensor infinity norm of
$\mathcal{X}$, denoted by $\|\mathcal{X}\|_{\infty}$, is defined as
$\|\mathcal{X}\|_{\infty}=\max_{i_{1},i_{2},\ldots,i_{d}}|\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}|$.

Let $p_{x_{1}}(y)$ and $p_{x_{2}}(y)$ denote the probability density functions
or probability mass functions of a real scalar random variable $y$ with
corresponding parameters $x_{1}$ and $x_{2}$, respectively. The Kullback-
Leibler (KL) divergence of $p_{x_{1}}(y)$ from $p_{x_{2}}(y)$ is denoted and
defined as follows:

$K(p_{x_{1}}(y)||p_{x_{2}}(y))=\mathbb{E}_{p_{x_{1}}(y)}\left[\log\frac{p_{x_{1}}(y)}{p_{x_{2}}(y)}\right].$

The Hellinger affinity between $p_{x_{1}}(y)$ and $p_{x_{2}}(y)$ is denoted
and defined as

$H(p_{x_{1}}(y)||p_{x_{2}}(y))=\mathbb{E}_{p_{x_{1}}}\left[\sqrt{\frac{p_{x_{2}}(y)}{p_{x_{1}}(y)}}\right]=\mathbb{E}_{p_{x_{2}}}\left[\sqrt{\frac{p_{x_{1}}(y)}{p_{x_{2}}(y)}}\right].$

The joint probability distributions of random tensors, denoted by
$p_{\mathcal{X}_{1}}(\mathcal{Y}),p_{\mathcal{X}_{2}}(\mathcal{Y})$, are the
joint probability distributions of the vectorization of the tensors. Then the
KL divergence of $p_{\mathcal{X}_{1}}(\mathcal{Y})$ from
$p_{\mathcal{X}_{2}}(\mathcal{Y})$ is defined as

$K(p_{\mathcal{X}_{1}}(\mathcal{Y})||p_{\mathcal{X}_{2}}(\mathcal{Y})):=\sum_{i_{1},i_{2},\ldots,i_{d}}K(p_{(\mathcal{X}_{1})_{i_{1}i_{2}\cdots
i_{d}}}(\mathcal{Y}_{i_{1}i_{2}\cdots
i_{d}})||p_{(\mathcal{X}_{2})_{i_{1}i_{2}\cdots
i_{d}}}(\mathcal{Y}_{i_{1}i_{2}\cdots i_{d}})),$

and its Hellinger affinity is defined as

$H(p_{\mathcal{X}_{1}}(\mathcal{Y})||p_{\mathcal{X}_{2}}(\mathcal{Y})):=\prod_{i_{1},i_{2},\ldots,i_{d}}H(p_{(\mathcal{X}_{1})_{i_{1}i_{2}\cdots
i_{d}}}(\mathcal{Y}_{i_{1}i_{2}\cdots
i_{d}}),p_{(\mathcal{X}_{2})_{i_{1}i_{2}\cdots
i_{d}}}(\mathcal{Y}_{i_{1}i_{2}\cdots i_{d}})).$

Let $\mathfrak{R}$ be a Euclidean space endowed with the Euclidean norm
$\|\cdot\|$. For an arbitrary closed proper function
$g:\mathfrak{R}\rightarrow(-\infty,+\infty]$, the proximal mapping of $g$ at
$y$ is the operator given by

$\mbox{Prox}_{g}(y)=\arg\min_{x\in\mathfrak{R}}\left\\{g(x)+\frac{1}{2}\|x-y\|^{2}\right\\}.$

##  3 Sparse NTD and Completion

Let $\mathcal{X}^{*}\in\mathbb{R}_{+}^{n_{1}\times n_{2}\times\cdots\times
n_{d}}$ be an unknown nonnegative tensor with Tucker rank
$(r_{1},r_{2},\ldots,r_{d})$. Assume that $\mathcal{X}^{*}$ admits a following
sparse nonnegative Tucker decomposition:

$\mathcal{X}^{*}=\mathcal{C}^{*}\times_{1}A_{1}^{*}\times_{2}A_{2}^{*}\cdots\times_{d}A_{d}^{*},$
(3)

where $\mathcal{C}^{*}\in\mathbb{R}_{+}^{r_{1}\times r_{2}\times\cdots\times
r_{d}}$ is the core tensor and $A_{i}^{*}\in\mathbb{R}_{+}^{n_{i}\times
r_{i}}$ are sparse factor matrices, $i\in[d]$. Suppose that $r_{i}\leq
n_{i},i\in[d]$. We also assume that each entries of
$\mathcal{X}^{*},\mathcal{C}^{*},A_{i}^{*}$ are bounded, i.e.,

$\begin{split}&0\leq\mathcal{X}_{i_{1}i_{2}\cdots i_{d}}^{*}\leq\frac{c}{2},\
\
(i_{1},i_{2},\ldots,i_{d})\in[n_{1}]\times[n_{2}]\times\cdots\times[n_{d}],\\\
&0\leq\mathcal{C}_{i_{1}i_{2}\cdots i_{d}}^{*}\leq 1,\ \
(i_{1},i_{2},\ldots,i_{d})\in[r_{1}]\times[r_{2}]\times\cdots\times[r_{d}],\\\
&0\leq(A_{i}^{*})_{lm}\leq a_{i},\ \
(l,m)\in[n_{i}]\times[r_{i}],i\in[d],\end{split}$ (4)

where $c,a_{i}>0$ are given constants. Here we use $\frac{c}{2}$ for brevity
in the subsequence analysis.

However, only a noisy and incompleted version of the underlying tensor
$\mathcal{X}^{*}$ is available in practice. Let
$\Omega\subseteq[n_{1}]\times[n_{2}]\times\cdots\times[n_{d}]$ be a subset at
which the entries of the observations $\mathcal{Y}$ are collected. Denote
$\mathcal{Y}_{\Omega}\in\mathbb{R}^{m}$ to be a vector such that the entries
of $\mathcal{Y}$ in the index $\Omega$ are vectorized by lexicographic order,
where $m$ is the number of the observed entries, i.e., $m=|\Omega|$. Assume
that $n_{1},n_{2},\ldots,n_{d}\geq 2$ and $d\geq 3$ throughout this paper. Let
the location $\Omega$ set be generated according to an independent Bernoulli
model with probability $p=\frac{m}{n_{1}n_{2}\cdots n_{d}}$ (denoted by
Bern($p$)), i.e., each index $(i_{1},i_{2},\ldots,i_{d})$ belongs to $\Omega$
with probability $p$, which is denoted as $\Omega\sim\text{Bern}(p)$. Suppose
that the entries of observations are conditionally independent, which implies
that the overall likelihood is just the product of the likelihoods. Therefore,
the joint probability density function or probability mass function of
observations $\mathcal{Y}_{\Omega}$ is given by

$p_{\mathcal{X}_{\Omega}^{*}}(\mathcal{Y}_{\Omega}):=\prod_{(i_{1},i_{2},\ldots,i_{d})\in\Omega}p_{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*}}(\mathcal{Y}_{i_{1}i_{2}\cdots i_{d}}).$ (5)

By taking the negative logarithm of the probability distribution, we propose
the following sparse NTD and completion model with nonnegative and bounded
constraints:

$\min_{\mathcal{X}\in\Upsilon}\left\\{-\log
p_{\mathcal{X}_{\Omega}}(\mathcal{Y}_{\Omega})+\sum_{i=1}^{d}\lambda_{i}\|A_{i}\|_{0}\right\\},$
(6)

where $\lambda_{i}>0$ are the regularization parameters, $i\in[d]$,
$\|A_{i}\|_{0}$ is employed to characterize the sparsity of the factor matrix
$A_{i}$ in Tucker decomposition, and $\Upsilon$ is defined as

$\begin{split}\Upsilon:=\Big{\\{}&\mathcal{X}=\mathcal{C}\times_{1}A_{1}\cdots\times_{d}A_{d}:\
\mathcal{C}\in\mathfrak{C},\ A_{i}\in\mathfrak{B}_{i},i\in[d],\\\
&0\leq\mathcal{X}_{i_{1}i_{2}\cdots i_{d}}\leq c,\
(i_{1},i_{2},\ldots,i_{d})\in[n_{1}]\times[n_{2}]\times\cdots\times[n_{d}]\Big{\\}}.\end{split}$
(7)

Here $\Pi$ is a countable set of candidate estimates, and $\mathfrak{C}$ and
$\mathfrak{B}_{i}$ are the nonnegative and bounded sets constructed as
follows: First, let

$\tau:=2^{\lceil\log_{2}(n_{m})^{\beta}\rceil}$ (8)

for a specified $\beta\geq 1,$ we construct $\mathfrak{C}$ to be the set of
all nonnegative tensors $\mathcal{C}\in\mathbb{R}_{+}^{r_{1}\times
r_{2}\times\cdots\times r_{d}}$ whose entries are discretized to one of $\tau$
uniformly sized bins in the range $[0,1]$, and $\mathfrak{B}_{i}$ to be the
set of all nonnegative matrices $A_{i}\in\mathbb{R}_{+}^{n_{i}\times r_{i}}$
whose entries either take the value $0$, or are discretized to one of $\tau$
uniformly sized bins in the range $[0,a_{i}]$, $i\in[d]$.

###### Remark 3.1

The model in (6) can address the observations with general noise
distributions, where one just needs to know the probability mass function or
probability density function of the observations $\mathcal{Y}_{\Omega}$. In
particular, the noise models including additive Gaussian noise, additive
Laplace noise, and Poisson observations will be discussed in detail in the
next section.

###### Remark 3.2

For the observations with additive Gaussian noise, Liu et al. [33] proposed a
sparse NTD model with special core tensor $\mathcal{C}$, which is effective
for image compression. Moreover, Xu [51] proposed a gradient descent type
algorithm for sparse NTD and completion with additive Gaussian noise, which
used the relaxation of $\ell_{0}$ norm for each factor matrix. Besides, Mørup
et al. [36] proposed two multiplicative update algorithms for sparse NTD,
where additive Gaussian noise and Poisson observations were considered with
full observations. However, there is no theoretical analysis of error bounds
about these models in [33, 36, 51]. Only efficient algorithms were proposed
and studied to solve their resulting sparse NTD models. Moreover, a general
loss function framework has been proposed in (6) for sparse NTD and
completion, while the special loss function with commonly used noise types is
considered in the existing literature, such as additive Gaussian noise [33,
51], Poisson observations [36].

###### Remark 3.3

Recently, Zhang et al. [54] proposed a sparse nonnegative tensor factorization
and completion model with general noise distributions based on the algebra
framework of tensor-tensor product for third-order tensor, where the
underlying tensor was decomposed into the tensor-tensor product of one sparse
nonnegative tensor and one nonnegative tensor. The difference between model
(6) and [54] is the factorization of the underlying tensor, where the tensor-
tensor product with block circulant structure [25] is used in [54] and the
Tucker decomposition of the underlying tensor is utilized in model (6). And
the nonnegative Tucker decomposition has more widely applications and
explanatory information than the nonnegative tensor factorization based on
tensor-tensor product, where the last one is mainly applied in X-ray CT
imaging, image compression and deblurring based on tensor patch dictionary
learning [37]. Besides, the model in [54] is only effective for third-order
tensors, while the model in (6) can be utilized to any order tensors and can
extract more physically meaningful latent components by the sparse and
nonnegative factors [56].

###### Remark 3.4

Jain et al. [24] proposed a noisy tensor completion method based on CP
decomposition, where one factor is sparse. The model (6) can reduce to sparse
nonnegative CP decomposition and completion when the core tensor is diagonal.
However, the core tensor of Tucker decomposition is not diagonal in general in
many real-world applications [56]. Moreover, Hong et al. [20] proposed a
general loss function framework for CP decomposition under various scenarios
and designed a gradient descent algorithm to solve the resulting model, while
there are no nonnegative and sparse constraints for the factor matrices, and
there is no theoretical analysis of error bounds about the model in [20].

##  4 Error Bounds

In this section, an upper error bound of the estimator of the sparse NTD and
completion model in (6) is established under a general class of noise
distributions, and then the upper error bounds of the estimators for the
observations with special noise models are derived, including additive
Gaussian noise, additive Laplace noise, and Poisson observations.

Let $\mathcal{X}^{\lambda}$ be a minimizer of (6). Next we state the main
theorem about the error bound of the estimator in (6).

###### Theorem 4.1

Let the sampling set $\Omega$ be drawn from the independent Bernoulli model
with probability $p=\frac{m}{n_{1}n_{2}\cdots n_{d}}$, i.e.,
$\Omega\sim\textup{Bern}(p)$, and the joint probability density/mass function
of $\mathcal{Y}_{\Omega}$ be defined as (5). For any

$\lambda_{i}\geq 4(\beta+2)\left(1+\frac{2\gamma}{3}\right)\log(n_{m}),\ \
i\in[d],$

where $\gamma$ is a constant satisfying

$\gamma\geq\max_{\mathcal{X}\in\Upsilon}\max_{i_{1},i_{2},\ldots,i_{d}}K\left(p_{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*}}(\mathcal{Y}_{i_{1}i_{2}\cdots
i_{d}})||p_{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}}(\mathcal{Y}_{i_{1}i_{2}\cdots i_{d}})\right),$ (9)

then

$\begin{split}&\ \frac{\mathbb{E}_{\Omega,\mathcal{Y}_{\Omega}}\left[-2\log
H(p_{{\mathcal{X}}^{\lambda}}(\mathcal{Y}),p_{\mathcal{X}^{*}}(\mathcal{Y}))\right]}{n_{1}n_{2}\cdots
n_{d}}\leq\frac{8\gamma\log(m)}{m}+3\cdot\min_{\mathcal{X}\in\Upsilon}\Biggl{\\{}\frac{K(p_{\mathcal{X}^{*}}(\mathcal{Y})||p_{\mathcal{X}}(\mathcal{Y}))}{n_{1}n_{2}\cdots
n_{d}}\\\ &\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\left(\max_{i}\\{\lambda_{i}\\}+\frac{8\gamma(\beta+2)\log(n_{m})}{3}\right)\frac{r_{1}r_{2}\cdots
r_{d}+\sum_{i=1}^{d}\|A_{i}\|_{0}}{m}\Biggl{\\}},\end{split}$

where the expectation is taken with respect to the joint distribution of
$\Omega$ and $\mathcal{Y}_{\Omega}$.

Theorem 4.1 is the upper error bound of the estimator in (6) under a general
class of noise distributions, which can be specified to some special noise
distributions, such as additive Gaussian noise, additive Laplace noise, and
Poisson observations. The key point for showing the upper error bound with
special noise distributions is to specify the logarithmic Hellinger affinity
$\log
H(p_{{\mathcal{X}}^{\lambda}}(\mathcal{Y}),p_{\mathcal{X}^{*}}(\mathcal{Y}))$
and the minimum KL divergence
$\min_{\mathcal{X}\in\Upsilon}K(p_{\mathcal{X}^{*}}(\mathcal{Y})||p_{\mathcal{X}}(\mathcal{Y}))$.

From Theorem 4.1, we can see that the negative logarithmic Hellinger affinity
with respect to the underlying tensor and the recovered tensor is on the order
of $O(\frac{r_{1}r_{2}\cdots
r_{d}+\sum_{i=1}^{d}\|A_{i}\|_{0}}{m}\log(n_{m}))$ if the KL divergence
$\min_{\mathcal{X}\in\Upsilon}K(p_{\mathcal{X}^{*}}(\mathcal{Y})||p_{\mathcal{X}}(\mathcal{Y}))$
is not too large and $\lambda_{i}$ are fixed, where $A_{i}$ are the factor
matrices of $\mathcal{X}\in\Upsilon$ in the Tucker decomposition. Therefore,
the Tucker rank is lower, the error bound of the estimator in (6) is smaller.

###### Remark 4.1

For the problem of sparse NTD with partial observations, the upper error bound
is related to the degree of freedoms of the tensor, i.e., $r_{1}r_{2}\cdots
r_{d}+\sum_{i=1}^{d}\|A_{i}\|_{0}$. Moreover, it follows from Theorem 4.1 that
the upper error bound decreases as the number of observed samples increases.

In the following subsections, we specify the noise models and establish the
detailed upper error bounds based on Theorem 4.1. In the discretization levels
(8), we choose

$\beta=1+\frac{\log\left(\frac{(2^{d+1}-1)\sqrt{d}r_{1}r_{2}\cdots
r_{d}a_{1}a_{2}\cdots a_{d}}{c\sqrt{n_{m}}}+1\right)}{\log(n_{m})},$ (10)

which implies that $\beta\geq 1$. For the regularization parameters
$\lambda_{i}$, we consider the specific choice

$\lambda_{i}=4(\beta+2)\left(1+\frac{2\gamma}{3}\right)\log(n_{m}),\ \
i\in[d],$ (11)

where $\gamma$ satisfies (9) and is related to the maximum KL divergence
between $p_{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*}}(\mathcal{Y}_{i_{1}i_{2}\cdots i_{d}})$ and
$p_{\mathcal{X}_{i_{1}i_{2}\cdots i_{d}}}(\mathcal{Y}_{i_{1}i_{2}\cdots
i_{d}})$ for any $\mathcal{X}\in\Upsilon$.

###  4.1 Additive Gaussian Noise

In this subsection, we establish the upper error bound of the estimator in (6)
for the observations corrupted by additive Gaussian noise, which has been
widely used in a variety of applications. In this case, the observation
$\mathcal{Y}_{\Omega}$ is generated by

$\mathcal{Y}_{i_{1}i_{2}\cdots i_{d}}=\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*}+\mathcal{N}_{i_{1}i_{2}\cdots i_{d}},\ \ \
(i_{1},i_{2},\ldots,i_{d})\in\Omega,$

where $\mathcal{N}_{i_{1}i_{2}\cdots i_{d}}$ are independent zero-mean
Gaussian noise with standard deviation $\sigma>0$. The distribution of
$\mathcal{Y}_{\Omega}$ obeys a multivariate Gaussian density with dimension
$|\Omega|$ whose mean and covariance matrix are $\mathcal{X}^{*}_{\Omega}$ and
$\sigma^{2}I_{\Omega}$, respectively, where $I_{\Omega}$ denotes the identity
matrix with dimension $|\Omega|\times|\Omega|$. Then the joint probability
density function of $\mathcal{Y}_{\Omega}$ is given by

$p_{\mathcal{X}_{\Omega}^{*}}(\mathcal{Y}_{\Omega})=\frac{1}{(2\pi\sigma^{2})^{|\Omega|/2}}\exp\left(-\frac{\|\mathcal{Y}_{\Omega}-\mathcal{X}_{\Omega}^{*}\|^{2}}{2\sigma^{2}}\right).$
(12)

Now we specify the upper error bound in Theorem 4.1 for the observations with
additive Gaussian noise in the following theorem, where the joint probability
density function of the observations is given by (12).

###### Theorem 4.2

Let $\beta$ and $\lambda_{i}$ be defined as (10) and (11), $i\in[d]$, where
$\gamma=\frac{c^{2}}{2\sigma^{2}}$ in (11). Assume that the sampling set
$\Omega\sim\textup{Bern}(p)$ with $p=\frac{m}{n_{1}n_{2}\cdots n_{d}}$ and the
joint probability density function of $\mathcal{Y}_{\Omega}$ is given by (12).
Then the estimator in (6) satisfies

$\begin{split}\frac{\mathbb{E}_{\Omega,\mathcal{Y}_{\Omega}}\left[\|\mathcal{X}^{\lambda}-\mathcal{X}^{*}\|_{F}^{2}\right]}{n_{1}n_{2}\cdots
n_{d}}\leq&\ \frac{22c^{2}\log(m)}{m}\\\ &\
+16(\beta+2)(2c^{2}+3\sigma^{2})\left(\frac{r_{1}r_{2}\cdots
r_{d}+\sum_{i=1}^{d}\|A_{i}^{*}\|_{0}}{m}\right)\log(n_{m}).\end{split}$

By Theorem 4.2, we know that the upper error bound of the estimator of (6) for
the observations with additive Gaussian noise is on the order of

$O\left(\frac{r_{1}r_{2}\cdots
r_{d}+\sum_{i=1}^{d}\|A_{i}^{*}\|_{0}}{m}\log(n_{m})\right).$

Consequently, if the factor matrices are sparser, the error upper bound of the
estimator is lower. Moreover, when $m=n_{1}\cdots n_{d}$ and
$\|A_{i}\|_{0}=r_{i}n_{i},i\in[d]$, which implies that $A_{i}$ are dense, the
error bound in Theorem 4.2 is just the error bound of the estimator of NTD
with additive Gaussian noise, which has been widely used and studied, see [10,
55, 36] and references therein.

###### Remark 4.2

For sparse NTD with additive Gaussian noise, Liu et al. [33] proposed a novel
model with a special core tensor in the Tucker decomposition, where the
$\ell_{0}$ norm is replaced by the $\ell_{1}$ norm for the factor matrices. In
this case, Theorem 4.2 establishes the error bound of the model in [33] by
employing the $\ell_{0}$ norm to characterize the sparsity of the factor
matrices. Moreover, Xu et al. [51] proposed a sparse NTD and completion
approach for the observations with additive Gaussian noise, which replaced the
$\ell_{0}$ norm by the $\ell_{1}$ norm in (6). This implies that Theorem 4.2
established the error bound of the estimator of the model in [51] when the
$\ell_{0}$ norm is utilized to measure the sparsity of the factor matrices.

###### Remark 4.3

If we ignore the internal structure of a tensor and unfold the tensor into a
matrix, we can compare the error bound in Theorem 4.2 with that in [46,
Corollary 3], where the matrix based method in [46] can reduce to sparse
nonnegative matrix factorization and completion with nonnegative and bounded
constraints. For a given tensor $\mathcal{X}^{*}$ in the form of (3), without
loss of generality, we unfold it along the first mode, i.e.,
$\mathcal{X}_{(1)}^{*}=A_{1}^{*}(\mathcal{C}_{(1)}^{*}(\otimes_{i=d}^{2}A_{i}^{*})^{T}).$
In this case, one is capable of applying the matrix based method in [46] to
$(\mathcal{X}_{(1)}^{*})^{T}\in\mathbb{R}^{(n_{2}\cdots n_{d})\times
n_{1}}_{+}$, where the sparse factor is $A_{1}^{*}$. Then the resulting error
bound of the estimator of the matrix based method [46, Corollary 3] is on the
order of

$O\left(\frac{r_{1}n_{2}\cdots
n_{d}+\|A_{1}^{*}\|_{0}}{m}\log\left(n_{1}\vee(n_{2}\cdots
n_{d})\right)\right).$

Now we compare it with the case that only nonnegative tensor decomposition is
considered. In this case, the error bound of the matrix based method in [46]
for the tensor with mode-$1$ unfolding is the order of
$O(\frac{n_{1}r_{1}+r_{1}n_{2}n_{3}\cdots n_{d}}{m}\log(n_{1}\vee(n_{2}\cdots
n_{d}))$. And the error bound in Theorem 4.2 is the order of
$O(\frac{r_{1}r_{2}\cdots r_{d}+\sum_{i=1}^{d}r_{i}n_{i}}{m}\log(n_{m}))$.
Therefore, the error bound in Theorem 4.2 is smaller than that in [46] if
$r_{i}$ is much smaller than $n_{i},i\in[d]$ and $d$ is not too large. In
particular, if $r_{1}=\cdots=r_{d}=r$ and $n_{1}=\cdots=n_{d}=n$, the error
bounds in Theorem 4.2 and in [46] are $O(\frac{r^{d}+drn}{m}\log(n))$ and
$O(\frac{rn+rn^{d-1}}{m}\log(n(d-1))$, respectively, where the error bound of
the matrix based method is larger than that in Theorem 4.2 if $r,d$ are much
smaller than $n$.

###### Remark 4.4

We compare the error bound in Theorem 4.2 with that in [54, Proposition 4.1],
which is based on the algebra framework of tensor-tensor product and is only
effective for third-order tensors. In general, the two error bound results are
not comparable directly since the sparse factors are different for the two
factorizations. However, when the factors matrices in Tucker decomposition and
the factor tensor in the algebra framework of tensor-tensor product are dense,
the two error bounds can be comparable. For a third-order tensor with size
$n_{1}\times n_{2}\times n_{3}$, the error bound of the estimator in Theorem
4.2 is $O(\frac{r_{1}r_{2}r_{3}+\sum_{i=1}^{3}r_{i}n_{i}}{m}\log(n_{m})),$
while it is $O(\frac{rn_{1}n_{3}+rn_{2}n_{3}}{m}\log(n_{1}\vee n_{2}))$ in
[54], where $r$ is the tubal rank of the underlying tensor. In this case, the
error bound in Theorem 4.2 is smaller than that in [54] if $r_{i}$ is much
smaller than $n_{i}$ and $r$ is close to $r_{i}$. In particular, if
$n_{1}=n_{2}=n_{3}=n$ and $r_{1}=r_{2}=r_{3}=r$, the error bound in Theorem
4.2 is $O(\frac{r^{3}+3rn}{m}\log(n)),$ which is smaller than that (i.e.,
$O(\frac{2rn^{2}}{m}\log(n))$) in [54] if $r\leq 2n-3$ and $n\geq 3$. These
conditions can be satisfied easily in real-world applications since $n$ is
generally large.

###### Remark 4.5

Jain et al. [24] proposed a noisy tensor completion model based on CP
decomposition with a special sparse factor for an $n_{1}\times n_{2}\times
n_{3}$ tensor, where the third factor matrix is sparse in CP decomposition. In
particular, when the constraints are nonnegative, the model in [24] reduces to
sparse nonnegative CP decomposition and completion. In this case, the error
bound of the estimator of their model for the observations with additive
Gaussian noise was

$\frac{\mathbb{E}_{\Omega,\mathcal{Y}_{\Omega}}[\|\mathcal{X}_{cp}-\mathcal{X}^{*}\|_{F}^{2}]}{n_{1}n_{2}n_{3}}=O\left(\frac{(n_{1}+n_{2})r+\|C^{*}\|_{0}}{m}\log(\max\\{n_{1},n_{2},n_{3}\\})\right),$
(13)

where $\mathcal{X}_{cp}$ is the estimator by [24], $r$ is the CP rank of
$\mathcal{X}^{*}$, and $C^{*}$ is the third factor matrix in the CP
decomposition of $\mathcal{X}^{*}$. However, (13) is hard to compare with the
error bound in Theorem 4.2 directly since the CP rank of a tensor is not
comparable to the Tucker rank of a tensor in general. Also the factor matrices
are typically different between the CP and Tucker decompositions of a tensor.
In particular, if the core tensor in Tucker decomposition is diagonal and
$r_{1}=r_{2}=r_{3}=r$, the error bound in Theorem 4.2 is comparable with (13).
In fact, the error bound in Theorem 4.2 is smaller than (13) if $r$ is much
smaller than $n_{i}$.

###  4.2 Additive Laplace Noise

In this subsection, we establish the error bound of the estimator of model (6)
for the observations with additive Laplace noise. In this case, the
observation $\mathcal{Y}_{\Omega}$ is modeled by

$\mathcal{Y}_{i_{1}\cdots i_{d}}=\mathcal{X}_{i_{1}\cdots
i_{d}}^{*}+\mathcal{N}_{i_{1}\cdots i_{d}},\ \ \
(i_{1},\ldots,i_{d})\in\Omega,$ (14)

where $\mathcal{N}_{i_{1}\cdots i_{d}}$ are independent Laplace distribution
with parameters $(0,\tau)$, $\tau>0$, which is denoted by
$\textup{Laplace}(0,\tau)$. The joint probability density function of
$\mathcal{Y}_{\Omega}$ is given by

$p_{\mathcal{X}_{\Omega}^{*}}(\mathcal{Y}_{\Omega})=\left(\frac{\tau}{2}\right)^{|\Omega|}\exp\left(-\frac{\|\mathcal{Y}_{\Omega}-\mathcal{X}_{\Omega}^{*}\|_{1}}{\tau}\right).$
(15)

Next we establish an upper error bound of the estimator of model (6) for the
observations satisfying the noisy model (14).

###### Theorem 4.3

Let $\beta$ and $\lambda_{i}$ be defined as (10) and (11), $i\in[d]$, where
$\gamma=\frac{c^{2}}{2\tau^{2}}$. Assume that the sampling set
$\Omega\sim\textup{Bern}(p)$ with $p=\frac{m}{n_{1}n_{2}\cdots n_{d}}$ and the
joint probability density function of $\mathcal{Y}_{\Omega}$ is given by (15).
Then the estimator of (6) satisfies

$\begin{split}\frac{\mathbb{E}_{\Omega,\mathcal{Y}_{\Omega}}\left[\|\mathcal{X}^{\lambda}-\mathcal{X}^{*}\|_{F}^{2}\right]}{n_{1}n_{2}\cdots
n_{d}}\leq&\ \frac{11c^{2}(2\tau+c)^{2}\log(m)}{2\tau^{2}m}\\\ &\
+12\left(1+\frac{2c^{2}}{3\tau^{2}}\right)(2\tau+c)^{2}(\beta+2)\log(n_{m})\frac{r_{1}r_{2}\cdots
r_{d}+\sum_{i=1}^{d}\|A_{i}^{*}\|_{0}}{m}.\end{split}$

From Theorem 4.3, we know that the upper error bound of the estimator of model
(6) is on the order of $O(\frac{r_{1}r_{2}\cdots
r_{d}+\sum_{i=1}^{d}\|A_{i}^{*}\|_{0}}{m}\log(n_{m}))$. Similar to the case of
additive Gaussian noise, the error bound of the estimator of model (6)
decreases as the number of observed samples increases. Moreover, if the factor
matrices of Tucker decomposition are sparser, the error bound obtained by (6)
is lower.

###### Remark 4.6

For the observations with additive Laplace noise, we compared the error bound
in Theorem 4.3 with that of the matrix based method in [46, Corollary 5]. For
a $d$th-order tensor
$\mathcal{X}^{*}\in\mathbb{R}_{+}^{n_{1}\times\cdots\times n_{d}}$ with Tucker
decomposition in (3), we unfold it into a matrix
$\mathcal{X}^{*}_{(1)}=A_{1}^{*}(\mathcal{C}_{(1)}^{*}(\otimes_{i=d}^{2}A_{i}^{*})^{T})\in\mathbb{R}_{+}^{n_{1}\times(n_{2}\cdots
n_{d})}$ along the first mode, where $A_{1}^{*}$ is sparse. The matrix based
method in [46] is then applied to
$(\mathcal{X}^{*}_{(1)})^{T}\in\mathbb{R}_{+}^{(n_{2}\cdots n_{d})\times
n_{1}}$. And the error bound of the estimator in [46, Corollary 5] is on the
order of $O(\frac{r_{1}n_{2}n_{3}\cdots
n_{d}+\|A_{1}^{*}\|_{0}}{m}\log(n_{1}\vee(n_{2}\cdots n_{d})))$. Therefore, if
$r_{i}$ is much smaller than $n_{i},i\in[d]$, and $d$ is not too large, the
error bound of the estimator by model (6) is smaller than that of the matrix
based method in [46]. In particular, if $r_{1}=\cdots=r_{d}=r$ and
$n_{1}=\cdots=n_{d}=n$, the error bound of the estimator in Theorem 4.3 is the
order of $O(\frac{r^{d}+\sum_{i=1}^{d}\|A_{i}^{*}\|_{0}}{m}\log(n))$, while
the error bound of the estimator of the matrix based method in [46] is
$O(\frac{rn^{d-1}+\|A_{1}^{*}\|_{0}}{m}(d-1)\log(n))$. This implies that the
error bound of the estimator in Theorem 4.3 is smaller than that of the matrix
based method in [46] as long as $r$ is much smaller than $n$ and
$\|A_{i}^{*}\|_{0}$ is close to $\|A_{1}^{*}\|_{0},i=2,\ldots,d$.

###### Remark 4.7

Similar to the case for the observations with additive Gaussian noise, we
compare the error bound in Theorem 4.3 with that of the sparse nonnegative
tensor factorization method under the tensor-tensor product framework in [54],
which is only effective for third-order tensors. Since the spare factors are
different between the sparse NTD and sparse nonnegative tensor factorization
with tensor-tensor product in [54], it is difficult to compare the error
bounds of the two methods directly. However, if the factor tensor is not
sparse, the error bound of the estimator in [54, Proposition 4.2] is on the
order of $O(\frac{rn_{1}n_{3}+rn_{2}n_{3}}{m}\log(n_{1}\vee n_{2}))$, where
$r$ is the tubal rank of the underlying tensor. The error upper bound of the
estimator by model (6) is on the order of
$O(\frac{r_{1}r_{2}r_{3}+\sum_{i=1}^{3}r_{i}n_{i}}{m}\log(n_{m}))$, which is
smaller than that of [54] if $r$ is close to $r_{i}$ and $r_{i}$ is much
smaller than $n_{i}$, $i=1,2,3$.

###  4.3 Poisson Observations

In this subsection, we establish an upper error bound of the estimator for
Poisson observations, where the observations are modeled as

$\mathcal{Y}_{i_{1}\cdots i_{d}}=\textup{Poisson}(\mathcal{X}_{i_{1}\cdots
i_{d}}^{*}),\ \ \ (i_{1},\ldots,i_{d})\in\Omega.$

Here $y=\textup{Poisson}(x)$ denotes that $y$ obeys a Poisson distribution
with parameter $x>0$. The joint probability mass function of
$\mathcal{Y}_{\Omega}$ is given by

$p_{\mathcal{X}_{\Omega}^{*}}(\mathcal{Y}_{\Omega})=\prod_{(i_{1},i_{2},\ldots,i_{d})\in\Omega}\frac{(\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*})^{\mathcal{Y}_{i_{1}i_{2}\cdots
i_{d}}}\exp({-\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*}})}{(\mathcal{Y}_{i_{1}i_{2}\cdots i_{d}})!}.$ (16)

Now we establish an upper error bound of the estimator of model (6), where the
joint probability mass function of $\mathcal{Y}_{\Omega}$ satisfies (16).

###### Theorem 4.4

Suppose that $\min_{i_{1},i_{2},\ldots,i_{d}}\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*}\geq\varrho$ and each entry of the tensor in the candidate set
$\Upsilon$ is also not smaller than $\varrho$, where $\varrho>0$ is a given
constant. Let $\beta$ and $\lambda_{i}$ be defined as (10) and (11),
$i\in[d]$, where $\gamma=\frac{c^{2}}{\varrho}$ in (11). Assume that the
sampling set $\Omega\sim\textup{Bern}(p)$ with $p=\frac{m}{n_{1}n_{2}\cdots
n_{d}}$ and the joint probability mass function of $\mathcal{Y}_{\Omega}$ is
given by (16). Then the estimator in (6) satisfies

$\begin{split}\frac{\mathbb{E}_{\Omega,\mathcal{Y}_{\Omega}}\left[\|\mathcal{X}^{\lambda}-\mathcal{X}^{*}\|_{F}^{2}\right]}{n_{1}n_{2}\cdots
n_{d}}\leq&\ \frac{44c^{3}\log(m)}{\varrho
m}+48c\left(1+\frac{4c^{2}}{3\varrho}\right)(\beta+2)\frac{r_{1}r_{2}\cdots
r_{d}+\sum_{i=1}^{d}\|A_{i}^{*}\|_{0}}{m}\log(n_{m}).\end{split}$

By Theorem 4.4, one can obtain that the upper error bound of the estimator of
model (6) with Poisson observations is on the order of
$O(\frac{r_{1}r_{2}\cdots
r_{d}+\sum_{i=1}^{d}\|A_{i}^{*}\|_{0}}{m}\log(n_{m}))$. Moreover, the error
bound will decrease as the number of observations increases. The lower bound
of $\mathcal{X}^{*}$ will also influence the error bound of the estimator, and
it will be difficult to recover $\mathcal{X}^{*}$ as $\varrho$ is close to
zero.

###### Remark 4.8

In order to compare the error bound in Theorem 4.4 with that of the matrix
based method in [46, Corollary 6], we need to unfold the underlying
nonnegative tensor $\mathcal{X}^{*}$ $(n_{1}\times\cdots\times n_{d})$ into a
matrix, where the multi-linear structure of a tensor is ignored. Without loss
of generality, assume that $\mathcal{X}^{*}$ is unfolded along the first mode,
in which the resulting matrix is
$\mathcal{X}_{(1)}^{*}=A_{1}^{*}(\mathcal{C}_{(1)}^{*}(\otimes_{i=d}^{2}A_{i}^{*})^{T})\in\mathbb{R}_{+}^{n_{1}\times(n_{2}\cdots
n_{d})}$. Then, by applying the matrix based method in [46] to
$(\mathcal{X}_{(1)}^{*})^{T}$, the error bound of the resulting estimator is
on the order of $O(\frac{r_{1}n_{2}n_{3}\cdots
n_{d}+\|A_{1}^{*}\|_{0}}{m}\log(n_{1}\vee(n_{2}\cdots n_{d})))$, which is
larger than that of Theorem 4.4 if $r_{i}$ is much smaller than
$n_{i},i\in[d]$, and $d$ is not too large. In particular, when
$r_{1}=\cdots=r_{d}=r,$ $n_{1}=\cdots=n_{d}=n$, and
$\|A_{1}^{*}\|_{0}=\cdots=\|A_{d}^{*}\|_{0}$, the error bound of the estimator
in [46, Corollary 6] is on the order of
$O(\frac{rn^{d-1}+\|A_{1}^{*}\|_{0}}{m}(d-1)\log(n))$, and the error bound in
Theorem 4.4 is $O(\frac{r^{d}+d\|A_{1}^{*}\|_{0}}{m}\log(n))$, which is
smaller than that of [46] as long as $d\geq 3$ and $r$ is much smaller than
$n$.

###### Remark 4.9

Cao et al. [7] proposed a matrix completion method with the nuclear norm
constraint for Poisson observations and established the error bound of the
estimator of their proposed model. For higher-order tensor completion, without
loss of generality, we unfolding the underlying tensor
$\mathcal{X}^{*}\in\mathbb{R}_{+}^{n_{1}\times n_{2}\times\cdots\times n_{d}}$
into a matrix along the first mode, i.e., $\mathcal{X}_{(1)}^{*}$. In this
case, when $m\geq(n_{1}+n_{2}\cdots n_{d})\log(n_{1}\cdots n_{d})$, the error
bound of the estimator obtained by [7, Theorem 2] is

$\frac{\|\mathcal{X}_{mt}-\mathcal{X}^{*}\|_{F}^{2}}{n_{1}n_{2}\cdots
n_{d}}\leq C\sqrt{\frac{r_{1}(n_{1}+n_{2}\cdots n_{d})}{m}}\log(n_{1}\cdots
n_{d})$ (17)

with high probability, where $\mathcal{X}_{mt}$ is the estimator obtained by
[7] and $C>0$ is a constant. In particular, when $n_{1}=\cdots=n_{d}=n$ and
$r_{1}=\cdots=r_{d}$, the upper bound in (17) reduces to
$O(d\sqrt{\frac{r_{1}(n+n^{d-1})}{m}}\log(n))$, which is larger than
$O(\frac{r_{1}^{d}+\sum_{i=1}^{d}\|A_{i}^{*}\|_{0}}{m}\log(n))$ if $r_{1}$ is
much smaller than $n$ (e.g., $r_{1}<(\frac{n^{d-1}}{4})^{\frac{1}{2d-1}}$).
Moreover, the factor matrix $A_{i}$ is sparser, the error bound in Theorem 4.4
is lower, while the factor matrix $A_{i}$ has not influence on the error bound
in [7, Theorem 2].

###### Remark 4.10

Similar to the observations with additive Laplace noise, we can only compare
the error bound in Theorem 4.4 with that in [54, Proposition 4.3] for third-
order tensors, which has smaller error bound than that of [53, Theorem 3.1]
for low-rank tensor completion with Poisson observations. More detailed
comparisons about the error bounds of the two methods can be referred to [53,
Section IV.C]. We do not consider the sparse case since the factor tensor
based on the algebra framework of tensor-tensor product is not comparable with
the factor matrices in Tucker decomposition. In this case, the error bound of
the estimator in [54, Proposition 4.3] is
$O(\frac{rn_{1}n_{3}+rn_{2}n_{3}}{m}\log(n_{1}\vee n_{2}))$, where $r$ is the
tubal rank of the underlying tensor. By Theorem 4.4, we obtain that the error
bound of the estimator of model (6) is on the order of
$O(\frac{r_{1}r_{2}r_{3}+\sum_{i=1}^{3}r_{i}n_{i}}{m}\log(n_{1}\vee n_{2}\vee
n_{3}))$, which is smaller than that of [54] if $r$ is close to $r_{i}$ and
$r_{i}$ is much smaller than $n_{i}$, $i=1,2,3$.

##  5 Minimax Lower Bounds

In this section, we establish the minimax lower bound of an estimator for the
observations satisfying (5). The accuracy of an estimator
$\widetilde{\mathcal{X}}$ for estimating the true tensor $\mathcal{X}^{*}$ can
be measured in terms of its risk [47], which is defined as

$\frac{\mathbb{E}_{\Omega,\mathcal{Y}_{\Omega}}[\|\widetilde{\mathcal{X}}-\mathcal{X}^{*}\|_{F}^{2}]}{n_{1}n_{2}\cdots
n_{d}}.$

Now we consider a class of tensors with Tucker rank
$\mathbf{r}=(r_{1},\ldots,r_{d})$ satisfying (3), where each factor matrix in
Tucker decomposition obeys $\|A_{i}\|_{0}\leq s_{i},\ i\in[d]$, and the
amplitudes of each entry of the core tensor and factor matrices are not larger
than $1$ and $a_{i}$, respectively. In this case, we define the following set

$\begin{split}\mathfrak{L}(\mathbf{s},\mathbf{r},\mathbf{a}):=\Big{\\{}&\mathcal{X}=\mathcal{C}\times_{1}A_{1}\cdots\times_{d}A_{d}:\mathcal{C}\in\mathbb{R}_{+}^{r_{1}\times\cdots\times
r_{d}},A_{i}\in\mathbb{R}_{+}^{n_{i}\times
r_{i}},0\leq\mathcal{C}_{i_{1}\cdots i_{d}}\leq 1,\|A_{i}\|_{0}\leq s_{i},\\\
&~{}~{}~{}0\leq(A_{i})_{lm}\leq
a_{i},(i_{1},\ldots,i_{d})\in[r_{1}]\times\cdots\times[r_{d}],(l,m)\in[n_{i}]\times[r_{i}],i\in[d]\Big{\\}}.\end{split}$
(18)

The worst-case performance of an estimator $\widetilde{\mathcal{X}}$ of
$\mathcal{X}^{*}$ is measured by the maximum risk on the set
$\mathfrak{L}(\mathbf{s},\mathbf{r},\mathbf{a})$, which is defined as

$\sup_{\mathcal{X}^{*}\in\mathfrak{L}(\mathbf{s},\mathbf{r},\mathbf{a})}\frac{\mathbb{E}_{\Omega,\mathcal{Y}_{\Omega}}[\|\widetilde{\mathcal{X}}-\mathcal{X}^{*}\|_{F}^{2}]}{n_{1}n_{2}\cdots
n_{d}}.$

The estimator, which has the smallest maximum risk among all possible
estimators, is said to achieve the minimax risk [47]. In this case, for the
sparse NTD and completion problem, the minimax risk is defined as

$\inf_{\widetilde{\mathcal{X}}}\sup_{\mathcal{X}^{*}\in\mathfrak{L}(\mathbf{s},\mathbf{r},\mathbf{a})}\frac{\mathbb{E}_{\Omega,\mathcal{Y}_{\Omega}}[\|\widetilde{\mathcal{X}}-\mathcal{X}^{*}\|_{F}^{2}]}{n_{1}n_{2}\cdots
n_{d}}.$ (19)

Next we will estimate the lower bound of (19), which is stated in the
following theorem.

###### Theorem 5.1

Suppose that the probability density function or probability mass function of
any two entries of the noisy observations satisfies

$K(\mathbb{P}_{x},\mathbb{P}_{y})\leq\frac{(x-y)^{2}}{2\mu},$ (20)

where $\mu>0$ depends on the noise distribution. Assume that $r_{i}\leq s_{i}$
in the set $\mathfrak{L}(\mathbf{s},\mathbf{r},\mathbf{a})$ in (18),
$i\in[d]$. Then, for the joint probability density function or probability
mass function of the observations $\mathcal{Y}_{\Omega}$ obeying (5), there
exist two constants $\widetilde{\alpha}_{m},\gamma_{m}>0$ such that

$\begin{split}&\inf_{\widetilde{\mathcal{X}}}\sup_{\mathcal{X}^{*}\in\mathfrak{L}(\mathbf{s},\mathbf{r},\mathbf{a})}\frac{\mathbb{E}_{\Omega,\mathcal{Y}_{\Omega}}[\|\widetilde{\mathcal{X}}-\mathcal{X}^{*}\|_{F}^{2}]}{n_{1}n_{2}\cdots
n_{d}}\\\ \geq\
&\frac{\widetilde{\alpha}_{m}}{2^{d+5}(d+1)}\min\left\\{\prod_{i=1}^{d}\Delta_{i}(s_{i},n_{i})a_{i}^{2},\gamma_{m}^{2}\mu^{2}\left(\frac{r_{1}r_{2}\cdots
r_{d}+\sum_{i=1}^{d}s_{i}}{m}\right)\right\\},\end{split}$

where

$\Delta_{i}(s_{i},n_{i})=\min\left\\{1,\frac{s_{i}}{n_{i}}\right\\},\
i\in[d].$ (21)

For a fixed $d$, Theorem 5.1 shows that the lower bound of the estimator
satisfying (18) is on the order of $O(\frac{r_{1}r_{2}\cdots
r_{d}+\sum_{i=1}^{d}s_{i}}{m})$, which matches to the upper bound in Theorem
4.1 up to a logarithmic factor. This demonstrates the upper bound in Theorem
4.1 is nearly optimal.

###### Remark 5.1

The proof of the minimax risk bound in Theorem 5.1 utilized the technique of
the tools for matrix completion in [26, 43], see also [54]. The key issue of
the proof is to construct the packing sets for the core tensors and sparse
factor matrices, respectively. Then the Varshamov-Gilbert bound [47, Lemma
2.9] is used to determine the Frobenius norm of any two tensors in the set.
Finally, one can get the minimax lower bound by the minimax analysis
techniques in [47, Theorem 2.5].

###### Remark 5.2

Based on Theorem 5.1, for the explicit lower bounds of the minimax risk with
special noise models in Section 4, we just need to specify $\mu$ in (20)
except for Poisson observations. In particular, for the observations of sparse
NTD and completion with additive Gaussian noise and additive Laplace noise, we
can set $\mu=\sigma$ and $\mu=\tau$, respetively. For Poisson observations, we
need each entry of the underlying tensor to be strictly positive in the
construction of the packing set. Similar results can be found in [43, 54]. For
brevity, we omit the detailed proof here for sparse NTD and completion with
Poisson observations.

##  6 Optimization Algorithm

In this section, we design an ADMM based algorithm [6] to solve problem (6).
Notice that the feasible set $\Pi$ in (7) is discrete, which leads to the
difficulty for the algorithm design. Similar to [24, 54], we drop the discrete
assumption of the core tensor and factor matrices in order to use continuous
optimization techniques, which may be justified by choosing a very large value
of $\tau$ and by noting that continuous optimization algorithms use finite
precision arithmetic when executed on a computer. Hence, we consider to solve
the following problem:

$\begin{split}\min_{\mathcal{X},\mathcal{C},A_{i}}\ &-\log
p_{\mathcal{X}_{\Omega}}(\mathcal{Y}_{\Omega})+\sum_{i=1}^{d}\lambda_{i}\|A_{i}\|_{0}\\\
{\textup{s.t.}}\
&\mathcal{X}=\mathcal{C}\times_{1}A_{1}\cdots\times_{d}A_{d},0\leq\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}\leq
c,(i_{1},i_{2},\ldots,i_{d})\in[n_{1}]\times[n_{2}]\times\cdots\times[n_{d}],\\\
\ &0\leq\mathcal{C}_{i_{1}i_{2}\cdots i_{d}}\leq
1,(i_{1},i_{2},\ldots,i_{d})\in[r_{1}]\times[r_{2}]\times\cdots\times[r_{d}],0\leq(A_{i})_{ml}\leq
a_{i},\\\ \ &(m,l)\in[n_{i}]\times[r_{i}],i\in[d].\end{split}$ (22)

Let

$\begin{split}&\Xi_{1}:=\left\\{\mathcal{X}\in\mathbb{R}_{+}^{n_{1}\times
n_{2}\times\cdots\times n_{d}}:0\leq\mathcal{X}_{i_{1}i_{2}\cdots i_{d}}\leq
c,(i_{1},i_{2},\ldots,i_{d})\in[n_{1}]\times[n_{2}]\times\cdots\times[n_{d}]\right\\},\\\
&\Xi_{2}:=\left\\{\mathcal{C}\in\mathbb{R}_{+}^{r_{1}\times
r_{2}\times\cdots\times r_{d}}:0\leq\mathcal{C}_{i_{1}i_{2}\cdots i_{d}}\leq
1,(i_{1},i_{2},\ldots,i_{d})\in[r_{1}]\times[r_{2}]\times\cdots\times[r_{d}]\right\\},\\\
&\Psi_{i}:=\left\\{A_{i}\in\mathbb{R}_{+}^{n_{i}\times
r_{i}}:0\leq(A_{i})_{ml}\leq
a_{i},(m,l)\in[n_{i}]\times[r_{i}]\right\\},i\in[d].\end{split}$

Then problem (22) can be rewritten as

$\begin{split}\min_{\mathcal{X},\mathcal{C},A_{i}}\ &-\log
p_{\mathcal{X}_{\Omega}}(\mathcal{Y}_{\Omega})+\sum_{i=1}^{d}\lambda_{i}\|A_{i}\|_{0}+\delta_{\Xi_{1}}(\mathcal{X})+\delta_{\Xi_{1}}(\mathcal{C})+\sum_{i=1}^{d}\delta_{\Psi_{i}}(A_{i})\\\
{\textup{s.t.}}\
&\mathcal{X}=\mathcal{C}\times_{1}A_{1}\cdots\times_{d}A_{d},\end{split}$ (23)

where $\delta(\cdot)$ denotes the indicator function over a set. Let
$\mathcal{X}=\mathcal{Z},\mathcal{C}=\mathcal{B},A_{i}=H_{i},A_{i}=S_{i},i\in[d]$.
Then problem (23) is equivalent to

$\begin{split}\min\ &-\log
p_{\mathcal{X}_{\Omega}}(\mathcal{Y}_{\Omega})+\sum_{i=1}^{d}\lambda_{i}\|H_{i}\|_{0}+\delta_{\Xi_{1}}(\mathcal{Z})+\delta_{\Xi_{2}}(\mathcal{B})+\sum_{i=1}^{d}\delta_{\Psi_{i}}(S_{i})\\\
{\textup{s.t.}}\
&\mathcal{X}=\mathcal{C}\times_{1}A_{1}\cdots\times_{d}A_{d},\mathcal{X}=\mathcal{Z},\mathcal{C}=\mathcal{B},A_{i}=H_{i},A_{i}=S_{i},i\in[d].\end{split}$
(24)

The augmented Lagrangian function associated with problem (24) is given by

$\begin{split}&L(\mathcal{X},\mathcal{C},A_{i},\mathcal{Z},\mathcal{B},H_{i},S_{i},\mathcal{T}_{i},M_{i},N_{i})\\\
=&-\log
p_{\mathcal{X}_{\Omega}}(\mathcal{Y}_{\Omega})+\sum_{i=1}^{d}\lambda_{i}\|H_{i}\|_{0}+\delta_{\Xi_{1}}(\mathcal{Z})+\delta_{\Xi_{2}}(\mathcal{B})+\sum_{i=1}^{d}\delta_{\Psi_{i}}(S_{i})\\\
&+\langle\mathcal{T}_{1},\mathcal{X}-\mathcal{C}\times_{1}A_{1}\cdots\times_{d}A_{d}\rangle+\langle\mathcal{T}_{2},\mathcal{X}-\mathcal{Z}\rangle+\langle\mathcal{T}_{3},\mathcal{C}-\mathcal{B}\rangle\\\
&+\sum_{i=1}^{d}\langle M_{i},A_{i}-H_{i}\rangle+\sum_{i=1}^{d}\langle
N_{i},A_{i}-S_{i}\rangle\\\
&+\frac{\beta_{1}}{2}\|\mathcal{X}-\mathcal{C}\times_{1}A_{1}\cdots\times_{d}A_{d}\|_{F}^{2}+\frac{\beta_{2}}{2}\|\mathcal{X}-\mathcal{Z}\|_{F}^{2}+\frac{\beta_{3}}{2}\|\mathcal{C}-\mathcal{B}\|_{F}^{2}\\\
&+\sum_{i=1}^{d}\left(\frac{\rho_{i}}{2}\|A_{i}-H_{i}\|_{F}^{2}+\frac{\alpha_{i}}{2}\|A_{i}-S_{i}\|_{F}^{2}\right),\end{split}$

where $\mathcal{T}_{1},\mathcal{T}_{2},\mathcal{T}_{3},M_{i},N_{i}$ are the
Lagrangian multipliers and
$\beta_{1},\beta_{2},\beta_{3},\rho_{i},\alpha_{i}>0$ are penalty parameters,
$i\in[d]$. The iteration template of ADMM is given as follows:

$\displaystyle\mathcal{X}^{k+1}=\mathop{\text{argmin}}_{\mathcal{X}}L(\mathcal{X},\mathcal{C}^{k},A_{i}^{k},\mathcal{Z}^{k},\mathcal{B}^{k},H_{i}^{k},S_{i}^{k},\mathcal{T}_{i}^{k},M_{i}^{k},N_{i}^{k}),$
(25)
$\displaystyle\mathcal{C}^{k+1}=\mathop{\text{argmin}}_{\mathcal{C}}L(\mathcal{X}^{k+1},\mathcal{C},A_{i}^{k},\mathcal{Z}^{k},\mathcal{B}^{k},H_{i}^{k},S_{i}^{k},\mathcal{T}_{i}^{k},M_{i}^{k},N_{i}^{k}),$
(26) $\displaystyle
A_{i}^{k+1}=\mathop{\text{argmin}}\limits_{A_{i}}L(\mathcal{X}^{k+1},\mathcal{C}^{k+1},A_{i},\mathcal{Z}^{k},\mathcal{B}^{k},H_{i}^{k},S_{i}^{k},\mathcal{T}_{i}^{k},M_{i}^{k},N_{i}^{k}),$
(27)
$\displaystyle\mathcal{Z}^{k+1}=\mathop{\text{argmin}}\limits_{\mathcal{Z}}L(\mathcal{X}^{k+1},\mathcal{C}^{k+1},A_{i}^{k+1},\mathcal{Z},\mathcal{B}^{k},H_{i}^{k},S_{i}^{k},\mathcal{T}_{i}^{k},M_{i}^{k},N_{i}^{k}),$
(28)
$\displaystyle\mathcal{B}^{k+1}=\mathop{\text{argmin}}\limits_{\mathcal{B}}L(\mathcal{X}^{k+1},\mathcal{C}^{k+1},A_{i}^{k+1},\mathcal{Z}^{k+1},\mathcal{B},H_{i}^{k},S_{i}^{k},\mathcal{T}_{i}^{k},M_{i}^{k},N_{i}^{k}),$
(29) $\displaystyle
H_{i}^{k+1}=\mathop{\text{argmin}}\limits_{H_{i}}L(\mathcal{X}^{k+1},\mathcal{C}^{k+1},A_{i}^{k+1},\mathcal{Z}^{k+1},\mathcal{B}^{k+1},H_{i},S_{i}^{k},\mathcal{T}_{i}^{k},M_{i}^{k},N_{i}^{k}),$
(30) $\displaystyle
S_{i}^{k+1}=\mathop{\text{argmin}}\limits_{S_{i}}L(\mathcal{X}^{k+1},\mathcal{C}^{k+1},A_{i}^{k+1},\mathcal{Z}^{k+1},\mathcal{B}^{k+1},H_{i}^{k+1},S_{i},\mathcal{T}_{i}^{k},M_{i}^{k},N_{i}^{k}),$
(31)
$\displaystyle\mathcal{T}_{1}^{k+1}=\mathcal{T}_{1}^{k}+\beta_{1}(\mathcal{X}^{k+1}-\mathcal{C}^{k+1}\times_{1}A_{1}^{k+1}\cdots\times_{d}A_{d}^{k+1}),$
(32)
$\displaystyle\mathcal{T}_{2}^{k+1}=\mathcal{T}_{2}^{k}+\beta_{2}(\mathcal{X}^{k+1}-\mathcal{Z}^{k+1}),\
\mathcal{T}_{3}^{k+1}=\mathcal{T}_{3}^{k}+\beta_{3}(\mathcal{C}^{k+1}-\mathcal{B}^{k+1}),$
(33) $\displaystyle M_{i}^{k+1}=M_{i}^{k}+\rho_{i}(A_{i}^{k+1}-H_{i}^{k+1}),\
N_{i}^{k+1}=N_{i}^{k}+\alpha_{i}(A_{i}^{k+1}-S_{i}^{k+1}),i\in[d].$ (34)

Next we give the explicit solutions of each subproblem in the ADMM. Let

$f(\mathcal{X}):=-\log p_{\mathcal{X}_{\Omega}}(\mathcal{Y}_{\Omega}).$

The optimal solution of (25) with respect to $\mathcal{X}$ is given by

$\mathcal{X}^{k+1}=\textup{Prox}_{\frac{1}{\beta_{1}+\beta_{2}}f}\left(\frac{1}{\beta_{1}+\beta_{2}}\left(\beta_{1}\mathcal{C}^{k}\times_{1}A_{1}^{k}\cdots\times_{d}A_{d}^{k}-\mathcal{T}_{1}^{k}+\beta_{2}\mathcal{Z}^{k}-\mathcal{T}_{2}^{k}\right)\right).$
(35)

According to (1), problem (26) is equivalent to

$\min_{\mathcal{C}}\frac{\beta_{1}}{2}\left\|\textup{vec}\left(\mathcal{X}^{k+1}+\frac{1}{\beta_{1}}\mathcal{T}_{1}^{k}\right)-(\otimes_{i=d}^{1}A_{i}^{k})\textup{vec}(\mathcal{C})\right\|_{F}^{2}+\frac{\beta_{3}}{2}\left\|\textup{vec}(\mathcal{C})-\textup{vec}\left(\mathcal{B}^{k}-\frac{1}{\beta_{3}}\mathcal{T}_{3}^{k}\right)\right\|_{F}^{2}.$
(36)

The optimal solution of (36) is given by

$\begin{split}\textup{vec}(\mathcal{C}^{k+1})&=\left(\beta_{1}(\otimes_{i=d}^{1}A_{i}^{k})^{T}(\otimes_{i=d}^{1}A_{i}^{k})+\beta_{2}I\right)^{-1}\mathbf{q}^{k}\\\
&=\left(\beta_{1}(\otimes_{i=d}^{1}(A_{i}^{k})^{T}A_{i}^{k})+\beta_{2}I\right)^{-1}\mathbf{q}^{k},\end{split}$
(37)

where
$\mathbf{q}^{k}:=(\otimes_{i=d}^{1}A_{i}^{k})^{T}\textup{vec}\left(\beta_{1}\mathcal{X}^{k+1}+\mathcal{T}_{1}^{k}\right)+\textup{vec}\left(\beta_{3}\mathcal{B}^{k}-\mathcal{T}_{3}^{k}\right)$
and the second equality holds by [22, Lemma 4.2.10].

The problem (27) can be reformulated as

$\min_{A_{i}}\frac{\beta_{1}}{2}\left\|\mathcal{X}^{k+1}_{(i)}-A_{i}R_{i}^{k}+\frac{1}{\beta_{1}}(\mathcal{T}_{1}^{k})_{(i)}\right\|_{F}^{2}+\frac{\rho_{i}}{2}\left\|A_{i}-H_{i}^{k}+\frac{1}{\rho_{i}}M_{i}^{k}\right\|_{F}^{2}$
(38)

where

$R_{i}^{k}:=\mathcal{C}_{(i)}^{k+1}(A_{d}^{k}\otimes\cdots\otimes
A_{i+1}^{k}\otimes A_{i-1}^{k+1}\otimes\cdots\otimes A_{1}^{k+1})^{T}.$ (39)

The optimal solution of (38) is represented as

$A_{i}^{k+1}=\left(\left(\beta_{1}\mathcal{X}^{k+1}_{(i)}+(\mathcal{T}_{1}^{k})_{(i)}\right)(R_{i}^{k})^{T}+\rho_{i}H_{i}^{k}-M_{i}^{k}\right)\left(\beta_{1}R_{i}^{k}(R_{i}^{k})^{T}+\rho_{i}I\right)^{-1}.$
(40)

The optimal solution of (28) with respect to $\mathcal{Z}$ is given by

$\mathcal{Z}^{k+1}=\arg\min_{\mathcal{Z}}\delta_{\Xi_{1}}(\mathcal{Z})+\frac{\beta_{2}}{2}\left\|\mathcal{Z}-\left(\mathcal{X}^{k+1}+\frac{1}{\beta_{2}}\mathcal{T}_{2}^{k}\right)\right\|_{F}^{2}=\mathbf{P}_{\Xi_{1}}\left(\mathcal{X}^{k+1}+\frac{1}{\beta_{2}}\mathcal{T}_{2}^{k}\right),$
(41)

where $\mathbf{P}_{\Xi_{1}}(\cdot)$ denotes the projection onto a given set
$\Xi_{1}$. Similarly, the optimal solution of (29) with respect to
$\mathcal{B}$ is given by

$\mathcal{B}^{k+1}=\arg\min_{\mathcal{B}}\delta_{\Xi_{2}}(\mathcal{B})+\frac{\beta_{3}}{2}\left\|\mathcal{B}-\left(\mathcal{C}^{k+1}+\frac{1}{\beta_{3}}\mathcal{T}_{3}^{k}\right)\right\|_{F}^{2}=\mathbf{P}_{\Xi_{2}}\left(\mathcal{C}^{k+1}+\frac{1}{\beta_{3}}\mathcal{T}_{3}^{k}\right).$
(42)

The optimal solutions of (30) and (31) are given by

$H_{i}^{k+1}=\arg\min_{H_{i}}\lambda_{i}\|H_{i}\|_{0}+\frac{\rho_{i}}{2}\left\|H_{i}-\left(A_{i}^{k+1}+\frac{1}{\rho_{i}}M_{i}^{k}\right)\right\|_{F}^{2}=\textup{Prox}_{\frac{\lambda_{i}}{\rho_{i}}\|\cdot\|_{0}}\left(A_{i}^{k+1}+\frac{1}{\rho_{i}}M_{i}^{k}\right)$
(43)

and

$S_{i}^{k+1}=\arg\min_{S_{i}}\delta_{\Upsilon_{i}}(S_{i})+\frac{\alpha_{i}}{2}\left\|S_{i}-\left(A_{i}^{k+1}+\frac{1}{\alpha_{i}}N_{i}^{k}\right)\right\|_{F}^{2}=\mathbf{P}_{\Psi_{i}}\left(A_{i}^{k+1}+\frac{1}{\alpha_{i}}N_{i}^{k}\right),$
(44)

where the proximal mapping of the $\ell_{0}$ norm can be implemented in a
point-wise manner, and for any $\lambda>0$, its scalar form is given by (see
[4, Example 6.10], for example)

$\textup{Prox}_{\lambda\|\cdot\|_{0}}(y)=\begin{cases}0,&\text{if
}|y|<\sqrt{2\lambda},\\\ \\{0,y\\},&\text{if }|y|=\sqrt{2\lambda},\\\
y,&\text{if }|y|>\sqrt{2\lambda}.\end{cases}$

The template of the ADMM for solving (24) is presented in Algorithm 1.

Algorithm 1 Alternating Direction Method of Multipliers for Solving Problem
(24).

1:Initialization:
$\mathcal{C}^{0},A_{i}^{0},\mathcal{Z}^{0},\mathcal{B}^{0},\mathcal{T}_{1}^{0},\mathcal{T}_{2}^{0},\mathcal{T}_{3}^{0},M_{i}^{0},i\in[d]$.

2:repeat

3: Step 1. Compute ${\mathcal{X}}^{k+1}$ by (35).

4: Step 2. Update ${\mathcal{C}}^{k+1}$ by (37).

5: Step 3. Compute $A_{i}^{k+1}$ via (40), $i\in[d]$.

6: Step 4. Compute $\mathcal{Z}^{k+1}$ by (41).

7: Step 5. Update $\mathcal{B}^{k+1}$ according to (42).

8: Step 6. Compute $H_{i}^{k+1}$ and $S_{i}^{k+1}$ by (43) and (44),
respectively, $i\in[d]$.

9: Step 7. Update $\mathcal{T}_{1}^{k+1}$, $\mathcal{T}_{2}^{k+1}$,
$\mathcal{T}_{3}^{k+1}$, $M_{i}^{k+1}$, $N_{i}^{k+1}$ by (32), (33) and (34),
respectively.

10:until A stopping condition is satisfied.

In the implementable of Algorithm 1, one needs to compute the proximal mapping
of $\frac{1}{\beta_{1}+\beta_{2}}f$ at $\mathcal{H}$, where

$\mathcal{H}:=\frac{1}{\beta_{1}+\beta_{2}}\left(\beta_{1}\mathcal{C}^{k}\times_{1}A_{1}^{k}\cdots\times_{d}A_{d}^{k}-\mathcal{T}_{1}^{k}+\beta_{2}\mathcal{Z}^{k}-\mathcal{T}_{2}^{k}\right).$
(45)

In particular, for the special noise observation models including additive
Gaussian noise, additive Laplace noise, and Poisson observations, the proximal
mapping in (35) are given in detail in [54].

  * •

Additive Gaussian noise: The joint probability density function satisfies (12)
and then

$\mathcal{X}^{k+1}=\frac{1}{1+\sigma^{2}(\beta_{1}+\beta_{2})}\mathcal{P}_{\Omega}\left(\mathcal{Y}+\sigma^{2}(\beta_{1}+\beta_{2})\mathcal{H}\right)+\mathcal{P}_{\overline{\Omega}}(\mathcal{H}),$

where $\overline{\Omega}$ is the complementary set of $\Omega$ and
$\mathcal{P}_{\Omega}$ is the projection operator onto the index $\Omega$ such
that

$\left(\mathcal{P}_{\Omega}\left(\mathcal{A}\right)\right)_{i_{1}i_{2}\cdots
i_{d}}=\begin{cases}\mathcal{A}_{i_{1}i_{2}\cdots i_{d}},&\textup{if}\
(i_{1},i_{2},\cdots,i_{d})\in\Omega,\\\ 0,&\text{otherwise.}\end{cases}$

  * •

Additive Laplace noise: The joint probability density function satisfies (15)
and then

$\mathcal{X}^{k+1}=\mathcal{P}_{\Omega}\left(\mathcal{Y}+\textup{sign}(\mathcal{H}-\mathcal{Y})\circ\max\left\\{|\mathcal{H}-\mathcal{Y}|-\frac{1}{\tau(\beta_{1}+\beta_{2})},0\right\\}\right)+\mathcal{P}_{\overline{\Omega}}(\mathcal{H}),$

where $\textup{sign}(\cdot)$ and $\circ$ represents the signum function and
the Hadamard point-wise product, respectively.

  * •

Poisson observations: The joint probability density function satisfies (16)
and then

$\mathcal{X}^{k+1}=\frac{1}{2(\beta_{1}+\beta_{2})}\mathcal{P}_{\Omega}\left((\beta_{1}+\beta_{2})\mathcal{H}-\mathbf{1}+\sqrt{((\beta_{1}+\beta_{2})\mathcal{H}-\mathbf{1})^{2}+4(\beta_{1}+\beta_{2})\mathcal{Y}}\right)+\mathcal{P}_{\overline{\Omega}}(\mathcal{H}),$

where $\mathbf{1}$ represents an $n_{1}\times n_{2}\times\cdots\times n_{d}$
tensor with all entries being $1$, the square root and square are performed in
the point-wise manner.

Now the computational complexity of ADMM for solving (24) is given as follows.
If the observations are chosen as additive Gaussian noise, additive Laplace
noise, or Poisson observations, the computation cost of $\mathcal{X}^{k+1}$ is
on the order of
$O(\sum_{j=1}^{d}(\prod_{i=1}^{j}n_{i})(\prod_{i=j}^{d}r_{i}))$, which mainly
follows from the computational cost of (45). The computational cost of
$\mathcal{C}^{k+1}$ is
$O(\sum_{i=1}^{d}r_{i}^{2}n_{i}+\prod_{i=1}^{d}r_{i}^{3}+\sum_{j=1}^{d}(\prod_{i=1}^{j}r_{i})(\prod_{i=j}^{d}n_{i}))$.
In the implementation of $R_{i}^{k}$ in (39), we do not compute the Kronecker
product directly. Let
$\mathcal{K}_{i}:=\mathcal{C}^{k+1}\times_{1}A_{1}^{k+1}\cdots\times_{i-1}A_{i-1}^{k+1}\times_{i+1}A_{i+1}^{k}\times\cdots\times_{d}A_{d}^{k}.$
Then $(\mathcal{K}_{i})_{(i)}=R_{i}^{k}$, which leads to the computation of
$R_{i}^{k}$ via $(\mathcal{K}_{i})_{(i)}$. The computational cost of
$\mathcal{K}_{i}$ is [51, Appendix B]

$O\left(\sum_{j=1}^{i-1}\left(\prod_{t=1}^{j}n_{t}\right)\left(\prod_{t=j}^{d}r_{t}\right)+r_{i}\left(\prod_{t=1}^{i-1}n_{t}\right)\sum_{j=i+1}^{d}\left(\prod_{t={i+1}}^{j}n_{t}\right)\left(\prod_{t=j}^{d}r_{t}\right)\right).$
(46)

For any $i\in[d]$, the computational cost of each factor matrix $A_{i}$ is the
sum of (46) and $O(r_{i}n_{1}\cdots n_{d}+r_{i}^{3}+n_{i}r_{i}^{2}).$ The
computational costs of $\mathcal{Z}^{k+1}$ and $\mathcal{B}^{k+1}$ are
$O(\prod_{i=1}^{d}n_{i})$ and $O(\prod_{i=1}^{d}r_{i})$, respectively. The
computational cost of $H_{i}^{k+1}$ and $S_{i}^{k+1}$ are both
$O(n_{i}r_{i})$. The computational costs of the multipliers
$\mathcal{T}_{1},\mathcal{T}_{2},\mathcal{T}_{3},M_{i},N_{i}$ are
$O(\sum_{j=1}^{d}(\prod_{i=1}^{j}n_{i})(\prod_{i=j}^{d}r_{i}))$. Note that
$r_{i}\leq n_{i},i\in[d]$, and the quantities of (46) and
$O(\sum_{j=1}^{d}(\prod_{i=1}^{j}r_{i})(\prod_{i=j}^{d}n_{i}))$ are similar,
where the detailed discussions can be found in [51, Appendix B]. Therefore,
the computational cost of ADMM in Algorithm 1 is

$O\left(d\sum_{j=1}^{d}\left(\prod_{i=1}^{j}n_{i}\right)\left(\prod_{i=j}^{d}r_{i}\right)+\sum_{j=1}^{d}\left(\prod_{i=1}^{j}r_{i}\right)\left(\prod_{i=j}^{d}n_{i}\right)+\sum_{i=1}^{d}r_{i}^{2}n_{i}+\prod_{i=1}^{d}r_{i}^{3}+\left(\prod_{i=1}^{d}n_{i}\right)\sum_{i=1}^{d}r_{i}\right).$

###### Remark 6.1

Since the constraint and the $\ell_{0}$ norm in the objective of (24) are both
nonconvex, the convergence of ADMM cannot be guaranteed in general. Although
great efforts have been made about the convergence of ADMM for solving
nonconvex problems, the problem with both constraints and objective are
nonconvex, where the objective is not Lipschitz continuous for general loss
function mentioned in Section 4. Therefore, the existing literature about the
convergence of ADMM for nonconvex problems (e.g., see [5, 21, 50]) cannot be
applied to our model directly.

##  7 Numerical Experiments

In this section, we evaluate the effectiveness the proposed sparse NTD and
completion model (SNTDC) using both synthetic data and real image data. We
compare SNTDC with a matrix based method in [46] (denoted by Matrix) and
sparse nonnegative tensor factorization and completion with tensor-tensor
product (SNTFTP) [54]. All experiments are performed in MATLAB 2020a on a
computer with Intel Xeon W-2133 and 32 GB of RAM.

Algorithm 1 will be stopped if the maximum number of iterations $300$ is
reached or the following condition is satisfied

$\frac{\|\mathcal{X}^{k+1}-\mathcal{X}^{k}\|_{F}}{\|\mathcal{X}^{k}\|_{F}}\leq
10^{-4},$

where $\mathcal{X}^{k}$ is the $k$th iteration of Algorithm 1. We use a
sequentially truncated high order singular value decomposition [49] on
$\mathcal{P}_{\Omega}(\mathcal{Y})$ to get the initial values
$\mathcal{C}^{0},A_{i}^{0},i\in[d]$. Other initial values in Algorithm 1 are
the same as $\mathcal{P}_{\Omega}(\mathcal{Y})$, or $\mathcal{C}^{0}$, or
$A_{i}^{0},i\in[d]$ as long as their sizes are the same. For the
regularization parameters $\lambda_{i}$, we set $\lambda_{i}$ to the same for
any $i\in[d]$ and choose them from the set $\\{5,10,50\\}$ to get the best
recovery performance. Moreover, $\beta_{i}$ are set the same for different $i$
and chosen from the set $\\{50,70,100,250,350,450,550\\}$. For $\alpha_{i}$
and $\rho_{i}$, we set them the same in each case for any $i\in[d]$, and
choose them from the set $\\{0.1,0.01,0.001\\}$ to get the best recovery
performance.

###  7.1 Synthetic Data

In this subsection, we test the synthetic data to demonstrate the
effectiveness of the proposed method. The synthetic tensor is generated as
follows: The underlying tensor has the Tucker decomposition
$\mathcal{X}=\mathcal{C}\times_{1}A_{1}\times_{2}A_{2}\cdots\times_{d}A_{d},$
where the core tensor $\mathcal{C}$ is generated by MATLAB command
$\textup{rand}([n_{1},\ldots,n_{d}])$, and the factor matrices $A_{i}$ are
generated by $a_{i}\cdot\textup{sptenrand}([n_{i},r_{i}],\upsilon_{i})$,
$i\in[d]$. Here $\upsilon_{i}$ is the sparse ratio of $A_{i}$, i.e.,
$\upsilon_{i}=\frac{\|A_{i}\|_{0}}{n_{i}r_{i}},i\in[d]$. We set
$c=2\|\mathcal{X}\|_{\infty}$. For these synthetic data, we test 10 trials in
each case and the average result of these trials is the last result.

We first test the third-order tensors with $n_{1}=n_{2}=n_{3}=100$, whose
Tucker rank is $(5,5,5)$. For the observations with additive Gaussian noise
and additive Laplace noise, we set $\sigma^{2}=0.01$ and $\tau=0.01$,
respectively. In Figure 1, we show the relative error versus sampling ratios
of different methods for the observations with additive Gaussian noise,
additive Laplace noise, and Poisson observations, respectively, where
$\upsilon_{i}=0.3,i=1,2,3$. It can be observed that the relative errors of the
matrix based method, SNTFTP, and SNTDC decrease as the sampling ratio
increases. Moreover, the relative errors obtained by SNTDC are smaller than
those obtained by the matrix based method and SNTFTP for the three kinds of
observations.

(a) Gaussian noise (b) Laplace noise (c) Poisson observations

Figure 1: Relative error versus sampling ratio for synthetic tensors with size
$100\times 100\times 100$ and Tucker rank $(5,5,5)$.

(a) Gaussian noise (b) Laplace noise (c) Poisson observations

Figure 2: Relative error versus sampling ratio for synthetic tensors with size
$50\times 50\times 50\times 50$ and Tucker rank $(5,5,5,5)$.

Now we test a fourth-order tensor $50\times 50\times 50\times 50$ with Tucker
rank $(5,5,5,5)$, where $\sigma^{2}=0.01$ for additive Gaussian noise,
$\tau=0.01$ for additive Laplace noise, and $\upsilon_{i}=0.3,i=1,2,3$. Since
the SNTFTP is only effective for third-order tensors, we do not compare with
it for fourth-order tensors in this experiment. In Figure 2, we show the
relative error versus sampling ratio for the observations with additive
Gaussian noise, additive Laplace noise, and Poisson observations,
respectively. We can see that the relative errors of the SNTDC are much
smaller than those of the matrix based method for the observations with
different noise distributions. Similarly, the relative errors decrease as the
number of samples increase for both the matrix based method and SNTDC.

###  7.2 Image Data

In this subsection, we test the Swimmer
dataset111https://stodden.net/Papers.html ($32\times 32\times 256$) [14] and
the Columbia Object Image Library
(COIL-100)222https://www.cs.columbia.edu/CAVE/software/softlib/coil-100.php
for sparse nonnegative Tucker decomposition and completion. For the Swimmer
dataset, the size of each image is $32\times 32$ and there are 256 swimmer
images. Similar to [51], the Tucker rank of the Swimmer dataset is set to
$(24,20,20)$ in Algorithm 1. For the COIL-100, which contains 100 objects, we
resize each image to $64\times 64\times 3$ and choose first 50 images for the
fifty-sixth object, In this case, the size of the resulting tensor is
$64\times 64\times 3\times 50$, whose Tucker rank is set to $(15,15,3,10)$ in
Algorithm 1.

In Figure 3, we show the relative error versus sampling ratio of different
methods for the observations with additive Gaussian noise, additive Laplace
noise, and Poisson observations, respectively, where $\sigma^{2}=0.01$ for
additive Gaussian noise and $\tau=0.01$ for additive Laplace noise. It can be
seen from this figure that the relative errors of the SNTDC are smaller than
those of the matrix based method and SNTFC for different sampling ratios.
Moreover, the performance of SNTF is better than that of the matrix based
method in terms of relative errors for the three kinds of observations. And
the relative errors of the matrix based method, SNTFC, and SNTDC decrease as
the number of samples increases for the observations with additive Gaussian
noise, additive Laplace noise, and Poisson observations.

(a) Gaussian noise (b) Laplace noise (c) Poisson observations

Figure 3: Relative error versus sampling ratio for the Swimmer dataset.

(a) Gaussian noise (b) Laplace noise (c) Poisson observations

Figure 4: Relative error versus sampling ratio for the COIL-100 dataset.

Since the SNTFTP is only effective for third-order tensors, we stack all
images of the COIL-100 based on the color channels into a third-order tensors
along the frontal slices for the SNTFTP, where the size of the resulting
tensor is $64\times 64\times 150$. In Figure 4, we display the relative error
versus sampling ratio of different methods for the COIL-100 dataset, where
$\sigma^{2}=0.05$ for additive Gaussian noise and $\tau=0.05$ for additive
Laplace noise. We can see that the relative errors obtained by the matrix
based method, SNTFTP, and SNTDC decrease as the sampling ratio increases for
different noise observations. And the relative errors obtained by SNTDC are
smaller than those obtained by the matrix based method and SNTFTP. Moreover,
the SNTFTP performs better than the matrix based method in terms of relative
error for the observations with additive Gaussian noise, additive Laplace
noise, and Poisson observations.

##  8 Conclusions and Future Work

In this paper, we have proposed and studied the problem of sparse NTD and
completion with a general class of noise observations, where the underlying
tensor is decomposed into a core tensor and several factor matrices with
nonnegativity and sparsity constraints. Moreover, the loss function is derived
by the maximum likelihood estimation of the noisy observations. The error
bound of the estimator of the SNTDC model is established under a class of
noise distributions. Then the error bounds are specified to some widely used
noise models including additive Gaussian noise, additive Laplace noise, and
Poisson observations. Besides, the minimax lower bound of the observed model
is derived, which matches the upper bound up to a logarithmic factor.
Numerical experiments demonstrate the effectiveness of the proposed sparse NTD
and completion model compared with other methods.

In the experiments, we need to set the Tucker rank of the underlying tensor in
advance, which is unknown in general for real-world datasets. An interesting
direction for future work on this problem is to update the Tucker rank
adaptively. Besides, in the sparse NTD and completion model, the $\ell_{0}$
norm is employed to characterise the sparsity of the factor matrices. However,
the $\ell_{0}$ norm is nonconvex and it is challenge to get the unique and
global optimal solution. Therefore, it is also of great interest to replace
the $\ell_{0}$ norm by the $\ell_{1}$ norm for the factor matrices and then
analyze the error bound of the corresponding model.

The convergence of Algorithm 1 cannot be guaranteed in general since our model
in (24) is nonconvex and multi-block. Future work is devoted to establishing
the convergence of ADMM for our nonconvex multi-block model. Moreover, due to
effectiveness and explanatory information of orthogonal constraints on the
factor matrices for nonnegative Tucker decomposition (c.f. [40]), it will also
be of great interest to extend the error bound of our model to that of
orthogonal nonnegative Tucker decomposition under general loss function.

## Appendix A. Proof of Theorem 4.1

At the beginning, the following lemma is first recorded, which plays a vital
role for the proof of Theorem 4.1.

###### Lemma 8.1

Let $\Upsilon$ be a countable collection of candidate reconstructions
$\mathcal{X}$ of $\mathcal{X}^{*}$ and its penalty
$\emph{pen}(\mathcal{X})\geq 1$ satisfying
$\sum_{\mathcal{X}\in\Upsilon}2^{-\textup{pen}(\mathcal{X})}\leq 1$. For any
integer $2^{d}\leq m\leq n_{1}n_{2}\cdots n_{d}$, let
$\Omega\sim\textup{Bern}(p)$. Moreover, the joint probability density/mass
function of the corresponding observations follows
$p_{\mathcal{X}_{\Omega}^{*}}(\mathcal{Y}_{\Omega})=\prod_{(i,j,k)\in\Omega}p_{\mathcal{X}_{ijk}^{*}}(\mathcal{Y}_{ijk})$,
which are assumed to be conditionally independent on the given $\Omega$. Let
$\gamma$ be a constant satisfying

$\gamma\geq\max_{\mathcal{X}\in\Upsilon}\max_{i_{1},i_{2},\ldots,i_{d}}K\left(p_{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*}}(\mathcal{Y}_{i_{1}i_{2}\cdots
i_{d}})||p_{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}}(\mathcal{Y}_{i_{1}i_{2}\cdots i_{d}})\right).$ (47)

Consider the following complexity penalized maximum likelihood estimator

${\mathcal{X}}^{\mu}\in\arg\min_{\mathcal{X}\in\Upsilon}\left\\{-\log
p_{\mathcal{X}_{\Omega}}(\mathcal{Y}_{\Omega})+\mu\cdot\emph{pen}(\mathcal{X})\right\\}.$
(48)

Then for any $\mu\geq 2\left(1+\frac{2\gamma}{3}\right)\log(2)$, one has

$\begin{split}&\ \frac{\mathbb{E}_{\Omega,\mathcal{Y}_{\Omega}}\left[-2\log
H(p_{{\mathcal{X}}^{\mu}}(\mathcal{Y}),p_{\mathcal{X}^{*}}(\mathcal{Y}))\right]}{n_{1}n_{2}\cdots
n_{d}}\\\ \leq&\
3\cdot\min_{\mathcal{X}\in\Upsilon}\left\\{\frac{K(p_{\mathcal{X}^{*}}(\mathcal{Y})||p_{\mathcal{X}}(\mathcal{Y})}{n_{1}n_{2}\cdots
n_{d}}+\left(\mu+\frac{4\gamma\log(2)}{3}\right)\frac{\emph{pen}(\mathcal{X})}{m}\right\\}+\frac{8\gamma\log(m)}{m},\end{split}$

where the expectation is taken with respect to the joint distribution of
$\Omega$ and $\mathcal{Y}_{\Omega}$.

The proof of Lemma 8.1 can be obtained easily based on the matrix case [46,
Lemma 8] and [31], see also for the cases of CP decomposition in [24, Lemma 1]
and tensor factorization via tensor-tensor product in [54, Lemma 1.1]. The
three steps of proof in [46, Lemma 8] are giving the “good” sample set
characteristics, a conditional error guarantee and some simple conditioning
arguments, which are in point-wise manners in fact for the KL divergence,
negative logarithmic Hellinger affinity, and maximum likelihood estimation. As
a consequence, we can extend them to the tensor case with Tucker decomposition
easily. Here we omit the details.

The proof of Theorem 4.1 follows the line of the proof of [46, Theorem 1], see
also the proofs of [41, Theorem 3] and [54, Theorem 4.1]. The main technique
of this proof is the well-known Kraft-McMillan inequality [35, 28]. The
penalty of the underlying tensor $\mathcal{X}$ is constructed by the codes of
the core tensors and sparse factor matrices, where the underlying tensor has
the Tucker decomposition form with the core tensor being nonnegative and the
factor matrices being nonnegative and sparse. Next we return to the proof of
Theorem 4.1.

Based on the result in Lemma 8.1, one should define the penalties
$\textup{pen}(\mathcal{X})\geq 1$ in the candidate reconstructions
$\mathcal{X}$ of $\mathcal{X}^{*}$ such that the penalties $\mathcal{X}$ in
$\Upsilon$ satisfy

$\sum_{\mathcal{X}\in\Upsilon}2^{-\textup{pen}(\mathcal{X})}\leq 1.$ (49)

The condition (49) is the well-known Kraft-McMillan inequality for coding
elements of $\Upsilon$ with an alphabet of size $2$, which is satisfied
automatically if we choose the penalties to be code lengths for the unique
decodable binary code for the elements $\mathcal{X}\in\Upsilon$ [28, 35], see
also [11, Section 5]. This also gives the constructions of the penalties.

Let
$\Upsilon_{1}:=\\{\mathcal{X}=\mathcal{C}\times_{1}A_{1}\cdots\times_{d}A_{d}:\
\mathcal{C}\in\mathfrak{C},\ A_{i}\in\mathfrak{B}_{i}\\},$ where
$\mathfrak{C}$ and $\mathfrak{B}_{i},i\in[d]$ are the same as those
constructed in Section 3. Note that $\Upsilon_{1}\subseteq\Upsilon$ by the
construction of $\Upsilon_{1}$. Now we consider the discretized core tensor
$\mathcal{C}\in\mathfrak{C}$ and sparse factor matrices
$A_{i}\in\mathfrak{B}_{i},i\in[d]$.

  * (1)

We encode the amplitude of each element of $\mathcal{C}$ using
$\log_{2}(\tau)$ bits, where $\tau$ is defined as (8). Then a total of
$r_{1}r_{2}\cdots r_{d}\log_{2}({\tau})$ bits are used to encode the core
tensor $\mathcal{C}$.

  * (2)

Let $\zeta_{i}:=2^{\lceil\log_{2}(r_{i}n_{i})\rceil},i\in[d]$. We encode each
nonzero element of $A_{i}$ using $\log_{2}(\zeta_{i})$ to denote the location
and $\log_{2}(\tau)$ bits for its amplitude. In this case, a total of
$\|A_{i}\|_{0}(\log_{2}(\tau)+\log_{2}(\zeta_{i}))$ bits is used to encode the
sparse factor matrix $A_{i}$.

  * (3)

Finally, we can assign each $\mathcal{X}\in\Upsilon_{1}$ whose code length
satisfies

$\textup{pen}(\mathcal{X})=r_{1}r_{2}\cdots
r_{d}\log_{2}({\tau})+\sum_{i=1}^{d}\|A_{i}\|_{0}(\log_{2}(\tau)+\log_{2}(\zeta_{i})).$

By the constructions, we know that such codes are uniquely decodable, which
implies that
$\sum_{\mathcal{X}\in\Upsilon_{1}}2^{-\textup{pen}(\mathcal{X})}\leq 1$ [28,
35]. Since $\Upsilon\subseteq\Upsilon_{1}$, one has

$\sum_{\mathcal{X}\in\Upsilon}2^{-\textup{pen}(\mathcal{X})}\leq\sum_{\mathcal{X}\in\Upsilon_{1}}2^{-\textup{pen}(\mathcal{X})}\leq
1.$

Consequently, by Lemma 8.1, we get that the estimator $\mathcal{X}^{\mu}$ in
the following

$\begin{split}\mathcal{X}^{\mu}&\in\arg\min_{\mathcal{X}\in\Upsilon}\left\\{-\log
p_{\mathcal{X}_{\Omega}}(\mathcal{Y}_{\Omega})+\mu\cdot\textup{pen}(\mathcal{X})\right\\}\\\
&=\arg\min_{\mathcal{X}\in\Upsilon}\left\\{-\log
p_{\mathcal{X}_{\Omega}}(\mathcal{Y}_{\Omega})+\mu\sum_{i=1}^{d}\|A_{i}\|_{0}(\log_{2}(\tau)+\log_{2}(\zeta_{i}))\right\\},\end{split}$
(50)

satisfies

$\begin{split}&\frac{\mathbb{E}_{\Omega,\mathcal{Y}_{\Omega}}\left[-2\log
H(p_{{\mathcal{X}}^{\mu}},p_{\mathcal{X}^{*}})\right]}{n_{1}n_{2}\cdots
n_{d}}\leq\frac{8\gamma\log(m)}{m}+3\cdot\min_{\mathcal{X}\in\Upsilon}\Biggl{\\{}\frac{K(p_{\mathcal{X}^{*}}(\mathcal{Y})||p_{\mathcal{X}}(\mathcal{Y})}{n_{1}n_{2}\cdots
n_{d}}+\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left(\mu+\frac{4\gamma\log(2)}{3}\right)\frac{r_{1}r_{2}\cdots
r_{d}\log_{2}(\tau)+\sum_{i=1}^{d}\|A_{i}\|_{0}(\log_{2}(\tau)+\log_{2}(\zeta_{i}))}{m}\Biggl{\\}},\end{split}$

where $\mu\geq 2\left(1+\frac{2\gamma}{3}\right)\log(2)$ and $\gamma$
satisfies (47).

Let $\lambda_{i}=\mu(\log_{2}(\tau)+\log_{2}(\zeta_{i})),i\in[d]$. Note that
$r_{i}\leq n_{i}$ and

$\log_{2}(\tau)+\log_{2}(\zeta_{i})\leq
2(\beta+2)\frac{\log(n_{m})}{\log(2)}.$ (51)

Therefore, for any $\lambda_{i}\geq
4(\beta+2)(1+\frac{2\gamma}{3})\log(n_{m}),i\in[d],$ the estimator of

$\begin{split}{\mathcal{X}}^{\lambda}\in\arg\min_{\mathcal{X}\in\Upsilon}\left\\{-\log
p_{\mathcal{X}_{\Omega}}(\mathcal{Y}_{\Omega})+\sum_{i=1}^{d}\lambda_{i}\|A_{i}\|_{0}\right\\}\end{split}$
(52)

satisfies

$\begin{split}&\ \frac{\mathbb{E}_{\Omega,\mathcal{Y}_{\Omega}}\left[-2\log
H(p_{{\mathcal{X}}^{\lambda}}(\mathcal{Y}),p_{\mathcal{X}^{*}}(\mathcal{Y}))\right]}{n_{1}n_{2}\cdots
n_{d}}\leq\frac{8\gamma\log(m)}{m}+3\cdot\min_{\mathcal{X}\in\Upsilon}\Biggl{\\{}\frac{K(p_{\mathcal{X}^{*}}(\mathcal{Y})||p_{\mathcal{X}}(\mathcal{Y}))}{n_{1}n_{2}\cdots
n_{d}}\\\ &\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\left(\max_{i}\\{\lambda_{i}\\}+\frac{8\gamma(\beta+2)\log(n_{m})}{3}\right)\frac{r_{1}r_{2}\cdots
r_{d}+\sum_{i=1}^{d}\|A_{i}\|_{0}}{m}\Biggl{\\}},\end{split}$

where the inequality holds by (51). This completes the proof. $\Box$

## Appendix B. Proof of Theorem 4.2

First, we establish an upper bound of the tensor infinity norm between the
underlying tensor $\mathcal{X}^{*}$ and its closest surrogate in $\Upsilon$,
where
$\mathcal{X}^{*}=\mathcal{C}^{*}\times_{1}A_{1}^{*}\times_{2}A_{2}^{*}\cdots\times_{d}A_{d}^{*}$,
and $\mathcal{C}^{*},A_{i}^{*}\in[d]$ are defined as (4). The following
estimate will be useful in the sequel.

###### Lemma 8.2

Let
$\mathcal{X}^{*}_{s}=\mathcal{C}^{s}\times_{1}A_{1}^{s}\cdots\times_{d}A_{d}^{s}$,
where the entries of $\mathcal{C}^{s}$ are the closest discretized surrogates
of the entries of $\mathcal{C}^{*}$ in $\Upsilon$, and the entries of
$A_{i}^{s}$ are the closest discretized surrogates of the entries of
$A_{i}^{*}$ in $\Upsilon$, $i\in[d]$. Then it holds that

$\|\mathcal{X}^{*}_{s}-\mathcal{X}^{*}\|_{\infty}\leq\frac{2^{d+1}-1}{\tau-1}\prod_{i=1}^{d}(a_{i}r_{i}).$

where $\tau$ is defined as (8).

Proof. Let $\mathcal{C}^{s}=\mathcal{C}^{*}+\triangle\mathcal{C}$ and
$A_{i}^{s}=A_{i}^{*}+\triangle A_{i},i\in[d]$. By the definition of
$\mathcal{C}^{*}$ and $A_{i}^{*},i\in[d]$ in (4), we get that
$\|\triangle\mathcal{C}\|_{\infty}\leq\frac{1}{\tau-1}$ and $\|\triangle
A_{i}\|_{\infty}\leq\frac{a_{i}}{\tau-1}$. It follows from (1) that

$\begin{split}\mathcal{X}^{*}_{s}-\mathcal{X}^{*}=&\
(\otimes_{i=d}^{1}(A_{i}^{*}+\triangle
A_{i}))\textup{vec}(\mathcal{C}^{*}+\triangle\mathcal{C})-(\otimes_{i=d}^{1}A_{i}^{*})\textup{vec}(\mathcal{C}^{*})\\\
=&\ (\otimes_{i=d}^{1}(A_{i}^{*}+\triangle
A_{i}))\textup{vec}(\mathcal{C}^{*})+(\otimes_{i=d}^{1}(A_{i}^{*}+\triangle
A_{i}))(\triangle\mathcal{C})-(\otimes_{i=d}^{1}A_{i}^{*})\textup{vec}(\mathcal{C}^{*}).\end{split}$
(53)

Furthermore, we have

$\begin{split}\ &\otimes_{i=d}^{1}(A_{i}^{*}+\triangle A_{i})\\\ =\
&A_{d}^{*}\otimes\cdots\otimes A_{1}^{*}+\triangle A_{d}\otimes
A_{d-1}^{*}\otimes\cdots\otimes A_{1}^{*}+\cdots+A_{d}^{*}\otimes\cdots\otimes
A_{2}^{*}\otimes\triangle A_{1}\\\ \ &+\triangle A_{d}\otimes\triangle
A_{d-1}\otimes A_{d-3}^{*}\otimes\cdots\otimes
A_{1}^{*}+\cdots+A_{d}^{*}\otimes\cdots\otimes A_{3}^{*}\otimes\triangle
A_{2}\otimes\triangle A_{1}\\\ \ &+\cdots\\\ \ &+\triangle
A_{d}\otimes\cdots\otimes\triangle A_{2}\otimes
A_{1}^{*}+\cdots+A_{d}^{*}\otimes\triangle
A_{d-1}\otimes\cdots\otimes\triangle A_{1}\\\ \ &+\triangle
A_{d}\otimes\cdots\otimes\triangle A_{1}.\end{split}$

Since $0\leq(A_{i}^{*})_{lm}\leq a_{i}$ and
$0\leq\mathcal{C}^{*}_{i_{1}i_{2}\cdots i_{d}}\leq 1$,
$\forall(l,m)\in[n_{i}]\times[r_{i}],(i_{1},i_{2},\ldots,i_{d})\in[r_{1}]\times[r_{2}]\times\cdots\times[r_{d}],i\in[d]$,
we have that

$\begin{split}&\|(A_{d}^{*}\otimes\cdots\otimes
A_{1}^{*})\textup{vec}(\mathcal{C}^{*})\|_{\infty}\leq\prod_{i=1}^{d}(a_{i}r_{i}),\\\
&\|(\triangle A_{d}\otimes A_{d-1}^{*}\otimes\cdots\otimes
A_{1}^{*}+\cdots+A_{d}^{*}\otimes\cdots\otimes A_{2}^{*}\otimes\triangle
A_{1})\textup{vec}(\mathcal{C}^{*})\|_{\infty}\leq\frac{C_{d}^{1}}{\tau-1}\prod_{i=1}^{d}(a_{i}r_{i}),\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\vdots\\\
&\|(A_{d}^{*}\otimes\triangle A_{d-1}\otimes\cdots\otimes\triangle
A_{1})\textup{vec}(\mathcal{C}^{*})\|_{\infty}\leq\frac{C_{d}^{d-1}}{(\tau-1)^{d-1}}\prod_{i=1}^{d}(a_{i}r_{i}),\\\
&\|(\triangle A_{d}\otimes\cdots\otimes\triangle
A_{1})\textup{vec}(\mathcal{C}^{*})\|_{\infty}\leq\frac{1}{(\tau-1)^{d}}\prod_{i=1}^{d}(a_{i}r_{i}),\end{split}$
(54)

which yields that

$\|(\otimes_{i=d}^{1}(A_{i}^{*}+\triangle
A_{i}))\textup{vec}(\mathcal{C}^{*})-(\otimes_{i=d}^{1}A_{i}^{*})\textup{vec}(\mathcal{C}^{*})\|_{\infty}\leq\left(1+\frac{1}{\tau-1}\right)^{d}\prod_{i=1}^{d}(a_{i}r_{i})-\prod_{i=1}^{d}(a_{i}r_{i}).$
(55)

Similarly, by $\|\triangle\mathcal{C}\|_{\infty}\leq\frac{1}{\tau-1}$, we can
easily get that

$\|(\otimes_{i=d}^{1}(A_{i}^{*}+\triangle
A_{i}))(\triangle\mathcal{C})\|_{\infty}\leq\frac{1}{\tau-1}\left(1+\frac{1}{\tau-1}\right)^{d}\prod_{i=1}^{d}(a_{i}r_{i}).$
(56)

Note that

$\begin{split}\left(1+\frac{1}{\tau-1}\right)^{d+1}-1&=C_{d+1}^{1}\frac{1}{\tau-1}+C_{d+1}^{2}\frac{1}{(\tau-1)^{2}}+\cdots+C_{d+1}^{d+1}\frac{1}{(\tau-1)^{d+1}}\\\
&=\frac{1}{\tau-1}\left(C_{d+1}^{1}+C_{d+1}^{2}\frac{1}{\tau-1}+\cdots+C_{d+1}^{d+1}\frac{1}{(\tau-1)^{d}}\right)\\\
&\leq\frac{1}{\tau-1}\left(C_{d+1}^{1}+C_{d+1}^{2}+\cdots+C_{d+1}^{d+1}\right)=\frac{2^{d+1}-1}{\tau-1},\end{split}$
(57)

where the inequality holds by $\tau\geq 2$. Therefore, combining (53), (54),
(55) and (56), we obtain that

$\begin{split}\|\mathcal{X}^{*}_{s}-\mathcal{X}^{*}\|_{\infty}&\leq\left(1+\frac{1}{\tau-1}\right)^{d}\prod_{i=1}^{d}(a_{i}r_{i})-\prod_{i=1}^{d}(a_{i}r_{i})+\frac{1}{\tau-1}\left(1+\frac{1}{\tau-1}\right)^{d}\prod_{i=1}^{d}(a_{i}r_{i})\\\
&=\left(\left(1+\frac{1}{\tau-1}\right)^{d+1}-1\right)\prod_{i=1}^{d}(a_{i}r_{i})\\\
&\leq\frac{2^{d+1}-1}{\tau-1}\prod_{i=1}^{d}(a_{i}r_{i}),\end{split}$ (58)

where the last inequality follows form (57). This furnishes the desired
statement. $\Box$

Now we return to the proof of Theorem 4.2. Recall that the negative
logarithmic Hellinger affinity and KL divergence for additive Gaussian noise
are (see, e.g., [46, 54])

$-\log(H(p_{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}}(\mathcal{Y}_{i_{1}i_{2}\cdots
i_{d}}),p_{(\mathcal{X}^{*})_{i_{1}i_{2}\cdots
i_{d}}}(\mathcal{Y}_{i_{1}i_{2}\cdots
i_{d}})))=\frac{(\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}-(\mathcal{X}^{*})_{i_{1}i_{2}\cdots i_{d}})^{2}}{8\sigma^{2}}$

and

$K\left(p_{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*}}(\mathcal{Y}_{i_{1}i_{2}\cdots
i_{d}})||p_{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}}(\mathcal{Y}_{i_{1}i_{2}\cdots
i_{d}})\right)=\frac{(\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}-\mathcal{X}_{i_{1}i_{2}\cdots i_{d}}^{*})^{2}}{2\sigma^{2}}.$

Therefore, we obtain that

$-2\log
H(p_{{\mathcal{X}}^{\lambda}}(\mathcal{Y}),p_{\mathcal{X}^{*}}(\mathcal{Y}))=\frac{\|\mathcal{X}^{\lambda}-\mathcal{X}^{*}\|_{F}^{2}}{4\sigma^{2}}.$
(59)

Since $\mathcal{X}\in\Upsilon$ and $0\leq\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*}\leq\frac{c}{2}$, we can take $\gamma=\frac{c^{2}}{2\sigma^{2}}$ in
(47). Therefore, it follows from Theorem 4.1 that

$\begin{split}&\
\frac{\mathbb{E}_{\Omega,\mathcal{Y}_{\Omega}}\left[\|\mathcal{X}^{\lambda}-\mathcal{X}^{*}\|_{F}^{2}\right]}{n_{1}n_{2}\cdots
n_{d}}\leq\frac{16c^{2}\log(m)}{m}+12\sigma^{2}\cdot\min_{\mathcal{X}\in\Upsilon}\Biggl{\\{}\frac{K(p_{\mathcal{X}^{*}}(\mathcal{Y})||p_{\mathcal{X}}(\mathcal{Y}))}{n_{1}n_{2}\cdots
n_{d}}\\\ &\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\left(\max_{i}\\{\lambda_{i}\\}+\frac{4c^{2}(\beta+2)\log(n_{m})}{3\sigma^{2}}\right)\frac{r_{1}r_{2}\cdots
r_{d}+\sum_{i=1}^{d}\|A_{i}\|_{0}}{m}\Biggl{\\}}.\end{split}$ (60)

By the definition of $\beta$ in (10), we know that

$\begin{split}\tau=2^{\lceil\log_{2}(n_{m})^{\beta}\rceil}&\geq
2^{\beta\log_{2}(n_{m})}\geq
2^{\log_{2}(n_{m})+\log_{2}\left(\frac{(2^{d+1}-1)\sqrt{d}r_{1}r_{2}\cdots
r_{d}a_{1}a_{2}\cdots a_{d}}{c\sqrt{n_{m}}}+1\right)}\\\
&\geq\frac{(2^{d+1}-1)\sqrt{dn_{m}}r_{1}r_{2}\cdots r_{d}a_{1}a_{2}\cdots
a_{d}}{c}+1.\end{split}$ (61)

For any $\mathcal{X}^{s}$ constructed in Lemma 8.2, we have

$\|\mathcal{X}^{s}-\mathcal{X}^{*}\|_{\infty}\leq\frac{c}{\sqrt{dn_{m}}}\leq\frac{c}{2},$

where the last inequality follows from $d,n_{m}\geq 2$. Therefore,
$\|\mathcal{X}^{s}\|_{\infty}\leq\|\mathcal{X}^{*}\|_{\infty}+\frac{c}{2}\leq
c$, which implies $\mathcal{X}^{s}\in\Upsilon$. As a consequence, we get that

$\begin{split}\min_{\mathcal{X}\in\Upsilon}\Biggl{\\{}\frac{K(p_{\mathcal{X}^{*}}(\mathcal{Y})||p_{\mathcal{X}}(\mathcal{Y}))}{n_{1}n_{2}\cdots
n_{d}}\Biggl{\\}}&=\min_{\mathcal{X}\in\Upsilon}\Biggl{\\{}\frac{\|\mathcal{X}-\mathcal{X}^{*}\|_{F}^{2}}{2\sigma^{2}n_{1}n_{2}\cdots
n_{d}}\Biggl{\\}}\\\
&\leq\frac{\|\mathcal{X}^{s}-\mathcal{X}^{*}\|_{F}^{2}}{2\sigma^{2}n_{1}n_{2}\cdots
n_{d}}\leq\frac{\|\mathcal{X}^{s}-\mathcal{X}^{*}\|_{\infty}^{2}}{2\sigma^{2}}\\\
&\leq\frac{c^{2}}{2\sigma^{2}dn_{m}}\leq\frac{c^{2}}{2\sigma^{2}m},\end{split}$
(62)

where the third inequality holds by (61) and Lemma 8.2, and the last
inequality holds by the fact that $m\leq dn_{m}$.

Moreover, by the construction of $A_{i}^{s}$ in Lemma 8.2, we know that
$\|A_{i}^{s}\|_{0}=\|A_{i}^{*}\|_{0},\forall i\in[d]$. Therefore, plugging
(62) and (11) into (60), we have

$\begin{split}\frac{\mathbb{E}_{\Omega,\mathcal{Y}_{\Omega}}\left[\|\mathcal{X}^{\lambda}-\mathcal{X}^{*}\|_{F}^{2}\right]}{n_{1}n_{2}\cdots
n_{d}}\leq&\ \frac{22c^{2}\log(m)}{m}\\\ &\
+16(\beta+2)(2c^{2}+3\sigma^{2})\left(\frac{r_{1}r_{2}\cdots
r_{d}+\sum_{i=1}^{d}\|A_{i}^{*}\|_{0}}{m}\right)\log(n_{m}),\end{split}$

which leads to the desired conclusion. $\Box$

## Appendix C. Proof of Theorem 4.3

Proof. For the observations with additive Laplace noise, by [54, Lemma 3.1],
we know that

$\begin{split}&\ K\left(p_{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*}}(\mathcal{Y}_{i_{1}i_{2}\cdots
i_{d}})||p_{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}}(\mathcal{Y}_{i_{1}i_{2}\cdots i_{d}})\right)\\\ =&\
\frac{|\mathcal{X}_{i_{1}i_{2}\cdots i_{d}}-\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*}|}{\tau}-1+\exp\left(-\frac{|\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}-\mathcal{X}_{i_{1}i_{2}\cdots i_{d}}^{*}|}{\tau}\right)\\\ \leq&\
\frac{(\mathcal{X}_{i_{1}i_{2}\cdots i_{d}}-\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*})^{2}}{2\tau^{2}},\end{split}$

where the inequality follows from the fact that $e^{-x}\leq
1-x+\frac{x^{2}}{2}$ for any $x>0$. Therefore, we can choose
$\gamma=\frac{c^{2}}{2\tau^{2}}$ in (47).

Let $f(t):=\log(1+\frac{t}{2\tau})$, where $t\in[0,c]$. Then by Taylor’s
expansion, we get that

$f(t)=f(0)+f^{\prime}(0)t+\frac{f^{\prime\prime}(\xi)}{2}t^{2},$ (63)

where $\xi\in[0,c]$. By [54, Lemma 3.1], we obtain that

$\begin{split}&\ -2\log(H(p_{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}}(\mathcal{Y}_{i_{1}i_{2}\cdots
i_{d}}),p_{(\mathcal{X}^{*})_{i_{1}i_{2}\cdots
i_{d}}}(\mathcal{Y}_{i_{1}i_{2}\cdots i_{d}})))\\\ =&\
\frac{|\mathcal{X}_{i_{1}i_{2}\cdots i_{d}}-\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*}|}{\tau}-2\log\left(1+\frac{|\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}-\mathcal{X}_{i_{1}i_{2}\cdots i_{d}}^{*}|}{2\tau}\right)\\\ \geq&\
\frac{(\mathcal{X}_{i_{1}i_{2}\cdots i_{d}}-\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*})^{2}}{(2\tau+c)^{2}}\end{split}$

where the inequality holds by letting $t=|\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}-\mathcal{X}_{i_{1}i_{2}\cdots i_{d}}^{*}|$ and the fact that
$f^{\prime\prime}(\xi)\leq-\frac{1}{(2\tau+c)^{2}}$ in (63). Therefore, by
Theorem 4.1, we can deduce

$\begin{split}&\
\frac{\mathbb{E}_{\Omega,\mathcal{Y}_{\Omega}}\left[\|\mathcal{X}^{\lambda}-\mathcal{X}^{*}\|_{F}^{2}\right]}{n_{1}n_{2}\cdots
n_{d}}\leq\frac{8c^{2}(2\tau+c)^{2}\log(m)}{2\tau^{2}m}+3(2\tau+c)^{2}\cdot\min_{\mathcal{X}\in\Upsilon}\Biggl{\\{}\frac{K(p_{\mathcal{X}^{*}}(\mathcal{Y})||p_{\mathcal{X}}(\mathcal{Y}))}{n_{1}n_{2}\cdots
n_{d}}\\\ &\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\left(\max_{i}\\{\lambda_{i}\\}+\frac{4c^{2}(\beta+2)\log(n_{m})}{3\tau^{2}}\right)\frac{r_{1}r_{2}\cdots
r_{d}+\sum_{i=1}^{d}\|A_{i}\|_{0}}{m}\Biggl{\\}}.\end{split}$ (64)

Then by a similar argument of (62) in the proof of Theorem 4.2, we know that

$\begin{split}\min_{\mathcal{X}\in\Upsilon}\Biggl{\\{}\frac{K(p_{\mathcal{X}^{*}}(\mathcal{Y})||p_{\mathcal{X}}(\mathcal{Y}))}{n_{1}n_{2}\cdots
n_{d}}\Biggl{\\}}&\leq\min_{\mathcal{X}\in\Upsilon}\Biggl{\\{}\frac{\|\mathcal{X}-\mathcal{X}^{*}\|_{F}^{2}}{2\tau^{2}n_{1}n_{2}\cdots
n_{d}}\Biggl{\\}}\\\
&\leq\frac{\|\mathcal{X}^{s}-\mathcal{X}^{*}\|_{F}^{2}}{2\tau^{2}n_{1}n_{2}\cdots
n_{d}}\leq\frac{\|\mathcal{X}^{s}-\mathcal{X}^{*}\|_{\infty}^{2}}{2\tau^{2}}\leq\frac{c^{2}}{2\tau^{2}m},\end{split}$
(65)

where $\mathcal{X}^{s}$ is defined in Lemma 8.2. By the construction of
$A_{i}^{s}$ in Lemma 8.2, we know that
$\|A_{i}^{s}\|_{0}=\|A_{i}^{*}\|_{0},\forall i\in[d]$. As a consequence,
plugging (65) and (11) into (64), the estimator of (6) satisfies

$\begin{split}\frac{\mathbb{E}_{\Omega,\mathcal{Y}_{\Omega}}\left[|\mathcal{X}^{\lambda}-\mathcal{X}^{*}\|_{F}^{2}\right]}{n_{1}n_{2}\cdots
n_{d}}\leq&\ \frac{11c^{2}(2\tau+c)^{2}\log(m)}{2\tau^{2}m}\\\ &\
+12\left(1+\frac{2c^{2}}{3\tau^{2}}\right)(2\tau+c)^{2}(\beta+2)\log(n_{m})\frac{r_{1}r_{2}\cdots
r_{d}+\sum_{i=1}^{d}\|A_{i}^{*}\|_{0}}{m},\end{split}$

which is the desired statement. $\Box$

## Appendix D. Proof of Theorem 4.4

Proof. By [7, Lemma 8], we obtain that the KL divergence of Poisson
observations is

$\begin{split}K\left(p_{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*}}(\mathcal{Y}_{i_{1}i_{2}\cdots
i_{d}})||p_{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}}(\mathcal{Y}_{i_{1}i_{2}\cdots
i_{d}})\right)&\leq\frac{1}{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}}(\mathcal{X}_{i_{1}i_{2}\cdots i_{d}}-\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*})^{2}\\\ &\leq\frac{1}{\varrho}(\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}-\mathcal{X}_{i_{1}i_{2}\cdots i_{d}}^{*})^{2}.\end{split}$

Consequently, we can choose $\gamma=\frac{c^{2}}{\varrho}$ in (11). Moreover,
according to the proof of [41, Appendix IV], we have

$\begin{split}-2\log(H(p_{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}}(\mathcal{Y}_{i_{1}i_{2}\cdots
i_{d}}),p_{(\mathcal{X}^{*})_{i_{1}i_{2}\cdots
i_{d}}}(\mathcal{Y}_{i_{1}i_{2}\cdots
i_{d}})))&=\left(\sqrt{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*}}-\sqrt{\mathcal{X}_{i_{1}i_{2}\cdots i_{d}}}\right)^{2}\\\
&\geq\frac{1}{4c}\left(\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}-\mathcal{X}_{i_{1}i_{2}\cdots i_{d}}^{*}\right)^{2},\end{split}$

where the inequality holds by

$\begin{split}(\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}-\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*})^{2}&=\left(\left(\sqrt{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*}}-\sqrt{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}}\right)\left(\sqrt{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*}}+\sqrt{\mathcal{X}_{i_{1}i_{2}\cdots i_{d}}}\right)\right)^{2}\\\
&\leq 4c\left(\sqrt{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}^{*}}-\sqrt{\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}}\right)^{2}.\end{split}$

Hence, by Theorem 4.1, we get that

$\begin{split}&\
\frac{\mathbb{E}_{\Omega,\mathcal{Y}_{\Omega}}\left[\|\mathcal{X}^{\lambda}-\mathcal{X}^{*}\|_{F}^{2}\right]}{n_{1}n_{2}\cdots
n_{d}}\leq\frac{32c^{3}\log(m)}{m\varrho}+12c\cdot\min_{\mathcal{X}\in\Upsilon}\Biggl{\\{}\frac{K(p_{\mathcal{X}^{*}}(\mathcal{Y})||p_{\mathcal{X}}(\mathcal{Y}))}{n_{1}n_{2}\cdots
n_{d}}\\\ &\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\left(\max_{i}\\{\lambda_{i}\\}+\frac{8c^{2}(\beta+2)\log(n_{m})}{3\varrho}\right)\frac{r_{1}r_{2}\cdots
r_{d}+\sum_{i=1}^{d}\|A_{i}\|_{0}}{m}\Biggl{\\}}.\end{split}$ (66)

Let $\mathcal{X}^{s}$ be defined as in Lemma 8.2. Then by a similar argument
of (62), we get that

$\begin{split}\min_{\mathcal{X}\in\Upsilon}\Biggl{\\{}\frac{K(p_{\mathcal{X}^{*}}(\mathcal{Y})||p_{\mathcal{X}}(\mathcal{Y}))}{n_{1}n_{2}\cdots
n_{d}}\Biggl{\\}}&\leq\min_{\mathcal{X}\in\Upsilon}\Biggl{\\{}\frac{\|\mathcal{X}-\mathcal{X}^{*}\|_{F}^{2}}{\varrho
n_{1}n_{2}\cdots n_{d}}\Biggl{\\}}\\\
&\leq\frac{\|\mathcal{X}^{s}-\mathcal{X}^{*}\|_{F}^{2}}{\varrho
n_{1}n_{2}\cdots
n_{d}}\leq\frac{\|\mathcal{X}^{s}-\mathcal{X}^{*}\|_{\infty}^{2}}{\varrho}\leq\frac{c^{2}}{\varrho
m}.\end{split}$ (67)

Notice that $\|A_{i}^{s}\|_{0}=\|A_{i}^{*}\|_{0},\forall i\in[d]$. By
combining (11), (66), and (67), after some rearrangements, we obtain that the
estimator of (6) satisfies

$\begin{split}\frac{\mathbb{E}_{\Omega,\mathcal{Y}_{\Omega}}\left[|\mathcal{X}^{\lambda}-\mathcal{X}^{*}\|_{F}^{2}\right]}{n_{1}n_{2}\cdots
n_{d}}\leq&\ \frac{44c^{3}\log(m)}{\varrho
m}+48c\left(1+\frac{4c^{2}}{3\varrho}\right)(\beta+2)\log(n_{m})\frac{r_{1}r_{2}\cdots
r_{d}+\sum_{i=1}^{d}\|A_{i}^{*}\|_{0}}{m}.\end{split}$

This concludes the proof. $\Box$

## Appendix E. Proof of Theorem 5.1

Let

$\mathfrak{L}=\left\\{\mathcal{X}=\mathcal{C}\times_{1}A_{1}\cdots\times_{d}A_{d}:\mathcal{C}\in\mathfrak{D},A_{i}\in\mathfrak{X}_{i},i\in[d]\right\\},$

where

$\mathfrak{D}=\left\\{\mathcal{C}\in\mathbb{R}_{+}^{r_{1}\times
r_{2}\times\cdots\times r_{d}}:\mathcal{C}_{i_{1}\cdots
i_{d}}\in\\{0,1,c_{0}\\},(i_{1},\ldots,i_{d})\in[n_{1}]\times\cdots\times[n_{d}]\right\\}$

with

$c_{0}=\min\left\\{1,\frac{\gamma_{c}\mu}{\prod_{i=1}^{d}(a_{i}\sqrt{\Delta_{i}(s_{i},n_{i})})}\sqrt{\frac{r_{1}r_{2}\cdots
r_{d}}{m}}\right\\},$ (68)

and

$\mathfrak{X}_{i}=\\{A_{i}\in\mathbb{R}_{+}^{n_{i}\times
r_{i}}:(A_{i})_{jk}\in\\{0,a_{i},b_{i}\\},\|A_{i}\|_{0}\leq
s_{i},(j,k)\in[n_{i}]\times[r_{i}]\\}$

with

$b_{i}=\min\left\\{a_{i},\frac{\gamma_{i}\mu}{\left(\prod_{j\neq
i}a_{j}\right)\sqrt{\prod_{j=1}^{d}\Delta(s_{j},n_{j})}}\sqrt{\frac{s_{i}}{m}}\right\\},\
i\in[d].$

We will specify $\gamma_{c},\gamma_{i},i\in[d]$ in detail later. Therefore, by
the construction of $\mathfrak{L}$, we know that
$\mathfrak{L}\subseteq\mathfrak{L}(\mathbf{s},\mathbf{r},\mathbf{a})$. Next we
consider to construct the packing sets of the core tensor and factor matrices,
respectively.

_Case I_. We construct the following set

$\mathfrak{L}_{\mathcal{C}}=\\{\mathcal{X}=\mathcal{C}\times_{1}A_{1}\cdots\times_{d}A_{d}:\mathcal{C}\in\mathfrak{D}_{0}\\},$

where $\mathfrak{D}_{0}$ is defined as

$\mathfrak{D}_{0}=\left\\{\mathcal{C}\in\mathbb{R}_{+}^{r_{1}\times
r_{2}\times\cdots\times r_{d}}:\mathcal{C}_{i_{1}i_{2}\cdots
i_{d}}\in\\{0,c_{0}\\},(i_{1},\ldots,i_{d})\in[n_{1}]\times\cdots\times[n_{d}]\right\\}$
(69)

and

$A_{i}=a_{i}\begin{pmatrix}I_{r_{i}}\\\ \vdots\\\ I_{r_{i}}\\\
\mathbf{0}_{r_{i}}\end{pmatrix}\in\mathbb{R}_{+}^{n_{i}\times r_{i}}\
\textup{with}\ I_{r_{i}}\in\mathbb{R}^{r_{i}\times r_{i}},\ i\in[d].$

There are $\lfloor\frac{s_{i}\wedge n_{i}}{r_{i}}\rfloor$ blocks identity
matrix $I_{r_{i}}$ in $A_{i}$ and
$\mathbf{0}_{r_{i}}\in\mathbb{R}^{(n_{i}-r_{i}\lfloor\frac{s_{i}\wedge
n_{i}}{r_{i}}\rfloor)\times r_{i}}$ is a zero matrix. Consequently,
$\mathfrak{L}_{\mathcal{C}}\subseteq\mathfrak{L}$.

For any $\mathcal{X}\in\mathfrak{L}_{\mathcal{C}}$, by (2), we have

$\begin{split}\mathcal{X}_{(1)}&=A_{1}\mathcal{C}_{(1)}(A_{d}\otimes
A_{d-1}\otimes\cdots\otimes A_{2})^{T}\\\
&=a_{1}\begin{pmatrix}\mathcal{C}_{(1)}\\\ \vdots\\\ \mathcal{C}_{(1)}\\\
\mathbf{0}_{r_{12}}\end{pmatrix}(A_{d}\otimes A_{d-1}\otimes\cdots\otimes
A_{2})^{T}\end{split}$

where $\mathbf{0}_{r_{12}}\in\mathbb{R}^{(n_{1}-r_{1}\lfloor\frac{s_{1}\wedge
n_{1}}{r_{1}}\rfloor)\times(r_{2}\cdots r_{d})}$ is a zero matrix. By some
simple calculations, we get

$\|\mathcal{X}\|_{F}^{2}=\|\mathcal{X}_{(1)}\|_{F}^{2}=(a_{1}\cdots
a_{d})^{2}\|\mathcal{C}_{(1)}\|_{F}^{2}\prod_{i=1}^{d}\left\lfloor\frac{s_{i}\wedge
n_{i}}{r_{i}}\right\rfloor.$ (70)

Note that the entries of $\mathcal{C}$ only takes $0$ or $c_{0}$ by the
construction of $\mathfrak{D}_{0}$ in (69). By the Varshamov-Gilbert bound
[47, Lemma 2.9], we know that there exists a subset
$\mathfrak{L}_{\mathcal{C}}^{0}\subseteq\mathfrak{L}_{\mathcal{C}}$ such that

$|\mathfrak{L}_{\mathcal{C}}^{0}|\geq 2^{r_{1}\cdots r_{d}/8}+1,$ (71)

and for any
$\mathcal{X}_{1},\mathcal{X}_{2}\in\mathfrak{L}_{\mathcal{C}}^{0}$,

$\begin{split}\|\mathcal{X}_{1}-\mathcal{X}_{2}\|_{F}^{2}&=\|(\mathcal{X}_{1})_{(1)}-(\mathcal{X}_{2})_{(1)}\|_{F}^{2}\\\
&\geq\frac{r_{1}\cdots
r_{d}}{8}\left(\prod_{i=1}^{d}\left\lfloor\frac{s_{i}\wedge
n_{i}}{r_{i}}\right\rfloor\right)(c_{0}a_{1}\cdots a_{d})^{2}\\\
&\geq\frac{n_{1}n_{2}\cdots
n_{d}}{2^{d+3}}\left(\prod_{i=1}^{d}\Delta_{i}(s_{i},n_{i})\right)(c_{0}a_{1}\cdots
a_{d})^{2}\\\ &=\frac{n_{1}n_{2}\cdots
n_{d}}{2^{d+3}}\min\left\\{\prod_{i=1}^{d}a_{i}^{2}\Delta_{i}(s_{i},n_{i}),\gamma_{c}^{2}\mu^{2}\frac{r_{1}r_{2}\cdots
r_{d}}{m}\right\\},\end{split}$

where the second inequality follows from the fact that $\lfloor
x\rfloor\geq\frac{x}{2}$ for any $x\geq 1$, $\Delta_{i}(s_{i},n_{i})$ is
defined as (21), and the last equality follows from the definition of $c_{0}$
in (68).

For an arbitrary tensor $\mathcal{X}\in\mathfrak{L}_{\mathcal{C}}^{0}$, we
have that

$\begin{split}K(\mathbb{P}_{\mathcal{X}},\mathbb{P}_{\mathbf{0}})&\leq\frac{m}{n_{1}n_{2}\cdots
n_{d}}\left(\frac{1}{2\mu^{2}}\right)\sum_{i_{1},i_{2},\ldots,i_{d}}|\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}|^{2}\\\ &\leq\frac{m}{n_{1}n_{2}\cdots
n_{d}}\left(\frac{1}{2\mu^{2}}\right)c_{0}^{2}\left(\prod_{i=1}^{d}r_{i}a_{i}^{2}\left\lfloor\frac{s_{i}\wedge
n_{i}}{r_{i}}\right\rfloor\right)\\\
&\leq\frac{m}{2\mu^{2}}\min\left\\{\prod_{i=1}^{d}a_{i}^{2}\Delta_{i}(s_{i},n_{i}),\gamma_{c}^{2}\mu^{2}\frac{r_{1}r_{2}\cdots
r_{d}}{m}\right\\}\\\ &\leq\frac{\gamma_{c}^{2}r_{1}\cdots r_{d}}{2}\leq
4\gamma_{c}^{2}\log_{2}(|\mathfrak{L}_{\mathcal{C}}^{0}|-1),\end{split}$ (72)

where the first inequality holds by (20) and $\Omega\sim\text{Bern}(p)$ with
$p=\frac{m}{n_{1}n_{2}\cdots n_{d}}$, the second inequality holds by (70) and
$\mathfrak{L}_{\mathcal{C}}^{0}\subseteq\mathfrak{L}_{\mathcal{C}}$, the third
inequality holds by $\lfloor x\rfloor\leq x$ for any $x>0$ and the definition
of $c_{0}$ in (68), and the last inequality holds by (71). Consequently,
summing up the inequality in (72) for all
$\mathcal{X}\in\mathfrak{L}_{\mathcal{C}}^{0}$ except for
$\mathcal{X}=\mathbf{0}$, we deduce

$\frac{1}{|\mathfrak{L}_{\mathcal{C}}^{0}|-1}\sum_{\mathcal{X}\in\mathfrak{L}_{\mathcal{C}}^{0}}K(\mathbb{P}_{\mathcal{X}},\mathbb{P}_{\mathbf{0}})\leq\alpha\log(|\mathfrak{L}_{\mathcal{C}}^{0}|-1)$

where $\gamma_{c}:=\frac{\sqrt{\alpha\log(2)}}{2}$ with $\alpha\in(0,1/8)$.
Hence, by [47, Theorem 2.5], we have

$\begin{split}&\
\inf_{\widetilde{\mathcal{X}}}\sup_{\mathcal{X}^{*}\in\mathfrak{L}(\mathbf{s},\mathbf{r},\mathbf{a})}\mathbb{P}\left(\frac{\|\widetilde{\mathcal{X}}-\mathcal{X}^{*}\|_{F}^{2}}{n_{1}\cdots
n_{d}}\geq\frac{1}{2^{d+4}}\min\left\\{\prod_{i=1}^{d}a_{i}^{2}\Delta_{i}(s_{i},n_{i}),\gamma_{c}^{2}\mu^{2}\frac{r_{1}r_{2}\cdots
r_{d}}{m}\right\\}\right)\\\ \geq&\
\inf_{\widetilde{\mathcal{X}}}\sup_{\mathcal{X}^{*}\in\mathfrak{L}_{\mathcal{C}}^{0}}\mathbb{P}\left(\frac{\|\widetilde{\mathcal{X}}-\mathcal{X}^{*}\|_{F}^{2}}{n_{1}\cdots
n_{d}}\geq\frac{1}{2^{d+4}}\min\left\\{\prod_{i=1}^{d}a_{i}^{2}\Delta_{i}(s_{i},n_{i}),\gamma_{c}^{2}\mu^{2}\frac{r_{1}r_{2}\cdots
r_{d}}{m}\right\\}\right)\\\ \geq&\ \widetilde{\alpha},\end{split}$ (73)

where $\widetilde{\alpha}$ is defined as

$\widetilde{\alpha}:=\frac{\sqrt{|\mathfrak{L}_{\mathcal{C}}^{0}|-1}}{1+\sqrt{|\mathfrak{L}_{\mathcal{C}}^{0}|-1}}\left(1-2\alpha-\sqrt{\frac{2\alpha}{\log(|\mathfrak{L}_{\mathcal{C}}^{0}|-1)}}\right).$

_Case II_. Now we consider the packing set about the factor matrices. For a
fixed $i$, then for any $j\neq i$, we let

$A_{j}=a_{j}\begin{pmatrix}I_{r_{j}}\\\ \vdots\\\ I_{r_{j}}\\\
\mathbf{0}_{j}\end{pmatrix}\in\mathbb{R}_{+}^{n_{j}\times r_{j}},$

where there are $\lfloor\frac{s_{j}\wedge n_{j}}{r_{j}}\rfloor$ blocks
identity matrix $I_{r_{j}}$ in $A_{j}$ and $\mathbf{0}_{j}$ is an
$\big{(}n_{j}-r_{j}\lfloor\frac{s_{j}\wedge n_{j}}{r_{j}}\rfloor\big{)}\times
r_{j}$ zero matrix. Let

$\widetilde{\mathbf{A}}_{i}:=\left\\{A_{i}=(A_{r_{i}^{\prime}}\
\mathbf{0}_{r_{i}^{\prime}})\in\mathbb{R}_{+}^{n_{i}\times
r_{i}}:A_{r_{i}^{\prime}}\in\mathbb{R}_{+}^{n_{i}\times r_{i}^{\prime}}\
\textup{with}\
(A_{r_{i}^{\prime}})_{kj}\in\\{0,b_{i}\\},\mathbf{0}_{r_{i}^{\prime}}\in\mathbb{R}^{n_{i}\times(r_{i}-r_{i}^{\prime})}\right\\},$

where $r_{i}^{\prime}:=\lceil\frac{s_{i}}{n_{i}}\rceil$, there are at most
$s_{i}$ nonzero entries in $A_{r_{i}^{\prime}}$, and
$\mathbf{0}_{r_{i}^{\prime}}$ is an $n_{i}\times(r_{i}-r_{i}^{\prime})$ zero
matrix. Now we construct the following set

$\mathfrak{L}_{A_{i}}=\left\\{\mathcal{X}=\mathcal{C}\times_{1}A_{1}\cdots\times_{i-1}A_{i-1}\times_{i}A_{i}\times_{i+1}A_{i+1}\cdots\times_{d}A_{d}:A_{i}\in\widetilde{\mathbf{A}}_{i}\right\\},$

where the mode-$i$ unfolding of $\mathcal{C}$ is defined as

$\mathcal{C}_{(i)}=\begin{pmatrix}I_{r_{i}^{\prime}}\ \ \ \
I_{r_{i}^{\prime}}\ \ \ \ \cdots\ \ \ \ I_{r_{i}^{\prime}}\ \ \ \ \ \ \ \
\mathbf{0}_{1r_{i}^{\prime}}\\\ \mathbf{0}_{2r_{i}^{\prime}}\ \ \ \
\mathbf{0}_{2r_{i}^{\prime}}\ \ \ \ \cdots\ \ \ \
\mathbf{0}_{2r_{i}^{\prime}}\ \ \ \
\mathbf{0}_{3r_{i}^{\prime}}\end{pmatrix}.$

Here there are $\lfloor\frac{\prod_{j\neq i}r_{j}}{r_{i}^{\prime}}\rfloor$
blocks identity matrices
$I_{r_{i}^{\prime}}\in\mathbb{R}^{r_{i}^{\prime}\times r_{i}^{\prime}}$ in
$\mathcal{C}_{(i)}$, and

$\mathbf{0}_{1r_{i}^{\prime}}\in\mathbb{R}^{r_{i}^{\prime}\times\big{(}\prod_{j\neq
i}r_{j}-r_{i}^{\prime}\lfloor\frac{\prod_{j\neq
i}r_{j}}{r_{i}^{\prime}}\rfloor\big{)}},\
\mathbf{0}_{2r_{i}^{\prime}}\in\mathbb{R}^{(r_{i}-r_{i}^{\prime})\times
r_{i}^{\prime}},\
\mathbf{0}_{3r_{i}^{\prime}}\in\mathbb{R}^{(r_{i}-r_{i}^{\prime})\times\big{(}\prod_{j\neq
i}r_{j}-r_{i}^{\prime}\lfloor\frac{\prod_{j\neq
i}r_{j}}{r_{i}^{\prime}}\rfloor\big{)}}$

are zero matrices. From this construction, we know that
$\mathfrak{L}_{A_{i}}\subseteq\mathfrak{L}$. Therefore, for any
$\mathcal{X}\in\mathfrak{L}_{A_{i}}$, we have

$\begin{split}\mathcal{X}_{(i)}&=A_{i}\mathcal{C}_{(i)}(A_{d}\otimes\cdots\otimes
A_{i+1}\otimes A_{i-1}\otimes\cdots\otimes A_{1})^{T}\\\
&=(A_{r_{i}^{\prime}}\ \cdots\ A_{r_{i}^{\prime}}\
\mathbf{0}_{n_{i}})(A_{d}\otimes\cdots\otimes A_{i+1}\otimes
A_{i-1}\otimes\cdots\otimes A_{1})^{T},\end{split}$

where the first equality follows from (2) and
$\mathbf{0}_{n_{i}}\in\mathbb{R}^{n_{i}\times\big{(}\prod_{j\neq
i}r_{j}-r_{i}^{\prime}\lfloor\frac{\prod_{j\neq
i}r_{j}}{r_{i}^{\prime}}\rfloor\big{)}}$ is a zero matrix. Hence, we can
obtain through some simple calculations that

$\|\mathcal{X}\|_{F}^{2}=\|\mathcal{X}_{(i)}\|_{F}^{2}=\left\lfloor\frac{\prod_{j\neq
i}r_{j}}{r_{i}^{\prime}}\right\rfloor\left(\prod_{j\neq
i}\left(\left\lfloor\frac{s_{j}\wedge n_{j}}{r_{j}}\right\rfloor
a_{j}^{2}\right)\right)\|A_{r_{i}^{\prime}}\|_{F}^{2}.$ (74)

Notice that the entries of $A_{r_{i}^{\prime}}$ only takes $0$ or $b_{i}$. By
the Varshamov-Gilbert bound [47, Lemma 2.9], we know that there exists
$\mathfrak{L}_{A_{i}}^{0}\subseteq\mathfrak{L}_{A_{i}}$ such that

$|\mathfrak{L}_{A_{i}}^{0}|\geq 2^{s_{i}/8}+1,$ (75)

and for any $\mathcal{X}_{1},\mathcal{X}_{2}\in\mathfrak{L}_{A_{i}}^{0}$,

$\begin{split}\|\mathcal{X}_{1}-\mathcal{X}_{2}\|_{F}^{2}&\geq\frac{s_{i}}{8}\left\lfloor\frac{\prod_{j\neq
i}r_{j}}{r_{i}^{\prime}}\right\rfloor\prod_{j\neq
i}\left(\left\lfloor\frac{s_{j}\wedge n_{j}}{r_{j}}\right\rfloor
a_{j}^{2}\right)b_{i}^{2}\\\
&\geq\frac{s_{i}b_{i}^{2}}{2^{d+3}r_{i}^{\prime}}\prod_{j\neq
i}\left({s_{j}\wedge n_{j}}\right)a_{j}^{2}\\\ &\geq\frac{n_{1}\cdots
n_{d}}{2^{d+4}}\left(\prod_{j=1}^{d}\Delta_{j}(s_{j},n_{j})\right)\left(\prod_{j\neq
i}a_{j}^{2}\right)\min\left\\{a_{i}^{2},\frac{\gamma_{i}^{2}\mu^{2}}{(\prod_{j\neq
i}a_{j}^{2})\prod_{j=1}^{d}\Delta_{j}(s_{j},n_{j})}\left(\frac{s_{i}}{m}\right)\right\\}\\\
&=\frac{n_{1}\cdots
n_{d}}{2^{d+4}}\min\left\\{\prod_{j=1}^{d}\Delta_{j}(s_{j},n_{j})a_{j}^{2},\gamma_{i}^{2}\mu^{2}\left(\frac{s_{i}}{m}\right)\right\\},\end{split}$
(76)

where the second inequality holds by $\lfloor x\rfloor\geq\frac{x}{2}$ for any
$x\geq 1$ and the third inequality holds by $x/\lceil
x\rceil\geq\frac{1}{2}\min\\{x,1\\}$ for any $x>0$.

In addition, for any $\mathcal{X}\in\mathfrak{L}_{A_{i}}^{0}$, we have that

$\begin{split}K(\mathbb{P}_{\mathcal{X}},\mathbb{P}_{\mathbf{0}})&\leq\frac{m}{n_{1}n_{2}\cdots
n_{d}}\left(\frac{1}{2\mu^{2}}\right)\sum_{i_{1},i_{2},\ldots,i_{d}}|\mathcal{X}_{i_{1}i_{2}\cdots
i_{d}}|^{2}\\\ &\leq\frac{m}{n_{1}n_{2}\cdots
n_{d}}\left(\frac{1}{2\mu^{2}}\right)\left\lfloor\frac{\prod_{j\neq
i}r_{j}}{r_{i}^{\prime}}\right\rfloor\prod_{j\neq
i}\left(\left\lfloor\frac{s_{j}\wedge n_{j}}{r_{j}}\right\rfloor
a_{j}^{2}\right)(n_{i}r_{i}^{\prime}\wedge s_{i})b_{i}^{2}\\\
&\leq\frac{m}{2\mu^{2}}\left(\prod_{j=1}^{d}\Delta_{j}(s_{j},n_{j})\right)\left(\prod_{j\neq
i}a_{j}\right)b_{i}^{2}\\\
&=\frac{m}{2\mu^{2}}\min\left\\{\prod_{j=1}^{d}a_{j}^{2}\Delta_{j}(s_{j},n_{j}),\gamma_{i}^{2}\mu^{2}\frac{s_{i}}{m}\right\\}\\\
&\leq\frac{\gamma_{i}^{2}s_{i}}{2}\leq
4\gamma_{i}^{2}\log_{2}(|\mathfrak{L}_{A_{i}}^{0}|-1),\end{split}$ (77)

where the first inequality holds by (20) and $\Omega\sim\text{Bern}(p)$ with
$p=\frac{m}{n_{1}n_{2}\cdots n_{d}}$, the second inequality holds by (74) and
$\|A_{r_{i}^{\prime}}\|_{F}^{2}\leq(n_{i}r_{i}^{\prime}\wedge
s_{i})b_{i}^{2}$, the third inequality holds by the fact that $\lfloor
x\rfloor\leq x$ for any $x>0$, and the last inequality holds by (75). Let
$\gamma_{i}={\sqrt{\alpha_{i}\log(2)}}/{2}$ with
$\alpha_{i}\in(0,\frac{1}{8})$. Summing up the inequality in (77) for all
$0\neq\mathcal{X}\in\mathfrak{L}_{A_{i}}$ yields

$\frac{1}{|\mathfrak{L}_{A_{i}}^{0}|-1}\sum_{\mathcal{X}\in\mathfrak{L}_{A_{i}}^{0}}K(\mathbb{P}_{\mathcal{X}},\mathbb{P}_{\mathbf{0}})\leq\alpha_{i}\log(|\mathfrak{L}_{A_{i}}^{0}|-1).$

It follows from [47, Theorem 2.5] that

$\begin{split}&\
\inf_{\widetilde{\mathcal{X}}}\sup_{\mathcal{X}^{*}\in\mathfrak{L}(\mathbf{s},\mathbf{r},\mathbf{a})}\mathbb{P}\left(\frac{\|\widetilde{\mathcal{X}}-\mathcal{X}^{*}\|_{F}^{2}}{n_{1}\cdots
n_{d}}\geq\frac{1}{2^{d+5}}\min\left\\{\prod_{j=1}^{d}\Delta_{j}(s_{j},n_{j})a_{j}^{2},\gamma_{i}^{2}\mu^{2}\left(\frac{s_{i}}{m}\right)\right\\}\right)\\\
&\
\inf_{\widetilde{\mathcal{X}}}\sup_{\mathcal{X}^{*}\in\mathfrak{L}_{A_{i}}^{0}}\mathbb{P}\left(\frac{\|\widetilde{\mathcal{X}}-\mathcal{X}^{*}\|_{F}^{2}}{n_{1}\cdots
n_{d}}\geq\frac{1}{2^{d+5}}\min\left\\{\prod_{j=1}^{d}\Delta_{j}(s_{j},n_{j})a_{j}^{2},\gamma_{i}^{2}\mu^{2}\left(\frac{s_{i}}{m}\right)\right\\}\right)\\\
\geq&\ \widetilde{\alpha}_{i},\end{split}$ (78)

where

$\widetilde{\alpha}_{i}=\frac{\sqrt{|\mathfrak{L}_{A_{i}}^{0}|-1}}{1+\sqrt{|\mathfrak{L}_{A_{i}}^{0}|-1}}\left(1-2\alpha_{i}-\sqrt{\frac{2\alpha_{i}}{\log(|\mathfrak{L}_{A_{i}}^{0}|-1)}}\right).$

Let

$\eta:=\min\left\\{\prod_{j=1}^{d}\Delta_{j}(s_{j},n_{j})a_{j}^{2},\gamma_{m}^{2}\mu^{2}\left(\frac{r_{1}r_{2}\cdots
r_{d}+\sum_{i=1}^{d}s_{i}}{m}\right)\right\\},$

where $\gamma_{m}=\min\\{\gamma_{c},\gamma_{1},\ldots,\gamma_{d}\\}$.
Combining (73) with (78), we deduce that

$\begin{split}&\inf_{\widetilde{\mathcal{X}}}\sup_{\mathcal{X}^{*}\in\mathfrak{L}(\mathbf{s},\mathbf{r},\mathbf{a})}\mathbb{P}\left(\frac{\|\widetilde{\mathcal{X}}-\mathcal{X}^{*}\|_{F}^{2}}{n_{1}\cdots
n_{d}}\geq\frac{1}{2^{d+5}(d+1)}\eta\right)\geq\widetilde{\alpha}_{m},\end{split}$

where
$\widetilde{\alpha}_{m}=\min\\{\widetilde{\alpha},\widetilde{\alpha}_{1},\ldots,\widetilde{\alpha}_{d}\\}\in(0,1).$
Therefore, by Markov’s inequality, we get that

$\begin{split}&\inf_{\widetilde{\mathcal{X}}}\sup_{\mathcal{X}^{*}\in\mathfrak{L}(\mathbf{s},\mathbf{r},\mathbf{a})}\frac{\mathbb{E}_{\Omega,\mathcal{Y}_{\Omega}}\|\widetilde{\mathcal{X}}-\mathcal{X}^{*}\|_{F}^{2}}{n_{1}\cdots
n_{d}}\\\ \geq\
&\frac{\widetilde{\alpha}_{m}}{2^{d+5}(d+1)}\min\left\\{\prod_{i=1}^{d}\Delta_{i}(s_{i},n_{i})a_{i}^{2},\gamma_{m}^{2}\mu^{2}\left(\frac{r_{1}r_{2}\cdots
r_{d}+\sum_{i=1}^{d}s_{i}}{m}\right)\right\\},\end{split}$

which concludes the proof. $\Box$

## References

  * [1] E. Acar, D. M. Dunlavy, T. G. Kolda, and M. Mørup.  Scalable tensor factorizations for incomplete data.  Chemometrics Intell. Lab. Syst., 106(1):41–56, 2011.
  * [2] K. Allab, L. Labiod, and M. Nadif.  A semi-NMF-PCA unified framework for data clustering.  IEEE Trans. Knowl. Data Eng., 29(1):2–16, 2017.
  * [3] X. Bai, F. Xu, L. Zhou, Y. Xing, L. Bai, and J. Zhou.  Nonlocal similarity based nonnegative Tucker decomposition for hyperspectral image denoising.  IEEE J. Select. Topics in Appl. Earth Obs. Remote Sens., 11(3):701–712, 2018.
  * [4] A. Beck.  First-Order Methods in Optimization.  SIAM: Philadelphia, 2017.
  * [5] R. I. Boţ and D.-K. Nguyen.  The proximal alternating direction method of multipliers in the nonconvex setting: convergence analysis and rates.  Math. Oper. Res., 45(2):682–712, 2020.
  * [6] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Found. Trends. Mach. Learn., 3(1):1–122, 2011.
  * [7] Y. Cao and Y. Xie.  Poisson matrix recovery and completion.  IEEE Trans. Signal Process., 64(6):1609–1620, 2016.
  * [8] Y. Chen, T.-Z. Huang, and X.-L. Zhao.  Destriping of multispectral remote sensing image using low-rank tensor decomposition.  IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens., 11(12):4950–4967, 2018.
  * [9] A. Cichocki, D. Mandic, L. De Lathauwer, G. Zhou, Q. Zhao, C. Caiafa, and H. A. Phan.  Tensor decompositions for signal processing applications: From two-way to multiway component analysis.  IEEE Signal Process. Mag., 32(2):145–163, 2015.
  * [10] A. Cichocki, R. Zdunek, A. H. Phan, and S.-i. Amari.  Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation.  U.K..: John Wiley & Sons, 2009.
  * [11] T. M. Cover and J. A. Thomas.  Elements of Information Theory.  2nd Edition, John Wiley & Sons, New Jersey, 2006.
  * [12] L. De Lathauwer, B. De Moor, and J. Vandewalle.  A multilinear singular value decomposition.  SIAM. J. Matrix Anal. Appl., 21(4):1253–1278, 2000.
  * [13] L. De Lathauwer, B. De Moor, and J. Vandewalle.  On the best rank-1 and rank-($r_{1},r_{2},\ldots,r_{n}$) approximation of higher-order tensors.  SIAM J. Matrix Anal. Appl., 21(4):1324–1342, 2000.
  * [14] D. Donoho and V. Stodden.  When does non-negative matrix factorization give a correct decomposition into parts?  In Advances in Neural Information Processing Systems, volume 16, 2003.
  * [15] N. Gillis.  Nonnegative Matrix Factorization.  Philadelphia, SIAM, 2020.
  * [16] N. Gillis and F. Glineur.  Using underapproximations for sparse nonnegative matrix factorization.  Pattern Recognit., 43(4):1676–1687, 2010.
  * [17] N. Gillis and R. Luce.  Robust near-separable nonnegative matrix factorization using linear optimization.  J. Mach. Learn. Res., 15(1):1249–1280, 2014.
  * [18] C. J. Hillar and L.-H. Lim.  Most tensor problems are NP-hard.  J. ACM, 60(6):45, 2013.
  * [19] F. L. Hitchcock.  The expression of a tensor or a polyadic as a sum of products.  J. Math. Physics, 6(1-4):164–189, 1927.
  * [20] D. Hong, T. G. Kolda, and J. A. Duersch.  Generalized canonical polyadic tensor decomposition.  SIAM Rev., 62(1):133–163, 2020.
  * [21] M. Hong, Z.-Q. Luo, and M. Razaviyayn.  Convergence analysis of alternating direction method of multipliers for a family of nonconvex problems.  SIAM J. Optim., 26(1):337–364, 2016.
  * [22] R. A. Horn and C. R. Johnson.  Topics in Matrix analysis.  Cambridge: Cambridge University Press, 1991.
  * [23] P. O. Hoyer.  Nonnegative matrix factorization with sparseness constraints.  J. Mach. Learn. Res., 5:1457–1469, 2004.
  * [24] S. Jain, A. Gutierrez, and J. Haupt.  Noisy tensor completion for tensors with a sparse canonical polyadic factor.  In 2017 IEEE International Symposium on Inf. Theory, pages 2153–2157. IEEE, 2017.
  * [25] M. E. Kilmer and C. D. Martin.  Factorization strategies for third-order tensors.  Linear Algebra Appl., 435(3):641–658, 2011.
  * [26] O. Klopp, K. Lounici, and A. B. Tsybakov.  Robust matrix completion.  Probab. Theory Relat. Fields, 169(1-2):523–564, 2017.
  * [27] T. G. Kolda and B. W. Bader.  Tensor decompositions and applications.  SIAM Rev., 51(3):455–500, 2009.
  * [28] L. G. Kraft.  A device for quantizing, grouping, and coding amplitude-modulated pulses.  Master’s thesis, Massachusetts Institute of Technology, 1949.
  * [29] D. D. Lee and H. S. Seung.  Learning the parts of objects by non-negative matrix factorization.  Nature, 401(6755):788–791, 1999.
  * [30] N. Lee, A.-H. Phan, F. Cong, and A. Cichocki.  Nonnegative tensor train decompositions for multi-domain feature extraction and clustering.  In Int. Conf. Neural Inform. Process., pages 87–95, 2016.
  * [31] Q. J. Li.  Estimation of Mixture Models.  PhD thesis, Yale University, 1999.
  * [32] X. Li, M. K. Ng, G. Cong, Y. Ye, and Q. Wu.  MR-NTD: Manifold regularization nonnegative Tucker decomposition for tensor data dimension reduction and representation.  IEEE Trans. Neural Netw. Learn. Syst., 28(8):1787–1800, 2017.
  * [33] J. Liu, J. Liu, P. Wonka, and J. Ye.  Sparse non-negative tensor factorization using columnwise.  Pattern Recognit., 45(1):649–656, 2012.
  * [34] X. Lu, H. Wu, Y. Yuan, P. Yan, and X. Li.  Manifold regularized sparse NMF for hyperspectral unmixing.  IEEE Trans. Geosci. Remote Sens., 51(5):2815–2826, 2013.
  * [35] B. McMillan.  Two inequalities implied by unique decipherability.  IEEE Trans. Inf. Theory, 2(4):115–116, 1956.
  * [36] M. Mørup, L. K. Hansen, and S. M. Arnfred.  Algorithms for sparse nonnegative Tucker decompositions.  Neural Comput., 20(8):2112–2131, 2008.
  * [37] E. Newman and M. E. Kilmer.  Non-negative tensor patch dictionary approaches for image compression and deblurring applications.  SIAM J. Imaging Sci., 13(3):1084–1112, 2020.
  * [38] I. V. Oseledets.  Tensor-train decomposition.  SIAM J. Sci. Comput., 33(5):2295–2317, 2011.
  * [39] J. Pan and N. Gillis.  Generalized separable nonnegative matrix factorization.  IEEE Trans. Pattern Anal. Mach. Intell., 43(5):1546–1561, 2021\.
  * [40] J. Pan, M. K. Ng, Y. Liu, X. Zhang, and H. Yan.  Orthogonal nonnegative Tucker decomposition.  SIAM J. Sci. Comput., 43(1):B55–B81, 2021.
  * [41] M. Raginsky, R. M. Willett, Z. T. Harmany, and R. F. Marcia.  Compressed sensing performance bounds under Poisson noise.  IEEE Trans. Signal Process., 58(8):3990–4002, 2010.
  * [42] M. Rajapakse, J. Tan, and J. Rajapakse.  Color channel encoding with NMF for face recognition.  In 2004 International Conference on Image Processing, volume 3, pages 2007–2010. IEEE, 2004.
  * [43] A. V. Sambasivan and J. D. Haupt.  Minimax lower bounds for noisy matrix completion under sparse factor models.  IEEE Trans. Inf. Theory, 64(5):3274–3285, 2018.
  * [44] N. D. Sidiropoulos, L. De Lathauwer, X. Fu, K. Huang, E. E. Papalexakis, and C. Faloutsos.  Tensor decomposition for signal processing and machine learning.  IEEE Trans. Signal Process., 65(13):3551–3582, 2017.
  * [45] S. Soltani, M. E. Kilmer, and P. C. Hansen.  A tensor-based dictionary learning approach to tomographic image reconstruction.  BIT Numer. Math., 56(4):1425–1454, 2016.
  * [46] A. Soni, S. Jain, J. Haupt, and S. Gonella.  Noisy matrix completion under sparse factor models.  IEEE Trans. Inf. Theory, 62(6):3636–3661, 2016.
  * [47] A. B. Tsybakov.  Introduction to Nonparametric Estimation.  Springer, New York, 2009.
  * [48] L. R. Tucker.  Some mathematical notes on three-mode factor analysis.  Psychometrika, 31(3):279–311, 1966.
  * [49] N. Vannieuwenhoven, R. Vandebril, and K. Meerbergen.  A new truncation strategy for the higher-order singular value decomposition.  SIAM J. Sci. Comput., 34(2):A1027–A1052, 2012.
  * [50] Y. Wang, W. Yin, and J. Zeng.  Global convergence of ADMM in nonconvex nonsmooth optimization.  J. Sci. Comput., 78(1):29–63, 2019.
  * [51] Y. Xu.  Alternating proximal gradient method for sparse nonnegative Tucker decomposition.  Math. Program.Comput., 7(1):39–70, 2015.
  * [52] Q. Zhang, H. Wang, R. J. Plemmons, and V. P. Pauca.  Tensor methods for hyperspectral data analysis: a space object material identification study.  J. Opt. Soc. Amer. A, 25(12):3001–3012, 2008.
  * [53] X. Zhang and M. K. Ng.  Low rank tensor completion with Poisson observations.  IEEE Trans. Pattern Anal. Mach. Intell., 44(8):4239–4251, 2022\.