# Population Graph Cross-Network Node Classification for Autism Detection
Across Sample Groups

Anna Stephens, Francisco Santos, Pang-Ning Tan, and Abdol-Hossein Esfahanian
Dept. of Computer Science and Engineering
Michigan State University
East Lansing, United States
{steph496, santosf3, ptan<EMAIL_ADDRESS>
###### Abstract

Graph neural networks (GNN) are a powerful tool for combining imaging and non-
imaging medical information for node classification tasks. Cross-network node
classification extends GNN techniques to account for domain drift, allowing
for node classification on an unlabeled target network. In this paper we
present OTGCN, a powerful, novel approach to cross-network node
classification. This approach leans on concepts from graph convolutional
networks to harness insights from graph data structures while simultaneously
applying strategies rooted in optimal transport to correct for the domain
drift that can occur between samples from different data collection sites.
This blended approach provides a practical solution for scenarios with many
distinct forms of data collected across different locations and equipment. We
demonstrate the effectiveness of this approach at classifying Autism Spectrum
Disorder subjects using a blend of imaging and non-imaging data.

###### Index Terms:

autism, graph, transfer learning, cross-network node classification

##  I Introduction

Autism Spectrum Disorder (ASD) refers to a condition characterised by specific
communication impairments, restricted interests, and repetitive behaviours
[1]. As the disorder typically presents itself early in life, early and
accurate detection can help reduce the severity of many lifelong symptoms.
Towards this end, automated techniques based on machine learning and deep
learning have been developed for the early detection of ASD. Specifically,
these techniques have been applied to a wide range of subject data including
detailed subject screening data [2], videos of subject movements [3], and
Magnetic Resonance Imaging (MRI) data [4]. Given the diverse modalities of
data available, it is likely that the best detection results can be achieved
by combining the data from different sources.

Previous works have established an effective method to combine various
modalities of data using graph structures because of their flexible yet
powerful representation [5] [6]. These approaches make use of population
graphs, where the nodes represent individuals and edges are generally defined
with some similarity measure. In several prior ASD research, the node features
are generated from image information while the edges consider a combination of
subject information (sex, age, etc) in addition to image similarity [7] [8]
[9] [10]. Once these graphs are constructed, detecting ASD becomes a _node
classification_ problem.

Figure 1: An example illustrating the perils of failing to account for concept
drift in multi-source data. Assume a binary node classification task where the
two classes are represented as x’s and o’s, respectively. The blue points
denote the dataset for the source domain while the orange points denote the
dataset for the target domain. Observe that a decision boundary constructed
from the blue dataset will have reduced accuracy when applied to the orange
dataset due to concept drift between the source and target networks.

Similar to other diagnosis problems, one of the major challenges in applying
machine learning to ASD detection is the limited amount of labeled data
available. This has led to growing interest in utilizing labeled data from
multiple sites to train the machine learning model. One limitation of these
approaches is that they fail to consider a more realistic scenario in which
there may be little or no labeled data associated with the population of
interest. For example, a model may be trained on a research dataset but needs
to be applied to another location with a different imaging equipment or
patient demographics that were not well represented in the training data. This
problem can be addressed using an approach generally known as domain
adaptation or transfer learning [11]. Within the context of node
classification tasks, it is also referred to as _cross-network node
classification_ (CNNC) approach. CNNC assumes the availability of a
sufficiently large number of labeled nodes in a source network and an
unlabeled target network, whose node labels are to be predicted accurately.

This paper focuses on addressing the CNNC task for ASD detection. We use the
popular ABIDE [12] dataset, which contains both resting-state functional
Magnetic Resonance Imaging (fMRI) information and phenotypic data such as age,
sex, and screening results. Samples in the ABIDE dataset were collected from
several different collection sites, which were then divided into 2 groups to
form the source and target networks for our experimental studies.

CNNC has two major challenges when applied to the ABIDE dataset. The first is
extracting relevant information from the networks for ASD detection. In this
work we will learn a graph embedding of the source and target networks via a
graph convolutional network (GCN) [13]. A GCN layer is capable of learning a
node embedding which contains information about a node and its immediate
neighbors. This method relies on the homophily principle [14], which states
that similar individuals tend to form neighborhoods within a graph. Therefore,
in a graph with high homophily we can improve a node embedding by adding
information from it’s neighborhood.

The second major challenge is to handle potential concept drift between the
source and target networks. As the separate data collection sites may have
different fMRI imaging equipment and procedural differences, this may
introduce some discrepancies or ”drift” between the two networks. The presence
of concept drift often leads to poor results if a classifier is simply trained
on the source dataset and applied as it is to the target dataset. Fig. 1
demonstrates the challenge of CNNC when the learned node embedding does not
account for such distributional shift. A decision boundary learned from the
source dataset is likely to incorrectly classify a significant portion of the
target dataset.

In this paper we introduce OTGCN, a novel approach to address the CNNC task
for ASD detection. Our proposed approach combines graph neural network with
optimal transport (OT) to handle the drift between the source and target
networks. OT is a method for mapping a transportation between two
distributions. We will use OT to map the learned source representation to the
target representation. This strategy allows us to train an accurate model for
classifying the target nodes. Experimental results on the ABIDE dataset
demonstrate the effectiveness of our approach at diagnosing ASD across
different collection sites compared to state of the art CNNC methods.

##  II Related Works

###  II-A Machine Learning Approaches to ASD Detection

ASD detection is a task that lends itself to a wide variety of approaches.
Zunino et al [3] employed computer vision approaches for automatic early
detection of ASD by evaluating videos of subject movements when grasping a
bottle. They then applied recurrent neural networks to distinguish subjects
who have ASD from those who do not. Erkan and Thanh [2] analyzed detailed
screening data collected from mobile app surveys using traditional machine
learning methods such as k-nearest neighbor, support vector machines, and
random forests to diagnose subjects with ASD. Yuan et al. [15] applied natural
language processing (NLP) techniques to analyze hand written medical forms of
potential ASD patients while Carette et al. [16] performed ASD detection on
eye tracking data using long-short term memory (LSTM) networks.

The majority of the recent works in this area, however, focuses on using MRI
data [4] [17] along with other subject information such as sex, age, and
screening results [7] [8] [10] [9] for screening potential ASD patients. For
example, Li et al. [4] presented a graph neural network approach to find
biomarkers that can be used to detect ASD while Parisot et al. [7] employed a
graph convolutional network (GCN) to perform the detection using both fMRI
imaging and non-imaging phenotypic data. Similar to these works, our paper
focuses on using a combination of fMRI and other subject data, though our
approach is also applicable to blend other forms of data given the
representation power of graph neural networks.

###  II-B Machine Learning on Multi-site fMRI Data

Previous works on diagnosing brain-related problems using fMRI data from
multiple sites generally fall into two major categories—transfer learning and
federated learning.

Transfer learning [11] is a machine learning approach that enables prediction
models trained from a given data domain (known as the _source_ domain) to be
applied to another domain (known as the _target_ domain). Such a domain
adaptation approach can be used even if the target domain has no labeled
training data. Previous works on applying transfer learning to fMRI data can
be found in [18] [19] [20]. However, these approaches are mostly designed for
using only the image information and do not consider more complex data
structures or the use of additional non-imaging information.

In contrast, federated learning [21] is designed for training prediction
models in a collaborative fashion on decentralized datasets. The approach
assumes restricted or indirect access to the source dataset and direct access
to a partially labeled target dataset. There has been a few works applying
federated learning approaches to multi-site fMRI data [22] [23], but none of
them consider non-imaging information.

###  II-C Cross-Network Node Classification (CNNC)

There are a handful of studies focusing on applying domain adaptation
approaches to multi-network data using graph neural networks. For node
classification, these approaches are also collectively known as cross-network
node classification (CNNC). For example, Shen et al. [24] presents an approach
called CDNE that uses maximum mean discrepancy (MMD) loss to learn graph
embeddings for the source and target networks before sending those embedding
to a shared classifier. The authors subsequently extended the approach to
ACDNE in [25] by adding an adversarial domain adaptation component to the
framework. AdaGCN [26] combines conventional graph convolutional network (GCN)
layers with adversarial domain adaptation, while UdaGCN [27] uses multiple
GCNs with both graph adjacency and positive point-wise mutual information
(PPMI) matrices to learn an improved combined feature embedding of the nodes.
A more recent work, ASN [28] combines the node embedding from UdaGCN with the
adversarial domain adaptation approach from AdaGCN to address the CNNC
problem. While all of these approaches were designed for CNNC tasks, they have
not been tested on the medical dataset nor with a combination involving
imaging and non-imaging data.

##  III Preliminaries

###  III-A Problem Statement

Let $\mathcal{G}(V,E,X,Y)$ be an attributed graph, with a node set $V$, edge
set $E\subseteq V\times V$, node feature matrix $X$, and node label $Y$. Let
$A$ denote the adjacency matrix representation of $E$, where $A_{ij}>0$ if
$(v_{i},v_{j})\in E$ and 0 otherwise. In a domain adaptation setting, we
assume there exists a source graph,
$\mathcal{G}_{s}(V_{s},E_{s},X_{s},Y_{s})$, where the node labels in $Y_{s}$
are known, and a target graph, $\mathcal{G}_{t}(V_{t},E_{t},X_{t},Y_{t})$,
where the node labels in $Y_{t}$ are unknown. The adjacency matrices
corresponding to the source and target graphs are denoted as
$A_{s}\in\mathbb{R}^{n_{s}\times n_{s}}$ and $A_{t}\in\mathbb{R}^{n_{t}\times
n_{t}}$, respectively, where $n_{s}$ and $n_{t}$ are the corresponding number
of source and target nodes.

We further assume that both graphs have identical features, i.e.,
$X_{s}\in\mathbb{R}^{n_{s}\times m}$ and $X_{t}\in\mathbb{R}^{n_{t}\times m}$.
Both graphs are also assumed to have the same set of class labels, i.e.,
$Y_{s}\in\\{0,1,\cdots,k-1\\}^{n_{s}}$ and
$Y_{t}\in\\{0,1,\cdots,k-1\\}^{n_{t}}$, where $k$ is the number of classes.
For ASD detection, we consider a binary node classification problem with
$k=2$. At times, we also consider the combined graph
$\mathcal{G}_{c}(V_{c},E_{c},X_{c},Y_{c})$, where $V_{c}=V_{s}\cup V_{t}$,
$E_{c}=E_{s}\cup E_{t}$, $X_{c}=[X_{s};X_{t}]\in\mathbb{R}^{n\times m}$,
$A_{c}=[A_{s}\ \mathbf{0};\ \mathbf{0}^{T}\ A_{t}]\in\mathbb{R}^{n\times n}$,
$Y_{c}\in\\{0,1,\cdots,k-1\\}^{n}$, $n=n_{s}+n_{t}$, and $\mathbf{0}$ is an
$n_{s}\times n_{t}$ null matrix.

######  Definition 1 (Cross Network Node Classification (CNNC))

Given a source graph, $\mathcal{G}_{s}=(V_{s},E_{s},X_{s},Y_{s})$ and target
graph, $\mathcal{G}_{t}=(V_{t},E_{t},X_{t},Y_{t})$, our goal is to learn a
target function, $f:V\rightarrow\\{0,1,\cdots,k-1\\}$, that accurately
classifies the labeled nodes in $Y_{s}$ as well as the unlabeled nodes in
$Y_{t}$.

###  III-B Graph Convolutional Network (GCN)

Graph convolutional networks [13] employ a message passing strategy to
succinctly capture both the node features and graph structure information when
learning the feature embedding of a node in the graph. Specifically, each
“message” corresponds to the feature embedding information of a node, which
will be passed to all of its immediate neighbors. By aggregating the features
gathered from the neighbors, a new embedding for the node will be generated.

The message passing strategy can be formally stated as follows. Given an
adjacency matrix $A$ and node feature matrix $X$, the feature embedding of the
nodes in layer $l+1$, $H^{(l+1)}$, can be computed based on its feature
embedding at its previous layer, $H^{(l)}$, as follows:

$H^{(l+1)}=\sigma(\tilde{D}^{-\frac{1}{2}}\tilde{A}\tilde{D}^{-\frac{1}{2}}H^{(l)}W^{(l)})$
(1) $\tilde{A}=A+I$ $\tilde{D}_{ii}=\sum_{j}\tilde{A}_{ij}$

where $H^{(0)}=X$ (i.e., original node features), $\sigma(\cdot)$ is the ReLU
activation function, and $\tilde{D}$ is a diagonal matrix containing the sum
of the edge weights associated with each node in the graph.

###  III-C Optimal Transport

Figure 2: Illustration of the cost and transport plan matrices of optimal
transport

Figure 3: High level architecture of our model. Combined source and target
datasets are fed to two GCN layers and two NFT layers to learn two distinct
embeddings. In pretraining (left) these embeddings are concatenated and sent
directly to a fully connected classifier. Model weights are then updated using
just cross entropy loss. In OT training (right) concatenated embeddings are
routed to an OT layer before being sent to the classifier. The OT layer
replaces the source embedding with a transported version and then sends the
target embedding and the transported source embedding to the classifier. At
this point the model weights are updated with a combination of cross entropy
and OT losses.

Optimal Transport (OT) provides a principled approach for comparing two
probability distributions by finding the least costly way to reshape one of
the distributions into the other while incorporating their geometric
information. The original OT problem was proposed by Monge [29], with its
modern formulation being attributed to Kantorovich [30]. Given two marginal
distributions $\mu_{s}\in\mathbb{R}^{n_{s}}$ and
$\mu_{t}\in\mathbb{R}^{n_{t}}$ let $C\in\mathbb{R}^{n_{s}\times n_{t}}$ be a
cost matrix, where $C_{ij}$ is the cost of transporting one unit from
$\mu^{s}_{i}$ to $\mu^{t}_{j}$. Consider a transportation plan matrix
$P\in\mathbb{R}^{n_{s}\times n_{t}}$, where $p_{ij}$ is the proportion of
probability mass $\mu^{s}_{i}$ that should be moved to $\mu^{t}_{j}$.

The optimal transport (OT) problem seeks to find a transportation plan $P$
that minimizes the total transportation cost. The minimum transportation cost
$W(\mu_{s},\mu_{t})$ is also called the Wasserstein distance or the earth-
mover distance.

$W(\mu_{s},\mu_{t}):=\min_{P\in U(\mu_{s},\mu_{t})}\langle P,C\rangle_{F}$ (2)
$U(\mu_{s},\mu_{t}):=\\{P\in\mathbb{R}_{+}^{n_{s}\times
n_{t}}|P1_{n_{t}}=\mu_{s},P^{T}1_{n_{s}}=\mu_{t}\\}$

The discrete OT formulation [31] shown above can be solved as a linear
programming problem. However, it is computationally expensive
$\mathcal{O}((n+m)nm\log(n+m))$, unstable, and is not guaranteed to have
unique solution. Fortunately these problems can be addressed by employing the
Sinkhorn distance $W_{\lambda}$ [32] shown below, which utilizes an entropy
regularization function $H(P)$ to accelerate the OT computation. Specifically,
$W_{\lambda}(\mu_{s},\mu_{t})$ and $P^{\lambda}$ can be solved using the well-
known Sinkhorn algorithm as established by Cuturi [32].

$\displaystyle W_{\lambda}(\mu_{s},\mu_{t})$ $\displaystyle=$
$\displaystyle\langle P^{\lambda},C\rangle$ (3) $\displaystyle P^{\lambda}$
$\displaystyle=$ $\displaystyle\underset{P\in
U(\mu_{s},\mu_{t})}{\mathrm{argmin}}\langle P,C\rangle-\lambda H(P)$
$\displaystyle U(\mu_{s},\mu_{t})$ $\displaystyle:=$
$\displaystyle\left\\{P\in\mathbb{R}_{+}^{n_{s}\times
n_{t}}\big{|}P1_{n_{t}}=\mu_{s},P^{T}1_{n_{s}}=\mu_{t}\right\\}$

where $H(P)=-\sum_{ij}P_{ij}\log P_{ij}$ is the entropy regularization and
$\lambda$ is a user-specified regularization hyperparameter.

##  IV Methodology

We employed a combination of graph convolutional networks (GCN) with node
feature transformation layers to learn the feature embedding of the nodes in
the source and target graphs. We then performed optimal transport on the
learned embedding of the source nodes to match the distribution of the learned
embedding of the target nodes. Figure 3 provides a high-level illustration of
our proposed deep neural network architecture along with its training
procedure. Details of the proposed architecture are described in the
subsections below.

###  IV-A Initial Node Embedding Construction & Pretraining

As previously noted, conventional GCN employs a message passing strategy to
transmit information about the feature embedding of a node to its immediate
neighbors. This allows each node to aggregate the feature information of its
neighbors when constructing its own latent embedding. Unfortunately, in graphs
with low homophily, it is possible that the neighborhood information obtained
via message passing is of less value than the original feature information of
the node itself, which was the case with the population graph constructed from
the ABIDE dataset. However, due to its current formulation (see Eqn. 1),
conventional GCN may not be able to attenuate the influence of the graph
structure in comparison to the influence of the original features.

To overcome this challenge, we propose a modification to the graph
convolutional network architecture to independently learn a nonlinear
embedding of the original node features. We call this the _node feature
transformation_ (NFT) layer in Fig 3. The NFT layer consists of a combination
of linear layer plus ReLU activation functions to transform the original node
features into their corresponding nonlinear embedding. The transformed
features are then concatenated with the structural embedding of the nodes
obtained from GCN for subsequent node classification. This strategy increases
the flexibility of the model to capture the relative influence of the node
features and homophily (i.e., neighborhood features) on the classification
task. As shown in the architecture diagram, the GCN layers are trained on
combined adjacency and feature matrices of the source and target graphs, while
the NFT layers are trained on the combined feature matrices. The weights of
the GCN and NFT layers are jointly updated during backpropagation.

The pretraining process of the deep neural network can be seen on the left
side of Fig 3. The purpose of pretraining is to ensure that the full network
with optimal transport (right side of Fig 3) can be seeded with a good set of
initial weights. During pretraining, we train the NFT and GCN layers to each
produce their own feature embeddings, $H^{L}$ (output of NFT layer) and
$H^{G}$ (output of GCN layer). The two embeddings are then concatenated and
sent to a fully connected network for node classification. The network is
trained to minimize the following cross-entropy loss function:

$\textrm{Pretraining:}\ \ \mathcal{L}_{CE}=-\frac{1}{n_{s}}\sum_{x_{i}\in
X_{s}}\sum^{l}_{j=0}Y_{ij}\log(\hat{Y}_{ij})$ (4)

###  IV-B Domain Adaptation via Optimal Transport GCN

We address the domain adaptation problem for cross network node classification
(CNNC) by combining the pretrained network with an optimal transport layer as
shown on the right side of Fig. 3. The optimal transport layer utilizes the
Sinkhorn algorithm [32] to learn the relevant transportation plan matrix $P$
that will transform the learned embedding of the source nodes to match the
distribution of the target node embedding.

The OT layer takes the source and target node embeddings along with the
entropy regularization hyperparameter $\lambda$ as inputs and returns a new
embedding of the source nodes, $\hat{H}_{s}$, which matches the distribution
of the target nodes, $H_{t}$. This can be accomplished by solving the OT
problem given in Eqn. (3) to learn the transportation plan $P$ and using the
barycentric mapping approach in [33] [34] to transform the source node
features. Specifically, for each source node $i$, its latent features will be
transported to a new embedding as follows:

$\hat{H}_{s,i}=\underset{H\in\mathbb{R}^{d}}{\mathrm{argmin}}\sum_{j}P^{\lambda}_{i,j}\
C(H,H^{t}_{j})$

If the cost function $C$ corresponds to the squared $l_{2}$ distance, then the
barycentric mapping reduces to the following form.

$\hat{H}_{s}=\text{diag}(P^{\lambda}1_{n_{t}})^{-1}P^{\lambda}H_{t}$

For domain adaptation, we typically choose the marginal distributions of the
source and target node embeddings, $\mu_{s}$ and $\mu_{t}$, to be a uniform
distribution. This allows us to further simplify the equation as follows:

$\hat{H}_{s}=n_{s}P^{\lambda}H_{t}$ (5)

where $n_{s}$ is the number of source nodes.

Both $\hat{H}_{s}$ and $H_{t}$ are then fed to a fully connected layer to
perform the node classification. By replacing $H_{s}$ with $\hat{H}_{s}$ we
train the model to work on data that has been transported to the target
domain. This addresses the domain adaptation problem and allows the classifier
to accurately classify the target nodes.

Finally, the full OTGCN network is trained to minimize the following joint
objective function:

$\displaystyle\mathcal{L}=\mathcal{L}_{CE}+\theta\mathcal{L}_{OT}$ (6)

where $\mathcal{L}_{CE}$ is the cross entropy loss given in Eqn. (4) and
$\mathcal{L}_{OT}$ is the optimal transport loss defined as follows:

$\mathcal{L}_{OT}=\langle P^{\lambda},C\rangle-\lambda H(P^{\lambda})$ (7)

with the hyperparameter $\theta$ controlling how much emphasis is placed on
the optimal transport term.

##  V Experimental Evaluation

This section presents the experiments performed to evaluate the effectiveness
of our proposed OTGCN framework. The source code for the framework and other
aspects of our experiments can be found at
https://github.com/ajoystephens/otgcn

###  V-A Data Preparation

TABLE I: The data collection sites for ABIDE and its selection as source/target graph. Site | Samples | Autism | Dataset
---|---|---|---
YALE | 41 | 22 | target
CALTECH | 15 | 5 | target
CMU | 11 | 6 | target
NYU | 172 | 74 | source
UM_1 | 86 | 34 | source
USM | 67 | 43 | source
UCLA_1 | 64 | 37 | source
PITT | 50 | 24 | source
MAX_MUN | 46 | 19 | source
TRINITY | 44 | 19 | source
UM_2 | 34 | 13 | source
KKI | 33 | 12 | source
LEUVEN_2 | 28 | 12 | source
LEUVEN_1 | 28 | 14 | source
OLIN | 28 | 14 | source
SDSU | 27 | 8 | source
SBL | 26 | 12 | source
STANFORD | 25 | 12 | source
OHSU | 25 | 12 | source
UCLA_2 | 21 | 11 | source
TABLE II: ABIDE Totals by Label and Dataset | Autism | Control | Total
---|---|---|---
Source | 370 (46.0%) | 434 (54.0%) | 804
Target | 33 (49.3%) | 34 (50.7%) | 67

The ABIDE [12] dataset is a combination of samples obtained from 20 different
collection sites as shown in Table I for ASD classification. We split these
collection sites into source and target sites for a total of 804 source
samples and 67 target samples. Following the terminology used in other
previous works in this area [7] [19] [4] we refer to the two classes in the
dataset as _autism_ and _control_. Autism refers to a subject who developed
ASD while control refers to a subject without ASD.

For the ABIDE dataset, the image data was prepared following the methodology
of [7]. A population graph is constructed by considering each subject as a
node and computing the similarity of their fMRI images for edge construction.
Specifically, an edge was placed between a pair of nodes only if their image
similarity exceeds certain threshold value. The edges are also weighted
according to their computed similarity. The node features were derived from
phenotype information. After removing all columns with a known direct
relationship to the class label, the features were one-hot-encoded as
applicable. They were then normalized and any features with an abnormally high
correlation to the label were also removed.

Figure 4: TSNE plot of combined embeddings before and after transport in the
first pass of our OT layer. Points represent subjects and the lines are
potential decision boundaries based off from each source dataset. The light
colors represent the source dataset, before transport on the left and after
transport on the right. The dark colors represent the target datset, which is
the same in both cases. The green circle points out a portion of the target
group which will likely be incorrectly classified before transportation
occurs. The green arrow points to that same group on the right.

###  V-B Experimental Setup

We compared the performance of OTGCN against the following baselines. Aside
from GCN, the rest of the baseline methods are designed for cross network node
classification problems.

  * •

GCN [13]: A graph neural network that uses a message-passing technique to
learn the feature embedding of the nodes. The architecture has been
incorporated into OTGCN as well as other graph neural networks [8] [9] [10].
The GCN implementation was written using code from the torch geometric library
[35] and involved two GCN layers followed by a fully connected classification
layer.

  * •

AdaGCN [26]: A CNNC technique which uses a series of GCN layers to learn
separate embeddings for the source and target networks. It then performs
domain adaptation by using an adversarial discriminator to force the two
embedding into a shared domain. The AdaGCN implementation used in this paper
was taken directly from the author’s provide source
code111https://github.com/daiquanyu/AdaGCN_TKDE. No significant changes were
made to the author’s implementation, which consists of three GCN layers and a
single layer discriminator.

  * •

ASN [28]: Another CNNC approach which consists of two separate GCN variational
autoencoders (VAE) for the source and target, a shared GCN encoder which looks
at both the source and target, and a adversarial discriminator. Here we used
source code provide by the author222https://github.com/yuntaodu/ASN with a
small modification to correct for NaNs produced by the VAE reconstruction
loss. The VAE reconstruction loss equation implemented in ASN includes a large
exponential term followed by an aggregation step which can lead to NaNs if
there is a large dataset or large values in VAE output. We addressed this
issue by restricting values in VAE-generated embedding to between -10 and 10.
Each of the three encoders in ASN were implemented with two GCN layers and the
two decoders had a single layer.

  * •

ACDNE [25]: A CNNC approach which uses special feature extraction layers to
learn separate embeddings for the features and structure of each network. It
then concatenates the two embeddings into a single source embedding and a
single target embedding before sending both to an adversarial discriminator.
We followed the authors source code333https://github.com/shenxiaocam/ACDNE
with this method as well and made no significant changes. Each of the four
feature extractors consisted of two layers.

For fair comparison we implemented OTGCN with two NFT layers and two GCN
layers, similar to most of the baseline methods. We use the macro- and
micro-F1 scores[36] as our evaluation metrics when comparing the performance
of different classification methods. We perform 10-fold cross validation to
select the hyperparameters of OTGCN as well as the reported baselines. For
OTGCN, we tune the hyperparameters for $\lambda$ and $\theta$ as well as
standard hyperparameters such as learning rate. Possible hyperparameter values
were [0.01, 0.03, 0.05] for $\lambda$ and [5, 10, 15] for $\theta$. Since the
target graph is completely unlabeled and has potentially different
distribution than the source graph, this poses a significant problem for
hyperparameter tuning. To our knowledge there is no ideal solution to this
problem. For this paper, the source dataset was split into 10 folds and the
target dataset was _excluded_ from the hyperparameter tuning process. In each
pass the fold selected for validation was removed from the remainder of the
source dataset and treated as a distinct target dataset. We followed this
process to choose the best hyperparameters for both our method and the
reported baselines.

For OTGCN the first and second NFT layers produced hidden embeddings of 32 and
64 units respectively. Similarly, the two GCN embeddings were also 32 and 64.
The selected hyperparameters were $\lambda=0.01$ and $\theta=10$.

Each of ACDNE’s feature extrators is build with two layer models which contain
64 and 32 hidden dimensions. We chose hyperparameter candidate values by
referencing the paper associated with the work. Our tuning script selected
$1\times 10^{-5}$ for the weight of the pairwise constraint loss and $1\times
10^{-5}$ for the weight of the l2 regularization term.

ADAGCN’s three GCN layers have 64, 32 and 16 hidden dimensions and the
discriminator has a single 16 unit hidden layer. We chose hyperparameter
candidate values by looking at existing values within their source code for
other datasets. There were three model-specific hyperparmeters for controlling
various portion of their loss functions. Our tuning script selected $5\times
10^{-5}$ for the weight of the l2 loss, 1.0 for the weight of the Wasserstein
loss, and 5.0 for the weight of the penalty loss.

ASN’s encoders were all two layers with 64 and 32 hidden units. Once again we
chose hyperparameter candidate values by looking at existing information
within the source code. Our tuning script chose 0.0001 for the weight of the
difference loss, 1.0 for the weight of the domain loss, and 0.5 for the weight
of the reconstruction loss.

After hyperparameter tuning each method was trained with the chosen
hyperparameters ten times with distinct random seeds. Each time the model was
trained on the entire labeled source dataset and evaluated on the target
dataset. The mean and standard deviation of the resulting micro and macro F1
results are recorded in Table III.

TABLE III: Performance results on target dataset Method | Macro F1 | Micro F1
---|---|---
GCN | 0.50265 +/- 0.03876 | 0.55821 +/- 0.02234
AdaGCN | 0.27853 +/- 0.18966 | 0.37112 +/- 0.22752
ASN | 0.38712 +/- 0.09156 | 0.52537 +/- 0.04418
ACDNE | 0.94310 +/- 0.04065 | 0.94328 +/- 0.04049
OTGCN | 0.97907 +/- 0.00733 | 0.97910 +/- 0.00731

###  V-C Results

Table III shows the micro- and macro-F1 results for OTGCN and all baselines on
the prepared ABIDE target dataset. These results demonstrate the improved
performance of OTGCN over existing CNNC baselines. The next best performing
result was from ACDNE, which did not rely on GCN for it’s graph embedding, but
rather extracted structural and feature information separately into separate
embeddings. OTGCN and ACDNE perfomed significantly better than all other
baselines, likely because these two methods do not rely strictly upon GCN and
the homophily principle.

Next we endeavor to establish the effectiveness of our optimal transport
layer. To do this we save off a GCN embedding just before we send it to the
very first OT transportation. We then transport the source embedding, use TSNE
to reduce the embeddings to two dimensions and plot the results.

Fig. 4 is an example of a plot after this process. Here the points are colored
according to label and dataset. On the left we see source before
transportation in lighter colors, autism is light blue and control is light
orange. The target dataset is plotted with it in darker colors, autism is dark
blue and control is dark orange. A dark green line illustrates a potential
decision boundary according to the source dataset, and it is clear that it
will likely misclassify a significant group of the target samples. On the
right we see a similar plot, but in this case the source samples have been
transported to the target domain via our OT layer. Once again a dark green
line illustrates a potential decision boundary based off from the source
samples. It is clear from this demonstration that a decision boundary derived
from the transported source embedding will more accurately classify subjects
in the target network.

##  VI Conclusion

This paper presents a deep learning framework called OTGCN to address the ASD
detection problem using imaging and non-imaging information from multiple
sites. OTGCN leverages ideas from graph neural networks to learn a feature
representation of the nodes and optimal transport to effectively tackle the
domain adaptation problem between source and target graphs. Our framework also
incorporates a nonlinear feature transformation layer to alleviate the issue
of graphs with low homophily. We experimentally compared the performance of
OTGCN against several state of the art CNNC baselines using the multi-site
ABIDE dataset. These experiments demonstrated the superior performance of
OTGCN over the baseline methods.

## Acknowledgment

This material is based upon work supported by the NSF Program on Fairness in
AI in collaboration with Amazon under grant IIS-1939368. Any opinion,
findings, and conclusions or recommendations expressed in this material are
those of the author(s) and do not necessarily reflect the views of the
National Science Foundation or Amazon.

## References

  * [1] C. Lord, T. S. Brugha, T. Charman, J. Cusack, G. Dumas, T. Frazier, E. J. Jones, R. M. Jones, A. Pickles, M. W. State _et al._ , “Autism spectrum disorder,” _Nature reviews Disease primers_ , vol. 6, no. 1, pp. 1–23, 2020\.
  * [2] U. Erkan and D. N. Thanh, “Autism spectrum disorder detection with machine learning methods,” _Current Psychiatry Research and Reviews Formerly: Current Psychiatry Reviews_ , vol. 15, no. 4, pp. 297–308, 2019.
  * [3] A. Zunino, P. Morerio, A. Cavallo, C. Ansuini, J. Podda, F. Battaglia, E. Veneselli, C. Becchio, and V. Murino, “Video gesture analysis for autism spectrum disorder detection,” in _2018 24th international conference on pattern recognition (ICPR)_. IEEE, 2018, pp. 3421–3426.
  * [4] X. Li, Y. Zhou, N. Dvornek, M. Zhang, S. Gao, J. Zhuang, D. Scheinost, L. H. Staib, P. Ventola, and J. S. Duncan, “Braingnn: Interpretable brain graph neural network for fmri analysis,” _Medical Image Analysis_ , vol. 74, p. 102233, 2021.
  * [5] T. Tong, K. Gray, Q. Gao, L. Chen, D. Rueckert, A. D. N. Initiative _et al._ , “Multi-modal classification of alzheimer’s disease using nonlinear graph fusion,” _Pattern recognition_ , vol. 63, pp. 171–181, 2017\.
  * [6] Y. Zeng, D. Cao, X. Wei, M. Liu, Z. Zhao, and Z. Qin, “Multi-modal relational graph for cross-modal video moment retrieval,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2021, pp. 2215–2224.
  * [7] S. Parisot, S. I. Ktena, E. Ferrante, M. Lee, R. G. Moreno, B. Glocker, and D. Rueckert, “Spectral graph convolutions for population-based disease prediction,” in _Medical Image Computing and Computer Assisted Intervention- MICCAI 2017: 20th International Conference, Quebec City, QC, Canada, September 11-13, 2017, Proceedings, Part III 20_. Springer, 2017, pp. 177–185.
  * [8] S. Parisot, S. I. Ktena, E. Ferrante, M. Lee, R. Guerrero, B. Glocker, and D. Rueckert, “Disease prediction using graph convolutional networks: application to autism spectrum disorder and alzheimer’s disease,” _Medical image analysis_ , vol. 48, pp. 117–130, 2018.
  * [9] A. Kazi, S. Shekarforoush, S. Arvind Krishna, H. Burwinkel, G. Vivar, K. Kortüm, S.-A. Ahmadi, S. Albarqouni, and N. Navab, “Inceptiongcn: receptive field aware graph convolutional network for disease prediction,” in _Information Processing in Medical Imaging: 26th International Conference, IPMI 2019, Hong Kong, China, June 2–7, 2019, Proceedings 26_. Springer, 2019, pp. 73–85.
  * [10] H. Jiang, P. Cao, M. Xu, J. Yang, and O. Zaiane, “Hi-gcn: A hierarchical graph convolution network for graph embedding learning of brain network and brain disorders prediction,” _Computers in Biology and Medicine_ , vol. 127, p. 104096, 2020.
  * [11] C. Tan, F. Sun, T. Kong, W. Zhang, C. Yang, and C. Liu, “A survey on deep transfer learning,” in _Artificial Neural Networks and Machine Learning–ICANN 2018: 27th International Conference on Artificial Neural Networks_. Springer, 2018, pp. 270–279.
  * [12] C. Craddock, S. Sikka, B. Cheung, R. Khanuja, S. S. Ghosh, C. Yan, Q. Li, D. Lurie, J. Vogelstein, R. Burns _et al._ , “Towards automated analysis of connectomes: The configurable pipeline for the analysis of connectomes (c-pac),” _Front Neuroinform_ , vol. 42, pp. 10–3389, 2013.
  * [13] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” _arXiv preprint arXiv:1609.02907_ , 2016.
  * [14] M. McPherson, L. Smith-Lovin, and J. M. Cook, “Birds of a feather: Homophily in social networks,” _Annual review of sociology_ , vol. 27, no. 1, pp. 415–444, 2001.
  * [15] J. Yuan, C. Holtz, T. Smith, and J. Luo, “Autism spectrum disorder detection from semi-structured and unstructured medical data,” _EURASIP Journal on Bioinformatics and Systems Biology_ , vol. 2017, pp. 1–9, 2016.
  * [16] R. Carette, F. Cilia, G. Dequen, J. Bosche, J.-L. Guerin, and L. Vandromme, “Automatic autism spectrum disorder detection thanks to eye-tracking and neural network-based approach,” in _Internet of Things (IoT) Technologies for HealthCare: 4th International Conference, HealthyIoT 2017, Angers, France, October 24-25, 2017, Proceedings 4_. Springer, 2018, pp. 75–81.
  * [17] H. S. Nogay and H. Adeli, “Machine learning (ml) for the diagnosis of autism spectrum disorder (asd) using brain imaging,” _Reviews in the Neurosciences_ , vol. 31, no. 8, pp. 825–841, 2020.
  * [18] B. Cheng, M. Liu, D. Shen, Z. Li, D. Zhang, and A. D. N. Initiative., “Multi-domain transfer learning for early diagnosis of alzheimer’s disease,” _Neuroinformatics_ , vol. 15, pp. 115–132, 2017.
  * [19] J. Wang, L. Zhang, Q. Wang, L. Chen, J. Shi, X. Chen, Z. Li, and D. Shen, “Multi-class asd classification based on functional connectivity and functional correlation tensor via multi-source domain adaptation and multi-view sparse representation,” _IEEE transactions on medical imaging_ , vol. 39, no. 10, pp. 3137–3147, 2020.
  * [20] H. Guan, L. Wang, and M. Liu, “Multi-source domain adaptation via optimal transport for brain dementia identification,” in _2021 IEEE 18th International Symposium on Biomedical Imaging (ISBI)_. IEEE, 2021, pp. 1514–1517.
  * [21] C. Zhang, Y. Xie, H. Bai, B. Yu, W. Li, and Y. Gao, “A survey on federated learning,” _Knowledge-Based Systems_ , vol. 216, p. 106775, 2021.
  * [22] X. Li, Y. Gu, N. Dvornek, L. H. Staib, P. Ventola, and J. S. Duncan, “Multi-site fmri analysis using privacy-preserving federated learning and domain adaptation: Abide results,” _Medical Image Analysis_ , vol. 65, p. 101765, 2020.
  * [23] Z. Wang, M. Ye, X. Zhu, L. Peng, L. Tian, and Y. Zhu, “Metateacher: Coordinating multi-model domain adaptation for medical image classification,” _Advances in Neural Information Processing Systems_ , vol. 35, pp. 20 823–20 837, 2022.
  * [24] X. Shen, Q. Dai, S. Mao, F.-l. Chung, and K.-S. Choi, “Network together: Node classification via cross-network deep network embedding,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 32, no. 5, pp. 1935–1948, 2020.
  * [25] X. Shen, Q. Dai, F.-l. Chung, W. Lu, and K.-S. Choi, “Adversarial deep network embedding for cross-network node classification,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 34, no. 03, 2020, pp. 2991–2999.
  * [26] Q. Dai, X. Shen, X.-M. Wu, and D. Wang, “Network transfer learning via adversarial domain adaptation with graph convolution,” _arXiv preprint arXiv:1909.01541_ , 2019.
  * [27] M. Wu, S. Pan, C. Zhou, X. Chang, and X. Zhu, “Unsupervised domain adaptive graph convolutional networks,” in _Proceedings of The Web Conference 2020_ , 2020, pp. 1457–1467.
  * [28] X. Zhang, Y. Du, R. Xie, and C. Wang, “Adversarial separation network for cross-network node classification,” in _Proceedings of the 30th ACM International Conference on Information & Knowledge Management_, 2021, pp. 2618–2626.
  * [29] G. Monge, “Mémoire sur la théorie des déblais et des remblais,” _Mem. Math. Phys. Acad. Royale Sci._ , pp. 666–704, 1781.
  * [30] L. V. Kantorovich, “On the translocation of masses,” _Journal of mathematical sciences_ , vol. 133, no. 4, pp. 1381–1382, 2006.
  * [31] L. Kantorovitch, “On the translocation of masses,” _Management Science_ , vol. 5, no. 1, pp. 1–4, 1958. [Online]. Available: https://doi.org/10.1287/mnsc.5.1.1
  * [32] M. Cuturi, “Sinkhorn distances: Lightspeed computation of optimal transport,” in _Advances in Neural Information Processing Systems_ , C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Weinberger, Eds., vol. 26. Curran Associates, Inc., 2013. [Online]. Available: https://proceedings.neurips.cc/paper/2013/file/af21d0c97db2e27e13572cbf59eb343d-Paper.pdf
  * [33] N. Courty, R. Flamary, and D. Tuia, “Domain adaptation with regularized optimal transport,” in _Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2014, Nancy, France, September 15-19, 2014. Proceedings, Part I 14_. Springer, 2014, pp. 274–289.
  * [34] S. Ferradans, N. Papadakis, G. Peyré, and J.-F. Aujol, “Regularized discrete optimal transport,” _SIAM Journal on Imaging Sciences_ , vol. 7, no. 3, pp. 1853–1882, 2014.
  * [35] M. Fey and J. E. Lenssen, “Fast graph representation learning with PyTorch Geometric,” in _ICLR Workshop on Representation Learning on Graphs and Manifolds_ , 2019.
  * [36] M. Grandini, E. Bagli, and G. Visani, “Metrics for multi-class classification: an overview,” _arXiv preprint arXiv:2008.05756_ , 2020.


# Multi-layer Clustering-based Residual Sparsifying Transform for Low-dose CT
Image Reconstruction

Xikai Yang , Zhishen Huang , Yong Long11footnotemark: 1 , Saiprasad
Ravishankar22footnotemark: 2 University of Michigan - Shanghai Jiao Tong
University Joint Institute, Shanghai Jiao Tong University, Shanghai 200240,
China<EMAIL_ADDRESS>of Computational Mathematics, Science
and Engineering (CMSE), Michigan State University, East Lansing, MI 48824,
USADepartment of Biomedical Engineering, Michigan State University, East
Lansing, MI 48824, USA

###### Abstract

Purpose: The recently proposed sparsifying transform models incur low
computational cost and have been applied to medical imaging. Meanwhile, deep
models with nested network structure reveal great potential for learning
features in different layers. In this study, we propose a network-structured
sparsifying transform learning approach for X-ray computed tomography (CT),
which we refer to as multi-layer clustering-based residual sparsifying
transform (MCST) learning. The proposed MCST scheme learns multiple different
unitary transforms in each layer by dividing each layer’s input into several
classes. We apply the MCST model to low-dose CT (LDCT) reconstruction by
deploying the learned MCST model into the regularizer in penalized weighted
least squares (PWLS) reconstruction.

Methods: The proposed MCST model combines a multi-layer sparse representation
structure with multiple clusters for the features in each layer that are
modeled by a rich collection of transforms. We train the MCST model in an
unsupervised manner via a block coordinate descent algorithm. Since our method
is patch-based, the training can be performed with a limited set of images.
For CT image reconstruction, we devise a novel algorithm called PWLS-MCST by
integrating the pre-learned MCST signal model with PWLS optimization.

Results: We conducted LDCT reconstruction experiments on XCAT phantom data and
Mayo Clinic data. We trained the MCST model with 2 (or 3) layers and with 5
clusters in each layer. The learned transforms in the same layer showed rich
features while additional information is extracted from representation
residuals. Our simulation results demonstrate that PWLS-MCST achieves better
image reconstruction quality than the conventional FBP method and PWLS with
edge-preserving (EP) regularizer. It also outperformed recent advanced methods
like PWLS with a learned multi-layer residual sparsifying transform prior
(MARS) and PWLS with a union of learned transforms (ULTRA), especially for
displaying clear edges and preserving subtle details.

Conclusions: In this work, a multi-layer sparse signal model with a nested
network structure is proposed. We dub this novel model as the MCST model that
exploits multi-layer residual maps to sparsify the underlying image and
clusters the inputs in each layer for accurate sparsification. We presented a
new PWLS framework with a learned MCST regularizer for LDCT reconstruction.
Experimental results show that the proposed PWLS-MCST provides clearer
reconstructions than several baseline methods.

##  1 Introduction

Natural signals such as images, audio and video often demonstrate sparsity in
some representation domain such as wavelets 17, Fourier transform, etc.
Dictionary learning methods can extract features from natural signal data and
represent the natural signal in a sparse format. Such sparse representation
methods have been used in many areas like image denoising 6, super-resolution
33, inpainting 12, medical image reconstruction 34, etc.

There are two typical approaches to find a sparse representation of a signal:
the synthesis model and the analysis model. Let $\mathbf{y}$ be the
measurement outcome of signal $\mathbf{x}$ with measurement operator $\Phi$.
For synthesis models, a signal $\mathbf{x}\in\mathbb{R}^{n}$ is encoded as a
combination of column vectors in a synthesis dictionary
$\mathbf{\Omega}_{S}\in\mathbb{R}^{n\times m}$ by assuming
$\mathbf{x}=\mathbf{\Omega}_{S}\mathbf{z}$, where
$\mathbf{z}\in\mathbb{R}^{m}$ contains only a few nonzero elements, and then
the data consistency is enforced simultaneously in the form of
$\min_{\mathbf{x}}\|\mathbf{y}-\Phi\mathbf{x}\|_{2}$. For analysis models, a
signal is assumed to be sparse in an (overcomplete) transform domain with
$\mathbf{\Omega}_{A}\mathbf{x}=\mathbf{z}$, where
$\mathbf{\Omega}_{A}\in\mathbb{R}^{m\times n}$ is referred to as the analysis
dictionary, and meanwhile one in parallel enforces the data consistency
$\min_{\mathbf{x}}\|\mathbf{y}-\Phi\mathbf{x}\|_{2}$.

The sparse coding step is NP-hard, and the algorithms to learn synthesis or
analysis dictionaries are computationally expensive such as K-SVD 1, Analysis
K-SVD 22, NAAOL 13, etc. The sparsifying transform model manifests advantage
in terms of low computational cost for computing the sparse coding of signals
21. By modeling $\mathbf{\Omega}\mathbf{x}=\mathbf{z}+\eta$ where $\eta$
indicates an error term, $\mathbf{x}$ is approximately sparse in the transform
domain. Such a sparsifying transform model 19 can be regarded as a generalized
analysis model. For instance, Ravishankar et al. proposed to learn doubly
sparse transforms for signals 18, where the sparsifying transform is a product
of two different transforms, a sparse transform and a fast, analytical
transform. Other works also developed efficient blind compressed sensing
algorithms 20 for joint sparsifying transform learning and image
reconstruction for MRI.

Recent years have witnessed the growing development of deep learning methods.
A typical method is to employ deep neural networks (DNN) and train an end-to-
end model between the input and the output. For example, Jin et al. proposed
the FBPConvNet based on the U-net framework 10, which constructs a mapping
from noisy low-dose CT reconstruction images to high-quality images. Chen et
al. combined the autoencoder and deconvolution network and proposed the RED-
CNN framework for low-dose CT imaging 3. Despite that DNN-based methods show
potential for capturing underlying features of big data, a large amount of
(often paired) training data is indispensable to train a feasible network. In
addition, an explicit understanding of the modeling or representation captured
in the network is currently lacking.

In the past few years, researchers imitated the nested network structure and
created more complex sparse signal models. For example, Tang et al. proposed a
Deep Micro-Dictionary Learning and Coding Network, where classic convolutional
layers are replaced with dictionary learning and feature encoding operations
8. Singhal et al. 29 presented a deep dictionary learning network (DDL) for
image classification by minimizing the reconstruction error after going
through multiple DL feature encoding steps. Ravishankar et al. exploited
multi-layer extension for sparsifying transform model, which enables
sparsifying an input image layer by layer 23. Other works proposed to improve
signal models by combining multiple signal modules in a parallel structure 30,
35. Wen et al. extended the single sparsifying transform model to a union of
sparsifying transforms model, where different inputs are assigned to different
clusters and transforms are learned with respect to each cluster 30.

Several applications of sparsifying transform (ST) learning and various
extensions have been demonstrated in the field of medical image reconstruction
such as magnetic resonance image reconstruction and computed tomography (CT)
image reconstruction 35, 32. The conventional method for regular-dose CT image
reconstruction is the analytical filtered back-projection (FBP) 7. However,
with reduced/low X-ray dose, severe artifacts and noise can degrade the
quality of the reconstructed image. A class of reconstruction methods takes
the physics model of the imaging system into account along with statistical
models of measurements and noise and often simple object priors. These methods
are referred to as model-based image reconstruction methods (MBIR) 9.

The penalized weighted least squares (PWLS) approach is a common MBIR method,
whose loss function includes a weighted quadratic data-fidelity term and a
regularizer term. The regularizer is conventionally based on hand-crafted
models and priors 24, 27, 28, 2. Previous studies manifested a promising
result by introducing learned sparsifying transform learning into the PWLS
scheme 15, 16, 14, 4. The sparsifying transform learned in an unsupervised
manner is incorporated into a prior term, which is equivalent to involving
effective information extracted from the training slices into the process of
image reconstruction. Moreover, using multiple sparsifying transforms
increases the number of obtained features 35, 32 and further boosts the
performance of reconstruction algorithms.

In this study, we propose a network-structured sparsifying transform learning
module, which we refer to as the Multi-layer Clustering-based residual
Sparsifying Transform (MCST) model. We extended the original sparsifying
transform model in two aspects. On the one hand, the multi-layer structure of
MCST model enables gaining potential features by exploiting residual maps
layer by layer. On the other hand, the clustering operation in each layer
implements an attention-like mechanism by grouping input patches. We propose a
pipeline for training the MCST model and show how to update the variables in
this pipeline. In addition, we have applied the proposed learned MCST model to
LDCT reconstruction through a PWLS optimization scheme. We named the new LDCT
reconstruction algorithm PWLS-MCST. Experimental results demonstrate that
PWLS-MCST not only provides better image reconstruction quality than
conventional methods like FBP and PWLS method with the non-adaptive edge-
preserving (EP) regularizer, but also outperforms several advanced methods
such as PWLS-ULTRA 35 and PWLS-MARS 32 in terms of RMSE and SSIM.

The rest of this paper is organized as follows. In Section 2, we introduce the
MCST model and describe its application to image reconstruction in the LDCT
setting. The experimental results are presented in Section 3. We summarize our
conclusions in Section 4.

##  2 Methods

###  2.1 Model description

In this section, we give a mathematical description of the MCST model and
introduce its framework. Figure 1 shows the structure of the proposed multi-
layer clustering-based residual sparsifying transform (MCST) model. Table 1 in
Appendix gives a detailed description of the notation used in our methodology.
Matrix $\mathbf{\Omega}_{l,k}$ denotes the transform in class $k$ in layer
$l$. $\mathbf{R}_{l,{C_{l,k}}}$ and $\mathbf{Z}_{l,{C_{l,k}}}$ represent
residual map and sparse code submatrices respectively, in layer $l$ with
column indices $i\in C_{l,k}$, where $C_{l,k}$ is a set containing indices of
clustered variables of class $k$ in the $l$-th layer. In particular, each
column in $\mathbf{R}_{1}$ is a vectorized overlapping image patch extracted
from the training dataset (of images) with an appropriate patch stride. Figure
2 gives an illustration of the patch extraction and patch clustering processes
on an image across multiple layers. For patches with the same index, even
though they are in different layers, the location of the included pixels in
the image is the same. Each patch is vectorized for further clustering
processing after extraction. The proposed MCST model extends the original
sparsifying transform model in two dimensions to strengthen the model’s
representation ability. The first type of extension is increasing the depth of
the model. We borrow the design idea of the recent MARS 32 model. Multi-layer
structure enforces the filtered residuals (the difference between the
transformed input signals and their sparse representation) to be further
sparsified layer by layer. Each layer keeps partial information or captures
effective features from the input. The second type of improvement is changing
the width of the MCST model. Similar to the popular attention mechanism, for
different parts of the input signal in each layer, we learn sparsifying
transform modules separately. Since our method is patch-based, a clustering
operation is necessary to assign each patch into a specific class. Moreover,
the relative order of the patches cannot be disrupted by the clustering
operation, which requires us to remember the index of each patch at each layer
and combine them in the same order before the next clustering.

Figure 1: Illustration of multi-layer clustering-based sparsifying transform
model. $\mathbf{\Omega}_{l,k}$ denotes the transform corresponding to class
$k$ in layer $l$. $\mathbf{R}_{l,{C_{l,k}}}$, $\mathbf{Z}_{l,{C_{l,k}}}$
represent residual map and sparsecode submatrices in layer $l$ with column
indices $i\in C_{l,k}$ respectively, where $C_{l,k}$ is a set containing
indices of clustered variables from class $k$ in the $l$-th layer.

The objective function corresponding to the MCST model’s training is as
follows. Due to the complexity of the proposed method, we present the notation
system used in (2.1) in Table 1.

$\displaystyle\min_{\\{\mathbf{\Omega}_{l},\mathbf{Z}_{l},C_{l,k}\\}}\sum_{l=1}^{L}\sum_{k=1}^{K_{l}}\sum_{i\in
C_{l,k}}\bigg{\\{}\|\mathbf{\Omega}_{l,k}\mathbf{r}_{l,i}-\mathbf{z}_{l,i}\|_{2}^{2}+\eta_{l}^{2}\|\mathbf{z}_{l,i}\|_{0}\bigg{\\}}$
$\displaystyle\mathrm{s.t.}\quad\begin{cases}\mathbf{r}_{l+1,i}=\mathbf{\Omega}_{l,k}\mathbf{r}_{l,i}-\mathbf{z}_{l,i},&1\leq
l\leq L-1,\\\
\mathbf{\Omega}_{l,k}\mathbf{\Omega}_{l,k}^{\top}=\mathbf{\Omega}_{l,k}^{\top}\mathbf{\Omega}_{l,k}=\mathbf{I},&\forall
l,k.\end{cases}$ (P0)

Figure 2: Illustration of the patch extraction and patch clustering processes
on an image across multiple layers. For patches with the same index, even
though they are in different layers, the location of the included pixels in
the image is the same. Each patch is vectorized for further clustering
processing after extraction.

###  2.2 Algorithm for model training and image reconstruction

Figure 3 gives an overview of the algorithms for model training and low-dose
CT reconstruction. The whole process is divided into two stages. In the
process of model training, (2.1) is solved by employing a block coordinate
descent (BCD) method. We propose an iterative algorithm involving cluster
assignment update, transform update, and sparse coding steps. Note that we use
unlabelled images to train the proposed MCST signal model. A similar cyclical
updating strategy is applied to the image reconstruction stage as well.
Instead of updating the transform in each iteration, for the reconstruction
stage, an image update step is included into the BCD algorithm. We designed
the prior function (regularizer) incorporating the learned MCST model to
harness the effective information obtained in the training stage.

Figure 3: Overview of the proposed PWLS-MCST algorithm for low-dose CT
reconstruction.

####  2.2.1 Algorithm for MCST model training

The formulation for training the MCST model is given as in (2.1). Since (2.1)
is nonconvex, similar to the recent MARS work 32, we apply the BCD algorithm
to solve this problem, which takes an iterative updating strategy among
different variables to be optimized. Algorithm 1 shows the full MCST learning
pipeline. We update variables in MCST by looping over layers. The specific
solving procedure is presented in the following steps:

  * a)

Cluster assignment update in $l$-th layer

In this step, we find the cluster assignment for every patch in the $l$th
layer by solving the following optimization problem (1). Random initialization
is employed for cluster assignment in the first iteration. When updating the
cluster assignment in the $l$-th layer, the patch assignments in the deeper
layers $\\{k(i,s),s>l\\}$ are fixed.

$\displaystyle\widehat{k}(i,l)=\arg\min_{1\leq k\leq
K_{l}}\|\mathbf{\Omega}_{l,k(i,l)}\mathbf{r}_{l,i}-\mathbf{z}_{l,i}\|_{2}^{2}+\sum_{s=l+1}^{L}\|\mathbf{\Omega}_{s,k(i,s)}\mathbf{r}_{s,i}-\mathbf{z}_{s,i}\|_{2}^{2}+\eta_{l}^{2}\|\mathbf{z}_{l,i}\|_{0}$
(1)
$\displaystyle\quad\mathrm{s.t.}\quad\mathbf{r}_{s,i}=\mathbf{\Omega}_{s-1,k(i,s-1)}\mathbf{r}_{s-1,i}-\mathbf{z}_{s-1,i},\quad
l+1\leq s\leq L-1.$

##### Interpretation

The objective function of optimization problem (1) consists of two parts:
sparse encoding residual and an $\ell_{0}$ penalty to enforce sparsity on the
encoded representation. When computing the cluster assignment of a patch $i$
in the $l$-th layer, which is denoted as $k(i,l)$, the transforms in later
layers ($\mathbf{\Omega}_{s,\cdot}$ where $s>l$) still play a role. This is
because the residuals $\mathbf{r}_{s,i}$ of the patch $i$ at later layers
($s>l$) depend on transforms ($\mathbf{\Omega}_{l,\cdot}$) and the residual
($\mathbf{r}_{l,i}$) in the $l$-th layer, which is reflected in the
constraint.

To solve problem (1), we start at layer $l$ and proceed to later layers. To
determine $\widehat{k}(i,l)$, one only needs to enumerate the values of the
objective function with $k\in[K_{l}]$, and choose the index $k$ that gives the
minimal objective function value. We repeat the solution process across all
layers until the cluster assignment no longer changes.

  * b)

Sparse coding step for $\mathbf{Z}_{l}$

The objective to update $\mathbf{Z}_{l}$ in each iteration is shown in (2).
With the cluster assignment for every patch $\\{C_{l,k},1\leq l\leq L,1\leq
k\leq K_{l}\\}$ and all transforms $\\{\mathbf{\Omega}_{l,k},1\leq l\leq
L,1\leq k\leq K_{l}\\}$ and sparse coefficient maps except the $l$th layer
$\\{\mathbf{Z}_{k},\forall k\neq l\\}$ fixed, the objective function in (2.1)
becomes

$\displaystyle\min_{\mathbf{Z}_{l}}\sum_{s=l}^{L}\sum_{k=1}^{K_{s}}\sum_{i\in
C_{s,k}}\bigg{\\{}\|\mathbf{\Omega}_{s,k}\mathbf{r}_{s,i}-\mathbf{z}_{s,i}\|_{2}^{2}\bigg{\\}}+\eta_{l}^{2}\|\mathbf{Z}_{l}\|_{0}.$
(2)

We decompose (2) into (3) which concerns only one patch, and solve for all
$\\{\mathbf{z}_{l,i},1\leq i\leq N\\}$ in parallel.

$\displaystyle\min_{\mathbf{z}_{l,i}}\sum_{s=l}^{L}\bigg{\\{}\|\mathbf{\Omega}_{s,k(i,s)}\mathbf{r}_{s,i}-\mathbf{z}_{s,i}\|_{2}^{2}\bigg{\\}}+\eta_{l}^{2}\|\mathbf{z}_{l,i}\|_{0}$
(3)

To solve (3), we fix the base layer index as $l$. Recall that the residuals
$\mathbf{r}_{s}$ in later layers $(s>l)$ depends on the residual
$\mathbf{r}_{l}$ and the encoding vector $\mathbf{z}_{l}$ in layer $l$.
Therefore, we need to peel off this dependence to find the solution for the
encoding vector in layer $l$. We use the unitary property of the transforms to
disentangle the dependence. Consider the encoding discrepancy term in a layer
$s$ where $s=l+1,l+2,\cdots,L$:

$\displaystyle\|\mathbf{\Omega}_{s,k(i,s)}\mathbf{r}_{s,i_{l}}-\mathbf{z}_{s,i}\|_{2}$
$\displaystyle=\|\underbrace{\mathbf{\Omega}_{s,k(i,s)}^{\top}\mathbf{\Omega}_{s,k(i,s)}}_{\mathbf{I}}\mathbf{r}_{s,i}-\mathbf{\Omega}_{s,k(i,s)}^{\top}\mathbf{z}_{s,i}\|_{2}$
$\displaystyle=\|\mathbf{r}_{s,i}-\mathbf{\Omega}_{s,k(i,s)}^{\top}\mathbf{z}_{s,i}\|_{2}$
$\displaystyle=\|\underbrace{\mathbf{\Omega}_{s-1,k(i,s-1)}\mathbf{r}_{s-1,i}-\mathbf{z}_{s-1,i}}_{\mathbf{r}_{s,i}}-\mathbf{b}_{i}^{(s-1)\leftarrow
s}\|_{2}$
$\displaystyle=\|\mathbf{r}_{s-1,i}-\underbrace{\mathbf{\Omega}_{s-1,k(i,s-1)}^{\top}\mathbf{z}_{s-1,i}-\mathbf{\Omega}_{s-1,k(i,s-1)}^{\top}\mathbf{b}_{i}^{(s-1)\leftarrow
s}}_{\mathbf{b}^{(s-2)\leftarrow s}_{i}}\|_{2}$
$\displaystyle=\|\mathbf{\Omega}_{s-2,k(i,s-1)}\mathbf{r}_{s-2,i}-\mathbf{z}_{s-2,i}-\mathbf{b}_{i}^{(s-2)\leftarrow
s}\|_{2}$ $\displaystyle=\cdots\cdots$
$\displaystyle=\|\mathbf{\Omega}_{l,k(i,l)}\mathbf{r}_{l,i}-\mathbf{z}_{l,i}-\mathbf{b}_{i}^{l\leftarrow
s}\|_{2},$ (4)

where we define the back-propagation information vector from the $q$-th to the
$p$-th layer for the patch $i$ as

$\displaystyle\mathbf{b}^{p\leftarrow q}_{i}$
$\displaystyle=\mathbf{\Omega}_{p+1,k(i,p+1)}^{\top}\mathbf{z}_{p+1,i}+\mathbf{\Omega}_{p+1,k(i,p+1)}^{\top}\mathbf{\Omega}_{p+2,k(i,p+1)}^{\top}\mathbf{z}_{p+2,i}+\cdots+$
(5)
$\displaystyle\quad\quad\mathbf{\Omega}_{p+1,k(i,p+1)}^{\top}\mathbf{\Omega}_{p+2,k(i,p+2)}^{\top}\cdots\mathbf{\Omega}_{q,k(i,q)}^{\top}\mathbf{z}_{q,i}$
$\displaystyle=\sum_{s=p+1}^{q}\bigg{(}\prod_{m=p+1}^{s}\mathbf{\Omega}_{m,k(i,m)}^{\top}\bigg{)}\mathbf{z}_{s,i}.$

Consequently, the objective function in the optimization problem (3) can be
re-written as

$\min_{\mathbf{z}_{l,i}}\|\mathbf{\Omega}_{l,k(i,l)}\mathbf{r}_{l,i}-\mathbf{z}_{l,i}\|^{2}_{2}+\sum_{s=l+1}^{L}\|\mathbf{\Omega}_{l,k(i,s)}\mathbf{r}_{l,i}-\mathbf{z}_{l,i}-\mathbf{b}_{i}^{l\leftarrow
s}\|^{2}_{2}+\eta_{l}^{2}\|\mathbf{z}_{l,i}\|_{0}.$

The solution of optimal $\mathbf{z}_{l,i}$ to the optimization problem above
is thus

$\widehat{\mathbf{z}}_{l,i}=\begin{cases}H_{\eta_{l}/\sqrt{L-l+1}}\bigg{(}\mathbf{\Omega}_{l,k(i,l)}\mathbf{r}_{l,i}-\frac{1}{L-l+1}\sum_{s=l+1}^{L}\mathbf{b}^{l\leftarrow
s}_{i}\bigg{)},&1\leq l\leq L-1,\\\
H_{\eta_{L}}(\mathbf{\Omega}_{L,k(i,L)}\mathbf{r}_{L,i}),&l=L,\end{cases}$ (6)

The solution form depends on the specific layer $l$, where optimization is
done in (2).

  * c)

Transform update for $\\{\mathbf{\Omega}_{l,k},1\leq k\leq K_{l}\\}$

To update the transforms in the $l$-th layer, we leave out quadratic terms
with layer depth less than $l$ and all $\|\cdot\|_{0}$ terms which are
unrelated to the target optimization variables. By assembling the column
vectors $\mathbf{r}_{s,i}$, $\mathbf{z}_{s,i}$ into $\mathbf{R}_{s,C_{s,k}}$
and $\mathbf{Z}_{s,C_{s,k}}$ respectively where $i\in C_{s,k}$, (2.1) can be
rewritten as the following optimization problem. With cluster assignment
$C_{l,k}$ and sparse code $\mathbf{Z}_{l}$ fixed, we solve (7) to update
$\mathbf{\Omega}_{l,k}$.

$\displaystyle\min_{\\{\mathbf{\Omega}_{l,k}\\}}\sum_{s=l}^{L}\sum_{k=1}^{K_{s}}\bigg{\\{}\|\mathbf{\Omega}_{s,k}\mathbf{R}_{s,C_{s,k}}-\mathbf{Z}_{s,C_{s,k}}\|_{\mathrm{F}}^{2}\bigg{\\}}\;\;\quad\mathrm{s.t.}\,\,\mathbf{\Omega}_{l,k}\mathbf{\Omega}_{l,k}^{\top}=\mathbf{\Omega}_{l,k}^{\top}\mathbf{\Omega}_{l,k}=\mathbf{I},\forall
l,k.$ (7)

Define $\mathbf{B}^{l\leftarrow s}_{C_{l,k}}$ as the matrix consisting of
back-propagation information vectors from the $s$-th layer to the $l$-th layer
for all patches in layer $l$ belonging to cluster $k$. We need to handle the
same entanglement between transforms in the former layers and later layers.
With the same argument to untwine such dependence as in equation (b)), the
objective function in problem (7) can be re-written as

$\displaystyle\min_{\\{\mathbf{\Omega}_{l,k}\\}}\sum_{k=1}^{K_{l}}\bigg{\\{}\|\mathbf{\Omega}_{l,k}\mathbf{R}_{l,C_{l,k}}-\mathbf{Z}_{l,C_{l,k}}\|_{\mathrm{F}}^{2}\bigg{\\}}+\sum_{s=l+1}^{L}\sum_{k=1}^{K_{l}}\bigg{\\{}\|\mathbf{\Omega}_{l,k}\mathbf{R}_{l,C_{l,k}}-\mathbf{Z}_{l,C_{l,k}}-\mathbf{B}^{l\leftarrow
s}_{C_{l,k}}\|_{\mathrm{F}}^{2}\bigg{\\}}$ (8)
$\displaystyle\sim\min_{\\{\mathbf{\Omega}_{l,k}\\}}\sum_{k=1}^{K_{l}}\bigg{\\{}\|\mathbf{\Omega}_{l,k}\mathbf{R}_{l,C_{l,k}}-\mathbf{Z}_{l,C_{l,k}}-\frac{1}{L-l+1}\sum_{s=l+1}^{L}\mathbf{B}^{l\leftarrow
s}_{C_{l,k}}\|_{\mathrm{F}}^{2}\bigg{\\}},$

where the symbol $\sim$ means equal up to some additive terms that do not
depend on the variables being optimized for.

From the rewritten objective function, the unitary constraint on transforms
reduces the optimization over each $\mathbf{\Omega}_{l,k}$ to the orthogonal
procrustes problem, which thus can be solved separately. Define
$\mathbf{G}_{l,k}$ as

$\mathbf{G}_{l,k}=\begin{cases}\mathbf{R}_{l,C_{l,k}}\big{(}\mathbf{Z}_{l,C_{l,k}}+\frac{1}{L-l+1}\sum_{s=l+1}^{L}\mathbf{B}_{C_{l,k}}^{l\leftarrow
s}\big{)}^{\top},&1\leq l\leq L-1,\\\
\mathbf{R}_{L,C_{L,k}}\mathbf{Z}_{L,C_{L,k}}^{\top},&l=L.\end{cases}$ (9)

In the $l$-th layer, the solution to each orthogonal procustes problem (i.e.
the solution to (7)) with respect to each transform $\mathbf{\Omega}_{l,k}$ is

$\widehat{\mathbf{\Omega}}_{l,k}=\mathbf{V}_{l,k}\mathbf{U}_{l,k}^{\top},$
(10)

where $\mathbf{U}_{l,k}$ and $\mathbf{V}_{l,k}$ are the singular vector
matrices in the singular value decomposition of the matrix
$\mathbf{G}_{l,k}=\mathbf{U}_{l,k}\mathbf{\Sigma}_{l,k}\mathbf{V}_{l,k}^{\top}$
25.

Algorithm 1 MARS Learning Algorithm 32

Input: training data $\mathbf{R}_{1}$, all-zero initial
$\\{\widetilde{\mathbf{Z}}_{l}^{(0)}\\}$, initial
$\\{\widetilde{\mathbf{\Omega}}_{1,k}^{(0)},1\leq k\leq K_{1}\\}=$ 2D DCT,
random matrices for initial $\\{\widetilde{\mathbf{\Omega}}_{l,k}^{(0)},2\leq
l\leq L,1\leq k\leq K_{l}\\}$, k-means initialization for the $1$st layer
cluster assignment $\\{\widetilde{C}_{1,k}^{(0)},1\leq k\leq K_{1}\\}$ and
random initialization for cluster assignment in the rest of layers, thresholds
$\\{\eta_{l}\\}$, number of iterations $T$.

Output: learned transforms $\\{\widetilde{\mathbf{\Omega}}_{l,k}^{(T)},1\leq
l\leq L,1\leq k\leq K_{l}\\}$.

for $t=1,2,\cdots,{T}$ do

for $l=1,2,\cdots,{L}$ do

1) Cluster assignment update for $\\{\widetilde{C}_{l,k}^{(t)},1\leq k\leq
K_{l}\\}$ by solving (1).

2) Sparse Coding for $\widetilde{\mathbf{Z}}_{l}^{(t)}$ via (6).

3) Updating $\\{\widetilde{\mathbf{\Omega}}_{l,k}^{(t)},1\leq k\leq K_{l}\\}$
via (10).

end for

end for

####  2.2.2 PWLS-MCST image reconstruction algorithm

We adopt the Penalized Weighted Least Squares (PWLS) approach to reconstruct
an image from its noisy sinogram. The learned MCST model is incorporated into
the cost function as the regularizer term. The specific PWLS problem is shown
in (P1), where $\mathbf{x}\in\mathbb{R}^{N_{p}}$ denotes the reconstructed
image and $N_{p}$ is the number of pixels in the image. Vector
$\mathbf{y}\in\mathbb{R}^{N_{d}}$ represents the noisy sinogram data and
$\mathbf{A}\in\mathbb{R}^{N_{d}\times N_{p}}$ is the system matrix of the CT
scan. $\mathbf{W}=\mathsf{diag}\\{w_{i}\\}\in\mathbb{R}^{N_{d}\times N_{d}}$
is the diagonal weighting matrix with $w_{i}$ indicating the inverse variance
of $y_{i}$. Matrix $\mathbf{H}_{\mathbf{A}}$ is the diagonal majorizing matrix
of $\mathbf{A}^{T}\mathbf{W}\mathbf{A}$. Operator $\mathbf{P}_{i}$ extracts
the $i$-th overlapping patch from the original image with patch stride of $1$
pixel. Parameter $\beta$ controls the trade-off between the data-fidelity term
and the regularizer term. $\\{\gamma_{l}\\}$ are tunable thresholds which
control the sparsity of the MCST model in each layer.

$\min_{\mathbf{x}\geq\mathbf{0}}\frac{1}{2}\|\mathbf{y}-\mathbf{A}\mathbf{x}\|^{2}_{\mathbf{W}}+\beta\mathsf{S}(\mathbf{x}),$
(P1)

where $\mathsf{S}(\mathbf{x})$ is defined as

$\displaystyle\mathsf{S}(\mathbf{x})\triangleq\min_{\\{\mathbf{Z}_{l},C_{l,k}\\}}\sum_{l=1}^{L}\sum_{k=1}^{K_{l}}\sum_{i\in
C_{l,k}}\bigg{\\{}\|\mathbf{\Omega}_{l,k}\mathbf{r}_{l,i}(\mathbf{x})-\mathbf{z}_{l,i}\|_{2}^{2}+\gamma_{l}^{2}\|\mathbf{z}_{l,i}\|_{0}\bigg{\\}}$
$\displaystyle\mathrm{s.t.}\quad\begin{cases}\mathbf{r}_{1,i}(\mathbf{x})=\mathbf{P}_{i}\mathbf{x},&\forall\,i,\\\
\mathbf{r}_{l+1,i}(\mathbf{x})=\mathbf{\Omega}_{l,k}\mathbf{r}_{l,i}(\mathbf{x})-\mathbf{z}_{l,i},&\forall
i,\,\,1\leq l\leq L-1.\end{cases}$

To solve (P1), we decompose (P1) into several subproblems (image update for
$\mathbf{x}$, cluster assignment update for $\\{C_{l,k},1\leq l\leq L,1\leq
k\leq K_{l}\\}$, sparse coding for $\\{\mathbf{Z}_{l},1\leq l\leq L\\}$) and
solve these subproblems sequentially in each iteration of the proposed
algorithm. The complete algorithm corresponding to the image reconstruction
stage is shown in Algorithm 2.

  * a)

Image update step

Here, we fix the cluster assignments $\\{C_{l,k},1\leq l\leq L,1\leq k\leq
K_{l}\\}$ and sparse representations $\\{\mathbf{Z}_{l},1\leq l\leq L\\}$. The
resulting subproblem can be rewritten as follows.

$\min_{\mathbf{x}\geq\mathbf{0}}\frac{1}{2}\|\mathbf{y}-\mathbf{A}\mathbf{x}\|^{2}_{\mathbf{W}}+\beta\mathsf{S}_{2}(\mathbf{x})$
(11)
$\displaystyle\mathsf{S}_{2}(\mathbf{x})\triangleq\sum_{l=1}^{L}\sum_{k=1}^{K_{l}}\sum_{i\in
C_{l,k}}\bigg{\\{}\|\mathbf{\Omega}_{l,k}\mathbf{r}_{l,i}(\mathbf{x})-\mathbf{z}_{l,i}\|_{2}^{2}\bigg{\\}}$
$\displaystyle\mathrm{s.t.}\quad\begin{cases}\mathbf{r}_{1,i}(\mathbf{x})=\mathbf{P}_{i}\mathbf{x},&\forall\,i\\\
\mathbf{r}_{l+1,i}(\mathbf{x})=\mathbf{\Omega}_{l,k}\mathbf{r}_{l,i}(\mathbf{x})-\mathbf{z}_{l,i},&\forall
i,\,1\leq l\leq L-1,.\end{cases}$

The solving process for (11) is identical to the image update step in Ref. 32
For each outer iteration, we update the image $T_{i}$ times with parameter
$\rho$ decreasing as in (12). The formulae to compute
$\nabla\mathsf{S}_{2}(\mathbf{x})$ and the Hessian matrix
$\mathbf{H}_{\mathsf{S}_{2}}$ for image update are given in (13) and (14),
respectively.

$\rho_{r}(\alpha)=\begin{cases}1,&r=0,\\\
\frac{\pi}{\alpha(r+1)}\sqrt{1-(\frac{\pi}{2\alpha(r+1)})^{2}},&\text{otherwise},\end{cases}\\\
$ (12)

$\nabla\mathsf{S}_{2}(\mathbf{x})=2\beta\sum_{i=1}^{N_{p}}(\mathbf{P}_{i})^{\top}\bigg{\\{}L\mathbf{P}_{i}\mathbf{x}-\sum_{k=1}^{L}\mathbf{b}^{0\leftarrow
k}_{i}\bigg{\\}},$ (13)
$\mathbf{H}_{\mathsf{S}_{2}}\triangleq\nabla^{2}\mathsf{S}_{2}(\mathbf{x})=2L\beta\sum_{i=1}^{N_{p}}(\mathbf{P}_{i})^{\top}\mathbf{P}_{i},\\\
$ (14)

  * b)

Cluster assignment update for $\\{C_{l,k},1\leq l\leq L,1\leq k\leq K_{l}\\}$

In the process of image reconstruction, the cluster assignment for every
reconstructed image patch is updated dynamically in the same manner as in the
learning stage. We find the optimal class for the $i$-th patch in the $l$-th
layer by solving the following problem as in Section 2.2.1.

$\displaystyle\widehat{k}(i,l)=\arg\min_{1\leq k\leq
K_{l}}\|\mathbf{\Omega}_{l,k}\mathbf{r}_{l,i}-\mathbf{z}_{l,i}\|_{2}^{2}+\sum_{s=l+1}^{L}\|\mathbf{\Omega}_{s,k(i,s)}\mathbf{r}_{s,i}-\mathbf{z}_{s,i}\|_{2}^{2}+\eta_{l}^{2}\|\mathbf{z}_{l,i}\|_{0}$
(15)
$\displaystyle\quad\mathrm{s.t.}\quad\begin{cases}\mathbf{r}_{1,i}&=\mathbf{P}_{i}\mathbf{x},\\\
\mathbf{r}_{s,i}&=\mathbf{\Omega}_{s-1,k(i,s-1)}\mathbf{r}_{s-1,i}-\mathbf{z}_{s-1,i},\,\,l+1\leq
s\leq L-1.\end{cases}$

  * c)

Sparse coding step for $\\{\mathbf{Z}_{l},1\leq l\leq L\\}$

In this step, by eliminating variables except the target sparsecode
$\mathbf{Z}_{l}$, (P1) gets simplified to (16) which is the same as the sparse
coding problem in (2).

$\displaystyle\min_{\mathbf{Z}_{l}}\sum_{s=l}^{L}\sum_{k=1}^{K_{s}}\sum_{i\in
C_{s,k}}\bigg{\\{}\|\mathbf{\Omega}_{s,k}\mathbf{r}_{s,i}-\mathbf{z}_{s,i}\|_{2}^{2}\bigg{\\}}+\gamma_{l}^{2}\|\mathbf{Z}_{l}\|_{0}$
(16)

The solution to (16) is given in (6).

Algorithm 2 Image Reconstruction Algorithm

Input: initial image $\widetilde{\mathbf{x}}^{(0)}$, all-zero initial
$\\{\widetilde{\mathbf{Z}}_{l}^{(0)}\\}$, pre-learned
$\\{\mathbf{\Omega}_{l}\\}$, thresholds $\\{\gamma_{l}\\}$,

$\alpha=1.999$, $\mathbf{H}_{\mathbf{A}}$, $\mathbf{H}_{\mathsf{S}_{2}}$,
number of outer iterations $T_{O}$, number of inner iterations $T_{i}$.

Output: reconstructed image $\widetilde{\mathbf{x}}^{(T_{O})}$.

for $t=0,1,2,\cdots,{T_{O}-1}$ do

1) Image Update: With $\\{\widetilde{\mathbf{Z}}_{l}^{(t)}\\}$ fixed,
Initialization: $\rho=1$, $\mathbf{x}^{(0)}=\widetilde{\mathbf{x}}^{(t)}$,
$\mathbf{g}^{(0)}=\bm{\zeta}^{(0)}=\mathbf{A}^{\top}\mathbf{W}(\mathbf{A}\mathbf{x}^{(0)}-\mathbf{y})$
and
$\mathbf{h}^{(0)}=\mathbf{H}_{\mathbf{A}}\mathbf{x}^{(0)}-\bm{\zeta}^{(0)}$.

for $r=0,1,2,\cdots,T_{i}-1,$ do $\left\\{\begin{aligned}
\mathbf{s}^{(r+1)}&=\rho(\mathbf{H}_{\mathbf{A}}\mathbf{x}^{(r)}-\mathbf{h}^{(r)})+(1-\rho)\mathbf{g}^{(r)}\\\
\mathbf{x}^{(r+1)}&=[\mathbf{x}^{(r)}-(\rho\mathbf{H}_{\mathbf{A}}+\mathbf{H}_{\mathsf{S}_{2}})^{-1}(\mathbf{s}^{(r+1)}+\nabla\mathsf{S}_{2}(\mathbf{x}^{(r)}))]_{+}\\\
\bm{\zeta}^{(r+1)}&\triangleq\mathbf{A}^{\top}\mathbf{W}(\mathbf{A}\mathbf{x}^{(r+1)}-\mathbf{y})\\\
\mathbf{g}^{(r+1)}&=\frac{\rho}{\rho+1}(\alpha\bm{\zeta}^{(r+1)}+(1-\alpha)\mathbf{g}^{(r)})+\frac{1}{\rho+1}\mathbf{g}^{(r)}\\\
\mathbf{h}^{(r+1)}&=\alpha(\mathbf{H}_{\mathbf{A}}\mathbf{x}^{(r+1)}-\bm{\zeta}^{(r+1)})+(1-\alpha)\mathbf{h}^{(r)}\end{aligned}\right.$

decreasing $\rho$ using (12).

end for

$\widetilde{\mathbf{x}}^{(t+1)}=\mathbf{x}^{(T_{i})}$.

2) Cluster assignment update: with $\widetilde{\mathbf{x}}^{(t+1)}$ and sparse
code $\\{\mathbf{Z}_{l},1\leq l\leq L\\}$ fixed, for each $1\leq l\leq L$,
update $\\{C_{l,k},1\leq k\leq K_{l}\\}$ sequentially by (15).

3) Sparse Coding: with $\widetilde{\mathbf{x}}^{(t+1)}$ and cluster assignment
$\\{C_{l,k},1\leq l\leq L,1\leq k\leq K_{l}\\}$ fixed, for each $1\leq l\leq
L$, update $\widetilde{\mathbf{Z}}_{l}^{(t+1)}$ sequentially by (6).

end for

##  3 Experiments

###  3.1 Experiment setup

In this study, we evaluate the performance of the proposed MCST sparse signal
model and the PWLS-MCST algorithm for LDCT reconstruction. We investigate the
PWLS-MCST algorithm with different numbers of layers (2, 3 layers) with each
layer containing 5 clusters to classify input patches into various groups. We
compare the proposed methods against several baseline methods including FBP 7,
PWLS with EP regularization 2, and the recent MARS 32 and ULTRA 35
regularizations in PWLS, respectively. The details of these methods are
described as follows.

  1. 1.

FBP: We used a Hanning window with $0.4$ width to accomplish the back-
projection reconstruction.

  2. 2.

PWLS-EP: An edge-preserving (EP) regularizer is employed in the PWLS
reconstruction scheme. The mathematical representation of the EP regularizer
is $\mathsf{R}(\mathbf{x})=\sum_{j=1}^{N_{p}}\sum_{k\in
N_{j}}\kappa_{j}\kappa_{k}\phi(x_{j}-x_{k})$, where $N_{p}$ counts the overall
number of pixels, $N_{j}$ denotes the neighborhood of the $j$th pixel , and
$\kappa_{j}$ and $\kappa_{k}$ denote the analytically determined parameters to
encourage uniform resolution.
$\phi(t)\triangleq\delta^{2}(\sqrt{1+|t/\delta|^{2}}-1)$ is a potential
function.

  3. 3.

PWLS-MARS: A multi-layer residual sparsifying transform (no clustering
involved) is adopted as the prior in PWLS reconstruction. For a fair
comparison, we choose MARS models with $2$, $3$ layers (MARS2, MARS3), which
have the same depth as the MCST models in our experiments.

  4. 4.

PWLS-ULTRA: PWLS with a union of transforms model, which is equivalent to the
MCST signal model with a single layer.

We used two metrics, RMSE and SSIM 31, 34, to evalute the quality of
reconstructions. We compute the root mean square error (RMSE) as RMSE
$=\sqrt{\Sigma_{i\in\text{ROI}}(\widehat{\mathbf{x}}_{i}-\mathbf{x}^{\star}_{i})^{2}/{N_{\text{ROI}}}}$,
where $\widehat{\mathbf{x}}$ and $\mathbf{x}^{\star}$ denote the reconstructed
image and ground truth image, respectively, and $N_{ROI}$ calculates the pixel
number inside the region of interest (ROI), which is a circular region
containing all structures and tissues.

Our experiments include two main parts: learning transforms for constructing
regularizers, and using them for low-dose CT image reconstruction. We work
with two datasets: XCAT phantom simulated data 26 and Mayo clinic data 11. We
generate the low-dose measurements with the “Poisson + Gaussian” noise model,
i.e.,
$\widehat{\bf{y}}_{i}=\mathrm{Poisson}\\{I_{0}\mathrm{e}^{-[\bf{Ax}]_{i}}\\}+\mathcal{N}\\{0,\sigma^{2}\\}$
5, where $I_{0}$ denotes the incident photon intensity of the X-ray beam and
we set $I_{0}=1\times 10^{4}$ per ray and with no scatter. $\sigma^{2}=5^{2}$
is the variance of electronic noise. For XCAT phantom data, we use 5 slices to
train the model. We simulate the low-dose measurements with GE 2D LightSpeed
fan-beam geometry. The size of measurements is $888\times 984$. The pixel size
for the XCAT phantom is $\Delta_{x}=\Delta_{y}=0.4883$mm. For the Mayo Clinic
data, 7 regular-dose images collected from three patients were used to train
the MCST model. We synthesized the low-dose measurements using a fan-beam CT
geometry with a monoenergetic source. The sinograms are of size $736\times
1152$. Specifically, the width of the detector column is $1.2858$mm. The
distances from the source to detector and to rotation center are $1085.6$mm
and $595$mm, respectively. The size of the reconstructed images is $512\times
512$ with $\Delta_{x}=\Delta_{y}=0.69$mm.

###  3.2 Low-dose experiments with XCAT phantom dataset

####  3.2.1 Behavior of the learned MCST models

First, we display the $64\times 64$ learned transforms from the XCAT phantom
training dataset. Figure 4 shows the image set with all training slices. The
number below each subfigure indicates its location in the volume. We extracted
overlapping $8\times 8$ patches with a patch stride $1\times 1$. The total
number of training patches is approximately $8.5\times 10^{5}$.

Figure 4: XCAT phantom slices for training MCST models. The number below each
subfigure indicates its location in the volume.

We learn the transforms in the MCST method with different model depth: ULTRA
(single-layer MCST), MCST2 (2-layer), and MCST3 (3-layer), where each layer
contains 5 clusters to enable learning rich features. To verify the
effectiveness of the multi-layer and multi-cluster extension, MARS2 (2-layer)
and MARS3 (3-layer) are included for comparison. We ran 1000 iterations to
train each model and vary the parameters as for the MARS model in our previous
work 32. Specifically, we set parameters $\eta=80$ for ULTRA,
$(\eta_{1},\eta_{2})=(80,60)$ for MCST2 and MARS2,
$(\eta_{1},\eta_{2},\eta_{3})=(90,80,60)$ for MARS3,
$(\eta_{1},\eta_{2},\eta_{3})=(80,60,40)$ for MCST3. We initialize the
transforms for MCST in the first layer with 2D DCT matrices and the rest of
the transforms are initialized with random matrices. Figure 5 shows the pre-
learned $64\times 64$ transforms for ULTRA with 5 clusters (shown in the
orange box), MARS with two layers (shown in the gray box), and MCST2 with 5
clusters in both layers (shown in the blue box). Each row of the learned
transforms is reorganized as an $8\times 8$ square matrix for display
purposes. The results show that the MCST model integrates the advantages of
ULTRA model and MARS model. In particular, the transforms in the first layer
of MCST show rich features (like ULTRA), while the transforms in the second
layer capture finer features by further sparsifying the representation of
residuals.

Figure 5: The pre-learned $64\times 64$ transforms of ULTRA with 5 clusters
(shown in the orange box), MARS with two layers (shown in the gray box), and
MCST2 with 5 clusters in both layers (shown in the blue box). Each row of the
transforms is reorganized as an $8\times 8$ square matrix for simplicity.

####  3.2.2 LDCT reconstruction results with PWLS-MCST algorithm

The reconstructed slice $20$ and slice $60$ of the XCAT phantom dataset are
illustrated in Figures 6 and 7, respectively. We compare the results of PWLS-
MCST with the conventional method FBP, PWLS-EP, PWLS-ULTRA, and PWLS-MARS. We
used the FBP reconstruction as initialization for PWLS-EP and set the
regularization parameter as $\beta=2^{16}$. We ran 1000 iterations of the
relaxed LALM algorithm for PWLS-EP and ran 1500 outer iterations for other ST-
based iterative methods to ensure their convergence. We take the same training
dataset and validate image slice (slice 48 of the XCAT phantom) as in Ref. 32.
The reconstruction parameters $\\{\beta,\gamma_{l}\\}$ are varied based on
their values in Ref. 32. Specifically, we set the regularization parameter
$\beta$ and sparsity parameters $\\{\gamma_{l}\\}$ as $(\beta,\gamma)=(2\times
10^{5},30)$ for ULTRA, $(\beta,\gamma_{1},\gamma_{2})=(9\times 10^{4},30,10)$
for MARS2 and MCST2, $(\beta,\gamma_{1},\gamma_{2},\gamma_{3})=(9\times
10^{4},25,15,10)$ for MARS3, and
$(\beta,\gamma_{1},\gamma_{2},\gamma_{3})=(9\times 10^{4},30,12,10)$ for
MCST3, respectively.

Figure 6: Comparison of reconstructions of slice 20 of the XCAT phantom with
FBP, PWLS-EP, PWLS-ULTRA, PWLS-MARS2, PWLS-MARS3, PWLS-MCST2, and PWLS-MCST3,
respectively. The display window is $[800,1200]$ HU. Figure 7: Comparison of
reconstructed slice 60 of the XCAT phantom with FBP, PWLS-EP, PWLS-ULTRA,
PWLS-MARS2, PWLS-MARS3, PWLS-MCST2, and PWLS-MCST3, respectively. The display
window is $[800,1200]$ HU.

The reconstruction results show that PWLS-MCST attains better performance than
all baseline methods. Compared with PWLS-MARS, the images reconstructed by
PWLS-MCST demonstrate greater contrast due to flexible image modeling with
multi-class transforms learned from different groups of patches. The
performance difference is exemplified in the zoomed-in portions and details
pointed by red arrows. Also, using multiple layers enables PWLS-MCST to
preserve more key features compared to PWLS-ULTRA. For example, in the
reconstructed images by PWLS-ULTRA, some subtle bone structures are missing
(zoomed-in region at the bottom-left corner of subfigures), whereas PWLS-MCST
mitigates this effect; the zoom-in regions located at the bottom-right corner
of subfigures for PWLS-MCST are closer to the ground truth than those with
PWLS-ULTRA. We point out that PWLS-MCST provides the best performance in terms
of RMSE and SSIM.

However, compared with the 2-layer MCST model, the improvement with deeper
MCST model (MCST3) is limited, which is in accordance with the performance of
PWLS-MARS. Since XCAT phantom images have relatively simple structures, the
complexity of MCST3 models may be too high to demonstrate significant
usefulness with this dataset. A larger number of tunable parameters in deeper
MCST models might also create potential sub-optimality in the training
process.

###  3.3 Low-dose experiments with Mayo Clinic dataset

####  3.3.1 MCST model training

Next, we train the MCST model on the Mayo Clinic dataset. We use 7 regular-
dose slices collected from three different patients (L096, L067, L143) for
training. Figure 8 displays all training slices used in the experiment.
Similar to the XCAT phantom experiment, we extracted $8\times 8$ overlapping
patches from these training slices. The total number of patches is $\approx
1.8\times 10^{6}$.

Figure 8: Training dataset selected from Mayo Clinic data. Seven regular-dose
slices collected from patients L096, L067, L143 are used to train the MCST
model.

We trained the ST-based signal model with similar parameters as for the XCAT
phantom experiment. We set iteration number for training as 1000 and set
$\eta=80$ for ULTRA, $(\eta_{1},\eta_{2})=(80,60)$ for MARS2 and MCST2,
$(\eta_{1},\eta_{2},\eta_{3})=(60,60,40)$ for MARS3, and
$(\eta_{1},\eta_{2},\eta_{3})=(80,60,40)$ for MCST3. The initialization of the
transforms in the MCST model are DCT matrices for the first layer and random
matrices for successive layers, respectively.

As Figure 9 shows, the learned transforms from the ULTRA learning algorithm
contain some compound features (e.g. the third transform in the first row).
Similar to the XCAT phantom experiment, MCST2 blends the advantages of MARS
with the ULTRA model and captures richer features from the training patch set.
Furthermore, considering the complexity of the Mayo Clinic dataset, MCST2 is
able to learn effective features in every transform of the first layer.

Figure 9: Learned transforms for the Mayo Clinic data. The transforms in the
orange box are obtained from the ULTRA learning algorithm, while the
transforms in the gray and blue boxes correspond to the MARS2 and MCST2 model,
respectively.

####  3.3.2 Numerical and visual results on Mayo Clinic data

After obtaining a series of transforms, we incorporate those pre-learned
transforms into the PWLS scheme and compare various reconstruction methods.
For PWLS-EP, we iterate the relaxed LALM algorithm with 1000 iterations with
the FBP reconstruction as initialization. The hyperparameter $\beta$ is set as
$2^{15.5}$. For PWLS-ULTRA, PWLS-MARS2, and PWLS-MARS3, we set the maximum
outer iterations $T_{O}$ as 1500 and inner image update times $T_{i}=2$.
Similar to the MARS work 32, we tuned the parameters on one validation image
(slice 100 of patient L506) and tested the performance of our proposed
algorithm on other slices (slice 140 of patient L333, slice 90 of patient
L109, and slice 120 of patient L067). We varied the reconstruction parameters
based on the fine-tuned parameters in the PWLS-MARS work. Specifically, we
chose $(\beta,\gamma)=(2.5\times 10^{4},30)$ for PWLS-ULTRA,
$(\beta,\gamma_{1},\gamma_{2})=(1.8\times 10^{4},30,10)$ for PWLS-MARS2 and
PWLS-MCST2, and $(\beta,\gamma_{1},\gamma_{2},\gamma_{3})=(1.8\times
10^{4},30,12,10)$ for PWLS-MARS3 and PWLS-MCST3, respectively.

Figures 10, 11, and 12 give the image reconstruction results of slice 140 of
patient L333, slice 90 of patient L109, and slice 120 of patient L067,
respectively. The proposed PWLS-MCST algorithm gives better quality of image
reconstructions compared to other recent unsupervised learning-based methods
(PWLS-ULTRA and PWLS-MARS), especially for recovering subtle details. Although
the result of PWLS-ULTRA is clearer than the initial PWLS-EP reconstruction,
it can be improved by incorporating the rich multi-layer structure from MCST.
Numerical results indicate that PWLS-MCST performs the best in terms of RMSE
and SSIM criteria.

Figure 10: Comparison of the reconstructed slice 140 of patient L333 in the
Mayo Clinic data with FBP, PWLS-EP, PWLS-ULTRA, PWLS-MARS2, PWLS-MARS3, PWLS-
MCST2, and PWLS-MCST3 schemes, respectively. The display window is
$[800,1200]$ HU. Figure 11: Comparison of the reconstructed slice 90 of
patient L109 in the Mayo Clinic data with FBP, PWLS-EP, PWLS-ULTRA, PWLS-
MARS2, PWLS-MARS3, PWLS-MCST2, and PWLS-MCST3 schemes, respectively. The
display window is $[800,1200]$ HU. Figure 12: Comparison of the reconstructed
slice 120 of patient L067 in the Mayo Clinic data with FBP, PWLS-EP, PWLS-
ULTRA, PWLS-MARS2, PWLS-MARS3, PWLS-MCST2, and PWLS-MCST3 schemes,
respectively. The display window is $[800,1200]$ HU.

##  4 Conclusion

In this work, we proposed a multi-layer clustering-based residual sparsifying
transform (MCST) network comprised of multiple sparsifying transform modules
to extract rich information from training images. The MCST model is different
from the previous single sparsifying transform model in two aspects. First,
learning multiple transforms (in each layer) by clustering the data allows
capturing rich features from patches. Second, the multi-layer structure aims
to better capture the sparsity of the signal representation across several
layers. We applied the learned MCST model to LDCT reconstruction and presented
a novel image reconstruction algorithm dubbed PWLS-MCST. Experimental results
indicate that PWLS-MCST provides better image reconstruction quality and image
features than conventional methods such as FBP and PWLS-EP. Furthermore, PWLS-
MCST shows advantages over several recent sophisticated algorithms for LDCT
reconstruction including PWLS-MARS and PWLS-ULTRA, especially for displaying
clearer edges and preserving subtle details. In the future, we plan to involve
more sophisticated architectures like skipped connections among different
layers and multi-scale transformation through down-sampling and up-sampling
within the MCST model and formulations. Additional imaging applications can be
explored to validate the general applicability of the proposed model.

## Appendix: Notations

We use lower case bold font to denote vectors and upper case bold font to
denote matrices. $\|\cdot\|_{2}$ is the 2-norm for vectors,
$\|\cdot\|_{\mathrm{F}}$ is the Frobenius norm for matrices, and
$\|\cdot\|_{0}$ counts the non-zero entries. Operator $H_{\eta}(\cdot)$ is the
hard-thresholding function which zeros out values less than $\eta$. We list
the notation used in the description of our methodology in Table 1.

Table 1: Notation for MCST problem formulation $K_{l}$ | Number of classes in layer $l$
---|---
$C_{l,k}$ | A set containing indices of patches in class $k$ in the $l$-th layer
$\mathbf{r}_{l,i},\mathbf{z}_{l,i}$ |  $i$-th column of the matrix $\mathbf{R}_{l}$ and $\mathbf{Z}_{l}$
$\mathbf{R}_{l,C},\mathbf{Z}_{l,C}$ | Submatrix of $\mathbf{R}_{l}$ and $\mathbf{Z}_{l}$ respectively, with column indices in the set $C$
$\mathbf{\Omega}_{l,k}$ | Transform of the class $k$ in $l$-th layer
$k(i,s)$ | Cluster assignment of the patch $i$ in $s$-th layer
$\mathbf{b}^{p\leftarrow q}_{i}$ | information back-propagation vector for the patch $i$ from $q$-th layer to $p$-th layer
$\mathbf{B}_{C}^{p\leftarrow q}$ |  Information back-propagation matrix from $q$-th layer to $p$-th layer for patches with indices in the set $C$

## References

  * AEB 06 M. Aharon, M. Elad, and A. Bruckstein.  K-SVD: an algorithm for designing overcomplete dictionaries for sparse representation.  IEEE Trans. Sig. Proc., 54(11):4311–4322, November 2006.
  * CF 15 J. H. Cho and J. A. Fessler.  Regularization designs for uniform spatial resolution and noise properties in statistical image reconstruction for 3D X-ray CT.  IEEE Trans. Med. Imag., 34(2):678–689, February 2015.
  * CZK+ 17 H. Chen, Y. Zhang, M. K. Kalra, F. Lin, Y. Chen, P. Liao, J. Zhou, and G. Wang.  Low-dose CT with a residual encoder-decoder convolutional neural network.  IEEE Trans. Med. Imag., 36(12):2524–2535, 2017.
  * CZLF 17 I. Y. Chun, X. Zheng, Y. Long, and J. A. Fessler.  Efficient sparse-view X-ray CT reconstruction using $\ell_{1}$ regularization with learned sparsifying transform.  In Proc. Intl. Mtg. on Fully 3D Image Recon. in Rad. and Nuc. Med, pages 115–119, June 2017.
  * DLZF 16 Q. Ding, Y. Long, X. Zhang, and J. A. Fessler.  Modeling mixed Poisson-Gaussian noise in statistical image reconstruction for X-ray CT.  In Proc. 4th Intl. Mtg. on image formation in X-ray CT, pages 399–402, 2016.
  * EA 06 M. Elad and M. Aharon.  Image denoising via sparse and redundant representations over learned dictionaries.  IEEE Trans. Im. Proc., 15(12):3736–3745, December 2006.
  * FDK 84 L. A. Feldkamp, L. C. Davis, and J. W. Kress.  Practical cone beam algorithm.  J. Opt. Soc. Am. A, 1(6):612–619, June 1984.
  * HHW+ 19 T. Hao, W. Heng, X. Wei, W. Wei, X. Dan, Y. Yan, and S. Nicu.  Deep micro-dictionary learning and coding network.  In 2019 IEEE Winter Conference on Applications of Computer Vision (WACV), pages 386–395, 2019.
  * HYMR 21 Z. Huang, S. Ye, M.T. McCann, and S. Ravishankar.  Model-based Reconstruction with Learning: From Unsupervised to Supervised and Beyond.  arXiv e-prints, March 2021.
  * JMFU 17 K. H. Jin, M. T. McCann, E. Froustey, and M. Unser.  Deep convolutional neural network for inverse problems in imaging.  IEEE Trans. Im. Proc., 26(9):4509–4522, September 2017.
  * McC 16 C. McCollough.  TU-FG-207A-04: Overview of the low dose CT grand challenge.  Med. Phys., 43(2):3759–60, 2016.
  * MES 08 J. Mairal, M. Elad, and G. Sapiro.  Sparse representation for color image restoration.  IEEE Trans. Im. Proc., 17(1):53–69, 2008.
  * MSRM 13 Y. Mehrdad, N. Sangnam, G. Rémi, and D. E. Mike.  Constrained overcomplete analysis operator learning for cosparse signal modelling.  IEEE Trans. Sig. Proc., 61(9):2341–2355, 2013.
  * 14 L. Pfister and Y. Bresler.  Adaptive sparsifying transforms for iterative tomographic reconstruction.  In Proc. 3rd Intl. Mtg. on image formation in X-ray CT, pages 107–110, 2014.
  * 15 L. Pfister and Y. Bresler.  Model-based iterative tomographic reconstruction with adaptive sparsifying transforms.  In Proc. SPIE, volume 9020, pages 90200H–1–90200H–11, 2014\.
  * 16 L. Pfister and Y. Bresler.  Tomographic reconstruction with adaptive sparsifying transforms.  In Proc. IEEE Conf. Acoust. Speech Sig. Proc., pages 6914–6918, 2014.
  * PRK 93 Y. Pati, R. Rezaiifar, and P. Krishnaprasad.  Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition.  In Asilomar Conf. on Signals, Systems and Computers, pages 40–44 vol.1, 1993.
  * 18 S. Ravishankar and Y. Bresler.  Learning doubly sparse transforms for images.  IEEE Trans. Im. Proc., 22(12):4598–4612, December 2013.
  * 19 S. Ravishankar and Y. Bresler.  Learning sparsifying transforms.  IEEE Trans. Sig. Proc., 61(5):1072–1086, March 2013.
  * 20 S. Ravishankar and Y. Bresler.  Efficient blind compressed sensing using sparsifying transforms with convergence guarantees and application to magnetic resonance imaging.  SIAM Journal on Imaging Sciences, 8(4):2519–2557, 2015.
  * 21 S. Ravishankar and Y. Bresler.  $\ell_{0}$ sparsifying transform learning with efficient optimal updates and convergence guarantees.  IEEE Trans. Sig. Proc., 63(9):2389–2404, May 2015.
  * RPE 13 R. Rubinstein, T. Peleg, and M. Elad.  Analysis K-SVD: A dictionary-learning algorithm for the analysis sparse model.  IEEE Trans. Sig. Proc., 61(3):661–677, February 2013.
  * RW 18 S. Ravishankar and B. Wohlberg.  Learning multi-layer transform models.  In Allerton Conf. on Comm., Control, and Computing, pages 160–165, 2018.
  * SB 93 K. Sauer and C. Bouman.  A local update strategy for iterative reconstruction from projections.  IEEE Trans. Sig. Proc., 41(2):534–548, February 1993.
  * Sch 66 P. H. Schönemann.  A generalized solution of the orthogonal procrustes problem.  Psychometrika, 31:1–10, 1966.
  * SMB+ 08 W. P. Segars, M. Mahesh, T. J. Beck, E. C. Frey, and B. M. W. Tsui.  Realistic CT simulation using the 4D XCAT phantom.  Med. Phys., 35(8):3800–3808, 2008.
  * TBSH 06 J-B. Thibault, C. A. Bouman, K. D. Sauer, and J. Hsieh.  A recursive filter for noise reduction in statistical iterative tomographic imaging.  In Proc. SPIE, volume 6065, pages 60650X–1–60650X–10, 2006\.
  * TSBH 07 J-B. Thibault, K. Sauer, C. Bouman, and J. Hsieh.  A three-dimensional statistical approach to improved image quality for multi-slice helical CT.  Med. Phys., 34(11):4526–4544, November 2007.
  * VA 18 S. Vanika and M. Angshul.  Majorization minimization technique for optimally solving deep dictionary learning.  Neural Processing Letters, 47:799–814, 2018.
  * WRB 15 B. Wen, S. Ravishankar, and Y. Bresler.  Structured overcomplete sparsifying transform learning with convergence guarantees and applications.  Intl. J. Comp. Vision, 114(2-3):137–167, 2015.
  * XYM+ 12 Q. Xu, H. Yu, X. Mou, L. Zhang, J. Hsieh, and G. Wang.  Low-dose x-ray CT reconstruction via dictionary learning.  IEEE Trans. Med. Imag., 31(9):1682–1697, 2012.
  * YLR 21 X. Yang, Y. Long, and S. Ravishankar.  Multilayer residual sparsifying transform (MARS) model for low-dose ct image reconstruction.  Med. Phys., 48(10):6388–6400, 2021.
  * YWHM 10 J. Yang, J. Wright, T. S. Huang, and Y. Ma.  Image super-resolution via sparse representation.  IEEE Trans. Im. Proc., 19(11):2861–2873, 2010.
  * ZMWY 17 Y. Zhang, X. Mou, G. Wang, and H. Yu.  Tensor-based dictionary learning for spectral CT reconstruction.  IEEE Trans. Med. Imag., 36(1):142–154, 2017.
  * ZRLF 18 X. Zheng, S. Ravishankar, Y. Long, and J. A. Fessler.  PWLS-ULTRA: An efficient clustering and learning-based approach for low-dose 3D CT image reconstruction.  IEEE Trans. Med. Imag., 37(6):1498–1510, June 2018.