# Hybrid sample size calculations for cluster randomised trials using
assurance

S. Faye Williamson1111Corresponding author. Email:
mailto:faye.williamson<EMAIL_ADDRESS>Svetlana V. Tishkovskaya2, Kevin J. Wilson3
1 Biostatistics Research Group, Population Health Sciences Institute,
Newcastle University, UK
2 Faculty of Health and Care, Lancashire Clinical Trials Unit, University of
Central Lancashire, UK
3 School of Mathematics, Statistics & Physics, Newcastle University, UK

###### Abstract

Background/Aims: Sample size determination for cluster randomised trials
(CRTs) is challenging because it requires robust estimation of the intra-
cluster correlation coefficient (ICC). Typically, the sample size is chosen to
provide a certain level of power to reject the null hypothesis in a two sample
hypothesis test. This relies on the minimal clinically important difference
(MCID) and estimates for the overall standard deviation, the ICC and, if
cluster sizes are assumed to be unequal, the coefficient of variation of the
cluster size. Varying any of these parameters can have a strong effect on the
required sample size. In particular, it is very sensitive to small differences
in the ICC. A relevant ICC estimate is often not available, or the available
estimate is imprecise due to being based on studies with low numbers of
clusters. If the ICC value used in the power calculation is far from the
unknown true value, this could lead to trials which are substantially over- or
under-powered.

Methods: In this paper, we propose a hybrid approach using Bayesian assurance
to determine the sample size for a CRT in combination with a frequentist
analysis. Assurance is an alternative to traditional power, which incorporates
the uncertainty on key parameters through a prior distribution. We suggest
specifying prior distributions for the overall standard deviation, ICC and
coefficient of variation of the cluster size, while still utilising the MCID.
We illustrate the approach through the design of a CRT in post-stroke
incontinence and compare the results to those obtained from a standard power
calculation.

Results: We show that assurance can be used to calculate a sample size based
on an elicited prior distribution for the ICC, whereas a power calculation
discards all of the information in the prior except for a single point
estimate. Results show that this approach can avoid misspecifying sample sizes
when the prior medians for the ICC are very similar, but the underlying prior
distributions exhibit quite different behaviour. Incorporating uncertainty on
all three of the nuisance parameters, rather than only on the ICC, does not
notably increase the required sample size.

Conclusion: Assurance provides a better understanding of the probability of
success of a trial given a particular MCID and can be used instead of power to
produce sample sizes which are more robust to parameter uncertainty. This is
especially useful when there is difficulty obtaining reliable parameter
estimates.

Keywords: Assurance, Bayesian design, cluster randomised trials, expected
power, hybrid approach, intra-cluster correlation, minimal clinically
important difference, sample size determination.

##  1 Background

Cluster randomised trials (CRTs) are a type of randomised controlled trial
(RCT) in which randomisation is at the cluster-level, rather than the
individual-level as in standard RCTs. This means that groups of individuals
(e.g. general practices, schools or communities) are randomly allocated to
different interventions (e.g. vaccination programmes or behavioural
interventions). A common reason for implementing this design is to mitigate
the risk of contamination or where individual randomisation is not feasible.
Other justifications are detailed in Eldridge and Kerry, (2012).

Individuals within a cluster are likely to share similar characteristics (e.g.
demographics), as well as be exposed to extraneous factors unique to the
cluster (e.g. delivery of the intervention by the same healthcare
professional). Consequently, outcomes from members of the same cluster are
often correlated, which can be quantified by the intra-cluster correlation
coefficient (ICC). This lack of independence reduces the statistical power
compared to a standard RCT of the same size, meaning that the sample size
needs to be inflated to allow for the clustering effect.

Various methods for sample size determination in CRTs exist (see Rutterford et
al.,, 2015; Gao et al.,, 2015), which all rely on estimation of the ICC. In
practice, ICC estimates are typically based on pilot studies, but these are
often too small to provide precise and reliable estimates (Eldridge et al.,,
2016). An alternative simple approach is to use a conservative estimate of the
ICC (e.g. the upper confidence interval limit) in the sample size calculation
(Browne,, 1995). However, this can lead to over-powered and unnecessarily
large trials. A more reliable method is to combine ICC estimates from multiple
sources, such as previous trials or databases listing ICC estimates (e.g.
Moerbeek and Teerenstra,, 2015), and use information on patterns in ICCs (see
Korevaar et al.,, 2021). This raises further issues such as how to effectively
combine the ICC estimates, how to adequately reflect their varying degrees of
relevance to the planned trial and how to capture the uncertainty in the
individual ICC estimates. Lewis and Julious, (2021) suggest integrating over a
range of possible ICC values, determined by confidence intervals obtained
using methods in Ukoumunne, (2002), to provide an ‘average’ sample size with
respect to the ICC. However, this does not consider the uncertainty present in
other design parameters, such as the treatment effect and variability of the
outcome measures. Further, it assumes that each value of the ICC is equally
likely.

Utilising a Bayesian approach for the trial design, in which prior
distributions are assigned to the unknown design parameters such as the ICC,
could further circumvent these issues and is particularly useful in settings
where ICC estimates are not readily available. In the CRT literature, prior
distributions for the ICC have been proposed based on subjective beliefs
(Spiegelhalter,, 2001) and single or multiple ICC estimates (Turner et al.,,
2004), which may be weighted by relevance of outcomes and patient population
(Turner et al.,, 2005). These are used to estimate a distribution for the
power of the planned trial for a given sample size. Within the Bayesian
framework, uncertainty in other design parameters can be incorporated into the
sample size calculation in a similar way, and the relative likelihood of
different parameter values is encompassed through specification of the prior
distribution. For example, Sarkodie et al., (2023) assign a prior to the
overall standard deviation, in addition to the ICC, then describe a ‘hybrid’
approach to determine the sample size required to attain a desired ‘expected
power’, defined as a weighted average of the probability that the null
hypothesis is rejected (with weights determined by the priors).

Hybrid approaches, which combine a Bayesian design with a frequentist analysis
of the final trial data, have gained increasing popularity, particularly with
respect to standard RCTs (Kunzmann et al.,, 2021). In this paper, we adopt a
hybrid approach by using the Bayesian concept of assurance to determine the
sample size for a two-arm parallel-group CRT with a Wald test for the
analysis. In contrast to traditional frequentist power, which represents a
conditional probability that the trial is a success, given the values chosen
for the design parameters and the hypothesised treatment effect, assurance
typically refers to the unconditional probability that the trial will be
‘successful’ (Chen and Ho,, 2017). We modify this definition by conditioning
on the minimal clinically important difference (MCID) instead of assigning a
prior distribution to, and integrating over, the treatment effect (as in
Kunzmann et al.,, 2021; Ciarleglio et al.,, 2016). This is more representative
of the design stage of a trial, in which the treatment effect is typically
fixed a priori by investigators. Moreover, this ensures that the assurance
will tend to one as the sample size increases so can be used analogously to
traditional power, thus aiding interpretation.

A key consideration when applying a Bayesian design is how to specify suitable
prior distributions. In contrast to the paper by Sarkodie et al., (2023),
which assumes independent priors on the ICC and standard deviation, we suggest
a joint prior distribution for these parameters, as described in Sections 2.3
and 3.2.2. In addition, we account for the fact that many CRTs have unequal
cluster sizes by defining a prior distribution on the coefficient of variation
of cluster size. This is often overlooked in standard sample size calculations
for CRTs (Eldridge et al.,, 2006; Zhan et al.,, 2021).

Our approach is motivated by a parallel-group CRT, Identifying Continence
OptioNs after Stroke (ICONS), outlined in Section 3.1. In Section 3.2, we
illustrate the effects of redesigning this trial using the entire ICC prior
distribution to inform sample size determination via an assurance calculation,
rather than relying on a single point estimate from this distribution as in
Tishkovskaya et al., (2023). The impacts of varying the ICC prior
distributions on the chosen sample size are evaluated in Section 3.3. We
perform sensitivity analyses on other design parameters in an additional
simulation study provided in the Appendix.

Jones et al., (2021) summarise the current state of play regarding the use of
Bayesian methods in CRTs. In doing so, they highlight the “need for further
Bayesian methodological development in the design and analysis of CRTs …in
order to increase the accessibility, availability and, ultimately, use of the
approach.” This paper is therefore a timely contribution.

##  2 Methods

###  2.1 Analysis for CRTs

Suppose that we are designing a two-arm, parallel-group CRT assuming 1:1
randomisation of clusters and normally distributed outcomes. A common analysis
following the trial is to use a linear mixed-effects model. That is, if
$Y_{ij}$ is the response for individual $i=1,\ldots,n_{j}$ in cluster
$j=1,\ldots,C$, then

$Y_{ij}=\alpha+X_{j}\delta+c_{j}+e_{ij},$ (1)

where $\alpha$ is an intercept term; $X_{j}$ is a binary variable which takes
the value 1 if cluster $j$ is allocated to the treatment arm and 0 if it is
allocated to the control arm, so that $\delta$ represents the treatment
effect; $c_{j}\sim\textrm{N}(0,\sigma^{2}_{b})$ is a random cluster effect
with $\sigma^{2}_{b}$ denoting the between cluster variation and
$e_{ij}\sim\textrm{N}(0,\sigma^{2}_{w})$ is the individual-level error with
$\sigma^{2}_{w}$ denoting the within cluster variation. This can be extended
for stepped-wedge CRTs by following the model in Hussey and Hughes, (2007).

The ratio of the variability between clusters $\sigma^{2}_{b}$ to the total
variability $\sigma^{2}=\sigma^{2}_{b}+\sigma^{2}_{w}$ determines the extent
to which clustering induces correlations between outcomes for individuals in
the same cluster. This is referred to as the intra-cluster correlation
coefficient (ICC), $\rho=\sigma^{2}_{b}/\sigma^{2}$ (Kerry and Bland,, 1998).

The superiority of the treatment is assessed via a hypothesis test of
$H_{0}:\delta\leq 0$ versus $H_{1}:\delta>0$. Using a Wald test, assuming
asymptotic normality, the test statistic is
$Z=\hat{\delta}/\sqrt{\textrm{Var}(\hat{\delta})}$, where $\hat{\delta}$ is
the estimate of $\delta$ and
$\textrm{Var}(\hat{\delta})=4\sigma^{2}[1+\\{(\nu^{2}+1)\bar{n}-1\\}\rho]/C\bar{n}$,
where $\bar{n}$ is the average sample size per cluster and $\nu$ is the
coefficient of variation of cluster size, i.e. the ratio of the standard
deviation of cluster sizes to the mean cluster size.

###  2.2 Choosing a sample size using assurance

The power of the one-sided Wald test for significance level $\alpha$ can be
approximated (Eldridge et al.,, 2016) by

$P(n\mid\delta,\bm{\psi})=\Phi\left(\delta\sqrt{\dfrac{C\bar{n}}{4\sigma^{2}[1+\\{(\nu^{2}+1)\bar{n}-1\\}\rho]}}-z_{1-\alpha}\right),$
(2)

where $z_{1-\alpha}$ is the $100(1-\alpha)\%$ percentile of the standard
normal distribution and $\bm{\psi}=(\sigma,\rho,\nu)$ is the vector of
“nuisance” parameters, excluding the treatment effect. For a two-sided Wald
test, $z_{1-\alpha}$ would be replaced by $z_{1-\alpha/2}$. For equal cluster
sizes, the power function would take the same form as equation (2), with
$\nu=0$ and $\bar{n}=n_{j}=n$.

In a standard power calculation, the sample size would be chosen as the
smallest value which gives 80% or 90% power, based on values for
$\bm{\theta}=(\delta,\bm{\psi})$. The treatment effect $\delta$ could be
specified as the MCID or an estimate based on a pilot study, similar
historical trials or expert knowledge. The values used for $\bm{\psi}$ are
typically estimates.

Alternatively, we can use assurance to choose the sample size. Whereas the
power is conditioned on the chosen estimates for $\bm{\psi}$ and possibly
$\delta$, the assurance represents the unconditional probability that an RCT
will achieve a successful outcome (O’Hagan and Stevens,, 2001). Assurance has
been used almost exclusively when the value to be used for $\delta$ is an
estimate. In this case, suppose that the CRT is a success if the null
hypothesis is rejected by the Wald test for $\delta$. Rather than using point
estimates for $\bm{\theta}$, we could assign a prior distribution
$\pi(\bm{\theta})$ to it and define the assurance $A(n)$ as the power,
averaged over the uncertainty in $\bm{\theta}$:

$\displaystyle A(n)$ $\displaystyle=$
$\displaystyle\int_{\bm{\theta}}\Pr(H_{0}\textrm{
rejected}\mid\bm{\theta})\pi(\bm{\theta})d\bm{\theta},$ (3) $\displaystyle=$
$\displaystyle\int_{\bm{\theta}}P(n\mid\bm{\theta})\pi(\bm{\theta})d\bm{\theta}.$

One disadvantage of the assurance is that it tends to $\Pr(\delta>0)$ under
$\pi(\delta)$ as the sample size increases. That is, unlike power, there may
be no sample size for which the assurance is above the typical thresholds of
80% or 90%. Kunzmann et al., (2021) avoid this by conditioning the prior
distribution for $\delta$ on $\delta>0$ in the assurance calculation. In this
paper, we consider the following alternative approach.

The assurance in (3) assumes that we choose $\delta$ in the sample size
calculation based on a priori considerations of the likelihood of the
treatment effect. Instead, we consider the assurance in conjunction with a
trial planned using the relevance argument, that is, using the MCID for
$\delta$, $\delta_{M}$. In this case, there is no need to define a prior
distribution for $\delta$, and the assurance reduces to

$\displaystyle A(n\mid\delta_{M})$ $\displaystyle=$
$\displaystyle\int_{\bm{\psi}}P(n\mid\delta_{M},\bm{\psi})\pi(\bm{\psi})d\bm{\psi}.$

The advantage of this is that the assurance will now tend to one as the sample
size increases.

To evaluate the assurance in practice, we sample values of
$(\bm{\psi}_{j})_{j=1,\ldots,S}$ from the prior distribution $\pi(\bm{\psi})$
for some large number of samples $S$, and use Monte Carlo simulation to
approximate the assurance as

$\displaystyle\tilde{A}(n\mid\delta_{M})$ $\displaystyle=$
$\displaystyle\dfrac{1}{S}\sum_{j=1}^{S}P(n\mid\delta_{M},\bm{\psi}_{j}),$ (4)
$\displaystyle\approx$
$\displaystyle\dfrac{1}{S}\sum_{j=1}^{S}\Phi\left(\delta_{M}\sqrt{\dfrac{C\bar{n}}{4\sigma_{j}^{2}[1+\\{(\nu_{j}^{2}+1)\bar{n}-1\\}\rho_{j}]}}-z_{1-\alpha}\right).$

###  2.3 Specification of priors

To evaluate the assurance, we are required to specify a prior distribution for
$\bm{\psi}$. This simplifies to specifying marginal prior distributions for
each parameter if they can be assumed independent. Given that $\sigma^{2}$ and
$\rho$ are both functions of $\sigma^{2}_{w}$ and $\sigma^{2}_{b}$, it is
unlikely that $\sigma$ and $\rho$ can be assumed independent. Therefore, we
consider a joint prior distribution for $(\sigma,\rho)$ and a marginal prior
distribution for $\nu$. In order for the assurance to be a meaningful
representation of the probability that the null hypothesis is rejected, these
prior distributions should be informative, representing the current state of
knowledge about the possible parameter values. This is an elicitation problem,
and information to specify the priors can be obtained from relevant past data,
expert knowledge or a combination (an example of this is given in Section 3).

Since the coefficient of variation can only take positive values, a gamma
distribution $\nu\sim\textrm{Gamma}(a_{\nu},b_{\nu})$ is a sensible choice for
its prior distribution. The hyperparameters $a_{\nu}$ and $b_{\nu}$ could be
chosen based on previous studies, via modelling or by eliciting expert
knowledge (Eldridge et al.,, 2016).

One way to specify a joint prior distribution for $(\sigma,\rho)$ is to assign
independent priors to $\sigma^{2}_{b}$ and $\sigma^{2}_{w}$, which will induce
a correlation between $\rho$ and $\sigma^{2}$. If we sample values of
$\sigma^{2}_{b}$ and $\sigma^{2}_{w}$ from their priors, we can obtain samples
from the joint prior of $(\sigma,\rho)$. Typical choices of prior
distributions for $\sigma^{2}_{b}$ and $\sigma^{2}_{w}$ are (inverse) gamma
distributions because they provide conjugacy.

An alternative approach, relevant to our application, is to specify the joint
distribution between $\rho$ and $\sigma$ directly. For example, we can utilise
a bivariate copula to encode the dependence between the parameters. A
bivariate copula is a joint distribution function on $[0,1]^{2}$ with standard
uniform marginal distributions (Nelson,, 2006). It can be used to construct a
joint prior for $\rho$ and $\sigma$ via

$\pi_{\rho,\sigma}(\rho,\sigma)=\pi_{\rho}(\rho)\pi_{\sigma}(\sigma)c(u,v),$

where $\pi_{\rho}$ and $\pi_{\sigma}$ are marginal prior distributions,
$c(u,v)$ is the bivariate copula density function evaluated at
$u=F_{\rho}(\rho)$ and $v=F_{\sigma}(\sigma)$ for prior cumulative
distribution functions (CDFs) $F_{\rho}$ and $F_{\sigma}$. One simple choice
is the Gaussian copula:

$c(u,v)=\dfrac{\partial^{2}}{\partial u\partial
v}\Phi_{\gamma}(\Phi^{-1}(u),\Phi^{-1}(v)),$

where $\Phi_{\gamma}$ is the CDF of the bivariate standard normal distribution
with correlation $\gamma$, and $\Phi^{-1}$ is the inverse univariate standard
normal CDF. The advantage of this structure is that it allows specification of
the marginal prior distributions for $\rho$ and $\sigma$ separately to their
dependence, which is given by $\gamma$.

##  3 Application

###  3.1 The ICONS post-stroke incontinence CRT

The approach developed in this paper is motivated by a planned parallel-group
CRT, “Identifying Continence OptioNs after Stroke” (ICONS), which investigates
the effectiveness of a systematic voiding programme in secondary care versus
usual care on post-stroke urinary incontinence for people admitted to NHS
stroke units (Thomas et al.,, 2015). The primary outcome is the severity of
urinary incontinence at three months post-randomisation, measured using the
International Consultation on Incontinence Questionnaire (Avery et al.,,
2004). Although a feasibility trial, ICONS-I (Thomas et al.,, 2014), was
conducted, the resulting ICC estimate was of low precision and could not be
used as a reliable single source to inform the planning of the proposed trial.

ICONS therefore considered a Bayesian approach to combine multiple ICC
estimates from 16 previous related trials. The opinions of eight experts
regarding the relevance of the previous ICC estimates were elicited (O’Hagan,,
2019) and used to assign weights to each study and each outcome within a
study. The elicited study and outcome weights were combined using mathematical
aggregation (O’Hagan et al.,, 2006) and incorporated into a Bayesian
hierarchical model following the method in Turner et al., (2005). The
resulting constructed ICC distribution had a posterior median of
$\hat{\rho}=0.0296$ with a $95\%$ credible interval of $(0.00131,0.330)$.
Details of the expert elicitation process and modelling are described in
Tishkovskaya et al., (2023).

For the ICONS CRT, the sample size was chosen to give 80% power with a 5%
significance level to detect $\delta_{M}=2.52$ using a two-tailed independent-
samples $t$-test and a common standard deviation $\sigma$ of 8.32 obtained
from the ICONS-I feasibility trial. The ICC was assumed to be less than or
equal to $\hat{\rho}=0.0296$. It was assessed as realistic to recruit between
40 and 50 stroke units, which required total sample sizes of $N=480$ and
$N=450$, respectively, and an average sample size per cluster of $n=12$ and
$n=9$, respectively. The original sample size calculation assumed equally
sized clusters. However, if we consider unequal cluster sizes with $\nu=0.49$
(obtained from ICONS-I) and apply the Wald test, the required sample sizes
remain the same.

###  3.2 Redesigning the ICONS CRT using assurance

We consider assurance as an alternative to power to determine the sample size
for the ICONS CRT. This seems like a more natural approach given the
uncertainty in the ICC and the extensive elicitation and modelling that was
conducted to construct the ICC posterior distribution (which forms the prior
distribution for the assurance-based sample size calculation). Moreover,
assurance incorporates the full ICC distribution into the sample size
calculation, rather than relying on a single point estimate from it as in the
power calculation.

We consider the following two forms of assurance.

####  3.2.1 Assurance based on the ICC prior only

In the first case, we fix $(\sigma,\nu)$ using the point estimates obtained
from ICONS-I and only consider the assurance with respect to the ICC. We
sample $S=10,000$ values of $\rho$ from its distribution (see Figure 1) and
approximate the assurance using (4).

Figure 1: Histogram of 10,000 samples of the ICC, $\rho$

To obtain an assurance of 80%, the resulting average sample sizes per cluster
are $\bar{n}=17$ for $C=40$ clusters ($C=20$ per arm) and $\bar{n}=11$ for
$C=50$ clusters ($C=25$ per arm), requiring total sample sizes of $N=680$ and
$N=550$, respectively (Table 1). Thus, the inclusion of uncertainty in the ICC
results in a larger sample size than when using the posterior median ICC, but
provides a more realistic and robust study design. Compared to the classical
approach, the total sample size attained is smaller for the smaller number of
clusters.

The left-hand side plot of Figure 2 illustrates the trade-off between cluster
size and assurance/power, for $C=40$ clusters ($C=20$ per arm). The power
calculation based on the median from the elicited prior distribution of $\rho$
is represented by the red curve and the assurance with a prior on $\rho$ only
by the black curve. We see that the assurance requires a larger sample size
than power when the target lies above 0.5. We also include the power curve
corresponding to the commonly used approach of taking the median of the 34 ICC
estimates (blue line). For a target power of 0.8 (horizontal line), Figure 2
shows that this method requires a larger sample size per cluster than the
aforementioned methods.

Figure 2: Power and assurance curves for the ICONS CRT (left). The power using
the posterior median ICC is red, the power using the median ICC from the 34
ICC estimates is light blue, the assurance with a prior only on the ICC is
black and the assurance with a prior on all of the nuisance parameters
$\bm{\psi}$ is green. The effect of varying the coefficient of variation $\nu$
on the assurance (right). $\nu$ varies between 0 (top curve) and 1 (bottom
curve), with the red line at $\nu=0.5$. The horizontal line indicates the
desired power/assurance. Each plot corresponds to $C=40$ ($C=20$ per arm).

We illustrate the effect of changing $\nu=(0,0.1,\ldots,1)$ on the assurance
in the right-hand side plot of Figure 2. The red curve corresponds to
$\nu=0.5$, the top curve to $\nu=0$ and the bottom curve to $\nu=1$. As $\nu$
increases, the assurance decreases for a given cluster size. We see that the
estimate of $\nu$ has a relatively strong effect on the assurance, and hence
the required sample size. This implies that it needs to be estimated
accurately, or its uncertainty should be taken into account in the assurance
calculation.

####  3.2.2 Assurance based on the prior for $\psi$

In the second case, we obtain the sample size required using an assurance
calculation which averages over a prior distribution on $\sigma$ and $\nu$, as
well as the ICC. Using the data from ICONS-I, we give $\sigma$ and $\nu$ gamma
marginal prior distributions, centred at their estimated values of 8.32 and
0.49, respectively. The standard deviations of the prior distributions are
chosen to represent a belief that $\sigma$ is very likely to be in the range
$[5,11]$ and $\nu$ is very likely to be in the range $[0.3,0.7]$.
Specifically, $\sigma\sim\textrm{Gamma}(a_{\sigma},b_{\sigma})$ and
$\nu\sim\textrm{Gamma}(a_{\nu},b_{\nu})$, where
$a_{\cdot}=m_{\cdot}^{2}/v_{\cdot}$ and $b_{\cdot}=m_{\cdot}/v_{\cdot}$,
$m_{\sigma}=8.32,v_{\sigma}=1^{2}$ and $m_{\nu}=0.49,v_{\nu}=0.066^{2}$.

To incorporate the dependence between $\rho$ and $\sigma$, we utilise the
Gaussian copula with $\gamma=0.44$. This is chosen to be consistent with the
correlation between $\rho$ and $\sigma$ that would result from independent
prior distributions on the between and within group variances of
$\sigma^{2}_{b}\sim\textrm{Gamma}(0.18,0.04)$ and
$\sigma^{2}_{w}\sim\textrm{Gamma}(21.06,0.32)$, respectively. The
hyperparameters of these two gamma prior distributions are chosen to provide
the correct marginal means and variances for $\rho$ and $\sigma$. To sample
values of $\rho$ and $\sigma$ from their joint prior distribution, we repeat
the following steps:

  1. 1.

Sample $(x_{i},y_{i})$, $i=1,\ldots,S$ from $\textrm{N}_{2}(\bm{0},R)$, where
$R$ is the prior correlation matrix with diagonal elements 1 and off-diagonal
elements $\gamma=0.44$.

  2. 2.

Calculate $(\rho_{i},\sigma_{i})$ as
$\left(F^{-1}_{\rho}(\Phi(x)),F^{-1}_{\sigma}(\Phi(y))\right)$.

The quantile function $F^{-1}_{\sigma}$ is that of the relevant normal
distribution. The empirical quantile function $F^{-1}_{\rho}$ is used for
$\rho$, based on the 10,000 prior samples.

The resulting joint prior distribution for $(\sigma,\rho)$ and marginal prior
distribution for $\nu$ are illustrated in Figure 3. We see that the marginal
prior for $\rho$ remains as in Figure 1, but the samples are positively
correlated with the values of $\sigma$.

Figure 3: The joint prior distribution between $\rho$ and $\sigma$ (left) and
the marginal prior distribution for $\nu$ (right), based on 10,000 samples.

The resulting average cluster sample sizes for an assurance of 80% are
$\bar{n}=18$ for $C=40$ clusters ($C=20$ per arm) and $\bar{n}=12$ for $C=50$
clusters ($C=25$ per arm), requiring total sample sizes of $N=720$ and
$N=600$, respectively. By incorporating uncertainty on $\sigma$ and $\nu$, as
well as $\rho$, the sample size increases only slightly, as illustrated in the
left-hand side plot of Figure 2 (green line). To achieve a target assurance of
80% (dashed horizontal line), the average sample size required per cluster
increases from 17 to 18 when $C=40$; an increase in total sample size of
approximately 5%.

Table 1 summarises the sample sizes required to attain a target
power/assurance of 80% for the various approaches applied to the ICONS trial.
“Classical approach” refers to the multiple-estimate method of taking the
median of the ICC estimates without taking the relevance of the different
studies into account. Relative to the classical approach that is often used in
practice, the total sample size required when using the assurance-based method
remains the same whilst incorporating uncertainty on all three parameters.

Method | Priors |  | Total number
---
of clusters, $C$
| Mean cluster
---
size, $\bar{n}$
| Total sample
---
size, $N$
| Power
---
(classical approach)
NA |  | 50
---
40
| 10
---
18
| 500
---
720
| Power
---
(based on posterior median)
NA |  | 50
---
40
| 9
---
12
| 450
---
480
Assurance | $\rho$ |  | 50
---
40
| 11
---
17
| 550
---
680
Assurance | $\bm{\psi}=(\sigma,\rho,\nu)$ |  | 50
---
40
| 12
---
18
| 600
---
720
Table 1: Summary of sample sizes obtained for the ICONS CRT based on power and
assurance calculations.

###  3.3 Sensitivity analysis for the ICC prior

In the above, we consider the ICC prior distribution based on all eight
reviewers and all 16 relevant studies. In this section, we investigate the
sensitivity of the assurance-based sample size (with priors on $\bm{\psi}$) to
varying assumptions on the reviewers and relevant studies, and compare this to
the sensitivity of the sample sizes from power calculations (using the
posterior median ICC).

To recognise uncertainty in the individual reviewers’ responses, and in how
these responses were pooled, the mathematical aggregation was refitted with
alternative reviewer importance weights: equal weights of 0.125 for all
reviewers and using a rank sum approach (see Tishkovskaya et al.,, 2023). For
the rank sum approach, we use Cronbach’s alpha score and assign ranks to each
reviewer according to this score. In addition, we rerun the Bayesian
hierarchical model for only the top 4 (25%), 8 (50%) and 12 (75%) most
relevant studies. We refer to the five variations of the original ICC prior
distribution as: equal weights, differentiated weights, top 4, top 8 and top
12.

The differentiated weights prior (red) and equal weights prior (green) are
provided alongside the original prior (black) in the left-hand side plot of
Figure 4. The top 4 prior (red), top 8 prior (green) and top 12 prior (blue)
are given alongside the original prior (black) in the right-hand side plot of
Figure 4. In both plots, the prior medians are given by vertical dashed lines.

Figure 4: Left: The densities of the differentiated weights (red), equal
weights (green) and original ICC prior (black). Right: The densities of the
top 4 (red), top 8 (green), top 12 (blue) and original ICC prior (black). The
prior medians are represented by vertical dashed lines.

We see that the ICC prior remains similar to the original prior whether
differentiated weights or equal weights are used, although both alternative
weightings assign more probability to the ICC taking larger values. There is a
larger change when using the top 4, top 8 or top 12 studies. In each case, the
alternative prior is more diffuse than the original prior. Relatively large
changes in the prior can cause only small changes in the prior median (e.g.
the original prior compared to the top 12 prior). The effects of the
alternative ICC priors on the sample sizes are shown in Table 2.

| $C=50$ | $C=40$
---|---|---
| Assurance | Power | Assurance | Power
ICC Estimate/Prior | $\bar{n}$ | $N$ | $\bar{n}$ | $N$ | $\bar{n}$ | $N$ | $\bar{n}$ | $N$
Original | 12 | 600 | 9 | 450 | 18 | 720 | 12 | 480
Differentiated weights | 13 | 650 | 10 | 500 | 20 | 800 | 13 | 520
Equal weights | 13 | 650 | 9 | 450 | 19 | 760 | 13 | 520
Top 4 | 12 | 600 | 8 | 400 | 18 | 720 | 11 | 440
Top 8 | 14 | 700 | 9 | 450 | 23 | 920 | 12 | 480
Top 12 | 15 | 750 | 9 | 450 | 24 | 960 | 12 | 480
Table 2: The average sample size per cluster $\bar{n}$ and the total sample
size $N$ required for the ICONS CRT using assurance (with priors on
$\bm{\psi}$) and power based on the original ICC estimate/prior and five
alternative estimates/priors when $C=50$ and $C=40$. The power is based on the
posterior median of the ICC.

We see smaller changes in sample sizes for $C=50$ than $C=40$ using assurance.
Overall, we observe larger changes in sample size using assurance than power
based on the prior median of the ICC. This illustrates the risk with using
just the median; it takes no account of the prior probability that the ICC
could be relatively large, so has the potential to systematically
underestimate the required sample size. In contrast, the assurance-based
sample size is sensitive to the entire ICC prior distribution, particularly
the upper tail.

To illustrate this point, compare the original ICC prior (black) to the top 12
prior (blue) in the right hand side of Figure 4. They have substantially
different priors, resulting in large differences in sample sizes required
under assurance (600 versus 750 when $C=50$, respectively). However, their
prior medians are almost identical, resulting in identical sample size
requirements under power (450 when $C=50$).

In the Appendix, we further evaluate the properties of the hybrid approach
compared to power via a simulation study.

##  4 Conclusions

A standard sample size calculation requires pre-specification of parameters
that are unknown at the design stage of a trial. Unique to sample size
calculations for typical CRTs is the ICC, which requires robust estimation to
avoid over- or under-powering the trial. Unnecessarily high ICC values, for
example, lead to inefficient trials, increasing the number of clusters and/or
participants, and overall trial costs. In practice, parameter uncertainty is
typically not considered, which can be problematic given the sensitivity of
the sample size to small differences in the ICC.

This paper proposes an alternative approach to sample size determination for
CRTs using the Bayesian concept of assurance to incorporate parameter
uncertainty into the design. The advantage of this approach is that it yields
designs that provide adequate power across the likely range of parameter
values and is therefore more robust to parameter misspecification. This is
particularly important when there is difficulty obtaining a reliable ICC
estimate, such as in the ICONS post-stroke incontinence CRT used to motivate
this work.

We assign prior distributions to the ICC, overall standard deviation and
coefficient of variation of the cluster size, whilst setting the treatment
effect equal to the MCID in line with standard practice. We consider a joint
prior for the ICC and standard deviation to model the dependency between these
parameters. In the motivating case-study, we use the entire ICC prior
distribution elicited from expert opinion and data from previous studies to
inform the sample size. Further work could consider using a commensurate prior
to synthesise multiple sources of pre-trial information on the ICC, as in
Zheng et al., (2023).

Sensitivity analyses of the assurance-based sample size to different ICC
priors showed that different behaviour of the prior, particularly in the upper
tail, can have quite a strong effect on the required sample size. Using a
point estimate from this prior, for example the median, can miss this overall
behaviour and result in sample sizes which are systematically too small, based
on current knowledge about the ICC. Additional sensitivity analyses conducted
on the overall standard deviation showed that the greater the uncertainty
expressed in the prior, the more robust the assurance-based sample size is
(see Appendix).

Uncertainty in the treatment effect can also be incorporated into the
assurance calculation in a similar way. This may be appropriate for non-
inferiority trials, for example, where the non-inferiority margin is fixed in
advance and the treatment difference can be considered a nuisance parameter.

In line with regulatory requirements, we have maintained a frequentist
analysis to present a hybrid framework. Further work could consider a fully
Bayesian approach by using assurance when the success criterion is based on
the posterior distribution of the treatment effect (e.g. Spiegelhalter,, 2001)

The hybrid approach presented in this paper can be applied to avoid
incorrectly powered studies resulting from ill-estimated model parameters, to
mitigate the impact of uncertainty in the ICC and other nuisance parameters,
and to incorporate expert opinion or historical data when designing a CRT.

## References

  * Avery et al., (2004) Avery, K., Donovan, J., Peters, T. J., Shaw, C., Gotoh, M., and Abrams, P. (2004).  ICIQ: a brief and robust measure for evaluating the symptoms and impact of urinary incontinence.  Neurourology and Urodynamics: Official Journal of the International Continence Society, 23(4):322–330.
  * Browne, (1995) Browne, R. H. (1995).  On the use of a pilot sample for sample size determination.  Statistics in Medicine, 14:1933–1940.
  * Chen and Ho, (2017) Chen, D.-G. and Ho, S. (2017).  From statistical power to statistical assurance: It’s time for a paradigm change in clinical trial design.  Communications in Statistics - Simulation and Computation, 46(10):7957–7971.
  * Ciarleglio et al., (2016) Ciarleglio, M. M., Arendt, C. D., and Peduzzi, P. N. (2016).  Selection of the effect size for sample size determination for a continuous response in a superiority clinical trial using a hybrid classical and Bayesian procedure.  Clinical Trials, 13(3):275–285.
  * Eldridge and Kerry, (2012) Eldridge, S. and Kerry, S. (2012).  A practical guide to cluster randomised trials in health services research.  John Wiley & Sons.
  * Eldridge et al., (2006) Eldridge, S. M., Ashby, D., and Kerry, S. (2006).  Sample size for cluster randomized trials: effect of coefficient of variation of cluster size and analysis method.  International journal of epidemiology, 35(5):1292–1300.
  * Eldridge et al., (2016) Eldridge, S. M., Costelloe, C. E., Kahan, B. C., Lancaster, G. A., and Kerry, S. M. (2016).  How big should the pilot study for my cluster randomised trial be?  Statistical methods in medical research, 25(3):1039–1056.
  * Gao et al., (2015) Gao, F., Earnest, A., Matchar, D. B., Campbell, M. J., and Machin, D. (2015).  Sample size calculations for the design of cluster randomized trials: a summary of methodology.  Contemporary clinical trials, 42:41–50.
  * Hussey and Hughes, (2007) Hussey, M. A. and Hughes, J. P. (2007).  Design and analysis of stepped wedge cluster randomized trials.  Contemporary clinical trials, 28(2):182–191.
  * Jones et al., (2021) Jones, B. G., Streeter, A. J., Baker, A., Moyeed, R., and Creanor, S. (2021).  Bayesian statistics in the design and analysis of cluster randomised controlled trials and their reporting quality: a methodological systematic review.  Systematic Reviews, 10(91):1–14.
  * Kerry and Bland, (1998) Kerry, S. M. and Bland, J. M. (1998).  The intracluster correlation coefficient in cluster randomisation.  BMJ, 316(7142):1455–1460.
  * Korevaar et al., (2021) Korevaar, E., Kasza, J., Taljaard, M., Hemming, K., Haines, T., Turner, E. L., Thompson, J. A., Hughes, J. P., and Forbes, A. B. (2021).  Intra-cluster correlations from the CLustered OUtcome Dataset bank to inform the design of longitudinal cluster trials.  Clinical Trials, 18(5):529–540.
  * Kunzmann et al., (2021) Kunzmann, K., Grayling, M., Lee, K., Robertson, D., Rufibach, K., and Wason, J. (2021).  A review of Bayesian perspectives on sample size derivation for confirmatory trials.  The American Statistician, 75:424–432.
  * Lewis and Julious, (2021) Lewis, J. and Julious, S. A. (2021).  Sample sizes for cluster-randomised trials with continuous outcomes: Accounting for uncertainty in a single intra-cluster correlation estimate.  Statistical methods in medical research, 30(11):2459–2470.
  * Moerbeek and Teerenstra, (2015) Moerbeek, M. and Teerenstra, S. (2015).  Power analysis of trials with multilevel data.  CRC Press.
  * Nelson, (2006) Nelson, R. (2006).  An Introduction to Copulas.  Springer-Verlag, New York.
  * O’Hagan, (2019) O’Hagan, A. (2019).  Expert knowledge elicitation: Subjective but scientific.  The American Statistician, 73(sup1):69–81.
  * O’Hagan et al., (2006) O’Hagan, A., Buck, C. E., Daneshkhah, A., Eiser, J. R., Garthwaite, P. H., Jenkinson, D. J., Oakley, J. E., and Rakow, T. (2006).  Uncertain judgements: eliciting experts’ probabilities.
  * O’Hagan and Stevens, (2001) O’Hagan, A. and Stevens, J. (2001).  Bayesian assessment of sample size for clinical trials of cost-effectiveness.  Medical Decision Making, 21:219–230.
  * Rutterford et al., (2015) Rutterford, C., Copas, A., and Eldridge, S. (2015).  Methods for sample size determination in cluster randomized trials.  International journal of epidemiology, 44(3):1051–1067.
  * Sarkodie et al., (2023) Sarkodie, S. K., Wason, J. M., and Grayling, M. J. (2023).  A hybrid approach to comparing parallel-group and stepped-wedge cluster-randomized trials with a continuous primary outcome when there is uncertainty in the intra-cluster correlation.  Clinical Trials, 20(1):59–70.
  * Spiegelhalter, (2001) Spiegelhalter, D. J. (2001).  Bayesian methods for cluster randomized trials with continuous responses.  Statistics in medicine, 20(3):435–452.
  * Thomas et al., (2015) Thomas, L. H., French, B., Sutton, C. J., Forshaw, D., Leathley, M. J., Burton, C. R., Roe, B., Cheater, F., Booth, J., McColl, E., et al. (2015).  Identifying Continence OptioNs after Stroke (ICONS): an evidence synthesis, case study and exploratory cluster randomised controlled trial of the introduction of a systematic voiding programme for patients with urinary incontinence after stroke in secondary care.  Programme Grants for Applied Research, 3(1).
  * Thomas et al., (2014) Thomas, L. H., Watkins, C. L., Sutton, C. J., Forshaw, D., Leathley, M. J., French, B., Burton, C. R., Cheater, F., Roe, B., Britt, D., Booth, J., McColl, E., et al. (2014).  Identifying continence options after stroke (ICONS): a cluster randomised controlled feasibility trial.  Trials, 15(509):1–15.
  * Tishkovskaya et al., (2023) Tishkovskaya, S. V., Sutton, C. J., Thomas, L. H., and Watkins, C. L. (2023).  Determining the sample size for a cluster-randomised trial: Bayesian hierarchical modelling of the intracluster correlation coefficient.  Clinical Trials, 20(3):293–306.
  * Turner et al., (2005) Turner, R. M., Thompson, S. G., and Spiegelhalter, D. J. (2005).  Prior distributions for the intracluster correlation coefficient, based on multiple previous estimates, and their application in cluster randomized trials.  Clinical trials, 2(2):108–118.
  * Turner et al., (2004) Turner, R. M., Toby Prevost, A., and Thompson, S. G. (2004).  Allowing for imprecision of the intracluster correlation coefficient in the design of cluster randomized trials.  Statistics in medicine, 23(8):1195–1214.
  * Ukoumunne, (2002) Ukoumunne, O. C. (2002).  A comparison of confidence interval methods for the intraclass correlation coefficient in cluster randomized trials.  Statistics in medicine, 21(24):3757–3774.
  * Zhan et al., (2021) Zhan, D., Xu, L., Ouyang, Y., Sawatzky, R., and Wong, H. (2021).  Methods for dealing with unequal cluster sizes in cluster randomized trials: A scoping review.  PLOS ONE, 16(7).
  * Zheng et al., (2023) Zheng, H., Jaki, T., and Wason, J. M. (2023).  Bayesian sample size determination using commensurate priors to leverage preexperimental data.  Biometrics, 79(2):669–683.


# SALSA: A Sequential Alternating Least Squares Approximation Method For MIMO
Channel Estimation

Sepideh Gherekhloo, Khaled Ardah, Martin Haardt S. Gherekhloo and M. Haardt
are with the Communications Research Laboratory (CRL), TU Ilmenau, Ilmenau,
Germany (e-mail: {sepideh.gherekhloo, martin.haardt}@tu-ilmenau.de). K. Ardah
is with Lenovo (Deutschland) GmbH (e-mail: kardah@lenovo.com).

###### Abstract

In this paper, we consider the channel estimation problem in sub-6 GHz uplink
wideband MIMO-OFDM communication systems, where a user equipment with a fully-
digital beamforming structure is communicating with a base station having a
hybrid analog-digital beamforming structure. A novel channel estimation method
called Sequential Alternating Least Squares Approximation (SALSA) is proposed
by exploiting a hidden tensor structure in the uplink measurement matrix.
Specifically, by showing that any MIMO channel matrix can be approximately
decomposed into a summation of $R$ factor matrices having a Kronecker
structure, the uplink measurement matrix can be reshaped into a 3-way tensor
admitting a Tucker decomposition. Exploiting the tensor structure, the MIMO
channel matrix is estimated sequentially using an alternating least squares
method. Detailed simulation results are provided showing the effectiveness of
the proposed SALSA method as compared to the classical least squares method.

###### Index Terms:

Channel estimation, massive MIMO, Tucker tensor decomposition, alternating
least squares

##  I Introduction

Massive MIMO [1] is one of the key enabling technologies of 5G-NR mobile
communications [2] and it shall remain relevant in future 6G wireless systems.
By employing a large number of antennas at the base station (BS) relative to
the number of scheduled users, massive MIMO systems increase the data
throughput relative to legacy systems by providing a large beamforming gain
and an improved multi-user interference suppression owing to its high spatial
resolution [3]. Recently, massive MIMO communications have received a special
attention with the introduction of millimeter-wave (mm-wave)-based wireless
communications [4], since the use of massive MIMO in such systems becomes a
requirement rather than an option to compensate the high pathloss encountered
in the wireless communication systems at higher frequencies. However, it is
well-known that the promised theoretical massive MIMO gains heavily rely on
the availability of accurate channel state information (CSI) and the
considered beamforming structure.

On the one hand, classical fully-digital (FD) beamforming structures, which
generally provide the maximum beamforming gain, require a dedicated radio
frequency (RF) chain for each antenna element. This increases not only the
implementation cost and complexity of massive MIMO systems, but also the
circuit energy consumption. A promising solution to these issues relies on the
recently introduced hybrid analog-digital (HAD) beamforming structures [5, 6,
4, 7, 8], which use a combination of analog beamforming in the RF domain and
digital beamforming in the baseband domain to reduce the number of RF chains
as compared to FD beamforming structures, e.g., the number of RF chains can be
as small as the number of transmitted data streams.

On the other hand, in 5G-NR systems, for example, the BS estimates the CSI
from uplink sounding reference signals (SRS) emitted by the user terminals
(UEs). In mm-wave systems, the CSI estimation problem is often transformed
into a multi-dimensional direction-of-arrival (DoA) estimation problem [9, 10,
11], thanks to the low-rank (sparse) nature of mm-wave MIMO channels [4],
where several techniques, e.g., compressed sensing [9, 10] and ESPRIT [11] can
be readily employed to obtain a high CSI estimation accuracy while requiring a
small number of training overhead. Differently, in sub-6 GHz-based systems,
the MIMO channels often experience a high-rank nature, which makes most, if
not all, mm-wave-based MIMO channel estimation methods unfeasible. To this
end, classical channel estimation techniques, e.g., least-squares (LS) and
minimum mean squared-error (MMSE) methods [12, 13] can be used to estimate
sub-6 GHz-based MIMO channels. However, these methods were originally
developed for single-antenna and small-scale MIMO systems and suffer from a
severe performance degradation in difficult scenarios, e.g., with small number
of training snapshots and/or a low signal-to-noise ratio (SNR). Since sub-6
GHz massive MIMO communications are, and will remain, an integral part of
current and future wireless communication systems, more efficient channel
estimation techniques than the classical methods are required.

In this paper, we consider the channel estimation problem in sub-6 GHz uplink
wideband MIMO-OFDM communication systems, where a single-user with a FD
beamforming structure communicates with a BS having a HAD beamforming
structure. By exploiting a hidden tensor structure in the uplink measurement
matrix, we propose a novel channel estimation method called Sequential
Alternating Least Squares Approximation (SALSA). Specifically, by showing that
any MIMO channel matrix can be approximately decomposed into a summation of
$R$ factor matrices having a Kronecker structure, the uplink measurement
matrix can be reshaped into a 3-way tensor admitting a Tucker decomposition
[14]. Exploiting such a tensor representation, the MIMO channel matrix can be
estimated sequentially using the classical ALS method [15]. Detailed
simulation results are provided showing that the SALSA-based approach can
achieve a more accurate channel estimation in difficult scenarios as compared
to the classical LS-based approach.

Notation: The transpose, the complex conjugate, the conjugate transpose
(Hermitian), and the Kronecker product are denoted as ${\bm{A}}^{\mathsf{T}}$,
${\bm{A}}^{*}$, ${\bm{A}}^{\mathsf{H}}$, and $\otimes$, respectively.
Moreover, ${\bm{I}}_{N}$ is the $N\times N$ identity matrix,
$\text{vec}\\{{\bm{A}}\\}$ forms a vector by staking the columns of ${\bm{A}}$
over each other, and the $n$-mode product of a tensor
$\bm{\mathcal{A}}\in\mathbb{C}^{I_{1}\times I_{2}\times\dots,\times I_{N}}$
with a matrix ${\bm{B}}\in\mathbb{C}^{J\times I_{n}}$ is denoted as
$\bm{\mathcal{A}}\times_{n}{\bm{B}}$.

##  II System Model

We consider an uplink single-user wideband MIMO-OFDM communication system, as
depicted in Fig. 1, where a UE with $N_{\text{UE}}$ antennas is communicating
with a BS with $N_{\text{BS}}$ antennas over $N_{\text{SC}}$ subcarriers. The
UE has a FD beamforming structure while the BS has a HAD beamforming structure
with $N_{\text{RF}}\leq N_{\text{BS}}$ radio-frequency (RF) chains. We assume
that the $N_{\text{BS}}$ antennas and the $N_{\text{RF}}$ RF chains are
divided equally111To simplify the exposition, we assume that $N_{\text{BS}}$,
$N_{\text{RF}}$, and $N_{\text{G}}$ are selected so that
$\mathring{N}_{\text{BS}}$ and $\mathring{N}_{\text{RF}}$ are integer numbers,
without loss of generality. into $N_{\text{G}}\geq 1$ groups, where each group
has $\mathring{N}_{\text{BS}}=\frac{N_{\text{BS}}}{N_{\text{G}}}$ antennas and
$\mathring{N}_{\text{RF}}=\frac{N_{\text{RF}}}{N_{\text{G}}}$ RF chains (i.e.,
$N_{\text{BS}}=N_{\text{G}}\cdot\mathring{N}_{\text{BS}}$ and
$N_{\text{RF}}=N_{\text{G}}\cdot\mathring{N}_{\text{RF}}$) and the RF chains
in every group are connected with every antenna element in the same group.
Moreover, we assume a block-fading channel model as shown in Fig. 2, where the
channel coherence-time $T_{\text{C}}$ is divided into
$T_{\text{BS}}T_{\text{UE}}$ transmission time intervals (TTIs), i.e., every
block has $T_{\text{UE}}$ snapshots.

Let $\bm{\bar{A}}_{i}\in\mathbb{C}^{N_{\text{BS}}\times N_{\text{RF}}}$ denote
the analog combining matrix at the $i$th block at the BS. Then, according to
our above assumptions, $\bm{\bar{A}}_{i}$ has a block-diagonal structure given
as222Note that if $N_{\text{G}}=1$, the above analog structure coincides with
the known fully-connected analog structure [5], where every RF chain is
connected to every antenna element. On the other hand, if
$N_{\text{G}}=N_{\text{RF}}$, the above analog structure coincides with the
known partially-connected analog structure [5], where every RF chain is
connected to a unique subset of antenna elements.

$\displaystyle\bm{\bar{A}}_{i}=\frac{1}{\sqrt{\mathring{N}_{\text{BS}}}}\cdot\begin{bmatrix}\bm{\bar{A}}_{i,1}&\dots&{\bm{0}}\\\
\vdots&\ddots&\vdots\\\
{\bm{0}}&\dots&\bm{\bar{A}}_{i,N_{\text{G}}}\end{bmatrix}\in\mathbb{C}^{N_{\text{BS}}\times
N_{\text{RF}}},$ (1)

where
$\bm{\bar{A}}_{i,g}\in\mathbb{C}^{\mathring{N}_{\text{BS}}\times\mathring{N}_{\text{RF}}}$
is the $g$th block-matrix with constant modulus entries, i.e.,
$\big{|}[\bm{\bar{A}}_{i,g}]_{[r,c]}\big{|}=1$, where
$[\bm{\bar{A}}_{i,g}]_{[r,c]}$ is the $(r,c)$th entry of $\bm{\bar{A}}_{i,g}$.

The received signal by the BS in the $(i,j)$th TTI over the $k$th subcarrier,
with $i\in\\{1,\dots,T_{\text{BS}}\\}$, $j\in\\{1,\dots,T_{\text{UE}}\\}$,
$k\in\\{1,\dots,N_{\text{SC}}\\}$, can be expressed as

$\displaystyle\bm{\bar{y}}_{k,i,j}=\bm{\bar{A}}^{\mathsf{H}}_{i}{\bm{H}}_{k}{\bm{f}}_{k,j}s_{k,j}+\bm{\bar{A}}^{\mathsf{H}}_{i}\bm{\bar{z}}_{k,i,j}\in\mathbb{C}^{N_{\text{RF}}},$
(2)

where ${\bm{f}}_{k,j}\in\mathbb{C}^{{N}_{\text{UE}}}$ is the $(k,j)$th
precoding vector, $s_{k,j}\in\mathbb{C}$ is the corresponding training symbol,
$\bm{\bar{z}}_{i,j}\in\mathbb{C}^{{N}_{\text{BS}}}$ is the BS additive white
Gaussian noise with zero mean and variance $\sigma_{n}^{2}$, and
${\bm{H}}_{k}\in\mathbb{C}^{{N}_{\text{BS}}\times{N}_{\text{UE}}}$ is the
$k$th subcarrier frequency-domain MIMO channel matrix.

Initially, we collect the measurement vectors
$\\{\bm{\bar{y}}_{k,i,j}\\}_{j=1}^{T_{\text{UE}}}$ next to each other as
$\bm{\bar{Y}}_{k,i}=[\bm{\bar{y}}_{k,i,1},\dots,\bm{\bar{y}}_{k,i,T_{\text{UE}}}]$,
which can be written as

$\displaystyle\bm{\bar{Y}}_{k,i}=\bm{\bar{A}}^{\mathsf{H}}_{i}{\bm{H}}_{k}{\bm{F}}_{k}+\bm{\bar{A}}^{\mathsf{H}}_{i}\bm{\bar{Z}}_{k,i}\in\mathbb{C}^{{N}_{\text{RF}}\times
T_{\text{UE}}},$ (3)

where
${\bm{F}}_{k}=[{\bm{f}}_{k,1}s_{k,1},\dots,{\bm{f}}_{k,T_{\text{UE}}}s_{k,T_{\text{UE}}}]\in\mathbb{C}^{{N}_{\text{UE}}\times
T_{\text{UE}}}$ and
$\bm{\bar{Z}}_{k,i}=[\bm{\bar{z}}_{k,i,1},\dots,\bm{\bar{z}}_{k,i,T_{\text{UE}}}]$.
We assume that ${\bm{F}}_{k},\forall k$, are designed with orthonormal rows,
i.e., ${\bm{F}}_{k}{\bm{F}}^{\mathsf{H}}_{k}={\bm{I}}_{N_{\text{UE}}},\forall
k$, and $T_{\text{UE}}\geq N_{\text{UE}}$. After applying the right-filtering
to (3) we obtain

$\displaystyle{\bm{Y}}_{k,i}=\bm{\bar{Y}}_{k,i}{\bm{F}}^{\mathsf{H}}=\bm{\bar{A}}^{\mathsf{H}}_{i}{\bm{H}}_{k}+{\bm{Z}}_{k,i}\in\mathbb{C}^{N_{\text{RF}}\times
N_{\text{UE}}},$ (4)

where
${\bm{Z}}_{k,i}=\bm{\bar{A}}^{\mathsf{H}}_{i}\bm{\bar{Z}}_{k,i}{\bm{F}}^{\mathsf{H}}_{k}$.
Next, we collect the measurement matrices
$\\{{\bm{Y}}_{k,i}\\}_{i=1}^{T_{\text{BS}}}$ on the top of each other as
${\bm{Y}}_{k}=\big{[}{\bm{Y}}^{\mathsf{T}}_{k,1},\dots,{\bm{Y}}^{\mathsf{T}}_{k,T_{\text{BS}}}\big{]}^{\mathsf{T}}$,
which can be written as

$\displaystyle{\bm{Y}}_{k}={\bm{A}}{\bm{H}}_{k}+{\bm{Z}}_{k}\in\mathbb{C}^{L\times
N_{\text{UE}}},$ (5)

where $L=T_{\text{BS}}{N}_{\text{RF}}$,
${\bm{A}}=\big{[}\bm{\bar{A}}_{1},\dots,\bm{\bar{A}}_{T_{\text{BS}}}\big{]}^{\mathsf{H}}\in\mathbb{C}^{L\times
N_{\text{BS}}}$, and
${\bm{Z}}_{k}=\big{[}{\bm{Z}}^{\mathsf{T}}_{k,1},\dots,{\bm{Z}}^{\mathsf{T}}_{k,T_{\text{BS}}}\big{]}^{\mathsf{T}}$.
After that, we collect the measurement matrices
$\\{{\bm{Y}}_{k}\\}_{k=1}^{N_{\text{SC}}}$ next to each other as
$\bm{{Y}}=[{\bm{Y}}_{1},\dots,{\bm{Y}}_{N_{\text{SC}}}]$, which can be written
as

$\displaystyle{\bm{Y}}={\bm{A}}{\bm{H}}+{\bm{Z}}\in\mathbb{C}^{L\times
N_{\text{UE}}N_{\text{SC}}},$ (6)

where ${\bm{Z}}=\big{[}{\bm{Z}}_{1},\dots,{\bm{Z}}_{N_{\text{SC}}}\big{]}$ and
${\bm{H}}=\big{[}\bm{H}_{1},\dots,\bm{H}_{N_{\text{SC}}}\big{]}\in\mathbb{C}^{N_{\text{BS}}\times
N_{\text{UE}}N_{\text{SC}}}$ is the total MIMO channel matrix.

Figure 1: The considered uplink MIMO-OFDM communication system. Figure 2: The
channel coherence time $T_{\text{C}}$ division.

The baseline LS-based channel estimation method: Given the measurement matrix
in (6), a least-squares (LS)-based method can be used to obtain an estimate of
the total MIMO channel matrix as

$\displaystyle\bm{\hat{H}}_{\text{LS}}=[{\bm{A}}]^{+}{\bm{Y}}=\big{[}\bm{\hat{H}}_{1},\dots,\bm{\hat{H}}_{N_{\text{SC}}}\big{]}\in\mathbb{C}^{N_{\text{BS}}\times
N_{\text{UE}}N_{\text{SC}}},$ (7)

where $[\cdot]^{+}$ denotes the Moore-Penrose pseudo-inverse. Note that, due
to the left filtering, the LS-based method requires that $L\geq
N_{\text{BS}}$, i.e., $T_{\text{BS}}\geq\frac{N_{\text{BS}}}{N_{\text{RF}}}$
to provide an accurate channel estimate.

##  III The proposed SALSA method

To obtain a more accurate channel estimate while reducing the training
overhead, we propose in this section a novel channel estimation method called
SALSA, which is derived by exploiting a hidden tensor structure in the
measurement matrix in (6). To show this, we first recall the following
propositions from [16, 17, 18].

Proposition 1: Let ${\bm{X}}$ be a matrix given as

$\displaystyle{\bm{X}}={\bm{X}}_{1}\otimes{\bm{X}}_{2}=\begin{bmatrix}{\bm{X}}_{1,1}&\dots&{\bm{X}}_{1,J_{1}}\\\
\dots\\\ {\bm{X}}_{I_{1},1}&\dots&{\bm{X}}_{I_{1},J_{1}}\\\
\end{bmatrix}\in\mathbb{C}^{I\times J},$ (8)

where ${\bm{X}}_{1}\in\mathbb{C}^{I_{1}\times J_{1}}$,
${\bm{X}}_{2}\in\mathbb{C}^{I_{2}\times J_{2}}$, $I=I_{1}I_{2}$,
$J=J_{1}J_{2}$, and ${\bm{X}}_{n,m}=[{\bm{X}}_{1}]_{[n,m]}{\bm{X}}_{2}$ is the
$(n,m)$th block-matrix of ${\bm{X}}$. Let
${\bm{K}}\in\mathbb{C}^{I_{1}J_{1}\times I_{2}J_{2}}$ be a rank-one matrix
given as

$\displaystyle{\bm{K}}=\begin{bmatrix}\text{vec}\\{{\bm{X}}_{1,1}\\}^{\text{T}}\\\
\dots\\\ \text{vec}\\{{\bm{X}}_{I_{1},1}\\}^{\text{T}}\\\ \dots\\\
\text{vec}\\{{\bm{X}}_{1,J_{1}}\\}^{\text{T}}\\\ \dots\\\
\text{vec}\\{{\bm{X}}_{I_{1},J_{1}}\\}^{\text{T}}\\\
\end{bmatrix}=\text{vec}\\{{\bm{X}}_{1}\\}\text{vec}\\{{\bm{X}}_{2}\\}^{\text{T}},$
(9)

with the rank-one truncated-SVD given as
${\bm{K}}=\sigma{\bm{u}}{\bm{v}}^{\text{H}}$, where
${\bm{u}}\in\mathbb{C}^{I_{1}J_{1}}$ and ${\bm{v}}\in\mathbb{C}^{I_{2}J_{2}}$
are the left and right singular vectors of ${\bm{K}}$, respectively, and
$\sigma$ is the associated singular value. Then, the optimal solution to

$\displaystyle\underset{{\bm{X}}_{1},{\bm{X}}_{2}}{\text{minimize
}}\|{\bm{X}}-\big{(}{\bm{X}}_{1}\otimes{\bm{X}}_{2}\big{)}\|^{2}_{\text{F}}$
(10)

can be obtained as

$\displaystyle{\bm{X}}_{1}$
$\displaystyle=\text{reshape}\\{\sqrt{\sigma}{\bm{u}},I_{1},J_{1}\\}$ (11)
$\displaystyle{\bm{X}}_{2}$
$\displaystyle=\text{reshape}\\{\sqrt{\sigma}{\bm{v}}^{*},I_{2},J_{2}\\}.$
(12)

Proof: Please refer to [17] for more details.

Proposition 2: For any given $I\times J$ matrix ${\bm{X}}$, it can be
approximately written as a summation of $R\geq 1$ factor matrices as

$\displaystyle{\bm{X}}=\sum_{r=1}^{R}{\bm{X}}_{r}=\sum_{r=1}^{R}{\bm{X}}_{1,r}\otimes{\bm{X}}_{2,r},$
(13)

where ${\bm{X}}_{r}={\bm{X}}_{1,r}\otimes{\bm{X}}_{2,r}$,
${\bm{X}}_{1,r}\in\mathbb{C}^{I_{1}\times J_{1}}$, and
${\bm{X}}_{2,r}\in\mathbb{C}^{I_{2}\times J_{2}}$, $I=I_{1}I_{2}$, and
$J=J_{1}J_{2}$.

Proof: The proof follows directly by applying Proposition 1 sequentially [18].
The corresponding Proposition is summarized in Algorithm 1.

Let $I={N}_{\text{BS}}$ and $J={N}_{\text{UE}}N_{\text{SC}}$. Then, from
Proposition 2, the total frequency-domain MIMO channel matrix
${\bm{H}}\in\mathbb{C}^{I\times J}$ in (6) can be approximately written as

$\displaystyle{\bm{H}}\approx\sum_{r=1}^{R}{\bm{C}}_{r}\otimes{\bm{B}}_{r}\in\mathbb{C}^{I\times
J},$ (14)

where ${\bm{B}}_{r}\in\mathbb{C}^{I_{1}\times J_{1}}$,
${\bm{C}}_{r}\in\mathbb{C}^{I_{2}\times J_{2}}$, $I=I_{1}I_{2}$, and
$J=J_{1}J_{2}$. As shown in Fig. 3, the approximation becomes tighter as the
number of channel factor matrices $R$ increases. More importantly, we can see
that in case of full rank channels, the optimal value of $R$, denoted in the
figure by $R_{\text{opt}}$, is dependent on the division scenario of $I$ and
$J$, where $R_{\text{opt}}\approx\min\\{I_{1}J_{1},I_{2}J_{2}\\}$. In other
words, reducing the dimension of one of the channel factor matrices, i.e.,
${\bm{B}}_{r}\in\mathbb{C}^{I_{1}\times J_{1}}$ or
${\bm{C}}_{r}\in\mathbb{C}^{I_{2}\times J_{2}}$, reduces the value of
$R_{\text{opt}}$.

Algorithm 1 Sequential Kronecker Factorization

1:Input: A matrix ${\bm{X}}\in\mathbb{C}^{I\times J}$

2:Select $R$, $I_{1}$, $J_{1}$, $I_{2}$, $J_{2}$ such that $I=I_{1}I_{2}$ and
$J=J_{1}J_{2}$

3:for $r=1$ to $R$ do

4: Get
${\bm{X}}_{r}={\bm{X}}-\sum_{r^{\prime}=1}^{r-1}{\bm{X}}_{1,r^{\prime}}\otimes{\bm{X}}_{2,r^{\prime}}$

5: Given ${\bm{X}}_{r}$, get ${\bm{X}}_{1,r}$ and ${\bm{X}}_{2,r}$ using (11)
and (12), respectively

6:end for

7:Output:
${\bm{\hat{}}X}=\sum_{r=1}^{R}{\bm{X}}_{1,r}\otimes{\bm{X}}_{2,r}\in\mathbb{C}^{I\times
J}$

Figure 3: MSE vs. the number of channel factor matrices $R$ assuming
$N_{\text{BS}}=64$, $N_{\text{UE}}=4$, and $N_{\text{SC}}=16$, where
$\text{MSE}=\|{\bm{H}}-\sum_{r=1}^{R}{\bm{C}}_{r}\otimes{\bm{B}}_{r}\|^{2}_{\text{F}}$.
Here, the total MIMO channel matrix
${\bm{H}}=\big{[}\bm{H}_{1},\dots,\bm{H}_{N_{\text{SC}}}\big{]}\in\mathbb{C}^{I\times
J}$ is generated following the 3GPP CDL channel model [19, 20] with the main
system parameters outlined in Table I. Please refer to Section IV for more
details.

Let $L=T_{\text{BS}}{N}_{\text{RF}}$. Then, by substituting (14) into (6), and
assuming $R$ is sufficiently large, we can write

$\displaystyle{\bm{Y}}$
$\displaystyle={\bm{A}}\Bigg{(}\sum_{r=1}^{R}{\bm{C}}_{r}\otimes{\bm{B}}_{r}\Bigg{)}+{\bm{Z}}=\sum_{r=1}^{R}{\bm{A}}({\bm{C}}_{r}\otimes{\bm{B}}_{r})+{\bm{Z}}$
$\displaystyle=\sum_{r=1}^{R}{\bm{Y}}_{r}+{\bm{Z}}\in\mathbb{C}^{L\times J},$
(15)

where
${\bm{Y}}_{r}={\bm{A}}({\bm{C}}_{r}\otimes{\bm{B}}_{r})\in\mathbb{C}^{L\times
J}$. From (III), we note that ${\bm{Y}}_{r}$ can be seen as the 1-mode
unfolding of a 3-way Tucker tensor given as [14]

$\displaystyle\bm{\mathcal{Y}}_{r}=\bm{\mathcal{S}}\times_{1}{\bm{A}}\times_{2}{\bm{B}}^{\mathsf{T}}_{r}\times_{3}{\bm{C}}^{\mathsf{T}}_{r}\in\mathbb{C}^{L\times
J_{1}\times J_{2}},$ (16)

where $\bm{\mathcal{S}}\in\mathbb{Z}^{I\times I_{1}\times I_{2}}$ is the core-
tensor with the 1-mode unfolding given as
$[\bm{\mathcal{S}}]_{(1)}\overset{\text{def}}{=}{\bm{I}}_{I}$. The $\ell$-mode
unfolding of $\bm{\mathcal{Y}}_{r}$, $\ell=\\{1,2,3\\}$, can be expressed as

$\displaystyle[\bm{\mathcal{Y}}_{r}]_{(1)}$
$\displaystyle={\bm{A}}[\bm{\mathcal{S}}]_{(1)}({\bm{C}}_{r}\otimes{\bm{B}}_{r})\in\mathbb{C}^{L\times
J},$ (17) $\displaystyle[\bm{\mathcal{Y}}_{r}]_{(2)}$
$\displaystyle={\bm{B}}^{\mathsf{T}}_{r}[\bm{\mathcal{S}}]_{(2)}({\bm{C}}_{r}\otimes{\bm{A}}^{\mathsf{T}})\in\mathbb{C}^{J_{1}\times
LJ_{2}},$ (18) $\displaystyle[\bm{\mathcal{Y}}_{r}]_{(3)}$
$\displaystyle={\bm{C}}^{\mathsf{T}}_{r}[\bm{\mathcal{S}}]_{(3)}({\bm{B}}_{r}\otimes{\bm{A}}^{\mathsf{T}})\in\mathbb{C}^{J_{2}\times
LJ_{1}}.$ (19)

From (16), the 3-way Tucker tensor form of (III) can be expressed as

$\displaystyle\bm{\mathcal{Y}}=\sum_{r=1}^{R}\bm{\mathcal{Y}}_{r}+\bm{\mathcal{Z}}\in\mathbb{C}^{L\times
J_{1}\times J_{2}},$ (20)

where $\bm{{{\mathcal{Z}}}}$ is the 3-way tensor representation of the noise
matrix ${\bm{Z}}$. This latter formulation suggests that the factor matrices
$\\{{\bm{B}}_{r},{\bm{C}}_{r}\\}_{r=1}^{R}$ can be estimated sequentially as
follows. Let $\bm{\mathcal{Y}}_{r}$ be the tensor obtained at the $r$th
sequential step as

$\displaystyle\bm{\mathcal{Y}}_{r}=\bm{\mathcal{Y}}-\sum_{r^{\prime}=1}^{r-1}\bm{\mathcal{Y}}_{r^{\prime}}\in\mathbb{C}^{L\times
J_{1}\times J_{2}}.$ (21)

Then, by exploiting the 2-mode and the 3-mode unfoldings, the $r$th factor
matrices ${\bm{B}}_{r}$ and ${\bm{C}}_{r}$ can be obtained using, e.g., the
ALS method [15], where one factor matrix is assumed to be fixed when solving
for the other. Specifically, ${\bm{B}}_{r}$ and ${\bm{C}}_{r}$ can be obtained
as

$\displaystyle{\bm{B}}^{\mathsf{T}}_{r}$
$\displaystyle=[\bm{\mathcal{Y}}_{r}]_{(2)}\big{[}\bm{\Psi}_{2}\big{]}^{+}=[\bm{\mathcal{Y}}_{r}]_{(2)}{\bm{\Psi}}_{2}^{\text{H}}[{\bm{\Psi}}_{2}{\bm{\Psi}}_{2}^{\text{H}}]^{-1}$
(22) $\displaystyle{\bm{C}}^{\mathsf{T}}_{r}$
$\displaystyle=[\bm{\mathcal{Y}}_{r}]_{(3)}\big{[}{\bm{\Psi}}_{3}\big{]}^{+}=[\bm{\mathcal{Y}}_{r}]_{(3)}{\bm{\Psi}}_{3}^{\text{H}}[{\bm{\Psi}}_{3}{\bm{\Psi}}_{3}^{\text{H}}]^{-1},$
(23)

where ${\bm{\Psi}}_{2}$ and ${\bm{\Psi}}_{3}$ are given as

$\displaystyle{\bm{\Psi}}_{2}$
$\displaystyle=[\bm{\mathcal{S}}]_{(2)}({\bm{C}}_{r}\otimes{\bm{A}}^{\mathsf{T}})\in\mathbb{C}^{I_{1}\times
LJ_{2}}$ (24) $\displaystyle{\bm{\Psi}}_{3}$
$\displaystyle=[\bm{\mathcal{S}}]_{(3)}({\bm{B}}_{r}\otimes{\bm{A}}^{\mathsf{T}})\in\mathbb{C}^{I_{2}\times
LJ_{1}}.$ (25)

Algorithm 2 summarizes the proposed SALSA method for estimating the total MIMO
channel matrix ${\bm{H}}\in\mathbb{C}^{I\times J}$, which is guaranteed to
converge monotonically to, at least, a local optimum solution [15].

Algorithm 2 SALSA For MIMO-OFDM Channel Estimation

1:Input: Measurement matrix ${\bm{Y}}\in\mathbb{C}^{L\times J}$ as in (6)

2:Select $R\geq 1$, $N_{\text{max-iter}}\geq 1$, $I_{1},I_{2},J_{1}$, and
$J_{2}$ such that $I=I_{1}I_{2}={N}_{\text{BS}}$ and
$J=J_{1}J_{2}={N}_{\text{UE}}N_{\text{SC}}$

3:Obtain the 3-way Tucker tensor $\bm{\mathcal{Y}}$ in (20) from ${\bm{Y}}$

4:for $r=1$ to $R$ do

5: Get
$\bm{\mathcal{Y}}_{r}=\bm{\mathcal{Y}}-\sum_{r^{\prime}=1}^{r-1}\bm{\mathcal{\hat{Y}}}_{r^{\prime}}$

6: Initialize ${\bm{C}}^{(0)}_{r}\in\mathbb{C}^{I_{2}\times J_{2}}$, e.g.,
randomly

7: for $n=1$ to $N_{\text{max-iter}}$ do

8: Get ${\bm{B}}^{(n)}_{r}$ using (22) for given ${\bm{C}}^{(n-1)}_{r}$

9: Get ${\bm{C}}^{(n)}_{r}$ using (23) for given ${\bm{B}}^{(n)}_{r}$

10: end for

11: Set $\bm{\hat{B}}_{r}={\bm{B}}^{(N_{\text{max-iter}})}_{r}$ and
$\bm{\hat{C}}_{r}={\bm{C}}^{(N_{\text{max-iter}})}_{r}$

12: Get
$\bm{\mathcal{\hat{Y}}}_{r}=\bm{\mathcal{S}}\times_{1}{\bm{A}}\times_{2}\bm{\hat{B}}^{\mathsf{T}}_{r}\times_{3}\bm{\hat{C}}^{\mathsf{T}}_{r}$,
go back to Step (5)

13:end for

14:Output:
$\bm{\hat{H}}_{\text{SALSA}}=\sum_{r=1}^{R}\bm{\hat{C}}_{r}\otimes\bm{\hat{B}}_{r}\in\mathbb{C}^{I\times
J}$

Note that, due to the right filtering, the SALSA method in Algorithm 2
requires that (C1) $I_{1}\leq LJ_{2}$ and (C2) $I_{2}\leq LJ_{1}$, i.e.,
$T_{\text{BS}}\geq\min\Big{\\{}\frac{I_{1}}{{N}_{\text{RF}}J_{2}},\frac{I_{2}}{{N}_{\text{RF}}J_{1}}\Big{\\}}$
to provide an accurate channel estimation. Therefore, under practical
settings, the SALSA method in Algorithm 2 requires less training overhead than
the LS method in (7). On the other hand, assuming that the complexity of
calculating the Moore-Penrose pseudo-inverse of an $n\times m$ matrix is on
the order of $\mathcal{O}(\min\\{n,m\\}^{3})$, then the complexity of the LS
method in (7) is on the order of $\mathcal{O}(\min\\{L,J\\})^{3}$, while for
the SALSA method in Algorithm 2 the complexity is on the order of
$\mathcal{O}(R\cdot N_{\text{max-iter}}\cdot I_{1}^{3}\cdot I_{2}^{3})$,
assuming that the (C1) and (C2) conditions are satisfied.

##  IV Simulation Results

We adopt the 3GPP clustered delay line (CDL) channel model described in TR
38.901 [19], where a step-by-step tutorial of it along the MATLAB scripts for
channel generation is presented in [20]. Specifically, in our simulation, we
first generate a time-domain channel tensor
$\bm{\mathcal{H}}\in\mathbb{C}^{N_{\text{BS}}\times N_{\text{UE}}\times
N_{\text{taps}}}$, where $N_{\text{taps}}$ represents the number of time-
domain channel taps calculated according to [20, Eqn. (64)] and using the
system parameters shown in Table I. Then, we perform a $N_{\text{SC}}$-point
FFT operation along the third dimension for each receive-transmit antenna pair
to obtain the frequency-domain channel tensor
$\bm{\mathcal{H}}\in\mathbb{C}^{N_{\text{BS}}\times N_{\text{UE}}\times
N_{\text{SC}}}$, where the $k$th slice matrix
$\bm{H}_{k}=\bm{\mathcal{H}}_{[:,:,k]}\in\mathbb{C}^{N_{\text{BS}}\times
N_{\text{UE}}}$ represents the the $k$th subcarrier frequency-domain MIMO
channel matrix.

We show the simulation results in terms of the normalized mean-square-error
(NMSE) that is defined as
$\text{NMSE}=\mathbb{E}\\{\|{\bm{H}}-\bm{\hat{H}}_{\text{X}}\|_{\text{F}}^{2}\\}/\mathbb{E}\|{\bm{H}}\|_{\text{F}}^{2}\\}$,
where $\text{X}\in\\{\text{LS, SALSA}\\}$. The signal-to-noise ratio (SNR) is
defined as
$\text{SNR}=\mathbb{E}\\{\|\bm{\mathcal{Y}}-\bm{\mathcal{Z}}\|_{\text{F}}^{2}\\}/\mathbb{E}\\{\|\bm{\mathcal{Z}}\|_{\text{F}}^{2}\\}$.
In all simulation scenarios, we set $N_{\text{BS}}=64$, $N_{\text{UE}}=4$,
$N_{\text{SC}}=16$, $T_{\text{UE}}=N_{\text{UE}}$, $N_{\text{RF}}=4$,
$N_{\text{G}}=2$, and assume a random generation of the analog decoding matrix
$\bm{A}\in\mathbb{C}^{T_{\text{BS}}{N}_{\text{RF}}\times N_{\text{BS}}}$,
where every nonzero entry is obtained as
$a={1}/{\sqrt{\mathring{N}_{\text{BS}}}}\cdot e^{j\phi}$, where
$\phi\in[0,2\pi]$.

TABLE I: System Parameters Parameter | Value
---|---
Scenario | UMi
Cell radius | 100 m
BS (UE) height | 10 (1.5) m
Carrier frequency $f_{c}$ | 4 GHz
Sampling frequency $f_{s}$ | 30.72 MSamples/s
No. of subcarriers $N_{\text{SC}}$ | $16$
No. of antennas at BS $N_{\text{BS}}$ | 64 $(8\times 8)$
No. of antennas at UE $N_{\text{UE}}$ | 4 $(2\times 2)$
Polarization | Single

Figure 4: NMSE vs. SNR for different $I$ and $J$ division scenarios.

Initially, we show simulation results investigating the best division scenario
of $I$ and $J$ with the constraints of $I=I_{1}I_{2}$, $J=J_{1}J_{2}$,
$I_{d}\geq 1$, $J_{d}\geq 1$, and $I_{d},J_{d}$ are Natural numbers, where
$d\in\\{1,2\\}$. Recall that $I={N}_{\text{BS}}$ and
$J={N}_{\text{UE}}N_{\text{SC}}$. Therefore, we have $I=J=64$ and the
candidate numbers of $I_{d}$ and $J_{d}$ are $1,2,4,8,16,32$, and $64$.
Therefore, we have in total 49 different division scenarios as illustrated in
Table II. We have simulated the SALSA algorithm using all the 49 possible
scenarios. In Fig. 4, we show the NMSE versus SNR results for some selected
$I$ and $J$ division scenarios. The other scenarios are not shown, due to
space limitations, but we note that their NMSE performance are inferior
compared to the shown scenarios.

From Fig. 4, when $T_{\text{BS}}=12$, i.e.,
$L=T_{\text{BS}}{N}_{\text{RF}}=48<N_{\text{BS}}$, the analog training matrix
$\bm{A}\in\mathbb{C}^{L\times N_{\text{BS}}}$, i.e., the 1st factor matrix of
the measurement tensor in (16), is left non-invertible, i.e.,
$[{\bm{A}}]^{+}{\bm{A}}\neq\bm{I}$. Therefore, the LS-based method has a very
bad channel estimation accuracy NMSE. On the other hand, we can see that the
best NMSE of SALSA method is achieved when $I_{1}=8$, $I_{2}=8$, $J_{1}=64$,
and $J_{2}=1$, i.e., when ${\bm{B}}_{r}\in\mathbb{C}^{8\times 64}$ and
${\bm{C}}_{r}\in\mathbb{C}^{8\times 1},\forall r$. The main reason is that by
dividing $I=64$ equally between $I_{1}$ and $I_{2}$, i.e., $I_{1}=I_{2}=8$,
SALSA reduces the impact of the non-invertibility of $\bm{A}$ by distributing
it between the 2nd (i.e., ${\bm{B}}_{r}$) and the 3rd (i.e., ${\bm{C}}_{r}$)
factor matrices of the measurement tensor, which leads to a better channel
estimation accuracy. On the other hand, by setting $J_{1}=64$ and $J_{2}=1$,
the required number of channel factor matrices $R$ reduces as compared to the
other division scenario, as we have illustrated above in Fig. 3.

Differently, when $T_{\text{BS}}=16$, i.e., $L=N_{\text{BS}}$, the analog
training matrix $\bm{A}$ is left invertible, i.e.,
$[{\bm{A}}]^{+}{\bm{A}}=\bm{I}$. Therefore, the LS-based method has an
accurate channel estimation accuracy. For SALSA method, on the other hand, we
can see that when $I_{1}=8$ and $I_{2}=8$, the estimation accuracy of SALSA
improves as we increase $J_{1}$ and decrease $J_{2}$, where the best result is
obtained when we have $J_{1}=64$ and $J_{2}=1$, i.e., similar to the case
above when $T_{\text{BS}}=12$. Nonetheless, we can see that the SALSA method
can obtain a more accurate channel estimation, compered to the LS-based
method, by setting $I_{1}=1$, $I_{2}=64$, $J_{1}=64$, and $J_{2}=1$, i.e.,
${\bm{B}}_{r}\in\mathbb{C}^{1\times 64}$ and
${\bm{C}}_{r}\in\mathbb{C}^{64\times 1}$ (or, not shown in the figure, by
setting $I_{1}=64$, $I_{2}=1$, $J_{1}=1$, and $J_{2}=64$, i.e.,
${\bm{B}}_{r}\in\mathbb{C}^{64\times 1}$ and
${\bm{C}}_{r}\in\mathbb{C}^{1\times 64}$). In the both these scenarios, the
channel matrix $\bm{H}\in\mathbb{C}^{64\times 64}$ in (14) is decomposed into
a summation of $R$ factor matrices
${\bm{B}}_{r}\otimes{\bm{C}}_{r}\in\mathbb{C}^{64\times 64}$, each having a
rank-one, i.e., rank$\\{{\bm{B}}_{r}\otimes{\bm{C}}_{r}\\}=1,\forall r$, which
leads to a better channel estimation accuracy.

TABLE II: Division scenarios of $I_{1},I_{2},J_{1}$, and $J_{2}$ Scenario No. |  $I_{1}$ and $I_{2}$ values |  $J_{1}$ and $J_{2}$ values
---|---|---
Scenario 1 | $[I_{1},I_{2}]=[64,1]$ | $[J_{1},J_{2}]=[64,1]$
| $\vdots$
Scenario 7 | $[J_{1},J_{2}]=[1,64]$
$\vdots$ | $\vdots$ | $\vdots$
Scenario 43 | $[I_{1},I_{2}]=[1,64]$ | $[J_{1},J_{2}]=[64,1]$
| $\vdots$
Scenario 49 | $[J_{1},J_{2}]=[1,64]$

In Figs. 5 and 6 we show NMSE versus SNR simulation results with varying the
number of channel training overhead, i.e., $T_{\text{BS}}$ and the number of
channel factor matrices, i.e., $R$, respectively. From Fig. 5, we can see that
the channel estimation accuracy of both methods, i.e., LS-based and SALSA
improves as $T_{\text{BS}}$ increases. However, SALSA significantly
outperforms LS-based with all $T_{\text{BS}}<16$, i.e., $L<N_{\text{BS}}$
scenarios, wherein the analog training matrix $\bm{A}\in\mathbb{C}^{L\times
N_{\text{BS}}}$ is left non-invertible.

On the other hand, we can see from Fig. 6 that the SALSA channel estimation
accuracy increases with the increasing $R$, in the high SNR regime, while it
decreases with the increasing $R$, in the low SNR regime. The main reason is
that, in the high SNR regime, the noise impact is minimal and by increasing
$R$, the channel estimation accuracy increases, as we have illustrated above
in Fig. 3. On the other hand, in the low SNR regime, the channel measurement
tensor is noise-limited and, therefore, the impact of noise increases by
increasing $R$, i.e., after a certain $R$, the estimated channel factor
matrices are very noisy that decreases the overall estimation accuracy.
Clearly, for every SNR regime/level, there is an optimal $R$ value, wherein
the channel estimation accuracy is maximized, which we leave for a follow up
future work.

Figure 5: NMSE vs. SNR with varying $T_{\text{BS}}$.
Figure 6: NMSE vs. SNR with varying $R$.

##  V Conclusion

In this paper, we have proposed a novel channel estimation method for MIMO-
OFDM sub-6 GHz communication systems called SALSA. We have shown that an
accurate channel estimation can be obtained with a small training overhead by
exploiting a hidden tensor structure in the received measurement matrix, which
estimates the channel matrix sequentially using an ALS-based method. Our
results show that the SALSA method outperforms the conventional LS-based
method, especially in the low training overhead, which makes it more appealing
for practical implementations.

## References

  * [1] T. L. Marzetta, E. G. Larsson, H. Yang, and H. Q. Ngo, _Fundamentals of Massive MIMO_. Cambridge University Press, 2016.
  * [2] E. Dahlman, S. Parkvall, and J. Skold, _5G NR: The Next Generation Wireless Access Technology_ , 1st ed. USA: Academic Press, Inc., 2018.
  * [3] E. Björnson, L. Sanguinetti, J. Hoydis, and M. Debbah, “Optimal design of energy-efficient multi-user MIMO systems: Is massive MIMO the answer?” _IEEE Trans. Wireless Commun._ , vol. 14, no. 6, pp. 3059–3075, Jun. 2015\.
  * [4] R. W. Heath, N. González-Prelcic, S. Rangan, W. Roh, and A. M. Sayeed, “An overview of signal processing techniques for Millimeter Wave MIMO systems,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 10, no. 3, pp. 436–453, 2016.
  * [5] K. Ardah, G. Fodor, Y. C. B. Silva, W. Cruz, and F. R. Cavalcanti, “A unifying design of hybrid beamforming architectures employing phase-shifters or switches,” _IEEE Trans. Veh. Technol._ , pp. 1–1, 2018.
  * [6] X. Gao, L. Dai, S. Han, C.-L. I, and R. W. Heath, “Energy-efficient hybrid analog and digital precoding for MmWave MIMO systems with large antenna arrays,” _IEEE Journal on Selected Areas in Communications_ , vol. 34, no. 4, pp. 998–1009, 2016.
  * [7] H. Shokri-Ghadikolaei, C. Fischione, G. Fodor, P. Popovski, and M. Zorzi, “Millimeter wave cellular networks: A MAC layer perspective,” _IEEE Transactions on Communications_ , vol. 63, no. 10, pp. 3437–3458, 2015.
  * [8] S. Gherekhloo, K. Ardah, and M. Haardt, “Hybrid beamforming design for downlink MU-MIMO-OFDM millimeter-wave systems,” in _Proc. IEEE 11th Sensor Array and Multichannel Signal Processing Workshop (SAM)_ , 2020, pp. 1–5.
  * [9] A. Alkhateeb, G. Leus, and R. W. Heath, “Compressed sensing based multi-user millimeter wave systems: How many measurements are needed?” in _Proc. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2015, pp. 2909–2913.
  * [10] K. Ardah, B. Sokal, A. L. F. de Almeida, and M. Haardt, “Compressed sensing based channel estimation and open-loop training design for hybrid analog-digital massive MIMO systems,” in _Proc. 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2020, pp. 4597–4601.
  * [11] J. Zhang and M. Haardt, “Channel estimation and training design for hybrid multi-carrier mmwave massive MIMO systems: The beamspace ESPRIT approach,” in _Proc. 25th European Signal Processing Conference (EUSIPCO)_ , 2017, pp. 385–389.
  * [12] I. Barhumi, G. Leus, and M. Moonen, “Optimal training design for MIMO OFDM systems in mobile wireless channels,” _IEEE Transactions on Signal Processing_ , vol. 51, no. 6, pp. 1615–1624, 2003.
  * [13] M. Biguesh and A. Gershman, “Training-based MIMO channel estimation: a study of estimator tradeoffs and optimal training signals,” _IEEE Transactions on Signal Processing_ , vol. 54, no. 3, pp. 884–893, 2006.
  * [14] M. Haardt, F. Roemer, and G. Del Galdo, “Higher-order SVD-based subspace estimation to improve the parameter estimation accuracy in multidimensional harmonic retrieval problems,” _IEEE Transactions on Signal Processing_ , vol. 56, no. 7, pp. 3198–3213, 2008.
  * [15] P. Comon, X. Luciani, and A. L. F. de Almeida, “Tensor decompositions, alternating least squares and other tales,” _Journal of Chemometrics_ , vol. 23, no. 7-8, pp. 393–405, 2009.
  * [16] C. F. Van Loan and N. Pitsianis, _Approximation with Kronecker Products_ , M. S. Moonen, G. H. Golub, and B. L. R. De Moor, Eds. Dordrecht: Springer Netherlands, 1993.
  * [17] K. K. Wu, Y. Yam, H. Meng, and M. Mesbahi, “Kronecker product approximation with multiple factor matrices via the tensor product algorithm,” in _Proc. IEEE International Conference on Systems, Man, and Cybernetics (SMC)_ , 2016, pp. 004 277–004 282.
  * [18] C. Garvey, C. Meng, and J. G. Nagy, “Singular value decomposition approximation via kronecker summations for imaging applications,” 2018. [Online]. Available: https://arxiv.org/abs/1803.11525
  * [19] 3GPP, “Study on Channel Model for Frequencies from 0.5 to 100 GHz,” 3rd Generation Partnership Project (3GPP), Technical Report (TR) 38.901, 2020, v16.1.0.
  * [20] D. G. Riviello, F. Di Stasio, and R. Tuninato, “Performance analysis of multi-user MIMO schemes under realistic 3GPP 3-D channel model for 5G mmwave cellular networks,” _Electronics_ , vol. 11, no. 3, 2022. [Online]. Available: https://www.mdpi.com/2079-9292/11/3/330


(2\pi)^2\frac{\rmd I^{\rm{in-out}}}{\rmd z \rmd^2 \p} &= -\frac{2\alpha_s}{z(1-z)E}P_{ba}(z)
\rmR \,i \int_0^L\rmd \tau \,\rme^{-i \frac{\p^2}{2 \omega}\tau}\rme^{-\frac{1}{12} \frac{\qhat \p^2}{\omega^2}  \tau^3 }\,.
\end{align}

The four-point function is not trivial to calculate in the eikonal approximation either. However, it can now be calculated through an ordinary differential equation, instead of a more complicated Schrödinger equation

\begin{equation}
\frac{\rmd}{\rmd t}\Cc^{(4)}_i(t,t_2) = \M_{ij}(t)\Cc_j(t)\,.
\end{equation}

The matrix in Eq. (<ref>) now has a simple time dependence

\begin{equation}\label{eq:potential-matrix-pair-eikonal}
\mathbb{M}(t) =
-\frac{\qhat \theta^2}{4C_F}
\begin{bmatrix}
& -\frac{1}{N_c}(t-t_1)(t-t_2) \\
N_c z(1-z)(t_2-t_1)^2
& [C_F-N_c z(1-z)](t_2-t_1)^2
\end{bmatrix}\,.
\end{equation}

This can be solved numerically.

In the large-$N_c$ approximation the above matrix simplifies

\begin{equation}\label{eq:potential-matrix-pair-eikonal-Nc}
\mathbb{M}(t) =
-\frac{\qhat \theta^2}{4}
\begin{bmatrix}
& 0 \\
& [1-2 z(1-z)](t_2-t_1)^2
\end{bmatrix}\,.
\end{equation}

The zero in the upper right entry means that $\Cc^{(4)}_1$ can be solved analytically

\begin{equation}
\Cc^{(4)}_1(L,t_2) = \rme^{-\frac{1}{12}\frac{\qhat \p^2}{\omega^2} \left[(L-t_2)^3+(L-t_1)^3-(t_2-t_1)^3\right]}\,.
\end{equation}

The physical solution that we need in the emission spectrum is $\Cc^{(4)}_2$. This can be solved using the solution of
$\Cc^{(4)}_1$ through

\begin{equation}
\Cc^{(4)}_2(L,t_2) = \Cc^{(4)}_{2,\rm{fac}}(L,t_2) -\frac12 \frac{\qhat \p^2}{\omega^2}z(1-z) (t_2-t_1)^2
\int_{t_2}^L \rmd s \,\Cc^{(4)}_{2,\rm{fac}}(L,s) \Cc^{(4)}_1(s,t_2)\,,
\end{equation}

where the factorizable solution is

\begin{equation}
\Cc^{(4)}_{2,\rm{fac}}(L,t_2) = \rme^{-\frac{1}{4}\frac{\qhat \p^2}{\omega^2} [1-2z(1-z)](L-t_2)(t_2-t_1)^2}\,.
\end{equation}

This leads to the in-in spectrum

\begin{align}\label{eq:spectrum-hom-eikonal}
(2\pi)^2\frac{\rmd I^{\rm{in-in}}}{\rmd z \rmd^2 \p} &=
\frac{4\alpha_s }{\qhat[1-2z(1-z)]} P_{ba}(z)\nn
&\times\rmR \, \int_0^L\rmd \tau \,\frac{1}{\tau^2}\rme^{-i \frac{\p^2}{2 \omega}\tau}\rme^{-\frac{1}{12} \frac{\qhat \p^2}{\omega^2}  \tau^3 }  \left(1-\rme^{-\frac{1}{4}\frac{\qhat \p^2}{\omega^2} [1-2z(1-z)](L-\tau)\tau^2}\right)\,.
\end{align}