# Toward Clinically Trustworthy Deep Learning

Cooper Gamble+
Radiology Informatics Lab
Mayo Clinic
Rochester, MN 55904
<EMAIL_ADDRESS>
Shahriar Faghani+
Radiology Informatics Lab
Mayo Clinic
Rochester, MN 55904
<EMAIL_ADDRESS>
Bradley J. Erickson∗
Radiology Informatics Lab
Mayo Clinic
Rochester, MN 55904
<EMAIL_ADDRESS>

###### Abstract

As deep learning (DL) continues to demonstrate its ability in radiological
tasks, it is critical that we optimize clinical DL solutions to include
safety. One of the principal concerns in the clinical adoption of DL tools is
trust. This study aims to apply conformal prediction as a step toward
trustworthiness for DL in radiology. This is a retrospective study of 491 non-
contrast head CTs from the CQ500 dataset, in which three senior radiologists
annotated slices containing intracranial hemorrhage (ICH). The dataset was
split into definite and challenging subsets, where challenging images were
defined to those in which there was disagreement among readers. A DL model was
trained on 146 patients (10,815 slices) from the definite data (training
dataset) to perform ICH localization and classification for five classes of
ICH. To develop an uncertainty-aware DL model, 1,546 cases of the definite
data (calibration dataset) was used for Mondrian conformal prediction (MCP).
The uncertainty-aware DL model was tested on 8,401 definite and challenging
cases to assess its ability to identify challenging cases. After the MCP
procedure, the model achieved an F1 score of 0.920 for ICH classification on
the test dataset. Additionally, it correctly identified 6,837 of the 6,856
total challenging cases as challenging (99.7% accuracy). It did not
incorrectly label any definite cases as challenging. The uncertainty-aware ICH
detector performs on par with state-of-the-art models. MCP’s performance in
detecting challenging cases demonstrates that it is useful in automated ICH
detection and promising for trustworthiness in radiological DL.

+Co-first author, ∗Corresponding author

_Keywords_ deep learning $\cdot$ conformal prediction $\cdot$ intracranial
hemorrhage

##  1 Introduction

Intracranial hemorrhage (ICH) is a widespread and potentially life-threatening
condition, affecting more than 50,000 adults annually in the United States
(1). Deep learning (DL) has demonstrated promising results as an automated
solution to ICH detection (2). However, because DL generally do not provide a
confidence score for their predictions, it is difficult for both clinicians
and researchers to trust them (3).

In general, a radiological DL model provides some predictions (diagnoses) in
the form of numerical outputs for the possible diagnoses that were included in
the training dataset. These outputs possess the mathematical properties of
measures of probability, but they are not representative of the true
probability of a given outcome because they are uncalibrated (3,4).
Calibration is the process by which a model’s predictions are made to more
accurately approximate the underlying distribution of the population on which
it will be applied. Conventionally, calibration maps a model’s output to the
observed frequency of an event. In this study, however, calibration refers to
the process of learning the test distribution from a small sample of the
population called the calibration dataset. Conformal prediction (i.e.
distribution-free uncertainty quantification) is a family of statistical
methods, the premise of which is to generate prediction sets (i.e.,
differential diagnoses) which are guaranteed to contain the true value
(diagnosis) at a user-specified confidence level (5). Though conformal methods
vary in nature and implementation, the general process begins by defining a
heuristic notion of uncertainty (HNU) for the DL model. Then, for every
instance in the calibration dataset, this uncertainty metric is calculated and
recorded. Finally, the uncertainty values obtained during calibration are then
applied to some computation during inference to yield sets of
predictions—instead of point predictions—for each test instance. The nature of
this computation depends on the nature of the selected conformal method.

Mondrian conformal prediction (MCP) is a subset of conformal prediction which
relies on a grouping function to split the calibration data based on some
shared characteristic (6). MCP operates by ranking the calibration HNU values
across a given group, which allows for the calculation of conformal scores for
predictions on the test dataset. For each test sample, a conformal score is
obtained for each group by dividing the would-be rank (i.e., the insertion
index) of the sample’s HNU value within a given group by the size of that
group. This conformal score is then compared against a user-specified
threshold, and groups for which the sample’s score exceeds the threshold are
included in its final prediction set (see Figure 1).

Figure 1: High-level overview of Mondrian conformal prediction as it is
applied in this study. IPH = intraparenchymal hemorrhage, IVH =
intraventricular hemorrhage, HNU = heuristic notion of uncertainty, i =
insertion index in sorted value array, p = conformal score.

In addition to uncalibrated models, another source of mistrust in DL models is
poor handling of out-of-domain cases. Traditional training methods restrict
the apparent knowledge of a DL model to the limits of its training domain (7),
such that as inputs trend away from being in-domain, the model’s predictions
are increasingly likely to be incorrect. Other artifacts like patient motion
or implanted devices can also significantly alter images, causing them to be
out-of-domain. The out-of-domain problem is exacerbated when a model’s outputs
are not calibrated.

In this study, we utilize MCP as a method to improve trustworthiness in DL
models. We train an object detection algorithm to classify and localize ICH in
non-contrast head CT slices. We demonstrate that MCP achieves state-of-the-art
performance in detection of five classes of ICH, and we further show that an
MCP-augmented pipeline enables highly accurate identification of challenging
cases. We define challenging cases as those in which radiologists disagreed on
whether ICH was present. We propose that if a DL tool can generate
statistically guaranteed differential diagnoses and reliably identify when it
encounters a challenging input, then physicians will be increasingly able to
trust this tool and therefore better positioned to adopt it as an automated
solution in their clinical or research workflow.

##  2 Materials and Methods

###  2.1 Data Collection

This is a retrospective study of 491 anonymous, non-contrast head CT scans
collected at various centers in New Delhi, India and compiled as the publicly
available CQ500 dataset (8). The multi-institutional 2019 Radiological Society
of North America (RSNA) Brain CT Hemorrhage Challenge dataset was used to
externally validate the performance of the DL model (9). See Supplementary
Descriptions 1 and 2 for more information about the datasets used in this
study.

###  2.2 Data Organization

Out of the CQ500 dataset, scans for which all three readers indicated that
none of the five hemorrhage classes were present were excluded and used as an
evaluative control group. After narrowing the dataset to contain only scans
positive for at least one class of hemorrhage, the scans were split into
definite and challenging subsets. The definite dataset was defined to contain
scans on whose labels readers unanimously agreed. Any scan for which there was
not unanimity among the readers for the five hemorrhage classes was placed
into the challenging dataset. The definite dataset was then further divided
into training, tuning, calibration, and test subsets with a 70/10/10/10 split,
assuring that images were grouped by their series and patient identifiers to
prevent information leakage between subsets (10). Slices were also stratified
by their ground truth labels for all five hemorrhage types to maintain a class
distribution as close to the true distribution as possible. The data
inclusion, exclusion, and splitting process is depicted in Supplementary
Figure 1.

The external validation dataset was split into positive (hemorrhagic) and
negative (non-hemorrhagic) subsets. Because the negative subset was
significantly larger than the positive subset, and because there already
existed a negative control group from the CQ500 dataset, a random sample of
the negative subset was used during external validation. The sample size was
equal to the total number of positive slices.

See Supplementary Notes 1 and 2 for additional remarks about the data
organization process.

###  2.3 Image Processing

Windowing was applied to all scans in the final datasets at a window level of
50 and a window width of 80 to enhance hemorrhage prominence (11). The
windowed scans were then converted from their original Digital Imaging and
Communications in Medicine (DICOM) format into two-dimensional Portable
Network Graphics (PNGs) of 512x512 pixels. 21,653 out of 674,254 images in the
RSNA dataset could not be converted from DICOM to PNG format, so these images
were excluded from the final external validation dataset. During training, the
DL pipeline applied standard augmentations to reduce the effects of
overfitting (12). Augmentations were not applied during other phases of the
study. Additional information regarding the augmentations used can be found in
Supplementary Table 1.

###  2.4 Model Development

Ultralytics’ You Only Look Once (YOLO) v8 computer vision model was selected
as the DL tool for this task because it is well-regarded as state-of-the-art
in object detection (13). The YOLO pipeline outputs both the model’s
confidence score for a bounding box containing an object of interest, as well
as its confidence score for a detected object belonging to the given class.
The final confidence score for a given prediction is determined by a non-
linear combination of these two confidence values whose weighting is optimized
as an auxiliary component of the cost function during training. The model was
trained for 300 epochs with a batch size of six using the AdamW optimization
algorithm with an initial learning rate of 0.01 and YOLO’s cosine annealing
learning rate scheduler (14). At the end of each training epoch, the model was
applied to the tuning dataset to benchmark its performance. Additional
training parameters are listed in Supplementary Table 2.

###  2.5 Model Calibration

When the model completed training, inference was performed on the calibration
dataset. The model was allowed to predict up to 10,000 instances of hemorrhage
per slice because to populate class-wise Mondrian groups, a bounding box and
confidence score for each hemorrhage class was needed for every slice. These
predictions were then filtered by identifying the highest confidence bounding
box for each hemorrhage class, resulting in five confidence scores—one for
each hemorrhage class—for each calibration sample. Five Mondrian groups were
then formed by sorting the calibration confidence scores across all samples
for each hemorrhage class. An additional five groups were formed for the
absence of each hemorrhage class by sorting the additive complements of all
calibration confidence scores within each hemorrhage class. Figure 1
demonstrates the flow of MCP from calibration to inference.

###  2.6 Model Performance Evaluation

After generating Mondrian groups, the model was evaluated on the test,
challenging, negative control, and RSNA external validation datasets. During
inference for each dataset, the model was once again allowed to predict up to
10,000 instances of hemorrhage per slice (Figure 2, Step 1). The highest
confidence bounding boxes for each hemorrhage class were then selected,
yielding five bounding boxes—each of a different class label—for each sample
(Figure 2, Step 2). Next, clusters were created by calculating the
intersection over union (IoU) between each bounding box and all other bounding
boxes (Figure 2, Step 3). A cluster was defined by any bounding box whose IoU
with all other bounding boxes was less than the given threshold. Then, each
cluster was filled by reinspecting all initial predictions, excluding those
that had already been selected, and including any bounding box whose IoU with
the cluster-delimiting bounding box was greater than or equal to the given IoU
threshold (Figure 2, Step 4). Penultimately, each cluster was condensed to
contain only five bounding boxes by retaining the cluster-delimiting bounding
box and the highest confidence bounding box for the remaining four hemorrhage
classes (Figure 2, Step 5). Finally, conformal scores were calculated for each
bounding box, and conformalized prediction sets were generated for each
cluster (Figure 2, Step 6).

Figure 2: Algorithm for clustering predictions with IoU thresholding and non-
maximum suppression. IPH = intraparenchymal hemorrhage, IVH = intraventricular
hemorrhage, SDH = subdural hemorrhage, EDH = epidural hemorrhage, SAH =
subarachnoid hemorrhage, C1 = Cluster 1, C2 = Cluster 2, C3 = Cluster 3, C4 =
Cluster 4, IoU = intersection over union.

Because this algorithm relied on two threshold values—one for the IoU and one
for the conformal score—an optimization process was conducted with the test
dataset to determine which thresholds to apply during evaluation on the
challenging, negative, and external validation datasets. Three confusion
matrices were defined. Each matrix tested different criteria for a prediction
to be counted as a true positive (TP), false positive (FP), true negative
(TN), or false negative (FN). Precise textual definitions of TP, FP, TN, and
FN for each confusion matrix can be found in Figure 3, and visual examples of
each case are depicted in Figure 4.

Figure 3: Textual descriptions of TP, FP, TN, and FN for the three confusion
matrices used during model evaluation. TP = true positive, FP = false
positive, TN = true negative, FN = false negative, IoU = intersection over
union. Figure 4: Visual depictions of examples of TP, FP, TN, and FN for the
three confusion matrices used during model evaluation. TP = true positive, FP
= false positive, TN = true negative, FN = false negative, IPH =
intraparenchymal hemorrhage, IVH = intraventricular hemorrhage, SDH = subdural
hemorrhage, EDH = epidural hemorrhage, SAH = subarachnoid hemorrhage.

During optimization, IoU and conformal score thresholds were varied between
0.00 and 1.00 at intervals of 0.05. The clustering algorithm described above
and depicted in Figure 2 was then applied using each combination of IoU and
Conformal score thresholds, yielding 441 conformalized prediction sets
(differential diagnoses) for each test instance. Then, for Matrices A and B,
the IoU and Conformal score thresholds which yielded the maximum the area
under the receiver operating characteristic curve (AUROCC) were selected. For
Matrix C, because only TP and FP were defined, the threshold values which
yielded the maximum positive predictive value (PPV) were selected. To avoid
confusion in the Results section, note that the threshold values obtained
during optimization for both Matrix A and Matrix B were identical.

These threshold values were then applied to the challenging, negative, and
external validation datasets. For the challenging dataset, all four counts
(TP, FP, TN, FN), as well as sensitivity, specificity, PPV, negative
predictive value (NPV), and the F1-score were reported for the thresholds
obtained from Matrices A and B. For Matrix C, TN and FN were not well-defined,
so only PPV was reported. For the control group, only the TN and FP counts and
specificity were reported. Other metrics were not defined because this dataset
only contained non-hemorrhagic slices. No bounding boxes were available for
the external validation dataset, so the evaluation was agnostic to
localization performance and assessed only the model’s classification ability.
As such, the definitions of TP, FP, TN, and FN were adjusted to yield
meaningful results (see Figure 5). For classification, sensitivity,
specificity, PPV, NPV, and the F1-score were reported. For exact
classification, because FN was not well-defined, only specificity and PPV were
reported.

Figure 5: Textual descriptions of TP, FP, TN, and FN for the two confusion
matrices used during external validation of the model. TP = true positive, FP
= false positive, TN = true negative, FN = false negative.

In addition to the metrics described above, mean average precision (mAP) was
reported for all datasets, except the control group because it contained no
positive samples and the external validation dataset because bounding boxes
were unavailable. mAP is a standard metric in object detection and though it
does not paint a full picture of a model’s performance, it can offer insight
into its positive predictive ability (14,15).

###  2.7 Trustworthiness Evaluation

This study also evaluated the performance of MCP as a mechanism for
trustworthiness in DL. The accuracy of the model in identifying challenging
cases was reported, where a sample was predicted as challenging if the final
prediction set for any cluster within the sample contained either of the
following: both the presence and absence of one hemorrhage class, or the
presence of more than one hemorrhage class. The number of samples in the
definite test dataset which the model incorrectly labeled challenging was also
reported to assess the model’s ability to distinguish between definite and
challenging cases.

###  2.8 Statistical Analysis

All calculations were performed using the NumPy (version 1.24.4) and scikit-
learn (version 1.2.2) libraries in Python (version 3.8.16).

##  3 Results

###  3.1 Dataset Characteristics

The definite dataset was split such that there were 10,815 training slices,
1,546 tuning slices, 1,546 calibration slices, and 1,545 test slices. The
challenging dataset comprised 6,856 slices. Instance frequencies for each
hemorrhage class within each dataset are shown in Table 1.

Class | Training | Tuning | Calibration | Test | Challenging
---|---|---|---|---|---
IPH | 4,743 | 664 | 664 | 683 | 732
IVH | 1,415 | 219 | 208 | 187 | 407
SDH | 3,183 | 432 | 466 | 451 | 4,251
EDH | 273 | 40 | 40 | 39 | 174
SAH | 5,237 | 754 | 765 | 745 | 3,557
Combined | 14,581 | 2,109 | 2,143 | 2,105 | 9,121

Table 1: Frequency of all hemorrhage classes for each dataset. IPH =
intraparenchymal hemorrhage, IVH = intraventricular hemorrhage, SDH = subdural
hemorrhage, EDH = epidural hemorrhage, SAH = subarachnoid hemorrhage.

Of the 652,601 slices from the RSNA dataset successfully converted to PNGs,
there were 32,564 instances of intraparenchymal hemorrhage (IPH), 23,766
instances of intraventricular hemorrhage (IVH), 42,496 instances of subdural
hemorrhage (SDH), 2,761 instances of epidural hemorrhage (EDH), and 32,122
instances of subarachnoid hemorrhage (SAH).

###  3.2 Model Development Cost

Training the model for 300 epochs required 10.867 hours on three NVIDIA
A100-SXM4-80GB cards. Inference required approximately 40 additional hours on
one of these cards.

###  3.3 Model Performance: Hemorrhage Detection

mAP values at an IoU threshold of 0.95<EMAIL_ADDRESS>for each hemorrhage class
are displayed for each dataset in Table 2.

Class | Training | Tuning | Calibration | Test | Challenging
---|---|---|---|---|---
IPH | 0.979 | 0.861 | 0.876 | 0.873 | 0.0492
IVH | 0.958 | 0.804 | 0.798 | 0.851 | 0.040
SDH | 0.976 | 0.891 | 0.880 | 0.884 | 0.0538
EDH | 0.964 | 0.889 | 0.877 | 0.867 | 0.0379
SAH | 0.980 | 0.837 | 0.842 | 0.864 | 0.00659
Combined | 0.972 | 0.856 | 0.855 | 0.868 | 0.0375

Table 2: mAP values at an IoU threshold of 0.95<EMAIL_ADDRESS>of all hemorrhage
classes for each dataset. IPH = intraparenchymal hemorrhage, IVH =
intraventricular hemorrhage, SDH = subdural hemorrhage, EDH = epidural
hemorrhage, SAH = subarachnoid hemorrhage.

For the test and challenging datasets, confusion matrices with IoU and
Conformal score thresholds optimized on the test dataset are shown in Table 3.

Dataset | Metric | Matrix A | Matrix B | Matrix C
---|---|---|---|---
Test | TP | 1,952 | 1,950 | 347
FP | 187 | 187 | 30
TN | 1,358 | 1,358 | —
FN | 153 | 155 | —
Sensitivity | 0.927 | 0.926 | —
Specificity | 0.879 | 0.879 | —
PPV | 0.913 | 0.912 | 0.920
NPV | 0.899 | 0.898 | —
F1 | 0.920 | 0.919 | —
Challenging | TP | 1,325 | 1,325 | 15
FP | 2,222 | 2,222 | 107
TN | 4,634 | 4,634 | —
FN | 7,796 | 7,796 | —
Sensitivity | 0.145 | 0.145 | —
Specificity | 0.676 | 0.676 | —
PPV | 0.374 | 0.374 | 0.123
NPV | 0.323 | 0.323 | —
F1 | 0.209 | 0.209 | —

Table 3: Confusion matrices and derived metrics for the test and challenging
datasets. Test metrics were calculated by optimizing the IoU and Conformal
score thresholds for AUROCC (Matrices A and B) and (Matrix C). Challenging
metrics were calculated using the optimal thresholds as determined by the test
dataset. TP = true positive, FP = false positive, TN = true negative, FN =
false negative, PPV = positive predictive value, NPV = negative predictive
value, IoU = intersection over union, AUROCC = area under the receiver
operating characteristic curve.

Table 4 depicts the model’s performance on the negative and external
validation datasets using the same optimal threshold values.

Dataset | Metric | Matrices A and B Thresholds | Matrix C Thresholds
---|---|---|---
Negative Control | TN | 61,637 | 82,074
FP | 20,684 | 247
Specificity | 0.749 | 0.997
External Validation (Classification) | TP | 131,994 | 94,225
FP | 60,088 | 839
TN | 127,748 | 186,997
FN | 55,842 | 93,611
Sensitivity | 0.703 | 0.502
Specificity | 0.680 | 0.995
PPV | 0.687 | 0.991
NPV | 0.696 | 0.666
F1 | 0.695 | 0.666
External Validation (Exact Classification) | TP | 27,295 | 302
FP | 43,934 | 668
TN | 60,284 | 93,624
Specificity | 0.578 | 0.993
PPV | 0.383 | 0.311

Table 4: Performance metrics for the negative control and external validation
datasets. Metrics are reported for two sets of threshold values: the first set
was determined by optimizing for AUROCC (Matrices A and B), and the second set
was determined by optimizing for PPV (Matrix C). TP = true positive, FP =
false positive, TN = true negative, FN = false negative, PPV = positive
predictive value, NPV = negative predictive value, AUROCC = area under the
receiver operating characteristic curve.

###  3.4 Model Performance: Identification of Challenging Cases

The model achieved an identification accuracy of 99.7% (6,837 out of 6,856) at
the optimized IoU and conformal score thresholds. The model labeled zero of
the 1,545 definite test dataset samples as challenging. This finding is
significant because the model performs markedly worse (p < 0.00001) on
challenging cases than definite ones, as shown in Table 2. This demonstrates
that the MCP pipeline can reduce the rate of false predictions by flagging
challenging cases and referring them to a human reviewer. An example of the
model’s failure on a challenging case is shown in Figure 6.

Figure 6: Example of challenging case prediction and ground truth. IPH =
intraparenchymal hemorrhage, SDH = subdural hemorrhage, SAH = subarachnoid
hemorrhage.

##  4 Discussion

This study aims to propose a method for increasing trust in medical DL tools
by demonstrating the efficacy of MCP in the task of ICH detection and
supplying a deployable toolkit for applications in other tasks. The high
accuracy of the proposed approach in identifying challenging cases indicates
that MCP is a promising method for enhancing trust in radiological AI systems.
Furthermore, the GitHub repository and usage tutorial released to accompany
this work offers a practical, scalable implementation of MCP to further
validate its performance.

DL-assisted ICH detection was first proposed in 2017 with Phong, et al.
demonstrating that a convolutional neural network could be trained to classify
the presence or absence of ICH with 98.9% average accuracy (16). Recently,
Teneggi, et al. assessed the performance differences between strongly- and
weakly- supervised models in ICH classification using the RSNA dataset as its
primary dataset and the CQ500 dataset as external validation (17). They
reported an AUROCC of 0.96 on their tuning subset of the RSNA dataset and of
0.95 on the CQ500 dataset. Additionally, Umapathy, et al. obtained an F1 score
of 0.97 on the CQ500 dataset with a convolutional ensemble trained on the RSNA
dataset (18). As for conformal prediction, Vovk, et al. introduced the
underlying theory in 2005 (19). More recently, Angelopoulos and Bates released
a comprehensive guide to conformal methods in modern DL (5), and Ding, et al.
proposed class-conditional conformal prediction (6), from which the
methodologies used in this study are derived.

This work presents a DL solution to ICH detection that performs on par with
state-of-the-art models and introduces MCP as a method of building trust with
radiological DL tools. The uncertainty-aware localization and classification
model generates statistically guaranteed differential diagnoses, and the
integration of MCP enables it to flag challenging cases with near-perfect
accuracy for review by radiologists. Moreover, the DL architecture used in
this solution makes it easily adaptable to a wide array of radiological and
medical tasks, especially considering that the IoU and conformal score
thresholds can be optimized uniquely for each task to which this pipeline is
applied. This opens the door for further investigation into the critical
frontier of trust in DL. The publicized code and deployment instructions
further encourage the continued research and development of MCP as a mechanism
for trust in radiological DL.

The DL model operates on two-dimensional CT slices, so it is not immediately
compatible with current Picture Archiving and Communication Systems (PACS).
Given the current state of radiological DL, however, this improvement is
certainly feasible without significant modification to the MCP backbone of the
trust component.

To deem a DL algorithm trustworthy, it must be able to reliably communicate
when it encounters cases about which it is uncertain. This study not only
demonstrates that an MCP-based approach can yield up to 99.7% accuracy in the
identification of challenging cases, but also that a trustworthy DL model
achieves performance metrics on par with state-of-the-art models. Further
efforts in trustworthy radiological DL will include increasing practicality by
expanding the application of MCP to three-dimensional models, as well as
continuous external validation to account for generalizability. Trustworthy DL
is the necessary next step toward the feasible clinical adoption of AI tools
in radiology.

## References

  * [1] Caceres JA, Goldstein JN. Intracranial hemorrhage. Emerg Med Clin North Am. 2012;30(3):771–794.
  * [2] Yeo M, Tahayori B, Kok HK, et al. Review of deep learning algorithms for the automatic detection of intracranial hemorrhages on computed tomography head imaging. J Neurointerv Surg. 2021;13(4):369–378.
  * [3] Faghani S, Moassefi M, Rouzrokh P, et al. Quantifying Uncertainty in Deep Learning of Radiologic Images. Radiology. 2023;308(2):e222217.
  * [4] Wang C. Calibration in Deep Learning: A Survey of the State-of-the-Art. arXiv [cs.LG]. 2023. doi: 10.48550/ARXIV.2308.01222.
  * [5] Angelopoulos AN, Bates S. Conformal Prediction: A Gentle Introduction. Foundations and Trends® in Machine Learning. Now Publishers; 2023;16(4):494–591.
  * [6] Ding T, Angelopoulos AN, Bates S, Jordan MI, Tibshirani RJ. Class-Conditional Conformal Prediction with Many Classes. arXiv [stat.ML]. 2023. doi: 10.48550/ARXIV.2306.09335.
  * [7] Janiesch C, Zschech P, Heinrich K. Machine learning and deep learning. Electronic Markets. 2021;31(3):685–695.
  * [8] Chilamkurthy S, Ghosh R, Tanamala S, et al. Deep learning algorithms for detection of critical findings in head CT scans: a retrospective study. Lancet. 2018;392(10162):2388–2396.
  * [9] Flanders AE, Prevedello LM, Shih G, et al. Construction of a Machine Learning Dataset through Collaboration: The RSNA 2019 Brain CT Hemorrhage Challenge. Radiol Artif Intell. 2020;2(3):e190211.
  * [10] Rouzrokh P, Khosravi B, Faghani S, et al. Mitigating Bias in Radiology Machine Learning: 1. Data Handling. Radiol Artif Intell. 2022;4(5):e210290.
  * [11] Baba Y, Murphy A. Windowing (CT). Radiopaedia.org. Radiopaedia.org; 2017. doi: 10.53347/rid-52108.
  * [12] Zhang K, Khosravi B, Vahdati S, et al. Mitigating Bias in Radiology Machine Learning: 2. Model Development. Radiol Artif Intell. 2022;4(5):e220010.
  * [13] ultralytics: NEW - YOLOv8 in PyTorch > ONNX > OpenVINO > CoreML > TFLite. Github; https://github.com/ultralytics/ultralytics. Accessed January 10, 2024.
  * [14] Henderson P, Ferrari V. End-to-end training of object class detectors for mean average precision. Computer Vision – ACCV 2016. Cham: Springer International Publishing; 2017. p. 198–213.
  * [15] Faghani S, Khosravi B, Zhang K, et al. Mitigating Bias in Radiology Machine Learning: 3. Performance Metrics. Radiol Artif Intell. 2022;4(5):e220061.
  * [16] Phong TD, Duong HN, Nguyen HT, et al. Brain Hemorrhage Diagnosis by Using Deep Learning. Proceedings of the 2017 International Conference on Machine Learning and Soft Computing. New York, NY, USA: Association for Computing Machinery; 2017. p. 34–39.
  * [17] Teneggi J, Yi PH, Sulam J. Examination-level Supervision for Deep Learning–based Intracranial Hemorrhage Detection at Head CT. Radiology: Artificial Intelligence. Radiological Society of North America; 2023;e230159.
  * [18] Umapathy S, Murugappan M, Bharathi D, Thakur M. Automated Computer-Aided Detection and Classification of Intracranial Hemorrhage Using Ensemble Deep Learning Techniques. Diagnostics (Basel). 2023;13(18). doi: 10.3390/diagnostics13182987.
  * [19] Vovk V, Gammerman A, Shafer G. Algorithmic Learning in a Random World. Springer Science & Business Media; 2005.

##  5 Supplementary Materials

Supplementary Description 1: CQ500 Dataset
The Centre for Advanced Research in Imaging, Neurosciences, and Genomics is
responsible for the compilation of the CQ500 dataset. Scans were collected in
two batches, B1 and B2. B1 comprised all head CTs conducted at a subset of the
collection centers from 20 November 2017 to 20 December 2017. B2 scans were
selectively included from a much larger dataset of scans obtained from the
entire set of collection centers over an unspecified date range. First,
radiology reports for scans in the larger dataset were processed with a
natural language processing algorithm to identify scans with IPH, IVH, SDH,
EDH, SAH, and calvarial fracture. Scans with reports containing mention of any
of these conditions were then randomly selected such that there were
approximately 80 scans for each intracranial condition. Scans were excluded if
they contained post-operative defects, did not contain a non-contrast series
with axial cuts and a soft reconstruction kernel covering the complete brain,
or the patient was less than seven years of age. The final CQ500 dataset was
read by three senior radiologists with eight, 12, and 20 years of experience
in cranial CT interpretation. Each reader provided labels for the presence or
absence of IPH, IVH, SDH, EDH, and SAH, as well as positive or negative
indicators for chronic SDH, calvarial fracture, midline shift, and mass effect
(8).

Supplementary Description 2: RSNA External Validation Dataset
The RSNA dataset is a multi-institutional effort with images contributed by
Stanford University, Universidade Federal de São Paulo, and Thomas Jefferson
University. Scans were labeled by a team of sixty junior and senior
radiologists who were members of the American Society of Neuroradiology. The
RSNA dataset was split between two stages for Kaggle competition; this study
incorporates the 674,254 training slices from Stage 1. Annotation was
performed on a commercial, browser-based annotation platform called MD.ai.
Readers used their personal computers to access the software, and each reader
was blind to all other readers’ annotations. Users were directed to go through
a series of axial brain CT images and choose from five hemorrhage subtype
image labels: IPH, IVH, SDH, EDH, and SAH. They could also label the initial
and final slices of a specific hemorrhage and let the software interpolate the
labels for the slices in between. In cases where multiple hemorrhage subtypes
were found, multiple labels were assigned to the appropriate slices. If no
hemorrhagic features were present in the entire examination, the annotator
could assign one of two labels: normal or abnormal/not hemorrhage (such as
stroke, atrophy, white matter disease, hydrocephalus, tumor). An extra flag
related to the examination was provided to indicate incomplete or incorrectly
identified body parts (9). Annotators underwent onboarding sessions in groups
of 10. Practice cases were utilized to acquaint readers with both the tool and
the annotation procedure. Each volunteer was given access to 24 training
tests, which they labeled using the annotation tool. Ratings for both
examination and image-based annotations were documented for every annotator.

Supplementary Note 1: Ambiguous Labeling
There were two small subsets of scans with ambiguous labeling: one subset
contained scans which readers identified as containing hemorrhage but for
which no bounding boxes were provided in the annotations, and the other subset
contained scans for which bounding boxes were provided but which readers had
identified as containing no hemorrhage. The first subset was excluded the from
the study and the second subset was placed into the challenging dataset.

Supplementary Note 2: Definite-Challenging Dataset Overlap
There was overlap between the definite and challenging datasets because
concordance was determined for each hemorrhage class, such that if for a
single scan readers agreed on one hemorrhage label but disagreed on another,
that scan was placed into the definite dataset with concordant label and into
the challenging dataset with the disagreed-upon label.

Supplementary Table 1: Augmentations Applied to Training Images

Augmentation Description | Value
---|---
Image Hue Saturation Value (HSV) Hue adjustment (fraction) | 0.015
Image HSV-Saturation adjustment (fraction) | 0.700
Image HSV-Value adjustment (fraction) | 0.400
Image rotation (+/- degrees) | 0.000
Image translation (+/- fraction) | 0.100
Image scale (+/- gain) | 0.500
Image shear (+/- degrees) | 0.000
Image perspective (+/- fraction) | 0.000
Image flip along y-axis (probability) | 0.000
Image flip along x-axis (probability) | 0.500
Image mosaic (probability) | 0.500
Image mix-up (probability) | 0.000

Supplementary Table 2: YOLOv8 Training Hyperparameters

Hyperparameter | Value
---|---
Seed | 829
Learning Rate | 0.01
Learning Rate Scheduler | Cosine Annealing
Non-Max Suppression Threshold | 0.45
Confidence Threshold | 0.25
Intersection over Union Threshold | 0.50
Convolutional Filters per Layer | 64
Convolutional Layers | 53
Activation Function | Leaky Rectified Linear Unit
Number of Anchors | 3
Network Architecture | YOLOv8x
Optimizer | AdamW
Cost Function | Varifocal Loss

Supplementary Figure 1: Flowchart of Data Inclusion, Exclusion, and Splitting

Supplementary Figure 1: Depiction of data inclusion, exclusion, and splitting.
RSNA = Radiological Society of North America, DICOM = Digital Imaging and
Communications in Medicine, PNG = Portable Network Graphics.


We consider the problem of learning low-rank tensors from partial observations with structural constraints, and propose a novel factorization of such tensors, which leads to a simpler optimization problem. The resulting problem is an optimization problem on manifolds. We develop first-order and second-order Riemannian optimization algorithms to solve it. The duality gap for the resulting problem is derived, and we experimentally verify the correctness of the proposed algorithm. We demonstrate the algorithm on nonnegative constraints and Hankel constraints.

§ INTRODUCTION

With the rise in the availability of multidimensional data such as colour images, video sequences, and hyperspectral images, tensor-based techniques have started gaining attention as traditional matrix-based methods cannot exploit the underlying structure present in higher-dimensional data. Many recent applications of tensor reconstruction techniques also enforce structural constraints such as nonnegativity ([32], [48]) or a Hankel structure ([45], [46]).

We propose a framework for dealing with tensor completion problems with general structural constraints. In particular, we consider the structured low-rank tensor learning problem of the following form:

\begin{equation}
\begin{aligned}
\underset{\W\in\Rtensor}{\text{min}}& &&C\|\W_\Omega-\Y_\Omega\|^2 + R(\W)\\
\text{subject to}&\hspace{0.5cm} &&A(\W) \ge 0, \label{eqn:primal_problem}
\end{aligned}
\end{equation}

where $\Y_\Omega\in \Rtensor$ is a partially observed tensor for indices given in the set $\Omega$, $(\W_\Omega)_{i_1,\ldots,i_K} = w_{i_1,\ldots,i_K}$ for $({i_1,\ldots,i_K})\in \Omega$, $C>0$ denotes the cost parameter, $\mapping{R}{\Rtensor}{\R}$ is a regularizer that induces low-rank constraint on the tensor, and $\mapping{A}{\Rtensor}{\R^n}$ is a linear map that induces structural constraint on the tensor.

Following [27], we learn the tensor as a sum of $K$ tensors, $\W = \W^{(1)} + \cdots + \W^{(K)}$, and use the following low-rank regularizer:

\begin{equation}
R(\W) = \sum_{k=1}^K \dfrac{1}{\lambda_k}\Norm{W_k^{(k)}}^2_*.
\end{equation}

The main contributions of the paper are given below.

* We propose a novel factorization for modeling structured low-rank tensors through a partial dual problem of (<ref>).

* We develop first-order and second-order Riemannian optimization algorithms that exploit proposed factorization's inherent geometric structure.

* We compute the expression for the duality gap and verify the correctness of the proposed algorithm through experiments.

* We apply the proposed algorithm to the nonnegative constraint and the Hankel constraint.

§ NOTATION

We follow [2] for our tensor notation. We present a few important notions here.

Tensors are denoted by uppercase calligraphic letters, e.g., $\W$. Matrices are denoted by uppercase letters, e.g., $X$. For a square matrix $X\in\R^{n\times n}$, we denote its trace by $tr(X)$. $\R_+$ denotes the interval $[0,\infty)$.
The inner product of two tensors is defined as follows:

\begin{equation*}
\innerproduct{\W}{\mathcal{U}} = \sum_{i_1=1}^{n_1}\sum_{i_2=1}^{n_2}\cdots \sum_{i_K=1}^{n_K} w_{i_1,\ldots,i_K}u_{i_1,\ldots,i_K}.
\end{equation*}

A mode-$k$ fiber of a tensor $\W\in\Rtensor$, denoted by $w_{i_1,\ldots,i_{k-1},:,i_{k+1},\ldots,i_K}$, is a vector obtained by fixing all but $k$-th index of $W$. The mode-$k$ unfolding of a tensor $\W\in\Rtensor$ is a matrix $unfold_k(\W) = W_k\in \R^{n_k\times n_1\cdots n_{k-1}n_{k+1}\cdots n_K}$ formed by arranging the mode-$k$ fibers to be the columns of the resulting matrix. Similarly, we can define the $k$-mode folding operation ($fold_k$) as the inverse of the unfolding operation - it converts a given matrix to a tensor of a suitable order.

The $k$-mode product of a tensor $\W \in \Rtensor$ with a matrix $U \in \R^{m \times n_k}$ is defined as follows:

\begin{equation*}
(\W\times_k U)_{i_1,\ldots,i_{k-1},j,i_{k+1},\ldots,i_K} = \sum_{i_k}^{n_k} w_{i_1,\ldots,i_K}u_{j,i_k}, \;\text{i.e.},\;\mathcal{X} = \W \times_k U \iff X_k = U W_k.
\end{equation*}

§ RELATED WORK

There are a number of well-known methods for tensor completion using different approaches for enforcing the low-rank constraint. Many methods use a generalization of the matrix trace norm regularizer to the tensor case. Two well-known regularizers are the overlapped trace norm regularizer ([29], [34] and [33]) and latent trace norm regularizer ([38]). The overlapped trace norm regularizer  uses the regularizer $R(\W) = \sum_{k=1}^{K} \|W_k\|_*$. In the latent trace norm regularizer model, the tensor is modelled as the sum of $K$-tensors, $\W = \sum_{k=1}^K \W^{(k)}\in\Rtensor$, and the regularizer is defined as

\begin{equation*}
R(\W) = \inf_{\substack{\W^{(k)},\\ k \in \{1,\ldots,K\}}} \sum_{k=1}^{K} \Norm{W^{(k)}_k}_*.
\end{equation*}

Other methods for tensor completion use tensor decompositions like the Tucker decomposition and CP decomposition to learn the tensors. [30] and [31] formulate the tensor completion problem as an optimization problem on the Riemannian manifold of fixed multi-linear rank tensors.

A general formulation for low-rank matrix completion problems with structural constraints was developed in [36]. It provided a unified framework for dealing with general linear inequality and equality constraints.

Recently, tensor completion problems with structural constraints have started attracting attention. [32] and [48] develop nonnegative tensor completion algorithms for image and video reconstruction tasks. [37] falls under the current framework which considers the special case of nonnegative constraints in detail with vast experimentation. [45] models image completion with missing slices as a higher-order Hankel tensor completion problem. [47] uses low-rank Hankel tensor model for estimating traffic states from partial observations.

§ DUAL FRAMEWORK

Following [27], we construct a partial dual problem to the primal problem (<ref>), incorporating the structural constraint into the formulation. We do this using the approach outlined in [36]. We generalize Theorem 1 of [27] to the case of tensors with structural constraints. The following lemma [28] is used in the development of the dual formulation.

For a matrix $X \in \R^{d \times T}$, the nuclear norm of $X$ satisfies the following relation:

\begin{equation*}
\|X\|_*^2 = \min_{\Theta \in P^d, \, \text{range}(X) \subseteq \text{range}(\Theta)} \innerproduct{\Theta^{\dagger}X}{X},
\end{equation*}
where $P^d = \{ S \in \R^{d \times d}:\, S \succeq 0, \text{tr}(S) = 1\} $, $\text{range}(\Theta) = \{\Theta z:\, z \in \R^d\}$, $\Theta^\dagger$ denotes the pseudo-inverse of $\Theta$. For a given $X$, the minimizer is $\bar{\Theta} = \sqrt{XX^T} / \text{tr}(\sqrt{XX^T}) $.

Using the above lemma, we can write (<ref>) as

\begin{equation}
\begin{aligned}
\underset{\substack{\Theta_k \in P^{n_k},\\ k \in \{1,\ldots,K\}}}{\text{min}}\;\underset{\substack{\W^{(k)},\\ k \in \{1,\ldots,K\}}}{\text{min}}& &&C\Norm{\bigg(\sum_{k=1}^K \W^{(k)}\bigg)_\Omega-\Y_\Omega}^2 + \sum_{k=1}^K \frac{1}{2\lambda_k} \innerproduct{\Theta_k^\dagger W^{(k)}_k}{W^{(k)}_k}\\
\text{subject to}&\hspace{0.5cm} &&A(\W) \ge 0. \label{eqn:primal_rewritten}
\end{aligned}
\end{equation}

The following minimax problem is equivalent to the problem (<ref>):

k ∈{1,…,K}min    ∈𝒞,   s∈^nmax            _Ω - 14C^2 -∑_k=1^K λ_k/2 Z_k + (A^*(s))_kΘ_k[Z_k + (A^*(s))_k],

where $\mathcal{C} = \{\Z\in\Rtensor:\,\Z = \Z_\Omega\}$.

The optimal solution $\bar{\W}$ of (<ref>) is related to the optimal solution $\{\bar{\Theta}_1,\ldots,\bar{\Theta}_K,\bar{\Z},\bar{s}\}$ of (<ref>) by

\begin{equation}\label{eqn:primal_dual_relation}
\bar{\W} = \sum_{k=1}^K \lambda_k(\bar{\Z}+A^*(\bar{s}))\times_k\bar{\Theta}_k.
\end{equation}

For proof, see Appendix <ref>.

From (<ref>) it can be seen that the structured and low-rank constraints on $\W$ can be decomposed into structured constraints on $s$ and low-rank constraints on $\Theta_k,$ which leads to a simpler optimization method.

§ PROPOSED ALGORITHM

Due to the low-rank constraint on $\W$, each $\Theta_k$ has a low rank. Therefore, a fixed-rank parameterized problem can be given by writing $\Theta_k = U_k U_k^T$. The problem (<ref>) can then be written as follows.

\begin{equation}\label{eqn:outer_manifold_problem}
\underset{U\in S^{n_1}_{r_1}\times\cdots\times S^{n_K}_{r_K}}{\min}\;g(U),
\end{equation}

where $U = (U_1,\ldots,U_K)$, $S^{n}_{r} = \{U\in\R^{n\times r}:\,\|U\|_F = 1\},$ and

\begin{equation}\label{eqn:inner_strong_convex_problem}
g(U) = \underset{\Z\in\mathcal{C},s\in\R^n}{\max}\; \innerproduct{\Z}{\Y_\Omega} - \dfrac{1}{4C}\|\Z\|^2 -\sum_{k=1}^K \frac{\lambda_k}{2} \Norm{U_k^T(Z_k + (A^*(s))_k)}^2.
\end{equation}

The optimization problem in (<ref>) is strongly convex for a given $U$, while problem (<ref>) is a non-convex problem in $U$.

The set $S^{n_1}_{r_1}\times\cdots\times S^{n_K}_{r_K}$ is a Riemannian manifold ([26], [27]), and thus problem (<ref>) is an optimization problem on a manifold. We solve it using either Riemannian conjugate gradient or Riemannian Trust-Region algorithm, depending on the structural constraint. The proposed algorithm is shown in Algo. <ref>. For more details on optimization on general manifolds, we refer the reader to [1] and [4].

The Riemannian optimization algorithm in Algo. <ref> requires computing the Euclidean gradient and its directional derivative, which are given in the following lemma.

Let $\{\hat{\Z},\hat{s}\}$ be the maximizer of the convex problem (<ref>) at $U$. Then, the Euclidean gradient $\nabla g(U)$ is given by

\begin{equation*}
\nabla g(U) = - (\lambda_1 P_1, \ldots, \lambda_K P_K),
\end{equation*}

where $P_k = (\hat{Z}_k + (A^*(\hat{s}))_k)(\hat{Z}_k+(A^*(\hat{s}))_k)^T U_k$. Let $V\in \R^{n_1\times r_1}\times\cdots\times \R^{n_K\times r_K}$ and $\Dot{Z}_k, \Dot{s}$ denote the directional derivatives of $Z_k$ and $s$ along $V$ respectively. Then, the directional derivative of $\nabla g$ at $U$ along $V$ is

\begin{equation*}
D\nabla g(U)[V] = - (\lambda_1 Q_1, \ldots, \lambda_K Q_K),
\end{equation*}

where $Q_k = (\hat{Z}_k + (A^*(\hat{s}))_k)(\hat{Z}_k+(A^*(\hat{s}))_k)^T V_k + 2 \, sym((\Dot{Z}_k + (A^*(\Dot{s}))_k)(\hat{Z}_k+(A^*(\hat{s}))_k)^T))U_k$ and $sym(X) = (X+X^T)/2$.

It can be seen that computing $D\nabla g(U)[V]$ requires the terms $\Dot{\Z}$ and $\Dot{s}$. These terms can be computed by applying directional derivative along $V$ on the first-order optimality conditions of the problem (<ref>) at $\{\hat{\Z},\hat{s}\}$.

Proposed Algorithm for Structured Low-Rank Tensor Completion
$\Y_\Omega$, rank=$(r_1,\ldots,r_K)$, $\varepsilon$, $(\lambda_1,\ldots,\lambda_K)$
$\hat{\W} = \sum_{k=1}^{K} \lambda_k(\hat{\Z}+A^*(\hat{s})) \times_k (U_k U_k^T)$

$t = 1,2,\cdots$
Check Termination: if $\|\nabla g(U^{(t)})\| \le \varepsilon$ then break
Solve for $\hat{\Z}^{(t)}$ and $\hat{s}^{(t)}$ in (<ref>)
Compute cost $g(U^{(t)})$, gradient $\nabla g(U^{(t)})$ and the directional derivative of $\nabla g(U^{(t)})$
Update $U$: $U^{(t+1)}$ = RiemannianCG-update$(U^{(t)})$ or
$U^{(t+1)}$ = RiemannianTR-update$(U^{(t)})$

We have the following result regarding the optimality of the proposed algorithm. It is a generalization of Theorem 3 in [36].

Let $\hat{U} = \left(U_1, \ldots, U_K \right)$ be a feasible solution of problem (<ref>) and $\left\{ \hat{\Z},\hat{s}\right\}$ be the maximizer of the convex problem (<ref>) at $U = \hat{U}$. Let $\hat{\A} = A^{*}(\hat{s})$, and $\sigma_k$ be the maximum singular value of $\hat{Z}_k + \hat{A}_k$. Let $\hat{\Theta} = (\hat{\Theta}_1,\ldots,\hat{\Theta}_K)$, where $\hat{\Theta}_k = \hat{U}_k \hat{U}_k^T$. Then,
$\left\{\hat{\Theta},\hat{\Z},\hat{s}\right\}$ is a candidate solution for the partial dual problem (<ref>) and we have the following expression for the duality gap $\Delta$:

\begin{equation}
\Delta = \sum_{k=1}^{K} \frac{\lambda_k}{2} \left( \sigma_k^2 -  \|\hat{U}_k^T(\hat{Z}_k+\hat{A}_k)\|^2 \right).\label{eqn:duality_gap}
\end{equation}

For proof see Appendix <ref>.

§ APPLICATIONS

We consider several popular applications for our proposed method. See [27] for the case where there are no structure constraints.

§.§ Nonnegative Tensor completion

The nonnegative tensor completion problem is

\begin{equation}
\begin{aligned}
\underset{\W\in\Rtensor}{\text{min}}& &&C\|\W_\Omega - \Y_\Omega\|^2 + \sum_{k=1}^K \dfrac{1}{\lambda_k}\Norm{W_k^{(k)}}^2_* \\
\text{subject to}&\hspace{0.5cm} &&\W \ge 0. \label{eqn:primal_problem_nn}
\end{aligned}
\end{equation}

The fixed-rank dual problem of (<ref>) is

\begin{equation}\label{eqn:outer_manifold_problem_nn}
\underset{U\in S^{n_1}_{r_1}\times\cdots\times S^{n_K}_{r_K}}{\min}\;g(U),
\end{equation}

$g(U)$ is given by

\begin{equation}\label{eqn:inner_strong_convex_problem_nn}
g(U) = \underset{\s\in\Rptensor}{\max}\bigg(\underset{\Z\in\mathcal{C}}{\max}\; \innerproduct{\Z}{\Y_\Omega} - \dfrac{1}{4C}\|\Z\|^2 -\sum_{k=1}^K \frac{\lambda_k}{2} \Norm{U_k^TZ_k + U_k^TS_k}^2\bigg).
\end{equation}
where $U = (U_1,\ldots,U_K)$ and $\Rptensor$ is the set of all tensors of size $n_1 \times \cdots n_K$ with non-negative entries.

The problem (<ref>) can be solved alternatively for $\Z$ and $\s$. Then the resulting problem in $\Z$ is unconstrained least-squares problem, which can be solved using linear conjugate gradient algorithm (For various preconditioned CG approaches, see [25, 5, 6, 7, 8, 9, 10, 14, 20, 24, 13, 15, 16, 17, 18, 19, 21, 22, 23]), and over $\s$ it is non-negative least-squares problem which can be solved using [43]. The cost $g(U)$ and the gradient $\nabla g(U)$ can be easily computed following Lemma <ref>, hence, we employ Riemannian conjugate gradient to solve the outer optimization problem over $U$. The per-iteration complexity of the proposed algorithm is $
O\bigg(T|\Omega|\sum_{k=1}^K{r_k}+\sum_{k=1}^K n_k r_k^2+\sum_{k=1}^K r_k^3\bigg)$, where $T$ is the sum of iterations of CG and NNLS solver for solving $\Z$ and $\s$ respectively.

§.§ Hankel Tensor Completion

The primal problem is given by

∈min          C_Ω-_Ω^2 + ∑_k=1^K 1λ_k_2k(^(k))^2_*,

where $\h_k = \unfold{k}\circ \h$ and $\h$ is the $k$-th order Hankel transform such that for $\W\in \Rtensor$, $\h(\W)\in \R^{\tau_1\times n_1-\tau_1+1\times \cdots \times \tau_K\times n_K-\tau_K+1}$ is a $2K$ order tensor for some $\tau = (\tau_1,\ldots,\tau_K),$ which is the duplication parameter of the Hankel transform (see [45]).

The fixed-rank dual problem of (<ref>) is

\begin{equation}\label{eqn:hankel_outer_manifold_problem}
\underset{U\in S^{n_1-\tau_1+1}_{r_1}\times\cdots\times S^{n_K-\tau_K+1}_{r_K}}{\min}\;g(U),
\end{equation}

where $U = (U_1,\ldots,U_K),$ and

\begin{equation}
\begin{aligned}
\label{eqn:hankel_inner_strong_convex_problem_constrained}
g(U) = \underset{\Z\in\mathcal{C},\s}{\max} \innerproduct{\Z}{\Y_\Omega} - \dfrac{1}{4C}\|\Z\|^2 -\sum_{k=1}^K \frac{\lambda_k}{2} \Norm{U_k^TS_k}^2, \quad \text{subject to} \quad\h^*(\s) = \Z.
\end{aligned}
\end{equation}

The problem (<ref>) can be solved using linear conjugate gradient algorithm over $\Z$ and $\s,$ and the equality constraint is ensured at each step by performing a projection step. The cost $g(U)$, the gradient $\nabla g(U),$ and the directional derivative $D\nabla g(U)[V]$ can be easily computed following Lemma <ref> and Remark <ref> above. We employ Riemannian trust-region algorithm to solve (<ref>). The per-iteration complexity of the proposed algorithm is $
O\bigg(TI|\Omega|\sum_{k=1}^K{r_k}+\sum_{k=1}^K n_k r_k^2+\sum_{k=1}^K r_k^3\bigg)$, where $T$ is the iterations of CG and $I|\Omega|$ is the number of non-zero entries of $\s$.

§ TOY EXPERIMENT

We consider a toy problem of size $100\times 100\times 3$ with 10% train data for both Nonnegative and Hankel Tensor Completion problems. The ranks $(r_1,r_2,r_3)$ are chosen as $(10,10,3)$ for  both the experiments, and the regularization constants $\lambda_k$'s are chosen according to [27]. The maximum iterations for Riemannian optimization problem was set to 200. The duplication parameter $\tau$ in the Hankel transform was set to $(10,10,1)$. The plots of the variation in gradient norm and relative duality gap with iterations for both problems is shown in Appendix <ref> in Figures <ref> and <ref>. We observe that in both the cases the plots of gradient norm, and the relative duality gap decreases rapidly with iterations.

[1]
Boumal, N.: An introduction to optimization on smooth manifolds, 2020. Accessed online from <http://www.nicolasboumal.net/book>

[2]
Kolda, T. G. and Bader, B. W.: Tensor decompositions and applications, SIAM Review, 51(3), pp. 455–500, 2009.

[4]
Absil, P.-A., Mahony, R. and Sepulchre, R.: Optimization Algorithms on Matrix Manifolds, Princeton University Press, Princeton, NJ, 2008.
Accessed online from <https://press.princeton.edu/absil>

[5]
M. Benzi. Preconditioning techniques for large linear systems: A
survey. Journal of Computational Physics, 182(2):418–477, 2002.

[6]
S. Das, S. Katyan, and P. Kumar. Domain de-
composition based preconditioned solver for bundle adjustment. In
R. Venkatesh Babu, Mahadeva Prasanna, and Vinay P. Namboodiri,
editors, Computer Vision, Pattern Recognition, Image Processing, and
Graphics, pages 64–75, Singapore, 2020. Springer Singapore.

[7]
S. Das, S. Katyan, and P. Kumar. A deflation based
fast and robust preconditioner for bundle adjustment. In Proceedings of
the IEEE/CVF Winter Conference on Applications of Computer Vision
(WACV), pages 1782–1789, January 2021.

[8]
S. Katyan, S. Das, and P. Kumar. Two-grid precon-
ditioned solver for bundle adjustment. In 2020 IEEE Winter Conference
on Applications of Computer Vision (WACV), pages 3588–3595, 2020

[9]
K. Mehta, A. Mahajan, P. Kumar Effects of spectral normal-
ization in multi-agent reinforcement learning. In IJCNN, 2023

[10]
P. Kumar. Aggregation based on graph matching and inexact coarse
grid solve for algebraic two grid. International Journal of Computer
Mathematics, 91(5):1061–1081, 2014

[11]
P Kumar, L Grigori, F Nataf, Q Niu, Combinative preconditioning based on relaxed nested factorization and tangential filtering preconditioner, INRIA, HAL Id: inria-00392881, 2009.

[12]
L. Grigori, P. Kumar, F. Nataf, K. Wang, A class of multilevel parallel preconditioning strategies, INRIA 7410, 2010.

[13]
Sahithi Rampalli, Natasha Sehgal, Ishita Bindlish, Tanya Tyagi, Pawan Kumar, Efficient fpga implementation of conjugate gradient methods for laplacian system using hls, arXiv preprint arXiv:1803.03797.

[14]
P. Kumar, L. Grigori, F. Nataf, and Q. Niu. On relaxed nested
factorization and combination preconditioning. International Journal of
Computer Mathematics, 93(1):179–199, 2016

[15]
P. Kumar, L. Grigori, Q. Niu, F. Nataf, Fourier Analysis of Modified Nested
Factorization Preconditioner for Three-Dimensional Isotropic Problems, HAL Id: inria-00448291, 2010.

[16]
A. Aggarwal, S. Kakkar, P. Kumar, Multithreaded Filtering Preconditioner for Diffusion Equation on Structured Grid, arXiv preprint arXiv:1909.09771, 2019.

[17]
P. Kumar, Fast solvers for nonsmooth optimization problems in phase separation, 2015 Federated Conference on Computer Science and Information Systems (FedCSIS), 2015.

[18]
P. Kumar, A class of preconditioning techniques suitable for partial differential equations of structured and unstructured mesh, PhD thesis, 2010.

[19]
P. Kumar, Fast Preconditioned Solver for Truncated Saddle Point Problem in Nonsmooth Cahn–Hilliard Model, Recent Advances in Computational Optimization: Results of the Workshop on Computational Optimization WCO 2015

[20]
P. Kumar, Stefano Markidis, Giovanni Lapenta, Karl Meerbergen,
and Dirk Roose. High performance solvers for implicit particle in cell
simulation. Procedia Computer Science, 18:2251–2258, 2013. 2013
International Conference on Computational Science

[21]
Q Niu, L Grigori, P Kumar, F Nataf, Modified tangential frequency filtering decomposition and its Fourier analysis, Numerische Mathematik 116 (1), 123-148, 2010.

[22]
P Kumar, Purely algebraic domain decomposition methods for the incompressible Navier-Stokes equations, arXiv preprint arXiv:1104.3349, 2011.

[23]
P. Kumar, Multithreaded direction preserving preconditioners, IEEE 13th International Symposium on Parallel and Distributed Computing, 2014.

[24]
P. Kumar, Karl Meerbergen, and Dirk Roose. Multi-threaded nested
filtering factorization preconditioner. PARA, Applied Parallel and Scientific Computing, pages 220–234,
Berlin, Heidelberg, 2013.

[25]
Y. Saad. Iterative Methods for Sparse Linear Systems. Society for
Industrial and Applied Mathematics, Philadelphia, PA, USA, 2 edition,

[26]
Journée, M., Bach, F., Absil, P.-A. and Sepulchre, R.: Low-rank optimization on the cone of positive semidefinite matrices, SIAM Journal on Optimization,  20(5), pp. 2327-2351, 2010.

[27]
Nimishakavi, M., Jawanpuria, P. and Mishra, B.: A Dual Framework for Trace Norm Regularized Low-rank Tensor Completion, Conference on Neural Information Processing Systems (NeurIPS), 2018.

[28]
Argyriou, A., Evgeniou, T. and Pontil, M.: Multi-task feature learning, NIPS, 2006.

[29]
Liu, J., Musialski, P., Wonka, P., and Ye, J.: Tensor completion for estimating missing values in visual data, IEEE transactions on pattern analysis and machine intelligence, 35(1), pp. 208-220, 2012.

[30]
Kressner, D., Steinlechner, M. and Vandereycken, B.: Low-rank tensor completion by Riemannian optimization, BIT Numerical Mathematics, 54(2), pp. 447-468, 2014.

[31]
Kasai, H. and Mishra, B.: Low-rank tensor completion: a Riemannian manifold preconditioning approach, International conference on machine learning (PMLR), pp. 1012-1021, 2016.

[32]
Xu, Y. and Yin, W.: A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion, SIAM Journal on imaging sciences, 6(3), pp. 1758-1789, 2013.

[33]
Signoretto, M., Dinh, Q. T., De Lathauwer, L., and Suykens, J. A.: Learning with tensors: a framework based on convex optimization and spectral regularization, Machine Learning, 94(3), pp. 303-351, 2014.

[34]
Cheng, H., Yu, Y., Zhang, X., Xing, E., and Schuurmans, D.: Scalable and sound low-rank tensor learning, Artificial Intelligence and Statistics (PMLR), pp. 1114-1123, 2016.

[36]
Jawanpuria, P. and Mishra, B.: A unified framework for structured low-rank matrix learning, ICML, 2018.

[37]
Sinha, T.K., Naram, J. and Kumar, P.: Nonnegative Low-Rank Tensor Completion via Dual Formulation with Applications to Image and Video Completion, Winter Conference on Applications of Computer Vision (WACV), 2022.

[38]
Guo, X., Yao, Q. and Kwok, J. T.: Efficient sparse low-rank tensor completion using the frank-wolfe algorithm, AAAI, 2017.

[43]
Kim, D., Sra, S. and Dhillon, I. S.: A non-monotonic method for large-scale non-negative least squares, Optimization Methods and Software, 28(5), pp. 1012–1039, 2013.

[45]
Yokota, T., Erem, B., Guler, S., Warfield, S. K. and Hontani, H.: Missing Slice Recovery for Tensors Using a Low-Rank Model in Embedded Space, IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 8251-8259, 2018.

[46]
Ying, J., Lu, H., Wei, Q., Cai, J.F., Guo, D., Wu, J., Chen, Z., Qu, X.: Hankel Matrix Nuclear Norm Regularized Tensor Completion for $N$-dimensional Exponential Signals, IEEE Transactions on Signal Processing, 65(14), pp. 3702-3717, 2017.

[47]
Wang, X., Wu, Y., Zhuang, D., Sun, L.: Low-Rank Hankel Tensor Completion for Traffic Speed Estimation, arXiv preprint arXiv:2105.11335, 2021.

[48]
Chen, B., Sun, T., Zhou, Z., Zeng, Y. and Cao, L.: Nonnegative tensor completion via low-rank Tucker decomposition: model and algorithm, IEEE Access, 7, pp. 95903-95914, 2019.

[54]
Boumal, N. and Absil, P.-A.: Low-rank matrix completion via preconditioned optimization on the Grassmann manifold, Linear Algebra and its Applications, vol. 475, pp. 200-239, 2015.

§ EXPERIMENTAL PLOTS

[Gradient norm v/s Iters]
[Relative duality gap v/s Iters]
Variation of gradient norm and relative duality gap  with iterations for Nonnegative Tensor Completion problem.

[Gradient norm v/s Iters]
[Relative duality gap v/s Iters]
Variation of gradient norm and relative duality gap  with iterations for Hankel Tensor Completion problem.

§ PROOF OF DUAL FORMULATION (PROOF OF THEOREM <REF>).

Consider the following subproblem of the above problem

\begin{equation}
\begin{aligned}
\underset{\substack{\W^{(k)},\\ k \in \{1,\ldots,K\}}}{\text{min}}& &&C\Norm{\bigg(\sum_{k=1}^K \W^{(k)}\bigg)_\Omega-\Y_\Omega}^2 + \sum_{k=1}^K \frac{1}{2\lambda_k} \innerproduct{\Theta_k^\dagger W^{(k)}_k}{W^{(k)}_k}. \\
\text{subject to}&\hspace{0.5cm} &&A(\W) \ge 0.
\end{aligned}
\end{equation}

We introduce auxiliary variables $U_k$ with the associated constraints $U_k = W_k^{(k)}$. Now, we construct the Lagrangian of the problem

\begin{equation}
\begin{aligned}\label{eqn:lagrangian}
\mathcal{L}(\W^{(1)},\ldots,\W^{(K)},U_1,&\ldots,U_K, \Lambda_1,\ldots, \Lambda_K, s) = C\Norm{\bigg(\sum_{k=1}^K \W^{(k)}\bigg)_\Omega-\Y_\Omega}^2 + \\
&\sum_{k=1}^K \frac{1}{2\lambda_k} \innerproduct{\Theta_k^\dagger U_k}{ U_k} + \sum_{k=1}^K \innerproduct{\Lambda_k}{W^{(k)}_k - U_k} - \innerproduct{s}{A(\W)},
\end{aligned}
\end{equation}
where $\Lambda_1,\ldots, \Lambda_K,s$ are the Lagrange multipliers.

The dual function of the above will be given by

q(Θ_1,…,Θ_K,Λ_1,…, Λ_K,s) =^(k), U_k,
k ∈{1,…,K}min ℒ(^(1),…,^(K),U_1,…,U_K, Λ_1, …, Λ_K, s).

Taking the derivative of the Lagrangian with respect to $\W^{(k)}$ and equating to 0, we get

\begin{align}
\frac{\partial \mathcal{L}}{\partial \W^{(k)}} = 0
&\implies 2C\bigg[\bigg(\sum_{k=1}^K \W^{(k)}\bigg)_\Omega-\Y_\Omega\bigg] + \fold{k} (\Lambda_k) - A^*(s) = 0 \nonumber\\
&\implies \fold{k} (\Lambda_k) = 2C\bigg[\Y_\Omega-\bigg(\sum_{k=1}^K \W^{(k)}\bigg)_\Omega\bigg] + A^{*}(s)\label{eqn:lambda_k},
\end{align}

where $A^*$ is the adjoint of $A$.

The right-hand side of (<ref>) is independent of $k$, so $\fold{k}(\Lambda_k)$ is also independent of $k$. Then we have,

\begin{equation}\label{eqn:lambda_Z}
\fold{k}(\Lambda_k) = \Z + A^*(s),\text{ where }\Z = 2C\bigg[\Y_\Omega-\bigg(\sum_{k=1}^K \W^{(k)}\bigg)_\Omega\bigg].
\end{equation}

Next, we minimize w.r.t. $U_k$, to get

\begin{equation}\label{eqn:U_k}
U_k = \lambda_k\Theta_k \Lambda_k.
\end{equation}

Now we compute (<ref>) term by term by using (<ref>) and (<ref>)

\begin{equation}\label{eqn:term1}
C\Norm{\bigg(\sum_{k=1}^K \W^{(k)}\bigg)_\Omega-\Y_\Omega}^2 = C\Big(\dfrac{1}{4C^2}\|\Z\|^2\Big) = \dfrac{1}{4C}\|\Z\|^2,
\end{equation}

\begin{align}\label{eqn:term2}
\sum_{k=1}^K \frac{1}{2\lambda_k} &\innerproduct{\Theta_k^\dagger U_k}{U_k}-\innerproduct{\Lambda_k}{U_k} = \frac{1}{2\lambda_k} \innerproduct{\Theta_k^\dagger U_k-2\lambda_k\Lambda_k}{U_k} \nonumber\\
&= \sum_{k=1}^K \frac{1}{2\lambda_k} \innerproduct{\Theta_k^\dagger \lambda_k\Theta_k[Z_k + (A^*(s))_k]-2\lambda_k[Z_k + (A^*(s))_k]}{\lambda_k\Theta_k[Z_k + (A^*(s))_k]} \nonumber\\
&= -\sum_{k=1}^K \frac{\lambda_k}{2} \innerproduct{Z_k + (A^*(s))_k}{\Theta_k[Z_k + (A^*(s))_k]},
\end{align}

\begin{align}\label{eqn:term3}
\sum_{k=1}^K \innerproduct{\Lambda_k}{W^{(k)}_k} - \innerproduct{A^*(s)}{\W} &= \sum_{k=1}^K \innerproduct{\fold{k}(\Lambda_k)}{\W^{(k)}} - \innerproduct{A^*(s)}{\W} \nonumber\\
&= \innerproduct{\Z+A^*(s)}{\sum_{k=1}^K \W^{(k)}} - \innerproduct{A^*(s)}{\W}\nonumber\\
&= \innerproduct{\Z+A^*(s)}{\W} - \innerproduct{A^*(s)}{\W} = \innerproduct{\Z}{\W} \overset{(*)}{=} \innerproduct{\Z}{\W_\Omega}\nonumber\\
&= \bigg\langle \Z,\bigg(\sum_{k=1}^K \W^{(k)}\bigg)_\Omega\bigg\rangle \nonumber \\
&= \innerproduct{\Z}{\Y_\Omega-(1/2C)\Z} = \innerproduct{\Z}{\Y_\Omega} - \dfrac{1}{2C}\|\Z\|^2,
\end{align}

where $(*)$ is obtained by noting that $\Z = \Z_\Omega$ by definition in (<ref>). Summing (<ref>), (<ref>), (<ref>), we obtain the expression for the dual function as

\begin{equation}
q(\Theta_1,\ldots,\Theta_K,s) = \innerproduct{\Z}{\Y_\Omega} - \dfrac{1}{4C}\|\Z\|^2 -\sum_{k=1}^K \frac{\lambda_k}{2} \innerproduct{Z_k + (A^*(s))_k}{\Theta_k[Z_k + (A^*(s))_k]}.\nonumber
\end{equation}

This gives the minimax problem (<ref>). From (<ref>) and (<ref>), we can deduce the said relation between optimal points of primal and minimax problems.

§ PROOF OF DUALITY GAP (PROOF OF THEOREM <REF>).

Following [27], we rewrite problem (<ref>) in the form

Θ∈P^n_1 ×⋯×P^n_Kmin          G_1(Θ),  where

G_1(Θ) = ∈𝒞, s∈^nmax            _Ω - 14C^2 -∑_k=1^K λ_k/2 Z_k + (A^*(s))_kΘ_k[Z_k + (A^*(s))_k].

Now, using the min-max interchange, the max-min problem equivalent to <ref> can be written as

∈𝒞, s∈^nmax            G_2(,s),

G_2(,s)    = Θ∈P^n_1 ×⋯×P^n_Kmin _Ω - 14C^2 -∑_k=1^K λ_k/2 Z_k + (A^*(s))_kΘ_k[Z_k + (A^*(s))_k]

=  _Ω - 14C^2 -∑_k=1^K  λ_k/2 Θ_k ∈P^n_kmax Z_k + (A^*(s))_kΘ_k[Z_k + (A^*(s))_k]

=  _Ω - 14C^2 -∑_k=1^K  λ_k/2 B_k(,s),

B_k(,s) = Θ_k ∈P^n_kmax            Z_k + (A^*(s))_kΘ_k[Z_k + (A^*(s))_k].

We know that the value of $B_k(\Z,s)$ is equal to the largest eigenvalue of $(Z_k + (A^*(s))_k)(Z_k + (A^*(s))_k)^T$. Equivalently, if $\sigma_k$ is the largest singular value of $Z_k + (A^*(s))_k$, then $B_k(\Z,s) = \sigma_k^2$. Now, the duality gap will be given by

\begin{equation*}
\Delta = G_1(\hat{\Theta}) - G_2(\hat{\Z},\hat{s}) = \sum_{k=1}^{K} \frac{\lambda_k}{2} \left( \sigma_k^2 -  \|\hat{U}_k^T(\hat{Z}_k+\hat{A}_k)\|^2 \right),
\end{equation*}
where $\hat{U}, \hat{Z},$ and $\hat{s}$ were defined in Theorem <ref>.