#  A Comparison of Deep Learning Architectures for Optical Galaxy Morphology
Classification

Ezra Fielding  Department of Computer Science
University of the Western Cape
Cape Town, South Africa
0000-0002-7936-0222  Clement N. Nyirenda  Department of Computer Science
University of the Western Cape
Cape Town, South Africa
0000-0002-4181-0478  Mattia Vaccari  Department of Physics & Astronomy
University of the Western Cape
Cape Town, South Africa
0000-0002-6748-0577

###### Abstract

The classification of galaxy morphology plays a crucial role in understanding
galaxy formation and evolution. Traditionally, this process is done manually.
The emergence of deep learning techniques has given room for the automation of
this process. As such, this paper offers a comparison of deep learning
architectures to determine which is best suited for optical galaxy morphology
classification. Adapting the model training method proposed by Walmsley et al
in 2021, the Zoobot Python library is used to train models to predict Galaxy
Zoo DECaLS decision tree responses, made by volunteers, using EfficientNet B0,
DenseNet121 and ResNet50 as core model architectures. The predicted results
are then used to generate accuracy metrics per decision tree question to
determine architecture performance. DenseNet121 was found to produce the best
results, in terms of accuracy, with a reasonable training time. In future,
further testing with more deep learning architectures could prove beneficial.

###### Index Terms:

Deep Learning, Optical Galaxy Morphology, Classification, Astronomy

##  I Introduction

Morphology plays a crucial role in the investigation of galaxy formation and
evolution [1]. As such, accurate and detailed morphological classifications
are required to advance research in this area. This process of classification
has historically been performed visually by scientists. However, modern
astronomical surveys have been producing more detailed images than scientists
can visually classify [1]. Projects, such as Galaxy Zoo, have attempted to
solve this problem by asking so called ‘citizen scientists’, volunteer members
of the public, to provide classifications through a web interface [2]. While
these efforts have been largely successful and have greatly benefited the
scientific community, it is only a matter of time until the output of modern
astronomical surveys exceeds efforts which rely solely on human intelligence
[3, 4]. In other words, projects which rely on humans to visually classify
morphology, will not be able to keep up with the sheer volume of data produced
by next generation surveys, such as the Square Kilometer Array (SKA),
resulting in missed discoveries and insights [5, 6].

Different approaches have been used to solve the problem of classifying
galaxies according to their morphology. In the original Galaxy Zoo study,
approximately one hundred thousand participants collectively provided more
than 4×$10^{7}$ individual classifications for nearly one million galaxies
from the Sloan Digital Sky Survey (SDSS) [2]. Volunteers submitted
classifications on the Galaxy Zoo website, where they were shown an image
cutout of an area of sky centered on a galaxy. These raw classifications were
then aggregated and used to form a final catalogue. These Galaxy Zoo results
were consistent with those for a subset of SDSS classified by professional
astronomers [2]. Subsequent projects such as Radio Galaxy Zoo and Galaxy Zoo
DECaLS have also found success in using ‘citizen scientists’, to visually
classify and inspect galaxies [7, 1]. The results of these Galaxy Zoo projects
have been used as the basis for a large number of studies on galaxy morphology
[3].

However, despite these successes, as the amount of survey data in need of
classification continues to increase, approaches relying mainly on human
classification and inspection will require too much time and effort to be
viable [4]. For example, the recent Galaxy Zoo DECaLS project would have taken
approximately 8 years to collect the standard 40 classifications for each of
the 311488 suitable galaxies with the current number of volunteers available
[1]. By that time, new surveys would have begun. This highlights the need for
classification methods which do not rely heavily on human input.

A fair amount of work has been done in implementing deep learning algorithms
for the purpose of classifying galaxy morphology [6]. A large portion of this
work makes use of artificial neural networks to predict galaxy morphology
based on images. These networks are trained using a labeled subset of existing
surveys. Convolutional neural networks (CNNs) have been shown to be
particularly useful and accurate in the classification of astronomical
objects, provided the sufficient availability of labeled datasets [8, 9, 10,
11, 12, 13]. Of particular interest is the work of Becker et al, which
comparatively investigates the performance of different CNN architectures
applied to radio galaxy classification. In their work, MCRGNET, RADIO GALAXY
ZOO, and CONVXPRESS (a novel classifier) were the architectures which had the
best balance of computational requirements and recognition performance [14].
Walmsley et al used Galaxy Zoo DECaLS response data to train an EfficientNet
B0 architecture to predict the responses to the GDZ-5 decision tree used by
the Galaxy Zoo project to classify galaxies [1].

Most existing works have focused on comparing state-of-the-art deep learning
architectures on radio-spectrum images. In terms of optical-spectrum images,
the work in [1], which is the leading one, only uses the EfficientNet B0 and
does not present a comparison with other prominent deep learning
architectures. Therefore, this paper aims to assess the performance of three
different, popular deep learning architectures in predicting volunteer vote
fractions for optical-spectrum galaxy images from the recent Galaxy Zoo DECaLS
project. EfficientNet B0 [15], DenseNet121 [16] and ResNet50 [17] will
compared in this work. Performance criteria will include training time and
accuracy metrics such as precision, recall and F1-score.

The rest of this paper is organised as follows. Section II presents a brief
overview of the deep learning architectures used in this work. Section III
presents the materials and methods used in this work, while Section IV
presents the results and their discussion. Section V draws the conclusion.

##  II Deep Learning Architectures Used in this Work

Three popular deep learning architectures were selected for this comparison.
DenseNet121 and ResNet50 were selected alongside EfficientNet B0, which is the
default architecture included with Walsmley’s Zoobot library [18]. These
architectures were chosen based on their size, popularity, and ease of
integration into the existing Zoobot library code base. These architectures
are briefly described as follows.

#### EfficientNet B0

EfficientNet is a CNN-based architecture that uses a scaling method that
uniformly scales all dimensions of depth, width, and resolution using a
compound coefficient. The EfficientNet family of models are, therefore,
optimised for both accuracy and efficiency, with the aim of performing better
than traditional CNNs [15]. Apart from the classification of galaxies,
EfficientNet has been applied to many other areas such as Covid-19 diagnosis
[19], automated fruit recognition [20] and multi-label classification of
fundus images [21].

#### DenseNet121

The DenseNet family of models differentiates itself from traditional CNN’s by
connecting each layer to every other layer in a feed-forward fashion [16].
This works to reduce the number of parameters and alleviate the vanishing
gradient problem. DenseNets achieve high performance while requiring less
memory and computation when compared to other CNNs. DenseNet has also been
applied to many areas such as remote sensing image classification [22],
malware detection [23], and classification of Alzheimer’s disease [24].

#### ResNet50

The ResNet Family of models reformulate the model layers as learning residual
functions with reference to the layer inputs. This results in models which are
easier to optimise and which show a gain in accuracy to increased model depth
[17]. Some of the recent applications include crime prediction [25, 26] and
skin cancer classification [27].

Table I presents the generic model that illustrates how the architectures were
implemented as layers in the Zoobot library’s default model. The top fully
connected layer for each architecture was removed, and replaced with by the
three augmentation layers which feed directly in to the network.

Table I: Generic Model Architecture Layer (type) | Output Shape
---|---
RandomRotation | (None, 300, 300, 1)
RandomFlip | (None, 300, 300, 1)
RandomCrop | (None, 224, 224, 1)
Core Architecturea | (None, 7, 7, )
GlobalAveragePooling2D | (None, )
Dense | (None, 34)
aThis layer will be swapped out to test each architecture.

With the top fully connected layer removed, EfficientNet B0 resulted in
approximately 4.0 million trainable parameters, while DenseNet121 and ResNet50
resulted in approximately 7.0 million and 23.6 million trainable parameters
respectively. The final dense layer of the core model architecture was
modified to give predictions smoothly, from 1 to 100 using softmax activation.
This was done in order to provide plausible Dirichlet posteriors for galaxy
morphology [1].

##  III Materials and Methods

The method used for this architecture comparison, follows closely with that
presented by Walmsley et al, 2021 for predicting Galaxy Zoo GZD-5 decision
tree volunteer responses [1]. However, instead of training three models per
architecture to form an ensemble network, only one model per architecture was
trained. The training was also conducted at a decreased batch size of 64. A
modified version of the Zoobot library was used for data pre-processing, model
creation, model training, and making predictions [18].

###  III-A Description of the Dataset

The Galaxy Zoo DECaLS dataset includes 253286 424x424 pixel color images (3
channels in the g, r, and z bands) [1]. Volunteer response data is provided
for each galaxy. This response data includes the total votes and vote
fractions for each question in the GZD-5 decision tree. Any galaxy with less
than 3 volunteer classifications was removed resulting in a 1.5% reduction in
the dataset size. This left a total of 249581 images for training and testing.
A standard random 80/20 train-test split was applied to the dataset. This
resulted in a training dataset containing 199664 galaxies and a test dataset
containing 49917 galaxies. The test dataset will serve as the validation
dataset.

The GZD-5 decision tree used for this dataset includes ten questions which
have a different number of options depending on the question. A chart
explaining the GZD-5 decision tree can be found in Fig. 1.

Figure 1: The GZD-5 Decision Tree [1].

###  III-B Data Processing

The color information of each image was removed by taking an average over its
3 channels to avoid biasing the morphology predictions [28]. The images were
then resized and saved as 300x300x1 matrices. Further augmentations were done
when images were loaded into the model. Random horizontal and vertical flips
were first applied, followed by a random rotation in the range (0, $\pi$). The
image was finally cropped to a size of 224x224 pixels about a random centroid.
These augmentations were applied as the images were fed into the model.

###  III-C Training and Testing

Tensorflow and Python 3.8 were used for model training. The three models were
identically trained using the Adam optimizer and loss was calculated by taking
the negative log likelihood of observing volunteer answers to all questions
from Dirichlet-Multinomial distributions [1]. Training was done with a batch
size of 64 and ended after the model loss failed to improve after 10
consecutive epochs. A batch size of 64 was chosen due to the memory
limitations of the training hardware. The test dataset was repeated once to
reduce variance from random dropouts and augmentations. Training was completed
on the South African Center for High Performance Computing’s (CHPC) Legau
Cluster on two NVIDIA V100 GPUs. Training prototyping and debugging, along
with model prediction and evaluation was performed on the ilifu cloud
computing facility.

To test and evaluate model performance, a new 80/20 train-test split was
applied to the dataset. The new test dataset which resulted from this split
was used for evaluation. A galaxy was only used to evaluate a question if 50%
or more of the volunteers shown that galaxy was asked that specific question.
Five forward passes per image were made with each model to account for random
dropout configurations and augmentations. This resulted in five sets of
predicted Dirichlet-Multinomial posteriors per image per model. The mean
posterior of the total votes for each model were then calculated and recorded,
to be used for calculating the predicted vote fraction for each question
option.

The highest predicted vote fraction for a question was rounded to 1, while the
remaining fractions were rounded to 0. This was done to produce a set of
discrete classifications for each of the GZD-5 decision tree questions. The
discrete classifications for each question were then used to generate accuracy
metrics, namely precision, recall and F1-score.

#### Precision

This is defined as the ratio of images that was correctly defined as belonging
to the positive class, to the total number of images that were classified as
belonging to the positive class.

#### Recall

This refers to the ratio of the number of images that were correctly
classified as belonging to the positive class, to the total number of images
in the positive class.

#### F1-Score

This is the weighted average of recall and precision.

The weighted average metrics are considered for each question. The weighted
average metrics are calculated by calculating the metrics for each question
option, and finding their average weighted by support. This allows the
accuracy metrics to account for label imbalance. Table II shows the support
for the metrics of each decision tree question.

Table II: Question Support Question | Support
---|---
smooth-or-featured | 49917
disk-edge-on | 15445
has-spiral-arms | 11380
bar | 11380
bulge-size | 11380
how-rounded | 32526
edge-on-bulge | 2475
spiral-winding | 7499
spiral-arm-count | 7499
merging | 49247

##  IV Results and Discussion

The following results were produced after running the experiment. Table III
shows the total number of epochs run during training and the total training
time.

Table III: Statistics from Training | EfficientNet B0 | DenseNet121 | ResNet50
---|---|---|---
Total Epochs | 52 | 37 | 37
Total Time (hours) | 19.795 | 11.723 | 9.642

DenseNet121 and ResNet50 both ran for 37 epochs, while EfficientNet B0 ran for
52 epochs. Despite being the largest model by number of parameters, ResNet50
trained for a total of 9.642 hours, while DenseNet121 completed training after
11.723 hours. EfficientNet B0 spent the longest time training, at 19.795
hours.

The weighted average metrics for each question for each of the three models
were placed into tables for comparison. Table IV presents a comparison of the
precision metric, Table V presents a comparison of the recall metric, and
Table VI presents a comparison of the F1-score metric. The weighted average
for all questions is displayed for each model at the bottom of each comparison
table.

Table IV: Precision Comparison Question | Core Model Architecture
---|---
Name | EfficientNet B0 | DenseNet121 | ResNet50
smooth-or-featured | 0.877 | 0.880 | 0.869
disk-edge-on | 0.954 | 0.957 | 0.955
has-spiral-arms | 0.891 | 0.893 | 0.868
bar | 0.697 | 0.698 | 0.673
bulge-size | 0.684 | 0.690 | 0.675
how-rounded | 0.870 | 0.876 | 0.867
edge-on-bulge | 0.794 | 0.826 | 0.785
spiral-winding | 0.685 | 0.695 | 0.677
spiral-arm-count | 0.664 | 0.695 | 0.652
merging | 0.849 | 0.851 | 0.845
Weighted Average | 0.838 | 0.843 | 0.831
Best model result highlighted in bold. Worst model result underlined.
Table V: Recall Comparison Question | Core Model Architecture
---|---
Name | EfficientNet B0 | DenseNet121 | ResNet50
smooth-or-featured | 0.876 | 0.881 | 0.863
disk-edge-on | 0.955 | 0.957 | 0.955
has-spiral-arms | 0.887 | 0.888 | 0.875
bar | 0.715 | 0.710 | 0.696
bulge-size | 0.694 | 0.696 | 0.691
how-rounded | 0.867 | 0.875 | 0.864
edge-on-bulge | 0.820 | 0.825 | 0.816
spiral-winding | 0.698 | 0.707 | 0.690
spiral-arm-count | 0.676 | 0.674 | 0.664
merging | 0.881 | 0.882 | 0.880
Weighted Average | 0.848 | 0.851 | 0.841
Best model result highlighted in bold. Worst model result underlined.
Table VI: F1-Score Comparison Question | Core Model Architecture
---|---
Name | EfficientNet B0 | DenseNet121 | ResNet50
smooth-or-featured | 0.876 | 0.881 | 0.865
disk-edge-on | 0.955 | 0.956 | 0.955
has-spiral-arms | 0.889 | 0.890 | 0.869
bar | 0.693 | 0.677 | 0.662
bulge-size | 0.669 | 0.674 | 0.672
how-rounded | 0.867 | 0.875 | 0.865
edge-on-bulge | 0.787 | 0.794 | 0.784
spiral-winding | 0.687 | 0.696 | 0.678
spiral-arm-count | 0.654 | 0.668 | 0.642
merging | 0.851 | 0.855 | 0.848
Weighted Average | 0.836 | 0.840 | 0.829
Best model result highlighted in bold. Worst model result underlined.

From the model metric comparison tables, as far as accuracy is concerned, it
is clear that DenseNet121 produces the best predictions for the given dataset.
ResNet50 generally performs the worst out of the three, and EfficientNet B0
lies in the middle between the two in terms of accuracy.

Taking into account accuracy and training time, DenseNet121 appears to be the
most optimal architecture for this dataset. While ResNet50 had the fastest
training time, it generally performed the worst for each question in each
metric. EfficientNet B0, on the other hand, performed second best in terms of
accuracy. However, it took significantly longer to train compared to the other
two architectures. DenseNet121 provides the best accuracy performance at a
reasonable model training time.

##  V Conclusion

In conclusion, three popular deep learning architectures (EfficientNet B0,
DenseNet121, and ResNet50) were evaluated and compared in terms of accuracy
and training performance with regards to classifying optical galaxy
morphology. Performance criteria included training time and accuracy metrics
such as precision, recall and F1-score. The method for training outlined by
Walmsley et al in [1] was followed with some minor changes. In terms of
accuracy, DenseNet121 was the best performing architecture with the weighted
averages of 0.843, 0.851 and 0.840 for precision, recall and F1-score
respectively. On the other hand, ResNet50 exhibited the lowest performance
with the weighted average scores of 0.831, 0.841, and 0.829 for precision,
recall and F1-score respectively. The training time for DenseNet121 was 11.723
hours, which was approximately 2 hours more than ResNet50. On the other hand,
EfficientNet B0 took 19.795 hours. Taking accuracy, training time and total
parameters into account, DenseNet121 is the best model to use for this use
case.

In future, further testing with a larger selection of deep learning
architectures could prove useful and beneficial for the morphology
classification of optical galaxies. Applying these architectures to other
Galaxy Zoo projects could also prove insightful.

## Acknowledgment

The authors would like to thank the Galaxy Zoo volunteers, whose efforts made
this dataset possible. Their efforts are individually acknowledged at
https://authors.galaxyzoo.org/. Special thanks to Mike Walmsley for granting
access to the Galaxy Zoo DECaLS dataset and making the Zoobot library publicly
available. The authors acknowledge the Centre for High Performance Computing
(CHPC), South Africa, and the ilifu cloud computing facility – www.ilifu.ac.za
- for providing computational resources to this research project.

## References

  * [1] M. Walmsley, C. Lintott, T. Geron, S. Kruk, C. Krawczyk, K. W. Willett, S. Bamford, W. Keel, L. S. Kelvin, L. Fortson _et al._ , “Galaxy zoo decals: Detailed visual morphology measurements from volunteers and deep learning for 314,000 galaxies,” _arXiv preprint arXiv:2102.08414_ , 2021\.
  * [2] C. J. Lintott, K. Schawinski, A. Slosar, K. Land, S. Bamford, D. Thomas, M. J. Raddick, R. C. Nichol, A. Szalay, D. Andreescu _et al._ , “Galaxy zoo: morphologies derived from visual inspection of galaxies from the sloan digital sky survey,” _Monthly Notices of the Royal Astronomical Society_ , vol. 389, no. 3, pp. 1179–1189, 2008.
  * [3] K. L. Masters and G. Z. Team, “Twelve years of galaxy zoo,” _Proceedings of the International Astronomical Union_ , vol. 14, no. S353, pp. 205–212, 2019\.
  * [4] T. J. Galvin, M. T. Huynh, R. P. Norris, X. Wang, E. Hopkins, K. Polsterer, N. O. Ralph, A. N. O’Brien, and G. H. Heald, “Cataloguing the radio-sky with unsupervised machine learning: a new approach for the ska era,” _Monthly Notices of the Royal Astronomical Society_ , vol. 497, no. 3, pp. 2730–2758, 2020.
  * [5] P. E. Dewdney, P. J. Hall, R. T. Schilizzi, and T. J. L. Lazio, “The square kilometre array,” _Proceedings of the IEEE_ , vol. 97, no. 8, pp. 1482–1496, 2009.
  * [6] C. J. Fluke and C. Jacobs, “Surveying the reach and maturity of machine learning and artificial intelligence in astronomy,” _Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery_ , vol. 10, no. 2, p. e1349, 2020.
  * [7] J. K. Banfield, O. Wong, K. W. Willett, R. P. Norris, L. Rudnick, S. S. Shabala, B. D. Simmons, C. Snyder, A. Garon, N. Seymour _et al._ , “Radio galaxy zoo: host galaxies and radio morphologies derived from visual inspection,” _Monthly Notices of the Royal Astronomical Society_ , vol. 453, no. 3, pp. 2326–2340, 2015.
  * [8] S. Dieleman, K. W. Willett, and J. Dambre, “Rotation-invariant convolutional neural networks for galaxy morphology prediction,” _Monthly notices of the royal astronomical society_ , vol. 450, no. 2, pp. 1441–1459, 2015.
  * [9] W. Alhassan, A. Taylor, and M. Vaccari, “The first classifier: compact and extended radio galaxy classification using deep convolutional neural networks,” _Monthly Notices of the Royal Astronomical Society_ , vol. 480, no. 2, pp. 2085–2093, 2018.
  * [10] H. Domínguez Sánchez, M. Huertas-Company, M. Bernardi, D. Tuccillo, and J. Fischer, “Improving galaxy morphologies for sdss with deep learning,” _Monthly Notices of the Royal Astronomical Society_ , vol. 476, no. 3, pp. 3661–3676, 2018.
  * [11] C. J. Burke, P. D. Aleo, Y.-C. Chen, X. Liu, J. R. Peterson, G. H. Sembroski, and J. Y.-Y. Lin, “Deblending and classifying astronomical sources with mask r-cnn deep learning,” _Monthly Notices of the Royal Astronomical Society_ , vol. 490, no. 3, pp. 3952–3965, 2019.
  * [12] V. Lukic, M. Brüggen, J. K. Banfield, O. Wong, L. Rudnick, R. P. Norris, and B. Simmons, “Radio galaxy zoo: compact and extended radio source classification with deep learning,” _Monthly Notices of the Royal Astronomical Society_ , vol. 476, no. 1, pp. 246–260, 2018.
  * [13] Z. Ma, H. Xu, J. Zhu, D. Hu, W. Li, C. Shan, Z. Zhu, L. Gu, J. Li, C. Liu _et al._ , “A machine learning based morphological classification of 14,245 radio agns selected from the best–heckman sample,” _The Astrophysical Journal Supplement Series_ , vol. 240, no. 2, p. 34, 2019.
  * [14] B. Becker, M. Vaccari, M. Prescott, and T. Grobler, “Cnn architecture comparison for radio galaxy classification,” _Monthly Notices of the Royal Astronomical Society_ , vol. 503, no. 2, pp. 1828–1846, 2021.
  * [15] M. Tan and Q. Le, “Efficientnet: Rethinking model scaling for convolutional neural networks,” in _International Conference on Machine Learning_. PMLR, 2019, pp. 6105–6114.
  * [16] G. Huang, Z. Liu, L. Van Der Maaten, and K. Q. Weinberger, “Densely connected convolutional networks,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 4700–4708.
  * [17] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [18] M. Walmsley, “Zoobot,” https://github.com/mwalmsley/zoobot, 2021.
  * [19] G. Marques, D. Agarwal, and I. de la Torre Díez, “Automated medical diagnosis of covid-19 through efficientnet convolutional neural network,” _Applied soft computing_ , vol. 96, p. 106691, 2020.
  * [20] L. T. Duong, P. T. Nguyen, C. Di Sipio, and D. Di Ruscio, “Automated fruit recognition using efficientnet and mixnet,” _Computers and Electronics in Agriculture_ , vol. 171, p. 105326, 2020.
  * [21] J. Wang, L. Yang, Z. Huo, W. He, and J. Luo, “Multi-label classification of fundus images with efficientnet,” _IEEE Access_ , vol. 8, pp. 212 499–212 508, 2020.
  * [22] W. Tong, W. Chen, W. Han, X. Li, and L. Wang, “Channel-attention-based densenet network for remote sensing image scene classification,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 13, pp. 4121–4132, 2020.
  * [23] J. Hemalatha, S. A. Roseline, S. Geetha, S. Kadry, and R. Damaševičius, “An efficient densenet-based deep learning model for malware detection,” _Entropy_ , vol. 23, no. 3, p. 344, 2021.
  * [24] J. Ruiz, M. Mahmud, M. Modasshir, M. S. Kaiser, f. t. Alzheimer’s Disease Neuroimaging Initiative _et al._ , “3d densenet ensemble in 4-way classification of alzheimer’s disease,” in _International Conference on Brain Informatics_. Springer, 2020, pp. 85–96.
  * [25] T. Matereke, C. Nyirenda, and M. Ghaziasgar, “A comparative evaluation of spatio temporal deep learning techniques for crime prediction,” in _Proceedings of the IEEE Africon Conference_. IEEE, 2021.
  * [26] B. Wang, P. Yin, A. L. Bertozzi, P. J. Brantingham, S. J. Osher, and J. Xin, “Deep learning for real-time crime forecasting and its ternarization,” _Chinese Annals of Mathematics, Series B_ , vol. 40, no. 6, pp. 949–966, 2019\.
  * [27] N. Gouda and J. Amudha, “Skin cancer classification using resnet,” in _2020 IEEE 5th International Conference on Computing Communication and Automation (ICCCA)_. IEEE, 2020, pp. 536–541.
  * [28] M. Walmsley, L. Smith, C. Lintott, Y. Gal, S. Bamford, H. Dickinson, L. Fortson, S. Kruk, K. Masters, C. Scarlata _et al._ , “Galaxy zoo: probabilistic morphology through bayesian cnns and active learning,” _Monthly Notices of the Royal Astronomical Society_ , vol. 491, no. 2, pp. 1554–1574, 2020.


# Unsupervised Instance and Subnetwork Selection for Network Data

Lin Zhang  International Digital
Economy Academy (IDEA)
Shenzhen, China
<EMAIL_ADDRESS>Nicholas Moskwa  Biological Sciences
University at Albany - SUNY
New York, USA
<EMAIL_ADDRESS>Melinda Larsen  Biological Sciences
University at Albany - SUNY
New York, USA
<EMAIL_ADDRESS>Petko Bogdanov  Computer Science
University at Albany - SUNY
New York, USA
<EMAIL_ADDRESS>
###### Abstract

Unlike tabular data, features in network data are interconnected within a
domain-specific graph. Examples of this setting include gene expression
overlaid on a protein interaction network (PPI) and user opinions in a social
network. Network data is typically high-dimensional (large number of nodes)
and often contains outlier snapshot instances and noise. In addition, it is
often non-trivial and time-consuming to annotate instances with global labels
(e.g., disease or normal). How can we jointly select discriminative
subnetworks and representative instances for network data without supervision?

We address these challenges within an unsupervised framework for joint
subnetwork and instance selection in network data, called UISS, via a convex
self-representation objective. Given an unlabeled network dataset, UISS
identifies representative instances while ignoring outliers. It outperforms
state-of-the-art baselines on both discriminative subnetwork selection and
representative instance selection, achieving up to $10\%$ accuracy improvement
on all real-world data sets we use for evaluation. When employed for
exploratory analysis in RNA-seq network samples from multiple studies it
produces interpretable and informative summaries.

##  I Introduction

Network data abounds in a wide range of application domains: from activation
snapshots of the human brain to the global state of sensor and online social
networks. Different from tabular data, network data includes a shared
structure associating features (nodes) in addition to their values. This
structured feature space enables robust and interpretable solutions in a host
of tasks, such as subnetwork selection [11, 36], graph signal decomposition
[33], network-aware distance measures [3] and summarization [39, 32].

A network sample (i.e., instance) consists of feature values (node weights)
and a structure shared by all instances. Consider the example in Fig. 1 of
gene expression for $6$ patients ($P1-P6$), where the expression levels for
each gene ($G1-G7$) are employed as features and the human interactome is the
structure associating the features. Our goal is to select a subspace of well-
connected subnetworks which also “distinguish” the natural instance clusters.
One such feature set in the example includes genes $\\{G3,G4,G5\\}$ forming a
connected PPI subnetwork which also elucidates well-pronounced clusters of
patients: $\\{P1,P3,P4\\}$ and $\\{P2,P6\\}$. Another goal is to identify
outliers, such as patient $P5$ who does not “align” well with the natural
clusters.

Figure 1:  An illustrative example for unsupervised instance and subnetwork
learning from a set of instances—patients’ gene expression, sharing a common
network structure—the human protein-protein interaction network (PPI)
associating genes of known interactions. Our goal is to jointly select
representative instances and discriminative subnetworks which can then be used
for various downstream tasks.

It is important to note that network instance data is different from network
structure data employed for structural subgraph mining and classification
since the latter employs subgraph occurrences as features as opposed to node
values within a fixed structure [21]. Our setting is also different from
collective classification [37] typically arising in social media mining [9],
where individual nodes are instances as opposed to the whole network snapshot.

One important task in network instance data analysis is the prediction of the
global state (label) of network instances [11, 48] with applications in
neuroscience [5], gene expression on protein-protein interaction (PPI)
networks [15] and analysis of the state of transportation and sensor networks
[6]. Existing work typically focuses on the supervised setting, where each
sample is characterized by a global label (e.g., disease/healthy). To tackle
the typical high dimensionality existing work exploits the network structure
to identify connected discriminative sub-spaces (subgraphs and associated
feature values) to construct classifiers which generalize and offer
interpretability [11, 36].

While it may be easy to obtain a network instance, obtaining reliable global
labels for supervision may be challenging due to expensive human-in-the-loop
tests (e.g., doctor assessment for cognitive impairments [1]), or complex
global states (e.g., regimes of operation of infrastructure networks such as
road networks and power grids). In such scenarios unsupervised methods can
both elucidate important subnetworks and representative instances to be
labeled (akin to active learning [38]), while avoiding outlier and redundant
samples.

We propose UISS, an unsupervised approach to jointly learn discriminative
subgraphs and representative instances for network data. We formalize the
problem as a self-representative factorization of the data which also
maximizes the connectivity of selected nodes and the discriminative power of
selected instances. In particular, we select subnetworks and instances that
can reconstruct the full data (self-representation) and discount loss due to
outlier instances. We enforce network locality for selected nodes by a
quadratic penalty based on the network Laplacian. Finally, the representative
power of selected instances is promoted via an orthogonality constraint.

The main contributions of this paper are as follows:
1\. Significance: We address an unexplored problem of unsupervised learning
for network samples by jointly detecting discriminative subgraphs and
representative instances.

2\. Robustness: UISS is robust to noise and outliers, outperforming baselines
in feature and instance selection and achieving up to $10\%$ accuracy
improvement in real-world data.

3\. Interpretability and applicability: UISS discovers interpretable subspaces
such as discriminative city-center localities in bike rental behavior, and
reconstructs the developmental timeline in RNA-seq samples from organoid
development.

##  II Related Work

Unsupervised feature selection: Methods in this category employ filtering [20]
and embedding [46], however, the resulting features are not guaranteed to be
discriminative for the actual (withheld) instance classes. Learning
discriminative features based on cluster (pseudo) labels have been proposed
when actual instance labels are absent [28]. Other discriminative feature
selection methods [35, 27, 14] directly exploit cluster information: UDFS [46]
maximizes the cluster separation and compactness; RUFS [35] learns jointly
clusters and feature subsets by non-negative matrix factorization; NFDS [28]
imposes orthogonality constraints to avoid cluster ambiguity; and MCFS [7]
preserves across-cluster structure within a spectral embedding. Gu et al. [17]
developed feature selection method using the structure among instances. While
all above methods optimize the discriminative power of features, (i) they
cannot select subnetworks in network samples,(ii) they are not robust to
outliers which have been shown to impact unsupervised learning [29]; and (iii)
they expect the number of clusters as an input parameter. We overcome these
shortcomings by combining network smoothness and cluster-representative
instance selection via an orthogonality constraint and by self-representation.

Feature selection in networked settings has also been proposed for cases when
the instances are associated within a known network [25, 44, 26, 45]. However,
this setting is different (and complementary) to ours in that the network
structure associates instances as opposed to features.

Instance and feature co-selection: Feature and instance selection are co-
dependent tasks since representative instances highlight relevant features and
vice versa. CoSelect [41] and gActive [23] learn features and instances
effectively but require global label information and are applicable in
supervised settings only. Li et al. [24] developed an alternative unsupervised
model, called ALFS which, however, does not consider a network structure among
features or the existence of outliers, rendering it less robust than our
solution as we demonstrate empirically. Our settings also differ from
supervised active learning [38], which has been widely employed in text [42]
and network data [41]. Active learning in a fully unsupervised setting has
been pursued via transductive experimental design (TED) [47] which performs
instance selection guided by reconstruction error. Following a similar
minimum-reconstruction-error approach, ALFS [24] was proposed to jointly
select representative instances and features.

Learning with network data: Closest to our setting are recent supervised
methods for labeled network data [11, 10, 36, 48]. They adopt sampling or
optimization to learn subgraphs with high discriminative power for supervised
network state prediction. Unsupervised learning for network data, however,
presents a new challenge which has not been studied in the literature to the
best of our knowledge.

##  III Preliminaries

A network instance ${\mathcal{S}}_{i}=\left\\{V_{i},E_{i},{{\bf
X}_{i}}\right\\}$, is the triplet containing nodes $V_{i}$, edges $E_{i}$ and
node/feature values ${\bf X}_{i}$. A dataset ${\mathcal{DS}}=\\{S_{1}\dots
S_{n}\\}$ is a set of $n$ network instances. All instances share the same
network structure. Edges in ${\mathcal{S}}_{i}$ are weighted by the fraction
of their occurrences in instances ${\bf M}_{pq}=1/n\sum_{i}E_{i}(p,q)$. The
graph Laplacian matrix of ${\mathcal{S}}_{i}$ is defined as ${\bf L}={\bf
D}-{\bf M}$, where ${\bf D}$ is the diagonal matrix of weighted node degrees,
i.e., ${\bf D}_{pp}=\sum_{q}{\bf M}_{pq}$. The data matrix
$\mathbf{X}=\left[{X}_{1},{X}_{2},...{X}_{n}\right]^{T}\in\mathbb{R}^{n\times
m}$ of a dataset ${\mathcal{DS}}$ is comprised of network instance values
${X}_{i}$ in its rows.

##  IV Problem Formulation

Given a dataset ${\mathcal{DS}}$ of $n$ unlabeled network instances with $m$
nodes, one of our goal is to select $p~{}(p\ll n)$ instances which are
representative of the underlying clusters in the dataset and also to discover
a connected subgraph subspace of $q~{}(q\ll m)$ nodes. Intuitively, the above
design objective assumes that some instances form (an unknown number of)
clusters, while others are outliers and do not belong to clusters.

Self-representative factorization: We model the joint instance and subgraph
selection as a sparse self-representative factorization of the observations
$\mathbf{X}$. In traditional matrix factorization
$\mathbf{X}=\mathbf{U}\mathbf{V}$ [13], the two factors can be thought of as
latent indicators of row and column clusters. Such representations are
efficient and compact, but are also shown to exchange points from separate
low-dimensional structures in the data, leading to sub-optimal clustering
performance [31]. Instead, we propose a self-representative reconstruction
based on the principle that a matrix can be represented by factors selected
from its own rows and columns. In particular, we introduce selectors for
features ${\bf P}\in\mathbb{R}^{m\times k}$ and instances ${\bf
Q}\in\mathbb{R}^{k\times n}$, leading to:

$\displaystyle\underset{{\bf
P,Q}}{\mathrm{argmin}}\left\|\mathbf{X}-\mathbf{UV}\right\|_{2,1}+\lambda_{1}\left\|\mathbf{P}\right\|_{2,1}+\lambda_{2}\left\|\mathbf{Q}^{T}\right\|_{2,1},$
(1)
$\displaystyle~{}~{}~{}~{}s.t.~{}\mathbf{U}=\mathbf{XP},\mathbf{V}=\mathbf{QX},\mathbf{P}\geq
0,\mathbf{Q}\geq 0,$

where the $L_{2,1}$ norm is employed to ensure robustness to outliers and
noise through balance parameters $\lambda_{1}$ and $\lambda_{2}$. The key
distinction from existing latent factorization is that factors are selected
from the data columns and rows themselves. One can also view our co-selection
strategy as imposing a reduced-rank structure on a CUR-like decomposition [31]
via a bi-linear factorization employing $\bf P$ and $\bf Q$. Such a low-rank
structure promotes a compact and interpretable latent space of the instances
and subnetworks. Note that subnetwork and instance selection is based on
thresholding $\bf P$’s rows and $\bf Q$’s columns rather than a fixed number
of components pre-specified by the user. In particular, the importance of
subnetworks and instances can be quantified based on the corresponding row
norm of $\bf P$ and column norm of $\bf Q$.

Discriminative power using an orthogonal instance subspace: So far our
objective ensures instance and subnetwork selection, but does not ensure that
the selected instances are representative of the underlying clusters as
opposed to outliers. Instead of explicitly modeling cluster memberships and
deciding on the number of clusters apriori, we enforce orthogonality on the
selected instances via a constraint of the form
$\mathbf{V}\mathbf{V}^{T}=\mathbf{I}$. Columns of $\bf V$ act as cluster
indicators, but also as a dictionary basis to reconstruct all instances, and
thus this orthogonal constraint within the factorization ensures separation of
selected instances. While the important role of orthogonality in traditional
factorization clustering has long been recognized [13], we impose
orthogonality within a self-representative factorization, leading to a
superior performance in both feature and instance selection for network data.

Graph connectivity: We also impose “smoothness” of the subnetwork selection
with respect to the common graph structure $\mathcal{S}$ to encourage
connectivity in the feature space via a trace norm:
$\textup{Tr}(\mathbf{P}^{T}\mathbf{L}\mathbf{P})$, where $\mathbf{L}$ is the
graph Laplacian matrix of the nodes’ network.

Combining all modeling goals we obtain:

$\displaystyle\underset{{\bf
P,Q}}{\mathrm{argmin}}\left\|\mathbf{X}-\mathbf{UV}\right\|_{2,1}+\lambda_{1}\left\|\mathbf{P}\right\|_{2,1}+\lambda_{2}\left\|\mathbf{Q}^{T}\right\|_{2,1}+\lambda_{3}\textup{Tr}(\mathbf{P}^{T}\mathbf{L}\mathbf{P}),$
(2)
$\displaystyle~{}s.t.~{}\mathbf{U}=\mathbf{XP},\mathbf{V}=\mathbf{QX},\mathbf{V}\mathbf{V}^{T}=\mathbf{I},\mathbf{P},\mathbf{Q}\geq
0.\vspace*{-0cm}$

The objective unifies our design goals: joint (i) subnetwork selection via
$\mathbf{P}$, and (ii) discriminative instance selection via $\mathbf{Q}$,
where discriminative power is promoted by $\mathbf{V}$’s orthogonality.

##  V Optimization Algorithm

Our objective from Eq. 2 is jointly convex, however, it is non-trivial to
develop gradient-based solutions for it directly because the $L_{2,1}$ norm is
non-smooth. Instead, we propose an alternating optimization algorithm which
efficiently updates the instance selector $\bf P$ and feature selector $\bf Q$
in turn until convergence.

Update of Q: To update $\bf Q$, we fix $\bf P$ and obtain:

$\displaystyle\underset{{\bf V,Q}}{\mathrm{argmin}}\hskip
1.42271pt\left\|{{\bf X}}-\mathbf{U}{\bf V}\right\|_{2,1}+\left\|{\bf V}-{\bf
QX}\right\|_{F}^{2}+\lambda_{2}\left\|\mathbf{Q}^{T}\right\|_{2,1},$ (3)
$\displaystyle~{}~{}~{}s.t.~{}\mathbf{V}\mathbf{V}^{T}=\mathbf{I},\mathbf{Q}\geq
0\vspace*{-0cm}$

We set the gradient w.r.t $\bf Q$ of Eq. 3 to zero:

$\displaystyle({\bf QX}-{\bf V}){\bf X}^{T}+\lambda_{2}{\bf
Q}{\bf\Pi}=\mathbf{0},$ (4)

where $\bf\Pi$ is a diagonal matrix:
${\bf\Pi}_{ii}=(2\left\|\mathbf{Q}^{i}\right\|_{2})^{-1}$. Thus, we obtain a
closed-form solution for Eq. 4: ${\bf{Q}}={\bf V}{\bf X}^{T}({\bf X}{\bf
X}^{T}+\lambda_{2}{\bf\Pi})^{-1}$. Note that $({\bf X}{\bf
X}^{T}+\lambda_{2}{\bf\Pi})$ is invertible since ${\bf X}^{T}{\bf X}$ is
positive definite and ${\bf\Pi}$ is a non-negative diagonal matrix.

To handle the orthogonality constraint in Eq. 3, we introduce an intermediate
variable $\bf{W}$ which approximates $\bf{V}$, leading to two subproblems:

$~{}\begin{aligned}
\begin{cases}\underset{{\mathbf{V}}}{\mathrm{argmin}}\hskip
1.42271pt\left\|{{\bf X}}-\mathbf{U}{\bf V}\right\|_{2,1}+\left\|{\bf V}-{\bf
QX}\right\|_{F}^{2}+\left\|{\bf V}-{\bf W}\right\|_{F}^{2}\textup{(a)}\\\
\underset{{\bf W}}{\mathrm{argmin}}\hskip 1.42271pt\left\|{\bf V}-{\bf
W}\right\|_{F}^{2},\hskip 1.42271pts.t\hskip
2.84544pt\mathbf{W}\mathbf{W}^{T}=\mathbf{I}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\textup{(b)}\end{cases}\end{aligned}$
(5)

We set the gradient w.r.t. $\bf V$ in Eq. 5(a) to $\mathbf{0}$:

$\displaystyle{\bf U}^{T}({\bf U}{\bf V}-{{\bf X}}){\bf\Theta}+{\bf V}-{\bf
QX}+{\bf V}-{\bf W}=\mathbf{0},$ (6)

where ${\bf\Theta}$ is diagonal with elements
${\bf\Theta}_{ii}=\frac{1}{2\left\|{{\bf X}}^{i}-{\bf U}{\bf
V}^{i}\right\|_{2}}$ when ${{\bf X}}^{i}-{\bf U}{\bf V}^{i}\neq\mathbf{0}$,
and $\mathbf{0}$ otherwise. Since ${\bf U}^{T}{\bf U}$ is symmetric and
positive semi-definite, its eigendecomposition ${\bf U}^{T}{\bf U}={\bf
A}{\bf\Sigma}{\bf A}^{T}$ has real eigenvalues (diagonal of ${\bf\Sigma}$) and
real orthonormal eigenvectors in A. Employing this decomposition, we can
rewrite (6) as

$~{}\begin{aligned} {\bf A}{\bf\Sigma}{\bf A}^{T}{\bf V}{\bf\Theta}+2{\bf
V}={{{\bf\Psi}}},\end{aligned}$ (7)

where ${{\bf\Psi}}={\bf U}^{T}{{\bf X}}{{\bf\Theta}}+{\bf QX}+{\bf W}$. Since
$\bf{A}$ is orthonormal, we can multiply Eq. 7 by ${\bf A}^{T}$ on the left
side, which after substituting ${{{\bf E}}}={\bf A}^{T}{\bf V}$ can be
simplified to ${\bf\Sigma}{{{\bf E}}}{\bf\Theta}+2{{{\bf E}}}={\bf
A}^{T}{{{\bf\Psi}}}$. Individual elements of ${{\bf E}}$ can then be updated
as follows : ${{{\bf E}}}_{ij}={\left[{\bf
A}^{T}{{{\bf\Psi}}}\right]}_{ij}/({\bf\Sigma}_{ii}{\bf\Theta}_{jj}+2)$ and
$\bf V$ can be obtained as ${\bf V}={\bf A}{{{\bf E}}}$.

There exists a closed-form solution
$\mathbf{W}=\mathbf{R}\mathbf{I}\mathbf{T}^{T}$ solution for Eq. 5(b) due to
Viklands et al. [43], where $\mathbf{R}$ and $\mathbf{T}$ are the left and
right singular vectors of $\bf V$ in an SVD decomposition.

Update of P: The subproblem w.r.t $\bf P$ becomes:

$~{}\begin{aligned} &\underset{{\bf P}}{\mathrm{argmin}}\hskip
1.42271pt\left\|{{\bf X}}-\mathbf{X}{\bf
P}\mathbf{V}\right\|_{2,1}+\lambda_{1}\left\|\mathbf{P}\right\|_{2,1}+\lambda_{3}\textup{Tr}(\mathbf{P}^{T}\mathbf{L}\mathbf{P})\end{aligned}$

We again employ an auxiliary ${\bf U}$ to approximate $\mathbf{X}{\bf P}$:

$~{}\begin{aligned} &\underset{{\bf P,U}}{\mathrm{argmin}}\hskip
1.42271pt\left\|{{\bf
X}}-\mathbf{U}\mathbf{X}\right\|_{2,1}+\lambda_{1}\left\|\mathbf{P}\right\|_{2,1}+\lambda_{3}\textup{Tr}(\mathbf{P}^{T}\mathbf{L}\mathbf{P})\\\
&+\left\|\mathbf{U}-\mathbf{XP}\right\|_{F}^{2},\end{aligned}$ (8)

and alternate between $\bf P$ and $\bf U$ updates. We set the gradient w.r.t.
$\bf U$ in Eq. 8 to $\bf 0$:

$\displaystyle\mathbf{K}(\mathbf{U}\mathbf{V}-{{\bf
X}})\mathbf{V}^{T}+\mathbf{U}-\mathbf{XP}=0,$ (9)

where $\bf{K}$ is a diagonal matrix with elements ${\bf
K}_{ii}=\left(2\left\|{{\bf X}}_{i}-{\bf U}_{i}{\bf
V}\right\|_{2}\right)^{-1}$ when ${{\bf X}}_{i}-{\bf U}_{i}{\bf X}\neq 0$, and
${\bf K}_{ii}=0$ otherwise. Similar to the solution for $\bf{V}$, we employ
the eigendecomposition of ${\bf V}{\bf V}^{T}={\bf B}{\bf\Lambda}{\bf B}^{T}$
and multiply by ${\bf B}$ on the right, thus simplifying Eq. 9 to
$\mathbf{K}\mathbf{H}{\bf\Lambda}+{\bf H}={\bf\Xi}{\bf B}$, where ${\bf
H}={\bf UB}$ and ${\bf\Xi}=\mathbf{XP}+\mathbf{K}{{\bf X}}{\bf V}^{T}$. After
updating the elements of ${\bf H}$ according to ${\bf H}_{ij}=\frac{[{\bf\Xi
B}]_{ij}}{{\bf K}_{ii}{\bf\Lambda}_{jj}+1}$, we get a closed-form solution for
$\bf{U}$: ${\bf U}={\bf HB}^{T}$. We set the gradient w.r.t. $\bf{P}$ in Eq. 8
to $\bf 0$:

$~{}\begin{aligned}
{\bf{X}}^{T}\left(\mathbf{XP}-\mathbf{U}\right)+\lambda_{1}{\bf{G}}\mathbf{P}+\lambda_{3}\mathbf{L}\mathbf{P}=0,\end{aligned}$
(10)

where $\bf{G}$ is diagonal with elements ${\bf
G}_{ii}=(2\left\|\mathbf{P}_{i}\right\|_{2})^{-1}$ when
$\left\|\mathbf{P}_{i}\right\|_{2}\neq 0$, and ${\bf G}_{ii}=0$ otherwise.
Thus, $\bf P$ has the following closed-form solution:

$~{}\begin{aligned} {\bf
P}=(\mathbf{X}^{T}\mathbf{X}+\lambda_{1}\mathbf{G}+\lambda_{3}\mathbf{L})^{-1}{\bf{X}}^{T}\mathbf{U},\end{aligned}$
(11)

where $\mathbf{X}^{T}\mathbf{X}+\lambda_{1}\mathbf{G}+\lambda_{3}\mathbf{L}$
is invertible due to all summands being positive semi-definite matrices.

Algorithm 1 UISS

1:Training data $\bf{X}$, and
parameters($\lambda_{1},\lambda_{2},\lambda_{3},k$)

2:Selection matrices $\bf P,Q$

3:Initialize $\bf P,Q$, ${\bf\Theta}$, ${\bf\Pi}$ and ${\bf K}$ to identity
matrices; $\bf{U}$, $\bf V$, to random matrices;

4:while $\bf P,Q$ not converged do

5: $\left({\bf A},{\bf\Sigma}\right)=evd\left({\bf U}^{T}{\bf U}\right)$

6: ${\bf\Theta}_{ii}=\left(2\left\|{{\bf X}}^{i}-{\bf U}{\bf
V}^{i}\right\|_{2}\right)^{-1}$ if ${{\bf X}}^{i}-{\bf U}{\bf
V}^{i}\neq\mathbf{0}$

7: while $\bf{W}$ and $\bf{V}$ have not converged do

8: ${{\bf\Psi}}={\bf U}^{T}{{\bf X}}{{\bf\Theta}}+{\bf QX}+{\bf W}$

9: ${{{\bf E}}}_{ij}={\left[{\bf
A}^{T}{{{\bf\Psi}}}\right]}_{ij}/({\bf\Sigma}_{ii}{\bf\Theta}_{jj}+2)$

10: ${\bf V}={\bf A}{\bf E}$

11: $({\mathbf{R}},{\mathbf{T}})=svd\left(\mathbf{V}\right)$

12: ${\bf W}={\bf R}{\bf I}{\bf T}^{T}$

13: ${\bf{Q}}=max[{\bf V}{\bf X}^{T}({\bf X}{\bf
X}^{T}+\lambda_{2}{\bf\Pi})^{-1},0]$

14: ${\bf\Pi}_{ii}=(2\left\|\mathbf{Q}^{i}\right\|_{2})^{-1}$ if ${{\bf
Q}}_{i}\neq 0$

15: $\left({\bf B},{\bf\Lambda}\right)=evd\left({\bf V}{\bf V}^{T}\right)$

16: ${\bf
P}=max[(\mathbf{X}^{T}\mathbf{X}+\lambda_{1}\mathbf{G}+\lambda_{3}\mathbf{L})^{-1}{\bf{X}}^{T}\mathbf{U},0]$

17: ${\bf\Xi}=\mathbf{XP}+\mathbf{K}{{\bf X}}{\bf V}^{T}$ ${\bf
H}_{ij}=[{\bf\Xi B}]_{ij}/({\bf K}_{ii}{\bf\Lambda}_{jj}+1)$

18: ${\bf U}={\bf H}{\bf B}^{T}$

19: ${\bf K}_{ii}=\left(2\left\|{{\bf X}}_{i}-{\bf U}_{i}{\bf
X}\right\|_{2}\right)^{-1}$ if ${{\bf X}}_{i}-{\bf U}_{i}{\bf X}\neq 0$

20:return $\bf P,Q$

Overall algorithm, complexity and convergence: The detailed steps of the UISS
are presented in Algorithm 1 and follow our update derivations above. We
considered both iteration to convergence for each of the sub-problems
independently and iterations of one update for each of the variables. The
former strategy results in faster convergence for $\bf{W}$ and $\bf{V}$, hence
this is the one we present in the Alg. 1. It is important to note that all
updates to diagonal matrices (Steps 5, 15, 21) are applied after they are re-
initialized to all-zeroes at each iteration.

The complexity of our algorithm is determined by the number of iterations to
convergence (in both loops) and the complexity of individual steps in the
loops. The most costly steps involve the SVD operation in (Step 10) and the
matrix inversions (Steps 13 and 17). While $svd$ has a super-quadratic
complexity in the worst case, efficient implementations enable good
scalability in practice. The matrix inversions can also be implemented
efficiently due to the low-rank updates in each iteration. Note that in Step
13 ${\bf XX}^{T}$ and in Step 17 ${\bf X}^{T}{\bf X}+\lambda_{3}{\bf L}$ are
constants and can be inverted once and re-used in subsequent iterations, thus
exploiting the sparse structure via the Sherman-Morrison formula for sparse
inverse updates:

$(\mathbf{A}+\gamma\mathbf{D})^{-1}=\mathbf{A}^{-1}-\frac{\mathbf{A}^{-1}\mathbf{d}\mathbf{d}^{T}\mathbf{A}^{-1}}{1+\mathbf{d}^{T}\mathbf{A}^{-1}\mathbf{d}},\vspace*{-0cm}$
(12)

where $A$ corresponds to the constant matrices in the two steps,
$\mathbf{d}=\sqrt{\gamma diag(\mathbf{D})}$ is a column vector and the square
root is applied element-wise thus $\gamma\mathbf{D}=d^{T}d$. Therefore, we can
reduce the complexity of inversion as $O\left(mn_{N}\right)$, where $n_{N}$ is
the number of non-zero elements in $\bf X$. The pair $({\bf D},\gamma)$ from
the above formula corresponds to $({\bf\Pi},\lambda_{2})$ in Step 12 and
$({\bf G},\lambda_{1})$ in Step 16 respectively. _We will publish the
implementation of UISS at the camera-ready paper._

All subproblems in the Alg. 1 are solved in an alternating manner and have
closed-form solutions, where $\left\\{{\bf P,Q},{\bf U},{\bf V}\right\\}$ are
optimal solution of Eq. 2. These closed-form solutions guarantee the overall
optimization sequence converges to the primal-dual optimal solutions.
Meanwhile, the constraints in Alg. 1 will converge to zero, i.e.
$\left\|\mathbf{P}^{iter+1}-\mathbf{P}^{iter}\right\|_{F}\rightarrow
0,\left\|\mathbf{Q}^{iter+1}-\mathbf{Q}^{iter}\right\|_{F}\rightarrow 0$.
Therefore, we can obtain the global convergence of UISS because it features
both sequence convergence and constraint convergence [19].

Dataset | $|\mathcal{V}|$ | $|E|$ | $|\mathcal{DS}|$ | Hidden Classes
---|---|---|---|---
Synthetic | 100 | 637 | 400 | positive/negative
Bike [2] | 142 | 1,723 | 299 | weekday/weekend
Embryo [15] | 1,321 | 5,227 | 34 | tissue layers
CCT [8] | 4,665 | 270,571 | 184 | weekday/weekend
ADNI [1] | 6,216 | 683,760 | 173 | AD/NC
Liver [22] | 7,383 | 251,916 | 123 | disease/control
TABLE I: Summary of evaluation datasets’ statistics.

##  VI Experimental Evaluation

###  VI-A Datasets.

We employ both synthetic and real-world datasets for evaluation and summarize
their statistics in Table. I. We synthesize geometric networks by uniformly
sampling nodes in a unit square and connecting nodes at distances smaller than
a threshold of $\tau=0.2$. We select well-connected ground truth subgraphs and
generate balanced set of instances labeled by two (hidden) global states.
Particularly, nodes from the ground truth are assigned a value between
$[50,100]$ for positive instances and $[-100,-50]$ for negative ones. Then we
add Gaussian noise to all nodes, where ground truth and non-ground-truth nodes
have different mean values set to $10$ and $70$ respectively. We also inject
outlier instances in the data for some of our experiments.

Nodes in the _Bike_ [2] dataset are bicycle rental stations in Boston and
edges connect stations based on a distance threshold. Nodes’ values correspond
to the number of check-outs in a day, where we employ the last $299$ days of
the data. We also employ two gene expression datasets: _Liver_ metastasis [22]
and Embryo [15]. Their network structures are PPI networks [12], while node
features correspond to gene expression values with hidden global labels:
healthy/normal subjects in Liver and tissue type in Embryo. The _ADNI_ [1]
dataset contains fMRI resting state measurements for subjects labeled by AD:
suffering Alzheimer’s disease and NC: healthy normal controls. The graph
structure associates functional links (nodes) with their level of coherence
(feature values). Nodes are connected if the corresponding functional links
share a brain region. _CCT_ contains city cellular HTTP traffic data records
[8] where nodes are stations, hourly requests are feature values, and node
pairs are connected based on a distance threshold. Hidden global labels
reflect if the snapshot occurred during workday hours (8am-16pm) or off hours.

Methods | FS | IS | Netw. | Description
---|---|---|---|---
ALFS [24] | ✓ | ✓ |  | Optimal reconstruction
DSLw/o [48] | ✓ |  | ✓ | Sparse self-representation
UFSOL [18] | ✓ |  |  | Locality preservation
RRSS [34] |  | ✓ |  | Subspace learning
UISS | ✓ | ✓ | ✓ | Our method
TABLE II: Summary of competing techniques. FS: feature selection; IS: instance
selection; Network: employs the inter-feature network.

###  VI-B Experimental setup.

Baselines: All competing techniques are summarized in Tbl. II. We employ
several recent feature selection methods to conduct comparison experiments,
ALFS [24], employs data reconstruction to find the instances and features that
minimize the reconstruction error; DSLw/o is a variant of the state-of-art
supervised discriminative subgraph learning method [48] which we restrict to
the unsupervised setting by removing the SVM-like supervision from its
objective. UFSOL [18] selects features by preserving the local structure in
the data. We also employ two recent unsupervised instance selectors: ALFS [24]
which selects instances and features for minimum-reconstruction error and RRSS
[34] which selects instance based on sparse subspace learning.

Metrics: Following the typical experimental setup in unsupervised feature [18]
and instance [24] selection, we first perform these tasks in an unsupervised
manner, then reveal the hidden class labels and measure the accuracy of the
selected features and/or instances. To allow for fair comparison, we employ
the same classifier SVM (linear kernel, $C=1$) for all methods’ selections and
report average accuracy of $50$ runs.For datasets with ground-truth features
(Liver and Synthetic) and synthetically injected outliers, we also measure the
ROC of recovering these ground truth elements respectively.

###  VI-C Experiments on synthetic network data.

Subnetwork selection: We first evaluate the ability of competing techniques to
detect injected ground truth (GT) feature subgraphs in synthetic data. The AUC
for ground truth subgraph detection with increasing noise variance is
presented in Fig. LABEL:fig:syn_noise_fs_roc. While all methods’ performance
decreases with increasing noise variance, UISS’s ability to recover the
injected GT subgraphs consistently dominates that of baselines. Though DSLw/o
uses the network structure for subnetwork selection, similar to UISS, it does
not optimize unsupervised discriminative power and is sensitive to noise and
outliers. The comparatively lower accuracy of UFSOL is due to the low
robustness to noise of its cluster-based feature selection strategy. Similarly
ALFS deteriorates quickly with the noise variance due to its preference for
high-variance features dominating its Frobenius reconstruction loss. UISS is
less sensitive to noise due to the $L_{2,1}$ reconstruction loss, the use of
the network structure among features and the imposed representative power via
orthogonality.

(a) Bike

(b) Embryo

(c) ADNI

(d) CCT

Figure 2: Comparison of classification accuracy on real world datasets with
increasing number of selected features.

We also study the robustness of subnetwork selection for increasing the number
of outlier instances and provide comparison in Fig.
LABEL:fig:syn_outlier_fs_roc. UISS dominates and degrades more gracefully with
the number of outlier compared to alternatives. This robustness can be
attributed to the combination of orthogonality of selected instances for self-
representation and the $L_{2,1}$ reconstruction cost, collectively discarding
outliers as unable to represent other instances and incurring high
reconstruction error. All baseline methods’ AUC suffers at all levels of
outliers since they treats all instances as cluster members and do not
explicitly model the existence of outliers.

Instance selection: Next we evaluate the utility of competing techniques as
discriminative instance selectors (Fig. LABEL:fig:syn_noise_is_acc). _Can
selected instances inform accurate classifiers if global labels are added
post-selection?_ Among the instance selector baselines, ALFS can also perform
feature selection while RRSS selects representative instances based on all
features. Note that while ALFS and UISS internally re-weigh features top
select instances, the final classifiers we train in this experiment are based
on all features in selected instances to enable a fair comparison with RRSS.

When all training instances are available, the accuracy of competitors is the
same as they employ the same SVM classifier and all features. As we constrain
the methods to select smaller number of instances for labeling, the baselines’
performance degrades, while that of UISS increases slightly and remains above
the all-instance accuracy even when we select only $10\%$ of the total
instances (Fig. LABEL:fig:syn_noise_is_acc). The advantage of our model stems
from the dependence of the instance selector $\bf Q$ on the feature selector
$\bf P$ which enforces smoothness in the network structure $\mathcal{S}$ as
well as the imposed orthogonality on the selected representative instances. In
other words, while we do not explicitly exclude features, UISS treats the
available features deferentially which in turn informs better instance
selection. This advantage is to some extent noticeable for ALFS as well,
although its inability to consider the network structure among features leads
to significantly degraded accuracy for decreasing number of selected
instances.

Outlier detection. We also employ instance selectors to detect outlier
instances injected in the dataset (Fig. LABEL:fig:syn_outlier_is_acc). We
invert the instance selection order reported by each method and plot the ROC
for predicting outlier instances (the dataset is balanced, i.e. the same
number of outliers and non-outlier instances). UISS is better at detecting
outliers than ALFS, since the latter employ Frobenius reconstruction
strategies, thus forcing outliers’ inclusion among important instances.
Although RRSS employs an $L_{2,1}$-norm fitting function, it does not to
consider the structure among features when informing the instance selection.
Our choice of $L_{2,1}$ norm for reconstruction and the orthogonality
regularization on the instance association matrix allows UISS to “spot”
outlier instances better.

(a) Classification accuracy

(b) Ground truth features

Figure 3: Classification accuracy (a) and ROC for selecting ground truth genes
(b) in the Liver dataset.

(a) Bike

(b) Liver

(c) ADNI

(d) CCT

Figure 4: Comparison of classification accuracy with increasing number of
selected instances on real world datasets.

(a) Outliers (Bike)

(b) Outliers (CCT)

(c) Joint FS and IS

(d) UISS

(e) UFSOL

Figure 5: (5(a)),(5(b)): Outlier detection on real-world datasets. (5(c))
Accuracy comparison of joint unsupervised feature and instance selection
between UISS (left) and ALFS (right) on the CCT dataset. The number of
selected features varies on the vertical axis. UISS achieves over $82\%$
accuracy with as few as $30$ features and $30$ instances, while ALFS requires
all instances to attain the same performances. (5(d)),(5(e)): Visualization of
the subgraphs selected by UISS and UFSOL on the Boston Bike dataset.

###  VI-D Experiments on real-world network data

Subnetwork selection:  Next we vary the number of features in the selected
subnetwork and then add the hidden labels to quantify the cross-validation
accuracy (based on selected features in all instances) in the network datasets
(Fig. 2). UISS’s accuracy dominates all baselines on all datasets and its
accuracy is relatively stable around its optimal value. A second pattern that
stands out is that there is no clear second best baseline. While DSLw/o
dominated other baselines in identifying GT features in synthetic (Figs.
LABEL:fig:syn_noise_fs_roc), its ability to select predictive features on real
datasets does not consistently dominate non-network baselines: better on CCT,
ADNI and Bike, but worse on Embryo where features are binary. The best
improvement of UISS over baselines in these datasets is about $10\%$ of
accuracy and it is retained over different number of selected features.

In addition, we quantify the subgraph selection quality for the Liver dataset
that has partial ground truth genes associated with the global disease state
labels and reported in [11] (Figs. 3). We observe that the methods which
consider the PPI network structure among features (genes), namely UISS and
DSLw/o, are better at recovering known genes associated with the disease than
alternative which do not consider this structure. In particular the AUC of
UISS and DSLw/o are $0.79$ and $0.65$ respectively, while those of non-network
alternatives are close to random (Fig. 3(b)). Without considering the network
structure, UFSOL and ALFS tend to select isolated nodes rather than a
connected subgraph and those are not part of potential target pathways which
form connected subgraphs in the PPI network. UFSOL achieves the second best
performance after UISS in terms of classification accuracy (3(a)), however,
its AUC for ground truth gene detection is the worst among competitors (Fig.
3(b)). It is important to note that the set of ground truth genes as discussed
in [11] is far from complete, i.e. there are potentially more genes associated
with the disease, and thus the competing techniques could be employed to
detect more target genes and pathways of biological significance.

Instance selection: Similar to synthetic data, we also compare instance
selectors on classification accuracy for increasing number of informative
instances selected in an unsupervised manner (Fig. 4). The observed behavior
in real-world dataset is similar to that in synthetic data. Namely, UISS
dominates alternatives for small number of selected instances and it achieves
its optimal accuracy without using all available instances. This is a very
promising result for applications in which acquiring annotations is expensive
or time consuming, since employing UISS would enable the creation of accurate
classifiers with minimal number of instances.

Outlier detection: We also evaluate the ability of competing techniques to
detect synthetically injected outliers in real-world datasets in Fig. 5. As in
synthetic, we invert the scores for instances to include, i.e. the lowest-
ranking instances are deemed highest-scoring outliers. Injected outlier
instances have random feature values with similar mean value to instances in
the respective real-world dataset. UISS achieves close-to-optimal performance
on the Bike and CCT datasets, while ALFS behaves close to random due to its
sensitivity to outliers. RRSS is the second best method on the Bike dataset as
it explicitly models outliers for instance selection. The behavior of
competitors on CCT is also close to random, although in this dataset their
initial FPR growth is steeper and later flattens due to ranking some outliers
among the most important instances.

###  VI-E Joint feature and instance selection.

To enable a fair comparison with non-network baselines in our feature and
instance selection experiments from the previous subsections, we sub-selected
one of those dimensions and used the other one fully. The performance of UISS
improves both when some number of features and a subset of all available
instances are excluded. The optimal subsets in both dimensions are likely
inter-independent, and hence, in what follows, we study this dependence by an
exhaustive sub-selection of both instances and features in the CCT dataset.
Note that this dataset is very high dimensional and learning reliable
predictive models on it is likely to benefit from both feature and instance
selection. The only baseline that performs feature and instance selection
jointly is ALFS, and hence, we focus on comparison between UISS and ALFS in
this experiment.

We vary the number of selected features from $10\%$ to $100\%$ (top to bottom
in Fig. 5(c)) and the number of instances between $5\%$ and $50\%$ (left to
right in Fig. 5(c)) in the CCT dataset and present the cross-validation
accuracy after revealing the hidden labels for the selected feature subsets by
the methods (Fig. 5(c)). UISS achieves its optimal accuracy ($81\%$) with as
few as $30\%$ of the features and $30\%$ of the instances. To get to a similar
performance ALFS employs $100\%$ of the instances and $20\%$ of the features,
i.e. three times the number of instances and $10\%$ less features. Overall,
both methods benefit from the joint selection of features and instances,
however, UISS utilizes the network structure and further employs
regularization making it more robust to outliers.

###  VI-F Performance on non-network data.

(a) Feature selection

(b) Instance selection

Figure 6: Feature and instance selection accuracy on a non-network dataset
USPS [4].

UISS is designed specifically as an unsupervised feature and instance selector
for network data. However, one important question about its performance is
whether its advantage is only due to the availability of network structure
among the data features. To test this, we turn to a common image (“non-
network”) dataset USPS [4]. This dataset was previously adopted for both
feature and instance selection by some of our baselines and also more broadly
in the machine learning literature. We compare a non-network version of UISS
($\lambda_{3}=0$) to all baselines for feature and instance selection cross-
validation accuracy (Fig. 6), following the same protocol as the one adopted
for network data. Interestingly, when employed for instance selection, UISS
once again dominates all baselines when limiting the number of instances
employed to train a classifier (Fig. 6(b)). The feature selection experiment
renders our method still remains competitive with state-of-art feature
selection methods when they are restricted to use between $1\%$ and $10\%$ of
the available features (Fig. 6(a)). While these experiments show promise for
the generality of UISS beyond network data, investigation of more non-network
datasets is necessary to confirm its utility in such settings. Extensive non-
network experiments are beyond the scope and space constraints of the current
manuscript and we plan to include such analysis in an extended journal
version.

###  VI-G UISS at work: salivary gland organoids and bike sharing.

Organoids: We next evaluate the ability of UISS to elucidate the organization
of RNA sequencing (RNA-seq) networked samples combined from two different
studies of mouse salivary gland development [40]. RNA-seq quantifies the
levels of RNAs present in a tissue sample and the Tanaka study [40] employed
RNA-seq to characterize the gland organoids derived from embryonic stem cells
which are grown ex vivo and engineered to mimic in vivo organ development.
Organoids offer rapid disease modeling with applications to infectious
diseases from Zika to SARS-CoV-V2.This is a fitting use case for UISS, since
while fusing multiple unlabeled, high-dimensional and noisy datasets may
enable important new insights, it also warrants examination for outliers as
well as artifacts of varying experimental protocols. The organoid data set
includes 38 instances: 19 time points (some with replicates), spanning from
embryonic day 12 (E12) to 12 weeks into adulthood (P84), as well as different
cell types (GE/GM/OE) and manipulations (iSG/T-iSG) labeled in Fig. 7. We
employed the fold changes in RNA-seq counts as input features after a
regularized logarithm transformation [30].

Figure 7: UISS employed to analyze networked RNA-seq samples from salivary
gland development. Samples are annotated post-analysis with (i) development
time (E12-P84), and (ii) sample types GM/GE/EB/OE/(T)iSG. The embedding
learned by UISS reconstructs the developmental timeline and isolates expected
stem cell outliers (EB) and differentiates between organoid (iSG) and in vivo
development. ($k=20,\lambda_{1}=0.1,\lambda_{2}=0.1,\lambda_{3}=10^{3}$).

We apply UISS on the data and obtain a 3D embedding of the instances by
employing PCA on the learned data projection $\bf XP$ presented in Fig. 7. It
is important to note that the labels in the figure were not used by UISS, i.e.
the analysis is fully unsupervised. Embryonic stem cell derived embryoid
bodies (EBs) are pluripotent and able to become all known cell types. Not
surprisingly, UISS renders them based on their RNA-seq as outliers compared to
other samples. After the EBs are induced to differentiate into salivary glands
(iSG) in vitro their transcriptional landscape changes to become more similar
to developing glands. Using the time course labels (E12 through P84), we can
place experimental data in frame with the known developmental (in vivo)
progression (black dotted line). Organoids (iSG) are transcriptionally similar
to an early salivary gland differentiation stage around E15-16, however, once
upon transplantation in vivo (samples T-iSG) the organoid overcome a
developmental wall furthering their maturation. The T-iSG transcriptomes
should be comparable to postnatal salivary glands from five days after birth
(P5) to 84 days after birth (P84) from the Gluck study [16].

Interestingly, we detect a wider range of variability in maturation of the
transplanted samples, which could represent either variability in the
transplanted organoids or differential responses of the individual organoids
to the in vivo environment. Although the representation by UISS cannot readily
explain the variable gene expression by the organoids, this example of staging
organoids illustrates its ability to map relationships between biological
datasets. As RNA-seq is broadly used in biomedical sciences, UISS will have
broader applicability in comparison of biological datasets. It is important to
note that, PCA on the UISS’ embedding of the data $\bf XP$ retains over $99\%$
of the variance in the first component, while PCA on the raw data $X$ requires
more than $12$ components for $99\%$ variance retention.

Bike sharing: We also investigate the interpretability of the feature
selection by UISS in the Bike dataset. Figs. 5(d), 5(e) visualize the selected
subgraphs (highlighted nodes and edges) by UISS and UFSOL—its closest
competitor in terms of accuracy when employing at least $6$ features (see Fig.
2(a)). As expected, UISS selects a mostly connected feature subgraph due to
the network smoothness regularization. This subgraph corresponds to a locality
of bike rental stations in downtown Boston which can differentiate between
weekday and weekend traffic (hidden classes in this dataset). UFSOL selects
features which form multiple connected components which are less obvious to
interpret and more importantly result in less predictive classifiers (as
demonstrated in Fig. 2(a)).

###  VI-H Parameter sensitivity and running time

We evaluate UISS’s scalability and sensitivity to various parameters. The
experiment shows that our method’s running time grows nearly linearly for both
increasing number of instances and nodes. Additionally, we also investigate
the sensitivity of UISS to its three hyper-parameters, i.e. $\lambda_{1}$,
$\lambda_{2}$ and $\lambda_{3}$. All combinations of hyperparamters exhibit
similar robust behavior. Please find more details at the supplement material.

We also investigate the sensitivity of UISS to its three hyper-parameters:
$\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$, which control the number of
selected subnetwork, instances, and level of graphs smoothness in feature
selection respectively. These parameters all have a physical meaning and can
provide useful control of the end user who want to enforce each of the
corresponding behaviors based expert knowledge about the domain in which she
employs UISS. However, in some cases the optimal parameter setting may be hard
to obtain apriori, a typical challenge with many unsupervised methods. Hence
we analyze the accuracy stability under small variations of these parameters.

Based on our analysis presented in Fig. 8 the quality of UISS does not vary
significantly with the parameter settings. For this analysis, we measure the
performance in terms of feature selection accuracy by fixing one parameter
($\lambda_{3}$) and updating the other two. The overall performance on
different combinations is stable and close to peak accuracy over large
parameter ranges. The performance drops when at one of the parameter becomes
significantly smaller (order of magnitude) than the other two. The reason for
this behavior is that the three regularization terms in our objective have
similar contribution to the cost, and thus keeping their importance balanced
results in optimal performance. In addition, small values of $\lambda_{1}$ and
$\lambda_{2}$ will not enforce sufficient sparsity for the representative
instance and feature selection and lead to poor performance due to negative
effect of outliers and noise. Other combinations of hyperparamters exhibit
similar robust behavior (in supplement material). Also, by default we set
$k=10$, but our analysis does not render UISS sensitive to third internal
dimension of the self-representative factorization.

(a) Running time

(b) $\lambda_{1}$ v.s. $\lambda_{2}$

(c) $\lambda_{1}$ v.s. $\lambda_{3}$

(d) $\lambda_{2}$ v.s. $\lambda_{3}$

Figure 8: Parameter sensitivity and running time on the synthetic dataset.

##  VII Conclusion

In this paper we proposed and evaluated UISS: a general unsupervised approach
for joint subnetwork and instance selection in network data. It performs
interpretable subnetwork and instance selections via a self-representative
factorization which enforces smoothness on the inter-feature network and
robustness to outlier instances and noise. UISS dominated alternatives in both
instance and subnetwork selection, often achieving its highest accuracy using
significantly fewer features and instances than those employed by
corresponding baselines. Meanwhile, our method is able to discover
interpretable subnetwork in the data such as city center localities of bike
rental behavior distinguishing between workdays and weekends.

## Acknowledgment

The work is supported in part by the NSF Smart and Connected Communities
(SC&C) CMMI grant 1831547 and in part by the National Institutes of Health
through NIH grant R01DE027953.

## References

  * [1] ADNI project.  http://www.adni-info.org/.
  * [2] Hubway data visualization challenge:.  http://hubwaydatachallenge.org.
  * [3] V. Amelkin, P. Bogdanov, and A. K. Singh.  A distance measure for the analysis of polar opinion dynamics in social networks.  In Proc. of the Intl. Conference on Data Engineering (ICDE). IEEE, 2017.
  * [4] P. N. Belhumeur, J. P. Hespanha, and D. J. Kriegman.  Eigenfaces vs. fisherfaces: Recognition using class specific linear projection.  1997\.
  * [5] P. Bogdanov, N. Dereli, X.-H. Dang, D. S. Bassett, N. F. Wymbs, S. T. Grafton, and A. K. Singh.  Learning about learning: Mining human brain sub-network biomarkers from fmri data.  PLoS One, 2017.
  * [6] P. Bogdanov, M. Mongiovi, and A. K. Singh.  Mining heavy subgraphs in time-evolving networks.  In Proceedings of the IEEE International Conference on Data Mining (ICDM), 2011.
  * [7] D. Cai, C. Zhang, and X. He.  Unsupervised feature selection for multi-cluster data.  In Proc. of ACM SIGKDD Intl. Conference on Knowledge Discovery and Data Mining, KDD ’10, 2010.
  * [8] X. Chen, Y. Jin, S. Qiang, W. Hu, and K. Jiang.  Analyzing and modeling spatio-temporal dependence of cellular traffic at city scale.  In 2015 IEEE international conference on communications (ICC), pages 3585–3591. IEEE, 2015.
  * [9] K. Cheng, J. Li, and H. Liu.  Unsupervised feature selection in signed social networks.  In Proc. of the SIGKDD Intl. Conference on Knowledge Discovery and Data Mining, KDD ’17, 2017.
  * [10] X. H. Dang, A. K. Singh, et al.  Discriminative subnetworks with regularized spectral learning for global-state network data.  In Proc. of the European Conference on Machine Learning and Knowledge Discovery in Databases - Volume 8724, ECML PKDD 2014, pages 290–306, 2014.
  * [11] X.-H. Dang, H. You, P. Bogdanov, and A. Singh.  Learning predictive substructures with regularization for network data.  In Proceedings of the IEEE International Conference on Data Mining (ICDM), 2015.
  * [12] R. Dannenfelser, N. R. Clark, and A. Ma’ayan.  Genes2fans: connecting genes through functional association networks.  BMC Bioinformatics, 13(1):156, Jul 2012.
  * [13] C. Ding, T. Li, W. Peng, and H. Park.  Orthogonal nonnegative matrix t-factorizations for clustering.  In Proc. of the 12th ACM SIGKDD Intl. Conference on Knowledge Discovery and Data Mining, KDD ’06, pages 126–135, New York, NY, USA, 2006. ACM.
  * [14] L. Du and Y.-D. Shen.  Unsupervised feature selection with adaptive structure learning.  In Proc. of the 21th ACM SIGKDD Intl. Conference on Knowledge Discovery and Data Mining, KDD ’15, pages 209–218, New York, NY, USA, 2015. ACM.
  * [15] J. Dutkowski and T. Ideker.  Protein networks as logic functions in development and cancer.  PLoS comp. biology, 7(9):e1002180, 2011.
  * [16] C. Gluck, S. Min, et al.  Rna-seq based transcriptomic map reveals new insights into mouse salivary gland development and maturation.  BMC genomics, 17(1):923, 2016.
  * [17] Q. Gu, Z. Li, and J. Han.  Joint feature selection and subspace learning.  In Proc. of the Intl. Joint Conf. on Artificial Intelligence, IJCAI’11, 2011.
  * [18] J. Guo, Y. Quo, X. Kong, and R. He.  Unsupervised feature selection with ordinal locality.  In IEEE Intl. Conference on Multimedia and Expo (ICME), pages 1213–1218, July 2017.
  * [19] B.-S. He, L.-Z. Liao, and H. Yang.  A new inexact alternating directions method for monotone variational inequalities.  Mathematical Programming, 92:103–118, 01 2002.
  * [20] X. He, D. Cai, and P. Niyogi.  Laplacian score for feature selection.  In Proc. of the Intl. Conference on Neural Information Processing Systems, NIPS’05, pages 507–514, Cambridge, MA, USA, 2005. MIT Press.
  * [21] C. Jiang, F. Coenen, and M. Zito.  A survey of frequent subgraph mining algorithms.  The Knowledge Engineering Review, 28(1):75–105, 2013.
  * [22] D. H. Ki, H.-C. Jeung, C. H. Park, S. H. Kang, G. Y. Lee, W. S. Lee, N. K. Kim, H. C. Chung, and S. Y. Rha.  Whole genome analysis for liver metastasis gene signatures in colorectal cancer.  International journal of cancer, 121(9):2005–2012, 2007.
  * [23] X. Kong, W. Fan, and P. S. Yu.  Dual active feature and sample selection for graph classification.  In Proceedings of the 17th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 654–662, 2011.
  * [24] C. Li, X. Wang, W. Dong, J. Yan, Q. Liu, and H. Zha.  Joint active learning with feature selection via cur matrix decomposition.  IEEE Trans. on Pattern Analysis and Machine Intelligence, 2018.
  * [25] J. Li, R. Guo, C. Liu, and H. Liu.  Adaptive unsupervised feature selection on attributed networks.  In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 92–100, 2019.
  * [26] J. Li, X. Hu, J. Tang, and H. Liu.  Unsupervised streaming feature selection in social media.  In Proc. of the 24th ACM Intl. on Conference on Information and Knowledge Management, CIKM ’15, pages 1041–1050, New York, NY, USA, 2015. ACM.
  * [27] J. Li, X. Hu, L. Wu, and H. Liu.  Robust unsupervised feature selection on networked data.  In Proceedings of the 2016 SIAM International Conference on Data Mining, pages 387–395. SIAM, 2016.
  * [28] Z. Li, Y. Yang, J. Liu, X. Zhou, and H. Lu.  Unsupervised feature selection using nonnegative spectral analysis.  In Proc. of the AAAI Conference on Artificial Intelligence, 2012\.
  * [29] H. Liu and S. Yan.  Robust graph mode seeking by graph shift.  In Proc. of the 27th Intl. Conference on Intl. Conference on Machine Learning, ICML’10, pages 671–678, USA, 2010. Omnipress.
  * [30] M. I. Love, W. Huber, and S. Anders.  Moderated estimation of fold change and dispersion for rna-seq data with deseq2.  Genome biology, 15(12):550, 2014.
  * [31] M. W. Mahoney and P. Drineas.  Cur matrix decompositions for improved data analysis.  Proceedings of the National Academy of Sciences, 106(3):697–702, 2009.
  * [32] M. McNeil, B. Ma, and P. Bogdanov.  SAGA: Signal-aware graph aggregation.  In Proceedings of SIAM International Conference on Data Mining (SDM), 2022.
  * [33] M. J. McNeil, L. Zhang, and P. Bogdanov.  Temporal graph signal decomposition.  In Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining, pages 1191–1201, 2021.
  * [34] F. Nie, H. Wang, H. Huang, and C. Ding.  Early active learning via robust representation and structured sparsity.  In Proc. of the Intl. Joint Conference on Artificial Intelligence, IJCAI ’13, 2013.
  * [35] M. Qian and C. Zhai.  Robust unsupervised feature selection.  In Proc. of the Intl. Joint Conf. on Artificial Intelligence, IJCAI ’13, 2013.
  * [36] S. Ranu, M. Hoang, and A. Singh.  Mining discriminative subgraphs from global-state networks.  In Proc. of the 19th ACM SIGKDD intl. conference on Knowledge discovery and data mining, pages 509–517. ACM, 2013.
  * [37] P. Sen, G. Namata, M. Bilgic, L. Getoor, B. Galligher, and T. Eliassi-Rad.  Collective classification in network data.  AI magazine, 29(3):93–93, 2008.
  * [38] B. Settles.  Active learning literature survey.  Technical report, 2010.
  * [39] A. Silva, P. Bogdanov, and A. Singh.  Hierarchical in-network attribute compression via importance sampling.  In Proc. of the 31st IEEE Intl. Conference on Data Engineering (ICDE), 2015.
  * [40] J. Tanaka, M. Ogawa, H. Hojo, Y. Kawashima, Y. Mabuchi, K. Hata, S. Nakamura, R. Yasuhara, K. Takamatsu, T. Irié, et al.  Generation of orthotopically functional salivary gland from embryonic stem cells.  Nature communications, 9(1):1–13, 2018.
  * [41] J. Tang and H. Liu.  CoSelect: feature selection with instance selection for social media data.  In Proc. of the SIAM Intl. Conference on Data Mining, 2013.
  * [42] S. Tong and D. Koller.  Support vector machine active learning with applications to text classification.  Journal of machine learning research, 2(Nov):45–66, 2001.
  * [43] T. Viklands.  Algorithms for the weighted orthogonal procrustes problem and other least squares problems.  2006\.
  * [44] X. Wei, B. Cao, and P. S. Yu.  Unsupervised feature selection on networks: A generative view.  In Thirtieth AAAI Conference on Artificial Intelligence, AAAI’16, page 2215–2221. AAAI Press, 2016.
  * [45] X. Wei, S. Xie, and P. S. Yu.  Efficient partial order preserving unsupervised feature selection on networks.  In Proceedings of the 2015 SIAM International Conference on Data Mining, pages 82–90. SIAM, 2015.
  * [46] Y. Yang, H. T. Shen, Z. Ma, Z. Huang, and X. Zhou.  L2,1-norm regularized discriminative feature selection for unsupervised learning.  In Proc. of the Twenty-second International Joint Conference on Artificial Intelligence (IJCAI), 2011.
  * [47] K. Yu, J. Bi, and V. Tresp.  Active learning via transductive experimental design.  In Proc. of the Intl. Conference on Machine Learning, ICML ’06, pages 1081–1088, New York, NY, USA, 2006. ACM.
  * [48] L. Zhang and P. Bogdanov.  DSL: discriminative subgraph learning via sparse self-representation.  In Proc. of the SIAM Intl. Conference on Data Mining (SDM), 2019\.