# HUMUS-Net: Hybrid Unrolled Multi-Scale Network Architecture for Accelerated
MRI Reconstruction

Zalan Fabian
Dept. of Electrical and Computer Engineering
University of Southern California
Los Angeles, CA
<EMAIL_ADDRESS>
&Berk Tinaz
Dept. of Electrical and Computer Engineering
University of Southern California
Los Angeles, CA
<EMAIL_ADDRESS>
&Mahdi Soltanolkotabi
Dept. of Electrical and Computer Engineering
University of Southern California
Los Angeles, CA
<EMAIL_ADDRESS>

###### Abstract

In accelerated MRI reconstruction, the anatomy of a patient is recovered from
a set of under-sampled and noisy measurements. Deep learning approaches have
been proven to be successful in solving this ill-posed inverse problem and are
capable of producing very high quality reconstructions. However, current
architectures heavily rely on convolutions, that are content-independent and
have difficulties modeling long-range dependencies in images. Recently,
Transformers, the workhorse of contemporary natural language processing, have
emerged as powerful building blocks for a multitude of vision tasks. These
models split input images into non-overlapping patches, embed the patches into
lower-dimensional tokens and utilize a self-attention mechanism that does not
suffer from the aforementioned weaknesses of convolutional architectures.
However, Transformers incur extremely high compute and memory cost when 1) the
input image resolution is high and 2) when the image needs to be split into a
large number of patches to preserve fine detail information, both of which are
typical in low-level vision problems such as MRI reconstruction, having a
compounding effect. To tackle these challenges, we propose HUMUS-Net, a hybrid
architecture that combines the beneficial implicit bias and efficiency of
convolutions with the power of Transformer blocks in an unrolled and multi-
scale network. HUMUS-Net extracts high-resolution features via convolutional
blocks and refines low-resolution features via a novel Transformer-based
multi-scale feature extractor. Features from both levels are then synthesized
into a high-resolution output reconstruction. Our network establishes new
state of the art on the largest publicly available MRI dataset, the fastMRI
dataset. We further demonstrate the performance of HUMUS-Net on two other
popular MRI datasets and perform fine-grained ablation studies to validate our
design.

††36th Conference on Neural Information Processing Systems 2022

##  1 Introduction

Magnetic resonance imaging (MRI) is a medical imaging technique that uses
strong magnetic fields to picture the anatomy and physiological processes of
the patient. MRI is one of the most popular imaging modalities as it is non-
invasive and doesn’t expose the patient to harmful ionizing radiation. The MRI
scanner obtains measurements of the body in the spatial frequency domain also
called k-space. The data acquisition process is MR is often time-consuming.
Accelerated MRI (Lustig et al., 2008) addresses this challenge by
undersampling in the k-space domain, thus reducing the time patients need to
spend in the scanner. However, recovering the underlying anatomy from
undersampled measurements is an ill-posed problem (less measurements than
unknowns) and thus incorporating some form of prior knowledge is crucial in
obtaining high quality reconstructions. Classical MRI reconstruction
algorithms rely on the assumption that the underlying signal is sparse in some
transform domain and attempt to recover a signal that best satisfies this
assumption in a technique known as compressed sensing (CS) (Candes et al.,
2006; Donoho, 2006). These classical CS techniques have slow reconstruction
speed and typically enforce limited forms of image priors.

With the emergence of deep learning (DL), data-driven reconstruction
algorithms have far surpassed CS techniques (see Ongie et al. (2020) for an
overview). DL models utilize large training datasets to extract flexible,
nuanced priors directly from data resulting in excellent reconstruction
quality. In recent years, there has been a flurry of activity aimed at
designing DL architectures tailored to the MRI reconstruction problem. The
most popular models are convolutional neural networks (CNNs) that typically
incorporate the physics of the MRI reconstruction problem, and utilize tools
from mainstream deep learning (residual learning, data augmentation, self-
supervised learning). Comparing the performance of such models has been
difficult mainly due to two reasons. First, there has been a large variation
in evaluation datasets spanning different scanners, anatomies, acquisition
models and undersampling patterns rendering direct comparison challenging.
Second, medical imaging datasets are often proprietary due to privacy
concerns, hindering reproducibility.

More recently, the fastMRI dataset (Zbontar et al., 2019), the largest
publicly available MRI dataset, has been gaining ground as a standard
benchmark to evaluate MRI reconstruction methods. An annual competition, the
fastMRI Challenge (Muckley et al., 2021), attracts significant attention from
the machine learning community and acts a driver of innovation in MRI
reconstruction. However, over the past years the public leaderboard has been
dominated by a single architecture, the End-to-End VarNet (Sriram et al.,
2020) 111At the time of writing this paper, the number one model on the
leaderboard is AIRS-Net, however the complete architecture/training details as
well as the code of this method are not available to the public and we were
unable to reproduce their results based on the limited available information.,
with most models concentrating very closely around the same performance
metrics, hinting at the saturation of current architectural choices.

In this work, we propose HUMUS-Net: a Hybrid, Unrolled, MUlti-Scale network
architecture for accelerated MRI reconstruction that combines the advantages
of well-established architectures in the field with the power of contemporary
Transformer-based models. We utilize the strong implicit bias of convolutions,
but also address their weaknesses, such as content-independence and inability
to model long-range dependencies, by incorporating a novel multi-scale feature
extractor that operates over embedded image patches via self-attention.
Moreover, we tackle the challenge of high input resolution typical in MRI by
performing the computationally most expensive operations on extracted low-
resolution features. HUMUS-Net establishes new state of the art in accelerated
MRI reconstruction on the largest available MRI knee dataset. At the time of
writing this paper, HUMUS-Net is the only Transformer-based architecture on
the highly competitive fastMRI Public Leaderboard. Our results are fully
reproducible and the source code is available online 222Code:
https://github.com/MathFLDS/HUMUS-Net.

##  2 Background

###  2.1 Inverse Problem Formulation of Accelerated MRI Reconstruction

An MR scanner obtains measurements of the patient anatomy in the frequency
domain, also referred to as k-space. Data acquisition is performed via various
receiver coils positioned around the anatomy being imaged, each with different
spatial sensitivity. Given a total number of $N$ receiver coils, the
measurements obtained by the $i$th coil can be written as

$\bm{k_{i}}={\mathcal{F}}\bm{S_{i}}\bm{x^{*}}+\bm{z_{i}},~{}~{}i=1,..,N,$

where $\bm{x^{*}}\in\mathbb{C}^{n}$ is the underlying patient anatomy of
interest, $\bm{S_{i}}$ is a diagonal matrix that represents the sensitivity
map of the $i$th coil, ${\mathcal{F}}$ is a multi-dimensional Fourier-
transform, and $z_{i}$ denotes the measurement noise corrupting the
observations obtained from coil $i$. We use
$\bm{k}=(\bm{k}_{1},...,\bm{k}_{N})$ as a shorthand for the concatenation of
individual coil measurements and $\bm{x}=(\bm{x}_{1},...,\bm{x}_{N})$ as the
corresponding image domain representation after inverse Fourier
transformation.

Since MR data acquisition time is proportional to the portion of k-space being
scanned, obtaining fully-sampled data is time-consuming. Therefore, in
accelerated MRI scan times are reduced by undersampling in k-space domain. The
undersampled k-space measurements from coil $i$ take the form

$\bm{\tilde{k}_{i}}={\bm{M}}\bm{k_{i}}~{}~{}i=1,..,N,$

where ${\bm{M}}$ is a diagonal matrix representing the binary undersampling
mask, that has $0$ values for all missing frequency components that have not
been sampled during accelerated acquisition.

The forward model that maps the underlying anatomy to coil measurements can be
written concisely as $\bm{\tilde{k}}={\mathcal{A}\left(\bm{x^{*}}\right)},$
where ${\mathcal{A}\left(\cdot\right)}$ is the linear forward mapping and
$\bm{\tilde{k}}$ is the stacked vector of all undersampled coil measurements.
Our target is to reconstruct the ground truth object $\bm{x^{*}}$ from the
noisy, undersampled measurements $\bm{\tilde{k}}$. Since we have fewer
observations than variables to recover, perfect reconstruction in general is
not possible. In order to make the problem solvable, prior knowledge on the
underlying object is typically incorporated in the form of sparsity in some
transform domain. This formulation, known as compressed sensing (Candes et
al., 2006; Donoho, 2006), provides a classical framework for accelerated MRI
reconstruction (Lustig et al., 2008). In particular, the above recovery
problem can be formulated as a regularized inverse problem

$\bm{\hat{x}}=\operatorname*{arg\,min}_{\bm{x}}\left\|{\mathcal{A}\left(\bm{x}\right)}-\bm{\tilde{k}}\right\|^{2}+\mathcal{R}(\bm{x}),$
(2.1)

where $\mathcal{R}(\cdot)$ is a regularizer that encapsulates prior knowledge
on the object, such as sparsity in some wavelet domain.

###  2.2 Deep Learning-based Accelerated MRI Reconstruction

More recently, data-driven deep learning-based algorithms tailored to the
accelerated MRI reconstruction problem have surpassed the classical compressed
sensing baselines. Convolutional neural networks trained on large datasets
have established new state of the art in many medical imaging tasks. The
highly popular U-Net (Ronneberger et al., 2015) and other similar encoder-
decoder architectures have proven to be successful in a range of medical image
reconstruction (Hyun et al., 2018; Han and Ye, 2018) and segmentation (Çiçek
et al., 2016; Zhou et al., 2018) problems. In the encoder path, the network
learns to extract a set of deep, low-dimensional features from images via a
series of convolutional and downsampling operations. These concise feature
representations are then gradually upsampled and filtered in the decoder to
the original image dimensions. Thus the network learns a hierarchical
representation over the input image distribution.

Unrolled networks constitute another line of work that has been inspired by
popular optimization algorithms used to solve compressed sensing
reconstruction problems. These deep learning models consist of a series of
sub-networks, also known as cascades, where each sub-network corresponds to an
unrolled iteration of popular algorithms such as gradient descent (Zhang and
Ghanem, 2018) or ADMM (Sun et al., 2016). In the context of MRI
reconstruction, one can view network unrolling as solving a sequence of
smaller denoising problems, instead of the complete recovery problem in one
step. Various convolutional neural networks have been employed in the
unrolling framework achieving excellent performance in accelerated MRI
reconstruction (Putzky et al., 2019; Hammernik et al., 2018, 2019). E2E-VarNet
(Sriram et al., 2020) is the current state-of-the-art convolutional model on
the fastMRI dataset. E2E-VarNet transforms the optimization problem in (2.1)
to the k-space domain and unrolls the gradient descent iterations into $T$
cascades, where the $t$th cascade represents the computation

$\bm{\hat{k}^{t+1}}=\bm{\hat{k}^{t}}-\mu^{t}{\bm{M}}(\bm{\hat{k}^{t}}-\bm{\tilde{k}})+\mathcal{G}(\bm{\hat{k}^{t}}),$
(2.2)

where $\bm{\hat{k}^{t}}$ is the estimated reconstruction in the k-space domain
at cascade $t$, $\mu^{t}$ is a learnable step size parameter and
$\mathcal{G}(\cdot)$ is a learned mapping representing the gradient of the
regularization term in (2.1). The first term is also known as data consistency
(DC) term as it enforces the consistency of the estimate with the available
measurements.

Figure 1: Overview of the HUMUS-Block architecture. First, we extract high-
resolution features $\bm{F_{H}}$ from the input noisy image through a
convolution layer $f_{H}$. Then, we apply a convolutional feature extractor
$f_{L}$ to obtain low-resolution features and process them using a
Transformer-convolutional hybrid multi-scale feature extractor. The shallow,
high-resolution and deep, low-resolution features are then synthesized into
the final high-resolution denoised image.

##  3 Related Work

Transformers in Vision – Vision Transformer (ViT) (Dosovitskiy et al., 2020),
a fully non-convolutional vision architecture, has demonstrated state-of-the-
art performance on image classification problems when pre-trained on large-
scale image datasets. The key idea of ViT is to split the input image into
non-overlapping patches, embed each patch via a learned linear mapping and
process the resulting tokens via stacked self-attention and multi-layer
perceptron (MLP) blocks. For more details we refer the reader to Appendix G
and (Dosovitskiy et al., 2020). The benefit of Transformers over convolutional
architectures in vision lies in their ability to capture long-range
dependencies in images via the self-attention mechanism.

Since the introduction of ViT, similar attention-based architectures have been
proposed for many other vision tasks such as object detection (Carion et al.,
2020), image segmentation (Wang et al., 2021b) and restoration (Cao et al.,
2021b; Chen et al., 2021; Liang et al., 2021; Zamir et al., 2021; Wang et al.,
2021c). A key challenge for Transformers in low-level vision problems is the
quadratic compute complexity of the global self-attention with respect to the
input dimension. In some works, this issue has been mitigated by splitting the
input image into fixed size patches and processing the patches independently
(Chen et al., 2021). Others focus on designing hierarchical Transformer
architectures (Heo et al., 2021; Wang et al., 2021b) similar to popular
ResNets (He et al., 2015). Authors in Zamir et al. (2021) propose applying
self-attention channel-wise rather than across the spatial dimension thus
reducing the compute overhead to linear complexity. Another successful
architecture, the Swin Transformer (Liu et al., 2021), tackles the quadratic
scaling issue by computing self-attention on smaller local windows. To
encourage cross-window interaction, windows in subsequent layers are shifted
relative to each other.

Transformers in Medical Imaging – Transformer architectures have been proposed
recently to tackle various medical imaging problems. Authors in Cao et al.
(2021a) design a U-Net-like architecture for medical image segmentation where
the traditional convolutional layers are replaced by Swin Transformer blocks.
They report strong results on multi-organ and cardiac image segmentation. In
Zhou et al. (2021) a hybrid convolutional and Transformer-based U-Net
architecture is proposed tailored to volumetric medical image segmentation
with excellent results on benchmark datasets. Similar encoder-decoder
architectures for various medical segmentation tasks have been investigated in
other works (Huang et al., 2021; Wu et al., 2022). However, these networks are
tailored for image segmentation, a task less sensitive to fine details in the
input, and thus larger patch sizes are often used (for instance $4$ in Cao et
al. (2021a) ). This allows the network to process larger input images, as the
number of token embedding is greatly reduced, but as we demonstrate in Section
5.2 embedding individual pixels as $1\times 1$ patches is crucial for MRI
reconstruction. Thus, compute and memory barriers stemming from large input
resolutions are more severe in the MRI reconstruction task and therefore novel
approaches are needed.

Promising results have been reported employing Transformers in medical image
denoising problems, such as low-dose CT denoising (Wang et al., 2021a; Luthra
et al., 2021) and low-count PET/MRI denoising (Zhang et al., 2021). However,
these studies fail to address the challenge of poor scaling to large input
resolutions, and only work on small images via either downsampling the
original dataset (Luthra et al., 2021) or by slicing the large input images
into smaller patches (Wang et al., 2021a). In contrast, our proposed
architecture works directly on the large resolution images that often arise in
MRI reconstruction. Even though some work exists on Transformer-based
architectures for supervised accelerated MRI reconstruction (Huang et al.,
2022; Lin and Heckel, 2022; Feng et al., 2021), and for unsupervised pre-
trained reconstruction (Korkmaz et al., 2022), to the best of our knowledge
ours is the first work to demonstrate state-of-the-art results on large-scale
MRI datasets such as the fastMRI dataset.

##  4 Method

HUMUS-Net combines the efficiency and beneficial implicit bias of
convolutional networks with the powerful general representations of
Transformers and their capability to capture long-range pixel dependencies.
The resulting hybrid network processes information both in image
representation (via convolutions) and in patch-embedded token representation
(via Transformer blocks). Our proposed architecture consists of a sequence of
sub-networks, also called cascades. Each cascade represents an unrolled
iteration of an underlying optimization algorithm in k-space (see (2.2)), with
an image-domain denoiser, the HUMUS-Block. First, we describe the architecture
of the HUMUS-Block, the core component in the sub-network. Then, we specify
the high-level, k-space unrolling architecture of HUMUS-Net in Section 4.3.

###  4.1 HUMUS-Block Architecture

The HUMUS-Block acts as an image-space denoiser that receives an intermediate
reconstruction from the previous cascade and performs a single step of
denoising to produce an improved reconstruction for the next cascade. It
extracts high-resolution, shallow features and low-resolution, deep features
through a novel multi-scale transformer-based block, and synthesizes high-
resolution features from those. The high-level overview of the HUMUS-Block is
depicted in Fig. 1.

High-resolution Feature Extraction– The input to our network is a noisy
complex-valued image $\bm{x}_{in}\in\mathbb{R}^{h\times w\times c_{in}}$
derived from undersampled k-space data, where the real and imaginary parts of
the image are concatenated along the channel dimension. First, we extract
high-resolution features $\bm{F_{H}}\in\mathbb{R}^{h\times w\times d_{H}}$
from the input noisy image through a convolution layer $f_{H}$ written as
$\bm{F_{H}}=f_{H}(\bm{x}_{in})$. This initial $3\times 3$ convolution layer
provides early visual processing at a relatively low cost and maps the input
to a higher, $d_{H}$ dimensional feature space. Important to note that the
resolution of the extracted features is the same as the spatial resolution of
the input image.

Low-resolution Feature Extraction– In case of MRI reconstruction, input
resolution is typically significantly higher than commonly used image datasets
($32\times 32$ \- $256\times 256$), posing a significant challenge to
contemporary Transformer-based models. Thus we apply a convolutional feature
extractor $f_{L}$ to obtain low-resolution features
$\bm{F_{L}}=f_{L}(\bm{F_{H}})$, with $\bm{F_{L}}\in\mathbb{R}^{h_{L}\times
w_{L}\times d_{L}}$ where $f_{L}$ consists of a sequence of convolutional
blocks and spatial downsampling operations. The specific architecture is
depicted in Figure 4. The purpose of this module is to perform deeper visual
processing and to provide manageable input size to the subsequent computation
and memory heavy hybrid processing module. In this work, we choose
$h_{L}=\frac{h}{2}$ and $w_{L}=\frac{w}{2}$ which strikes a balance between
preserving spatial information and resource demands. Furthermore, in order to
compensate for the reduced resolution we increase the feature dimension to
$d_{L}=2\cdot d_{H}:=d$.

Deep Feature Extraction– The most important part of our model is MUST, a
MUlti-scale residual Swin Transformer network. MUST is a multi-scale hybrid
feature extractor that takes the low-resolution image representations
$\bm{F_{L}}$ and performs hierarchical Transformer-convolutional hybrid
processing in an encoder-decoder fashion, producing deep features
$\bm{F_{D}}=f_{D}(\bm{F_{L}}),$ where the specific architecture behind $f_{D}$
is detailed in Subsection 4.2.

High-resolution Image Reconstruction– Finally, we combine information from
shallow, high-resolution features $\bm{F_{H}}$ and deep, low-resolution
features $\bm{F_{D}}$ to reconstruct the high-resolution residual image via a
convolutional reconstruction module $f_{R}$. The residual learning paradigm
allows us to learn the difference between noisy and clean images and helps
information flow within the network (He et al., 2015). Thus the final denoised
image $\bm{x}_{out}\in\mathbb{R}^{h\times w\times c_{in}}$ is obtained as
$\bm{x}_{out}=\bm{x}_{in}+f_{R}(\bm{F_{H}},\bm{F_{D}})$. The specific
architecture of the reconstruction network is depicted in Figure 4.

###  4.2 Multi-scale Hybrid Feature Extraction via MUST

The key component to our architecture is MUST, a multi-scale hybrid encoder-
decoder architecture that performs deep feature extraction in both image and
token representation (Figure 1, bottom). First, individual pixels of the input
representation of shape $\frac{h}{2}\times\frac{w}{2}\times{d}$ are flattened
and passed through a learned linear mapping to yield
$\frac{h}{2}\cdot\frac{w}{2}$ tokens of dimension $d$. Tokens corresponding to
different image patches are subsequently merged in the encoder path, resulting
in a concise latent representation. This highly descriptive representation is
passed through a bottleneck block and progressively expanded by combining
tokens from the encoder path via skip connections. The final output is
rearranged to match the exact shape of the input low-resolution features
$\bm{F_{L}}$, yielding a deep feature representation $\bm{F_{D}}$.

Our design is inspired by the success of Residual Swin Transformer Blocks
(RSTB) in image denoising and super-resolution (Liang et al., 2021). RSTB
features a stack of Swin Transformer layers (STL) that operate on tokens via a
windowed self-attention mechanism (Liu et al., 2021), followed by convolution
in image representation. However, RSTB blocks operate on a single scale,
therefore they cannot be readily applied in a hierarchical encoder-decoder
architecture. Therefore, we design three variations of RSTB to facilitate
multi-scale processing as depicted in Figure 2.

RSTB-B is the bottleneck block responsible for processing the encoded latent
representation while maintaining feature dimensions. Thus, we keep the default
RSTB architecture for our bottleneck block, which already operates on a single
scale.

RSTB-D has a similar function to convolutional downsampling blocks in U-Nets,
but it operates on embedded tokens. Given an input with size $h_{i}\cdot
w_{i}\times d$, we pass it through an RSTB-B block and apply PatchMerge
operation. PatchMerge linearly combines tokens corresponding to $2\times 2$
non-overlapping image patches, while simultaneously increasing the embedding
dimension (see Figure 3, top and Figure 10 in the appendix for more details)
resulting in an output of size $\frac{h_{i}}{2}\cdot\frac{w_{i}}{2}\times
2\cdot d$. Furthermore, RSTB-D outputs the higher dimensional representation
before patch merging to be subsequently used in the decoder path via skip
connection.

RSTB-U used in the decoder path is analogous to convolutional upsampling
blocks. An input with size $h_{i}\cdot w_{i}\times d$ is first expanded into a
larger number of lower dimensional tokens through a linear mapping via
PatchExpand (see Figure 3, bottom and Figure 10 in the appendix for more
details). PatchExpand reverses the effect of PatchMerge on feature size, thus
resulting in $2h_{i}\cdot 2w_{i}$ tokens of dimension $\frac{d}{2}$. Next, we
mix information from the obtained expanded tokens with skip embeddings from
higher scales via TokenMix. This operation linearly combines tokens from both
paths and normalizes the resulting vectors. Finally, the mixed tokens are
processed by an RSTB-B block.

Figure 2: Depiction of different RSTB modules used in the HUMUS-Block.
Figure 3: Patch merge and expand operations used in our multi-scale feature
extractor.

Figure 4: Architecture of convolutional blocks for feature extraction and
reconstruction.

###  4.3 Iterative Unrolling

Architectures derived from unrolling the iterations of various optimization
algorithms have proven to be successful in tackling various inverse problems
including MRI reconstruction. These architecture can be interpreted as a
cascade of simpler denoisers, each of which progressively refines the estimate
from the preceding unrolled iteration (see more details in Appendix F).

Following Sriram et al. (2020), we unroll the gradient descent iterations of
the inverse problem in (2.1) in k-space domain, yielding the iterative update
scheme in (2.2). We apply regularization in image domain via our proposed
HUMUS-Block, that is we have
$\mathcal{G}(\bm{k})={\mathcal{F}}\left(\mathcal{E}\left(\textbf{D}\left(\mathcal{R}\left({\mathcal{F}^{-1}}(\bm{k})\right)\right)\right)\right),$
where D denotes the HUMUS-Block,
$\mathcal{R}(\bm{x_{1}},...,\bm{x_{N}})=\sum_{i=1}^{N}\bm{S}_{i}^{*}x_{i}$ is
the reduce operator that combines coil images via the corresponding
sensitivity maps and
$\mathcal{E}(\bm{x})=\left(\bm{S}_{1}\bm{x},...,\bm{S}_{N}\bm{x}\right)$ is
the expand operator that maps the combined image back to individual coil
images. The sensitivity maps can be estimated a priori using methods such as
ESPIRiT (Uecker et al., 2014) or learned in an end-to-end fashion via a
Sensitivity Map Estimator (SME) network proposed in Sriram et al. (2020). In
this work we aspire to design an end-to-end approach and thus we use the
latter method and estimate the sensitivity maps from the low-frequency (ACS)
region of the undersampled input measurements during training using a standard
U-Net network. A visual overview of our unrolled architecture is shown in
Figure 9 in the appendix.

###  4.4 Adjacent Slice Reconstruction (ASR)

We observe improvements in reconstruction quality when instead of processing
the undersampled data slice-by-slice, we jointly reconstruct a set of adjacent
slices via HUMUS-Net (Figure 5). That is, if we have a volume of undersampled
data $\tilde{\bm{k}}^{vol}=(\tilde{\bm{k}}^{1},...,\tilde{\bm{k}}^{K})$ with
$K$ slices, when reconstructing slice $c$, we instead reconstruct the volume
$(\tilde{\bm{k}}^{c-a},...,\tilde{\bm{k}}^{c-1},\tilde{\bm{k}}^{c},\tilde{\bm{k}}^{c+1},...,\tilde{\bm{k}}^{c+a})$
by concatenating them along the coil channel dimension, where $a$ denotes the
number of adjacent slices added on each side. However, we only calculate and
backpropagate the loss on the center slice $c$ of the reconstructed volume.
The benefit of ASR is that the network can remove artifacts corrupting
individual slices as it sees a larger context of the slice by observing its
neighbors. Even though ASR increases compute cost, it is important to note
that it does not impact the number of token embeddings (spatial resolution is
unchanged) and thus can be combined favorably with Transformer-based methods.

Figure 5:  Adjacent slice reconstruction (depicted in image domain for visual
clarity): HUMUS-Net takes a volume of adjacent slices
$(\tilde{\bm{x}}^{c-a},...,\tilde{\bm{x}}^{c},...,\tilde{\bm{x}}^{c+a})$ and
jointly reconstructs a volume
$(\hat{\bm{x}}^{c-a},...,\hat{\bm{x}}^{c},...,\hat{\bm{x}}^{c+a})$. The
reconstruction loss $\mathcal{L}$ is calculated only on the center slice
$\hat{\bm{x}}^{c}$.

##  5 Experiments

In this section we provide experimental results on our proposed architecture,
HUMUS-Net. First, we demonstrate the reconstruction performance of our model
on various datasets, including the large-scale fastMRI dataset. Then, we
justify our design choices through a set of ablation studies.

###  5.1 Benchmark Experiments

We investigate the performance of HUMUS-Net on three different datasets. We
use the structural similarity index measure (SSIM) (Wang et al., 2003) as a
basis of our evaluation, which is the most common evaluation metric in medical
image reconstruction. In all of our experiments, we follow the setup of the
fastMRI multi-coil knee track with $8\times$ acceleration in order to provide
comparison with state-of-the-art networks of the Public Leaderboard. That is,
we perform retrospective undersampling of the fully-sampled k-space data by
randomly subsampling $12.5\%$ of whole k-space lines in the phase encoding
direction, keeping $4\%$ of lowest frequency adjacent k-space lines.
Experiments on other acceleration ratios can be found in Appendix C. During
training, we generate random masks following the above method, whereas for the
validation dataset we keep the masks fixed for each k-space volume. For HUMUS-
Net, we center crop and pad inputs to $384\times 384$. We compare the
reconstruction quality of our proposed model with the current best performing
network, E2E-VarNet. For details on HUMUS-Net hyperparameters and training, we
refer the reader to Appendix A. For E2E-VarNet, we use the hyperparameters
specified in Sriram et al. (2020).

fastMRI knee multi-coil $\times 8$ test data
---
Method |  # of params (approx.) |  SSIM($\uparrow$) |  PSNR($\uparrow$) |  NMSE($\downarrow$)
E2E-VarNet (Sriram et al., 2020) | 30M | $0.8900$ | $36.9$ | $0.0089$
E2E-VarNet†(Sriram et al., 2020) | $0.8920$ | $37.1$ | $0.0085$
XPDNet (Ramzi et al., 2020) | 155M | $0.8893$ | $37.2$ | $0.0083$
$\Sigma$-Net (Hammernik et al., 2019) | 676M | $0.8877$ | $36.7$ | $0.0091$
i-RIM (Putzky et al., 2019) | 300M | $0.8875$ | $36.7$ | $0.0091$
U-Net (Zbontar et al., 2019) | 214M | $0.8640$ | $34.7$ | $0.0132$
HUMUS-Net (ours) | 109M | $0.8936$ | $37.0$ | $0.0086$
HUMUS-Net (ours)† | $0.8945$ | $37.3$ | $0.0081$
HUMUS-Net-L (ours) | 228M | $0.8944$ | $37.3$ | $0.0083$
HUMUS-Net-L (ours)† | 0.8951 | 37.4 | 0.0080

Table 1:  Performance of state-of-the-art accelerated MRI reconstruction
techniques on the fastMRI knee test dataset. Most models are trained only on
the fastMRI training dataset, if available we show results of models trained
on the fastMRI combined training and validation dataset denoted by
(${\dagger}$).

Validation SSIM
---
Method | fastMRI knee | Stanford 2D | Stanford 3D
E2E-VarNet (Sriram et al., 2020) | $0.8908$ | $0.8928\pm 0.0168$ | $0.9432\pm 0.0063$
HUMUS-Net (ours) | 0.8934 | $\textbf{0.8954}\pm\textbf{0.0136}$ | $\textbf{0.9453}\pm\textbf{0.0065}$

Table 2: Validation SSIM of HUMUS-Net on various datasets. For datasets with
multiple train-validation split runs we show the mean and standard error of
the runs.

Ablation studies
---
Method | Unrolled? | Multi-scale? | Low-res features? | Patch size | Embed. dim. | SSIM
Un-SS | $\checkmark$ | ✗ | ✗ | $1$ | $12$ | $0.9319\pm 0.0080$
Un-MS | $\checkmark$ | $\checkmark$ | ✗ | $1$ | $12$ | $0.9357\pm 0.0038$
Un-MS-Patch2 | $\checkmark$ | $\checkmark$ | ✗ | $2$ | $36$ | $0.9171\pm 0.0075$
HUMUS-Net | $\checkmark$ | $\checkmark$ | $\checkmark$ | $1$ | $36$ | 0.9449 $\pm$ 0.0064
SwinIR | $0.9336\pm 0.0069$
E2E-VarNet | $0.9432\pm 0.0063$

Table 3: Results of ablation studies on HUMUS-Net, evaluated on the Stanford
3D MRI dataset.

fastMRI – The fastMRI dataset (Zbontar et al., 2019) is the largest publicly
available MRI dataset with competitive baseline models and a public
leaderboard, and thus provides an opportunity to directly compare different
algorithms. Specifically, we run experiments on the multi-coil knee dataset,
consisting of close to $35k$ slices in $973$ volumes. We use the default
HUMUS-Net model defined above with $3$ adjacent slices as input. Furthermore,
we design a large variant of our model, HUMUS-Net-L, which has increased
embedding dimension compared to the default model (see details in Appendix A).
We train models both only on the training split, and also on the training and
validation splits combined (additional $\approx 20\%$ data) for the
leaderboard. Table 3 demonstrates our results compared to the best published
models from the fastMRI Leaderboard333https://fastmri.org/leaderboards
evaluated on the test dataset. Our model establishes new state of the art in
terms of SSIM on this dataset by a large margin, and achieves comparable or
better performance than other methods in terms of PSNR and NMSE. Moreover, as
seen in the second column of Table 3, we evaluated our model on the fastMRI
validation dataset as well and compared our results to E2E-VarNet, the best
performing model from the leaderboard. We observe similar improvements in
terms of the reconstruction SSIM metric to the test dataset. Visual inspection
of reconstructions shows that HUMUS-Net recovers very fine details in images
that may be missed by other state-of-the-art reconstruction algorithms (see
attached figures in Figure 11 in the appendix). Comparison based on further
image quality metrics can be found in Appendix H. We point out that even
though our model has more parameters than E2E-VarNet, our proposed image-
domain denoiser is more efficient than the U-Net deployed in E2E-VarNet, even
when their number of parameters are matched as discussed in Section 5.3.
Overall, larger model size does not necessarily correlate with better
performance, as seen in other competitive models in Table 3. Finally, we note
that the additional training data (in the form of the validation split)
provides a consistent small boost to model performance. We refer the reader to
Klug and Heckel (2022) for an overview of scaling properties of reconstruction
models.

Stanford 2D – Next, we run experiments on the Stanford2D FSE (Cheng, )
dataset, a publicly available MRI dataset consisting of scans from various
anatomies (pelvis, lower extremity and more) in $89$ fully-sampled volumes. We
randomly sample $80\%$ of volumes as train data and use the rest for
validation. We randomly generate $3$ different train-validation splits this
way to reduce variations in the presented metrics. As slices in this dataset
have widely varying shapes across volumes, we center crop the target images to
keep spatial resolution within $384\times 384$. We use the default HUMUS-Net
defined above with single slices as input. Our results comparing the best
performing MRI reconstruction model with HUMUS-Net is shown in the third
column of Table 3. We present the mean SSIM of all runs along with the
standard error. We achieve improvements of similar magnitude as on the fastMRI
dataset. These results demonstrate the effectiveness of HUMUS-Net on a more
diverse dataset featuring multiple anatomies.

Stanford 3D – Finally, we evaluate our model on the Stanford Fullysampled 3D
FSE Knees dataset (Sawyer et al., 2013), a public MRI dataset including $20$
volumes of knee MRI scans. We generate train-validation splits using the
method described for Stanford 2D and perform $3$ runs. We use the default
HUMUS-Net network with single slices as input. The last column of Table 3
compares our results to E2E-VarNet, showing improvements of similar scales as
on other datasets we have investigated in this work. This experiment
demonstrates that HUMUS-Net performs well not only on large-scale MRI
datasets, but also on smaller problems.

###  5.2 Ablation Studies

In this section, we motivate our design choices through a set of ablation
studies. We start from SwinIR, a general image reconstruction network, and
highlight its weaknesses for MRI. Then, we demonstrate step-by-step how we
addressed these shortcomings and arrived at the HUMUS-Net architecture. We
train the models on the Stanford 3D dataset. More details can be found in
Appendix B. The results of our ablation studies are summarized in Table 3.

First, we investigate SwinIR, a state-of-the-art image denoiser and super-
resolution model that features a hybrid Transformer-convolutional
architecture. In order to handle the $10\times$ larger input sizes ($320\times
320$) in our MRI dataset compared to the input images this network has been
designed for ($128\times 128$), we reduce the embedding dimension of SwinIR to
fit into GPU memory (16 GB). We find that compared to models designed for MRI
reconstruction, such as E2E-VarNet (last row in Table 3) SwinIR performs
poorly. This is not only due to the reduced network size, but also due to the
fact that SwinIR is not specifically designed to take the MRI forward model
into consideration.

Next, we unroll SwinIR and add a sensitivity map estimator. We refer to this
model as Un-SS. Due to unrolling, we have to further reduce the embedding
dimension of the denoiser and also decrease the depth of the network in order
to fit into GPU memory. Un-SS, due to its small size, performs slightly worse
than vanilla SwinIR and significantly lags behind the E2E-VarNet
architectures. We note that SwinIR operates over a single, full-resolution
scale, whereas state-of-the-art MRI reconstruction models typically
incorporate multi-scale processing in the form of U-Net-like architectures.

Thus, we replace SwinIR by MUST, our proposed multi-scale hybrid processing
unit, but keep the embedding dimension in the largest-resolution scale fixed.
The obtained network, which we call Un-MS, has overall lower computational
cost when compared with Un-SS, however as Table 3 shows MRI reconstruction
performance has significantly improved compared to both Un-SS and vanilla
SwinIR, which highlights the efficiency of our proposed multi-scale feature
extractor. Reconstruction performance is limited by the low dimension of patch
embeddings, which we are unable to increase further due to our compute and
memory constraints originating in the high-resolution inputs.

The most straightforward approach to tackle the challenge of high input
resolution is to increase the patch size. To test this idea, we take Un-MS and
embed $2\times 2$ patches of the inputs, thus reducing the number of tokens
processed by the network by a factor of $4$. We refer to this model as Un-MS-
Patch2. This reduction in compute and memory load allows us to increase
network capacity by increasing the embedding dimension $3$-folds (to fill GPU
memory again). However, Un-MS-Patch2 performs much worse than previous models
using patch size of $1$. For classification problems, where the general image
context is more important than small details, patches of $16\times 16$ or
$8\times 8$ are considered typical (Dosovitskiy et al., 2020). Even for more
dense prediction tasks such as medical image segmentation, patch size of
$4\times 4$ has been used successfully (Cao et al., 2021a). However, our
experiments suggest that in low-level vision tasks such as MRI reconstruction
using patches larger than $1\times 1$ may be detrimental due to loss of
crucial high-frequency detail information.

Our approach to address the heavy computational load of Transformers for large
input resolutions where increasing the patch size is not an option is to
process lower resolution features extracted via convolutions. That is we
replace MUST in Un-MS by a HUMUS-Block, resulting in our proposed HUMUS-Net
architecture. We train a smaller version of HUMUS-Net with the same embedding
dimension as Un-MS. As seen in Table 3, our model achieves the best
performance across all other proposed solutions, even surpassing E2E-VarNet.
This series of incremental studies highlights the importance of each
architectural design choice leading to our proposed HUMUS-Net architecture.
Further ablation studies on the effect of the number of unrolled iterations
and adjacent slice reconstruction can be found in Appendix D and Appendix E
respectively.

###  5.3 Direct comparison of image-domain denoisers

In order to further demonstrate the advantage of HUMUS-Net over E2E-VarNet, we
provide direct comparison between the image-domain denoisers used in the above
methods. E2E-VarNet unrolls a fully convolutional U-Net, whereas we deploy our
hybrid HUMUS-Block architecture as a denoiser (Fig. 1). We scale down the
HUMUS-Block used in HUMUS-Net to match the size of the U-Net in E2E-Varnet in
terms of number of model parameters. To this end, we reduce the embedding
dimension from $66$ to $30$. We train both networks on magnitude images from
the Stanford3D dataset. The results are summarized in Table 4. We observe that
given a fixed parameter budget, our proposed denoiser outperforms the widely
used convolutional U-Net architecture in MRI reconstruction, further
demonstarting the efficiency of HUMUS-Net. This experiment suggests that our
HUMUS-Block could serve as an excellent denoiser, replacing convolutional
U-Nets, in a broad range of image restoration applications outside of MRI,
which we leave for future work.

Model | # of parameters | SSIM($\uparrow$) | PSNR($\uparrow$) | NMSE($\downarrow$)
---|---|---|---|---
U-Net | 2.5M | $0.9348\pm 0.0072$ | $39.0\pm 0.6$ | $0.0257\pm 0.0007$
HUMUS-Block | 2.4M | $\textbf{0.9378}\pm\textbf{0.0065}$ | $\textbf{39.2}\pm\textbf{0.5}$ | $\textbf{0.0246}\pm\textbf{0.0004}$
Table 4:  Direct comparison of denoisers on the Stanford 3D dataset. Mean and
standard error of $3$ random training-validation splits is shown.

##  6 Conclusion

In this paper, we introduce HUMUS-Net, an unrolled, Transformer-convolutional
hybrid network for accelerated MRI reconstruction. HUMUS-Net achieves state-
of-the-art performance on the fastMRI dataset and greatly outperforms all
previous published and reproducible methods. We demonstrate the performance of
our proposed method on two other MRI datasets and perform fine-grained
ablation studies to motivate our design choices and emphasize the compute and
memory challenges of Transformer-based architectures on low-level and dense
computer vision tasks such as MRI reconstruction. A limitation of our current
architecture is that it requires fixed-size inputs, which we intend to address
with a more flexible design in the future. This work opens the door for the
adoption of a multitude of promising techniques introduced recently in the
literature for Transformers, which we leave for future work.

## Acknowledgments

M. Soltanolkotabi is supported by the Packard Fellowship in Science and
Engineering, a Sloan Research Fellowship in Mathematics, an NSF-CAREER under
award #1846369, DARPA Learning with Less Labels (LwLL) and FastNICS programs,
and NSF-CIF awards #1813877 and #2008443. This research is also in part
supported by AWS credits through an Amazon Faculty research award.

## References

  * Candes et al. (2006) Emmanuel J. Candes, Justin K. Romberg, and Terence Tao.  Stable signal recovery from incomplete and inaccurate measurements.  _Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences_ , 59(8):1207–1223, 2006.
  * Cao et al. (2021a) Hu Cao, Yueyue Wang, Joy Chen, Dongsheng Jiang, Xiaopeng Zhang, Qi Tian, and Manning Wang.  Swin-Unet: Unet-like pure Transformer for medical image segmentation.  _arXiv:2105.05537_ , May 2021a.  arXiv: 2105.05537.
  * Cao et al. (2021b) Jiezhang Cao, Yawei Li, Kai Zhang, and Luc Van Gool.  Video super-resolution Transformer.  _arXiv preprint arXiv:2106.06847_ , 2021b.
  * Carion et al. (2020) Nicolas Carion, Francisco Massa, Gabriel Synnaeve, Nicolas Usunier, Alexander Kirillov, and Sergey Zagoruyko.  End-to-end object detection with Transformers.  In _European Conference on Computer Vision_ , pages 213–229. Springer, 2020.
  * Chen et al. (2021) Hanting Chen, Yunhe Wang, Tianyu Guo, Chang Xu, Yiping Deng, Zhenhua Liu, Siwei Ma, Chunjing Xu, Chao Xu, and Wen Gao.  Pre-trained image processing Transformer.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 12299–12310, 2021.
  * (6) Joseph Y. Cheng.  Stanford 2D FSE.  http://mridata.org/list?project=Stanford 2D FSE.
  * Çiçek et al. (2016) Özgün Çiçek, Ahmed Abdulkadir, Soeren S Lienkamp, Thomas Brox, and Olaf Ronneberger.  3D U-Net: Learning dense volumetric segmentation from sparse annotation.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , pages 424–432. Springer, 2016.
  * Donoho (2006) David L. Donoho.  Compressed sensing.  _IEEE Transactions on Information Theory_ , 52(4):1289–1306, 2006.
  * Dosovitskiy et al. (2020) Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit, and Neil Houlsby.  An image is worth 16x16 words: Transformers for image recognition at scale.  _arXiv:2010.11929 [cs]_ , 2020.
  * Feng et al. (2021) Chun-Mei Feng, Yunlu Yan, Huazhu Fu, Li Chen, and Yong Xu.  Task Transformer network for joint MRI reconstruction and super-resolution.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , pages 307–317, 2021.
  * Hammernik et al. (2018) Kerstin Hammernik, Teresa Klatzer, Erich Kobler, Michael P Recht, Daniel K Sodickson, Thomas Pock, and Florian Knoll.  Learning a variational network for reconstruction of accelerated MRI data.  _Magnetic Resonance in Medicine_ , 79(6):3055–3071, 2018.
  * Hammernik et al. (2019) Kerstin Hammernik, Jo Schlemper, Chen Qin, Jinming Duan, Ronald M Summers, and Daniel Rueckert.  $\Sigma$-net: Systematic evaluation of iterative deep neural networks for fast parallel MR image reconstruction.  _arXiv preprint arXiv:1912.09278_ , 2019.
  * Han and Ye (2018) Yoseob Han and Jong Chul Ye.  Framing U-Net via deep convolutional framelets: Application to sparse-view CT.  _IEEE Transactions on Medical Imaging_ , 37(6):1418–1429, 2018.
  * He et al. (2015) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  _arXiv:1512.03385_ , 2015.
  * Heo et al. (2021) Byeongho Heo, Sangdoo Yun, Dongyoon Han, Sanghyuk Chun, Junsuk Choe, and Seong Joon Oh.  Rethinking spatial dimensions of Vision Transformers.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 11936–11945, 2021.
  * Huang et al. (2022) Jiahao Huang, Yingying Fang, Yinzhe Wu, Huanjun Wu, Zhifan Gao, Yang Li, Javier Del Ser, Jun Xia, and Guang Yang.  Swin Transformer for fast MRI.  _arXiv preprint arXiv:2201.03230_ , 2022.
  * Huang et al. (2021) Xiaohong Huang, Zhifang Deng, Dandan Li, and Xueguang Yuan.  Missformer: An effective medical image segmentation Transformer.  _arXiv preprint arXiv:2109.07162_ , 2021.
  * Hyun et al. (2018) Chang Min Hyun, Hwa Pyung Kim, Sung Min Lee, Sungchul Lee, and Jin Keun Seo.  Deep learning for undersampled MRI reconstruction.  _Physics in Medicine & Biology_, 63(13):135007, 2018.
  * Klug and Heckel (2022) Tobit Klug and Reinhard Heckel.  Scaling laws for deep learning based image reconstruction.  _arXiv preprint arXiv:2209.13435_ , 2022.
  * Korkmaz et al. (2022) Yilmaz Korkmaz, Salman UH Dar, Mahmut Yurt, Muzaffer Özbey, and Tolga Cukur.  Unsupervised MRI reconstruction via zero-shot learned adversarial Transformers.  _IEEE Transactions on Medical Imaging_ , 2022.
  * Liang et al. (2021) Jingyun Liang, Jiezhang Cao, Guolei Sun, Kai Zhang, Luc Van Gool, and Radu Timofte.  SwinIR: Image restoration using Swin Transformer.  _arXiv:2108.10257_ , 2021.
  * Lin and Heckel (2022) Kang Lin and Reinhard Heckel.  Vision Transformers enable fast and robust accelerated MRI.  In _Medical Imaging with Deep Learning_ , 2022.
  * Liu et al. (2021) Ze Liu, Yutong Lin, Yue Cao, Han Hu, Yixuan Wei, Zheng Zhang, Stephen Lin, and Baining Guo.  Swin Transformer: Hierarchical Vision Transformer using shifted windows.  _arXiv:2103.14030_ , 2021.
  * Lustig et al. (2008) Michael Lustig, David L. Donoho, Juan M. Santos, and John M. Pauly.  Compressed sensing MRI.  _IEEE Signal Processing Magazine_ , 25(2):72–82, 2008.
  * Luthra et al. (2021) Achleshwar Luthra, Harsh Sulakhe, Tanish Mittal, Abhishek Iyer, and Santosh Yadav.  Eformer: Edge enhancement based Transformer for medical image denoising.  _arXiv preprint arXiv:2109.08044_ , 2021.
  * Muckley et al. (2021) Matthew J Muckley, Bruno Riemenschneider, Alireza Radmanesh, Sunwoo Kim, Geunu Jeong, Jingyu Ko, Yohan Jun, Hyungseob Shin, Dosik Hwang, Mahmoud Mostapha, et al.  Results of the 2020 fastMRI challenge for machine learning MR image reconstruction.  _IEEE transactions on Medical Imaging_ , 40(9):2306–2317, 2021.
  * Ongie et al. (2020) Gregory Ongie, Ajil Jalal, Christopher A. Metzler Richard G. Baraniuk, Alexandros G. Dimakis, and Rebecca Willett.  Deep learning techniques for inverse problems in imaging.  _IEEE Journal on Selected Areas in Information Theory_ , 2020.
  * Putzky et al. (2019) Patrick Putzky, Dimitrios Karkalousos, Jonas Teuwen, Nikita Miriakov, Bart Bakker, Matthan Caan, and Max Welling.  I-RIM applied to the fastMRI Challenge.  _arXiv:1910.08952_ , 2019.
  * Ramzi et al. (2020) Zaccharie Ramzi, Philippe Ciuciu, and Jean-Luc Starck.  XPDNet for MRI reconstruction: an application to the 2020 fastMRI Challenge.  _arXiv preprint arXiv:2010.07290_ , 2020.
  * Ronneberger et al. (2015) Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , pages 234–241, 2015.
  * Sawyer et al. (2013) Anne Marie Sawyer, Michael Lustig, Marcus Alley, Phdmartin Uecker, Patrick Virtue, Peng Lai, and Shreyas Vasanawala.  Creation of fully sampled MR data repository for compressed sensing of the knee.  2013\.
  * Sriram et al. (2020) Anuroop Sriram, Jure Zbontar, Tullie Murrell, Aaron Defazio, C. Lawrence Zitnick, Nafissa Yakubova, Florian Knoll, and Patricia Johnson.  End-to-end variational networks for accelerated MRI reconstruction.  In _Medical Image Computing and Computer Assisted Intervention_ , pages 64–73, 2020.
  * Sun et al. (2016) Jian Sun, Huibin Li, Zongben Xu, et al.  Deep ADMM-Net for compressive sensing MRI.  _Advances in Neural Information Processing Systems_ , 29, 2016.
  * Uecker et al. (2014) 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_ , 71(3):990–1001, 2014.
  * Wang et al. (2021a) Dayang Wang, Zhan Wu, and Hengyong Yu.  TED-net: Convolution-free T2T VisionTransformer-based encoder-decoder dilation network for low-dose CT denoising.  In _International Workshop on Machine Learning in Medical Imaging_ , pages 416–425, 2021a.
  * Wang et al. (2021b) Wenhai Wang, Enze Xie, Xiang Li, Deng-Ping Fan, Kaitao Song, Ding Liang, Tong Lu, Ping Luo, and Ling Shao.  Pyramid Vision Transformer: A versatile backbone for dense prediction without convolutions.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 568–578, 2021b.
  * Wang et al. (2021c) Zhendong Wang, Xiaodong Cun, Jianmin Bao, and Jianzhuang Liu.  Uformer: A general U-shaped Transformer for image restoration.  _arXiv preprint arXiv:2106.03106_ , 2021c.
  * Wang et al. (2003) Zhou Wang, Eero P Simoncelli, and Alan C Bovik.  Multiscale structural similarity for image quality assessment.  In _The Thrity-Seventh Asilomar Conference on Signals, Systems & Computers, 2003_, volume 2, pages 1398–1402. Ieee, 2003.
  * Wu et al. (2022) Yixuan Wu, Kuanlun Liao, Jintai Chen, Danny Z Chen, Jinhong Wang, Honghao Gao, and Jian Wu.  D-Former: A U-shaped dilated Transformer for 3D medical image segmentation.  _arXiv preprint arXiv:2201.00462_ , 2022.
  * Zamir et al. (2021) Syed Waqas Zamir, Aditya Arora, Salman Khan, Munawar Hayat, Fahad Shahbaz Khan, and Ming-Hsuan Yang.  Restormer: Efficient Transformer for high-resolution image restoration.  _arXiv:2111.09881_ , 2021.
  * Zbontar et al. (2019) Jure Zbontar, Florian Knoll, Anuroop Sriram, Tullie Murrell, Zhengnan Huang, Matthew J. Muckley, Aaron Defazio, Ruben Stern, Patricia Johnson, Mary Bruno, Marc Parente, Krzysztof J. Geras, Joe Katsnelson, Hersh Chandarana, Zizhao Zhang, Michal Drozdzal, Adriana Romero, Michael Rabbat, Pascal Vincent, Nafissa Yakubova, James Pinkerton, Duo Wang, Erich Owens, C. Lawrence Zitnick, Michael P. Recht, Daniel K. Sodickson, and Yvonne W. Lui.  fastMRI: An open dataset and benchmarks for accelerated MRI.  _arXiv:1811.08839_ , 2019.
  * Zhang and Ghanem (2018) Jian Zhang and Bernard Ghanem.  ISTA-Net: Interpretable optimization-inspired deep network for image compressive sensing.  In _2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 1828–1837, 2018.
  * Zhang et al. (2021) Lipei Zhang, Zizheng Xiao, Chao Zhou, Jianmin Yuan, Qiang He, Yongfeng Yang, Xin Liu, Dong Liang, Hairong Zheng, Wei Fan, et al.  Spatial adaptive and Transformer fusion network (STFNet) for low-count PET blind denoising with MRI.  _Medical Physics_ , 2021.
  * Zhou et al. (2021) Hong-Yu Zhou, Jiansen Guo, Yinghao Zhang, Lequan Yu, Liansheng Wang, and Yizhou Yu.  nnFormer: Interleaved Transformer for volumetric segmentation.  _arXiv preprint arXiv:2109.03201_ , 2021.
  * Zhou et al. (2018) Zongwei Zhou, Md Mahfuzur Rahman Siddiquee, Nima Tajbakhsh, and Jianming Liang.  Unet++: A nested U-net architecture for medical image segmentation.  In _Deep Learning in Medical Image Analysis and Multimodal Learning for Clinical Decision Support_ , pages 3–11. 2018.

## Appendix

##  Appendix A HUMUS-Net baseline details

Our default model has $3$ RSTB-D downsampling blocks, $2$ RSTB-B bottleneck
blocks and $3$ RSTB-U upsampling blocks with $3-6-12$ attention heads in the
$D/U$ blocks and $24$ attention heads in the bottleneck block. For Swin
Transformers layers, the window size is $8$ for all methods and MLP ratio
($hidden\\_dim/input\\_dim$) of $2$ is used. Each RSTB block consists of $2$
STLs with embedding dimension of $66$. For HUMUS-Net-L, we increase the
embedding dimension to $96$. We use $8$ cascades of unrolling with a U-Net as
sensitivity map estimator (same as in E2E-VarNet) with $16$ channels.

We center crop and reflection pad the input images to $384\times 384$
resolution for HUMUS-Net and use the complete images for VarNet. In all
experiments, we minimize the SSIM loss between the target image $\bm{x}^{*}$
and the reconstruction $\hat{\bm{x}}$ defined as

$\mathcal{L}_{SSIM}(\bm{x}^{*},\hat{\bm{x}})=1-SSIM(\bm{x}^{*},\hat{\bm{x}}).$

fastMRI experiments– We train HUMUS-Net using Adam for $50$ epochs with a
learning rate of $0.0001$, dropped by a factor of $10$ at epoch $40$ and apply
adjacent slice reconstruction with $3$ slices. We run the experiments on
$8\times$ Titan RTX 24GB GPUs with a per-GPU batch size of $1$. We train
HUMUS-Net-L using a learning rate of $0.00007$ (using inverse square root
scaling with embedding dimension) for $45$ epochs, dropped by a factor of $10$
for further $5$ epochs. We trained HUMUS-Net-L on $4\times$ A100 40GB GPUs on
Amazon SageMaker with per-GPU batch size of 1.

Stanford 3D experiments– We train for $25$ epochs with a learning rate of
$0.0001$, dropped by a factor of $10$ at epoch $20$ and reconstruct single
slices. We run the experiments on $4\times$ Quadro RTX 5000 16GB GPUs with a
per-GPU batch size of $1$.

Stanford 2D experiments– We train for $50$ epochs with a learning rate of
$0.0001$, dropped by a factor of $10$ at epoch $40$ and reconstruct single
slices. Moreover, we crop reconstruction targets to fit into a $384\times 384$
box. We run the experiments on $8\times$ Quadro RTX 5000 16GB GPUs with a per-
GPU batch size of $1$.

##  Appendix B Ablation study experimental details

Here we discuss the experimental setting and hyperparameters used in our
ablation study in Section 5.2. In all experiments, we train the models on the
Stanford 3D dataset for $3$ different train-validation splits and report the
mean and standard error of the results. We use Adam optimizer with learning
rate $0.0001$ and train for $25$ epochs, decaying the learning rate by a
factor of $10$ at $20$ epochs. For all unrolled methods, we unroll $12$
iterations in order to provide direct comparison with the best performing
E2E-VarNet, which uses the same number of cascades. We reconstruct single
slices and do not use the method of adjacent slice reconstruction discussed in
Section 4.4. In models with sensitivity map estimator,we used the default
U-Net with $8$ channels (same as default E2E-VarNet). For Swin Transformers
layers, the window size is $8$ for all methods and MLP ratio
($hidden\\_dim/input\\_dim$) of $2$ is used. Further hyperparameters are
summarized in Table 5.

Method | Embedding dim. | # of STLs in RSTB | # attention heads | Patch size
---|---|---|---|---
Un-SS | $12$ | $2-2-2$ | $6-6-6$ | 1
Un-MS | $12$ | $D/U:~{}2-2-2,~{}B:~{}2$ | $D/U:~{}3-6-12,~{}B:~{}24$ | 1
Un-MS-Patch2 | $36$ | $D/U:~{}2-2-2,~{}B:~{}2$ | $D/U:~{}3-6-12,~{}B:~{}24$ | 2
HUMUS-Net | $36$ | $D/U:~{}2-2-2,~{}B:~{}2$ | $D/U:~{}3-6-12,~{}B:~{}24$ | 1
SwinIR | $66$ | $6-6-6-6$ | $6-6-6-6$ | 1

Table 5: Ablation study experimental details. We show the number of STL layers
per RSTB blocks and number of attention heads for multi-scale networks in
downsampling (D) , bottleneck (B) and upsampling (U) paths separately.

We run the experiments on $4\times$ Quadro RTX 5000 16GB GPUs with a per-GPU
batch size of $1$.

##  Appendix C Results on additional accelerations

We observe that HUMUS-Net achieves state-of-the-art performance across a wide
range of acceleration ratios. We perform experiments on the Stanford 3D
dataset using a small HUMUS-Net model with only $6$ cascades, embedding
dimension of $66$, adjacent slice reconstruction with $3$ slices and for the
sake of simplicity we removed the residual path from the HUMUS-Block. We call
this model HUMUS-Net-S. We set the learning rate to $0.0002$ and train the
model with Adam optimizer. We generate $3$ random training-validation splits
on the Stanford 3D dataset (same splits as in other Stanford 3D experiments in
this paper), and show the mean and standard error of the results in Table 6.
In our experiments, HUMUS-Net achieved higher quality reconstructions in every
metric over E2E-VarNet.

Acceleration | Model | SSIM($\uparrow$) | PSNR($\uparrow$) | NMSE($\downarrow$)
---|---|---|---|---
4x | E2E-VarNet | $0.9623\pm 0.0038$ | $42.9\pm 0.5$ | $0.0103\pm 0.0001$
HUMUS-Net-S | $\textbf{0.9640}\pm\textbf{0.0040}$ | $\textbf{43.3}\pm\textbf{0.5}$ | $\textbf{0.0096}\pm\textbf{0.0002}$
8x | E2E-VarNet | $0.9432\pm 0.0063$ | $40.0\pm 0.6$ | $0.0203\pm 0.0006$
HUMUS-Net-S | $\textbf{0.9459}\pm\textbf{0.0065}$ | $\textbf{40.4}\pm\textbf{0.7}$ | $\textbf{0.0184}\pm\textbf{0.0008}$
12x | E2E-VarNet | $0.9259\pm 0.0084$ | $37.7\pm 0.7$ | $0.0347\pm 0.0018$
HUMUS-Net-S | $\textbf{0.9283}\pm\textbf{0.0052}$ | $\textbf{38.0}\pm\textbf{0.5}$ | $\textbf{0.0298}\pm\textbf{0.0019}$

Table 6:  Experiments on various acceleration factors on the Stanford 3D
dataset. Mean and standard error of $3$ random training-validation splits is
shown.

##  Appendix D Effect of the number of unrolled iterations

The number of cascades in unrolled networks has a fundamental impact on their
performance. Deeper networks typically perform better, but also incur heavy
computational and memory cost. In this experiment we investigate the scaling
of HUMUS-Net with respect to the number of iterative unrollings.

We perform an ablation study on the Stanford 3D dataset with $\times 8$
acceleration for a fixed training-validation split, where 20% of the training
set has been set aside for validation. We reconstruct single slices, without
applying adjacent slice reconstruction. We plot the highest SSIM throughout
training on the validation dataset for each network on Figure 6.

We observe improvements in reconstruction performance with increasing number
of cascades, and this improvement has not saturated yet in the range of model
sizes we have investigated. In our main experiments, we select 8 cascades due
to memory and compute limitations for deeper models. However, our experiment
suggests that HUMUS-Net can potentially obtain even better reconstruction
results given enough computational resources.

Figure 6: Validation SSIM as a function of number of cascades (unrolled
itarations) in HUMUS-Net on the Stanford 3D dataset. We observe a steady
increase in reconstruction performance with more cascades.

##  Appendix E Effect of Adjacent Slice Reconstruction

In this section, we demonstrate that even though adjacent slice reconstruction
is in general helpful, it is not the main reason why HUMUS-Net performs better
than E2E-VarNet.

To this end, we add adjacent slice reconstruction to E2E-VarNet and
investigate its effect on model performance. We run experiments on the
Stanford 3D dataset for $3$ different random train-validation splits. We use
adjacent slice reconstruction with $3$ slices. We have found that the default
learning rate for E2E-VarNet is not optimal with the increased input size,
therefore we tune the learning rate using grid search and set it to $0.0005$.
We compare the results with HUMUS-Net, where we match the number of cascades
in E2E-VarNet, and use the default embedding dimension of $66$. We do not use
adjacent slice reconstruction when training HUMUS-Net in order to ablate its
effect.

The results are summarized in Table 7. We observe that ASR boosts the
reconstruction quality of E2E-VarNet. However, ASR alone cannot close the gap
between the two models, as HUMUS-Net without ASR still outperforms the best
E2E-VarNet model with ASR.

Model | ASR? | SSIM | PSNR | NMSE
---|---|---|---|---
E2E-VarNet | ✗ | $0.9432\pm 0.0063$ | $40.0\pm 0.6$ | $0.0203\pm 0.0006$
E2E-VarNet | ✓ | $0.9457\pm 0.0064$ | $40.4\pm 0.7$ | $0.0186\pm 0.0007$
HUMUS-Net | ✗ | $\textbf{0.9467}\pm\textbf{0.0059}$ | $\textbf{40.6}\pm\textbf{0.6}$ | $\textbf{0.0178}\pm\textbf{0.0005}$
Table 7: Results of the adjacent slice reconstruction ablation study on the
Stanford 3D dataset. Mean and standard error over $3$ random train-validation
splits is shown. ASR improves the performance of E2E-VarNet. However, HUMUS-
Net outperforms E2E-VarNet in all cases even without ASR.

##  Appendix F Iterative denoising visualization

Here, we provide more discussion on the iterative denoising interpretation of
unrolled networks, such as HUMUS-Net. Consider the regularized inverse problem
formulation of MRI reconstruction as

$\bm{\hat{x}}=\operatorname*{arg\,min}_{\bm{x}}\left\|{\mathcal{A}\left(\bm{x}\right)}-\bm{\tilde{k}}\right\|^{2}+\mathcal{R}(\bm{x}),$
(F.1)

and design the architecture based on unrolling the gradient descent steps on
the above problem. Applying data consistency corresponds to the gradient of
the first loss term in (F.1), whereas we learn the gradient of the regularizer
parameterized by the HUMUS-Block architecture yielding the update rule in
k-space

$\bm{\hat{k}^{t+1}}=\bm{\hat{k}^{t}}-\mu^{t}{\bm{M}}(\bm{\hat{k}^{t}}-\bm{\tilde{k}})+\mathcal{G}(\bm{\hat{k}^{t}}).$

This can be conceptualized as alternating between enforcing consistency with
the measurements and applying a denoiser based on some prior, represented by
the HUMUS-Block in our architecture.

We show an example of a sequence of intermediate reconstructions in order to
support the denoising intuition in Figure 7. We plot the magnitude image of
reconstructions at the outputs of consecutive HUMUS-Blocks, along with the
model input zero-padded reconstruction. We point out that in general the
intermediate reconstruction quality (for instance in terms of SSIM) is not
necessarily an increasing function of cascades. This is due to the highly non-
linear nature of the mapping represented by the neural network, which might
not always be intuitive and interpretable.

Figure 7: Visualization of intermediate reconstructions in HUMUS-Net.

##  Appendix G Vision Transformer terminology overview

Here we provide a brief overview of the terms related to the Transformer
architecture, some of which has been carried over from the Natural Language
Processing (NLP) literature. We introduce the key concepts through the Vision
Transformer architecture, which has inspired most Transformer-based
architectures in computer vision.

Traditional Transformers for NLP receive a sequence of 1D token embeddings.
Each token may represent for instance a word in a sentence. To extend this
idea to 2D images, the input $x\in\mathbb{R}^{H\times W\times C}$ image ($H$
and $W$ stand for spatial dimensions, $C$ for channel dimension) is split into
$N$ patches of size $P\times P$, and each patch is flattened into a 1D vector
to produce an input of flattened patches $x_{f}$ with shape $N\times P^{2}C$:

$x_{f}=\text{PatchFlatten}(x).$

In order to set the latent dimension of the input entering the Vision
Transformer, a trainable linear mapping $E\in\mathbb{R}^{P^{2}C\times D}$ is
applied to the flattened patches, resulting in the so called patch embeddings
$z$ of shape $N\times D$:

$z=x_{f}E.$

In order to encode information with respect to the position of image patches,
a learnable positional embedding $E_{pos}\in\mathbb{R}^{N\times D}$ is added
to the patch embeddings:

$z_{0}=z+E_{pos}.$

The input to the Transformer encoder is this $N\times D$ representation, which
we also refer to in the paper as token representation, as each row in the
representation corresponds to a token (in our case an image patch) in the
original input.

##  Appendix H Detailed validation results

Dataset | Model | SSIM($\uparrow$) | PSNR($\uparrow$) | NMSE($\downarrow$)
---|---|---|---|---
fastMRI | E2E-VarNet | 0.8908 | 36.8 | 0.0092
HUMUS-Net | 0.8934 | 37.0 | 0.0090
Stanford 2D | E2E-VarNet | $0.8928\pm 0.0168$ | $\textbf{33.9}\pm\textbf{0.7}$ | $0.0339\pm 0.0037$
HUMUS-Net | $\textbf{0.8954}\pm\textbf{0.0136}$ | $33.7\pm 0.6$ | $\textbf{0.0337}\pm\textbf{0.0024}$
Stanford 3D | E2E-VarNet | $0.9432\pm 0.0063$ | $40.0\pm 0.6$ | $0.0203\pm 0.0006$
HUMUS-Net | $\textbf{0.9453}\pm\textbf{0.0065}$ | $\textbf{40.4}\pm\textbf{0.6}$ | $\textbf{0.0187}\pm\textbf{0.0009}$

Table 8:  Detailed validation results of HUMUS-Net on various datasets. For
datasets with multiple train-validation split runs we show the mean and
standard error of the runs.

##  Appendix I Additional figures

Figure 8: Swin Transformer Layer, the fundamental building block of the
Residual Swin Transformer Block.

Figure 9: Depiction of iterative unrolling with sensitivity map estimator
(SME). HUMUS-Net applies a highly efficient denoiser to progressively improve
reconstructions in a cascade of sub-networks.

Figure 10: Patch merge and expand operations used in MUST.

Figure 10:
Figure 10:

Figure 10:
Figure 10:

Figure 10:
Figure 10:

Figure 10:
Figure 10:

Figure 10:
Figure 10: Ground truth

Figure 10: E2E-VarNet
Figure 10: HUMUS-Net (ours)

Figure 11: Visual comparison of reconstructions from the fastMRI knee dataset
with $\times 8$ acceleration. HUMUS-Net reconstructs fine details on MRI
images that other state-of-the-art methods may miss.


# Some notes on Goodman’s marginal-free correspondence analysis

Choulakian V., Université de Moncton, Canada<EMAIL_ADDRESS>
(February 2021)

###### Abstract

In his seminal paper Goodman (1996) introduced marginal-free correspondence
analysis; where his principal aim was to reconcile Pearson correlation measure
with Yule’s association measure in the analysis of contingency tables. We show
that marginal-free correspondence analysis is a particular case of
correspondence analysis with prespecified weights studied in the beginning of
the 1980s by Benzécri and his students. Furthermore, we show that it is also a
particular first-order approximation of logratio analysis with uniform
weights.

Key words: Marginal-free correspondence analysis; logratio analysis;
interactions; scale invariance; taxicab singular value decomposition.

AMS 2010 subject classifications: 62H25, 62H30

##  1 Introduction

Correspondence analysis (CA) and logratio analysis (LRA) are two popular
methods for the analysis and visualization of a contingency table (two-way
frequency counts data having $I$ rows and $J$ columns) or a compositional data
set ($I$ individuals, also named samples, of $J$ compositional parts). The
reference book on CA is Benzécri (1973); Beh and Lombardo (2014) present a
panoramic review of CA and its variants.

LRA includes two independently well developed methods: RC association models
for the analysis of contingency tables by Goodman (1979, 1981a, 1981b, 1991,
1996) and compositional data analysis (CoDA) by Aitchison (1986). CA and LRA
are based on three different principles: CA on Benzécri’s distribututional
equivalence principle, RC association models on Yule’s scale invariance
principle, and CoDA on Aitchison’s subcompositional coherence principle.

From a statistical point of view there is a fundamental difference between the
structures of a two-way contingency table $\mathbf{N}=(n_{ij})$ and a
compositional data set $\mathbf{X}=(x_{ij})$ for $i=1,...,I$ and $j=1,...,J$;
while from a mathematical point of view the form of the resulting statistical
equations arising from different departure assumptions may be identical in
Goodman’s RC association models and Aitchison’s CoDA.

Goodman (1996, equation (46)) in his seminal paper introduced marginal-free
correspondence analysis (mfCA); where his principal aim was to reconcile
Pearson correlation measure with Yule’s association measure in the analysis of
contingency tables. In this paper, we show that mfCA is a particular case of
CA with prespecified weights, which has been studied in the beginning of 1980s
under the direction of Benzécri. In Benzécri’s edited journal Les Cahiers de
l’Analyse des Données, the following papers appeared [Madre (1980), Cholakian
(1980, 1984), Benzécri (1983a, 1983b), Benzécri et al. (1980) and Moussaoui
(1987)]. Furthermore, we show that mfCA is also a particular first-order
approximation of LRA analysis with uniform weights.

This paper is organized as follows: Section 2 presents three different basic
ways of representing the concept of interaction in a contingency table;
section 3 discusses the the important consequences of Yule’s scale invariance
association index; section 4 presents Goodman’s marginial-free CA; section 5
discusses an example; section 6 presents the R code to do the computations;
finally we conclude in section 7.

##  2 Preliminaries on analysis of contingency tables

Let $\mathbf{P=N/}n=(p_{ij})$ of size $I\times J$ be the associated
correspondence matrix (probability table) of a contingency table N. We define
as usual $p_{i+}=\sum_{j=1}^{J}p_{ij}$ , $p_{+j}=\sum_{i=1}^{I}p_{ij},$ the
vector $\mathbf{r=(}p_{i+})\in\mathbb{R}^{I},$ the vector
$\mathbf{c=(}p_{+j})\in\mathbb{R}^{J}$, and $\mathbf{M}_{I}=Diag(\mathbf{r})$
the diagonal matrix having diagonal elements $p_{i+},$ and similarly
$\mathbf{M}_{J}=Diag(\mathbf{c}).$ We suppose that $\mathbf{M}_{I}$ and
$\mathbf{M}_{J}$ are positive definite metric matrices of size $I\times I$ and
$J\times J$, respectively; this means that the diagonal elements of
$\mathbf{M}_{I}$ and $\mathbf{M}_{J}$ are strictly positive.

###  2.1 Independence of the row and column categories

a) The $I$ row categories are independent of the $J$ column categories,

$\sigma_{ij}=p_{ij}-p_{i+}p_{+j}=0,$ (1)

where $\sigma_{ij}$ is the residual matrix of $p_{ij}$ with respect to the
independence model $p_{i+}p_{+j}.$

Remark 1: The contingency table $\mathbf{N}=(n_{ij})$ can also be represented
(coded) as an indicator matrix $\mathbf{Z=}\left[\mathbf{Z}^{I}\ \
\mathbf{Z}^{J}\right]=\left[\mathbf{(}z_{\alpha i})\ \ \mathbf{(}z_{\alpha
j})\right]\ $of size $n$ by $I+J,$ where $z_{\alpha i}=0$ if individual
$\alpha$ does not have level $i$ of the row variable, $z_{\alpha i}=1$ if
individual $\alpha$ has level $i$ of the row variable; $z_{\alpha j}=0$ if
individual $\alpha$ does not have level $j$ of the column variable, $z_{\alpha
j}=1$ if individual $\alpha$ has level $j$ of the column variable. Note that
$\mathbf{N}=(\mathbf{Z}^{I})^{\prime}\mathbf{Z}^{J}\ $and
$\sigma_{ij}=p_{ij}-p_{i+}p_{+j}$ is the covariance between the $i$-th column
of $\mathbf{Z}^{I}$ and the $j$-th column of $\mathbf{Z}^{J}$.

b) The independence assumption $\sigma_{ij}=0$ can also be interpreted in
another way as

$\displaystyle\Delta_{ij}$ $\displaystyle=(\frac{p_{ij}}{p_{i+}p_{+j}}-1)=0$
(2) $\displaystyle=\frac{1}{p_{i+}}(\frac{p_{ij}}{p_{+j}}-p_{i+})=0$
$\displaystyle=\frac{1}{p_{+j}}(\frac{p_{ij}}{p_{i+}}-p_{+j});$

this is the column and row homogeneity models. Benzécri (1973, p.31) named the
conditional probability vector ($\frac{p_{ij}}{p_{+j}}$ for $i=1,...,I$ and
$j$ fixed) the profile of the $j$th column; and the element
$\frac{p_{ij}}{p_{i+}p_{+j}}$ the density function of the probability measure
$(p_{ij})$ with respect to the product measure $p_{i+}p_{+j}$. The element
$\frac{p_{ij}}{p_{i+}p_{+j}}$is named Pearson ratio in Goodman (1996) and Beh
and Lombardo (2014, p.123).

c) A third way to represent the indepence assumption $\sigma_{ij}=0$ and the
row and column homogeneity models $\Delta_{ij}=0$ is via the ($w_{i}^{R}$,
$w_{j}^{C})$ weighted loglinear formulation, equation (3), assuming $p_{ij}>0$
and defining $G_{ij}=\log(p_{ij}),$

$\displaystyle\lambda_{ij}$ $\displaystyle=0$ (3)
$\displaystyle=G_{ij}-G_{i+}-G_{+j}+G_{++},$

where $G_{i+}=\sum_{j=1}^{J}G_{ij}w_{j}^{C},$
$G_{+j}=\sum_{i=1}^{I}G_{ij}w_{i}^{R}$ and
$G_{++}=\sum_{j=1}^{J}\sum_{i=1}^{I}G_{ij}w_{j}^{C}w_{i}^{R}$; $w_{j}^{C}>0$
and $w_{i}^{R}>0,$ satisfying
$\sum_{j=1}^{J}w_{j}^{C}=\sum_{i=1}^{I}w_{i}^{R}\ =1,$ are a priori fixed
probability weights. Two popular weights are marginal ($w_{j}^{C}=p_{+j}$,
$w_{i}^{R}=p_{i+})$ and uniform ($w_{j}^{C}=1/J$,$w_{i}^{R}=1/I).$ This is
implicit in equation 7 in Goodman (1996) or equation 2.2.6 in Goodman (1991);
and explicit in Egozcue et al. (2015).

Equation (3) is equivalent to the logratios

$\log(\frac{p_{ij}p_{i_{1}j_{1}}}{p_{ij_{1}}p_{i_{1}j}})=0\text{ \ \ for
}i\neq i_{1}\text{ and }j=j_{1},$

which Goodman (1979, equation 2.2) names it ”null association” model.

Equation (3) is also equivalent to

$p_{ij}=\frac{\exp(G_{i+})\exp(G_{+j})}{\exp(G_{++})},$

from which we deduce that : under the independence assumption the marginal row
probability vector ($p_{i+})$ is proportional to the vector of weighted
geometric means ($\exp(G_{i+}));$ and a similar property is true also for the
columns; see for instance Egozcue et al. (2015).

###  2.2 Interaction factorization

Suppose the independence-homogeneity-null association models are not true,
then each of the three equivalent model formulations (1,2,3) can be
generalized to explain the nonindependence-nonhomogeneity-association, named
interaction, among the $I$ rows and the $J$ columns by adding $k$ bilinear
terms, where $k=rank(\mathbf{N)-}1$. We designate any one of the interaction
indices (1,2,3) by $\tau_{ij}.$

Benzécri (1973, Vol.1, p. 31-32) emphasized the importance of row and column
weights or metrics in multidimensional data analysis; this is the reason in
the french data analysis circles any study starts with a triplet
$\mathbf{(X,M}_{I},\mathbf{M}_{J}\mathbf{)}$, where X represents the data set,
$\mathbf{M}_{I}=(Diag(m_{i}^{r}))$ is the metric defined on the rows and
$\mathbf{M}_{J}=(Diag(m_{j}^{c}))$ the metric defined on the columns. We
follow the same procedure where:

a) In covariance analysis, $\mathbf{X}=(\tau_{ij}\mathbf{)=(}\sigma_{ij})$ and
$(\mathbf{M}_{I},\mathbf{M}_{J})=(Diag(1/I),Diag(1/J));$

b) In CA, $\mathbf{X}=(\tau_{ij}\mathbf{)=(}\Delta_{ij})$ and
$(\mathbf{M}_{I},\mathbf{M}_{J})=(Diag(p_{i+}),Diag(p_{+j}));$

c) In LRA, $\mathbf{X}=(\tau_{ij}\mathbf{)=(}\lambda_{ij})$ and
$(\mathbf{M}_{I},\mathbf{M}_{J})=(Diag(w_{i}^{R}),Diag(w_{j}^{C}))$ with
$\sum_{j=1}^{J}w_{j}^{C}=\sum_{i=1}^{I}w_{i}^{R}=1.$

We factorize the interactions in (1,2,3) by singular value decomposition (SVD)
or taxicab SVD (TSVD) as

$\tau_{ij}=\sum_{\alpha=1}^{k}f_{\alpha}(i)g_{\alpha}(j)/\delta_{\alpha}.$ (4)

In the SVD case the parameters $(f_{\alpha}(i),g_{\alpha}(j),\delta_{\alpha})$
satisfy the conditions: for $\alpha,\beta=1,...,k$

$\delta_{\alpha}^{2}=\sum_{\alpha=1}^{k}f_{\alpha}^{\ \
2}(i)m_{i}^{r}=\sum_{\alpha=1}^{k}g_{\alpha}^{2}(j)m_{j}^{c}$
$0=\sum_{\alpha=1}^{k}f_{\alpha}(i)m_{i}^{r}=\sum_{\alpha=1}^{k}g_{\alpha}(j)m_{j}^{c}$
$0=\sum_{\alpha=1}^{k}f_{\alpha}(i)f_{\beta}(i)m_{i}^{r}=\sum_{\alpha=1}^{k}g_{\alpha}(j)g_{\beta}(j)m_{j}^{c}$

In the TSVD case the parameters
$(f_{\alpha}(i),g_{\alpha}(j),\delta_{\alpha})$ satisfy the conditions: for
$\alpha,\beta=1,...,k$

$\delta_{\alpha}=\sum_{\alpha=1}^{k}|f_{\alpha}^{\
}|(i)m_{i}^{r}=\sum_{\alpha=1}^{k}|g_{\alpha}(j)|m_{j}^{c}$
$0=\sum_{\alpha=1}^{k}f_{\alpha}(i)m_{i}^{r}=\sum_{\alpha=1}^{k}g_{\alpha}(j)m_{j}^{c}$
$0=\sum_{\alpha=1}^{k}f_{\alpha}(i)\
sign(f_{\beta}(i))m_{i}^{r}=\sum_{\alpha=1}^{k}g_{\alpha}(j)\
sign(g_{\beta}(j))\ m_{j}^{c}\text{\ \ for }\alpha>\beta.$

A description of TSVD can be found, among others, in Choulakian (2006, 2016).

Remark 2

a) In the case $(\tau_{ij}\mathbf{)=(}\sigma_{ij})$, the bilinear
decomposition (4) is also named interbattery analysis first proposed by Tucker
(1958); later on, Tenenhaus and Augendre (1996) reintroduced it within
correspondence analysis circles, where they showed that the Tucker
decomposition by SVD produced on some correpodence tables more interesting
structure, more interpretable, than CA.

b) In the case $(\tau_{ij}\mathbf{)=(}\Delta_{ij})$, the CA decomposition has
many interpretations. Essentially, for data analysis purposes Benzécri (1973)
interpreted it as weighted principal components analysis of row and column
profiles. Another useful interpretation of CA, comparable to Tucker
interbattery analysis, is Hotelling(1936)’s canonical correlation analysis,
see Lancaster (1958) and Goodman (1991, 1996).

##  3 Yule’s principle of scale invariance

We start by quoting Goodman (1996, section 10) to really understand Yule’s
principle of scale invariance: ”Pearson’s approach to the analysis of cross-
classified data was based primarily on the bivariate normal. He assumed that
the row and column classifications arise from underlying continuous random
variables having a bivariate normal distribution, so that the sample
contingency table comes from a discretized bivariate normal; and he then was
concerned with the estimation of the correlation coefficient for the
underlying bivariate normal. On the other hand, Yule felt that, for many kinds
of contingency tables, it was not desirable in scientific work to introduce
assumptions about an underlying bivariate normal in the analysis of these
tables; and for such tables, he used, to a great extent, coefficients based on
the odds-ratios (for example, Yule’s Q and Y), coefficients that did not
require any assumptions about underlying distributions. The Pearson approach
and the Yule approach appear to be wholly different, but a kind of
reconciliation of the two perspectives was obtained in Goodman (1981a)”. An
elementary exposition of these ideas with examples can also be found in
Mosteller (1968).

In the notation of our paper, Goodman’s reconciliation is based on defining
the a priori weights in the association index (3),
$\lambda_{ij}=\lambda(p_{ij},w_{j}^{C},w_{i}^{R}),$ where by its decomposition
into bilinear terms, mwLRA will correspond to Pearson’s approach, while uwLRA
to Yule’s approach. Because log-odds

$\displaystyle\log(\frac{p_{ij}p_{i_{1}j_{1}}}{p_{ij_{1}}p_{i_{1}j}})$
$\displaystyle=\lambda_{ij}+\lambda_{{}_{i_{1}j_{1}}}-\lambda_{{}_{i_{1}j}}-\lambda_{ij_{1}}$
$\displaystyle=\sum_{\alpha=1}^{k}(f_{\alpha}(i)-f_{\alpha}(i_{1}))(g_{\alpha}(j)-g_{\alpha}(j_{1})/\delta_{\alpha}.$
(5)

To have a clear picture of LRA with general a priori prescribed weights
($w_{j}^{C},w_{i}^{R})$, we first study the properties of the association
index $\lambda_{ij},$ that distinguishes it from interaction indices (2,3).

###  3.1 Scale invariance of an interaction index

We are concerned with the property of scale dependence or independence of the
three interaction indices (1,2,3). We note that in (1,2,3), $p_{ij}$ depends
on $n_{ij},$ $p_{ij}=n_{ij}/\sum_{i,j}n_{ij}.$ To emphasize this dependence,
we express an interaction index by
$\tau_{ij}(n_{ij})=\tau(p_{ij},m_{i}^{R},m_{j}^{C})$ where: in the case of the
association index $\tau_{ij}(n_{ij})=\lambda_{ij}$ is defined in (3), in the
case of the nonhomogeneity index $\tau_{ij}(n_{ij})=\Delta_{ij}$ is defined in
(2), and in the case of the nonindependence index
$\tau_{ij}(n_{ij})=\sigma_{ij}$ is defined in (1). Following Yule (1912), we
state the following

Definition 1: An interaction index $\tau_{ij}(n_{ij})$ is scale invariant if
$\tau_{ij}(n_{ij})=\tau_{ij}(a_{i}n_{ij}b_{j})$ for scales $a_{i}>0$ and
$b_{j}>0$.

It is important to note that Yule’s principle of scale invariance concerns a
function of four interaction terms, see equation (5); while in Definition 1
the invariance concerns each interaction term.

It is evident that the interaction indices (1 and 2) are not scale invariant:
because they are marginal-dependent.

Concerning the association index (3) we have

Lemma 1: The association index (3) is scale invariant.

Proof: Let $n^{\ast}=\sum_{i,j}a_{i}n_{ij}b_{j},$ then

$\displaystyle\tau_{ij}(a_{i}n_{ij}b_{j})$
$\displaystyle=\lambda(a_{i}n_{ij}b_{j}/n^{\ast},w_{j}^{C},w_{i}^{R})$
$\displaystyle=\log(a_{i}n_{ij}b_{j}/n^{\ast})-\sum_{j=1}^{J}w_{j}^{C}\log(a_{i}n_{ij}b_{j}/n^{\ast})$
$\displaystyle-\sum_{i=1}^{I}w_{i}^{R}\log(a_{i}n_{ij}b_{j}/n^{\ast})+\sum_{j=1}^{J}\sum_{i=1}^{I}w_{j}^{C}w_{i}^{R}\log(a_{i}n_{ij}b_{j}/n^{\ast})$
$\displaystyle=\lambda(n_{ij},w_{j}^{C},w_{i}^{R})=\lambda(p_{ij},w_{j}^{C},w_{i}^{R})=\tau_{ij}(n_{ij})$
$\displaystyle=\lambda(a_{i}p_{ij}b_{j},w_{j}^{C},w_{i}^{R}).$ (6)

Lemma 2: To a first-order approximation,
$\lambda_{ij}\approx(\frac{p_{ij}}{w_{j}^{C}w_{i}^{R}}-\frac{p_{i+}}{w_{i}^{R}}-\frac{p_{+j}}{w_{j}^{C}}+1).$

Proof: The average value of the density function
$\frac{p_{ij}}{w_{j}^{C}w_{i}^{R}}$ with respect to the product measure
$w_{j}^{C}w_{i}^{R}$ is 1; so the $IJ$ values
$\frac{p_{ij}}{w_{j}^{C}w_{i}^{R}}$ are distributed around 1. By Taylor series
expansion of $\log x$ in the neighborhood of $x=1$, we have to a first-order
$\log x\approx x-1.$ Putting $a_{i}=1/w_{i}^{R}$ and $b_{j}=1/w_{j}^{C}$ in
(6), and by using
$\log(\frac{p_{ij}}{w_{j}^{C}w_{i}^{R}})\approx\frac{p_{ij}}{w_{j}^{C}w_{i}^{R}}-1$

$\displaystyle\lambda(p_{ij},w_{j}^{C},w_{i}^{R})$
$\displaystyle=\lambda(\frac{p_{ij}}{w_{j}^{C}w_{i}^{R}},w_{j}^{C},w_{i}^{R})$
$\displaystyle=\log(\frac{p_{ij}}{w_{j}^{C}w_{i}^{R}})-\sum_{j=1}^{J}w_{j}^{C}\log(\frac{p_{ij}}{w_{j}^{C}w_{i}^{R}})$
$\displaystyle-\sum_{i=1}^{I}w_{i}^{R}\log(\frac{p_{ij}}{w_{j}^{C}w_{i}^{R}})+\sum_{j=1}^{J}\sum_{i=1}^{I}w_{j}^{C}w_{i}^{R}\log(\frac{p_{ij}}{w_{j}^{C}w_{i}^{R}}$
$\displaystyle\approx\frac{p_{ij}}{w_{j}^{C}w_{i}^{R}}-1-(\frac{p_{i+}}{w_{i}^{R}}-1)-(\frac{p_{+j}}{w_{j}^{C}}-1)+0,$

which is the required result.

Remark 3: Lemma 2 provides a first order approximation to mwTLRA and uwTLRA,
where we see that both first-order approximations are marginal-dependent but
in different ways.

a) In the case $(a_{i},b_{j})=(1/p_{i+},1/p_{+j})$ and$\
(w_{j}^{C},w_{i}^{R})=(p_{+j},p_{i+})$ in Lemma 2,
$\lambda_{ij}=\lambda(p_{ij},p_{+j},p_{i+})\approx$
$\frac{p_{ij}}{p_{+j}p_{i+}}-1;$ which implies that CA (or TCA) is a first-
order approximation of mwLRA (or mwTLRA), a result stated in Cuadras et al.
(2006).

b) In the case $(a_{i},b_{j})=(I,J)$ and$\ (w_{j}^{C},w_{i}^{R})=(1/I,1/J)$ in
Lemma 2, $\lambda_{ij}=\lambda(p_{ij},1/J,1/I)\approx$
$IJp_{ij}-Ip_{i+}-Jp_{+j}+1;$ which implies that the bilinear expansion of the
right side by TSVD (or SVD) is a first-order approximation of uwTLRA (or
uwLRA).

In this subsection, we discussed the approximation of LRA to CA related
methods. Greenacre (2009) posed the reciprocal question: when CA related
methods converge to LRA? And he stated two results; which we provide a proof
in the following subsection.

###  3.2 Box-Cox transformation

Theoretically CA and LRA have been presented in a unified mathematical
framework via Box-Cox transformation by Goodman (1996), where the bilinear
terms have been estimated by SVD. Goodman’s framework was further considered,
among others, by Cuadras et al. (2006), Greenacre (2009, 2010), and Cuadras
and Cuadras (2015 ).

Consider the triplet (X, Q, D), where $\mathbf{X=(}x_{ij})$ with $x_{ij}>0$
represents the data set, and
$(\mathbf{Q,D})=(Diag(w_{i}^{R}),Diag((w_{j}^{C}))$ with
$\sum_{j=1}^{J}w_{j}^{C}=\sum_{i=1}^{I}w_{i}^{R}=1.$ Let $\alpha$ be a
nonnegative real number. Following Goodman (1996, equations (3,4,5)), we
define the interaction index,

$\displaystyle Int(\frac{x_{ij}^{\alpha}}{\alpha},w_{j}^{C},w_{i}^{R})$
$\displaystyle=\frac{x_{ij}^{\alpha}}{\alpha}-\sum_{j=1}^{J}w_{j}^{C}\frac{x_{ij}^{\alpha}}{\alpha}-\sum_{i=1}^{I}w_{i}^{R}\frac{x_{ij}^{\alpha}}{\alpha}+\sum_{j=1}^{J}\sum_{i=1}^{I}w_{j}^{C}w_{i}^{R}\frac{x_{ij}^{\alpha}}{\alpha}$
$\displaystyle=(\frac{x_{ij}^{\alpha}-1}{\alpha})-\sum_{j=1}^{J}w_{j}^{C}(\frac{x_{ij}^{\alpha}-1}{\alpha})$
$\displaystyle-\sum_{i=1}^{I}w_{i}^{R}(\frac{x_{ij}^{\alpha}-1}{\alpha})+\sum_{j=1}^{J}\sum_{i=1}^{I}w_{j}^{C}w_{i}^{R}(\frac{x_{ij}^{\alpha}-1}{\alpha}).$
(7)

Using the well-known result based on Hopital’s rule, lim${}_{\alpha\rightarrow
0}(\frac{x_{ij}^{\alpha}-1}{\alpha})=\log(x_{ij}),$ (7) converges to

$\lambda(x_{ij},w_{j}^{C},w_{i}^{R})=\log(x_{ij})-\sum_{j=1}^{J}w_{j}^{C}\log(x_{ij})-\sum_{i=1}^{I}w_{i}^{R}\log(x_{ij})+\sum_{j=1}^{J}\sum_{i=1}^{I}w_{j}^{C}w_{i}^{R}\log(x_{ij}).$
(8)

We consider two cases of (7, 8):

a)
$\lambda(x_{ij},w_{j}^{C},w_{i}^{R})=\lambda(p_{ij},p_{+j},p_{i+})=\lambda(p_{ij}/(p_{+j}p_{i+}),p_{+j},p_{i+}),$
which is the interaction term of mwLRA, and equivalent to Result 2 in
Greenacre (2010).

b)
$\lambda(x_{ij},w_{j}^{C},w_{i}^{R})=\lambda(p_{ij},1/J,1/I)=\lambda(IJp_{ij},1/J,1/I)=\lambda(p_{ij}/(p_{+j}p_{i+}),1/J,1/I),$
which is the interaction term of uwLRA; this is similar to Result 1 in
Greenacre (2010).

Equation (7) can also be applied differently, where:

In (7) we replace $w_{j}^{C}$ by$\
w_{j}^{C}(\alpha)=\frac{\sum_{i=1}^{I}x_{ij}^{\alpha}}{\sum_{i=1}^{I}\sum_{j=1}^{J}x_{ij}^{\alpha}},\
w_{i}^{R}$ by$\
w_{i}^{R}(\alpha)=\frac{\sum_{j=1}^{J}x_{ij}^{\alpha}}{\sum_{i=1}^{I}\sum_{j=1}^{J}x_{ij}^{\alpha}};$
we see that lim${}_{\alpha\rightarrow 0}w_{j}^{C}(\alpha)=\frac{I}{IJ}=1/J;$
similarly lim${}_{\alpha\rightarrow 0}w_{i}^{R}(\alpha)=1/I.$ Then we get

$\text{lim}_{\alpha\rightarrow
0}Int(x_{ij}^{\alpha},w_{j}^{C}(\alpha),w_{i}^{R}(\alpha))=\lambda(x_{ij},1/J,1/I).$

In particular lim${}_{\alpha\rightarrow
0}Int(x_{ij}^{\alpha},w_{j}^{C}(\alpha),w_{i}^{R}(\alpha))=\lambda(p_{ij},1/J,1/I)=\lambda(IJp_{ij},1/J,1/I)=\lambda(p_{ij}/(p_{+j}p_{i+}),1/J,1/I)$.

##  4 CA with prescribed weights and Goodman’s mfCA

CA with prescribed weights is done in two steps in the following way: We
observe a probability table $\mathbf{P}=(p_{ij})$ of size $I$ by $J$. Let
$\mathbf{Q}$ of size $I$ by $J$ be an unknown probability table with known
marginals $q_{i+}$ and $q_{+j}.$ The two steps are:

Step1: We construct Q which is in a sense ”nearest to $\mathbf{P}$”. Two
general criteria are:
$Int(q_{ij},q_{+j},q_{i+})=\lambda(q_{ij}=a_{i}p_{ij}b_{j},q_{+j},q_{i+})$
based on (3) and
min${}_{q_{ij}}\sum_{i,j}(\frac{q_{ij}}{q_{i+}q_{+j}}-\frac{p_{ij}}{p_{i+}p_{+j}})^{2}q_{i+}q_{+j}$
based on (2).

Step 2: We apply CA to the constructed probability Q

$\frac{q_{ij}-q_{i+}q_{+j}}{q_{i+}q_{+j}}=\sum_{\alpha=1}^{k}f_{\alpha}(i)g_{\alpha}(j)/\delta_{\alpha},$

which represents CA of $\mathbf{P}$ with prespecified weights
($q_{i+},q_{+j})$. Cholakian (1980) presents an example, where both criteria
have been applied and similar results have been obtained.

In the particular case, where
$Int(q_{ij},1/J,1/I)=\lambda(q_{ij}=a_{i}p_{ij}b_{j},1/J,1/I),$ we get
Goodman’s mfCA, see Goodman (1996, equation (46)). $\mathbf{Q=}(q_{ij})$ is
related to $\mathbf{P=(}p_{ij})$ via the strictly positive scales
($a_{i},b_{j}),$ that keeps the association between the $i$-th row and the
$j$-th column unchanged. The famous iterative proportional fitting algorithm
(IPF) is used to construct Q. That is, the constructed probability table
($q_{ij})$ has uniform marginals $q_{+j}=1/J$ and $q_{i+}=1/I.$ So in Step 2,
CA representation is

$IJq_{ij}-1=\sum_{\alpha=1}^{k}f_{\alpha}(i)g_{\alpha}(j)/\delta_{\alpha},$
(9)

which represents a first-order approximation to both uwLRA and mwLRA by Remark
3\. Furthermore, by Remark 2 we see that mfCA can be interpreted both as
Tucker and Hotelling decompositions.

##  5 Examples

We present the analysis of two datasets for comparative purposes.

###  5.1 Example 1

This dataset is contrived and taken from Goodman (1991, Table 10(1), that we
reproduce below

$\begin{bmatrix}4&10&1\\\ 10&50&10\\\ 1&10&4\end{bmatrix}$

According to Goodman, LRA has one principal dimension, while CA has 2
principal dimensions. Here we compare the dispersion results of the 4 methods:
CA, TCA, mfCA and mfTCA.

In CA: $\delta_{1}=corr(f_{1}(i),g_{1}(j)=0.20$ and
$\delta_{2}=corr(f_{2}(i),g_{2}(j)=0.048$

In mfCA: $\delta_{1}=corr(f_{1}(i),g_{1}(j)=0.41$ and
$\delta_{2}=corr(f_{2}(i),g_{2}(j)=0.050$

In TCA: $\ \ \ \ \delta_{1}=0.070$ and $\delta_{2}=0.034$

In mfTCA: $\delta_{1}=0.285$ and $\delta_{2}=0.043$

###  5.2 Example 2

We consider the rodent data set of size 28 by 9 found in TaxicabCA in R
package. This is an abundance data set of 9 kinds of rats in 28 cities in
California. It can be considered both a contingency table and a compositional
data set. Choulakian (2017) analyzed it by comparing the CA and TCA maps;
furthermore Choulakian (2021) showed that it has quasi-2-blocks structure.
Here we compare the dispersion results of the first 2 principal dimensions in
the 4 methods: CA, TCA, mfCA and mfTCA:

In CA: $\delta_{1}=corr(f_{1}(i),g_{1}(j)=0.864$ and
$\delta_{2}=corr(f_{2}(i),g_{2}(j)=0.678$

In mfCA:$\ \delta_{1}=corr(f_{1}(i),g_{1}(j)=0.827$ and
$\delta_{2}=corr(f_{2}(i),g_{2}(j)=0.679$

In TCA$:\ \ \ \ \delta_{1}=0.478$ and $\delta_{2}=0.422$

In mfTCA: $\delta_{1}=0.743$ and $\delta_{2}=0.541$

The curious reader can apply the R code below to campare the 4 maps: CA, mfCA,
TCA and mfTCA.

##  6 R code

#

# install packages

install.packages(c(”ipfr”, ”ca”, ”TaxicabCA”))

#

library(TaxicabCA)

dataMatrix = as.matrix(rodent)

nRow $<$\- nrow(dataMatrix)

nCol $<$\- ncol(dataMatrix)

ssize $<$\- sum(dataMatrix)

#

#Computation of Q matrix of rodent

library(ipfr)

mtx $<$\- dataMatrix

row_targets $<$\- rep(ssize/nRow, nRow)

column_targets $<$\- rep(ssize/nCol, nCol)

QMatrix $<$\- ipu_matrix(mtx, row_targets, column_targets)

rownames(QMatrix) $<$\- paste(””,1:nRow,sep=””)

colnames(QMatrix) $<$\- paste(”C”,1:nCol,sep=””)

#

#CA map of rodent dataset

library(ca)

plot(ca(dataMatrix))

#mfCA map of rodent

plot(ca(QMatrix))

#

# TCA map of rodent

tca.Data $<$\- tca(dataMatrix, nAxes=2,algorithm = ”exhaustive”)

plot(

tca.Data,

axes = c(1, 2),

labels.rc = c(1, 1),

col.rc = c(”blue”, ”red”),

pch.rc = c(5, 5, 0.3, 0.3),

mass.rc = c(F, F),

cex.rc = c(0.6, 0.6),

jitter = c(F, T)

)

#mfTCA map of rodent dimensions 1-2

tca.Data $<$\- tca(QMatrix, nAxes=2,algorithm = ”exhaustive”)

plot(

tca.Data,

axes = c(1, 2),

labels.rc = c(1, 1),

col.rc = c(”blue”, ”red”),

pch.rc = c(5, 5, 0.3, 0.3),

mass.rc = c(F, T),

cex.rc = c(0.6, 0.6),

jitter = c(F, T)

)

#mfTCA map of rodent dimensions 2-3

tca.Data $<$\- tca(QMatrix, nAxes=2,algorithm = ”exhaustive”)

plot(

tca.Data,

axes = c(2, 3),

labels.rc = c(1, 1),

col.rc = c(”blue”, ”red”),

pch.rc = c(5, 5, 0.3, 0.3),

mass.rc = c(F, T),

cex.rc = c(0.6, 0.6),

jitter = c(F, T)

)

##  7 Conclusion

In his seminal paper Goodman (1996) introduced marginal-free correspondence
analysis; where his principal aim was to reconcile Pearson correlation measure
with Yule’s association measure in the analysis of contingency tables. We
showed that marginal-free correspondence analysis is a particular case of
correspondence analysis with prespecified weights studied in the beginning of
the 1980s by Benzécri and his students. Furthermore, we showed that it is also
a particular first-order approximation of logratio analysis with uniform
weights.

    Acknowledgements.

Choulakian’s research has been supported by NSERC of Canada.

References

Aitchison J (1986) The Statistical Analysis of Compositional Data. London:
Chapman and Hall.

Beh E and Lombardo R (2014) Correspondence Analysis: Theory, Practice and New
Strategies. N.Y: Wiley.

Benzécri JP (1973) L’Analyse des Données: Vol. 2: L’Analyse des
Correspondances. Paris: Dunod.

Benzécri JP (1983a) Ajustement d’un tableau à des marges sous l’hypothèse
d’absence d’interaction ternaire. Les Cahiers de l’Analyse des Données, 8(2),
227-233

Benzécri JP (1983b) Sur une généralisation du problème de l’ajustement d’une
mesure à des marges. Les cahiers de l’analyse des données, 8(3), 359-370

Benzécri JP, Bourgarit C, Madre JL (980) Problème: ajustement d’un tableau à
ses marges d’après la formule de reconstitution. Les Cahiers de l’Analyse des
Données, 5(l), 163-172

Cholakian V (1980) Un exemple d’application de diverses méthodes d’ajustement
d’un tableau à des marges imposées. Les Cahiers de l’Analyse des Données,
5(2), 173-176

Cholakian V (1984) Méthodes et critères pour l’ajustement d’un tableau à des
marges imposées. Les Cahiers de l’Analyse des Données, 9(l), pp. 113-117

Choulakian V (2006) Taxicab correspondence analysis. Psychometrika, 71,
333-345.

Choulakian V (2016) Matrix factorizations based on induced norms. Statistics,
Optimization and Information Computing, 4, 1-14.

Choulakian V (2017) Taxicab correspondence analysis of sparse contingency
tables. Italian Journal of Applied Statistics, 29 (2-3), 153-179.

Choulakian V (2021) Quantification of intrinsic quality of a principal
dimension in correspondence analysis and taxicab correspondence analysis.
Available on arXiv:2108.10685.

Cuadras CM, Cuadras D, Greenacre M (2006) A comparison of different methods
for representing categorical data. Communications in Statistics-Simul. and
Comp, 35(2), 447-459.

Cuadras CM, Cuadras D (2015 ) A unified approach for the multivariate analysis
of contingency tables. Open Journal of Statistics, 5, 223-232

Egozcue JJ, Pawlowsky-Glahn V, Templ M, Hron K (2015) Independence in
contingency tables using simplicial geometry. Communications in Statistics -
Theory and Methods, 44:18, 3978-3996

Goodman LA (1979) Simple models for the analysis of association in cross-
classifications having ordered categories. Journal of the American Statistical
Association, 74,537-55

Goodman LA (1981a) Association models and the bivariate normal for contingency
tables with ordered categories. Biometrika, 68, 347-355

Goodman LA (1981b) Association models and canonical correlation in the
analysis of cross-classifications having ordered categories. Journal of the
American Statistical Association, 76, 320-334

Goodman, LA (1991) Measures, models, and graphical displays in the analysis of
cross-classified data. Journal of the American Statistical Association, 86
(4), 1085-1111

Goodman LA (1996) A single general method for the analysis of cross-classified
data: Reconciliation and synthesis of some methods of Pearson, Yule, and
Fisher, and also some methods of correspondence analysis and association
analysis. Journal of the American Statistical Association, 91, 408-428

Greenacre M (2009) Power transformations in correspondence analysis.
Computational Statistics & Data Analysis, 53(8), 3107-3116

Greenacre M (2010) Log-ratio analysis is a limiting case of correspondence
analysis. Mathematical Geosciences, 42, 129-134

Madre JL (1980) Méthodes d’ajustement d’un tableau à des marges. Les Cahiers
de l’Analyse des Données, 5 (1), 87-99

Mosteller F (1968) Association and estimation in contingency tables, Journal
of the American Statistical Association, 63, 1-28

Moussaoui AE (1987) Sur la reconstruction approchée d’un tableau de
correspondance à partir du tableau cumulé par blocs suivant suivant deux
partitions des ensembles I et J. Les Cahiers de l’Analyse des Données, 12(3),
365-370

Yule, G.U. (1912). On the methods of measuring association between two
attributes. JRSS, 75, 579-642