# Improving tensor regression by optimal model averaging

Qiushi Bu  Academy of Mathematics and Systems Science, Chinese Academy of
Sciences, Beijing, China  University of Chinese Academy of Sciences, Beijing,
China  Hua Liang  Department of Statistics, George Washington University,
Washington, DC  Xinyu Zhang  Academy of Mathematics and Systems Science,
Chinese Academy of Sciences, Beijing, China  Jiahui Zou  School of Statistics,
Capital University of Economics and Business, Beijing, China

###### Abstract

Tensors have broad applications in neuroimaging, data mining, digital
marketing, etc. CANDECOMP/PARAFAC (CP) tensor decomposition can effectively
reduce the number of parameters to gain dimensionality-reduction and thus
plays a key role in tensor regression. However, in CP decomposition, there is
uncertainty which rank to use. In this article, we develop a model averaging
method to handle this uncertainty by weighting the estimators from candidate
tensor regression models with different ranks. When all candidate models are
misspecified, we prove that the model averaging estimator is asymptotically
optimal. When correct models are included in the candidate models, we prove
the consistency of parameters and the convergence of the model averaging
weight. Simulations and empirical studies illustrate that the proposed method
has superiority over the competition methods and has promising applications.

Keywords: Tensor regression; Model averaging; Model misspecification; Cross-
validation; Asymptotic optimality

##  1 Introduction

In the information age, a huge amount of information can be collected, and the
dimension of the data to be processed is increasing. Tensors, also called
multidimensional arrays, refer to data with higher spatial dimensions, rather
than high-dimensional arrays with many variables in one dimension. For
example, neuroimaging data (Zhou et al.,, 2013; Li et al., 2018b, ),
longitudinal time series data (Wang et al.,, 2021; Si et al.,, 2022), and
digital marketing data (Bi et al.,, 2018) all have 3 or more dimensions. A
simple way to process a tensor is to ignore the spatial structure and
vectorizing it. However, simply vectorizing a tensor and ignoring its spatial
structure may result in a much larger parameter size than the sample size, and
the parameters solved by this method may be inefficient (Yuan and Zhang,,
2016). Taking the magnetic resonance imaging (MRI) data with dimension
$\mathbb{R}^{64\times 64\times 64}$ as an example, the number of free
parameters will be $64\times 64\times 64=262144$ if we use linear regression
after vectorization.

Tensor decomposition solves this problem very well because it can effectively
reduce the size of parameters. Various tensor decomposition methods have been
developed, and two widely used methods are CP (Harshman,, 1970) and Tucker
tensor decomposition (Tucker,, 1966). The CP tensor decomposition can
significantly compress a tensor by representing it as a linear combination of
rank-1 basis tensors. An $N$-way tensor $\mathcal{X}\in\mathbb{R}^{I_{1}\times
I_{2}\times\dots\times I_{N}}$ is rank-1 if it can be presented as the outer
product of $N$ vectors, and the rank of a tensor is the smallest number for
which a tensor can be written as the sum of rank-1 tensors (Kolda and Bader,,
2009). If we apply CP tensor decomposition with a rank-5 approximation to the
same MRI data, the number of free parameters will reduce to $(64+64+64)\times
5=960$, which is much smaller than the value used in the vectorization method.
Guo et al., (2011) first apply the CP tensor decomposition in the linear
regression models, and Zhou et al., (2013) extend it to the generalized linear
models (GLMs). Since then, there has been rapid development in the area of CP
tensor regression. For example, Lock, (2018) studies the case where the
response variable is also a tensor, and Ke et al., (2023) apply the CP tensor
decomposition to quantile regressions.

Algorithms for tensor decomposition all require a given rank, but the rank is
unknown in applications. Håstad, (1990)demonstrates that calculating the rank
is an NP-complete problem. The existing literature mainly uses model selection
methods to choose a rank. For example, Wang et al., (2022) use the Akaike
information criterion (AIC) to choose rank, while Zhou et al., (2013), Zhang
et al., (2019) and Ke et al., (2023) apply the Bayesian information criterion
(BIC) to select the model. Zhang et al., (2019) and Ke et al., (2023) also
prove the rank selected by BIC recovers the true rank for the Gaussian linear
model when the true model is included in the set of candidate models. However,
model selection methods are sometimes unstable since a small perturbation in
the data may cause the selected model to be completely different (Yuan and
Yang,, 2005). Simulations in the current paper also show the poor performance
of model selection methods under high noise levels. Differing from model
selection which chooses the “best” model based on certain criteria, model
averaging weights the estimators of different models and avoids putting all
eggs in one unevenly woven basket (Longford,, 2005).

There are two main streams of model averaging: Bayesian model averaging (BMA)
and frequentist model averaging (FMA). Although BMA is flexible and can work
for a wide range of models (Hoeting et al.,, 1999), it still has drawbacks.
Different prior distributions can have a profound impact on regression
results, but choosing an appropriate prior distribution is often challenging
and experiential. FMA can be further divided into different categories,
including smoothed information criteria (Buckland et al.,, 1997; Hjort and
Claeskens,, 2003), adaptive weighting (Yuan and Yang,, 2005; Yang,, 2001;
Zhang et al.,, 2013), and optimal weighting (Wan et al.,, 2010; Liu,, 2015; Lu
and Su,, 2015). Among these methods, the optimal weighting method is well
developed and widely used. Initially, it focuses on linear models such as
Mallows model averaging (Hansen,, 2007), jackknife model averaging (Hansen and
Racine,, 2012), and heteroskedasticity robust $C_{p}$ (Liu and Okui,, 2013).
In addition to linear regression models, researchers have studied model
averaging under other frameworks such as GLMs (Zhang et al.,, 2016), nonlinear
models (Feng et al.,, 2022) and semi-parametric models (Li et al., 2018a, ).
Model averaging for high-dimensional linear regression is also widely studied
(Ando and Li,, 2014; Yu and Feng,, 2014; Feng and Liu,, 2020). However, few
researchers have focused on model averaging on multidimensional arrays. This
article develops the optimal model averaging method for the CP tensor
regression for the first time.

Based on the block relaxation algorithm proposed by Zhou et al., (2013), we
develop a model averaging method based on Kullback–Leibler ($\mathrm{KL}$)
loss to weight the estimators from candidate models with different ranks under
the GLM framework. $\mathrm{KL}$ divergence is commonly used in evaluating the
performance of the generalized linear model (Ando and Li,, 2017; Zou et al.,,
2022) as a replacement for the square loss. To avoid overfitting, we use a
$J-$fold cross-validation (CV) to calculate model averaging weights. Hansen
and Racine, (2012) and Ando and Li, (2017) apply leave-one-out cross-
validation on linear models and generalized linear models, respectively. Once
the weights are obtained, the model averaging estimator can be represented as
a weighted average of estimators from different models.

One theoretical contribution of this article is to show that our proposed
estimator is asymptotically optimal when all candidate models are
misspecified. Although Zhang et al., (2019) and Ke et al., (2023) prove that
BIC model selection can choose the true rank, it only works when the candidate
model set contains the true model. The asymptotic optimality guarantees that
our method is at least better than BIC model selection in the case of
misspecification. The simulation results also confirm this. Furthermore, the
consistency of the model averaging estimator and the convergence of the weight
are studied when at least one correct model is included in the set of
candidate models. To the best of our knowledge, this is the first time the CP
tensor regression is involved in the FMA framework. This article provides a
new solution to the rank selection problem in tensor regressions and tensor
decompositions.

The remainder of the article is organized as follows. Section 2 begins with
the notation and methodologies of the CP tensor regression, and then develops
the model averaging procedure. Section 3 provides theoretical results of the
model averaging estimator, including asymptotic optimality under the
misspecified framework, root$-n$ consistency and the weight convergence when
the set of candidate models contains the correct model. Section 4 presents
experimental results on simulated data and human brain MRI data to illustrate
the advantage of the proposed method. Further discussions and conclusions are
given in Section 5. Finally, technical proofs are provided in the Appendix.

##  2 Model setup and estimation

###  2.1 Notation and preliminaries

We introduce the notations about tensors, i.e., we use uppercase calligraphic
letters, such as $\mathcal{X}$ and $\mathcal{Y}$, to represent tensors,
uppercase bold letters, such as $\mathbf{X}$ and $\mathbf{Y}$, to represent
matrices, and lowercase bold letters, such as $\mathbf{x}$ and $\mathbf{y}$,
to represent vectors. $\left\|\cdot\right\|$ denotes the $\ell_{2}$ norm of a
vector. The $\mathrm{vec}(\cdot)$ operator stacks the elements of a tensor
$\mathcal{B}$ into a column vector. Specifically, the $(i_{1},\ldots,i_{D})$
entry of $\mathcal{B}\in\mathbb{R}^{p_{1}\times\dots\times p_{D}}$ maps to the
$j$th entry of $\mathrm{vec}(\mathcal{B})$, where
$j=1+\sum_{d=1}^{D}\left(i_{d}-1\right)\prod_{d^{\prime}=1}^{d-1}p_{d^{\prime}}$.
The operator $\circ$ represents the outer product, which means each element of
the outer product is the product of the corresponding vector elements:
$(\boldsymbol{\beta}^{(r)}_{1}\circ\cdots\circ\boldsymbol{\beta}^{(r)}_{D})_{i_{1},\dots,i_{D}}=\beta_{1,i_{1}}^{(r)}\times\dots\times\beta_{D,i_{D}}^{(r)}$.
The inner product of two tensor is written as
$<\mathcal{X},\mathcal{B}>=\mathrm{vec}(\mathcal{X})^{\mathrm{T}}\mathrm{vec}(\mathcal{B})$.
Given a $D-$dimensional rank-$R$ tensor
$\mathcal{B}\in\mathbb{R}^{p_{1}\times\dots\times p_{D}}$, the CP tensor
decomposition factorizes it into a sum of $R$ rank-1 tensors:

$\mathcal{B}=\sum_{r=1}^{R}\boldsymbol{\beta}_{1}^{(r)}\circ\cdots\circ\boldsymbol{\beta}_{D}^{(r)},$

where
$\boldsymbol{\beta}_{d}^{(r)}=(\beta^{(r)}_{d,1},\dots,\beta^{(r)}_{d,p_{d}})^{\mathrm{T}}\in\mathbb{R}^{p_{d}}$
is a column vector for $d=1,\cdots,D$ and $r=1,\cdots,R$. Let
$\mathbf{B}_{d}=[\boldsymbol{\beta}_{d}^{(1)},\dots,\boldsymbol{\beta}_{d}^{(R)}]\in\mathbb{R}^{p_{d}\times
R}$ for $d=1,\dots,D$, then the decomposition can be concisely expressed as
$\mathcal{B}=\llbracket\mathbf{B}_{1},\ldots,\mathbf{B}_{D}\rrbracket$ for a
shorthand (Kolda,, 2006). The operator $\llbracket\cdot\rrbracket$ establishes
a connection between the $D$-dimensional tensor and these $D$ matrices. Figure
1 gives an intuitive explanation of the CP decomposition in the 3-dimensional
case.

Figure 1: The CP decomposition of a 3-dimensional tensor

###  2.2 The CP tensor regression model

Suppose there are $n$ independent observations
$\\{(y_{i},\mathcal{X}_{i}),i=1,\dots,n\\}$, where
$\mathcal{X}_{i}\in\mathbb{R}^{p_{1}\times\cdots\times p_{D}}$ is a
$D$-dimensional tensor. $p_{1},\dots,p_{D}$ are all finite numbers and
independent of $n$. We assume that the response variable $y_{i}$ belongs to an
exponential family:

$\displaystyle\text{Pr}(y_{i}|\theta_{i},\phi)=\mathrm{exp}\left\\{\frac{y_{i}\theta_{i}-b(\theta_{i})}{\phi}+c(y_{i},\phi)\right\\},$
(1)

where $b(\cdot)$ and $c(\cdot,\cdot)$ are known functions, and $\theta_{i}$
and $\phi>0$ denote the natural and dispersion parameters, respectively. Under
the canonical link function, $\theta_{i}$ relates the parameter
$\mathcal{B}\in\mathbb{R}^{p_{1}\times\cdots\times p_{D}}$ to the covariate
$\mathcal{X}_{i}$ in the form

$\displaystyle\theta_{i}=<\mathcal{B},\mathcal{X}_{i}>.$

However, the number of free parameters, $\prod_{i=1}^{D}p_{i}$, is often very
large. Therefore, directly using the maximum likelihood method to estimate the
parameters in (1) will result in significant computational difficulties. If we
know the rank of $\mathcal{B}$ is $R$, we can utilize the CP tensor
decomposition to reduce the parameters, i.e.,

$\displaystyle\mathcal{B}=\sum_{r=1}^{R}\boldsymbol{\beta}_{1}^{(r)}\circ\cdots\circ\boldsymbol{\beta}_{D}^{(r)}=\llbracket\mathbf{B}_{1},\ldots,\mathbf{B}_{D}\rrbracket,$

where $\boldsymbol{\beta}_{d}^{(r)}\in\mathbb{R}^{p_{d}}$ is a column vector
for $d=1,\cdots,D$, $r=1,\cdots,R$, and
$\mathbf{B}_{i}=[\boldsymbol{\beta}_{i}^{(1)},\dots,\boldsymbol{\beta}_{i}^{(R)}]$.
Under the canonical link function, we can estimate the parameters
$\mathcal{B}$ via maximizing the loglikelihood function:

$\begin{split}&\ell\left(\mathbf{B}_{1},\ldots,\mathbf{B}_{D}\right)\\\
&=\sum_{i=1}^{n}\frac{y_{i}\theta_{i}-b\left(\theta_{i}\right)}{\phi}+\sum_{i=1}^{n}c\left(y_{i},\phi\right)\\\
&=\sum_{i=1}^{n}\frac{y_{i}<\mathcal{B},\mathcal{X}_{i}>-b\left(<\mathcal{B},\mathcal{X}_{i}>\right)}{\phi}+\sum_{i=1}^{n}c\left(y_{i},\phi\right).\end{split}$
(2)

Zhou et al., (2013) introduce a block relaxation algorithm to estimate the
parameters, where $\mathbf{B}_{d}$ is sequentially updated for $d=1,\dots,D$.
To address the non-uniqueness issue in the solution, they incorporate some
constraints to control the scale and permutation indeterminacy for ensuring
the uniqueness of estimated parameters.

###  2.3 Model averaging estimation

As mentioned above, a non-negligible problem is that we do not know which $R$
is the best since the algorithm assumes a given rank. Model averaging inspires
us to average possible ranks rather than to just rely on one “good” $R$, which
helps reduce the risk of choosing a “poor” $R$ by mistake. Supposing we have
$S$ candidate models, the rank and the estimated parameter of the $s$th model
are $r_{s}$ and
$\widehat{\mathcal{B}}_{(s)}=\sum_{r=1}^{r_{s}}\widehat{\boldsymbol{\beta}}_{1}^{(r)}\circ\cdots\circ\widehat{\boldsymbol{\beta}}_{D}^{(r)}$,
respectively. $S$ is fixed and does not change with the sample size.
Furthermore, we assume that the rank of the candidate models is sorted in
ascending order, i.e., $r_{1}<r_{2}<\dots<r_{S}$. For notational simplicity,
we add 0 vectors to the low-rank decomposition so that all decompositions have
the same matrix size. Then

$\widehat{\mathcal{B}}_{(s)}=\sum_{r=1}^{r_{s}}\widehat{\boldsymbol{\beta}}_{1}^{(r)}\circ\cdots\circ\widehat{\boldsymbol{\beta}}_{D}^{(r)}+\sum_{r=r_{s}+1}^{r_{S}}\boldsymbol{0}\circ\cdots\circ\boldsymbol{0}=\llbracket\widehat{\mathbf{B}}_{(s),1},\ldots,\widehat{\mathbf{B}}_{(s),D}\rrbracket,$

where
$\widehat{\mathbf{B}}_{(s),d}=[\widehat{\boldsymbol{\beta}}_{d}^{(1)},\dots,\widehat{\boldsymbol{\beta}}_{d}^{(r_{s})},\boldsymbol{0},\dots,\boldsymbol{0}]$
for $d=1,\dots,D$. Let
$\mathbf{w}=(w_{1}\cdots,w_{S})^{\mathrm{T}}\in\mathcal{W}=\\{\mathbf{w}\in[0,1]^{S}:\sum_{s=1}^{S}w_{s}=1\\}$
be the weight vector, and
$\widehat{\mathbf{B}}_{d}(\mathbf{w})=\sum_{s=1}^{S}w_{s}\widehat{\mathbf{B}}_{(s),d}$
for $d=1,\dots,D$. From the perspective of model averaging, we define the
model averaging estimator as the combination of several coefficient arrays
with different ranks:

$\widehat{\mathcal{B}}(\mathbf{w})=\sum_{s=1}^{S}w_{s}\widehat{\mathcal{B}}_{(s)}=\sum_{s=1}^{S}w_{s}\llbracket\widehat{\mathbf{B}}_{(s),1},\ldots,\widehat{\mathbf{B}}_{(s),D}\rrbracket=\llbracket\widehat{\mathbf{B}}_{1}(\mathbf{w}),\ldots,\widehat{\mathbf{B}}_{D}(\mathbf{w})\rrbracket.$

Let $\mathcal{B}_{0}$ be the true parameter for $\mathcal{B}$,
$\mathbf{y}=(y_{1},\dots,y_{n})^{\mathrm{T}}$,
$\boldsymbol{\mu}=\mathbb{E}(\mathbf{y})$,
$\theta_{0i}=<\mathcal{B}_{0},\mathcal{X}_{i}>$ for $i=1,\dots,n$,
$B_{0}=\sum_{i=1}^{n}b(\theta_{0i})$, and
$\boldsymbol{\theta}_{0}=(\theta_{01},\cdots,\theta_{0n})^{\mathrm{T}}$.
Furthermore, let
$B\\{\widehat{\mathcal{B}}(\mathbf{w})\\}=\sum_{i=1}^{n}b[\theta_{i}\\{\widehat{\mathcal{B}}(\mathbf{w})\\}]$,
and
$\boldsymbol{\theta}\\{\widehat{\mathcal{B}}(\mathbf{w})\\}=(\theta_{1}\\{\widehat{\mathcal{B}}(\mathbf{w})\\},\cdots,\theta_{n}\\{\widehat{\mathcal{B}}(\mathbf{w})\\})^{\mathrm{T}}$.
After obtaining the estimators
${\widehat{\mathcal{B}}}_{(1)},\dots,{\widehat{\mathcal{B}}}_{(S)}$, we get
the $\mathrm{KL}$ loss of
$\boldsymbol{\theta}\\{\widehat{\mathcal{B}}(\mathbf{w})\\}$:

$\begin{split}\mathrm{KL}(\mathbf{w})&=2\mathbb{E}_{\mathbf{y}^{*}}\log\frac{\prod_{i=1}^{n}\text{Pr}(y_{i}^{*}|\theta_{0i},\phi)}{\prod_{i=1}^{n}\text{Pr}[y_{i}^{*}|\theta_{i}\\{\widehat{\mathcal{B}}(\mathbf{w})\\},\phi]}\\\
&=2\phi^{-1}B\\{\widehat{\mathcal{B}}(\mathbf{w})\\}-2\phi^{-1}\boldsymbol{\mu}^{\mathrm{T}}\boldsymbol{\theta}\\{\widehat{\mathcal{B}}(\mathbf{w})\\}-2\phi^{-1}B_{0}+2\phi^{-1}\boldsymbol{\mu}^{\mathrm{T}}\boldsymbol{\theta}_{0},\end{split}$
(3)

where $y_{i}^{*}$ is another realization from
$\text{Pr}(\cdot|\theta_{i},\phi)$ and independent of $y_{i}$. Note that the
last two terms are not related to $\mathbf{w}$, so we just need to minimize
$2\phi^{-1}B\\{\widehat{\mathcal{B}}(\mathbf{w})\\}-2\phi^{-1}\boldsymbol{\mu}^{\mathrm{T}}\boldsymbol{\theta}\\{\widehat{\mathcal{B}}(\mathbf{w})\\}$
to obtain the optimal weights. Considering $\boldsymbol{\mu}$ is unknown, we
tend to use $\mathbf{y}$ to replace $\boldsymbol{\mu}$, but this may lead to
overfitting. Hence, we utilize the $J-$fold cross-validation to relieve this
trouble and set the criterion to calculate model-averaging weights. For
simplicity of expression, we assume that $n/J$ is an integer and $n_{0}=n/J$.
We introduce additional notations for the weight selection criterion after
introducing cross-validation. Let
$\theta_{j,i}\\{\mathcal{B}\\}=\theta_{(j-1)n_{0}+i}\\{\mathcal{B}\\}$,
$\mathcal{X}_{j,i}=\mathcal{X}_{(j-1)n_{0}+i}$, $y_{j,i}=y_{(j-1)n_{0}+i}$,
and $\mathbf{y}^{(j)}=(y_{j,1},\dots,y_{j,n_{0}})^{\mathrm{T}}$ be the
elements in the $j$th fold, and $\widehat{\mathcal{B}}^{[-j]}_{(s)}$ be the
estimator of $\mathcal{B}$ of the $s$th model without the $j$th fold. Then we
have
$\widehat{\mathcal{B}}^{[-j]}(\mathbf{w})=\sum_{s=1}^{S}w_{s}\widehat{\mathcal{B}}^{[-j]}_{(s)}$,
$\boldsymbol{\theta}^{[-j]}\\{\widehat{\mathcal{B}}^{[-j]}(\mathbf{w})\\}=(\theta_{j,1}\\{\widehat{\mathcal{B}}^{[-j]}(\mathbf{w})\\},\cdots,\theta_{j,n_{0}}\\{\widehat{\mathcal{B}}^{[-j]}(\mathbf{w})\\})^{\mathrm{T}}$
and
$B^{[-j]}\\{\widehat{\mathcal{B}}^{[-j]}(\mathbf{w})\\}=\sum_{i=1}^{n_{0}}b[\theta_{j,i}\\{\widehat{\mathcal{B}}^{[-j]}(\mathbf{w})\\}]$.
$B\\{\widehat{\mathcal{B}}(\mathbf{w})\\}$ and
$\boldsymbol{\theta}\\{\widehat{\mathcal{B}}(\mathbf{w})\\}$ in (3) change
into

$\widetilde{B}\\{\widetilde{\mathcal{B}}(\mathbf{w})\\}=\sum_{j=1}^{J}B^{[-j]}\\{\widehat{\mathcal{B}}^{[-j]}(\mathbf{w})\\}=\sum_{j=1}^{J}\sum_{i=1}^{n_{0}}b[\theta_{j,i}\\{\widehat{\mathcal{B}}^{[-j]}(\mathbf{w})\\}],$

and

$\widetilde{\boldsymbol{\theta}}\\{\widetilde{\mathcal{B}}(\mathbf{w})\\}=(\boldsymbol{\theta}^{[-1]}\\{\widehat{\mathcal{B}}^{[-1]}(\mathbf{w})\\}^{\mathrm{T}},\dots,\boldsymbol{\theta}^{[-J]}\\{\widehat{\mathcal{B}}^{[-J]}(\mathbf{w})\\}^{\mathrm{T}})^{\mathrm{T}}.$

Then the weight selection criterion $\mathrm{CV}_{J}(\mathbf{w})$ is written
as

$\mathrm{CV}_{J}(\mathbf{w})=2\phi^{-1}\widetilde{B}\\{\widetilde{\mathcal{B}}(\mathbf{w})\\}-2\phi^{-1}\mathbf{y}^{\mathrm{T}}\widetilde{\boldsymbol{\theta}}\\{\widetilde{\mathcal{B}}(\mathbf{w})\\}.$

In general, the selection of $J$ is 5 or 10 (Ando and Li,, 2014; Zhang and
Liu,, 2023; Gao et al.,, 2023). Arlot and Lerasle, (2016) also argue that
choosing $J$ to be 5 or 10 is a good choice. In this article, we assume
$J=O(1)$ in the asymptotic theory part, and $J=5$ in the simulations and real
data experiments. Finally, the resultant weight vector is defined as

$\displaystyle\widehat{\mathbf{w}}=\operatorname{argmin}_{\mathbf{w}\in\mathcal{W}}\mathrm{CV}_{J}(\mathbf{w}).$
(4)

Once the estimated weights are available, we can calculate the model averaging
estimator, denoted by
$\widehat{\mathcal{B}}(\widehat{\mathbf{w}})=\sum_{s=1}^{S}\widehat{w}_{s}\widehat{\mathcal{B}}_{(s)}$.
We term the proposed method as Tensor Regression Model Averaging (TRMA).

##  3 Asymptotic theory

###  3.1 Asymptotic optimality

It should be emphasized that the conclusions in this paper are based on the
assumption that the number of parameters $\prod_{d=1}^{D}p_{d}$ and the number
of models $S$ are fixed, and the sample size $n$ tends to infinity. Before
presenting the theorems, we first specify some regularity conditions. Let
$\boldsymbol{\beta}_{\mathcal{B}}=\mathrm{vec}([{\mathbf{B}}_{1},\dots,{\mathbf{B}}_{D}])$,
$\epsilon_{i}=y_{i}-\mu_{i}$ for $i=1,\dots,n$ and
$\boldsymbol{\epsilon}=(\epsilon_{1},\dots,\epsilon_{n})^{\mathrm{T}}=\mathbf{y}-\boldsymbol{\mu}$.

######  Condition 1.

For any $s=1,\dots,S$, there exists
$\mathcal{B}^{*}_{(s)}=\llbracket\mathbf{B}^{*}_{(s),1},\dots,\mathbf{B}^{*}_{(s),D}\rrbracket$
belonging to a compact space $\mathbb{B}$, such that

$\left\|\boldsymbol{\beta}_{\widehat{\mathcal{B}}_{(s)}}-\boldsymbol{\beta}_{\mathcal{B}^{*}_{(s)}}\right\|=O_{p}(n^{-1/2}).$
(5)

Besides, $\mathcal{B}^{*}_{(s)}$ and the true value $\mathcal{B}_{0}$ for the
unknown parameter lie in the interior of the compact space $\mathbb{B}$.

######  Remark 1.

$\mathcal{B}^{*}_{(s)}$ is the quasi-true parameter that minimizes
$\mathrm{KL}$ loss between the $s$th candidate model and the true model. The
quasi-true parameter equals the true parameter if the working model is
correctly specified. Remark 5 in Zhou et al., (2013) guarantees that the
maximum likelihood estimator is the best rank-$R$ approximation of
$\mathcal{B}_{true}$ in the sense of Kullback–Leibler distance. Given
Condition 1, it is easy to verify that

$\left\|\mathrm{vec}(\widehat{\mathcal{B}}_{(s)})-\mathrm{vec}(\mathcal{B}^{*}_{(s)})\right\|=O_{p}(n^{-1/2}).$
(6)

The difference between (5) and (6) is that (5) requires the consistency of the
vectors after the CP decomposition, while (6) indicates the consistency of the
whole tensor. Equation (5) is stronger than (6).

######  Condition 2.

There is a constant $c_{0}$ such that

$\sup_{\mathcal{B}\in\mathbb{B}}\max_{1\leq i\leq
n}\left\|\frac{\partial\theta_{i}\\{\mathcal{B}\\}}{\partial\boldsymbol{\beta}_{\mathcal{B}}}\right\|\leq
c_{0}<\infty.$

######  Remark 2.

Condition 2 states that the gradients of the systematic part are well defined,
which is similar to Condition (C1) in Zhang et al., (2016).

######  Condition 3.

There exist positive constants $c_{1}$ and $c_{2}$ such that

$\frac{\left\|\boldsymbol{\mu}\right\|^{2}}{n}\leq c_{1}<\infty$

and

$\max_{1\leq i\leq n}\left\|\mathrm{vec}(\mathcal{X}_{i})\right\|\leq
c_{2}<\infty.$

######  Condition 4.

As $n\rightarrow\infty$, there exist positive constants $c_{\min}$ and
$c_{\max}$ independent of n such that

$0<c_{\min}\leq\lambda_{\min}\left\\{\frac{1}{n}\sum_{i=1}^{n}\mathrm{vec}(\mathcal{X}_{i})\mathrm{vec}(\mathcal{X}_{i})^{\mathrm{T}}\right\\},$

and

$\lambda_{\max}\left\\{\frac{1}{n}\sum_{i=1}^{n}\mathrm{vec}(\mathcal{X}_{i})\mathrm{vec}(\mathcal{X}_{i})^{\mathrm{T}}\right\\}\leq
c_{\max}<\infty,$

where $\lambda_{\min}(\cdot)$ and $\lambda_{\max}(\cdot)$ represent the
minimum and maximum eigenvalue of a matrix, respectively.

######  Remark 3.

Condition 3 places restrictions on the variability of covariates, which is
similar to Condition 8 in Ando and Li, (2014) and Conditions 2-3 in Zou et
al., (2022). Condition 4 requires that the maximum and minimum eigenvalues of
the matrix
$1/n\sum_{i=1}^{n}\mathrm{vec}(\mathcal{X}_{i})\mathrm{vec}(\mathcal{X}_{i})^{\mathrm{T}}$
are bounded, which is similar to Condition C.1 in Zhao et al., (2020).
Condition 4 also implies that the matrix
$1/n\sum_{i=1}^{n}\mathrm{vec}(\mathcal{X}_{i})\mathrm{vec}(\mathcal{X}_{i})^{\mathrm{T}}$
is of full rank, which is a necessary condition to ensure the identifiability
of the tensor regression.

######  Condition 5.

(i) There exists a positive constant $c_{3}$ such that for $k=1,2$,

$\sup_{\mathcal{B}\in\mathbb{B}}\max_{1\leq i\leq
n}\left|b^{(k)}\left(\theta_{i}\\{\mathcal{B}\\}\right)\right|\leq
c_{3}<\infty;$

(ii) there is a positive constant $c_{4}$ such that

$\inf_{\mathcal{B}\in\mathbb{B}}\min_{1\leq i\leq
n}b^{{}^{\prime\prime}}\left(\theta_{i}\\{\mathcal{B}\\}\right)\geq c_{4}>0;$

######  Remark 4.

Condition 5 regulates the behavior of function $b(\cdot)$. Common
distributions such as the Poisson distribution, binary distribution and normal
distribution all satisfy this condition. This condition is similar to
Condition R2 in Ando and Li, (2017).

Let

$\mathrm{KL}^{*}(\mathbf{w})=2\phi^{-1}B\\{{\mathcal{B}}^{*}(\mathbf{w})\\}-2\phi^{-1}\boldsymbol{\mu}^{\mathrm{T}}\boldsymbol{\theta}\\{{\mathcal{B}}^{*}(\mathbf{w})\\}-2\phi^{-1}B_{0}+2\phi^{-1}\boldsymbol{\mu}^{\mathrm{T}}\boldsymbol{\theta}_{0},$

and $\xi_{n}=\inf_{\mathbf{w}\in\mathcal{W}}\mathrm{KL}^{*}(\mathbf{w})$.

######  Condition 6.

$n\xi_{n}^{-2}=o(1)$.

######  Remark 5.

Condition 6 requires $\xi_{n}$ to grow faster than $n^{1/2}$, which means the
candidate models are not too close to the true model. It rules out the case
where the true model exists in the set of candidate models. In particular, if
the ${s^{*}}$th model is the correct model and included in the set of
candidate models, we have
$\xi_{n}\leq\mathrm{KL}^{*}(\mathbf{w}_{s^{*}})=2\phi^{-1}(B\\{{\mathcal{B}}^{*}_{(s^{*})}\\}-B_{0})-2\phi^{-1}\boldsymbol{\mu}^{\mathrm{T}}(\boldsymbol{\theta}\\{{\mathcal{B}}^{*}_{(s^{*})}\\}-\boldsymbol{\theta}_{0})=0$,
where $\mathbf{w}_{s^{*}}$ is a weight vector in which the $s^{*}$ element is
1 and other elements are 0. In this case, $\xi_{n}=0$ and $n\xi_{n}^{-2}$ will
not converge to 0. This is also a common condition in model averaging
articles, for example in Ando and Li, (2014), Liu et al., (2020) and Zhang et
al., (2023).

######  Theorem 1.

If Conditions 1-6 hold, we have

$\displaystyle\frac{\mathrm{KL}(\widehat{\mathbf{w}})}{\inf_{\mathbf{w}\in\mathcal{W}}\mathrm{KL}(\mathbf{w})}\to
1$ (7)

in probability as $n\to\infty$.

Theorem 1 shows that the TRMA estimator given by (4) is asymptotically optimal
in the sense of $\mathrm{KL}$ divergence when all candidate models are
misspecified. In other words, the weight $\widehat{\mathbf{w}}$ yields a KL
loss that is asymptotically equivalent to that of the infeasible optimal
weight. Model misspecification is the common case because real tensor data
rarely have low-rank decompositions. Theorem 1 guarantees that our proposed
method will not be worse than other methods as $n\to\infty$ under model
misspecification.

###  3.2 Consistency

Following Zhang and Liu, (2019), the true model refers to the just-fitted
model. Models $1,\dots,s_{0}$ are termed as the underfitted models and models
$s_{0}+2,\dots,S$ correspond to the overfitted models. Both just-fitted models
and overfitted models are called correct models.

Zhou et al., (2013) prove that the estimated parameter of the true model
$\boldsymbol{\beta}_{\widehat{\mathcal{B}}}$ converges to the true parameter
$\boldsymbol{\beta}_{\mathcal{B}_{0}}$. In this article, we do not require the
decomposed parameter to converge to $\boldsymbol{\beta}_{\mathcal{B}_{0}}$,
but are more concerned with the consistency of the original tensor. That is,
we focus on whether the model averaging estimator converges to
$\mathcal{B}_{0}$. The following theorem gives the answer.

######  Theorem 2.

If there exist correct models in the set of candidate models and Conditions
1-5 hold, then

$\displaystyle\left\|\mathrm{vec}(\widehat{\mathcal{B}}(\widehat{\mathbf{w}}))-\mathrm{vec}(\mathcal{B}_{0})\right\|=O_{p}(n^{-1/2}).$
(8)

Theorem 2 states that the TRMA estimator converges to the true parameter at
the rate of $\sqrt{n}$ when the correct models are contained in the candidate
model set. This theorem can be viewed as a supplement of Theorem 1, which
assumes all candidate models are misspecified.

###  3.3 Weights of misclassified models

Now we focus on the properties of weights when the correct models exist. Let

$\Delta_{j}=\left\\{\left(\mathrm{vec}(\widehat{\mathcal{B}}^{[-j]}_{(1)})-\mathrm{vec}(\widehat{\mathcal{B}}^{[-j]}_{(s_{0}+1)})\right),\ldots,\left(\mathrm{vec}(\widehat{\mathcal{B}}^{[-j]}_{(s_{0})})-\mathrm{vec}(\widehat{\mathcal{B}}^{[-j]}_{(s_{0}+1)})\right)\right\\}^{\mathrm{T}},$

and $\widehat{w}_{\Delta}=\sum_{s=1}^{s_{0}}\widehat{w}_{s}$ be the sum of the
optimal weights assigned to underfitted models.

######  Condition 7.

For sufficiently large $n$, there exists a positive constant $c_{0}$ such that

$\min_{1\leq j\leq
J}\lambda_{\min}\left({\Delta_{j}}{\Delta_{j}}^{\mathrm{T}}\right)\geq c_{0}$

holds almost surely.

######  Remark 6.

Condition 7 is similar to Assumption 5 in Liu and Zhang, (2022). This
condition requires that there is a certain gap between the parameters of the
underfitted models and those of the true model.

######  Theorem 3.

If Conditions 1-5 and 7 hold and there exist correct models in the set of
candidate models, we have

$\widehat{w}_{\Delta}=o_{p}(1).$

Theorem 3 means that the weights assigned to the underfitted models converge
to zero as $n\to\infty$. In other words, the TRMA estimator ignores the
“wrong” models and is determined only by the correct models. The three
theorems take into account both misspecified and correctly specified
frameworks, and in each case, our proposed method demonstrates a certain
degree of superiority.

##  4 Simulation Study

###  4.1 2-D simulation

Similar to Zhou et al., (2013), we first conduct two-dimensional shape
examples. We compare the performances of different methods under a variety of
signal shapes, sample sizes and noise levels. The response $y_{i}$ is
generated from

$y_{i}=\eta_{i}+\epsilon_{i}=<\mathcal{X}_{i},\mathcal{B}>+\epsilon_{i},$

where $\mathcal{X}_{i}$ is a $64\times 64$ random matrix with all elements
being independent standard normal distribution. The noise term
$\epsilon_{i}\sim N(0,\sigma^{2})$, and we set the noise level to be
$5\%,10\%$, and $25\%$ of the standard deviation of mean $\eta_{i}$ for
different noise levels. $\mathcal{B}$ is a $64\times 64$ matrix, with the
signal region being 1 and the rest 0. Figure 2 illustrates the six signal
shapes, where the black area is 1 and the rest is 0. For sample sizes, we set
$n$=500, 750 and 1000. There are five candidate models with the CP tensor
regression from rank-1 to rank-5, respectively. To evaluate the performance of
TRMA, the following methods are used as competitors:

  * •

AIC model selection (Akaike,, 1973).

  * •

BIC model selection (Schwarz,, 1978).

  * •

$\mathbf{w}=(w_{1},\dots,w_{S})^{\mathrm{T}}$, where
$w_{s}=\mathrm{exp}(-\text{AIC}_{s})/\sum_{i=1}^{S}\mathrm{exp}(-\text{AIC}_{i})$,
proposed by Buckland et al., (1997). We call it “SAIC”.

  * •

$\mathbf{w}=(w_{1},\dots,w_{S})^{\mathrm{T}}$, where
$w_{s}=\mathrm{exp}(-\text{BIC}_{s})/\sum_{i=1}^{S}\mathrm{exp}(-\text{BIC}_{i})$,
which is also proposed by Buckland et al., (1997). We call it “SBIC”.

  * •

$\mathbf{w}=(0,\dots,0,1)^{\mathrm{T}}$. This is the tensor regression model
corresponding to the maximal rank decomposition. We call it “MAX”.

  * •

$\mathbf{w}=(1/S,\dots,1/S)^{\mathrm{T}}$, where S is the number of models.
This is the simple average of all models and we call it “EQMA”.

  * •

Vectorize $\mathcal{X}$ and use LASSO to get the parameters. This model is
used to illustrate the necessity of tensor decomposition on tensor data. We
use 5-fold cross-validation to search for the best penalty parameter
$\lambda$. We call it “LASSO”.

  * •

5-fold TRMA, the proposed method. We call it “TRMA”.

We evaluate six different signals, where the first three signals have rank-1,
rank-2, and rank-3 CP decomposition and the rest have no low-rank CP
decomposition. Different experimental settings allow us to evaluate the
performances of different methods when the set of candidate models is
misspecified or correctly specified.

(a) Signal 1

(b) Signal 2

(c) Signal 3

(d) Signal 4

(e) Signal 5

(f) Signal 6

Figure 2: Six different signals, where the black area is 1 and the rest is 0

Figures 3(a)-3(c) show the restoration of six signals by different methods
when the noise level is 5% of the standard deviation of mean $\eta_{i}$, with
sample sizes of 500, 750, and 1000. For the LASSO method, the number of
parameters is $64\times 64=4096$, which is much larger than the sample size.
In comparison, the rank-5 CP tensor regression method has only $64\times
2\times 5-5^{2}+5=620$ parameters. The last two terms “$-5^{2}+5$” are
restrictions imposed to avoid indeterminacy. As a result, in Figure 3(a) we
can see that the LASSO method cannot recover the signal very well due to the
small sample size. Therefore, using tensor decomposition methods to reduce the
number of parameters is of great significance when the number of parameters
exceeds the sample size. In Figure 3(a), the recovery results of the TRMA and
EQMA methods show the contour of the original image on Signals 3-6, while the
results of other methods are mostly chaotic. This demonstrates the robustness
of our proposed method when the sample size is relatively small. As the sample
size increases from 500 to 1000, the AIC, BIC, SAIC, and SBIC methods show
clearer signal recovery for Signals 3-6. However, our method still outperforms
these methods by producing more distinct black and white areas with higher
contrast.

(a) $n=500$

(b) $n=750$

(c) $n=1000$

Figure 3: Recovered signal images by different methods for (a) $n=500$, (b)
$750$, and (c) $1000$ and the noise level is 5% of the standard deviation of
mean $\eta_{i}$.

The root mean squared errors (RMSEs) of $\mathcal{B}$ obtained from different
methods are listed below to quantitatively compare the performance of
different methods. Each experimental setup is replicated 100 times, and the
means and standard deviations of RMSE are shown in Table 1. Most of the time,
the TRMA method yields the smallest RMSE, with the second smallest result in
the remaining cases. The LASSO method does not perform well due to the large
number of parameters. When there exist correct models in the set of candidate
models, i.e., in Signals 1-3, BIC and TRMA have the smallest RMSE. When the
signals do not have low-rank decompositions, such as in Signals 4-6, AIC and
BIC tend to choose the maximal rank model, resulting in the same RMSE as MAX.
However, we can see our TRMA method outperforms even the maximal rank model.
Furthermore, as the noise level increases, AIC and BIC perform worse and
cannot find the true model in Signals 1-3. Under these circumstances, TRMA has
a smaller RMSE, which highlights the instability of model selection under high
noise levels. Although EQMA sometimes has the best performance when the
signals do not have low-rank decompositions, it is not robust because it has a
larger RMSE when the correct model is included in the set of candidate models.
EQMA simply averages the parameters of all models. So when the true model
exists, the misspecified models will greatly affect the results.

Table 1: RMSEs of $\widehat{\mathcal{B}}$ when $n=500,750,1000$ and the noise
level is 5%, 10% and 20% of the standard deviation of mean $\eta_{i}$.

$n=500$ and the noise level is 5% of the standard deviation of mean
$\eta_{i}$.
---
| Signal 1 | Signal 2 | Signal 3 | Signal 4 | Signal 5 | Signal 6
| Mean | Std | Mean | Std | Mean | Std | Mean | Std | Mean | Std | Mean | Std
AIC | $0.0094^{(2)}$ | 0.0022 | $0.0159^{(1)}$ | 0.0009 | 0.5596 | 0.1838 | 0.6102 | 0.0825 | 0.8796 | 0.0770 | 0.9078 | 0.0863
BIC | $0.0091^{(1)}$ | 0.0007 | $0.0159^{(1)}$ | 0.0009 | 0.5488 | 0.1957 | 0.6102 | 0.0825 | 0.8796 | 0.0770 | 0.9103 | 0.0869
SAIC | $0.0094^{(2)}$ | 0.0019 | $0.0159^{(1)}$ | 0.0009 | 0.5596 | 0.1838 | 0.6102 | 0.0825 | 0.8796 | 0.0770 | 0.9078 | 0.0863
SBIC | $0.0091^{(1)}$ | 0.0007 | $0.0159^{(1)}$ | 0.0009 | 0.5488 | 0.1957 | 0.6102 | 0.0825 | 0.8796 | 0.0770 | 0.9103 | 0.0869
MAX | 0.6083 | 0.1000 | 0.7558 | 0.0670 | 0.6392 | 0.0458 | 0.7317 | 0.0774 | 0.9057 | 0.0580 | 0.9425 | 0.0732
EQMA | 0.1313 | 0.0218 | 0.2267 | 0.0181 | $0.2502^{(2)}$ | 0.0170 | $0.2474^{(2)}$ | 0.0139 | $0.3893^{(2)}$ | 0.0159 | $0.3876^{(2)}$ | 0.0183
LASSO | 0.3008 | 0.0030 | 0.3545 | 0.0019 | 0.2862 | 0.0030 | 0.3457 | 0.0019 | 0.4083 | 0.0014 | 0.4304 | 0.0016
TRMA | $0.0091^{(1)}$ | 0.0007 | $0.0260^{(2)}$ | 0.0246 | $0.1975^{(1)}$ | 0.0177 | $0.1602^{(1)}$ | 0.0059 | $0.3094^{(1)}$ | 0.0079 | $0.2986^{(1)}$ | 0.0072
$n=750$ and the noise level is 5% of the standard deviation of mean
$\eta_{i}$.
| Signal 1 | Signal 2 | Signal 3 | Signal 4 | Signal 5 | Signal 6
| Mean | Std | Mean | Std | Mean | Std | Mean | Std | Mean | Std | Mean | Std
AIC | $0.0077^{(2)}$ | 0.0023 | $0.0119^{(2)}$ | 0.0027 | $0.0162^{(2)}$ | 0.0029 | 0.2915 | 0.0238 | 0.5895 | 0.0510 | 0.5629 | 0.0446
BIC | $0.0070^{(1)}$ | 0.0005 | $0.0112^{(1)}$ | 0.0006 | $0.0156^{(1)}$ | 0.0008 | 0.2915 | 0.0238 | 0.5895 | 0.0510 | 0.5629 | 0.0446
SAIC | $0.0077^{(2)}$ | 0.0023 | $0.0118^{(2)}$ | 0.0025 | $0.0162^{(2)}$ | 0.0029 | 0.2915 | 0.0238 | 0.5895 | 0.0510 | 0.5629 | 0.0446
SBIC | $0.0070^{(1)}$ | 0.0005 | $0.0112^{(1)}$ | 0.0006 | $0.0156^{(1)}$ | 0.0008 | 0.2915 | 0.0238 | 0.5895 | 0.0510 | 0.5629 | 0.0446
MAX | 0.0401 | 0.0026 | 0.0421 | 0.0030 | 0.0482 | 0.0307 | 0.2915 | 0.0238 | 0.5895 | 0.0510 | 0.5629 | 0.0446
EQMA | 0.0137 | 0.0006 | 0.0411 | 0.0010 | 0.0659 | 0.0042 | $0.1276^{(2)}$ | 0.0048 | $0.2621^{(2)}$ | 0.0083 | $0.2457^{(2)}$ | 0.0074
LASSO | 0.2889 | 0.0062 | 0.3463 | 0.0050 | 0.2718 | 0.0062 | 0.3356 | 0.0047 | 0.4023 | 0.0042 | 0.4250 | 0.0043
TRMA | $0.0070^{(1)}$ | 0.0005 | $0.0112^{(1)}$ | 0.0006 | $0.0156^{(1)}$ | 0.0008 | $0.1208^{(1)}$ | 0.0029 | $0.2465^{(1)}$ | 0.0059 | $0.2334^{(1)}$ | 0.0060
$n=1000$ and the noise level is 5% of the standard deviation of mean
$\eta_{i}$.
| Signal 1 | Signal 2 | Signal 3 | Signal 4 | Signal 5 | Signal 6
| Mean | Std | Mean | Std | Mean | Std | Mean | Std | Mean | Std | Mean | Std
AIC | 0.0068 | 0.0021 | $0.0100^{(2)}$ | 0.0023 | $0.0126^{(2)}$ | 0.0018 | 0.1403 | 0.0091 | 0.3119 | 0.0198 | 0.2694 | 0.0199
BIC | $0.0060^{(1)}$ | 0.0004 | $0.0091^{(1)}$ | 0.0005 | $0.0122^{(1)}$ | 0.0006 | 0.1403 | 0.0091 | 0.3119 | 0.0198 | 0.2694 | 0.0199
SAIC | $0.0067^{(2)}$ | 0.0020 | $0.0100^{(2)}$ | 0.0022 | $0.0126^{(2)}$ | 0.0018 | 0.1403 | 0.0091 | 0.3119 | 0.0198 | 0.2694 | 0.0199
SBIC | $0.0060^{(1)}$ | 0.0004 | $0.0091^{(1)}$ | 0.0005 | $0.0122^{(1)}$ | 0.0006 | 0.1403 | 0.0091 | 0.3119 | 0.0198 | 0.2694 | 0.0199
MAX | 0.0272 | 0.0010 | 0.0264 | 0.0010 | 0.0254 | 0.0012 | 0.1403 | 0.0091 | 0.3119 | 0.0198 | 0.2694 | 0.0199
EQMA | 0.0111 | 0.0006 | 0.0384 | 0.0007 | 0.0616 | 0.0012 | $0.0997^{(1)}$ | 0.0027 | $0.2112^{(1)}$ | 0.0057 | $0.1914^{(1)}$ | 0.0055
LASSO | 0.2617 | 0.0065 | 0.3275 | 0.0072 | 0.2391 | 0.0072 | 0.3170 | 0.0066 | 0.3877 | 0.0058 | 0.4130 | 0.0066
TRMA | $0.0060^{(1)}$ | 0.0004 | $0.0091^{(1)}$ | 0.0005 | $0.0122^{(1)}$ | 0.0006 | $0.1013^{(2)}$ | 0.0024 | $0.2153^{(2)}$ | 0.0047 | $0.1980^{(2)}$ | 0.0053
$n=1000$ and the noise level is 10% of the standard deviation of mean
$\eta_{i}$.
| Signal 1 | Signal 2 | Signal 3 | Signal 4 | Signal 5 | Signal 6
| Mean | Std | Mean | Std | Mean | Std | Mean | Std | Mean | Std | Mean | Std
AIC | 0.0537 | 0.0022 | 0.0527 | 0.0022 | 0.0518 | 0.0025 | 0.1533 | 0.0095 | 0.3156 | 0.0184 | 0.2766 | 0.0203
BIC | 0.0238 | 0.0061 | $0.0289^{(2)}$ | 0.0056 | $0.0350^{(2)}$ | 0.0061 | 0.1533 | 0.0095 | 0.3156 | 0.0184 | 0.2766 | 0.0203
SAIC | 0.0537 | 0.0022 | 0.0527 | 0.0022 | 0.0518 | 0.0025 | 0.1533 | 0.0095 | 0.3156 | 0.0184 | 0.2766 | 0.0203
SBIC | 0.0238 | 0.0060 | $0.0289^{(2)}$ | 0.0056 | $0.0349^{(2)}$ | 0.0061 | 0.1533 | 0.0095 | 0.3156 | 0.0184 | 0.2766 | 0.0203
MAX | 0.0537 | 0.0022 | 0.0527 | 0.0022 | 0.0518 | 0.0025 | 0.1533 | 0.0095 | 0.3156 | 0.0184 | 0.2766 | 0.0203
EQMA | $0.0219^{(2)}$ | 0.0011 | 0.0428 | 0.0013 | 0.0643 | 0.0013 | $0.1033^{(1)}$ | 0.0025 | $0.2113^{(1)}$ | 0.0062 | $0.1937^{(1)}$ | 0.0064
LASSO | 0.2636 | 0.0065 | 0.3285 | 0.0073 | 0.2417 | 0.0074 | 0.3179 | 0.0059 | 0.3892 | 0.0059 | 0.4140 | 0.0068
TRMA | $0.0119^{(1)}$ | 0.0009 | $0.0181^{(1)}$ | 0.0009 | $0.0244^{(1)}$ | 0.0012 | $0.1045^{(2)}$ | 0.0024 | $0.2152^{(2)}$ | 0.0048 | $0.1995^{(2)}$ | 0.0051
$n=1000$ and the noise level is 25% of the standard deviation of mean
$\eta_{i}$.
| Signal 1 | Signal 2 | Signal 3 | Signal 4 | Signal 5 | Signal 6
| Mean | Std | Mean | Std | Mean | Std | Mean | Std | Mean | Std | Mean | Std
AIC | 0.1360 | 0.0051 | 0.1346 | 0.0047 | 0.1326 | 0.0066 | 0.2154 | 0.0115 | 0.3562 | 0.0211 | 0.3269 | 0.0232
BIC | 0.1360 | 0.0051 | 0.1346 | 0.0047 | 0.1326 | 0.0066 | 0.2154 | 0.0115 | 0.3562 | 0.0211 | 0.3269 | 0.0232
SAIC | 0.1360 | 0.0051 | 0.1346 | 0.0047 | 0.1326 | 0.0066 | 0.2154 | 0.0115 | 0.3562 | 0.0211 | 0.3269 | 0.0232
SBIC | 0.1360 | 0.0051 | 0.1346 | 0.0047 | 0.1326 | 0.0066 | 0.2154 | 0.0115 | 0.3562 | 0.0211 | 0.3269 | 0.0232
MAX | 0.1360 | 0.0051 | 0.1346 | 0.0047 | 0.1326 | 0.0066 | 0.2154 | 0.0115 | 0.3562 | 0.0211 | 0.3269 | 0.0232
EQMA | $0.0550^{(2)}$ | 0.0021 | $0.0673^{(2)}$ | 0.0027 | $0.0839^{(2)}$ | 0.0030 | $0.1221^{(2)}$ | 0.0039 | $0.2224^{(1)}$ | 0.0064 | $0.2061^{(1)}$ | 0.0064
LASSO | 0.2719 | 0.0061 | 0.3328 | 0.0078 | 0.2521 | 0.0069 | 0.3219 | 0.0064 | 0.3912 | 0.0058 | 0.4150 | 0.0062
TRMA | $0.0301^{(1)}$ | 0.0021 | $0.0455^{(1)}$ | 0.0026 | $0.0617^{(1)}$ | 0.0031 | $0.1201^{(1)}$ | 0.0035 | $0.2246^{(2)}$ | 0.0055 | $0.2093^{(2)}$ | 0.0056

Notes: The means and the standard deviations are obtained from 100
replications. Signals 1-3 correspond to the case where the correct models
exist in the candidate model set, and Signals 4-6 correspond to the case where
all candidate models are misspecified. The smallest and the second smallest
results of each setting of experiments are flagged by (1) and (2),
respectively.

Next, we analyze the performance of the maximal rank model. Due to the
existence of the noise term, the maximal model is not the best as it will lead
to overfitting. To verify this, we use the same method to generate testing
data with 200 observations, and the prediction RMSE is shown in Table 2. The
maximal rank model has the smallest RMSE on the training data, but it has a
larger RMSE on the testing data. In contrast, TRMA has the smallest RMSE in
most cases. Even when it is not optimal, the difference between TRMA and the
optimal method is very small.

Table 2: RMSEs of $\widehat{\mathbf{y}}$ on training data and testing data
when $n=1000$ and the noise level is 25% of the standard deviation of mean
$\eta_{i}$.

| Signal 1 | Signal 2 | Signal 3
---|---|---|---
| Training | Testing | Training | Testing | Training | Testing
AIC | 0.9139 | 1.1002 | 0.8392 | 1.2097 | 0.7740 | 1.2954
BIC | 0.9329 | $1.0739^{(1)}$ | 0.8612 | $1.1712^{(1)}$ | 0.7841 | $1.2777^{(1)}$
SAIC | 0.9134 | 1.0987 | $0.8384^{(2)}$ | 1.2095 | $0.7739^{(2)}$ | 1.2954
SBIC | 0.9329 | $1.0739^{(1)}$ | 0.8612 | $1.1712^{(1)}$ | 0.7841 | $1.2777^{(1)}$
MAX | $0.3852^{(1)}$ | 2.0130 | $0.4226^{(1)}$ | 1.9758 | $0.4673^{(1)}$ | 1.9128
EQMA | $0.5454^{(2)}$ | 1.2309 | 2.1282 | 2.6530 | 3.2246 | 4.0480
LASSO | 7.5535 | 16.8616 | 11.5289 | 21.0827 | 6.2145 | 15.4437
TRMA | 0.9151 | $1.0747^{(2)}$ | 0.8472 | $1.1726^{(2)}$ | 0.7752 | $1.2779^{(2)}$
| Signal 4 | Signal 5 | Signal 6
| Training | Testing | Training | Testing | Training | Testing
AIC | $2.1414^{(1)}$ | 8.9572 | $4.4436^{(1)}$ | 19.7966 | $4.0468^{(1)}$ | 17.2221
BIC | $2.1414^{(1)}$ | 8.9572 | $4.4436^{(1)}$ | 19.7966 | $4.0468^{(1)}$ | 17.2221
SAIC | $2.1414^{(1)}$ | 8.9572 | $4.4436^{(1)}$ | 19.7966 | $4.0468^{(1)}$ | 17.2221
SBIC | $2.1414^{(1)}$ | 8.9572 | $4.4436^{(1)}$ | 19.7966 | $4.0468^{(1)}$ | 17.2221
MAX | $2.1414^{(1)}$ | 8.9572 | $4.4436^{(1)}$ | 19.7966 | $4.0468^{(1)}$ | 17.2221
EQMA | $3.7510^{(2)}$ | $6.4427^{(1)}$ | $7.3754^{(2)}$ | $13.4480^{(1)}$ | $7.0013^{(2)}$ | $12.3670^{(1)}$
LASSO | 10.9486 | 20.2374 | 14.9444 | 25.0418 | 17.0848 | 26.4597
TRMA | 4.2672 | $6.5670^{(2)}$ | 9.0910 | $13.7247^{(2)}$ | 8.4024 | $12.8091^{(2)}$

Notes: The table displays $\left\|\widehat{\mathbf{y}}-\mathbf{y}\right\|$ on
training data and testing data. The means and the standard deviations of the
RMSEs are obtained from 100 replications. The smallest and the second smallest
results of each setting of experiments are flagged by (1) and (2),
respectively.

Next, we demonstrate the performance of different models under the
$\mathrm{KL}$ loss. Table 3 presents the $\mathrm{KL}$ loss of different
methods under different signals and different sample sizes. As the sample size
$n$ increases, the $\mathrm{KL}$ loss of each method significantly decreases.
When the true model exists, as in Signals 1-3, the $\mathrm{KL}$ loss of each
method is relatively small. However, in the case of model misspecification, as
in Signals 4-6, the $\mathrm{KL}$ loss is relatively large. After comparing
different model selection and model averaging methods, it is clear that the
TRMA method is quite stable, consistently producing either the best or second-
best results.

Table 3: $\mathrm{KL}$ loss of different methods under different signals when
the noise level is 5% of the standard deviation of mean $\eta_{i}$.

| Signal 1 | Signal 2 | Signal 3
---|---|---|---
| 500 | 750 | 1000 | 500 | 750 | 1000 | 500 | 750 | 1000
AIC | 0.20 | $0.14^{(2)}$ | $0.10^{(2)}$ | $0.54^{(1)}$ | 0.31 | $0.23^{(2)}$ | 719.25 | $0.56^{(2)}$ | $0.34^{(2)}$
BIC | $0.18^{(1)}$ | $0.11^{(1)}$ | $0.07^{(1)}$ | $0.54^{(1)}$ | $0.25^{(1)}$ | $0.18^{(1)}$ | 704.99 | $0.50^{(1)}$ | $0.31^{(1)}$
SAIC | $0.20^{(2)}$ | $0.14^{(2)}$ | $0.10^{(2)}$ | $0.54^{(1)}$ | $0.30^{(2)}$ | $0.23^{(2)}$ | 719.25 | $0.56^{(2)}$ | $0.34^{(2)}$
SBIC | $0.18^{(1)}$ | $0.11^{(1)}$ | $0.07^{(1)}$ | $0.54^{(1)}$ | $0.25^{(1)}$ | $0.18^{(1)}$ | 704.99 | $0.50^{(1)}$ | $0.31^{(1)}$
MAX | 786.86 | 3.26 | 1.53 | 1163.78 | 3.68 | 1.45 | 834.88 | 6.35 | 1.33
EQMA | 36.80 | 0.39 | 0.25 | 106.27 | 3.51 | 3.02 | $129.52^{(2)}$ | 9.20 | 7.71
TRMA | $0.18^{(1)}$ | $0.11^{(1)}$ | $0.07^{(1)}$ | $2.71^{(2)}$ | $0.25^{(1)}$ | $0.18^{(1)}$ | $80.92^{(1)}$ | $0.50^{(1)}$ | $0.31^{(1)}$
| Signal 4 | Signal 5 | Signal 6
| 500 | 750 | 1000 | 500 | 750 | 1000 | 500 | 750 | 1000
AIC | 775.95 | 177.34 | 39.84 | 1616.23 | 698.80 | 196.70 | 1718.88 | 664.16 | 148.93
BIC | 775.95 | 177.34 | 39.84 | 1616.23 | 698.80 | 196.70 | 1723.48 | 664.16 | 148.93
SAIC | 775.95 | 177.34 | 39.84 | 1616.23 | 698.80 | 196.70 | 1718.88 | 664.16 | 148.93
SBIC | 775.95 | 177.34 | 39.84 | 1616.23 | 698.80 | 196.70 | 1723.48 | 664.16 | 148.93
MAX | 1111.10 | 177.34 | 39.84 | 1696.11 | 698.80 | 196.70 | 1856.59 | 664.16 | 148.93
EQMA | $126.39^{(2)}$ | $33.59^{(2)}$ | $20.32^{(1)}$ | $320.38^{(2)}$ | $136.97^{(2)}$ | $90.13^{(1)}$ | $313.32^{(2)}$ | $123.99^{(2)}$ | $76.18^{(1)}$
TRMA | $53.25^{(1)}$ | $29.59^{(1)}$ | $21.13^{(2)}$ | $200.74^{(1)}$ | $122.28^{(1)}$ | $93.89^{(2)}$ | $185.60^{(1)}$ | $112.13^{(1)}$ | $81.75^{(2)}$

Notes: The table displays the KL divergence of different methods under various
sample sizes and signals. The means and the standard deviations of the RMSEs
are obtained from 100 replications. Signals 1-3 correspond to the case where
the correct models exist in the candidate model set, and Signals 4-6
correspond to the case where all candidate models are misspecified. The
smallest and the second smallest results of each setting of experiments are
flagged by (1) and (2), respectively.

In Figure 4, the ratio of the $\mathrm{KL}$ loss corresponding to the model
averaging estimator with respect to the smallest $\mathrm{KL}$ loss, i.e.,
$\mathrm{KL}(\widehat{\mathbf{w}})/\inf_{\mathbf{w}}\mathrm{KL}(\mathbf{w})$,
is plotted under model misspecification. It can be seen that as $n$ increases,
the ratio monotonically converges to 1, which confirms the asymptotic
optimality of Theorem 1.

(a) Signal 4

(b) Signal 5

(c) Signal 6

Figure 4: The ratio of KL divergence of TRMA over the infeasible best possible
model averaging.

To check the root$-n$ consistency of the model averaging estimator, we set $n$
to 500, 750, 1000, 2000 and 4000 with the noise level being 5% of the standard
deviations of mean $\eta$. Table 4 shows the values of
$\left\|\mathrm{vec}(\widehat{\mathcal{B}}(\widehat{\mathbf{w}}))-\mathrm{vec}(\mathcal{B}_{0})\right\|$
when correct models are contained. As $n$ grows, we can see that the
difference between $\widehat{\mathcal{B}}(\widehat{\mathbf{w}})$ and
$\mathrm{vec}(\mathcal{B}_{0})$ monotonically converges to 0, which reflects
the consistency of the model averaging estimator.

Table 4: $\left\|\mathrm{vec}(\widehat{\mathcal{B}}(\widehat{\mathbf{w}}))-\mathrm{vec}(\mathcal{B}_{0})\right\|$ of the TRMA estimator. $n$ | 500 | 750 | 1000 | 2000 | 4000
---|---|---|---|---|---
Signal 1 | 0.0103 | 0.0072 | 0.0054 | 0.0041 | 0.0025
Signal 2 | 0.0342 | 0.0107 | 0.0097 | 0.0059 | 0.0038
Signal 3 | 0.2248 | 0.0169 | 0.0128 | 0.0073 | 0.0044

In the experiment, we also record the weights assigned to each model to verify
Theorem 3. It can be seen in Table 5 that the weights assigned to the
underfitted models decrease as $n$ grows. When $n=500$, $\widehat{w}_{\Delta}$
is very large, most likely because of the relatively small sample size, and
the parameters of each candidate model cannot be calculated well. This
phenomenon gradually disappears as $n$ gets larger.

Table 5: $\widehat{w}_{\Delta}$ of the underfitted models.

$n$ | 500 | 750 | 1000 | 2000 | 4000
---|---|---|---|---|---
Signal 2 | 0.0533 | 0.0023 | 0.0021 | 0.0007 | 0.0002
Signal 3 | 0.8602 | 0.0112 | 0.0084 | 0.0021 | 0.0010

Notes: this table does not include Signal 1, as the CP rank of Signal 1 is 1
and no candidate model is underfitted.

###  4.2 3-D simulation

Similar to Section 4.1, we evaluate normal, binomial and Poisson tensor
regression on 3-D images. First we set
$\eta_{i}=<\mathcal{X}_{i},\mathcal{B}>$, where $\mathcal{X}_{i}$ and
$\mathcal{B}\in\mathbb{R}^{32\times 32\times 32}$ are 3-D tensors. Then we
generate the response variable through different models: for the normal model,
$y_{i}\sim\text{Normal}(\eta_{i},1)$; for the binomial model,
$y_{i}\sim\text{Bernoulli}(p_{i})$ with
$p_{i}=\mathrm{exp}(0.1\eta_{i})/[1+\mathrm{exp}(0.1\eta_{i})]$; for the
Poisson model, $y_{i}\sim\text{Poisson}(\mu_{i})$ with
$\mu_{i}=\mathrm{exp}(0.01\eta_{i})$. We assess two 3-D signals, the first
with a rank-2 CP decomposition and the second without a low-rank CP
decomposition. Tables 6-7 show the means and standard deviations of the RMSE
of $\widehat{\mathcal{B}}$ over 100 replications. We do not test the LASSO
method because the previous experiments have shown its poor without a CP
decomposition.

Table 6: RMSEs of $\widehat{\mathcal{B}}$ of 3-D Signal 1.

| Normal | Binomial | Poisson
---|---|---|---
| Mean | Std | Mean | Std | Mean | Std
AIC | $0.0046^{(2)}$ | 0.0007 | 1.4596 | 0.3070 | 0.4623 | 0.0894
BIC | $0.0027^{(1)}$ | 0.0001 | $0.1174^{(1)}$ | 0.0153 | $0.3370^{(2)}$ | 0.0288
SAIC | $0.0046^{(2)}$ | 0.0006 | 1.4595 | 0.3072 | 0.4566 | 0.0892
SBIC | $0.0027^{(1)}$ | 0.0001 | $0.1174^{(1)}$ | 0.0153 | $0.3370^{(2)}$ | 0.0288
MAX | 0.0072 | 0.0003 | 0.8930 | 0.1283 | 1.1687 | 0.0642
EQMA | 0.0199 | 0.0025 | 0.4911 | 0.0402 | 0.3798 | 0.0163
TRMA | $0.0027^{(1)}$ | 0.0002 | $0.1176^{(2)}$ | 0.0149 | $0.2874^{(1)}$ | 0.0203

Notes: The table displays
$\left\|\widehat{\mathcal{B}}-\mathcal{B}_{0}\right\|$ on 3-D Signal 1. In
this case, there are correct models in the set of candidate models. The means
and the standard deviations of the RMSEs are obtained from 100 replications.
The smallest and the second smallest results of each setting of experiments
are flagged by (1) and (2), respectively.

Table 7: RMSEs of $\widehat{\mathcal{B}}$ of 3-D Signal 2.

| Normal | Binomial | Poisson
---|---|---|---
| Mean | Std | Mean | Std | Mean | Std
AIC | 0.1573 | 0.0275 | 1.7578 | 0.1849 | 0.5259 | 0.0932
BIC | 0.1567 | 0.0279 | $0.1806^{(2)}$ | 0.0256 | $0.3651^{(2)}$ | 0.0251
SAIC | 0.1573 | 0.0275 | 1.7578 | 0.1849 | 0.5143 | 0.0964
SBIC | 0.1567 | 0.0279 | $0.1806^{(2)}$ | 0.0256 | $0.3651^{(2)}$ | 0.0251
MAX | 0.1587 | 0.0268 | 0.9934 | 0.0910 | 1.1784 | 0.0665
EQMA | $0.1107^{(1)}$ | 0.0120 | 0.5240 | 0.0385 | 0.4018 | 0.0182
TRMA | $0.1170^{(2)}$ | 0.0104 | $0.1801^{(1)}$ | 0.0255 | $0.3156^{(1)}$ | 0.0187

Notes: The table displays
$\left\|\widehat{\mathcal{B}}-\mathcal{B}_{0}\right\|$ on 3-D Signal 2, where
all candidate models are misspecified. The means and the standard deviations
of the RMSEs are obtained from 100 replications. The smallest and the second
smallest results of each setting of experiments are flagged by (1) and (2),
respectively.

Similar to the 2-D simulation, when the signal has a low-rank CP
decomposition, BIC and SBIC have the best performance, while the difference
between BIC and TRMA is small. However, when the signal has no low-rank CP
decomposition, TRMA performs better. From the results obtained from the
binomial and Poisson models, we find that TRMA not only performs well in the
normal model, but also outperforms in other distribution families. In the
Poisson model, TRMA outperforms AIC and BIC when the correct models exist.
Table 8 also provides the $\mathrm{KL}$ loss of different methods. The
rankings of different methods are similar to the RMSE of
$\widehat{\mathcal{B}}$. Our method gives either the smallest or the second
smallest $\mathrm{KL}$ loss.

Table 8: $\mathrm{KL}$ loss of 3-D Signal 1 and 3-D Signal 2.

| Normal | Binomial | Poisson
---|---|---|---
| Signal 1 | Signal 2 | Signal 1 | Signal 2 | Signal 1 | Signal 2
AIC | 0.3659 | 424.4208 | 5.6251 | 9.3273 | 0.4251 | 0.6409
BIC | $0.1192^{(1)}$ | 420.2293 | $0.1854^{(1)}$ | $0.3039^{(1)}$ | $0.2036^{(2)}$ | $0.2377^{(2)}$
SAIC | 0.3612 | 424.4208 | 5.6341 | 9.3273 | 0.4132 | 0.5992
SBIC | $0.1192^{(1)}$ | 420.2293 | $0.1854^{(1)}$ | $0.3039^{(1)}$ | $0.2036^{(2)}$ | $0.2377^{(2)}$
MAX | 0.8320 | 430.3220 | 5.1145 | 5.9969 | 7.0758 | 7.8166
EQMA | 6.3901 | $207.6599^{(1)}$ | 1.5039 | 2.3838 | 0.2576 | 0.2963
TRMA | $0.1220^{(2)}$ | $229.8151^{(2)}$ | $0.1898^{(2)}$ | $0.3083^{(2)}$ | $0.1407^{(1)}$ | $0.1785^{(1)}$

Notes: The table displays the KL divergence of various methods under different
regression model frameworks and signals. The means and the standard deviations
are obtained from 100 replications. The smallest and the second smallest
results of each setting of experiments are flagged by (1) and (2),
respectively.

Overall, the simulation results from the above experiments demonstrate that
the TRMA method is effective and provides more accurate predictions compared
to methods such as AIC, BIC, and SBIC. This is especially true when the data
has high levels of noise and lacks low-rank decompositions.

###  4.3 Real data examples

In this section, we apply the TRMA method and other methods to analyse two
datasets. The first one is the skin cancer dataset (https://www.isic-
archive.com/). Moles are types of skin growths or lesions, and most of them
are benign and harmless. Environmental factors such as prolonged sun exposure,
or changes in hormone levels such as puberty, may cause moles to darken or
lesion. In this case, the benign mole can develop into melanoma, a serious
skin cancer. The data contain 1800 photos of benign skin moles and 1497 photos
of malignant skin moles (ISIC,, 2022). The purpose of this experiment is to
distinguish whether the moles are benign or malignant. We downsize the origin
photos to a size of $64\times 64\times 3$ for the convenience of calculation,
and use each image as the independent variable $\mathcal{X}_{i}$. The response
variable $y_{i}$ is binary, with 1 representing the malignant mole and 0
representing the benign mole. The set of candidate models consists of rank-1
to rank-5 CP binary tensor regression models. We then apply various model
averaging and model selection methods to these five models. We conduct three
sets of experiments, randomly choosing 50%, 75%, and 90% of the whole data as
the training set to estimate the parameters, and the remaining as the testing
set to compare the misclassification rates of different models. To avoid
errors caused by the number of training samples and uneven data partitioning,
we repeat the experiment 100 times. Figure 5 presents the boxplot of the
misclassification rates on the testing dataset over the 100 replications.

(a) 50% of the data is used as the training set

(b) 75% of the data is used as the training set

(c) 90% of the data is used as the training set

Figure 5: Misclassification rates of different methods on the skin cancer
dataset when 50%, 75% and 90% of the data is used as the training set. The box
in the middle of the plot represents the interquartile range (IQR), which
includes 50% of the data points between the 25th percentile (Q1) and the 75th
percentile (Q3). The median is represented by a horizontal line within the
box. The upper and lower whiskers extending from the box indicate the maximum
and minimum values within 1.5 times the IQR from the edge of the box. Any data
points outside the whiskers are considered outliers and are shown as
individual points beyond the whiskers. The dashed line represents the median
misclassification rate of TRMA.

Figure 5 shows that the TRMA method has the highest accuracy, followed by
EQMA, BIC and SBIC, and the worst performing method is the MAX method. The
performance of the TRMA method is very stable and outperforms other methods
under different proportions of the training set.

The second dataset is the attention deficit hyperactivity disorder (ADHD) data
(http://fcon_1000.projects.nitrc.org/indi/adhd200/). ADHD is one of the most
common neurodevelopmental disorders of childhood. Common ADHD symptoms include
inattention, poor impulse control, and emotional hyperactivity. The data are
obtained from the ADHD-200 Global Competition datasets (ADHD,, 2012), which
have been split into 774 training samples and 172 testing samples. The
training samples contain 285 ADHD subjects and 489 normal controls. Among
them, there are 484 males and 290 females aged from 7 to 21. The raw data have
been preprocessed by standard steps, including skull-stripping, segmentation
with SPM12 and CAT12 toolboxes in MATLAB, etc. The size of the processed MRI
images is $121\times 145\times 121$. We downsize each image to $10\times
12\times 10$ by Haar wavelet transform, and use it as covariate
$\mathcal{X}_{i}$. The response variable $y_{i}$ is binary, with 1
representing ADHD and 0 representing the control group. The set of candidate
models include rank-1 to rank-5 CP binary tensor regression models, and we
apply various model averaging and model selection methods to these five
models. Figure 6 presents the boxplot of the misclassification rates on the
testing dataset over the 100 data replications.

Figure 6: Misclassification rates of different methods on the ADHD-200
dataset. The box in the middle of the plot represents the interquartile range
(IQR), which includes 50% of the data points between the 25th percentile (Q1)
and the 75th percentile (Q3). The median is represented by a horizontal line
within the box. The upper and lower whiskers extending from the box indicate
the maximum and minimum values within 1.5 times the IQR from the edge of the
box. Any data points outside the whiskers are considered outliers and are
shown as individual points beyond the whiskers. The dashed line represents the
median misclassification rate of TRMA.

The TRMA method outperforms BIC and SBIC, while BIC and SBIC perform much
better than the other models. The misclassification rate of the maximal rank
model is the largest. This does not mean that the MRI image has a low-rank CP
decomposition. The high-rank model may overfit the noise part, resulting in
poor prediction results. The misclassification rate of 0.3609 is also better
than that of the Tucker decomposition, the CP decomposition and the Bayes
method mentioned in Li et al., 2018b .

##  5 Discussion

In this article, we propose a model averaging method to avoid selecting a
single poor model for the CP tensor regression by weighting estimators of
different models. The weights are determined by minimizing a CV-type
$\mathrm{KL}$ divergence. We prove that the proposed estimator is
asymptotically optimal when all candidate models are misspecified. When at
least one candidate model is correctly specified, we also prove the TRMA
estimator is root$-n$ consistent and the weights of the model averaging
estimator are assigned to the correct models. The proposed method is
demonstrated to be efficient and stable through both simulations and real data
examples.

However, there are still some areas worthy of further research. First, while
our method is based on the CP tensor regression, other tensor decomposition
methods such as Tucker decomposition and tensor train decomposition also have
their advantages. Combining model averaging with other tensor regression
methods will be a promising avenue for future research. Second, the root$-n$
consistency of the model averaging estimator requires the number of parameters
to be fixed. In practice, especially in clinical data, the sample sizes are
relatively small. It is a challenging question whether the statistic has the
same property after relaxing this restriction. Finally, the article only
considers the case where the data are independent and identically distributed.
For longitudinal data or time series data, the conclusion may be different,
and this is a field that requires a further investigation.

## References

  * ADHD, (2012) ADHD (2012).  The ADHD-200 Sample.  http://fcon_1000.projects.nitrc.org/indi/adhd200/.
  * Akaike, (1973) Akaike, H. (1973).  Maximum likelihood identification of Gaussian autoregressive moving average models.  Biometrika, 60(2):255–265.
  * Ando and Li, (2014) Ando, T. and Li, K.-C. (2014).  A model-averaging approach for high-dimensional regression.  Journal of the American Statistical Association, 109(505):254–265.
  * Ando and Li, (2017) Ando, T. and Li, K.-C. (2017).  A weight-relaxed model averaging approach for high-dimensional generalized linear models.  The Annals of Statistics, 45(6):2654–2679.
  * Arlot and Lerasle, (2016) Arlot, S. and Lerasle, M. (2016).  Choice of $V$ for $V$-fold cross-validation in least-squares density estimation.  The Journal of Machine Learning Research, 17(1):7256–7305.
  * Bi et al., (2018) Bi, X., Qu, A., and Shen, X. (2018).  Multilayer tensor factorization with applications to recommender systems.  The Annals of Statistics, 46(6B):3308–3333.
  * Buckland et al., (1997) Buckland, S. T., Burnham, K. P., and Augustin, N. H. (1997).  Model selection: an integral part of inference.  Biometrics, pages 603–618.
  * Feng and Liu, (2020) Feng, Y. and Liu, Q. (2020).  Nested model averaging on solution path for high-dimensional linear regression.  Stat, 9(1):e317.
  * Feng et al., (2022) Feng, Y., Liu, Q., Yao, Q., and Zhao, G. (2022).  Model averaging for nonlinear regression models.  Journal of Business & Economic Statistics, 40(2):785–798.
  * Gao et al., (2023) Gao, Z., Zou, J., Zhang, X., and Ma, Y. (2023).  Frequentist model averaging for envelope models.  Scandinavian Journal of Statistics.
  * Guo et al., (2011) Guo, W., Kotsia, I., and Patras, I. (2011).  Tensor learning for regression.  IEEE Transactions on Image Processing, 21(2):816–827.
  * Hansen, (2007) Hansen, B. E. (2007).  Least squares model averaging.  Econometrica, 75(4):1175–1189.
  * Hansen and Racine, (2012) Hansen, B. E. and Racine, J. S. (2012).  Jackknife model averaging.  Journal of Econometrics, 167(1):38–46.
  * Harshman, (1970) Harshman, R. (1970).  Foundations of the parafac procedure: Models and conditions for an” explanatory” multimodal factor analysis.  UCLA Working Papers in Phonetics, 16(1):84.
  * Håstad, (1990) Håstad, J. (1990).  Tensor rank is NP-complete.  Journal of Algorithms, 11(4):644–654.
  * Hjort and Claeskens, (2003) Hjort, N. L. and Claeskens, G. (2003).  Frequentist model average estimators.  Journal of the American Statistical Association, 98(464):879–899.
  * Hoeting et al., (1999) Hoeting, J. A., Madigan, D., Raftery, A. E., and Volinsky, C. T. (1999).  Bayesian model averaging: a tutorial.  Statistical Science, 14(4):382–417.
  * ISIC, (2022) ISIC (2022).  Skin cancer dataset.  https://www.isic-archive.com/.
  * Ke et al., (2023) Ke, B., Zhao, W., and Wang, L. (2023).  Smoothed tensor quantile regression estimation for longitudinal data.  Computational Statistics & Data Analysis, 178.
  * Kolda, (2006) Kolda, T. (2006).  Multilinear operators for higher-order decompositions.  Sandia National Laboratories, Albuquerque, NM and Livermore, CA.
  * Kolda and Bader, (2009) Kolda, T. G. and Bader, B. W. (2009).  Tensor decompositions and applications.  SIAM Review, 51(3):455–500.
  * (22) Li, J., Xia, X., Wong, W. K., and Nott, D. (2018a).  Varying-coefficient semiparametric model averaging prediction.  Biometrics, 74(4):1417–1426.
  * (23) Li, X., Xu, D., Zhou, H., and Li, L. (2018b).  Tucker tensor regression and neuroimaging analysis.  Statistics in Biosciences, 10(3):520–545.
  * Liu, (2015) Liu, C.-A. (2015).  Distribution theory of the least squares averaging estimator.  Journal of Econometrics, 186(1):142–159.
  * Liu and Zhang, (2022) Liu, H. and Zhang, X. (2022).  Frequentist model averaging for undirected Gaussian graphical models.  Biometrics, pages 1–13.
  * Liu and Okui, (2013) Liu, Q. and Okui, R. (2013).  Heteroscedasticity-robust $C_{p}$ model averaging.  The Econometrics Journal, 16(3):463–472.
  * Liu et al., (2020) Liu, Q., Yao, Q., and Zhao, G. (2020).  Model averaging estimation for conditional volatility models with an application to stock market volatility forecast.  Journal of Forecasting, 39(5):841–863.
  * Lock, (2018) Lock, E. F. (2018).  Tensor-on-tensor regression.  Journal of Computational and Graphical Statistics, 27(3):638–647.
  * Longford, (2005) Longford, N. T. (2005).  Model Selection and Efficiency: Is ‘Which Model…?’ the Right Question?  Journal of the Royal Statistical Society. Series A (Statistics in Society), 168(3):469–472.
  * Lu and Su, (2015) Lu, X. and Su, L. (2015).  Jackknife model averaging for quantile regressions.  Journal of Econometrics, 188(1):40–58.
  * Schwarz, (1978) Schwarz, G. (1978).  Estimating the dimension of a model.  The Annals of Statistics, 6(2):461–464.
  * Si et al., (2022) Si, Y., Zhang, Y., and Li, G. (2022).  An efficient tensor regression for high-dimensional data.  arXiv preprint arXiv:2205.13734.
  * Tucker, (1966) Tucker, L. R. (1966).  Some mathematical notes on three-mode factor analysis.  Psychometrika, 31(3):279–311.
  * Wan et al., (2010) Wan, A. T., Zhang, X., and Zou, G. (2010).  Least squares model averaging by Mallows criterion.  Journal of Econometrics, 156(2):277–283.
  * Wang et al., (2021) Wang, D., Zheng, Y., and Li, G. (2021).  High-dimensional low-rank tensor autoregressive time series modeling.  arXiv preprint arXiv:2101.04276.
  * Wang et al., (2022) Wang, J., Hou, J., and Eldar, Y. C. (2022).  Tensor robust principal component analysis from multilevel quantized observations.  IEEE Transactions on Information Theory, 69(1):383–406.
  * Yang, (2001) Yang, Y. (2001).  Adaptive regression by mixing.  Journal of the American Statistical Association, 96(454):574–588.
  * Yu and Feng, (2014) Yu, Y. and Feng, Y. (2014).  Modified cross-validation for penalized high-dimensional linear regression models.  Journal of Computational and Graphical Statistics, 23(4):1009–1027.
  * Yuan and Zhang, (2016) Yuan, M. and Zhang, C.-H. (2016).  On tensor completion via nuclear norm minimization.  Foundations of Computational Mathematics, 16(4):1031–1068.
  * Yuan and Yang, (2005) Yuan, Z. and Yang, Y. (2005).  Combining linear regression models: When and how?  Journal of the American Statistical Association, 100(472):1202–1214.
  * Zhang et al., (2019) Zhang, X., Li, L., Zhou, H., Zhou, Y., Shen, D., and ADNI (2019).  Tensor generalized estimating equations for longitudinal imaging analysis.  Statistica Sinica, 29(4):1977–2005.
  * Zhang and Liu, (2019) Zhang, X. and Liu, C.-A. (2019).  Inference after model averaging in linear regression models.  Econometric Theory, 35(4):816–841.
  * Zhang and Liu, (2023) Zhang, X. and Liu, C.-A. (2023).  Model averaging prediction by $K$-fold cross-validation.  Journal of Econometrics, 235(1):280–301.
  * Zhang et al., (2013) Zhang, X., Lu, Z., and Zou, G. (2013).  Adaptively combined forecasting for discrete response time series.  Journal of Econometrics, 176(1):80–91.
  * Zhang et al., (2016) Zhang, X., Yu, D., Zou, G., and Liang, H. (2016).  Optimal model averaging estimation for generalized linear models and generalized linear mixed-effects models.  Journal of the American Statistical Association, 111(516):1775–1790.
  * Zhang et al., (2023) Zhang, X., Zhang, X., and Ma, Y. (2023).  A model-averaging treatment of multiple instruments in Poisson models with errors.  Canadian Journal of Statistics, 51(1):173–198.
  * Zhao et al., (2020) Zhao, S., Liao, J., and Yu, D. (2020).  Model averaging estimator in ridge regression and its large sample properties.  Statistical Papers, 61(4):1719–1739.
  * Zhou et al., (2013) Zhou, H., Li, L., and Zhu, H. (2013).  Tensor regression with applications in neuroimaging data analysis.  Journal of the American Statistical Association, 108(502):540–552.
  * Zou et al., (2022) Zou, J., Wang, W., Zhang, X., and Zou, G. (2022).  Optimal model averaging for divergent-dimensional Poisson regressions.  Econometric Reviews, 41(7):775–805.


# Surface Coil Intensity Correction for MRI

###### Abstract

Modern MRI scanners utilize one or more arrays of small receive-only coils to
collect k-space data. The sensitivity maps of the coils, when estimated using
traditional methods, differ from the true sensitivity maps, which are
generally unknown. Consequently, the reconstructed MR images exhibit undesired
spatial variation in intensity. These intensity variations can be at least
partially corrected using pre-scan data. In this work, we propose an intensity
correction method that utilizes pre-scan data. For demonstration, we apply our
method to a digital phantom, as well as to cardiac MRI data collected from a
commercial scanner by Siemens Healthineers. The code is available at
https://github.com/OSU-MR/SCC.

Index Terms—  MRI, brightness correction, surface coil

##  1 Introduction

MRI is a well-established noninvasive diagnostic tool. To meet the demands of
shorter scan times, finer resolution, and higher signal-to-noise ratio, MRI
hardware as well as data collection and processing methods have evolved
rapidly in the last three decades. In particular, MRI scanners now use arrays
of receive-only surface coils for signal detection. The k-space data are
collected simultaneously across all receive coils, facilitating higher k-space
undersampling and thus shorter scan times. However, the sensitivity maps of
the receive coils are spatially varying and generally not known in advance.
The commonly used methods to find sensitivity maps from the k-space data do
not yield true sensitivity maps but estimates in which the voxel-wise square-
root of the sum-of-squares of the coil sensitivity maps equals one. Although
this constraint allows estimating a unique set of sensitivity maps, the
resulting estimated maps are spatially modulated relative to the true maps.
When these imperfect sensitivity-map estimates are used in any SENSE-based
reconstruction method [1], the reconstructed images also exhibit spatially
varying intensity, leading to incorrect contrast and overly bright and dark
regions. This can be specifically problematic for thoracic and abdominal
imaging in larger patients because the central area of the body can appear
excessively dark compared to the tissue closer to the received coils. Post-
processing methods based on filtering or histogram manipulation have been
proposed to suppress intensity variation but have not been widely adopted [2,
3]. Although modern commercial scanners utilize data from pre-scans to
partially correct the intensity variation, these corrections are not typically
employed in research settings. This leads to mismatch between clinical and
research imaging results. This work presents an open-source implementation for
a Surface Coil intensity Correction method (SCC) for data collected on Siemens
scanners.

##  2 Methods

###  2.1 Problem formulation

In MRI, the data are collected in the spatial frequency domain, called
k-space, simultaneously from $K\geq 1$ receive coils. The data measured from
the $k^{\text{th}}$ receive coil is modeled as

$\displaystyle\bm{y}^{(k)}=\underbrace{\bm{P}\bm{F}\bm{S}^{(k)}}_{\displaystyle\triangleq\bm{A}^{(k)}}\bm{x}+\bm{\eta}^{(k)},$
(1)

where $\bm{y}^{(k)}\in{\mathbb{C}}^{M}$ is measured k-space data,
$\bm{x}\in{\mathbb{C}}^{N}$ is a vectorized 2D or 3D image,
$\bm{S}^{(k)}\in{\mathbb{C}}^{N\times N}$ is a diagonal matrix representing
the sensitivity map, $\bm{F}\in{\mathbb{C}}^{N\times N}$ is the 2D or 3D
discrete Fourier transform, $\bm{P}$ is the sampling matrix, which contains
$M$ rows from an $N\times N$ identity matrix, and
$\bm{\eta}^{(k)}\in{\mathbb{C}}^{M}$ is additive circularly symmetric Gaussian
noise with variance $\sigma^{2}$. By vertically stacking $\bm{y}^{(k)}$,
$\bm{A}^{(k)}$, and $\bm{\eta}^{(k)}$ into $\bm{y}\in{\mathbb{C}}^{KM}$,
$\bm{A}\in{\mathbb{C}}^{KM\times N}$, and $\bm{\eta}\in{\mathbb{C}}^{KM}$,
respectively, the data measured from all $K$ coils can be expressed as

$\displaystyle\bm{y}=\bm{A}\bm{x}+\bm{\eta}.$ (2)

To estimate $\bm{x}$ in (2), the regularized least-square (RL) solution is a
common choice, which yields

$\displaystyle\bm{x}_{\text{RL}}=\operatorname*{arg\,min}_{\bm{x}}\frac{1}{\sigma^{2}}\|\bm{A}\bm{x}-\bm{y}\|_{2}^{2}+\mathcal{R}(\bm{x}),$
(3)

where the regularization term $\mathcal{R}(\bm{x})$ encourages
$\bm{x}_{\text{RL}}$ to adhere to a prior belief about the set of true
$\bm{x}$.

To construct $\bm{A}$ in (3), one must know the true sensitivity maps
$\\{\bm{S}^{(k)}\\}_{k=1}^{K}$. In practice, the true sensitivity maps are
unknown, and so one typically uses an estimate
$\\{\widehat{\bm{S}}^{(k)}\\}_{k=1}^{K}$ produced by ESPIRiT [4], by the
method of Walsh et al. [5], or by simply dividing individual coil images by
the square-root of the sum-of-squares (SSoS) [1], leading to

$\displaystyle\widehat{\bm{x}}=\operatorname*{arg\,min}_{\bm{x}}\frac{1}{\sigma^{2}}\|\widehat{\bm{A}}\bm{x}-\bm{y}\|_{2}^{2}+\mathcal{R}(\bm{x}),$
(4)

where
$\widehat{\bm{A}}={\big{[}{\widehat{\bm{A}}^{(1){\mathsf{T}}}},\dots,{\widehat{\bm{A}}^{(K){\mathsf{T}}}}\big{]}}^{\mathsf{T}}$
and $\widehat{\bm{A}}^{(k)}=\bm{P}\bm{F}\widehat{\bm{S}}^{(k)}$.

With ESPIRiT, SSoS, or similar methods, the estimates are constrained to
$\widehat{\bm{S}}^{\mathsf{H}}\widehat{\bm{S}}=\bm{I}$ for
$\widehat{\bm{S}}=\big{[}\widehat{\bm{S}}^{(1){\mathsf{T}}},\dots,\widehat{\bm{S}}^{(K){\mathsf{T}}}\big{]}^{\mathsf{T}}$.
Because the true sensitivity maps generally do not obey
$\bm{S}^{\mathsf{H}}\bm{S}=\bm{I}$, the estimated ones do not match the truth.

In SSoS, $\widehat{\bm{S}}$ are estimated by elementwise dividing each coil
image $\bm{S}^{(k)}\bm{x}$ with the square-root of the sum-of-squares image.
Therefore, there exists a diagonal matrix $\bm{G}\in{\mathbb{R}}^{N\times N}$
such that $\bm{S}^{(k)}=\widehat{\bm{S}}^{(k)}\bm{G}$ holds for all non-zero
voxels in the noise-less case. In the typical case that $\bm{G}\neq\bm{I}$,
the recovered $\widehat{\bm{x}}$ in (4) incurs unwanted spatial variation
compared to $\bm{x}$. This problem can be addressed by estimating $\bm{G}$
from pre-scan data and integrating that estimate, $\widehat{\bm{G}}$, into the
forward model. In particular, we propose to replace $\widehat{\bm{S}}^{(k)}$
in (4) with the corrected map $\widehat{\bm{S}}^{(k)}\widehat{\bm{G}}$ for
each $k\in\\{1,\dots,K\\}$, giving

$\displaystyle\widehat{\bm{x}}_{G}=\operatorname*{arg\,min}_{\bm{x}}\frac{1}{\sigma^{2}}\|\widehat{\bm{A}}\widehat{\bm{G}}\bm{x}-\bm{y}\|_{2}^{2}+\mathcal{R}(\bm{x}).$
(5)

Note that (5) corrects the sensitivity maps prior to the reconstruction of
$\bm{x}$. We also propose an alternative approach that first reconstructs
$\bm{x}$ using uncorrected sensitivity maps and then corrects the
reconstructed image with an estimate, $\widehat{\bm{H}}$, of the true
intensity correction map $\bm{H}\in{\mathbb{R}}^{N\times N}$, i.e.,

$\displaystyle\widehat{\bm{x}}_{H}=\widehat{\bm{H}}\widehat{\bm{x}},$ (6)

where $\bm{H}$ and $\widehat{\bm{H}}$ are both diagonal. In the following
section, we describe our proposed “SCC” method of estimating $\bm{G}$ and
$\bm{H}$ from a short pre-scan that is invariably performed right before the
primary data acquisition stage.

3D body coil images, $\underline{\bm{x}}_{\text{b}}$ coil combining correction
map using (15) interpolation, $\mathcal{I}(\widehat{\bm{g}}_{s};\bm{q})$
reconstruction using (5) 3D surface coil images,
$\underline{\bm{x}}_{\text{s}}$ coil combining MRI data, $\bm{y}$, slice
location, $\bm{q}$
$\bm{x}_{\text{bc}}$$\widehat{\bm{g}}_{\text{s}}$$\widehat{\bm{G}}$$\widehat{\bm{x}}_{G}$$\bm{x}_{\text{sc}}$$\bm{q}$$\bm{y}$pre-
scan(a)

3D body coil images, $\underline{\bm{x}}_{\text{b}}$ coil combining correction
map using (16) interpolation, $\mathcal{I}(\widehat{\bm{h}}_{s};\bm{q})$
correction using (6) 3D surface coil images, $\underline{\bm{x}}_{\text{s}}$
coil combining MRI data, $\bm{y}$, slice location, $\bm{q}$ an image
reconstruction
$\bm{x}_{\text{bc}}$$\widehat{\bm{h}}_{\text{s}}$$\widehat{\bm{H}}$$\widehat{\bm{x}}_{H}$$\bm{x}_{\text{sc}}$$\bm{q}$$\bm{y}$$\widehat{\bm{x}}$pre-
scan(b)

Fig. 1: Two flavors of the proposed SCC method. In (a), the correction is
applied to the sensitivity maps that are then used for image reconstruction.
In (b), the correction is applied to images that have already been
reconstructed. The content of the shaded boxes is available in the raw data
file.

###  2.2 Estimating correction maps from pre-scan

Most clinical MRI scanners have two sets of RF coils: smaller surface receive
coils that are placed close to the patient, and larger transmit/receive “body
coils” that are hidden behind the bore liner. Since the advent of surface
receive coils with high channel counts, body coils are primarily used for RF
transmission. However, during the pre-scan, which is performed immediately
before data acquisition, fully sampled low-resolution 3D images are
simultaneously collected from surface and body coils. The main goal of this
pre-scan is to facilitate the estimation of correction maps $\bm{G}$ or
$\bm{H}$.

Let $\bm{x}_{\text{3D}}\in{\mathbb{C}}^{L}$ denote the true, low-resolution 3D
image with $L$ voxels, and let the diagonal matrix
$\bm{S}_{\text{b}}^{(j)}\in{\mathbb{C}}^{L\times L}$ be the true sensitivity
map of the $j^{\text{th}}$ body coil, with $j=1,2,\dots,J$. Then, the stack of
true body-coil images $\underline{\bm{x}}_{\text{b}}\in{\mathbb{C}}^{JL}$ is
given by $\underline{\bm{x}}_{\text{b}}=\bm{S}_{\text{b}}\bm{x}_{\text{3D}}$,
where
$\bm{S}_{\text{b}}=\big{[}\bm{S}_{\text{b}}^{(1){\mathsf{T}}},\bm{S}_{\text{b}}^{(2){\mathsf{T}}},\dots,\bm{S}_{\text{b}}^{(J){\mathsf{T}}}\big{]}^{\mathsf{T}}\in{\mathbb{C}}^{JL\times
L}$.

Let $\widehat{\bm{S}}_{\text{b}}\in{\mathbb{C}}^{JL\times L}$ be the estimate
of $\bm{S}_{\text{b}}$ constructed using SSoS. We assume that
$\widehat{\bm{S}}_{\text{b}}\approx\bm{S}_{\text{b}}$ because the body coils
are large and have approximately uniform sensitivity maps. In this case, we
can rewrite
$\underline{\bm{x}}_{\text{b}}=\bm{S}_{\text{b}}\bm{x}_{\text{3D}}$ as

$\displaystyle\underline{\bm{x}}_{\text{b}}$
$\displaystyle\approx\widehat{\bm{S}}_{\text{b}}\bm{x}_{\text{3D}}.$ (7)

By multiplying both sides of (7) with
$\widehat{\bm{S}}_{\text{b}}^{\mathsf{H}}$, we get

$\displaystyle\underbrace{\widehat{\bm{S}}_{\text{b}}^{\mathsf{H}}\underline{\bm{x}}_{\text{b}}}_{\displaystyle\triangleq\bm{x}_{\text{bc}}}$
$\displaystyle\approx\underbrace{\widehat{\bm{S}}_{\text{b}}^{\mathsf{H}}\widehat{\bm{S}}_{\text{b}}}_{\displaystyle=\bm{I}}\bm{x}_{\text{3D}},$
(8)

which implies

$\displaystyle\bm{x}_{\text{bc}}\approx\bm{x}_{\text{3D}}.$ (9)

Let the diagonal matrix $\bm{S}_{\text{s}}^{(k)}\in{\mathbb{C}}^{L\times L}$
be the true sensitivity map of the $k^{\text{th}}$ surface coil, with
$k=1,2,\dots K$. Then the surface-coil image stack
$\underline{\bm{x}}_{\text{s}}\in{\mathbb{C}}^{KL}$ is given by
$\underline{\bm{x}}_{\text{s}}=\bm{S}_{\text{s}}\bm{x}_{\text{3D}}$, where
$\bm{S}_{\text{s}}=\big{[}\bm{S}_{\text{s}}^{(1){\mathsf{T}}},\bm{S}_{\text{s}}^{(2){\mathsf{T}}},\dots,\bm{S}_{\text{s}}^{(K){\mathsf{T}}}\big{]}^{\mathsf{T}}\in{\mathbb{C}}^{KL\times
L}$. Let $\widehat{\bm{S}}_{\text{s}}\in{\mathbb{C}}^{KL\times L}$ represent
the coil maps estimated using SSoS. Then, there exists a diagonal
$\bm{G}_{\text{s}}\in{\mathbb{C}}^{L\times L}$ such that
$\bm{S}_{\text{s}}=\widehat{\bm{S}}_{\text{s}}\bm{G}_{\text{s}}$ for all non-
zero voxels in the noise-less case, in which case

$\displaystyle\underline{\bm{x}}_{\text{s}}$
$\displaystyle=\widehat{\bm{S}}_{\text{s}}\bm{G}_{\text{s}}\bm{x}_{\text{3D}}.$
(10)

By multiplying both sides of (10) with
$\widehat{\bm{S}}_{\text{s}}^{\mathsf{H}}$, we get

$\displaystyle\underbrace{\widehat{\bm{S}}_{\text{s}}^{\mathsf{H}}\underline{\bm{x}}_{\text{s}}}_{\displaystyle\triangleq\bm{x}_{\text{sc}}}$
$\displaystyle=\underbrace{\widehat{\bm{S}}_{\text{s}}^{\mathsf{H}}\widehat{\bm{S}}_{\text{s}}}_{\displaystyle=\bm{I}}\bm{G}_{\text{s}}\bm{x}_{\text{3D}},$
(11)

which yields the coil-combined surface-coil image

$\displaystyle\bm{x}_{\text{sc}}$
$\displaystyle=\bm{G}_{\text{s}}\bm{x}_{\text{3D}}.$ (12)

Combining (9) and (12), we get

$\displaystyle\bm{x}_{\text{sc}}\approx\bm{G}_{\text{s}}\bm{x}_{\text{bc}}.$
(13)

Since $\bm{G}_{\text{s}}$ is a diagonal matrix, (13) can also be written as

$\displaystyle\bm{x}_{\text{sc}}\approx\bm{X}_{\text{bc}}\bm{g}_{\text{s}}$
(14)

with
$\bm{X}_{\text{bc}}\triangleq\operatorname{Diag}(\bm{x}_{\text{bc}})\in{\mathbb{R}}^{L\times
L}$ and
$\bm{g}_{\text{s}}\triangleq\operatorname{diag}(\bm{G}_{\text{s}})\in{\mathbb{R}}^{L}$.

We solve for $\bm{g}_{\text{s}}$ in (14) using regularized least-squares:

$\displaystyle\widehat{\bm{g}}_{\text{s}}=\operatorname*{arg\,min}_{\bm{g}_{\text{s}}}\|\bm{X}_{\text{bc}}\bm{g}_{\text{s}}-\bm{x}_{\text{sc}}\|_{2}^{2}+\phi(\bm{g}_{\text{s}}),$
(15)

where $\phi(\cdot)$ is a regularization term that promotes spatial smoothness.
The 3D correction map $\widehat{\bm{g}}_{\text{s}}$ can be interpolated along
a 2D imaging plane to yield a 2D correction map
$\widehat{\bm{g}}=\mathcal{I}\big{(}\widehat{\bm{g}}_{\text{s}};\bm{q}\big{)}\in{\mathbb{R}}^{N}$,
where $\bm{q}$ describes the orientation and position of the 2D imaging plane
and $\mathcal{I}(\cdot;\bm{q})$ represents interpolation along that plane.
Finally, $\widehat{\bm{g}}$ can be embedded along the diagonal of
$\widehat{\bm{G}}$ and used in (5) to obtain $\widehat{\bm{x}}_{G}$.

Alternatively, one can rewrite (13) as
$\bm{x}_{\text{bc}}\approx\bm{H}_{\text{s}}\bm{x}_{\text{bc}}$, with diagonal
$\bm{H}_{\text{s}}\in{\mathbb{R}}^{L\times L}$ and estimate
$\bm{h}_{\text{s}}\triangleq\operatorname{diag}(\bm{H}_{\text{s}})$ via

$\displaystyle\widehat{\bm{h}}_{\text{s}}=\operatorname*{arg\,min}_{\bm{h}_{\text{s}}}\|\bm{X}_{\text{sc}}\bm{h}_{\text{s}}-\bm{x}_{\text{bc}}\|_{2}^{2}+\phi(\bm{h}_{\text{s}}),$
(16)

where
$\bm{X}_{\text{sc}}\triangleq\operatorname{Diag}(\bm{x}_{\text{sc}})\in{\mathbb{R}}^{L\times
L}$. The 3D correction map $\widehat{\bm{h}}_{\text{s}}$ can be interpolated
along a 2D imaging plane to yield a 2D correction map
$\widehat{\bm{h}}=\mathcal{I}\big{(}\widehat{\bm{h}}_{\text{s}};\bm{q}\big{)}\in{\mathbb{R}}^{N}$.
Then $\widehat{\bm{h}}$ can be embedded along the diagonal of
$\widehat{\bm{H}}$ for use in (6) to yield $\widehat{\bm{x}}_{H}$.

Fig. 1 describes the processing pipeline for producing corrected images
$\widehat{\bm{x}}_{G}$ and $\widehat{\bm{x}}_{H}$ using SCC.

body coil

surface coil

Fig. 2: Layout of the digital phantom with four small surface coils and two
large body coils.

###  2.3 Implementation of SCC

SCC was implemented with the regularization term
$\phi(\cdot)\triangleq\lambda\|\left[\nabla_{x}^{\mathsf{T}}(\cdot),\nabla_{y}^{\mathsf{T}}(\cdot),\nabla_{z}^{\mathsf{T}}(\cdot)\right]^{\mathsf{T}}\|_{2}^{2}$,
where $\nabla_{x}(\cdot)$, $\nabla_{y}(\cdot)$, and $\nabla_{z}(\cdot)$
compute first-order finite differences along the three directions of the input
and $\lambda>0$ controls the regularization strength. Since both terms in (15)
and (16) are quadratic, these optimization problems can be solved using
closed-form expressions. Despite implementing $\bm{X}_{\text{bc}}$ and
$\bm{X}_{\text{sc}}$ as sparse matrices, we observed that executing the
closed-form solution was slow. In response, we solved the optimization problem
in (15) and (16) using the conjugate gradient (CG) method. For 3D low-
resolution images of size $64\times 64\times 64$, the CPU-only computation
time for either of the optimization problems was less than 5 seconds on a
standard workstation.

For real MRI data, $\underline{\bm{x}}_{\text{b}}$,
$\underline{\bm{x}}_{\text{s}}$, $\bm{y}$, and $\bm{q}$ were extracted from
the raw data file. After discarding oversampling in the readout direction, the
field-of-view (FOV) and matrix size for $\underline{\bm{x}}_{\text{b}}$ and
$\underline{\bm{x}}_{\text{s}}$ were $500\times 500\times 500~{}\text{mm}^{3}$
and $64\times 32\times 32$, respectively, along the frequency encoding and two
phase-encoding directions. The pre-scan data in the discrete Fourier domain
were apodized with the Tukey window function to reduce Gibbs ringing and then
zero-padded to deliver an isotropic resolution of $7.8\times 7.8\times
7.8~{}\text{mm}^{3}$ and a matrix size of $64\times 64\times 64$. The multi-
coil $\underline{\bm{x}}_{\text{s}}$ and $\underline{\bm{x}}_{\text{b}}$
images were coil-combined using SSoS, resulting in real-valued
$\bm{x}_{\text{sc}}$ and $\bm{x}_{\text{bc}}$, respectively, of size $64\times
64\times 64$. Both $\bm{x}_{\text{sc}}$ and $\bm{x}_{\text{bc}}$ were
normalized by dividing with $\max\\!\left[\bm{x}_{\text{bc}}\right]$ for (15)
and $\max\\!\left[\bm{x}_{\text{sc}}\right]$ for (16). The value of $\lambda$
was set at $5\times 10^{-2}$. A Python implementation of SCC is available at
https://github.com/OSU-MR/SCC.

(a)

(b)

(c)

(d)

(e)

Fig. 3: A simulation example showing the impact applying SCC. (a) True
sensitivity maps simulated using the Biot-Savart law. (b) Element-wise
multiplication of true sensitivity maps with a $256\times 256$ digital
phantom. (c) Sensitivity maps estimated using SSoS. (d) From left to right:
reference magnitude image $|\bm{x}|$, magnitude image $|\widehat{\bm{x}}|$
recovered using SENSE with SSoS maps, $12\times|\widehat{\bm{x}}|$, and
$2\times|\widehat{\bm{x}}|^{1/2}$. (e) From left to right: the correction map
$\widehat{\bm{g}}$, $|\widehat{\bm{x}}_{G}|$, the correction map
$\widehat{\bm{h}}$, and $|\widehat{\bm{x}}_{H}|$. All grayscale images are
displayed on a scale of 0 to 1.

$\widehat{\bm{g}}$$\widehat{\bm{h}}$$|\widehat{\bm{x}}|$$|\widehat{\bm{x}}_{G}|$$|\widehat{\bm{x}}_{H}|$

Fig. 4: From left to right, two correction maps $\widehat{\bm{g}}$ and
$\widehat{\bm{h}}$, the magnitude of the uncorrected image, and the magnitudes
of the two corrected images. A representative two-chamber view of the heart is
shown.

##  3 Experiment and Results

###  3.1 Simulation study

In this study, we illustrate how SCC can improve the fidelity of the
reconstructed image using a digital phantom. To this end, we simulated a
$256\times 256$ image, $\bm{x}$, depicted in Fig. 2. The true sensitivity maps
from four surface coils were generated using the Biot-Savart law. The coils
were assumed to be circular wire loops, with a radius that is 20% of the FOV,
with the cylindrical axes of the coils pointing towards the center of the
image. The true sensitivity maps were element-wise multiplied with $\bm{x}$ to
generate individual coil images. Using SSoS as sensitivity maps, the image,
$\widehat{\bm{x}}$, was recovered using SENSE. As shown in the Fig. 3d,
$\widehat{\bm{x}}$ exhibits significant intensity variation, with the center
of the image appearing dark. One could adjust the image intensity by histogram
manipulation, but such measures do not fully address the problem. The true
sensitivity maps (not shown) of two larger body coils, each with a radius
equal to the FOV, were also generated using the Biot-Savart law. To simulate
pre-scan, low-resolution data comprising of $32\times 32$ center of k-space
were generated from both four small surface coils
($\underline{\bm{x}}_{\text{s}}$) and two large body coils
($\underline{\bm{x}}_{\text{b}}$). Following the procedure described in the
previous section, two correction maps were constructed using SCC. However, the
interpolation step was skipped because the digital phantom is already 2D. The
last row in Fig. 3 shows $\widehat{\bm{g}}$, $\widehat{\bm{x}}_{G}$,
$\widehat{\bm{h}}$, and $\widehat{\bm{x}}_{H}$. Both $\widehat{\bm{x}}_{G}$
and $\widehat{\bm{x}}_{H}$ show marked improvement in the uniformity of image
intensity. The normalized mean squared error (NMSE) defined as
$\mathrm{NMSE}(\bm{x},\widehat{\bm{x}})=20\log_{10}\left(\|\bm{x}-\widehat{\bm{x}}\|/\|\bm{x}\|\right)$
was $-5.67$ dB, $-27.64$ dB, and $-27.63$ dB for $\widehat{\bm{x}}$,
$\widehat{\bm{x}}_{G}$, and $\widehat{\bm{x}}_{H}$, respectively. The NMSE of
$\widehat{\bm{x}}_{H}$ or $\widehat{\bm{x}}_{G}$ is primarily limited by the
uniformity of $\underline{\bm{x}}_{\text{b}}$, which, in turn, depends on the
body-coil size.

###  3.2 Cardiac MRI study

For measured MRI data, single-shot free-breathing T1-weighted images were
collected from a healthy volunteer on a 3T scanner (MAGNETOM Vida, Siemens
Healthcare, Erlangen, Germany) running on an XA31A platform. The data were
collected without contrast using phase-sensitive inversion recovery sequence
[6] and were prospectively undersampled at the acceleration rate of two. The
raw data file, containing the pre-scan as well as the k-space measurements
$\bm{y}$, was copied from the scanner for offline processing in Python.
Representative results are shown in Fig. 4. Both $\widehat{\bm{x}}_{G}$ and
$\widehat{\bm{x}}_{H}$ offer improvement over the uncorrected image
$\widehat{\bm{x}}$.

###  3.3 Discussion

The quality of SCC correction depends on the spatial homogeneity of
$\bm{x}_{\text{bc}}$. Due to the much larger size of the body coils, one
expects $\bm{x}_{\text{bc}}$ to be more uniform than
$\underline{\bm{x}}_{\text{s}}$. In practice, however, $\bm{x}_{\text{bc}}$
may not offer the perfect uniformity, leading to residual intensity variation
in $\widehat{\bm{x}}_{G}$ and $\widehat{\bm{x}}_{H}$. The results presented in
Fig. 3 and Fig. 4 show little difference between $\widehat{\bm{x}}_{G}$ and
$\widehat{\bm{x}}_{H}$. These reconstructions are based on SENSE without
regularization. However, for deep learning-based methods or methods that use
regularization, e.g., compressed sensing [7], we expect $\widehat{\bm{x}}_{G}$
and $\widehat{\bm{x}}_{H}$ to differ. This is because $\widehat{\bm{x}}_{H}$
has better spatial uniformity than $\widehat{\bm{x}}$, which can potentially
make $\widehat{\bm{x}}_{H}$ more amenable to sparsity-based or generative
priors.

##  4 Conclusions

MRI reconstructions performed offline typically do not utilize the pre-scan
data to normalize image intensity, leading to uneven image intensity, which
can be especially problematic for abdominal and thoracic imaging from larger
subjects. We propose an optimization-based method to correct surface coil
intensity and provide an open-source Python implementation. Future efforts
will include a demonstration of SCC in other MRI applications and with
scanners from other vendors.

##  5 Additional Information

Acknowledgment: This work was funded by NIH grants R01EB029957 and
R01HL151697.

Compliance with ethical standards: This study was performed in line with the
principles of the Declaration of Helsinki. Approval was granted by the
Institutional Review Board (IRB) of The Ohio State University (2019H0076).

## References

  * [1] Klaas P Pruessmann, Markus Weiger, Markus B Scheidegger, and Peter Boesiger,  “SENSE: Sensitivity encoding for fast MRI,”  Magnetic Resonance in Medicine, vol. 42, no. 5, pp. 952–962, 1999\.
  * [2] James W Murakami, Cecil E Hayes, and ED Weinberger,  “Intensity correction of phased-array surface coil images,”  Magnetic Resonance in Medicine, vol. 35, no. 4, pp. 585–590, 1996\.
  * [3] C DeCarli, DGM Murphy, D Teichberg, G Campbell, and GS Sobering,  “Local histogram correction of MRI spatially dependent image pixel intensity nonuniformity,”  Journal of Magnetic Resonance Imaging, vol. 6, no. 3, pp. 519–528, 1996.
  * [4] Martin Uecker, Peng Lai, Mark J Murphy, Patrick Virtue, Michael Elad, John M Pauly, Shreyas S Vasanawala, and Michael Lustig,  “ESPIRiT–An eigenvalue approach to autocalibrating parallel MRI: Where SENSE meets GRAPPA,”  Magnetic Resonance in Medicine, vol. 71, no. 3, pp. 990–1001, 2014\.
  * [5] David O Walsh, Arthur F Gmitro, and Michael W Marcellin,  “Adaptive reconstruction of phased array MR imagery,”  Magnetic Resonance in Medicine, vol. 43, no. 5, pp. 682–690, 2000\.
  * [6] Peter Kellman, Andrew E Arai, Elliot R McVeigh, and Anthony H Aletras,  “Phase-sensitive inversion recovery for detecting myocardial infarction using gadolinium-delayed hyperenhancement,”  Magnetic Resonance in Medicine, vol. 47, no. 2, pp. 372–383, 2002\.
  * [7] Michael Lustig, David Donoho, and John M Pauly,  “Sparse MRI: The application of compressed sensing for rapid MR imaging,”  Magnetic Resonance in Medicine, vol. 58, no. 6, pp. 1182–1195, 2007\.