# Density-based Isometric Mapping ††thanks:  This is the authors’ version.
B. Yousefi, M. Khansari, R. Trask, P. Tallon, C. Carino, L. Ma are with the
University of Maryland, College Park, MD 20742.
A. Afrasiyabi is with Yale University, New Haven, CT 06520.
V. Kundra and L. Ren are with the University of Maryland School of Medicine,
Baltimore, MD 21201.
K. Farahani is with the Center for Biomedical Informatics & Information
Technology, National Cancer Institute (NCI), National Institute of Health
(NIH), Bethesda, MD 20814.
M. Hershman is with the University of Pennsylvania, and Boise Radiology Group,
Philadelphia PA 19104.
The corresponding author is B. Yousefi, BSE 4117, 9636 Gudelsky Drive,
Rockville, MD 20850, Tel: 240-665-6529, e-mail<EMAIL_ADDRESS>

Bardia Yousefi, Mélina Khansari, Ryan Trask, Patrick Tallon, Carina Carino,
Arman Afrasiyabi,
Vikas Kundra, Lan Ma, Lei Ren, Keyvan Farahani, Michelle Hershman

###### Abstract

The isometric mapping method employs the shortest path algorithm to estimate
the Euclidean distance between points on High dimensional (HD) manifolds. This
may not be sufficient for weakly uniformed HD data as it could lead to
overestimating distances between far neighboring points, resulting in
inconsistencies between the intrinsic (local) and extrinsic (global) distances
during the projection. To address this issue, we modify the shortest path
algorithm by adding a novel constraint inspired by the Parzen-Rosenblatt (PR)
window, which helps to maintain the uniformity of the constructed shortest-
path graph in Isomap. Multiple imaging datasets overall of 72,236 cases,
70,000 MINST data, 1596 from multiple Chest-XRay pneumonia datasets, and three
NSCLC CT/PET datasets with a total of 640 lung cancer patients, were used to
benchmark and validate PR-Isomap. 431 imaging biomarkers were extracted from
each modality. Our results indicate that PR-Isomap projects HD attributes into
a lower-dimensional (LD) space while preserving information, which is
visualized by the MNIST dataset indicating the maintaining local and global
distances. PR-Isomap achieved the highest comparative accuracies of 80.9%
(±5.8) for pneumonia and 78.5% (±4.4), 88.4% (±1.4), and 61.4% (±11.4) for
three NSCLC datasets, with a confidence interval of 95% for outcome
prediction. Similarly, the multivariate Cox model showed higher overall
survival, measured with c-statistics and log-likelihood test, of PR-Isomap
compared to other dimensionality reduction methods. Kaplan Meier survival
curve also signifies the notable ability of PR-Isomap to distinguish between
high-risk and low-risk patients using multimodal imaging biomarkers preserving
HD imaging characteristics for precision medicine.

_K_ eywords High dimensional data $\cdot$ Manifold learning $\cdot$ Density-
based Isometric mapping $\cdot$ Parzen-Rosenblatt constraint $\cdot$ Precision
medicine.

##  1 Introduction

High dimensional data projection in lower space is always interesting since
the data becomes sparse (in an exponential manner) by the _curse of
dimension_. Hence, several different dimension reduction (DR) [1, 2, 3, 49]
have been proposed under the machine learning paradigm to visualize or build a
model with better generalization performance. Within this context, a viable
approach to streamline data analysis involves the assumption that the dataset
in question inherently resides within a lower-dimensional (LD) space. Data
characterized by lower dimensions can be conveniently visualized within this
reduced-dimensional framework, rendering them more amenable to analysis and
interpretation. DR techniques play a pivotal role in this process, as they
facilitate the transformation of data between the HD space and the LD
embedding, serving as a bridge between these two domains. Principal component
analysis (PCA) [1] \- also called Karhunen–Loève theorem – KLT is a
representation of a stochastic process as an infinite linear combination of
orthogonal functions, analogous to a Fourier series representation of a
function on a bounded interval. This is modified by singular value
decomposition (SVD) and factor analysis. In particular, PCA is a DR technique
that effectively preserves the pattern of data in a linear subspace. However,
in many instances, data exhibits a nonlinear manifold that cannot be
geometrically represented by PCA. Here, manifolds are geometrical structures
that can be resided with a few degrees of freedom in the latent space.
Manifold models can correctly symbolize the association of data points in
physical systems or be controlled by a small number of continuous projected
parameters. Any point on the manifold can be profoundly retained by LD
vectorized representation. The overarching goal of the manifold learning
methods is to acquire appropriate non-linear embedding methods in the LD space
that preserve HD data configuration to simplify data analysis and prevent
issues with HD data. Alas, most of the well-known DR methods such as
multidimensional scaling (MDS) [2], isometrically mapping (Isomap) [3],
locally linear embedding (LLE) [4], Laplacian eigenmaps [5], Hessian eigenmaps
[6], Sammon mapping [7], and graph embedding [8] are inherently predicated
outside of a stringent uniformity assumption concerning the underlying data
distribution. Therefore, these methods are vulnerable to nonuniform data
spread in the manifold. However, in real-world applications, we may encounter
non-uniform data, particularly in medical data, which may affect the
diagnostic/prognostic outcomes.

Unsupervised methods within the spectral approach are Sammon [7] and Isomap
[3], which are related to MDS [9] and attempt to nonlinearly preserve the
similarities and distances, respectively. Isomap uses geodesic distance
instead of Euclidean distance in the formulation of the kernel [3]. Isomap is
classified as a non-linear method owing to its reliance on geodesic distance.
An important point to notice is that Isomap uses geodesic distance to stay on
the manifold for the calculation of the distance between two points far on the
curved manifold [3]. This discriminates Isomap from classical MDS, where
Euclidean distance was calculated for the points, and creates problems once
points are distant from each other. With geodesic distance in the Isomap,
there will be many smaller distances selected and approximated with Euclidean
distances with the k-nearest neighbor path defined on the manifold. Standard
Isomap considers the k-nearest neighbor for finding the shortest path with
specific k neighbors using Dijkstra’s algorithm. There is no guarantee that
these neighbors are close enough to each other so an approximation of geodesic
with Euclidean distances stays on the manifold. In other words, if we select
large k, there will be points that are far enough on the manifold that
Euclidean distance would not be calculated on the manifold and imitates the
same issue that the MDS faces. To relax the discussed data uniformity
assumption and enforce the geodesic points to stay on the manifold we propose
a constraint Parzen–Rosenblatt (PR) window [13, 14] -based method which
constraint to modify local neighbors’ graph to highlight the uniformity of the
data distribution on the HD manifold for isometric mapping, PR-Isomap, which
preserves both distances and uniformity criteria while nonlinearly reducing
the dimensionality. In terms of distant points on curved manifolds, PR-Isomap
limits the number of neighbor points selected on the manifold’s surface. This
overcomes the issue of far points on the manifold that standard Isomap does
not have, which is considered the main novelty of our approach. With PR-Isomap
there is more guarantee of approximating geodesic distance between two far
points into many smaller Euclidean distances on the manifold. This modifies
and simultaneously maintains the original assumption of isometric mapping [3].

As Fig. 1 illustrates, our dimension reduction method can extend to medical
data analysis investigation, where we generate the HD manifold using
multimodal imaging from the lung cancer NSCLC Radiogenomics dataset [29, 30,
31]. We previously applied data reduction for NSCLC prognosis for liquid
biopsy data integration ( [40, 41]), NSCLC intraobserver variability [42], and
breast cancer diagnosis while reducing data hierarchically between temporally
and for features [43, 44]. Particularly, we apply the standard t-SNE, PCA,
PHATE (Potential of Heat-diffusion for Affinity-based Trajectory Embedding)
[49], standard Isomap, and our proposed PR-Isomap method to reduce the
dimensionality of the data. PR-Isomap preserves distance in LD embedding and
showed more promising results in predicting patients’ survival while it is
tested for 1,596 cases of Pneumonia and 70,000 MNIST data. In our PR-Isomap,
for preserving local distances for the formation of data, we adopt a similar
concept as in the original Isomap to keep the geodesic distance of pairwise
points on the HD manifold. This unfolds the manifolds and makes it easier to
project data points to the LD space. We adopt PR window formation and use it
as a constraint for local neighbors to emphasize the continuity of data points
on the manifold while calculating geodesic distance. Overall, this paper
presents the following contributions:

  * -

We propose a novel DR method, PR-Isomap, which can project HD data points into
LD space while adeptly conserving both local and global distances, ensuring
the preservation of the overall data structure.

  * -

We modify Isomap to use the PR window constraint in the local neighbors’ graph
to emphasize the data’s uniformity distribution while computing LD
similarities which is critical to maintaining uniformity and maintaining
approximation of geodesic distance to many smaller Euclidean distances on the
manifold.

  * -

We conduct a comparative analysis of known state-of-the-art DR methods and
reduce HD data, and multimodal radiogenomic dataset and use the latent space
representation of these data to predict the survival of patients under cancer
treatment, Pneumonia diagnosis, and visualizing the MNIST dataset. These
applications were conducted to better understand the strengths and weaknesses
of the proposed approach. Throughout all the circumstances, PR-Isomap delivers
better results than the baselines.

In the next section, we discuss related works in DR and Isomap. In the section
after that, our DR method to tackle HD imaging biomarkers and constraint
isometric mapping is presented (Section 3). The experimental and computational
results, as well as the discussion, are then presented in Sections 4 and 5,
respectively. Finally, the conclusions are presented in Section 6.

Figure 1: Our PR-Isomap method deflates a high-dimensional manifold onto a
low-dimensional representation while preserving their similarities. We
generate the HD-manifold using NSCLC Radiogenomics dataset [29, 30, 31]. We
apply the standard t-SNE [11], Isomap [3], PCA [1], and our proposed PR-Isomap
method to reduce the dimensionality of the data. PR-Isomap preserves
similarity in low-dimensional embeddings and showed more promising results in
predicting patients’ survival.

##  2 Related Work

Unsupervised DR techniques can generally be divided into two categories:

  * LD Feature extraction: The goal here is to transform the original HD data into the LD space by extracting new features that are a linear or nonlinear combination of the original features, e.g., principal component analysis (PCA) [1]. PCA effectively preserves the pattern of data in a linear subspace. However, in many instances, data exhibits a nonlinear manifold that cannot be geometrically represented by PCA due to the linearity of the approach.

  * Manifold learning: These methodologies endeavor to unveil an LD manifold that is embedded in HD data and captures the fundamental structure of the data. Examples of manifold learning techniques include MDS [2], Isomap [3], LLE [4], Hessian eigenmaps [6], Sammon mapping [7], and graph embedding [8]. Manifold learning can be broadly divided into two categories based on the type of manifold structure that is assumed to be present in the data.

  * Global manifold learning: The assumption here is that the data can be represented by a solitary, comprehensive LD manifold that encapsulates the fundamental structure of the data. Illustrations of global manifold learning methods consist of Isomap and Laplacian Eigenmaps.

  * Local manifold learning: Such techniques operate under the assumption that the data can be modeled by a set of local LD manifolds, each capturing the structure of the HD data, i.e., LLE [4] and Hessian Eigenmaps [6]. LLE aims to reconstruct each data point linearly using its neighbors in the original space, and it employs the same weights of linear reconstruction for embedding in the LD space. Meanwhile, Kernel LLE conducts the stages of finding neighbors and linear reconstruction in the feature space.

Notwithstanding their common objective of uncovering the intrinsic structure
of data in an LD space while preserving the critical characteristics of HD
data, manifold learning techniques are differentiated by their underlying
manifold assumptions, criteria, and methods of constructing the LD
representations.

Similarity criteria for DR methods: t-distribution stochastically embedding
(t-SNE) [11] use Kullback-Leibler (KL) divergence between the two joint
probability distribution and often used for visualization without emphasizing
on uniformity of data distribution. Uniform manifold approximation and
projection (UMAP) is the only method that focuses on the uniformity of data
distributions on manifolds applying Riemannian geometry and category-theoretic
approaches with fuzzy simplicial sets [20]. UMAP exhibits a lack of strong
interpretability despite significant performance.

Distance metric-based DR methods: Sammon mapping, PHATE (Potential of Heat-
diffusion for Affinity-based Trajectory Embedding) [49], and Isomap are
related to MDS [2], attempting to nonlinearly preserve the pairwise
similarities and distances, respectively. LLE [4] and Kernel LLE [18] also use
distance measures to find neighbors and linear reconstruction in the feature
space. Isomap uses geodesic distance rather than Euclidean distance in the
kernel formulation [3].

Isometric mapping and its variations: Isomap endeavors to conserve the
intrinsic structure of the data. However, its intricate computation presents a
formidable obstacle when it comes to processing large datasets. Several Isomap
variants exist to address its computational complexity for large datasets. One
such variant is C-Isomap [50], which assigns different weights to
neighborhoods based on the density of data points, providing a more accurate
representation of the underlying data structure. Similarly, C-Isomap has a
computational complexity of $O(N^{3})$, making it challenging to use for large
datasets. To address the issue, another variant of Isomap selects a small set
of $M$ landmark points ($M<<N$) to compute the geodesic distances between each
landmark and the rest of the data points, reducing the computational
complexity to $O(M^{3})$, making it more manageable for large datasets.

Angle constraints are applied to the shortest path to ensure the smoothness of
the paths, and to accurately estimate geodesic distances, even in situations
where the surface intersects itself [23], which does not focus on the
uniformity data points. Spatiotemporal Isomap (ST-Isomap) [17] analytically
modifies the local neighbors’ graph weights to highlight the similarities to
temporally related observations. Our PR-Isomap method, as proposed, can be
effectively employed in HD manifold embeddings. It enhances the emphasis on
uniformity, preserving both pairwise distances and the overall uniformity of
the data.

Figure 2: Workflow of the proposed approach. A constrained isometrically
mapping method is used to reduce feature dimensionality. Then, a subset of
radiomic signatures is used to predict the survival of the patients using the
Cox proportional hazard multivariate model. Afterwards, the parameters of the
system were frozen, and the same analyses in an independent NSCLC dataset.

##  3 Materials and Methods

In the following sections, we introduce the isometric mapping method and
discuss the innovative techniques used to extract maximum LD information for
predicting survival and classification tasks (Fig. 1).

###  3.1 Isometrically Mapping DR

Isomap is closely linked to the standard MDS, which seeks to identify
underlying dimensions that may help to describe differences among HD features.
To deal with nonlinear complicated data, there are solutions to be followed.
Here, one way is to use kernel methods to project data to the HD space where
we can linearly group data. Another approach is to use nonlinear methods to
unfold the nonlinearity in the data. Isomap uses geodesic distance to unfold
the HD manifolds to a simpler one and makes them easier to project to LD
space. Our data points can be shown as a combination of multiple vectorized
data points,
$X=\left[\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{n}\right]\in\mathbb{R}^{d\times
n}$, where $\big{\\{}{\mathbf{x}_{i}\in\mathbb{R}^{d}}\big{\\}}_{i=1}^{n}$. We
assume that there is a finite data point drawn from our X. We can build a
geometric graph $G=(V;E)$ that has our data points as vertices and connects
vertices that are close to each other. We construct a k-nearest neighbor graph
for the $\textbf{x}_{i}$ connected to $\textbf{x}_{j}$ correspondingly and for
$\textbf{x}_{j}$ to $\textbf{x}_{i}$, show them by
$\left(\mathbf{x}_{i}\right)\ $ and $\gamma\left(\mathbf{x}_{j}\right)$,
$\gamma\left(\mathbf{x}_{i}\right),\
\gamma\left(\mathbf{x}_{j}\right)\in\mathcal{M}$. $\gamma(.)$ represents a
vector point on the manifold. The length of the path is defined by summing the
edge weights along the path on the manifold. The shortest path
($\mathcal{SP}$) distance
$D_{\mathcal{SP}}(\gamma\left(\mathbf{x}_{i}\right),\
\gamma\left(\mathbf{x}_{j}\right))$ between two vertices
$\gamma\left(\mathbf{x}_{i}\right),\ \gamma\left(\mathbf{x}_{j}\right)\in V$
and defines based on the length of the $\mathcal{SP}$ connecting them.
$g\left(\cdot\right)$ is a positive continuous function for $X$ in the given
path $\gamma(\cdot)$ connecting $\mathbf{x}_{i}\ $ and $\mathbf{x}_{j}$, which
is parameterized by $g$-length of the path as

$D_{g,\gamma}(t)=\ \int_{\gamma}{g\left(\gamma\left(t\right)\right)\ \left|\
\gamma^{\prime}\left(t\right)\right|\ \ dt},$ (1)

where $\gamma^{\prime}(\cdot)$ defines a small element on the manifold, also
known as the line integral along the path $\gamma(\cdot)$ with respect to
$g(\cdot)$. $D_{g,\gamma}\ $ denotes as distance on the manifold through the
shortest path. The g-geodesic path connecting $\mathbf{x}_{i}\ $ and
$\mathbf{x}_{j}$ is the path with minimized g-length [21]. The shortest path
is used for the Isomap algorithm to find the path on the HD manifold for
connecting two points. We assume that the LD embedding of the training data
with respect to the vectorized data points of the observation can be shown as
$Y=\left[\mathbf{y}_{1},{\
\mathbf{y}}_{2},\ldots,\mathbf{y}_{n}\right]\in\mathbb{R}^{p\times n}$, where
$\big{\\{}{\mathbf{y}_{i}\in\mathbb{R}^{p}}\big{\\}}_{i=1}^{n}$, and $p\leq
d$, desirably and usually $p\ll d$. The MDS method aims to project similar
points of the data close to each other and dissimilar points as far as
possible. This is followed by the equation below:

$\displaystyle\underset{Y}{\text{min}}\|D_{g,\gamma}(X)-D_{g,\gamma}(Y)\|_{F}^{2},$
(2)

where $D_{g,\gamma}(\cdot)$ is the shortest path distance from equation (1).
$D_{g,\gamma}\left(X\right)\ $ and $D_{g,\gamma}(Y)$ are any pairwise
distances of the data in HD space and in the subspace and can be written by
similarity $X^{T}X$ and $Y^{T}Y$ in the projected space, respectively. In MDS,
the solution for this problem is $Y=\Delta^{\frac{1}{2}}{\ V}^{T}$, where
$V^{T}$ and $\Delta\ $ are the eigenvector and eigenvalue matrices of the
input data, $X$. This is equivalent to the PCA method. If we use a generalized
kernel with respect to the kernel distance matrix $D$ as:

$\displaystyle K=\ -\frac{1}{2}HDH$ (3)

with the centering matrix of $H$, $H=I-\frac{1}{n}1{1}^{T}$, where $I$ is
centering matrix and $1$ are $[1,\dots,1]^{T}$ identity matrix. In generalized
MDS embedding we use eigen decomposition of $K$. The MDS applies Euclidean
distances between data points as a similarity measurement. In HD space with
curves, i.e., non-Romanian manifolds, the Euclidean distance may not serve as
an accurate measure of the distance between the points on the manifold.
However, it’s worth noting that the intrinsic and extrinsic distances between
two points on the curved manifold may differ. In other words, Euclidean
distance may not be an accurate representation of the true distance between
points in such spaces. Isomap solved this problem [3] by using the geodesic
distance between the points on the HD manifold. The Isomap algorithm comprises
two main steps:

1\. Creating a k-nearest-neighbor graph in HD space of the points on the
manifold, and then calculating the geodesic distance between each pair of
points applying shortest-path distance.

2\. Employ MDS and make the distance matrix to locate the points in the
embedded LD space. Isomap follows the Cayton assumptions [22, 23] to retrieve
the manifold parametrization:

  * i. An isometric embedding of the manifold in $\mathbb{R}^{p}$ exists.

  * ii. The manifold in $\mathbb{R}^{d}$ is convex as well as its parameter space. In other words, there is no hole in the manifold. Also, the geodesic distance between any two points in the parameter space remains entirely in the parameter space.

  * iii. Everywhere on the manifold, points are sufficiently sampled.

  * iv. Maintain the compactness of the manifold topologically.

Increasing the distance between the paired points mitigates the similarity of
these two distances and the accuracy of such an approximation, which
ultimately leads to the original MDS problem.

###  3.2 Constrained Isomap

The Isomap algorithm creates a k-nearest neighbor graph and uses the shortest
path algorithm to estimate the geodesic distances between points. After
finding geodesic distances between points, it searches for geodesic distances
equivalent to Euclidean distances. Due to isometrically embedded manifold
conditions, some points satisfy these conditions along with translation and
rotation, e.g., MDS is a classic example of finding these points [24]. This
study focuses on the Isomap algorithm for nonlinear manifold learning, which
aims to maintain the pairwise geodesic distances between data points [3, 25,
26]. To calculate the geodesic distance and the shortest path between two data
points, we make the assumption that the data points belong to a
$d$-dimensional manifold that is embedded in an $\mathbb{R}^{d}$ ambient
space. The crucial assumption for the Isomap method is the presence of an
isometric chart, in which pair points such as $\mathbf{x}_{i},\mathbf{x}_{j}$
lie on the HD manifold.
$D_{g,\gamma}\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)\ $ denotes the
geodesic distance of the pairwise points $\gamma:\mathcal{M}\rightarrow\
\mathbb{R}^{d}$ such that:

$D_{Euc}\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)\
\propto\|\gamma\left(\mathbf{x}_{i}\right)-\gamma\left(\mathbf{x}_{j}\right)\|_{2}$
(4)

where $\mathbf{x}_{i}\ $ and $\mathbf{x}_{j}$ are near each other on the
manifold. Geodesic distance of $\gamma\left(\mathbf{x}_{i}\right)\ $ and
$\gamma\left(\mathbf{x}_{j}\right),(\gamma\left(\mathbf{x}_{i}\right),\
\gamma\left(\mathbf{x}_{j}\right))\in\mathcal{M}$ defines along the manifold
surface, $\mathcal{M}$, and is denoted by
$D_{Euc}\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)$ and can be approximate by
$D_{g,\gamma}\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)\ $.

By building a neighborhood graph where each point is connected exclusively to
its k nearest neighbors, Isomap first calculates an estimate of the geodesic
distances between every pair of data points located on the manifold
$\mathcal{M}$; the edge weights are equivalent to the corresponding pairwise
distances. The Euclidean distance serves as an estimate for the geodesic
distance for adjacent pairs of data points. This assumption allows finding an
approximation of $D_{g,\gamma}\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)$
applying many small Euclidean distances, e.g.,

$D_{g,\gamma}\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)\approx\|\gamma\left(\mathbf{x}_{i}\right)-\gamma\left(\mathbf{x}_{j}\right)\|_{2}\
\ \forall\ \
\gamma\left(\mathbf{x}_{j}\right)\in\mathcal{N}_{k}\left(\gamma\left(\mathbf{x}_{j}\right)\right),$

where $\mathcal{N}_{k}\left(\gamma\left(\mathbf{x}_{j}\right)\right)$ denotes
the collection of the point $\gamma\left(\mathbf{x}_{j}\right)\in\
\mathcal{M}$ for $k$ closest neighbors on the manifold surface, $\mathcal{M}$.
The geodesic distance is calculated for non-neighboring points as the shortest
path length along the neighborhood graph, which can be discovered using
Dijkstra’s algorithm. The generated geodesic distance matrix is then subjected
to MDS [2] to identify a group of LD points that most closely match such
distances. Forming a k-nearest neighbor graph for the data points might create
inconsistency between geodesic and Euclidean distances in opposition to the
aforementioned condition. Particularly, if the third Cayton assumption does
not or is weakly satisfied, it indirectly implies that points in the HD
manifold should be uniformly distributed. If otherwise, for large $k$ in the
k-nearest neighbor graph there will be cases that create a discrepancy in an
approximation of the geodesic and Euclidean distances. This will cause a
difference between intrinsic and extrinsic distances on the manifold. Here, we
tackle this problem using the Parzen–Rosenblatt window constraint on the
k-nearest neighbor for a weak-uniformly distributed manifold which enforces
the accuracy of distance approximation in the Isomap algorithm (Algorithm 1
and Fig. 3). In this respect, Parzen–Rosenblatt (PR) window is defined as

$p_{h}\left(\mathbf{x}\right)=\frac{1}{k}\sum_{i=1}^{k}{\frac{1}{h^{2}}\Phi\Big{(}\frac{\mathbf{x}_{i}-\mathbf{x}}{h}\Big{)}},$
(5)

where $k$ is the number of neighbors centered in the vector, $x$, $x$,
$p_{h}(\mathbf{x})$ denotes the probability density of $\mathbf{x}$, while $h$
is the diameter of the window that helps to satisfy the approximation of
geodesic and Euclidean distances, and $\Phi$ represents the window function,
This could be a rectangular function or a uniform distribution, similar to a
window of size $h$. This function is also called the density estimation
function. In addition, PR-Isomap alters the graph weights of a portion of each
point’s k-nearest neighbors with respect to the pairwise distances on the
surface of the manifold, $\mathcal{M}$, that meet the requirements for the PR
window. The set of a PR-limited points
$\gamma\left(\mathbf{x}_{i}\right),\gamma\left(\mathbf{x}_{i}\right)\in\
\mathcal{M}$ are considered for projection as neighboring elements. Suppose
that the point $\gamma\left(\mathbf{x}_{j}\right)$, where
$\gamma\left(\mathbf{x}_{i}\right)\in\
\mathcal{N}_{k,h}\left(\gamma\left(\mathbf{x}_{j}\right)\right)$ is the
closest matching point that can be used for measuring the distance
$D_{g,\gamma}\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)$, i.e.,

$\displaystyle\underset{\gamma\left(\mathbf{x}_{k}\right)\in\
\mathcal{N}_{k,h}\left(\gamma\left(\mathbf{x}_{k}\right)\right)}{\text{min}}$
$\displaystyle D_{g,\gamma}\left(\mathbf{x}_{i},\mathbf{x}_{j}\right).$ (6)

PR-based isometric mapping uses a subset of k-nearest neighbor with a distance
less than $h$ on the surface of the manifold, $\mathcal{M}$, which is
represented by uniform neighbors $(\mathcal{UN}(\mathbf{x}_{i}))$, i.e.,

$\mathcal{UN}\left(\mathbf{x}_{i}\right)=\\{\gamma\left(\mathbf{x}_{i}\right)\in\mathcal{M},\
\ \gamma\left(\mathbf{x}_{i}\right)\in\
\mathcal{N}_{k,h}\left(\gamma\left(\mathbf{x}_{i}\right)\right),D_{g,\gamma}\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)\leq\begin{aligned}
&\underset{\gamma\left(\mathbf{x}_{k}\right)\in\
\mathcal{N}_{k,h}\left(\gamma\left(\mathbf{x}_{k}\right)\right)}{\text{min}}&&D_{g,\gamma}\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)\\\
\end{aligned}\\}$

In other words, $\mathcal{UN}\left(\mathbf{x}_{i}\right)$ is a subset of the
nearest neighbor graph, which falls into the PR radius, containing points on
the manifold $\mathcal{M}$ and they are close to the targeted point out of the
entire trivial match. $\mathcal{UN}\ $ points are employed to detect data
points on the local neighborhood of every point
$\gamma\left(\mathbf{x}_{i}\right)\ $ with a higher expectation to be
analogous to $\gamma\left(\mathbf{x}_{i}\right)$. Using the MDS generalized
optimization:

$\displaystyle\underset{Y}{\text{min}}$
$\displaystyle\|D_{g,\gamma}(X)-D_{g,\gamma}(Y)\|_{F}^{2}.$ (7) subject to
$\displaystyle\gamma\left(\mathbf{x}_{k}\right)\in\
\mathcal{N}_{k,h}\left(\gamma\left(\mathbf{x}_{k}\right)\right)$

Data:
{$\mathbf{x}_{1},\mathbf{x}_{1},\dots,\mathbf{x}_{n}\\}\in\mathbb{R}^{D}$.

Result:
$\\{\mathbf{y}_{1},\mathbf{y}_{1},\dots,\mathbf{y}_{n}\\}\in\mathbb{R}^{d}$.

Estimate: MDS on $D$\- Geodesic distance pairwise point via shortest path &
neighborhood graph dissimilarity matrix $D$. $D\in\mathbb{R}^{n\times
n}(D_{ii}=0,D_{ij}\geq 0)$, $K=\ -\frac{1}{2}HDH$

$K=\ -\frac{1}{2}HDH\>$ centering matrix of $H$, $H=I-\frac{1}{n}{11}^{T}$

$Y=\Delta^{\frac{1}{2}}{\ V}^{T}$

Constrained Dijkstra’s Algorithm

Form a weighted undirected graph $\mathcal{G}_{e}=(\nu_{e},E_{e},w_{e})$ for
q-nearest neighbor, $q\leq k$.

For $\forall$ vertex $v$, $v\in\mathcal{G}:D_{\gamma}[v]\leftarrow\infty,$
$parent[v]\leftarrow[],D_{\gamma}[s]]\leftarrow 0$,

$Q:=\forall\>$ nodes,$n$, $n\in\mathcal{G}$

while  _$Q\neq[]$_ do

$u\leftarrow v\in Q\>\&\>\min D_{\gamma}[u]$

remove $u$ from $Q$

for  _$\forall\ >\mathcal{N}_{q}[u,v]\in Q$_ do

if  _$D_{\gamma}[u,v] <h,\>\>u,v\in\mathcal{N}_{k,h}$_ then

alt $\leftarrow D_{\gamma}[u]+\mathcal{G}_{e}edge(u,v)$

if  _alt $<$ $D_{\gamma}[v]$_ then

$D_{\gamma}[v]\leftarrow$ alt

$parent[v]\leftarrow u$

end if

end if

end for

end while

Algorithm 1 PR-Isometric Mapping.

_Definition 1:_ Consider a geometric graph based on a fixed set of points
$\mathbf{x}_{1},\ldots,\mathbf{x}_{2}\in\mathbb{R}^{d}$. Let $h$ be a real
number defined for
$\mathcal{N}_{k,h}\left(\gamma\left(\mathbf{x}_{k}\right)\right)$ such that
$D_{g,\gamma}\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)\ \leq h$ implies that
$\mathbf{x}_{i}$ is connected to $\mathbf{x}_{j}$ on the graph.

Using Definition 1, if we rewrite equation (7) instead of uniform neighboring
points on the HD and the LD space, but with a difference in pairwise distances
in the HD space as a constraint:

$\displaystyle\underset{Y}{\text{min}}$
$\displaystyle\|D_{g,\gamma}(X)-D_{g,\gamma}(Y)\|_{F}^{2}$ (8) subject to
$\displaystyle D_{g,\gamma}\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)\leq h$

where $D_{g,\gamma}(X)\ $ redefines as any pairwise distances on the surface
of the HD manifold, $\mathcal{M}$, while the $\mathcal{SP}$ distance using an
unweighted PR-based k-nearest neighbor graph, $p_{h}$. Fig. 3 presents the
geometrical interpretation of the shortest path and PR-Isomap under
$\mathcal{N}_{k,h}\left(\gamma\left(\mathbf{x}_{k}\right)\right)$, where PR-
window constraints the k-nearest neighbor path.
To evaluate the strength of LD features, we use multiple classification
scenarios i.e., Random forest, and the Cox proportional hazard model to
predict patients’ Pneumonia and survival, respectively.

###  3.3 Outcome prediction and statistical analysis

In the proposed model, statistical analysis of image intensity, texture,
contextual information corresponding to intratumor, and ROIs are obtained and
listed as HD imaging attributes. These HD radiomics spanned onto LD space
while their maximum information is preserved. To evaluate the performance of
LD radiomics in predicting overall survival (OS), we performed binary
classifications to predict treatment failure. Similar to the PR-Isomap method,
we applied other commonly used state-of-the-art unsupervised DR methods to the
same data as a comparison analysis. In addition, we tested these methods on
multiple datasets for lung cancer with diverse imaging parameters (Tables
1,2,3) and Pneumonia (Fig. 4) to ensure the performance of the approach. We
categorize the classification tasks into the following divisions:
1) We used various frequently used DR models to reduce HD radiomics as
baselines such as the standard Isomap, tSNE, PCA, and PHATE models.
2) We conducted binary classification to investigate the power of prediction
using LD radiomics of multiple machine learning methods, i.e., random forest,
logistic regression, Naïve Bayesian classifier, and kernel support vector
machine (SVM). We performed 10-fold cross-validation and fixed hyperparameters
to ensure the quality of comparison.
3) The Cox proportional hazard (CPH) model is used to predict OS. Kaplan Meier
survival curves were also used to distinguish between high- and low- risk
patients with respect to median hazard.

Figure 3: Geometric interpretation of the shortest path and PR-Isomap under
$\mathcal{N}_{k,h}\left(\gamma\left(\mathbf{x}_{k}\right)\right)$,
Parzen–Rosenblatt window constraint on $k$-nearest neighbor windows on HD
manifold (a) and in the LD space after PR-Isomap projection (b).
$d_{Euc}(\mathbf{x}_{a},\mathbf{x}_{b})$ and
$D_{g,\gamma}\left(\mathbf{x}_{a},\mathbf{x}_{b}\right)\ $are Euclidian and
geodesic distances between $\mathbf{x}_{a}$ and $\mathbf{x}_{b}$.

###  3.4 Computational Complexity

PR-Isomap changes the complexity of Dijkstra’s algorithm on the new graph (see
Algorithm 1), $G_{e}$, is $O(log(N_{e})E_{e})$, $N_{e}$ and $E_{e}$ are the
number of nodes and edges in $G_{e}$, respectively. For a $q$-nearest-neighbor
graph with the parameter $q$, we have $qN/2$ edges, resulting in $N_{e}=qN/2$
nodes in $G_{e}$. Connecting each node to all $q$ neighbors, edges in $G_{e}$
in the worst case scenario is $E_{e}=((qN/2)q)/2$, leading to a complexity of
$O(log(qN/2)((qN/2)q)/2)=O((q^{2}N/4)log(qN/2))$.

##  4 Results

The proposed method modifies Isomap and examines it by HD throughputs from
various datasets.

###  4.1 Study population

Diverse datasets comprising CT, PET, and ChestXRay images derived from
patients diagnosed with Non-Small Cell Lung Cancer (NSCLC) and Pneumonia were
examined to evaluate the proposed approach. There were 70,000 digit image data
from the MNIST dataset [48] also used for visualizing the separability of the
LD features.

NSCLC Radiogenomics (Multimodal PET/CT): The first, the NSCLC-Radiogenomics
dataset [51, 30, 31], is a publicly available dataset from the National
Institutes of Health (NIH) — The Cancer Imaging Archive (TCIA), providing
images of 211 patients; we analyzed a sub-cohort of 130 patients for whom
information on miRNA was available. These provided a subset data of 57836 CT
images which had corresponding PET images. Overall, the dataset contains CT
and PET images. We used our annotations of the tumors as observed on the
medical images using maps of tumors in the CT scans and from the PET scans.
This is combined with other clinical information and survival outcomes to
predict the survival of the patients. A board-certified thoracic radiologist
(M.H.) manually segmented the tumor area using the ITK-SNAP software (version
3.6.0).

The NSCLC Radiogenomics dataset included 20, 79, 27, and 7 patients with slice
thickness (ST) of $ST\leq 1,1<ST\leq 2,2<ST\leq 4,$ and $ST\leq 5mm$,
respectively; 16 sets of images were acquired with Siemens Healthineers
(Erlangen, Germany) scanners, 110 with General Electric Healthcare (Chicago,
IL) scanners, 3 with Philips Medical Systems Technologies Ltd. (Amsterdam,
Netherlands), and 1 with Toshiba Corporation (Minato, Tokyo, Japan). This
dataset had 39 contrast-enhanced and 54 non-contrast-enhanced images (Table
1).

Figure 4: Tumors cluster in the two phenotypes. Visualizations of the original
CT/PET images with tumors in the field of view for each dataset with ChestXRay
pneumonia data. Tumors’ morphology and shape along with other features affect
the quantitative attributes obtained from the selected contours.

NSCLC-Radiomics (LUNG01): The second dataset contains images from 422 non-
small cell lung cancer (NSCLC) patients treated at MAASTRO Clinic, The
Netherlands. The patients’ pretreatment was conducted using CT scans, and
manual delineation by a radiation oncologist of the 3D volume of the total
tumor volume with additive clinical outcome data. This dataset refers to the
Lung1 dataset and is used to extract high dimensional radiomics and their
effects on prediction survival [52, 53]. The LUNG1 dataset comprised 422
patients with $2<ST\leq 4$; 324 sets of images were acquired with Siemens
Healthineers (Erlangen, Germany) scanners, and 98 with CMS Inc. Healthcare
Solutions (Charleston, SC). This dataset had 157 contrast-enhanced and 264
non-contrast-enhanced images (Tables 2.

NSCLC-Radiomics (LUNG03): The third dataset contains images and mRNA sequence
data used to investigate the association of radiomic imaging features with
gene-expression profiles consisting of 88 NSCLC CT scans [33]. The dataset
consisted of 60 ($68.1\%$) males and 28 ($31.8\%$) females; 43 ($48.9\%$)
patients with adenocarcinoma, 33 ($37.5\%$) with squamous cell carcinoma, and
12($13.4\%$) with another type of NSCLC; 39 ($44.3\%$) patients had stage I
disease, 27 ($30.7\%$) had stage II disease, 12 ($13.6\%$) had stage III, and
10 ($11.4\%$) patients had an unknown stage. There were some variabilities in
the outcome of the approach while stratifying patients based on CT parameters.
The LUNG3 dataset included 3, 27, and 58 patients with ST of $ST\leq
2,2<ST\leq 4$, and $ST\leq 5mm$, respectively; 38 were acquired with Siemens
Healthineers (Erlangen, Germany) scanners, 1 with a General Electric
Healthcare (Chicago, IL) scanner, and 49 with Philips Medical Systems
Technologies Ltd. (Amsterdam, Netherlands). This dataset had 38 contrast-
enhanced and 4 non-contrast-enhanced images, where there was no information in
the DICOM header files for the rest of these images regarding contrast
information (Tables 3).

Pneumonia dataset: A combination of 1125 chest X-ray cases collected, which
contains 500 pneumonia cases, 500 no-findings, 125 cases of COVID-19 (Fig. 4).
COVID-19 cases were selected from [34] and 1000 images selected from the
ChestX-ray8 database [35] contained 500 normal and 500 pneumonia cases [36].
55 COVID-19 cases frontal chest X-ray images were selected from Figure 1,
DarwinAI and the University of Waterloo have launched an open-source project
as a part of the COVIDx dataset to develop our models for COVID-19 detection
(COVID-Net) and COVID-19 risk stratification (COVID-RiskNet) [37, 35, 38, 39].
All the aforementioned data were collected from multiple sources, which were
publicly available databases and online websites such as Figure1.com [39],
Radiopaedia.org, the Italian Society of Medical and Interventional Radiology,
and the Hannover Medical School, even the images from online publications,
websites, or directly from the PDF using the tool pdfimages [46]. All the
images were used for different parts of the proposed approach with the overall
size of frontal X-rays of $512\times 512$ pixels [47].

MNIST: The Modified National Institute of Standards and Technology (MNIST)
dataset [48] is used for the non-biological visualization benchmark of the
proposed DR approach. It gives an intuitive understanding of the model’s
functionality. MNIST contains 60,000 small square 28×28-pixel grayscale images
of handwritten single digits in different scripts between 0 and 9. The pursued
task is to inclusively categorize a given handwritten digit while using
different DR models into 10 classes indicating 0 to 9 integer values.

Table 1: The CT parameters for NSCLC Radiogenomics dataset.  CT imaging
parameters - NSCLC Radiogenomics dataset
---
Characteristic | SIEMENS H | GE HC | PHILIPS MST | TOSHIBA C
Contrast Enhanced (CE) | 5 | 34 | - | -
Non- CE | 4 | 50 | - | 1
Unknown | 7 | 26 | 3 |
Imaging | 16 | 110 | 3 | 1
Convolution Kernel |  |  |  |
Standard (soft-tissue) | 2 | 36 | 1 (B,YB) | -
Hard (lung, bone, boneplus) | 14 | 74 | 1 (C) | 1 (FC52)
Tube Voltage (kVp) | 120 (100,140) | 120 (100,140) | 120 (120,140) | 120
Slice Thickness (mm) |  |  |  |
<1 | 12 | 7 | 1 | -
1<. <2 | 1 | 77 | - | 1
2 <. <4 | 3 | 23 | 1 | -
<5 | - | 5 | 1 | -
Table 2: The CT parameters for NSCLC Radiomics-genomics (LUNG01) dataset. CT
imaging parameters - NSCLC Radiomics-genomics LUNG01 dataset
---
Characteristic | SIEMENS H. | CMS Inc.
|  | Healthcare Sol.
Contrast Enhanced (CE) | 157 | -
Non- CE | 167 | 97
Unknown | - | 1
Imaging | 324 | 98
Convolution Kernel |  |
Standard (soft-tissue) | 306 | -
Hard (lung, bone, boneplus) | 18 | -
Unknown | - | 98
Tube Voltage (kVp) | 120 (120,140) | 119(40,360)
Slice Thickness (mm) |  |
$\leq$ 1 | - | -
1<. $\leq$ 2 | - | -
2<. $\leq$ 4 | 324 | 98
$\leq$ 5 | - | -
Table 3: The CT parameters for LUNG03 dataset. CT imaging parameters - NSCLC
Radiogenomics LUNG03 dataset
---
Characteristic | SIEMENS H. | GE HC | PHILIPS MST
Contrast Enhanced (CE) | - | - | 4
Non- CE | 37 | 1 | -
Unknown | 1 | - | 45
Imaging | 38 | 1 | 49
Convolution Kernel |  |  |
Standard (soft-tissue) | 3 | - | 41 (A,B)
Hard (lung, bone, boneplus) | 35 | 1 | 8 (C,D)
Tube Voltage (kVp) | 120 (120,140) | 120 | 120 (120,140)
Slice Thickness (mm) |  |  |
$\leq$ 1 | - | - | -
1<. $\leq$ 2 | - | - | 3
2 <. $\leq$ 4 | 7 | 1 | 19
$\leq$ 5 | 31 | - | 27
Figure 5: Visualization of the data dimensionality tested on CT radiomics of
NSCLC Radiogenomics dataset. a) Representation of HD radiomics,
$R^{(}130×431)$, presented using correlative distance. b) LD demonstration of
the HD radiomics projected in the LD space using PR-Isomap. c) Representing
the application of HD radiomics in a lower-dimensional space using the
original Isomap method. d) Using PCA, projecting HD features onto LD space. e)
t-SNE heatmap to make the same LD imaging biomarkers.

###  4.2 Implementation Details

The proposed method modifies Isomap and is examined by HD throughputs from
handwritten digits and lung cancer datasets.
Baselines. We compare PR-constrained Isomap with four baselines,

  * •

The standard t-SNE [11]: This method is the standard t-SNE that uses
Kullback–Leibler divergence similarity and measures the statistical difference
of probability distribution of points in the HD space corresponding to the
probability distribution of the same points on the LD space after projection.

  * •

PCA [1]: PCA is a popularly employed method to reduce the dimensionality of
data. It captures the maximum variance direction in the data, which maintains
the interpretability and maximum amount of information. PCA is a linear method
to reduce dimensionality.

  * •

The standard Isomap [3]: Isomap is applied for computing a quasi-isometric and
is considered to be a non-linear dimensionality reduction approach, which
attempts to maintain the geodesic distances from HD space into the LD space
embedding. The standard Isomap does not apply any constraints to maintain the
distance between the points.

  * •

PHATE [49]: PHATE emerges as a recently developed DR tool tailored for the
visualization of HD data. This method leverages a conceptual framework using
the local affinity matrix, Markov transition matrix, and Von Neumann Entropy
prior to MDS to acquire and represent the manifold, preserving not only local
but also global inter-data point distances. Operating as a visualization
technique, PHATE adeptly encapsulates the intricate local and global nonlinear
structures within the dataset, achieving this through the application of an
information-geometric distance metric interlinking data points.

To visualize the LD embedded data in 3D, we use UMAP [20] to preserve the
data’s global structure.

|  |  |  |
---|---|---|---|---
(a) PR-Isomap | (b) Isomap | (c) tSNE | (d) PCA | (e) PHATE
Figure 6: Visualizing the MNIST HD data projected to LD space in 3D space
using various DR methods. Mapping with PR-Isomap (a) showed the separation of
different categories being closer despite having lesser between-classes
(extrinsic or global) distances due to the additive constraint to maintain the
uniformity of the classes while intrinsic distances are represented well. The
results of PHATE exhibit similar behavior (e) but with lesser strength in
grouping the within-class (intrinsic) categories. Isomap (b) and tSNE (c)
showed lower uniformity and within-class grouping behavior, respectively.
Table 4: Predictive power via random forest for five different datasets,
overall of 2,236 cases across two different diseases.

Datasets | Number of cases | Modality | PR-Isomap | Isomap | PCA | tSNE | PHATE
---|---|---|---|---|---|---|---
|  |  | (Ours) |  |  |  |
| 130 | CT | 78.5 (±4.4) | 76.9(±5.8) | 73.8 (±7.6) | 76.9(±4.1) | 76.2(±7.3)
NSCLC Radiogenomics |  |  |  |  |  |  |
|  | PET | 78.5 (±5.1) | 78.5(±6.3) | 73.2 (±5.6) | 76.9(±5.3) | 75.3(±9.9)
LUNG03 | 88 | CT | 61.4 (±11.4) | 52.7(±14.8) | 49.6 (±15.1) | 40.3(±15.5) | 45.5(±14.4)
LUNG01 | 422 | CT | 88.4 (±1.4) | 87.7(±2.8) | 86.9 (±2.5) | 87.7(±3.1) | 88.3 (±0.6)
Pneumonia | 1596 | Chest X-Ray | 80.9 (±5.8) | 80.8(±4.8) | 80.3 (±5.8) | 79.5(±7.9) | 78.2(±6.6)
Multiple datasets |  |  |  |  |  |  |

Tasks. We consider multiple datasets, overall of 72,236 cases, to confirm the
efficacy of PR-Isomap: MNIST dataset of handwritten digits, Pneumonia, and
NSCLC-Radiogenomics.
Implementation Details. To apply PCA, t-SNE, and PHATE, we use the
implementation from [55], and [49], respectively. To apply UMAP for LD data
visualization, we apply the implementation from [20]. For standard Isomap, we
use [3] implemented by Jake Vanderplas (2011). We implement PR-Isomap based on
the standard implementation of Isomap by restricting k-nearest neighbors and
the shortest path on the manifold and updating it with the PR constraint.
Hyperparameters. For PCA and UMAP, we used the default hyperparameters, while
we changed the tSNE iteration number to 500 and a random state of 3. For the
binary classifiers, we used untuned classification models but we froze the
hyperparameters throughout the comparison for different DR methods and all
models to keep the integrity of the analysis. In Random Forest, we used 20
estimators, a maximum depth of the forest with four, and a random state of 45.
The training follows the standard cross-validation score with k = 8. For PR-
Isomap, $\Phi$ was a window function. h in the PR function was 16, and 10 for
CT and PET radiomics, while it was 70000 in the MNIST dataset and selected
based on the lower-bound amount for converting LD features. We used Gap
statistics and Elbow to identify the optimum number of LD spaces.

###  4.3 LD Representation of HD Manifolds

Classification evaluation and outcome prediction for the patients screening or
under treatment is a popular way of showing the predictability power of LD
radiomics. A supervised algorithm, random forest classifier, was used to
perform binary prediction of the final survival of the patients to assess the
prediction power of LD radiomics generated by the baselines. The model was
untuned but constantly across the comparison to perverse the integrity of the
evaluations. Table 4 reports the quantitative accuracy produced by PR-Isomap,
Isomap, tSNE, PCA, and PHATE. PR-isomap produces better accuracy than other DR
techniques. Thus, PR-Isomap can be used to reduce HD radiomics and facilitate
outcome predictions. In the context of visualizing 70,000 observations of
MNIST HD manifold projected into the LD 3D space utilizing various DR baseline
techniques. Notably, when employing PR-Isomap (6-a), it becomes evident that
the different categories exhibit closer proximity, despite possessing
relatively smaller between-class (global or extrinsic) distances. This can be
attributed to the additional PR constraint applied to sustain the uniformity
of the class distributions.
Moving to the results obtained through PHATE (6-e), a somewhat analogous
pattern emerges, albeit with a slightly less pronounced ability to effectively
group the within-class categories (intrinsic or local). In contrast, when
considering Isomap (6-b) and t-SNE (6-c), we observe distinct characteristics.
Isomap appears to exhibit lower uniformity in the distribution of data points,
while t-SNE showcases a propensity for within-class categories to cluster more
distinctly.
These observations collectively shed light on the diverse behaviors exhibited
by various DR methods when projecting MNIST HD data into an LD 3D space,
emphasizing the utility of PR-Isomap in fostering closer inter-category
proximity while upholding the uniformity of class distributions.

###  4.4 Cox Proportional Hazard Model to Predict Survival

In Cox modeling of OS on the multimodal PET/CT dataset, the C-statistic of CT
and PET LD radiomics were 0.68 and 0.67 ($95\%$ CI), with the long-rank
likelihood test p-value of $<0.005$, and $0.007$, respectively. This result
generated by PR-Isomap DR reduction outperformed other DR methods used for the
CT biomarker analysis while maintaining the separation of high- and low- risk
patients (Fig. 7). C-statistics of original Isomap while having LD CT and PET
radiomics yield 0.66 (p-value = $0.02$) and 0.67 (p-value = $0.02$), then the
results of the CT biomarkers considered to be the second-best accuracy after
our proposed method. Despite slightly higher C-statistics of PHATE for the
NSCLC Radiogenomics dataset, $0.70$ ($0.06$) for CT and $0.72$ ($0.02$) for
PET, the separation of high- and low- risk patients under treatment was
noticeably lower than other approaches (Fig. 7). In PET imaging markers, PCA
generates the highest C statistic, which is inconsistent with binary survival
prediction. the lowest C statistic belonged to CT-PCA and PET-tSNE radiomics
yielding $0.63$ ($0.04$) and $0.6$ ($0.02$), respectively. A model including
the CT-tSNE radiomic avatars resulted in C-statistics of $0.64$ with a p-value
of $0.09$. For both Lung01 and Lung03 datasets, PR-Isomap generated LD
radiomics alone yielded the highest C-Statistics of $0.59$ ($<0.0005$) and
$0.66$ ($0.002$), respectively. Also, the Kaplan-Meier survival curve
indicates higher separation of the high- and low-risk patients based on the
median of the hazard (Fig. 7). The Cox models of OS all datasets using PR-
Isomap radiomic avatars had statistically significant separation of the Kaplan
Meier curves for patients above versus below the median hazard compared to
other DR methods.

Figure 7: LD Multivariant Survival prediction. The proposed constrained
isometrically mapping DR method exhibits impressive discriminatory power in
distinguishing between high-risk and low-risk patients based on the median of
diverse imaging biomarkers. The incorporation of clinical information and
demographic variables using LD radiomics significantly improves the
separability of patient groups with distinct risk profiles. This capability
holds tremendous potential in accurately predicting patient survival across
the course of treatment.

##  5 Discussion

The HD to LD manifold learning, supported by extensive studies [51, 30, 31,
52, 53], proves instrumental in precision medicine, diagnosis, and optimizing
treatment planning. This approach effectively addresses the challenges of HD
multivariate data overload, notably in ’omics disciplines, by employing DR
methods. Such breakthroughs enable the prediction of genetic information
directly from CT images and the identification of gene mutations.
Additionally, certain imaging features correlate with gene overexpression, a
crucial therapeutic target. Given the complexity and dimensionality issues
inherent in these models, efficient DR techniques are indispensable. These
techniques aid in extracting non-invasive imaging biomarkers and diverse
genetic information, ultimately enhancing the precision of targeted therapy
and delivering clinical benefits [added]. Both genetic and radiomic data often
grapple with HD parameters, rendering DR methods indispensable for the precise
processing of HD medical data.

Several methods have been employed in previous studies, including PCA, MRMR
[54], clustering, MDS [2], Isomap [3], LLE [4], tSNE [55], and PHATE [49].
These methods project HD data to LD space and share some similarities. We used
KL-divergence, which is similarly used in MRMR, for comparison with tSNE. We
also considered MDS-based projection by comparing it with standard Isomap and
more recently developed PHATE [49]. PR-Isomap can be used to reduce the
dimensionality of data points from HD space while preserving the local
pairwise distances and structure of the HD data. In addition, Isomap and
similarly PR-Isomap exhibit clustering behavior. To measure the strength of
the proposed approach, PR-Isomap was compared to unsupervised DR methods with
fixed hyperparameters acquired in the tuning phase. Furthermore, our
unsupervised model was developed and validated using different independent
NSCLC datasets, and comparative analyses were conducted to demonstrate its
strength with respect to the state-of-the-art unsupervised feature reduction
methods, thereby increasing the novelty of our proposed model.

In precision medicine analysis, LD-projected radiomics models, employing
various DR methods, especially PR-Isomap, exhibit the capability to predict
outcome events. PR-Isomap consistently outperforms other state-of-the-art
unsupervised DR methods with the highest prediction accuracy 4. Although the
computational complexity of our proposed method aligns with that of the
original Isomap, the influence of spatial resolution, contrast, and unbalanced
data introduces minor variations in binary classification results. These
challenges are addressed using the Cox model and the "time to event" variable,
revealing PR-Isomap’s superior performance in survival prediction and
enhancing the accuracy of survival analysis compared to classification. PR-
Isomap surpasses all baseline DR models when utilizing multivariate Cox
proportional hazards models. However, in the first NSCLC data, PHATE exhibits
slightly higher C-statistics in the NSCLC-Radigenomics dataset. Yet, it
struggles with distinguishing high- and low-risk patients using LD radiomic
signatures (Fig. 7). Furthermore, the Kaplan-Meier survival curve of PR-Isomap
across DR models and different NSCLC datasets demonstrates superior separation
between high- and low-risk patients (see Fig. 7).

For prognostic analysis, our validation datasets encompass 510 patients,
primarily comprising CT imaging data, with only 130 patients possessing PET
radiomics data. Expanding the dataset to include genomic analysis is a
potential avenue for future research. Notably, the image sets within the
dataset exhibit heterogeneity in acquisition characteristics, such as contrast
enhancement, device model and manufacturer, slice thickness, and convolutional
kernel. Despite this potential source of noise, we successfully model survival
using patient imaging data. Furthermore, the absence of associations between
these imaging factors and survival strengthens our results’ credibility.
Nonetheless, to attain more robust prognostic results, larger and more
homogeneous datasets would be advantageous. Additionally, incorporating gene
expression data into the DR analysis could enhance the reliability of our
proposed method when dealing with HD data.

For diagnostic prediction using LD radiomic signatures, PR-Isomap outperformed
other state-of-the-art models, despite a considerably large cohort of patients
for pneumonia, more analysis, involving other types of diseases, could help
better understand and ensure the performance of the PR-Isomap. Such
investigations would provide a more holistic understanding of the model’s
capabilities and robustness in different clinical scenarios.

It is important to note the limitations of our study. PR-Isomap still suffers
from the same weaknesses as the original Isomap, even with the homogeneity
constraint, despite the alleviation of this problem. This is due to the nature
of the data, which often requires different transformations before projecting
it to LD space to overcome the nonuniform structure of data points on the
manifold. Other methods like UMAP using differential geometry in HD projection
may be beneficial in improving the data processing and interpretability of the
results. For that, mapping the HD manifolds to other spaces to redefine
uniformity could be an interesting direction to explore.

##  6 Conclusion

This study proposed a novel dimensionality reduction and manifold learning
method, PR-Isomap, which preserves information while projecting HD data to LD
space using a Parzen-Rosenblatt constraint. PR-Isomap was applied to create LD
projection of multiple datasets, overall 72,236 cases, and experimental
results showed significant superiority over other commonly used state-of-the-
art unsupervised DR methods to predict outcomes or diagnose diseases. PR-
Isomap method may be valuable for processing HD medical data, particularly
radiogenomics studies. This method has potential applications in radiogenomic-
driven diagnosis, treatment planning, and outcome prediction in HD multimodal
imaging.

## References

  * [1] Jolliffe, I. Principal component analysis (pp. 1094-1096). Springer Berlin Heidelberg. 28243034 pp. 28243034 (2011)
  * [2] Torgerson, W. Multidimensional scaling of similarity. Psychometrika. 30, 379-393 (1965)
  * [3] Tenenbaum, J., Silva, V. & Langford, J. A global geometric framework for nonlinear dimensionality reduction. Science. 290, 2319-2323 (2000)
  * [4] Roweis, S. & Saul, L. Nonlinear dimensionality reduction by locally linear embedding. Science. 290, 2323-2326 (2000)
  * [5] Belkin, M. & Niyogi, P. Laplacian eigenmaps for dimensionality reduction and data representation. Neural Computation. 15, 1373-1396 (2003)
  * [6] Donoho, D. & Grimes, C. Hessian eigenmaps: Locally linear embedding techniques for high-dimensional data. Proceedings Of The National Academy Of Sciences. 100, 5591-5596 (2003)
  * [7] Sammon, J. A nonlinear mapping for data structure analysis. IEEE Transactions On Computers. 100, 401-409 (1969)
  * [8] Gou, J., Yang, Y., Yi, Z., Lv, J., Mao, Q. & Zhan, Y. Discriminative globality and locality preserving graph embedding for dimensionality reduction. Expert Systems With Applications. 144 pp. 113079 (2020)
  * [9] Lee, J., Verleysen, M. & Others Nonlinear dimensionality reduction. (Springer,2007)
  * [10] Zhao, X. & Zhang, S. Facial expression recognition using local binary patterns and discriminant kernel locally linear embedding. EURASIP Journal On Advances In Signal Processing. 2012 pp. 1-9 (2012)
  * [11] Maaten, L. & Hinton, G. Visualizing data using t-SNE.. Journal Of Machine Learning Research. 9 (2008)
  * [12] Hinton, G. & Roweis, S. Stochastic neighbor embedding. Advances In Neural Information Processing Systems. 15 (2002)
  * [13] Parzen, E. On estimation of a probability density function and mode. The Annals Of Mathematical Statistics. 33, 1065-1076 (1962)
  * [14] Rosenblatt, M. Remarks on some nonparametric estimates of a density function. The Annals Of Mathematical Statistics. pp. 832-837 (1956)
  * [15] Zhao, X. & Zhang, S. Facial expression recognition using local binary patterns and discriminant kernel locally linear embedding. EURASIP Journal On Advances In Signal Processing. 2012 pp. 1-9 (2012)
  * [16] Lewandowski, M., Rincon, J., Makris, D. & Nebel, J. Temporal extension of laplacian eigenmaps for unsupervised dimensionality reduction of time series. 2010 20th International Conference On Pattern Recognition. pp. 161-164 (2010)
  * [17] Jenkins, O. & Matarić, M. A spatio-temporal extension to isomap nonlinear dimension reduction. Proceedings Of The Twenty-first International Conference On Machine Learning. pp. 56 (2004)
  * [18] Lin, R., Liu, C., Yang, M., Ahuja, N. & Levinson, S. Learning nonlinear manifolds from time series. Computer Vision–ECCV 2006: 9th European Conference On Computer Vision, Graz, Austria, May 7-13, 2006. Proceedings, Part II 9. pp. 245-256 (2006)
  * [19] Guo, Y., Guo, H. & Yu, S. Co-sne: Dimensionality reduction and visualization for hyperbolic data. Proceedings Of The IEEE/CVF Conference On Computer Vision And Pattern Recognition. pp. 21-30 (2022)
  * [20] McInnes, L., Healy, J. & Melville, J. Umap: Uniform manifold approximation and projection for dimension reduction. ArXiv Preprint ArXiv:1802.03426. (2018)
  * [21] Alamgir, M. & Von Luxburg, U. Shortest path distance in random k-nearest neighbor graphs. ArXiv Preprint ArXiv:1206.6381. (2012)
  * [22] Cayton, L. Algorithms for manifold learning. Univ. Of California At San Diego Tech. Rep. 12, 1 (2005)
  * [23] Babaeian, A. Part I-Constrained Shortest-Path For Manifold Learning And Multiple Manifold Clustering Part II-Community Detection In Large Graphs; Analysis, Design And Implementation. (UC San Diego,2017)
  * [24] Cox, T. & Cox, M. Multidimensional scaling. (Chapman,2001)
  * [25] Dadkhahi, H., Duarte, M. & Marlin, B. Out-of-sample extension for dimensionality reduction of noisy time series. IEEE Transactions On Image Processing. 26, 5435-5446 (2017)
  * [26] Najafi, A., Joudaki, A. & Fatemizadeh, E. Nonlinear dimensionality reduction via path-based isometric mapping. IEEE Transactions On Pattern Analysis And Machine Intelligence. 38, 1452-1464 (2015)
  * [27] Minka, T. Automatic choice of dimensionality for PCA. Advances In Neural Information Processing Systems. 13 (2000)
  * [28] Hinton, G. & Maaten, L. Visualizing data using t-sne journal of machine learning research. (2008)
  * [29] Bakr, S., Gevaert, O., Echegaray, S., Ayers, K., Zhou, M., Shafiq, M., Zheng, H., Benson, J., Zhang, W., Leung, A. & Others A radiogenomic dataset of non-small cell lung cancer. Scientific Data. 5, 1-9 (2018)
  * [30] Gevaert, O., Xu, J., Hoang, C., Leung, A., Xu, Y., Quon, A., Rubin, D., Napel, S. & Plevritis, S. Non–small cell lung cancer: identifying prognostic imaging biomarkers by leveraging public gene expression microarray data—methods and preliminary results. Radiology. 264, 387-396 (2012)
  * [31] Clark, K., Vendt, B., Smith, K., Freymann, J., Kirby, J., Koppel, P., Moore, S., Phillips, S., Maffitt, D., Pringle, M. & Others The Cancer Imaging Archive (TCIA): maintaining and operating a public information repository. Journal Of Digital Imaging. 26 pp. 1045-1057 (2013)
  * [32] Wilson, R. & Devaraj, A. Radiomics of pulmonary nodules and lung cancer.. Pulmonary Nodules And Lung Cancer. pp. 128 (2017)
  * [33] Van Griethuysen, J., Fedorov, A., Parmar, C., Hosny, A., Aucoin, N., Narayan, V., Beets-Tan, R., Fillion-Robin, J., Pieper, S. & Aerts, H. Computational radiomics system to decode the radiographic phenotype. Cancer Research. 77, e104-e107 (2017)
  * [34] Cohen, J., Morrison, P. & Dao, L. COVID-19 image data collection. ArXiv Preprint ArXiv:2003.11597. (2020)
  * [35] Wang, X., Peng, Y., Lu, L., Lu, Z., Bagheri, M. & Summers, R. Chestx-ray8: Hospital-scale chest x-ray database and benchmarks on weakly-supervised classification and localization of common thorax diseases. Proceedings Of The IEEE Conference On Computer Vision And Pattern Recognition. pp. 2097-2106 (2017)
  * [36] Ozturk, T., Talo, M., Yildirim, E., Baloglu, U., Yildirim, O. & Acharya, U. Automated detection of COVID-19 cases using deep neural networks with X-ray images. Computers In Biology And Medicine. 121 pp. 103792 (2020)
  * [37] Gunraj, H., Wang, L. & Wong, A. Covidnet-ct: A tailored deep convolutional neural network design for detection of covid-19 cases from chest ct images. Frontiers In Medicine. 7 pp. 608525 (2020)
  * [38] Wang, L., Lin, Z. & Wong, A. Covid-net: A tailored deep convolutional neural network design for detection of covid-19 cases from chest x-ray images. Scientific Reports. 10, 19549 (2020)
  * [39] Wang, L., Wong, A., Lin, Z., Lee, J., McInnis, P., Chung, A. & A. Figure 1 COVID-19 Chest X-ray Dataset Initiative. (2020), https://github.com/agchung/Figure1-COVID-chestxray-dataset
  * [40] Yousefi, B., LaRiviere, M., Cohen, E., Buckingham, T., Yee, S., Black, T., Chien, A., Noël, P., Hwang, W., Katz, S. & Others Combining radiomic phenotypes of non-small cell lung cancer with liquid biopsy data may improve prediction of response to EGFR inhibitors. Scientific Reports. 11, 9984 (2021)
  * [41] Yousefi, B., Jahani, N., LaRiviere, M., Cohen, E., Hsieh, M., Luna, J., Chitalia, R., Thompson, J., Carpenter, E., Katz, S. & Others Correlative hierarchical clustering-based low-rank dimensionality reduction of radiomics-driven phenotype in non-small cell lung cancer. Medical Imaging 2019: Imaging Informatics For Healthcare, Research, And Applications. 10954 pp. 278-285 (2019)
  * [42] Hershman, M., Yousefi, B., Serletti, L., Galperin-Aizenberg, M., Roshkovan, L., Luna, J., Thompson, J., Aggarwal, C., Carpenter, E., Kontos, D. & Others Impact of interobserver variability in manual segmentation of non-small cell lung cancer (NSCLC) applying low-rank radiomic representation on computed tomography. Cancers. 13, 5985 (2021)
  * [43] Yousefi, B., Sharifipour, H. & Maldague, X. A diagnostic biomarker for breast cancer screening via Hilbert embedded deep low-rank matrix approximation. IEEE Transactions On Instrumentation And Measurement. 70 pp. 1-9 (2021)
  * [44] Vigil, N., Nouri, B., Fernandes, H., Ibarra-Castanedo, C., Maldague, X. & Yousefi, B. Convex Factorization Embedding Thermography for Breast Cancer Diagnostic. IEEE Open Journal Of Instrumentation And Measurement. 1 pp. 1-8 (2022)
  * [45] Darvish, M., Trask, R., Tallon, P., Khansari, M., Ren, L., Hershman, M. & Yousefi, B. AI-Enabled Lung Cancer Prognosis. . arXiv preprint arXiv:2402.09476.
  * [46] Winther, H., Laser, H., Gerbel, S., Maschke, S., Hinrichs, J., Vogel-Claussen, Jens, W., Frank K., H., M., M. & Meyer, B. COVID-19 Image Repository.
  * [47] Yousefi, B., Kawakita, S., Amini, A., Akbari, H., Advani, S., Akhloufi, M., Maldague, X. & Ahadian, S. Impartially validated multiple deep-chain models to detect COVID-19 in chest X-ray using latent space radiomics. Journal Of Clinical Medicine. 10, 3100 (2021)
  * [48] LeCun, Y. The MNIST database of handwritten digits. Http://yann. Lecun. Com/exdb/mnist/. (1998)
  * [49] Moon, K., Dijk, D., Wang, Z., Gigante, S., Burkhardt, D., Chen, W., Yim, K., Elzen, A., Hirn, M., Coifman, R. & Others Visualizing structure and transitions in high-dimensional biological data. Nature Biotechnology. 37, 1482-1492 (2019)
  * [50] Silva, V. & Tenenbaum, J. Global versus local methods in nonlinear dimensionality reduction. Advances In Neural Information Processing Systems. 15 (2002)
  * [51] Bakr, S., Gevaert, O., Echegaray, S., Ayers, K., Zhou, M., Shafiq, M., Zheng, H., Benson, J., Zhang, W., Leung, A. & Others Data for NSCLC Radiogenomics Collection. The Cancer Imaging Archive. 2017.
  * [52] Aerts, H., Velazquez, E., Leijenaar, R. & Others Data from NSCLC-radiomics-genomics. Cancer Imaging Arch. (2017)
  * [53] Aerts, H., Velazquez, E., Leijenaar, R., Parmar, C., Grossmann, P., Carvalho, S., Bussink, J., Monshouwer, R., Haibe-Kains, B., Rietveld, D. & Others Decoding tumour phenotype by noninvasive imaging using a quantitative radiomics approach. Nature Communications. 5, 4006 (2014)
  * [54] Radovic, M., Ghalwash, M., Filipovic, N. & Obradovic, Z. Minimum redundancy maximum relevance feature selection approach for temporal gene expression data. BMC Bioinformatics. 18, 1-14 (2017)
  * [55] Maaten, L. & Hinton, G. Visualizing data using t-SNE.. Journal Of Machine Learning Research. 9 (2008)