# RIFLE: Robust Inference from Low Order Marginals

Sina Baharlouei
<EMAIL_ADDRESS>
&Kelechi Ogudu
<EMAIL_ADDRESS>
&Sze-chuan Suen 11footnotemark: 1
<EMAIL_ADDRESS>
&Meisam Razaviyayn 11footnotemark: 1
<EMAIL_ADDRESS>
Industrial and Systems Engineering, University of Southern, California, Los
Angeles, CA 90089, USADepartment of Computer Science, University of Southern,
California, Los Angeles, CA 90089, USA

###### Abstract

The ubiquity of missing values in real-world datasets poses a challenge for
statistical inference and can prevent similar datasets from being analyzed in
the same study, precluding many existing datasets from being used for new
analyses. While an extensive collection of packages and algorithms have been
developed for data imputation, the overwhelming majority perform poorly if
there are many missing values and low sample size, which are unfortunately
common characteristics in empirical data. Such low-accuracy estimations
adversely affect the performance of downstream statistical models. We develop
a statistical inference framework for predicting the target variable without
imputing missing values. Our framework, RIFLE (Robust InFerence via Low-order
moment Estimations), estimates low-order moments with corresponding confidence
intervals to learn a distributionally robust model. We specialize our
framework to linear regression and normal discriminant analysis, and we
provide convergence and performance guarantees. This framework can also be
adapted to impute missing data. In numerical experiments, we compare RIFLE
with state-of-the-art approaches (including MICE, Amelia, MissForest, KNN-
imputer, MIDA, and Mean Imputer). Our experiments demonstrate that RIFLE
outperforms other benchmark algorithms when the percentage of missing values
is high and/or when the number of data points is relatively small. RIFLE is
publicly available111https://github.com/optimization-for-data-driven-
science/RIFLE.

##  1 Introduction

Machine learning algorithms have shown promise when applied to a wide range of
problems, including healthcare, finance, social data analysis, image
processing, and speech recognition. However, this success mainly relied on the
availability of large-scale, high-quality datasets, which may be scarce in
many practical problems, especially in medical and health applications
(Pedersen et al., ; Sterne et al., 2009; Beaulieu-Jones et al., 2018).
Moreover, in such applications, many experiments and datasets suffer from
small sample size. At the same time, an increasingly large number of datasets
are publicly available. To fully and effectively utilize information on
related research questions from diverse datasets, information across various
datasets needs to be combined in a reliable fashion. After appending these
datasets together, the obtained dataset can contain large blocks of missing
values as they may not share exactly the same features (Figure 1).

For example, similar data may be collected from different health centers or
across waves of a dataset, although questions may not be exactly the same and
some questions be asked while others omitted. Appending two or more of these
datasets may then lead to many missing values in the resulting dataset. Thus,
a pre-processing step should be performed to remove or impute missing values.
However, removing rows containing missing values is not an option when the
percentage of missingness in a dataset is high, or the distribution of missing
values is MNAR (Missing Not At Random). For instance, as it is demonstrated in
Figure 1, if we eliminate the rows that have at least one missing entry, the
entire dataset will be discarded.

Figure 1: Consider the problem of predicting the trait $y$ from feature vector
$(x_{1},\dots,x_{100})$. Suppose that we have access to three data sets: The
first dataset includes the measurements of $(x_{1},x_{2},\dots,x_{40},y)$ for
$n_{1}$ individuals. The second dataset collects data from another $n_{2}$
individuals by measuring $(x_{30},\dots,x_{80})$ with no measurements of the
target variable $y$ in it; and the third dataset contains the measurements
from the variables $(x_{70},\dots,x_{100},y)$ for $n_{3}$ number of
individuals; How one should learn the predictor
$\hat{y}=h(x_{1},\dots,x_{100})$ from these three datasets?

The most common approach for handling missing values in a learning task is to
impute them in a pre-processing stage. The general idea behind data imputation
approaches is that the missing values can be predicted using the other,
available data points and correlated features. The imputation of data allows
practitioners to run standard statistical algorithms requiring complete data.
However, the performance of the prediction model can be highly reliant on the
accuracy of the imputer. High error rates in prediction of missing values by
the imputer can lead to catastrophic performance of the statistical methods
executed on the imputed data.

An alternative approach is to learn the parameters of the model directly
(without imputation) based on the available data. In this paper, we propose
RIFLE (Robust InFerence via Low-order moment Estimations) for the direct
inference of a target variable based on a set of features containing missing
values. The proposed framework does not require the data to be imputed in a
pre-processing stage. However, if needed, it can be used as a pre-processing
tool for imputing data as well. The main idea of the proposed framework is to
estimate the first and second-order moments of the data and their confidence
intervals by bootstrapping on the available data matrix entries. Then, RIFLE
finds the optimal parameters of the statistical model for the worst-case
distribution with the low-order moments (mean and variance) within the
estimated confidence intervals (See Figure 2).

Figure 2: Prediction of the target variable without imputation. RIFLE
estimates confidence intervals for low-order (first and second order)
marginals from the input data containing missing values. Then, it solves a
distributionally robust problem over all distributions whose low-order
marginals are within the estimated confidence intervals.

###  1.1 Related Work

The fundamental idea behind many data imputation approaches is that the
missing values can be predicted based on the available data of other data
points and correlated features. One of the most straightforward imputation
techniques is to replace missing values by the mean (or median) of that
feature calculated from what data is available see (Little and Rubin, 2019,
Chapter 3). However, this naïve approach ignores the correlation between
features and does not preserve the variance of features. Another class of
imputers have been developed based on the least-square methods (Raghunathan et
al., 2001; Kim et al., 2005; Zhang et al., 2008; Cai et al., 2006).
Raghunathan et al. (2001) learns a linear model with multivariate Gaussian
noise for the feature with the least amount of missing values. It repeats the
same procedure on the updated data to impute the next feature with the least
amount of missing values until all features are completely imputed. One
drawback of this approach is that the error from the imputation of previous
features can be propagated to subsequent features. To impute entries of a
given feature in a dataset, Kim et al. (2005) learns several univariate
regression models that consider that feature as the response, and then it
takes the average of these predictions as the final value of imputation. This
approach fails to learn the correlations involving more than two features.

Many more complex algorithms have been developed for imputation, although many
are sensitive to initial assumptions and may not converge. For instance, KNN-
Imputer imputes a missing feature of a data point by taking the mean value of
the $K$ closest complete data points (Troyanskaya et al., 2001). MissForest,
on the other hand, imputes the missing values of each feature by learning a
random forest classifier using other features of the training data (Stekhoven
and Bühlmann, 2012). MissForest does not need to assume that all features are
continuous (Honaker et al., 2011), or all categorical (Schafer, 1997).
However, both KNN-imputer and MissForest do not provide statistical or
computational convergence guarantees for their algorithm. Moreover, when the
proportion of missing values is high, both are likely to have a severe drop in
performance as it is demonstrated in Section 5. The Expectation Maximization
(EM) algorithm is another popular approach that learns the parameters of a
prior distribution on the data using available values based on the EM
algorithm of (Dempster et al., 1977); see also Ghahramani and Jordan (1994)
and Honaker et al. (2011). The EM algorithm is also used in Amelia, which fits
a jointly normal distribution to the data using EM and bootstrap technique
(Honaker et al., 2011). While Amelia demonstrates a superior performance on
datasets following a normal distribution, it is highly sensitive to the
violation of the normality assumption (as discussed in Bertsimas et al.
(2017)). Ghahramani and Jordan (1994) adopt the EM algorithm to learn a joint
Bernoulli distribution for the categorical data and a joint Gaussian
distribution for the continuous variables independently. While those
algorithms can be viewed as inference methods based on low-order estimates of
moments, they do not consider uncertainty in such low-order moments estimates.
By contrast, our framework utilizes robust optimization to consider the
uncertainty around the estimated moments. Moreover, our optimization procedure
for imputation and prediction is guaranteed to converge unlike some of the
aforementioned algorithms.

Another popular method for data imputation is multiple imputations by chained
equations (MICE). MICE learns a parametric distribution for each feature
conditional on the remaining features. For instance, it assumes that the
current target variable is a linear function of other features with a zero-
mean Gaussian noise. Each feature can have its own distinct distribution and
parameters (e.g., Poisson regression, logistic regression, etc). Based on the
learned parameters of conditional distributions, MICE can generate one or more
imputed datasets (Buuren and Groothuis-Oudshoorn, 2010). More recently,
several neural network-based imputers have been proposed. GAIN (Generative
Adversarial Imputation Network) adopts the idea of learning a generative
adversarial network based on the available data, and then imputing the missing
values using the trained generator (Yoon et al., 2018). One advantage of GAIN
over other existing GAN imputers is that it does not need a complete dataset
during the training phase. MIDA (Multiple Imputation using Denoising
Autoencoders) is an auto-encoder based approach that trains a denoising auto-
encoder on the available data considering the missing entries as noise.
Similar to other neural network-based methods, these algorithms suffer from
their black-box nature. They are challenging to interpret/explain, which makes
them not popular in mission-critical approaches such as healthcare. In
addition, no statistical or computational guarantees are provided for these
algorithms.

Bertsimas et al. (2017) formulates the imputation task as a constrained
optimization problem where the constraints are determined by the underlying
classification model such as KNN ($k$-nearest neighbors), SVM (Support Vector
Machine), and Decision Trees. Their general framework is non-convex, and the
authors relax the optimization for each choice of the cost function using
first-order methods. The relaxed problem is then optimized by the block
coordinate descent algorithm. They show the convergence and accuracy of their
proposed algorithm numerically, while a theoretical analysis which guarantees
the convergence of the algorithm is absent in their work.

Beside these imputation methods, Shivaswamy et al. (2006) and Xu et al. (2009)
adopt robust optimization approaches to learn the parameters of a support
vector machine model. They consider uncertainty sets for the missing entries
in the dataset, and solve a min-max problem over that set. The obtained
classifiers are robust to the uncertainty of missing entries within that
region. In contrast to the imputation based approaches, the robust
classification formulation does not carry the imputation error to the
classification phase. However, finding appropriate intervals for each missing
entry is a challenging problem by itself, and it might not be available
apriori in many real datasets. Moreover, their proposed algorithms are limited
to the SVM classifier. On the other hand, as we will see in this work,
estimating uncertainty regions for estimated low-order marginals can be done
naturally. Furthermore, our framework is not restricted to any particular
machine learning module.

###  1.2 Contributions

We list our main contributions in this report:

  1. 1.

We present a general robust optimization framework for learning the parameters
of a learning model in the presence of missing values. The proposed framework
does not require data imputation as a pre-processing stage. We specialize the
framework to linear regression and normal discriminant analysis models in two
case studies. The proposed framework provides a novel strategy for inference
in the presence of missing values, especially for datasets with large
proportions of missing values.

  2. 2.

We provide theoretical convergence guarantees and analyze the iteration
complexity of the presented algorithms for robust formulations of linear
regression and normal discriminant analysis. Moreover, we analyze the
asymptotic statistical properties of the solutions found by the algorithms.

  3. 3.

While the robust inference framework is primarily designed for direct
statistical inference in the presence of missing values without performing
data imputation, it can be adopted as an imputation tool as well. To
demonstrate the quality of the proposed imputer, we compare its performance
with several widely used imputation packages such as MICE, Amelia, MissForest,
KNN-Imputer and MIDA on real and synthetic datasets. Generally speaking, our
method outperforms all of the mentioned packages when the number of missing
entries is large.

  4. 4.

Based on our developments in this work, we created a publicly available
package for inference in the presence of missing data. Our implementation and
the code is available at https://github.com/optimization-for-data-driven-
science/RIFLE.

###  1.3 Organization of the Paper

In Section 2, we present the general robust inference framework in the
presence of missing data used in this work. In Section 3, we specialize the
proposed framework for the robust linear regression problem. The algorithms in
Section 3 are analyzed from both optimization and statistical standpoints.
Section 4 is dedicated to the robust classification task. In particular, we
demonstrate how to derive an algorithm for the robust normal discriminant
analysis as a case study of robust classification framework. Finally, in
Section 5, we provide numerical experiments to evaluate the performance of the
proposed algorithms.

##  2 Robust Inference via Estimating Low-order Moments

RIFLE is a distributionally robust optimization (DRO) framework based on
estimated low-order marginals. Assume that
$({\mathbf{x}},y)\in\mathbb{R}^{d}\times\mathbb{R}$ follows a joint
probability distribution $P^{*}$. A standard objective in parametric
statistical learning is to find the parameter ${\bm{\theta}}$ that minimizes
the population risk with respect to a given loss function $\ell$:

$\min_{\theta}\quad\mathbb{E}_{({\mathbf{x}},y)\sim
P^{*}}\Big{[}\ell\Big{(}{\mathbf{x}},y;\bm{\theta}\Big{)}\Big{]}.$ (1)

Since the underlying distribution of data is not usually available in
practice, the above problem cannot be solved directly. The most common
approach for approximating the population risk problem (1) is to minimize the
empirical risk with respect to $n$ given i.i.d samples
$({\mathbf{x}}_{1},y_{1}),\dots,({\mathbf{x}}_{n},y_{n})$ drawn from the joint
distribution $P^{*}$:

$\min_{\bm{\theta}}\quad\frac{1}{n}\sum_{i=1}^{n}\ell({\mathbf{x}}_{i},y_{i};\bm{\theta}).$

In the presence of missing entries in the data points
$\\{({\mathbf{x}}_{1},y_{1}),\dots,({\mathbf{x}}_{n},y_{n})\\}$, the above
formulation cannot be used. A natural alternative approach in the presence of
uncertainties and missing data is training via robust optimization. For
example, (Shivaswamy et al., 2006; Xu et al., 2009) suggest the formulation:

$\min_{{\bm{\theta}}}\;\;\max_{\\{\delta_{i}\in\mathcal{N}_{i}\\}_{i=1}^{n}}\;\;\frac{1}{n}\sum_{i=1}^{n}\ell({\mathbf{x}}_{i}-\delta_{i},y_{i};{\bm{\theta}}),$
(2)

where $\mathcal{N}_{i}$ represents the uncertainty region of data point $i$.
(Shivaswamy et al., 2006) obtains such an uncertainty set by assuming a known
distribution on the missing entries of data point $i$. (Xu et al., 2009)
defines $\mathcal{N}_{i}$ as a “box” constraint around the missing entries of
data point $i$. However, in practice, it is difficult to find such “per data
point” constraint sets. In this work, we instead take distributionally robust
approach to estimate the uncertainty set for the data distribution. In
particular, we aim to fit the best parameters of a statistical learning model
for the worst distribution in a given uncertainty set by solving:

$\quad\underset{{\bm{\theta}}}{\min}\quad\underset{P\in\mathcal{P}}{\max}\quad\mathbb{E}_{({\mathbf{x}},y)\sim
P}[\ell({\mathbf{x}},y;{\bm{\theta}})],$ (3)

where $\mathcal{P}$ is an uncertainty set over the underlying distribution of
data. A key observation is that defining the uncertainty set $\mathcal{P}$ in
(3) is easier (and can be done more naturally) than defining the uncertainty
sets $\\{\mathcal{N}_{i}\\}_{i=1}^{n}$ in (2). In particular, the uncertainty
set $\mathcal{P}$ can be obtained naturally by estimating low-order moments of
data distribution using available data. To explain this idea and to simplify
the notations, let ${\mathbf{z}}=({\mathbf{x}},y)$,
$\bar{{\bm{\mu}}}^{{\mathbf{z}}}\triangleq\mathbb{E}[{\mathbf{z}}]$, and
$\bar{{\mathbf{C}}}^{{\mathbf{z}}}\triangleq\mathbb{E}[{\mathbf{z}}{\mathbf{z}}^{T}].$
While $\bar{{\bm{\mu}}}^{{\mathbf{z}}}$ and
$\bar{{\mathbf{C}}}^{{\mathbf{z}}}$ are typically not known exactly, one can
estimate them (within certain confidence intervals) from the available data
despite having missing data points (see Subsection 2.1). In particular, we can
estimate the values
${\bm{\mu}}_{\min}^{{\mathbf{z}}},{\bm{\mu}}_{\max}^{{\mathbf{z}}},{\mathbf{C}}^{{\mathbf{z}}}_{\min}$,
and ${\mathbf{C}}^{{\mathbf{z}}}_{\max}$ from data such that
${\bm{\mu}}^{{\mathbf{z}}}_{\min}\leq\bar{{\bm{\mu}}}^{{\mathbf{z}}}\leq{\bm{\mu}}^{{\mathbf{z}}}_{\max}$
and
${\mathbf{C}}^{{\mathbf{z}}}_{\min}\leq\bar{{\mathbf{C}}}^{{\mathbf{z}}}\leq{\mathbf{C}}_{\max}^{{\mathbf{z}}}$
with high probability (where the inequalities for matrices and vectores denote
component-wise relations). Given these estimated confidence intervals from
data, problem (3) can be reformulated as

$\begin{split}\min_{{\bm{\theta}}}\quad\max_{P}\;\;&\mathbb{E}_{P}[\ell({\mathbf{z}};{\bm{\theta}})]\\\
{\rm
s.t.}\quad&{\bm{\mu}}_{\min}^{{\mathbf{z}}}\leq\mathbb{E}_{P}[{\mathbf{z}}]\leq{\bm{\mu}}_{\max}^{{\mathbf{z}}},\\\
&{\mathbf{C}}_{\min}^{{\mathbf{z}}}\leq\mathbb{E}_{P}[{\mathbf{z}}{\mathbf{z}}^{T}]\leq{\mathbf{C}}_{\max}^{{\mathbf{z}}}.\end{split}$
(4)

In Section 3 and Section 4, we simplify this formulation and develop efficient
algorithms for special popular forms of the loss function $\ell(\cdot;\cdot)$.
However, before proceeding further, we will first explain how we estimate the
parameters ${\bm{\mu}}_{\min},{\bm{\mu}}_{\max},{\mathbf{C}}_{\min}$, and
${\mathbf{C}}_{\max}$ from data in the next subsection.

###  2.1 Estimating Confidence Intervals of Low-order Moments

In this section, we explain the methodology of estimating confidence intervals
for $\mathbb{E}[{\mathbf{z}}_{i}]$ and
$\mathbb{E}[{\mathbf{z}}_{i}{\mathbf{z}}_{j}]$. Let
$({\mathbf{x}}^{1},y^{1}),({\mathbf{x}}^{2},y^{2}),\dots,({\mathbf{x}}^{n},y^{n})$
be $n$ given training data points, where
${\mathbf{x}}^{i}\in\tilde{\mathbb{R}}^{n},y^{i}\in\tilde{\mathbb{R}}$. We
define $\tilde{\mathbb{R}}=\mathbb{R}\cup\\{*\\}$ where the $*$ symbol
represents a missing entry. Assume that ${\mathbf{X}}$ is the data matrix
whose i-th row is ${\mathbf{x}}^{i}$, and ${\mathbf{y}}$ is the vector of
labels. Moreover, assume that ${\mathbf{a}}_{i}$ represents the $i$-th column
(feature) of matrix ${\mathbf{X}}$. We define:

$\tilde{{\mathbf{a}}}_{i}(k)=\begin{cases}{\mathbf{a}}_{i}(k)&\text{if}\>{\mathbf{a}}_{i}(k)\neq\text{*}\\\
0&\text{if}\>{\mathbf{a}}_{i}(k)=\text{*}\end{cases}$

Thus, $\tilde{{\mathbf{a}}}$ is obtained by replacing the missing values with
$0$. We estimate the confidence intervals for mean and covariance of features
using multiple bootstrap samples on the available data. Let
${\mathbf{C}}_{0}[i][j]$ and $\bm{\Delta}_{0}[i][j]$ be the center and the
radius of the confidence interval for ${\mathbf{C}}[i][j]$, respectively. We
compute the center of the confidence interval for ${\mathbf{C}}[i][j]$ as
follows:

${\mathbf{C}}_{0}[i][j]=\frac{1}{m_{ij}}\tilde{{\mathbf{a}}}_{i}^{T}\tilde{{\mathbf{a}}}_{j}$
(5)

where $m_{i}=|\\{k:{\mathbf{a}}_{i}(k)\neq*\\}|$ and
$m_{ij}=|\\{k:{\mathbf{a}}_{i}(k)\neq*,{\mathbf{a}}_{j}(k)\neq*\\}|$. This
estimator is obtained from the rows where both features are available.

To find the radius of confidence intervals for each given pair $(i,j)$ of
features, we choose $k$ different bootstrap samples with length $n$ on the
rows where both features $i$ and $j$ are available. Then, we compute
${\mathbf{C}}_{0}[i][j]$ of two features in each bootstrap sample. The
standard deviation of these estimations determines the radius of the
corresponding confidence interval. Algorithm 1 summarizes the required steps
for computing the radius of confidence interval for the $ij$-th entry of
covariance matrix $\Delta$. Note that the confidence intervals for $\bm{\mu}$
can be computed in the same manner.

Algorithm 1 Estimating Confidence Interval Length $\bm{\Delta}_{ij}$ for
Feature $i$ and Feature $j$.

1:Input: $k:$ Number of estimations

2:for $t=1,\ldots,k$ do

3: Pick $n$ samples with replacement from the rows where both $i$-th and
$j$-th are available.

4: Let $(\hat{X}_{i1},\hat{X}_{j1}),\dots,(\hat{X}_{in},\hat{X}_{jn})$ be the
$i$-th and $j$-th features of the selected samples

5: $C_{t}=\frac{1}{n}\sum_{r=1}^{n}\hat{X}_{ir}$

6:$\bm{\Delta}_{ij}=\textrm{std}(C_{1},C_{2},\dots,C_{K})$

Having ${\mathbf{C}}_{0}$ and $\bm{\Delta}$, the confidence interval for the
matrix ${\mathbf{C}}$ is computed as follows:

$\displaystyle{\mathbf{C}}_{min}={\mathbf{C}}_{0}-c\bm{\Delta}$ (6)
$\displaystyle{\mathbf{C}}_{max}={\mathbf{C}}_{0}+c\bm{\Delta},$ (7)

where the constant $c$ is a hyper-parameter determining the level of
robustness. In Section 3 and Section (4) we specialize the general
distributionally robust problem described in (4) to the linear regression and
normal discriminant analysis models.

###### Remark 1

Since the computation of confidence intervals for different entries of the
covariance matrix are independent of each other, they can be computed in
parallel. In particular, if $\gamma$ cores are available, $\lceil
d/\gamma\rceil$ features (columns of covariance matrix) can be assigned to
each one of the available cores.

##  3 Robust Linear Regression in the Presence of Missing Values

In a linear regression model, we aim to minimize the following objective
function:

$\quad\underset{{\bm{\theta}}}{\min}\quad\|{\mathbf{X}}{\bm{\theta}}-{\mathbf{y}}\|_{2}^{2}$

The optimal ${\bm{\theta}}$ for the linear regression model can be written as

${\bm{\theta}}=({\mathbf{X}}^{T}{\mathbf{X}})^{-1}({\mathbf{X}}^{T}{\mathbf{y}}).$
(8)

Thus, for obtaining the optimal solution, we only need to compute the second-
order moments of the data ${\mathbf{C}}={\mathbf{X}}^{T}{\mathbf{X}}$ and
${\mathbf{b}}={\mathbf{X}}^{T}{\mathbf{y}}$. In other words, It suffices to
compute the inner product of any two column vectors ${\mathbf{a}}_{i}$,
${\mathbf{a}}_{j}$ of ${\mathbf{X}}$, and the inner product of any column
vector ${\mathbf{a}}_{i}$ of ${\mathbf{X}}$ with vector ${\mathbf{y}}$. Since
the matrix ${\mathbf{X}}$ and vector ${\mathbf{y}}$ are not fully observed due
to the existence of missing values, use the same approach as (5) to compute
the point estimators ${\mathbf{C}}_{0}$ and ${\mathbf{b}}_{0}$. These point
estimators can be highly inaccurate, especially when the number of non-missing
rows for two given columns is small. In addition, if the pattern of missing
entries does not follow the MCAR assumption, the point estimators are not
unbiased estimators of ${\mathbf{C}}$ and ${\mathbf{b}}$. To learn a linear
regression model which is robust against large proportion of missing values,
and MNAR missing patterns, we specialize the framework presented in Section 2
to the linear regression model.

###  3.1 A Distributionally Robust Formulation of Linear Regression

As we mentioned above, for solving the linear regression problem, we only need
to estimate the second-order moments of the data
(${\mathbf{X}}^{T}{\mathbf{X}}$ and ${\mathbf{X}}^{T}{\mathbf{y}}$). Thus, the
distributionally robust formulation described in (4) is equivalent to the
following optimization problem for the linear regression model:

$\begin{array}[]{lll}\displaystyle{\min_{{\bm{\theta}}}\max_{{\mathbf{C}},{\mathbf{b}}}}&\displaystyle{{\bm{\theta}}^{T}{\mathbf{C}}{\bm{\theta}}-2{\mathbf{b}}^{T}{\bm{\theta}}+\lambda\|{\bm{\theta}}\|_{2}^{2}}\\\
{\rm
s.t.}&{\mathbf{C}}_{0}-c\bm{\Delta}\leq{\mathbf{C}}\leq{\mathbf{C}}_{0}+c\bm{\Delta},\\\
&{\mathbf{b}}_{0}-c\bm{\delta}\leq{\mathbf{b}}_{0}+c\bm{\delta},\\\
&{\mathbf{C}}\succeq 0,\end{array}$ (9)

where the last constraint guarantees that the covariance matrix is positive
semi-definite. To compute confidence intervals in the above formulation we
specialize Algorithm 1 to the linear regression problem. Let ${\mathbf{M}}$ be
the mask of the input data matrix ${\mathbf{X}}$ defined as:

${\mathbf{M}}_{ij}=\begin{cases}0,&\textrm{if }{\mathbf{X}}_{ij}\textrm{ is
missing,}\\\ 1,&\textrm{otherwise}.\end{cases}$

Assume that $m_{ij}=({\mathbf{M}}^{T}{\mathbf{M}})_{ij}$, which is the number
of rows in the dataset where both features $i$ and $j$ are available. To
estimate the confidence intervals for ${\mathbf{C}}_{ij}$, we select multiple
($k$) samples of size $N=m_{ij}$ from the rows where both features are
available. Each one of these samples with size $m_{ij}$ are obtained by
applying a bootstrap sampler (sampling with replacement) on the $m_{ij}$ rows
where both features are available. Then, for each sample we compute the
second-order moment of two features. Let ${\mathbf{S}}^{2}_{ij}$ be the
variance of the these $k$ estimations. We set
$\bm{\Delta_{ij}}=c{\mathbf{S}}_{ij}$, where $c$ is a hyper-parameter
determining the level of robustness. A larger $c$ corresponds to bigger
confidence intervals and thus a more robust estimator. However, $c$ leads to
very large confidence intervals that can adversely affect the performance of
the trained model. Algorithm 2 summarizes this procedure of computing
confidence intervals.

Algorithm 2 Estimating Confidence Interval Length $\bm{\Delta}_{ij}$ for
Feature $i$ and Feature $j$.

1:Input: $k:$ Number of estimations, feature indices $i$ and $j$.

2:for $t=1,\ldots,k$ do

3: Pick $m_{ij}$ samples with replacement from the rows where both $i$-th and
$j$-th are available.

4: Let $(\hat{X}_{i1},\hat{X}_{j1}),\dots,(\hat{X}_{in},\hat{X}_{jn})$ be the
$i$-th and $j$-th features of the selected samples

5: $C_{t}=\frac{1}{m_{ij}}\sum_{r=1}^{m_{ij}}\hat{X}_{ir}\hat{X}_{jr}$

6:$\bm{\Delta_{ij}}=\textrm{std}(C_{1},C_{2},\dots,C_{K})$

Fixing ${\mathbf{b}}$ and ${\mathbf{C}}$, the minimization problem in (9) has
a closed form solution with respect to $\theta$. To solve the inner
maximization problem, we apply the projected gradient ascent to ${\mathbf{C}}$
and ${\mathbf{b}}$. Note that the projection to the box constraint
${\mathbf{C}}_{0}-c\bm{\Delta}\leq{\mathbf{C}}\leq{\mathbf{C}}_{0}+c\bm{\Delta}$
can be done entriwise and has the following closed form:

$\Pi_{\Delta}({\mathbf{C}}[i][j])=\begin{cases}{\mathbf{C}}[i][j]&\textrm{if
}\quad{\mathbf{C}}_{0}[i][j]-c\bm{\Delta}[i][j]\leq{\mathbf{C}}[i][j]\leq{\mathbf{C}}_{0}[i][j]+c\bm{\Delta}[i][j],\\\
{\mathbf{C}}_{0}[i][j]-c\bm{\Delta}[i][j]&\textrm{if
}\quad{\mathbf{C}}[i][j]\leq{\mathbf{C}}_{0}[i][j]-c\bm{\Delta}[i][j],\\\
{\mathbf{C}}[i][j]+c\bm{\Delta}[i][j]&\textrm{if
}\quad{\mathbf{C}}_{0}[i][j]+c\bm{\Delta}[i][j]\leq{\mathbf{C}}[i][j].\end{cases}$

However, the projection to the set of positive semi-definite matrices requires
a singular value decomposition (SVD) of ${\mathbf{C}}$ at each iteration.
Moreover, satisfying both constraints simultaneously is a challenging problem
and cannot be done in closed-form. To avoid this complexity and time-consuming
singular value decomposition at each iteration, we relax the problem by
removing the PSD constraint on ${\mathbf{C}}$. This relaxation does not
drastically change the performance of the algorithm as our experiments show in
Section 5. In Section 3.4 we propose an algorithm for solving the original
problem efficiently. Algorithm 3, demonstrates the steps of optimizing problem
(9). At each iteration of the algorithm, we first perform one step of
projected gradient ascent on matrix ${\mathbf{C}}$ and ${\mathbf{b}}$, then we
solve the minimization problem with respect to $\theta$ for the updated
${\mathbf{C}}$ and ${\mathbf{b}}$ exactly.

Algorithm 3 Robust Linear Regression

1:Input: ${\mathbf{C}}_{0},{\mathbf{b}}_{0},\bm{\Delta},\bm{\delta},T$

2:Initialize: ${\mathbf{C}}={\mathbf{C}}_{0},{\mathbf{b}}={\mathbf{b}}_{0}$.

3:for $i=1,\ldots,T$ do

4: Update
${\mathbf{C}}=\Pi_{\bm{\Delta}}({\mathbf{C}}+\alpha{\bm{\theta}}{\bm{\theta}}^{T})$

5: Update ${\mathbf{b}}=\Pi_{\bm{\delta}}({\mathbf{b}}-2\alpha{\bm{\theta}})$

6: Set ${\bm{\theta}}=({\mathbf{C}}+\lambda{\mathbf{I}})^{-1}{\mathbf{b}}$

Since, the objective function is convex in ${\bm{\theta}}$ and concave in
${\mathbf{b}}$ and ${\mathbf{C}}$, the minimization and maximization sub-
problems are interchangeable (Sion et al., 1958). Thus, we can equivalently
rewrite Problem (9) as:

$\begin{array}[]{lll}\displaystyle{\max_{{\mathbf{C}},{\mathbf{b}}}}&g({\mathbf{C}},{\mathbf{b}})\\\
{\rm
s.t.}&{\mathbf{C}}_{0}-c\bm{\Delta}\leq{\mathbf{C}}_{0}\leq{\mathbf{C}}_{0}+c\bm{\Delta},\\\
&{\mathbf{b}}_{0}-c\bm{\delta}_{0}\leq{\mathbf{b}}_{0}\leq{\mathbf{b}}_{0}+c\bm{\delta},\\\
&{\mathbf{C}}\succeq 0,\end{array}$ (10)

where
$g({\mathbf{b}},{\mathbf{C}})=\min_{{\bm{\theta}}}{\bm{\theta}}^{T}{\mathbf{C}}{\bm{\theta}}-2{\mathbf{b}}^{T}{\bm{\theta}}+\lambda\|{\bm{\theta}}\|^{2}$.
Algorithm 3 describes the procedure of applying projected gradient ascent on
function $g$ to find an optimal solution of Problem (10). We use
${\mathbf{C}}_{0}$ and ${\mathbf{b}}_{0}$ obtained from Equation (5) as the
initialization points for ${\mathbf{C}}$ and ${\mathbf{b}}$.

###  3.2 Convergence Guarantee of Robust Linear Regression Algorithm

Since the minimization problem in definition of function $g$ is strongly
convex with respect to ${\bm{\theta}}$, the minimizer is unique and can be
found in closed form. Moreover, function $g$ is concave with respect to
${\mathbf{b}}$ and ${\mathbf{C}}$. Thus, applying the projected gradient
ascent algorithm on $g$ with appropriate choice of step-size will converge to
the global maximizer of function $g$. Note that, the convergence rate of
algorithm 3 depends on the Lipschitz constant of function $g$. We use Lemma 1
in (Barazandeh and Razaviyayn, 2020) to compute the Lipschitz constant of
function $g$. Without loss of generality, by rescaling of the data, if it is
necessary, we can assume that $||{\bm{\theta}}||_{2}\leq 1$. Having the
Lipschitz constant of function $g$, Theorem 3 establishes the convergence rate
of Algorithm 3 to an $\epsilon$-optimal saddle point of Problem (10).

###### Lemma 2

The Lipschitz constant of the gradient of the function $g$ defined in Problem
(10) is equal to $L=\frac{2}{\lambda}$.

Proof  The proof is relegated to Appendix B.

###### Theorem 3

Assume that $D=\|{\mathbf{C}}_{0}-{\mathbf{C}}^{*}\|_{F}$. Then Algorithm 3
computes an $\epsilon$-stationary solution of the objective function in (10)
in ${\cal O}(\frac{LD^{2}}{\epsilon})$ iterations, where
$L=\frac{2}{\lambda}$.

Proof  The proof is relegated to Appendix B.

###  3.3 Obtaining Faster Convergence and Acceleration

The convergence rate of Algorithm 3 to a saddle point of optimization problem
(9) can be slow in practice, since the algorithm requires to do a matrix
inversion for updating ${\bm{\theta}}$ and applying the box constraint to
${\mathbf{C}}$ and ${\mathbf{b}}$ at each iteration. While we update the
minimization problem in closed-form with respect to ${\bm{\theta}}$, we can
speed up the convergence rate of the maximization problem by applying
Nesterov’s acceleration method to function $g({\mathbf{b}},{\mathbf{C}})$ in
(10). Note that, we can compute the gradient of function $g$ with respect to
${\mathbf{C}}$ and ${\mathbf{b}}$ using Danskin’s theorem.

Algorithm 4 describes the steps to optimize Problem (10) using Nesterov’s
acceleration method.

Algorithm 4 Applying the Nesterov’s Acceleration Method to Robust Linear
Regression

1:${\mathbf{C}}_{0},{\mathbf{b}}_{0},\bm{\Delta},\bm{\delta},T$

2:Initialize:
${\mathbf{C}}_{1}={\mathbf{C}}_{0},{\mathbf{b}}_{1}={\mathbf{b}}_{0},\gamma_{0}=0,\gamma_{1}=1$.

3:for $i=1,\ldots,T$ do

4: $\gamma_{i+1}=\frac{1+\sqrt{1+4\gamma_{i}^{2}}}{2}$

5:
$Y_{{\mathbf{C}}_{i}}={\mathbf{C}}_{i}+\frac{\gamma_{i}-1}{\gamma_{i+1}}({\mathbf{C}}_{i}-{\mathbf{C}}_{i-1})$

6:
${\mathbf{C}}_{i+1}=[\Pi_{\Delta}(Y_{{\mathbf{C}}_{i}}+\frac{1}{L}\theta\theta^{T})]_{+}$

7:
$Y_{{\mathbf{b}}_{i}}={\mathbf{b}}_{i}+\frac{\gamma_{i}-1}{\gamma_{i+1}}({\mathbf{b}}_{i}-{\mathbf{b}}_{i-1})$

8:
${\mathbf{b}}_{i+1}=[\Pi_{\bm{\Delta}}(Y_{{\mathbf{b}}_{i}}-\frac{2{\bm{\theta}}}{L})]_{+}$

9: Set ${\bm{\theta}}=({\mathbf{C}}_{i+1}+\lambda I)^{-1}{\mathbf{b}}_{i+1}$

###### Theorem 4

Assume that $D=|\hat{\mathbf{C}}-\mathbf{C}^{*}|_{F}$. Then Algorithm 4
computes an $\epsilon$-stationary solution of the objective function in (10)
in ${\cal O}\Big{(}\sqrt{\frac{LD^{2}}{\epsilon}}\Big{)}$ iterations, where
$L=\frac{2}{\lambda}$.

Proof  The proof is relegated to Appendix B.

###  3.4 Solving the Dual Problem of Robust Linear Regression via ADMM

Alternating Direction Method of Multipliers (ADMM) is a popular algorithm for
efficiently solving linearly constrained optimization problems (Gabay and
Mercier, 1976; Hong et al., 2016). It has been extensively applied to large-
scale optimization problems in machine learning and statistical inference in
recent years (Assländer et al., 2018; Zhang et al., 2018). Consider the
following optimization problem consisting of two blocks of variables
${\mathbf{x}}$ and ${\mathbf{y}}$ that are linearly coupled:

$\begin{array}[]{lll}\displaystyle{\min_{{\mathbf{w}},{\mathbf{z}}}}&\displaystyle{f({\mathbf{w}})+g({\mathbf{z}})}\\\
{\rm
s.t.}&{\mathbf{A}}{\mathbf{w}}+{\mathbf{B}}{\mathbf{z}}={\mathbf{c}},\end{array}$
(11)

The augmented Lagrangian of the above problem can be written as:

$\begin{array}[]{lll}\displaystyle{\min_{{\mathbf{w}},{\mathbf{z}}}}&\displaystyle{f({\mathbf{w}})+g({\mathbf{z}})+\langle{\mathbf{A}}{\mathbf{w}}+{\mathbf{B}}{\mathbf{z}}-{\mathbf{c}},\lambda\rangle+\frac{\rho}{2}\|{\mathbf{A}}{\mathbf{w}}+{\mathbf{B}}{\mathbf{z}}-{\mathbf{c}}\|^{2}}\\\
\end{array}$ (12)

ADMM schema updates the primal and dual variables iteratively as it is
depicted in Algorithm 5.

Algorithm 5 Schema of ADMM algorithm Applied to the Augmented Lagrangian
Problem

1:for $t=1,\ldots,T$ do

2:
${\mathbf{w}}^{t+1}=\operatorname*{argmin}_{{\mathbf{w}}}f({\mathbf{w}})+\langle
A{\mathbf{w}}+B{\mathbf{z}}^{t}-{\mathbf{c}},\lambda\rangle+\frac{\rho}{2}\|A{\mathbf{w}}+B{\mathbf{z}}^{t}-{\mathbf{c}}\|^{2}$

3:
${\mathbf{z}}^{t+1}=\operatorname*{argmin}_{{\mathbf{z}}}f({\mathbf{w}}^{t+1})+\langle{\mathbf{A}}{\mathbf{w}}^{t+1}+{\mathbf{B}}{\mathbf{z}}-{\mathbf{c}},\lambda\rangle+\frac{\rho}{2}\|{\mathbf{A}}{\mathbf{w}}^{t+1}+{\mathbf{B}}{\mathbf{z}}-{\mathbf{c}}\|^{2}$

4:
$\lambda^{t+1}=\lambda^{t}+\rho({\mathbf{A}}{\mathbf{w}}^{t+1}+{\mathbf{B}}{\mathbf{z}}^{t+1}-{\mathbf{c}})$

As we mentioned earlier, simultaneous projection of ${\mathbf{C}}$ to the set
of positive semi-definite matrices and to the box constraint
${\mathbf{C}}_{min}\leq{\mathbf{C}}\leq{\mathbf{C}}_{max}$ is difficult. To
deal with this problem, we eliminate the PSD constraint in Algorithm 3 and
Algorithm 4. Moreover, in those algorithms, we need to carefully tune the
step-size to avoid inconsistency and guarantee convergence.

An alternative approach for solving Problem (9) that avoids these issues is to
solve the dual of the inner maximization problem. Since the maximization
problem is concave with respect to ${\mathbf{C}}$ and ${\mathbf{b}}$, and the
relative interior of the feasible set of constraints is non-empty, the duality
gap is zero. Hence, instead of solving the inner maximization problem, we can
solve its dual which is a minimization problem. Theorem 5 describes the dual
problem of the inner maximization problem in (9). Thus, Problem (9) can be
alternatively formulated as a minimization problem rather than a min-max
problem. We can solve such a constrained minimization problem efficiently via
ADMM algorithm. As we will show, the ADMM algorithm applied to the dual
problem does not need tuning of step-size or doing simultaneous projections to
the box constraints and positive semi-definite (PSD) constraints.

###### Theorem 5

(Dual Problem) The inner maximization problem described in (9) can be
equivalently formulated as:

$\begin{array}[]{lll}\displaystyle{\min_{{\mathbf{A}},{\mathbf{B}},{\mathbf{d}},{\mathbf{e}},{\mathbf{H}}}}&\displaystyle{-\langle{\mathbf{b}}_{min},{\mathbf{d}}\rangle+\langle{\mathbf{b}}_{max},{\mathbf{e}}\rangle-\langle{\mathbf{C}}_{min},{\mathbf{A}}\rangle+\langle{\mathbf{C}}_{max},{\mathbf{B}}\rangle+\lambda\|{\bm{\theta}}\|^{2}}\\\
{\rm
s.t.}&-{\bm{\theta}}{\bm{\theta}}^{T}-{\mathbf{A}}+{\mathbf{B}}-{\mathbf{H}}=0,\\\
&2{\bm{\theta}}-{\mathbf{d}}+{\mathbf{e}}=0,\\\
&{\mathbf{A}},{\mathbf{B}},{\mathbf{d}},{\mathbf{e}}\geq 0,\\\
&{\mathbf{H}}\succeq 0.\end{array}$

Therefore, Problem (9) can be alternatively written as:

$\begin{array}[]{lll}\displaystyle{\min_{{\bm{\theta}},{\mathbf{A}},{\mathbf{B}},{\mathbf{d}},{\mathbf{e}},{\mathbf{H}}}}&\displaystyle{-\langle{\mathbf{b}}_{min},{\mathbf{d}}\rangle+\langle{\mathbf{b}}_{max},{\mathbf{e}}\rangle-\langle{\mathbf{C}}_{min},{\mathbf{A}}\rangle+\langle{\mathbf{C}}_{max},{\mathbf{B}}\rangle+\lambda\|{\bm{\theta}}\|^{2}}\\\
{\rm
s.t.}&-{\bm{\theta}}{\bm{\theta}}^{T}-{\mathbf{A}}+{\mathbf{B}}-{\mathbf{H}}=0,\\\
&2{\bm{\theta}}-{\mathbf{d}}+{\mathbf{e}}=0,\\\
&{\mathbf{A}},{\mathbf{B}},{\mathbf{d}},{\mathbf{e}}\geq 0,\\\
&{\mathbf{H}}\succeq 0.\end{array}$ (13)

Proof  The proof is relegated to Appendix B.

To apply the ADMM method to the dual problem, we require to divide the
optimization variables into two blocks as in (11) such that both sub-problems
in Algorithm 5 can be efficiently solved. To do so, first, we introduce the
auxiliary variables
${\mathbf{d}}^{{}^{\prime}},{\mathbf{e}}^{{}^{\prime}},{\bm{\theta}}^{{}^{\prime}},{\mathbf{A}}^{{}^{\prime}}$
and ${\mathbf{B}}^{{}^{\prime}}$ to the dual problem. Also, let
${\mathbf{G}}={\mathbf{H}}+{\bm{\theta}}^{{}^{\prime}}{\bm{\theta}}^{{}^{\prime}T}$.
Therefore, Problem (40) is equivalent to:

$\begin{array}[]{lll}\displaystyle{\min_{{\bm{\theta}},{\mathbf{A}},{\mathbf{B}},{\mathbf{d}},{\mathbf{e}},{\mathbf{H}}}}&\displaystyle{-\langle{\mathbf{b}}_{min},{\mathbf{d}}\rangle+\langle{\mathbf{b}}_{max},{\mathbf{e}}\rangle-\langle{\mathbf{C}}_{min},{\mathbf{A}}\rangle+\langle{\mathbf{C}}_{max},{\mathbf{B}}\rangle+\lambda\|{\bm{\theta}}\|^{2}}\\\
{\rm s.t.}&{\mathbf{B}}-{\mathbf{A}}={\mathbf{G}},\\\
&2{\bm{\theta}}-{\mathbf{d}}+{\mathbf{e}}=0,\\\
&{\mathbf{A}}={\mathbf{A}}^{{}^{\prime}},{\mathbf{B}}={\mathbf{B}}^{{}^{\prime}},\\\
&{\mathbf{d}}={\mathbf{d}}^{{}^{\prime}},{\mathbf{e}}={\mathbf{e}}^{{}^{\prime}},{\bm{\theta}}={\bm{\theta}}^{{}^{\prime}},\\\
&{\mathbf{A}}^{{}^{\prime}},{\mathbf{B}}^{{}^{\prime}},{\mathbf{d}}^{{}^{\prime}},{\mathbf{e}}^{{}^{\prime}}\geq
0,\\\
&{\mathbf{G}}\succeq{\bm{\theta}}^{{}^{\prime}}{\bm{\theta}}^{{}^{\prime}T}.\end{array}$
(14)

Since handling both constraints on ${\bm{\theta}}$ in Problem (13) is
difficult, we interchange ${\bm{\theta}}$ with ${\bm{\theta}}^{{}^{\prime}}$
in the first constraint. Moreover, the non-negativity constraints on
${\mathbf{A}},{\mathbf{B}},{\mathbf{d}}$ and ${\mathbf{e}}$ are exchanged with
non-negativity constraints on
${\mathbf{A}}^{{}^{\prime}},{\mathbf{B}}^{{}^{\prime}},{\mathbf{d}}^{{}^{\prime}}$
and ${\mathbf{e}}^{{}^{\prime}}$. For the simplicity of presentation, assume
that
${\mathbf{c}}_{1}^{t}={\mathbf{b}}_{min}-{\bm{\mu}}_{d}^{t}+\rho{\mathbf{d}}^{{}^{\prime}t}+{\bm{\eta}}^{t}$,
${\mathbf{c}}_{2}^{t}=-{\mathbf{b}}_{max}-{\bm{\mu}}_{e}^{t}+\rho{\mathbf{e}}^{{}^{\prime}t}-{\bm{\eta}}^{t}$,
${\mathbf{c}}_{3}^{t}=-{\bm{\mu}}_{{\bm{\theta}}}+\rho{\bm{\theta}}^{{}^{\prime}t}-2{\bm{\eta}}^{t}$,
${\mathbf{D}}_{1}^{t}=\rho{\mathbf{A}}^{{}^{\prime}t}-\rho{\mathbf{G}}^{t}+{\bm{\Gamma}}^{t}-{\mathbf{M}}_{A}^{t}+{\mathbf{C}}_{min}$,
and
${\mathbf{D}}_{2}^{t}=\rho{\mathbf{B}}^{{}^{\prime}t}+\rho{\mathbf{G}}^{t}-{\bm{\Gamma}}^{t}-{\mathbf{M}}_{B}^{t}-{\mathbf{C}}_{max}$.
Algorithm 6 describes the ADMM algorithm applied to Problem (14) (see Appendix
C for the derivation).

Algorithm 6 Applying ADMM to the Dual Problem of Robust Linear Regression

1:Given:
${\mathbf{b}}_{min},{\mathbf{b}}_{max},{\mathbf{C}}_{min},{\mathbf{C}}_{max},\lambda,\rho$

2:Initialize:
${\mathbf{C}}_{0}={\mathbf{C}}_{1}=\hat{{\mathbf{C}}},{\mathbf{b}}_{0}={\mathbf{b}}_{1}=\hat{{\mathbf{b}}},\gamma_{0}=0,\gamma_{1}=1$.

3:for $t=0,\ldots,T$ do

4:
${\bm{\theta}}^{t+1}=\frac{1}{6\lambda+7\rho}(2{\mathbf{c}}_{1}^{t}-2{\mathbf{c}}_{2}^{t}-3{\mathbf{c}}_{3}^{t})$

5:
${\mathbf{d}}^{t+1}=\frac{1}{6\lambda+7\rho}(\frac{6\rho+4\lambda}{\rho}{\mathbf{c}}_{1}^{t}+\frac{4\rho+4\lambda}{\rho}{\mathbf{c}}_{2}^{t}+2{\mathbf{c}}_{3}^{t})$

6:
${\mathbf{e}}^{t+1}=\frac{2}{6\lambda+7\rho}(\frac{\rho+2\lambda}{\rho}{\mathbf{c}}_{1}^{t}+\frac{3\rho+2\lambda}{\rho}{\mathbf{c}}_{2}^{t}-{\mathbf{c}}_{3}^{t})$

7:
${\mathbf{A}}^{{}^{\prime}t+1}=\max({\mathbf{A}}^{t}+\frac{{\mathbf{M}}_{A}^{t}}{\rho},0)$

8:
${\mathbf{B}}^{{}^{\prime}t+1}=\max({\mathbf{B}}^{t}+\frac{{\mathbf{M}}_{B}^{t}}{\rho},0)$

9:
${\mathbf{G}}^{t+1}=[{\mathbf{B}}^{t}-{\mathbf{A}}^{t}+\frac{{\bm{\Gamma}}^{t}}{\rho}-{\bm{\theta}}^{{}^{\prime}t}{\bm{\theta}}^{{}^{\prime}tT}]_{+}+{\bm{\theta}}^{{}^{\prime}t}{\bm{\theta}}^{{}^{\prime}tT}$

10:
${\mathbf{d}}^{{}^{\prime}t+1}=\max({\mathbf{d}}^{t}+\frac{{\bm{\mu}}_{d}^{t}}{\rho},0)$

11:
${\mathbf{e}}^{{}^{\prime}t+1}=\max({\mathbf{e}}^{t}+\frac{{\bm{\mu}}_{e}^{t}}{\rho},0)$

12:
${\bm{\theta}}^{{}^{\prime}t+1}=\operatorname*{argmin}_{{\bm{\theta}}^{{}^{\prime}}}\quad\|{\bm{\theta}}^{t+1}-{\bm{\theta}}^{{}^{\prime}}\|^{2}+\langle{\bm{\mu}}_{{\bm{\theta}}}^{t},{\bm{\theta}}^{t+1}-{\bm{\theta}}^{{}^{\prime}}\rangle\quad{\rm
s.t.}\>{\mathbf{G}}^{t+1}\succeq{\bm{\theta}}^{{}^{\prime}T}{\bm{\theta}}^{{}^{\prime}}$

13:
${\mathbf{A}}^{t+1}=\frac{1}{3\rho}(2{\mathbf{D}}_{1}^{t}+{\mathbf{D}}_{2}^{t})$

14:
${\mathbf{B}}^{t+1}=\frac{1}{3\rho}({\mathbf{D}}_{1}^{t}+2{\mathbf{D}}_{2}^{t})$

15:
${\mathbf{M}}_{A}^{t+1}={\mathbf{M}}_{A}^{t}+\rho({\mathbf{A}}^{t+1}-{\mathbf{A}}^{{}^{\prime}t+1})$

16:
${\mathbf{M}}_{B}^{t+1}={\mathbf{M}}_{B}^{t}+\rho({\mathbf{B}}^{t+1}-{\mathbf{B}}^{{}^{\prime}t+1})$

17:
${\bm{\mu}}_{d}^{t+1}={\bm{\mu}}_{d}^{t}+\rho({\mathbf{d}}^{t+1}-{\mathbf{d}}^{{}^{\prime}t+1})$

18:
${\bm{\mu}}_{e}^{t+1}={\bm{\mu}}_{e}^{t}+\rho({\mathbf{e}}^{t+1}-{\mathbf{e}}^{{}^{\prime}t+1})$

19:
${\bm{\mu}}_{{\bm{\theta}}}^{t+1}={\bm{\mu}}_{\theta}^{t}+\rho({{\bm{\theta}}}^{t+1}-{{\bm{\theta}}}^{{}^{\prime}t+1})$

20:
${\bm{\eta}}^{t+1}={\bm{\eta}}^{t}+\rho(2{\bm{\theta}}^{t+1}-{\mathbf{d}}^{t+1}+{\mathbf{e}}^{t+1})$

21:
${\bm{\Gamma}}^{t+1}={\bm{\Gamma}}^{t}+\rho({\mathbf{B}}^{t+1}-{\mathbf{A}}^{t+1}-{\mathbf{G}}^{t+1})$

###### Corollary 6

If the feasible set of Problem (9) has non-empty interior, then Algorithm 6
converges to an $\epsilon$-optimal solution of Problem (14) with a linear rate
$\mathcal{O}(\frac{1}{\epsilon})$.

Proof  Since the inner maximization problem in (9) is convex, and its interior
feasible set is not empty, the duality gap is zero by the Slater’s condition.
Thus, according to Theorem $6.1$ in (He and Yuan, 2015), Algorithm 6 converges
to an optimal solution of the primal-dual problem with a linear rate.
Moreover, the sequence of constraint residuals converges to zero with a linear
rate as well.

###### Remark 7

The optimal solution obtained from the ADMM algorithm can be different from
the one given by Algorithm 3 because we remove the positive semi-definite
constraint on ${\mathbf{C}}$ in the latter. We investigate the difference
between solutions of two algorithms in three cases: First, we generate a small
positive semi-definite matrix ${\mathbf{C}}^{*}$ and the matrix of confidence
intervals ($\bm{\Delta}$) as follows:

${\mathbf{C}}^{*}=\begin{bmatrix}97&40&92\\\ 40&17&38\\\
92&38&88\end{bmatrix},\quad\bm{\Delta}=\begin{bmatrix}0.2&0.3&0.2\\\
0.3&0.1&0.2\\\ 0.1&0.3&0.1\end{bmatrix}.$

Moreover, let ${\mathbf{b}}^{*}$ and $\bm{\delta}$ are generated as follows:

${\mathbf{b}}^{*}=\begin{bmatrix}6.65\\\ 8.97\\\
5.40\end{bmatrix},\quad\bm{\delta}=\begin{bmatrix}0.1\\\ 0.2\\\
0.2\end{bmatrix}.$

Initializing both algorithms with a random matrix within
${\mathbf{C}}_{min}={\mathbf{C}}^{*}-\bm{\Delta}$ and
${\mathbf{C}}_{max}={\mathbf{C}}^{*}+\bm{\Delta}$, and a random vector within
${\mathbf{b}}_{min}={\mathbf{b}}^{*}-\bm{\delta}$ and
${\mathbf{b}}_{max}={\mathbf{b}}^{*}+\bm{\delta}$, ADMM algorithm returns a
different solution from Algorithm 3.

Besides, the difference in the performance of algorithms during the test phase
can be observed in the experiments on synthetic datasets depicted in Figure 6
as well, especially when the number of samples is smaller.

###  3.5 Performance Guarantees for RIFLE

In the previous section, we discussed how to efficiently solve the robust
linear regression problem in the presence of missing values up to optimality.
A natural question in this context is the performance of the proposed
algorithms in the previous section on the unseen test data points. Theorem 8
investigates such question from two perspectives: Under the assumption that
the missing values are distributed completely at random (MCAR), the optimal
solution of Problem (3) is consistent, meaning as the number of samples goes
to infinity, the optimal solution converges to the ground truth parameter
${\bm{\theta}}^{*}$. Moreover, for the finite case, Theorem 8 part (b) states
that with the proper choice of confidence intervals, with high probability,
the population (test) loss is bounded by the training loss.

###### Theorem 8

(a) Consistency of the Covariance Estimator: If the missing pattern of the
data follows MCAR, the estimator proposed in Algorithm 2 converges to its
actual value. Thus, when $n$ goes to infinity, the estimated confidence
intervals for each entry of the covariance matrix converges to $0$.

$\lim_{n\rightarrow\infty}|\hat{{\mathbf{C}}}_{ij}-{\mathbf{C}}^{*}_{ij}|=0$
(15)

(b) Let ${\mathbf{X}}^{n\times d}$, ${\mathbf{y}}$ be the matrix of data, and
the corresponding assigned labels to each data point. Assume that each entry
of ${\mathbf{X}}$ and ${\mathbf{y}}$ is missing with probability of $p$.
Moreover, Let $({\mathbf{C}}^{*},{\mathbf{b}}^{*})$ be the ground-truth
parameters, and
$(\tilde{{\bm{\theta}}},\tilde{{\mathbf{C}}},\tilde{{\mathbf{b}}})$ be a
saddle point of Problem (8). Then with probability of at least
$1-\frac{d^{2}}{2nc}$, we have:

$L_{test}(\tilde{{\bm{\theta}}})\leq L_{train}(\tilde{{\bm{\theta}}})$ (16)

Proof  The proof is relegated to Appendix B.

###  3.6 How to Handle Missing Values in the Test Phase

Let $\hat{{\bm{\theta}}}$ be the optimal parameter of the robust linear
regression problem (3) obtained by running one of the algorithms described in
the previous section on the training dataset. If, a test data point
${\mathbf{x}}$ does not contain missing values, then $\hat{y}$ can be easily
computed by $\hat{y}=\hat{{\bm{\theta}}}^{T}{\mathbf{x}}$. There are several
options for handling data points with several missing entries. A naïve idea is
to make all of the missing entries zero. However, this approach can affect the
performance of the model drastically. On the other hand, it is possible to
restrict the training dataset to the features that are available in the given
test data point. The main problem of this approach is that for any test data
point we must solve the entire min-max problem that can be extremely
inefficient unless the set of test data points is small. In practice, we keep
the optimal ${\mathbf{C}}$ and ${\mathbf{b}}$ obtained from optimizing problem
(9), and find a ${\bm{\theta}}$ for each one of the testing points, by
restricting the optimal ${\mathbf{C}}$ and ${\mathbf{b}}$ to the features that
are available for the testing point. Thus, it suffices to compute
${\bm{\theta}}^{{}^{\prime}}=({\mathbf{C}}^{{}^{\prime}}+\lambda\mathbf{I})^{-1}{\mathbf{b}}^{{}^{\prime}}$,
where ${\mathbf{C}}^{{}^{\prime}}$ and ${\mathbf{b}}^{{}^{\prime}}$ are
obtained by eliminating the features that are not available for the current
test point from the optimal ${\mathbf{C}}$ and ${\mathbf{b}}$.

###  3.7 Imputation of Missing Values via Robust Linear Regression

We impute a given dataset containing missing values feature by feature. In
other words, at each iteration, we consider one feature as the target variable
$y$, and the rest as input data features (${\mathbf{X}}$). Then, we learn a
robust linear regression model by solving the min-max problem (8). Let the
obtained optimal solutions be ${\mathbf{C}}^{*},{\mathbf{b}}^{*},$ and
${\bm{\theta}}^{*}$. For a given row, we restrict ${\mathbf{C}}^{*}$ and
${\mathbf{b}}^{*}$ to available features on that row. Then, we find the
corresponding optimal ${\bm{\theta}}$ in closed-form for the restricted
${\mathbf{C}}^{*}$ and ${\mathbf{b}}^{*}$ (similar to what we explained for
handling test data points containing missing values). Thus, to impute each
feature of the dataset, we only solve the min-max problem (8) once, and for
each row we find its optimal ${\bm{\theta}}$ by restricting ${\mathbf{C}}^{*}$
and ${\mathbf{b}}^{*}$ to available entries of that row. Note that, if the
dataset only contains a few missing patterns for different rows (for instance
we have $3$ different missing patterns in Figure 1) we can find the optimal
${\bm{\theta}}$ only for the distinct patterns instead of solving the problem
with respect to ${\bm{\theta}}$ for each row separately. Since imputation of
each feature is completely independent of the others, features can be
distributed to multiple cores (or computers) without losing performance.

##  4 Robust Classification Framework

In this section we specialize the proposed framework in Section 2 to the
classification task in the presence of missing values. For simplicity of
presentation, assume that $y$, the target feature, is a binary variable. A
basic idea for learning a robust classifier in this context is to apply the
proposed robust linear regression algorithm presented in (3) to the given
dataset assuming that $y$ is continuous. Then, the assigned labels can be
obtained by thresholding on the output of the robust regression model. As it
is discussed earlier, since the linear regression loss depends on data through
${\mathbf{X}}^{T}{\mathbf{X}}$ and ${\mathbf{X}}^{T}{\mathbf{y}}$, it suffices
to estimate first and second order moments of features (means and covariances)
using the available data.

For the classification task, the most commonplace measure of accuracy is 0-1
loss. Since 0-1 loss is highly non-convex, and intractable in general,
different classification algorithms minimize an approximation of 0-1 loss. A
relevant question in this context is whether we can guarantee a bound on 0-1
loss by minimizing the least-square loss proposed in (8). The following
theorem addresses such a question.

Beside the idea of thresholding on the robust linear regression output, one
can adapt the general framework described in Equation (4) to a classifier
model such as logistic regression or support vector machines. Let the data
points ${\mathbf{z}}_{1},\dots,{\mathbf{z}}_{n}$ come from a parametric
distribution $P_{{\bm{\theta}}^{*}}({\mathbf{z}})$. Assume that
$\hat{{\bm{\theta}}}$ is the estimation of ${\bm{\theta}}^{*}$ based on the
given $n$ incomplete data points. The general robust classification
formulation problem in the presence of missing values can be described as
follows:

$\min_{\mathbf{w}}\;\max_{\bm{\delta}_{\min}\leq{\bm{\theta}}-\hat{{\bm{\theta}}}\leq\bm{\delta}_{\max}}\;\;\mathbb{E}_{P_{\theta}}[\ell({\mathbf{z}},\mathbf{w})],$
(17)

where $\bm{\delta}_{\min}$ and $\bm{\delta}_{\max}$ are the estimated
confidence intervals for the difference of ${\bm{\theta}}$ and
$\hat{{\bm{\theta}}}$. In the following section, we present the robust
logistic regression in the presence of missing values as a case study of the
above framework.

###  4.1 Robust Linear Discriminant Analysis

The problem of training a logistic regression model on datasets containing
missing values has been studied extensively in the literature. Beside the
deletion of missing values, and imputation based approaches, (Fung and Wrobel,
1989) formulates the logistic regression as an equivalent linear discriminant
analysis problem. Mathematically speaking, they assume that the data points
assigned to a specific label follow a Gaussian distribution
$({\mathbf{x}}|y=i)\sim N({\bm{\mu}}_{i},\bm{\Sigma}))$. They use the
available data to estimate the parameters of each Gaussian distribution.
Therefore, the parameters of the logistic regression model can be assigned
based on the estimated parameters of the Gaussian distributions for different
classes. Similar to the linear regression case, the estimations of means and
covariances are unbiased only when the data satisfies MCAR condition.
Moreover, when the number of data points in the dataset is very small
$\Big{(}\frac{n(1-p)}{d}=\mathcal{O}(1)\Big{)}$ the variance of the
estimations can be very high. Thus, to train a logistic regression model which
is robust to the percentage and different types of missing values, we
specialize the general robust classification framework formulated in (17) to
the logistic regression model. Assume that $\mathbf{z}=(\mathbf{x},y)\sim
N({\bm{\mu}},\bm{\Sigma})$. We aim to find the optimal solution of the
following problem:

$\begin{array}[]{lll}\displaystyle{\min_{\mathbf{w}}\max_{\bm{\mu}_{0},\bm{\mu}_{1},\bm{\Sigma}_{0},\bm{\Sigma}_{1}}}&\displaystyle{\mathbb{E}_{\mathbf{x}|y=1\sim
N(\bm{\mu}_{1},\bm{\Sigma}_{1})}\Big{[}-\log\Big{(}\sigma(\mathbf{w}^{T}\mathbf{x})\Big{)}\Big{]}\mathbb{P}(y=1)}\>+\\\
&\displaystyle{\mathbb{E}_{\mathbf{x}|y=0\sim
N(\bm{\mu}_{0},\bm{\Sigma}_{0})}\Big{[}-\log\Big{(}1-\sigma(\mathbf{w}^{T}\mathbf{x})\Big{)}\Big{]}}\mathbb{P}(y=0)\\\
{\rm s.t.}&\bm{\mu}_{\min_{0}}\leq\bm{\mu}_{0}\leq\bm{\mu}_{\max_{0}}\\\
&\bm{\mu}_{\min_{1}}\leq\bm{\mu}_{1}\leq\bm{\mu}_{\max_{1}}\\\
&\bm{\Sigma}_{\min_{0}}\leq\bm{\Sigma}_{0}\leq\bm{\Sigma}_{\max_{0}}\\\
&\bm{\Sigma}_{\min_{1}}\leq\bm{\Sigma}_{1}\leq\bm{\Sigma}_{\max_{1}}\\\
\end{array}$ (18)

###  4.2 No Missing Values in Target Variable

To solve Problem (18), first we focus on the case when the target variable has
no missing values. In this case, each data point contributes to the estimation
of either $({\bm{\mu}}_{1},\bm{\Sigma}_{1})$ or
$({\bm{\mu}}_{0},\bm{\Sigma}_{0})$. Similar to the robust linear regression
case, we can apply Algorithm 2 to estimate the confidence intervals for
$\bm{\mu}_{i},\bm{\Sigma}_{i}$ using data points whose target variable equals
$i$ ($y=i$).

Obviously, the objective function is convex in $\mathbf{w}$, since the
logistic regression loss is convex, and the expectation of loss can be seen as
a weighted summation, which is convex. Thus, fixing ${\bm{\mu}},\bm{\Sigma}$
the outer minimization problem can be solved with respect to ${\mathbf{w}}$
using standard first order methods such as the gradient descent.

Although, the robust reformulation of logistic regression stated in (18) is
convex in ${\mathbf{w}}$ and concave in ${\bm{\mu}}_{0}$ and ${\bm{\mu}}_{1}$,
the inner maximization problem is intractable with respect to $\bm{\Sigma}$.
We approximate Problem (18) in the following manner:

$\begin{array}[]{lll}\displaystyle{\min_{\mathbf{w}}\max_{\bm{\mu}_{0},\bm{\Sigma}_{0},\bm{\mu}_{1},\bm{\Sigma}_{1}}}&\displaystyle{\pi_{1}\mathbb{E}_{\mathbf{x}|y=1\sim
N({\bm{\mu}}_{1},\bm{\Sigma}_{1})}\Big{[}-\log\Big{(}\sigma(\mathbf{w}^{T}\mathbf{x})\Big{)}\Big{]}+\pi_{0}\mathbb{E}_{\mathbf{x}|y=0\sim
N({\bm{\mu}}_{0},\bm{\Sigma}_{0})}\Big{[}-\log\Big{(}1-\sigma(\mathbf{w}^{T}\mathbf{x})\Big{)}\Big{]}},\\\
{\rm s.t.}&\bm{\mu}_{\min_{0}}\leq\bm{\mu}_{0}\leq\bm{\mu}_{\max_{0}}\\\
&\bm{\mu}_{\min_{1}}\leq\bm{\mu}_{1}\leq\bm{\mu}_{\max_{1}}\\\
&\bm{\Sigma}_{0}\in\\{\bm{\Sigma}_{01},\bm{\Sigma}_{02},\dots,\bm{\Sigma}_{0k}\\}\\\
&\bm{\Sigma}_{1}\in\\{\bm{\Sigma}_{11},\bm{\Sigma}_{12},\dots,\bm{\Sigma}_{1k}\\},\end{array}$
(19)

where $\pi_{1}=\mathbb{P}(y=1)$ and $\pi_{0}=\mathbb{P}(y=0)$. To compute
optimal $\bm{\mu}_{0}$ and $\bm{\mu}_{1}$, we have:

$\max_{\mu_{1}}\quad\mathbb{E}_{\mathbf{x}\sim
N(\mu_{1},\Sigma_{1})}\Big{[}-\log\Big{(}\sigma(\mathbf{w}^{T}\mathbf{x})\Big{)}\Big{]}\quad{\rm
s.t.}\quad\bm{\mu}_{min}\leq\bm{\mu}_{1}\leq\bm{\mu}_{max}$ (20)

###### Theorem 9

Let ${\mathbf{a}}[i]$ be the $i$-th element of vector ${\mathbf{a}}$. The
optimal solution of Problem (20) has the following form:

${\bm{\mu}}^{*}_{1}[i]=\begin{cases}{\bm{\mu}}_{\max}[i],\quad\quad\mbox{if
}\>{\mathbf{w}}[i]\leq 0\\\ {\bm{\mu}}_{\min}[i],\quad\quad\mbox{if
}\>{\mathbf{w}}[i]>0.\end{cases}$ (21)

Note that we relaxed (18) by taking the maximization problem over a finite set
of $\bm{\Sigma}$ estimations. We estimate each $\bm{\Sigma}$ by bootstraping
on the available data using Algorithm 2. Define $f_{i}({\mathbf{w}})$ as:

$\begin{array}[]{lll}f_{i}({\mathbf{w}})=\displaystyle{\pi_{1}\mathbb{E}_{\mathbf{x}\sim
N({\bm{\mu}}^{*}_{1},\bm{\Sigma}_{i1})}\Big{[}-\log\Big{(}\sigma(\mathbf{w}^{T}\mathbf{x})\Big{)}\Big{]}}\end{array}$
(22)

Similarly we can define:

$\begin{array}[]{lll}g_{i}({\mathbf{w}})=\displaystyle{\pi_{0}\mathbb{E}_{\mathbf{x}\sim
N({\bm{\mu}}^{*}_{0},\bm{\Sigma}_{i0})}\Big{[}-\log\Big{(}1-\sigma(\mathbf{w}^{T}\mathbf{x})\Big{)}\Big{]}}\end{array}$
(23)

Since, the maximization problem is over a finite set, we can rewrite Problem
(19) as:

$\begin{array}[]{lll}\displaystyle{\min_{\mathbf{w}}\max_{i,j\in\\{1,\dots,k\\}}}&\displaystyle{f_{i}({\mathbf{w}})+g_{j}({\mathbf{w}})}=\displaystyle{\min_{\mathbf{w}}\max_{p_{1},\dots,p_{k},q_{1},\dots,q_{k}}}&\displaystyle{\sum_{i=1}^{k}p_{i}f_{i}({\mathbf{w}})+\sum_{j=1}^{k}p_{i}g_{j}({\mathbf{w}})}\\\
{\rm s.t.}&\sum_{i=1}^{k}p_{i}=1,\quad p_{i}\geq 0\\\
&\sum_{j=1}^{k}q_{j}=1,\quad q_{j}\geq 0\end{array}$ (24)

To solve this problem faster, we add a regularization term to the maximization
problem as follows:

$\begin{array}[]{lll}\displaystyle{\min_{\mathbf{w}}\max_{p_{1},\dots,p_{k},q_{1},\dots,q_{k}}}&\displaystyle{\sum_{i=1}^{k}p_{i}f_{i}({\mathbf{w}})}-\delta\sum_{i=1}^{k}p_{i}^{2}+\sum_{j=1}^{k}q_{j}g_{j}({\mathbf{w}})-\delta\sum_{j=1}^{k}q_{j}^{2}\\\
{\rm s.t.}&\sum_{i=1}^{k}p_{i}=1,\quad p_{i}\geq 0\\\
&\sum_{j=1}^{k}q_{j}=1,\quad q_{j}\geq 0\end{array}$ (25)

First, we show how to solve the inner maximization problem. Note that the
$p_{i}$’s and $q_{i}$’s are independent of each other. We show how to find
optimal $p_{i}$’s. Optimizing with respect to $q_{i}$’s is similar. Since, the
maximization problem is a constrained quadratic program, we can write the
Lagrangian function as follows:

$\begin{array}[]{lll}\displaystyle{\max_{p_{1},\dots,p_{k}}}&\displaystyle{\sum_{i=1}^{k}p_{i}f_{i}({\mathbf{w}})}-\delta\sum_{i=1}^{k}p_{i}^{2}-\lambda(\sum_{i=1}^{k}p_{i}-1)\\\
{\rm s.t.}\quad p_{i}\geq 0\end{array}$ (26)

Having the optimal $\lambda$, the above problem has the closed form solution
with respect to each $p_{i}$, where can be written as:

$p^{*}_{i}=\left[\frac{-\lambda+f_{i}}{2\delta}\right]_{+}$

Thus, it suffices to find the optimal $\lambda$ in (26). Since $p^{*}_{i}$ is
a non-increasing function with respect to $\lambda$, we can find the optimal
value of $\lambda$ using the following bisection algorithm. Algorithm 7
demonstrates how to find an $\epsilon$-optimal $\lambda$ and $p_{i}^{*}$’s
efficiently using the bisection idea.

Algorithm 7 Finding the optimal $\lambda$ and $p_{i}$’s using the bisection
idea

1:Initialize:
$\lambda_{\textrm{low}}=0,\lambda_{\textrm{high}}=\max_{i}f_{i},p_{i}=0\quad\forall
i\in\\{1,2,\dots,k\\}$.

2:while $|\sum_{i=1}^{n}p_{k}-1|>\epsilon$ do

3: $\lambda=\frac{\lambda_{\textrm{low}}+\lambda_{\textrm{high}}}{2}$

4: Set $p_{i}=[\frac{-\lambda+f_{i}}{2\delta}]_{+}\quad\forall
i\in\\{1,2,\dots,k\\}$

5: if $\sum_{i=1}^{k}p_{i}<1$ then

6: $\lambda_{\textrm{high}}=\lambda$

7: else

8: $\lambda_{\textrm{low}}=\lambda$

9:return $\lambda,p_{1},p_{2},\dots,p_{k}$.

###### Remark 10

An alternative method for finding optimal $\lambda$, and $p_{i}$’s is to sort
$f_{i}$ values in $\mathcal{O}(k\log k)$ first, and then finding the smallest
$f_{i}$ such that if we set $\lambda=f_{i}$, the sum of $p_{i}$’s is bigger
than $1$ (let $j$ be the index of that value). Without loss of generality
assume that $f_{1}\leq\dots\leq f_{k}$. Then,
$\sum_{i=j}^{k}\frac{-\lambda+f_{i}}{2\delta}=1$, which has a closed-form
solution with respect to $\lambda$.

To update ${\mathbf{w}}$, we need to solve the following optimization problem:

$\begin{array}[]{lll}\displaystyle{\min_{\mathbf{w}}}&\displaystyle{\sum_{i=1}^{k}p^{*}_{i}f_{i}({\mathbf{w}})+\sum_{j=1}^{k}q^{*}_{j}g_{i}({\mathbf{w}})},\\\
\end{array}$ (27)

Similar to the standard statistical learning framework, we solve the following
empirical risk minimization problem by applying the gradient descent to
${\mathbf{w}}$ on a finite data sample. Define $\hat{f}_{i}$ as follows:

$\hat{f}_{i}({\mathbf{w}})=\pi_{1}\sum_{t=1}^{n}\Big{[}-\log\Big{(}\sigma(\mathbf{w}^{T}\mathbf{x}_{t})\Big{)}\Big{]},$
(28)

where ${\mathbf{x}}_{1},\dots,{\mathbf{x}}_{n}$ are generated from the
distribution $\mathcal{N}(\mathbf{\mu}_{1}^{*},\mathbf{\Sigma_{1i}})$. the
empirical risk minimization problem can be written as:

$\begin{array}[]{lll}\displaystyle{\min_{\mathbf{w}}}&\displaystyle{\sum_{i=1}^{k}p^{*}_{i}\hat{f}_{i}({\mathbf{w}})+\sum_{j=1}^{k}q^{*}_{j}\hat{g}_{i}({\mathbf{w}})},\end{array}$
(29)

Algorithm 8 summarizes the robust linear discriminant analysis method for the
case where the label of all data points is available.

Algorithm 8 Robust Normal Discriminant Analysis

1:Input: ${\mathbf{X}}_{0},{\mathbf{X}}_{1}$: matrix of data points with
labels $0$ and $1$ respectively, $T:$ Number of iterations, $\alpha:$ Step-
size.

2:Estimate $\bm{\mu}_{\textrm{min}_{0}}$ and $\bm{\mu}_{\textrm{max}_{0}}$
using the available entries of ${\mathbf{X}}_{0}$.

3:Estimate $\bm{\mu}_{\textrm{min}_{1}}$ and $\bm{\mu}_{\textrm{max}_{1}}$
using the available entries of ${\mathbf{X}}_{1}$.

4:Estimate $\mathbf{\Sigma}_{01},\dots,\mathbf{\Sigma}_{0k}$ using bootstrap
estimator on the available data of ${\mathbf{X}}_{0}$.

5:Estimate $\bm{\Sigma}_{11},\dots,\bm{\Sigma}_{1k}$ using bootstrap estimator
on the available data of ${\mathbf{X}}_{1}$.

6:for $i=1,\ldots,T$ do

7: Compute ${\bm{\mu}}^{*}_{1}$ and ${\bm{\mu}}^{*}_{0}$ by Equation (21).

8: Find optimal $p_{1},\dots,p_{k}$, and $q_{1},\dots,q_{k}$ using Algorithm
7.

9:
${\mathbf{w}}={\mathbf{w}}-\alpha\Bigg{(}\sum_{i=1}^{k}p^{*}_{i}\nabla\hat{f}_{i}({\mathbf{w}})+\sum_{j=1}^{k}q^{*}_{j}\nabla\hat{g}_{i}({\mathbf{w}})\Bigg{)}$

Theorem 11 demonstrates the convergence of gradient descent algorithm applied
to (29) in
$\mathcal{O}\Big{(}\frac{k}{\epsilon}\log(\frac{M}{\epsilon})\Big{)}$
iterations to an $\epsilon$-optimal solution.

###### Theorem 11

Assume that $M=\max_{i}f_{i}$. Gradient descent algorithm requires
$\mathcal{O}\Big{(}\frac{k}{\epsilon}\log(\frac{M}{\epsilon})\Big{)}$ gradient
evaluations for converging to an $\epsilon$-optimal saddle point of the
optimization problem (29).

###  4.3 Presence of Missing Values in the Target Feature

If the target feature contains missing values, the proposed algorithm for
solving the optimization problem (19) does not exploit the data points whose
target feature is not available. However, such points can contain valuable
statistical information about the underlying distribution of data. Thus, we
apply an Expectation Maximization (EM) procedure on the dataset as follows:

Assume that a dataset consisting of $n+m$ samples. Let
$({\mathbf{X}}_{1},y_{1}),\dots,({\mathbf{X}}_{n},y_{n})$ be $n$ samples whose
target variable is available and
$({\mathbf{X}}_{n+1},z_{1}),\dots,({\mathbf{X}}_{n+m},z_{m})$ are $m$ samples
where their corresponding labels are missing. Similar to the previous case, we
assume:

${\mathbf{X}}_{i}|z_{i}=j\sim\mathcal{N}({\bm{\mu}}_{j},\bm{\Sigma}_{j}),\quad
j=0,1.$

Thus, the probability of observing a data point $x_{i}$ can be written as:

$\displaystyle
P({\mathbf{X}}_{i}={\mathbf{x}}_{i})=\pi_{0}P({\mathbf{X}}_{i}={\mathbf{x}}_{i}|z_{i}=0)+\pi_{1}P({\mathbf{X}}_{i}={\mathbf{x}}_{i}|{\mathbf{z}}_{i}=1)$
$\displaystyle=\pi_{0}\mathcal{N}({\mathbf{x}}_{i};\mu_{0},\bm{\Sigma}_{0})+\pi_{1}\mathcal{N}({\mathbf{x}}_{i};{\bm{\mu}}_{1},\bm{\Sigma}_{1})$

The log of likelihood function can be formulated as follows:

$\ell({\bm{\mu}}_{0},\bm{\Sigma}_{0},{\bm{\mu}}_{1},\bm{\Sigma}_{1})=\sum_{i=1}^{n+m}\log\Big{(}\pi_{0}\mathcal{N}({\mathbf{x}}_{i};{\bm{\mu}}_{0},\bm{\Sigma}_{0})+\pi_{1}\mathcal{N}(x_{i};{\bm{\mu}}_{1},\bm{\Sigma}_{1})\Big{)}$

We apply Expectation Maximization procedure to jointly update
$\bm{\Sigma}_{0},\bm{\Sigma}_{1},{\bm{\mu}}_{0},{\bm{\mu}}_{1}$ and $z_{i}$’s.
Note that, the posterior distribution of $z_{i}$ can be written as:

$P(Z_{i}=t|{\mathbf{X}}_{i}={\mathbf{x}}_{i})=\frac{P({\mathbf{X}}_{i}={\mathbf{x}}_{i}|Z_{i}=t)P(Z_{i}=t)}{P({\mathbf{X}}_{i}={\mathbf{x}}_{i})}=\frac{\pi_{t}\mathcal{N}({\mathbf{x}}_{i};{\bm{\mu}}_{t},\bm{\Sigma}_{t})}{P({\mathbf{X}}_{i}={\mathbf{x}}_{i})}$

In E-step, we update $z_{i}$ values by comparing the posterior probabilities
for two possible labels. Precisely, we assign label $1$ to $Z_{i}$ if and only
if:

$\pi_{1}\mathcal{N}({\mathbf{x}}_{i};{\bm{\mu}}_{1},\bm{\Sigma}_{1})>\pi_{0}\mathcal{N}({\mathbf{x}}_{i};{\bm{\mu}}_{0},\bm{\Sigma}_{0})$

In M-step, we estimate
$\bm{\Sigma}_{0},\bm{\Sigma}_{1},{\bm{\mu}}_{0},{\bm{\mu}}_{1},\pi_{0}$ and
$pi_{0}$ by fixing the $z_{i}$ values. Since in M-step all labels (both
already available $y_{i}$’s and estimated $z_{i}$’s in E-step) are assigned,
updating the aforementioned parameters can be done as follows:

$\displaystyle{\bm{\mu}}_{1}[j]=\frac{1}{|\mathcal{S}_{1}\cap\mathcal{T}_{j}|}\sum_{i\in\mathcal{S}_{1}\cap\mathcal{T}_{j}}x_{i}[j]$
(30)
$\displaystyle{\bm{\mu}}_{0}[j]=\frac{1}{|\mathcal{S}_{0}\cap\mathcal{T}_{j}|}\sum_{i\in\mathcal{S}_{0}\cap\mathcal{T}_{j}}x_{i}[j]$
(31)
$\displaystyle\bm{\Sigma}_{1}[i][j]=\frac{1}{|\mathcal{S}_{1}\cap\mathcal{T}_{i}\cap\mathcal{T}_{j}|}\sum_{t\in\mathcal{S}_{1}\cap\mathcal{T}_{i}\cap\mathcal{T}_{j}}x_{t}[i]x_{t}[j]$
(32)
$\displaystyle\bm{\Sigma}_{0}[i][j]=\frac{1}{|\mathcal{S}_{0}\cap\mathcal{T}_{i}\cap\mathcal{T}_{j}|}\sum_{t\in\mathcal{S}_{0}\cap\mathcal{T}_{i}\cap\mathcal{T}_{j}}x_{t}[i]x_{t}[j]$
(33) $\displaystyle\pi_{1}=\frac{|S_{1}|}{|S_{0}\cup S_{1}|}$ (34)
$\displaystyle\pi_{0}=\frac{|S_{0}|}{|S_{0}\cup S_{1}|}$ (35)

We apply the M-step and E-step algorithm iteratively to obtain $\Sigma_{1}$
and $\Sigma_{0}$. Based on the random initialization of $z_{i}$’s we can
obtain different values for $\mu_{0},\mu_{1},\Sigma_{0}$ and $\Sigma_{1}$.
Having these estimations, we apply Algorithm 8 to solve the robust normal
discriminant analysis formulated in (19).

Algorithm 9 Expectation Maximization Procedure for Learning a Robust Normal
Discriminant Analysis

1:Input: $T$: Number of EM iterations, $k$: Number of covariance estimations
at each iteration.

2:Initialize: Set each missing labels randomly to 0 or 1.

3:for $i=1,\ldots,T$ do

4: Estimate $k$ covariance matrices by sampling with replacement from the
available entries

5: Find an optimal ${\mathbf{w}}$ for Problem (24)

6: Update the missing labels using the new ${\mathbf{w}}$ obtained above.

##  5 Numerical Results

To evaluate the proposed algorithms for learning regression and classification
models in the presence of missing values, we compare their performance with
several state-of-the-art approaches for data imputation on synthetic and real-
world datasets. The experiments are designed in a manner that the sensitivity
of the model to factors such as the number of samples, dimension of data,
types and proportion of missing values can be evaluated. In Appendix A, we
introduce the datasets used in the experiments.

###  5.1 Generating MCAR and MNAR Missing Values

To evaluate RIFLE and other state-of-the-art imputation approaches, we require
to have access to the ground-truth values of the missing entries. Hence, we
artificially mask a proportion of available data entries and predict them with
different imputation methods. The closer the predicted values to the actual
values, the better the performance of the method. In all experiments,
generated missing entries follow either a missing completely at random (MCAR)
pattern of missing not at random (MNAR) pattern.

Formally, the distribution of missing values in a dataset follows a missing
completely at random (MCAR) pattern if the probability of having a missing
value for a given entry is constant, independent of other available and
missing entries. On the other hand, a dataset follows a Missing At Random
(MAR) pattern if the missingness of each entry only depends on the available
data of other features. Finally, if the distribution of missing values does
not follow a MCAR or MAR pattern, we call it missing not at random (MNAR).

To generate the MCAR pattern on a given dataset, we fix a constant probability
$0<p<1$ and make each data entry unavailable with probability of $p$. We also
generate MNAR patterns for testing. The generation of the MNAR pattern is
based on the idea that if the value of an entry is farther from the mean of
its corresponding feature, then the probability of missingness for that entry
is larger. Algorithm 10 describes the procedure of generating MNAR missing
values for a given column of a dataset:

Algorithm 10 Generating MNAR Pattern for a Given Column of a Dataset

1:Input: $x_{1},x_{2},\dots,x_{n}$: The entries of the current column in the
dataset, $a,b$: Hyper-parameters controlling the percentage of missing values

2:Initialize: Set $\mu=\frac{1}{n}\sum_{i=1}^{n}x_{i}$ and
$\sigma^{2}=\frac{1}{n}\sum_{i=1}^{n}x_{i}^{2}-\mu^{2}$.

3:for $i=1,\ldots,n$ do

4: $x^{{}^{\prime}}_{i}=\frac{x_{i}-\mu}{\sigma}$

5: $p_{i}=F(a|x^{{}^{\prime}}_{i}|+b)$

6: Set $x_{i}=*$ with probability of $p_{i}$

Note that $F$ in the above algorithm is the cumulative distribution function
of a standard Gaussian random variable. $a$ and $b$ control the percentage of
missing values in the given column. As $a$ and $b$ increase, the probability
of having more missing values is higher. Since the availability of each data
entry clearly depends on its value, the generated missing pattern is missing
not at random (MNAR).

###  5.2 Performance Evaluation

To measure the performance of RIFLE and the existing approaches on a
regression task for a given test dataset consisting of $N$ data points, we use
normalized root mean squared error (NRMSE) defined as:

$\textrm{NRMSE}=\frac{\sqrt{\frac{1}{N}\sum_{i=1}^{N}(y_{i}-\hat{y}_{i})^{2}}}{\sqrt{\frac{1}{N}\sum_{i=1}^{N}(y_{i}-\bar{y})^{2}}}$

where $y_{i}$, $\hat{y}_{i}$, and $\bar{y}$ represent the true value of the
$i$-th data point, the predicted value of the $i$-th data point, and the
average of true values of data points respectively.

###  5.3 Convergence of ADMM, Nesterov and Projected Gradient Ascent
Algorithms

We presented three algorithms for solving the robust linear regression problem
formulated in (9): Projected gradient ascent, Nesterov acceleration method,
and Alternating Direction Method of Multipliers (ADMM) applied on the dual
problem. We established the convergence rate of the gradient ascent and
Nesterov acceleration methods in Theorem 3 and Theorem 4. To investigate the
convergence of ADMM algorithm and its dependence to $\rho$, we perform
Algorithm 6 on Super Conductivity dataset (Description in Appendix A) with
$30\%$ MCAR missing values. Figure 3 demonstrates the convergence of ADMM
algorithm for multiple values of $\rho$ applied to the Super Conductivity
dataset as described above. Describe number of train and test samples, number
of features. As it can be observed, decreasing the value of $\rho$ accelerates
the ADMM convergence to the optimal value.

Figure 3: Convergence of ADMM algorithm to the optimal solution of Problem
(13) for different values of $\rho$. The left plot measures the objective
function of Problem (13) per iteration (without considering the constraints),
while the right plot demonstrates the constraint violation of the algorithm
per iteration. The constraint violation can be measured by adding all
regularization terms in the augmented Lagrangian function formulated in
Problem (42).

In the next experiment, we compare the three proposed algorithms in terms of
number of iterations required to reach a certain level of test accuracy on
Super Conductivity dataset. The number of training samples is $1000$, and it
contains $40\%$ of MCAR missing values on both input features and the target
variable. The test dataset contain $2000$ samples. As depicted in Figure 4,
ADMM and Nesterov algorithms require less number of iterations to reach to the
optimal level of accuracy compared to Algorithm 3. However, the cost per
iteration of ADMM algorithm (Algorithm 6) is higher than the Nesterov
acceleration and Algorithm 3.

Figure 4: Performance of Nesterov acceleration method, projected gradient
ascent, and ADMM on Super Conductivity dataset vs the number of iterations.

###  5.4 Consistency of RIFLE

An algorithm learning a parameterized model with the ground truth parameter
$\theta^{*}$ is consistent if it can recover $\theta^{*}$, as the number of
training samples goes to infinity. Besides, among different learning
algorithms, the one that can recover the true parameter faster and with fewer
number of samples is preferable.

Figure 5: Consistency of ADMM (Algorithm 6) and Projected Gradient Ascent on
function $g$ (Algorithm 3) on the synthetic datasets with $40\%$, $60\%$ and
$80\%$ missing values.

In the first experiment, we investigate the consistency of RIFLE on synthetic
datasets with different proportions of missing values. The synthetic datasets
have $50$ input features following a jointly normal distribution. The target
variable is a linear function of input features added to a mean zero normal
noise. As observed in Figure 5, by increasing the number of samples, the NRMSE
of our framework converges to $0.01$ which is the standard deviation of the
zero-mean Gaussian noise added to each target value (the dashed line).

Moreover, as depicted in Figure 6, RIFLE requires fewer samples to recover the
ground truth parameters of the model compared to MissForest and KNN Imputer.
We compare RIFLE with these two methods since other existing packages do not
support the imputation of the test data based on the learned parameters during
the training phase.

Figure 6: Comparing the consistency of RIFLE, MissForest, and KNN Imputer on a
syntetic dataset containing $40\%$ of missing values.

###  5.5 Data Imputation via RIFLE

As explained in Section 3, while the primary goal of RIFLE is to learn a
robust regression model in the presence of missing values, it can be used as
an imputation tool as well. We run RIFLE and several state-of-the-art
approaches on five datasets from UCI repository (Spam, Housing, Clouds, Breast
Cancer, and Parkinson) with different proportions of MCAR missing values.
Then, we compute the NMRSE of imputed entries. Table 1 shows the performance
of RIFLE compared to other approaches for the datasets where the proportion of
missing values is relatively high
$\Big{(}\frac{n(1-p)}{d}\approx\mathcal{O}(1)\Big{)}$. RIFLE outperforms these
methods in almost all cases, and has a slightly better performance compared to
MissForest, which uses a highly non-linear model (random forest) to impute
missing values.

| Dataset Properties | NRMSE of Methods
---|---|---
Dataset | Missing Percentage | $n$ | $d$ | $n(1-p)/d$ | RIFLE | MICE | Amelia | MissForest | KNN | MIDA | Mean
Spam | 30% | 4602 | 57 | 56.51 | 0.87 | 1.23 | 1.26 | 0.90 | 1.04 | 0.97 | 1.02
Housing | 70% | 100 | 14 | 2.14 | 0.92 | 1.22 | 1.18 | 0.93 | 1.04 | 1.02 | 1.01
Clouds | 70% | 100 | 10 | 3 | 0.87 | 1.16 | 1.19 | 0.81 | 1.00 | 0.93 | 1.01
Breast Cancer | 40% | 198 | 34 | 3.49 | 0.54 | 0.77 | 0.83 | 0.57 | 0.79 | 0.72 | 0.99
Parkinson | 60% | 197 | 23 | 3.42 | 0.63 | 0.81 | 0.78 | 0.69 | 0.77 | 0.73 | 0.99

Table 1: Performance comparison of RIFLE and state-of-the-art methods on
several UCI datasets.

###  5.6 Sensitivity of Robust Imputer to Number of Samples and Proportion of
Missing Values

In this section we analyze the sensitivity of RIFLE and other state-of-the-art
approaches to the number of samples and proportion of missing values.

Figure 7: Performance Comparison of RIFLE, MICE and MissForest on four UCI
datasets: Parkinson, Spam, Wave Energy Converter, and Breast Cancer. For each
dataset we count the number of features that each method has a better
performance over other two methods.

In the first experiment, we create $5$ datasets containing $40\%$, $50\%$,
$60\%$, $70\%$, and $80\%$ of MCAR missing values respectively for four real
datasets (Spam, Parkinson, Wave Energy Converter, and Breast Cancer) from UCI
Repository. Given a feature of a dataset containing missing values, we say an
imputer wins that feature if the imputation error in terms of NRMSE for that
imputer is less than the error of other two imputers. Figure 7 reports the
number of features won by each one of the imputers on the created datasets
described above. As we can observe, the number of wins for RIFLE increases as
we increase the proportion of missing values in each one of four datasets.
This observation shows that sensitivity of RIFLE as an imputer to the
proportion of missing values is less than MissForest and MICE in general.

Figure 7 does not show how the NRMSE of one imputer is changed when the
proportion of missing values is increased. In the next experiment, we analyze
the sensitivity of RIFLE and several imputers in terms of NMRSE to change of
missing value proportions. Fixing the proportion of missing values, we
generate $10$ random datasets containing missing values in random locations on
Drive dataset (description of datasets is available in Appendix A). For each
dataset we impute the missing values with RIFLE, MissForest, Mean Imputation,
and MICE. Figure 8 shows the average and the standard deviation of these $4$
imputers’ performances for different proportions of missing values ($10\%$ to
$90\%$). Figure 8 depicts the sensitivity of MissForest and RIFLE to the
proportion of missing value in Drive dataset. For each experiment, we select
400 data points with different proportion of missing values (from $10\%$ to
$90\%$) and report the average NRMSE of imputed entries.

Figure 8: Sensitivity of RIFLE, MissForest, , Amelia, KNN Imputer, MIDA and
Mean Imputer to the percentage of missing values on Drive dataset. Increasing
the percentage of missing value entries degrades the performance of the
benchmarks as compared to RIFLE. KNN-imputer implementation cannot be executed
on datasets containing $80\%$ (or more) missing entries. Moreover, Amelia and
MIDA do not converge to a solution when the percentage of missing value
entries is higher than $70\%$.

###  5.7 Performance Comparison on Real Datasets

In this section, we compare the performance of RIFLE to several state-of-the-
art approaches including MICE (Buuren and Groothuis-Oudshoorn, 2010), Amelia
(Honaker et al., 2011), MissForest (Stekhoven and Bühlmann, 2012), KNN Imputer
(Raghunathan et al., 2001) and MIDA (Gondara and Wang, 2018).

There are two primary ways to do this. One method to predict a continuous
target variable in a dataset with many missing values is to first impute the
missing data with a state-of-the-art package, then run linear regression. An
alternative approach is to directly learn the target variable as we discussed
in Section 3.

Table 2 compares the performance of mean imputation, MICE, MIDAS, MissForest,
and KNN to that of RIFLE on three datasets: NHANES, Blog Feedback, and super
conductivity. Both Blog Feedback and Super Conductivity datasets contain
$30\%$ of MNAR missing values generated by Algorithm 10, with $10000$ and
$20000$ training samples respectively. The description of the NHANES data and
its distribution of missing values can be found in Appendix A.

Figure 9: Sensitivity of RIFLE, MissForest, MICE, Amelia, Mean Imputer, KNN
Imputer, and MIDA to the number of samples for the imputations of Drive
dataset containing $40\%$ of MCAR missing values. When the number of samples
is limited, RIFLE demonstrates better performance compared to other methods
and its performance is very close to the non-linear imputer MissForest for
larger samples.

Methods | Datasets
---|---
Super Conductivity | Blog Feedback | NHANES
Regression on Complete Data | 0.4601 | 0.7432 | 0.6287
RIFLE | 0.4894 | 0.8311 | 0.6262
Mean Imputation + Regression | 0.6114 | 0.9235 | 0.6287
MICE Imputation + Regression | 0.5094 | 0.8507 | 0.6612
MIDA Imputation + Regression | 0.5213 | 0.8394 | 0.6513
MissForest | 0.4911 | 0.8193 | 0.6365
KNN Imputer | 0.5436 | 0.8842 | 0.6401

Table 2: Normlaized RMSE of RIFLE and several state-of-the-art Methods on
Super Conductivity, Blog Feedback and NHANES datasets. The first two datasets
contain $30\%$ Missing Not At Random (MNAR) missing values in the training
phase generated by Algorithm 10.

Since MICE and MIDA are not capable of predicting values during the test phase
without data imputation, we use them as an imputer and then we apply the
linear regression to the imputed dataset. On the other hand, RIFLE, KNN
imputer and MissForest can predict the target variable without imputing the
train dataset. As we can observe in Table 2 RIFLE outperforms all other state-
of-the-art approaches executed on the three mentioned datasets. In particular,
the performance of RIFLE is higher than MissForest while the underlying model
used by RIFLE is simpler (linear) compared to the model (random forest)
utilized by Missforest.

###  5.8 Robust Normal Discriminant Analysis

In Section 4 we discussed how to specialize RIFLE to robust normal
discriminant analysis in presence of missing values. Since the maximization
problem with respect to the second-moment of the data ($\Sigma$) is
intractable, we solved the maximization problem over a set of $k$ covariance
matrices estimated by bootstrap sampling. To investigate the effect of
choosing $k$ on the performance of the robust classifier, we train robust
normal discriminant analysis models for different values of $k$ on two
training datasets (Avila and Magic) containing $40\%$ MCAR missing values. The
description of datasets can be found in Appendix A. For $k=1$ there is no
maximization problem, and thus, it is equivalent to the classifier proposed in
(Fung and Wrobel, 1989). As it can be observed in Figure 10, increasing the
number of covariance estimations, generally enhance the accuracy of the
classifier in the test phase. However, as it is shown in Theorem 11, the
required time for completing the training phase grows linearly in terms of the
number of covariance estimation.

Figure 10: Effect of the number of covariance estimations on the performance
(left) and run time (right) of robust LDA on Avila and Magic datasets.
Increasing the number of covariance estimations ($k$) improves the accuracy of
the model on the test data, however, it takes longer training time.

Number of Training Data Points | Method
---|---
LDA | Robust LDA (Common Covariance) | Robust LDA
50 | 52.38% | 62.14% | 61.36%
100 | 61.24% | 68.46% | 70.07%
200 | 73.49% | 73.35% | 73.51%

Table 3: Sensitivity of Linear Discriminant Analysis, Robust LDA with Common
Covariance Matrices for both groups, and Robust LDA (Different Covariance
matrices for two groups) to the number of training samples.

###  5.9 Comparison of Robust Linear Regression and Robust LDA

We presented two different approaches for learning a classification task in
the presence of missing values in Section 4: Thresholding on robust linear
regression and robust linear discriminant analysis (LDA). Table 4 compares the
performance of two classifiers on three different datasets. As it is
demonstrated in the table, when all features are continuous, linear
discriminant analysis has a better performance. It shows, the LDA model is
highly reliant on the normality assumption, while robust linear regression
handles the categorical features better than robust LDA.

|  | Accuracy of Methods
---|---|---
Dataset | Feature Type | RIFLE | Robust LDA | MissForest | MICE | KNN Imputer
Glass Identification | Continuous | 67.12% | 69.54% | 65.76% | 62.48% | 60.37%
Annealing | Mixed | 63.41% | 59.51% | 64.91% | 60.66% | 57.44%
Abalone | Mixed | 68.41% | 63.27% | 69.40% | 63.12% | 62.43%
Lymphography | Discrete | 65.32% | 58.15% | 66.11% | 55.73% | 57.39%
Adult | Discrete | 72.42% | 60.36% | 70.34% | 63.30% | 60.14%

Table 4: Accuracy of different approaches. When the input features are
continuous and the target variable is discrete, Robust LDA can perform better
than other methods. However, in the mixed and discrete settings, RIFLE results
in higher accuracy.

##  6 Conclusions

In this work, we have developed a novel robust framework for statistical
inference on incomplete datasets using the first and second-order moments of
the available data. RIFLE fits a statistical model for the worst-case
distribution such that the mean and variance fall within bootstrapped
confidence intervals. This framework can be widely applied to a variety of
conventional learning methods, and we provided theoretical convergence
guarantees, insights into the iteration complexity, and statistical properties
of the associated algorithms for linear regression and normal discriminant
analysis.

While we have designed the robust inference framework for direct statistical
inference in the presence of missing values without performing data
imputation, it can be adopted as an imputation tool as well. We have
demonstrated the performance of this method in direct inference and imputation
against several widely used packages. Our method generally outperforms these
existing methods when the number of missing entries is large.

While much work remains in the area of imputation and inference in the
presence of missing data, we believe that this framework will substantially
contribute to the field. In particular, we hope that this framework can be
used to help synthesize data to increase statistical power and more fully take
advantage of an increasingly large body of publicly available, if difficult to
combine, information.

## Acknowledgement

We would like to acknowledge support for this project from the NSF/NIH Joint
Initiative on Generalizable Data Science Methods for Biomedical Research
(1R01LM013315-01).

## References

  * Assländer et al. [2018] J. Assländer, M. A. Cloos, F. Knoll, D. K. Sodickson, J. Hennig, and R. Lattanzi.  Low rank alternating direction method of multipliers reconstruction for mr fingerprinting.  _Magnetic resonance in medicine_ , 79(1):83–96, 2018.
  * Barazandeh and Razaviyayn [2020] B. Barazandeh and M. Razaviyayn.  Solving non-convex non-differentiable min-max games using proximal gradient method.  _arXiv preprint arXiv:2003.08093_ , 2020.
  * Beaulieu-Jones et al. [2018] B. K. Beaulieu-Jones, D. R. Lavage, J. W. Snyder, J. H. Moore, S. A. Pendergrass, and C. R. Bauer.  Characterizing and managing missing structured data in electronic health records: data analysis.  _JMIR medical informatics_ , 6(1):e11, 2018.
  * Bertsimas et al. [2017] D. Bertsimas, C. Pawlowski, and Y. D. Zhuo.  From predictive methods to missing data imputation: an optimization approach.  _The Journal of Machine Learning Research_ , 18(1):7133–7171, 2017.
  * Bubeck [2014] S. Bubeck.  Convex optimization: Algorithms and complexity.  _arXiv preprint arXiv:1405.4980_ , 2014.
  * Buuren and Groothuis-Oudshoorn [2010] S. v. Buuren and K. Groothuis-Oudshoorn.  mice: Multivariate imputation by chained equations in r.  _Journal of statistical software_ , pages 1–68, 2010.
  * Cai et al. [2006] Z. Cai, M. Heydari, and G. Lin.  Iterated local least squares microarray missing value imputation.  _Journal of bioinformatics and computational biology_ , 4(05):935–957, 2006.
  * Dempster et al. [1977] A. P. Dempster, N. M. Laird, and D. B. Rubin.  Maximum likelihood from incomplete data via the em algorithm.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 39(1):1–22, 1977.
  * Fung and Wrobel [1989] K. Y. Fung and B. A. Wrobel.  The treatment of missing values in logistic regression.  _Biometrical Journal_ , 31(1):35–47, 1989.
  * Gabay and Mercier [1976] D. Gabay and B. Mercier.  A dual algorithm for the solution of nonlinear variational problems via finite element approximation.  _Computers & mathematics with applications_, 2(1):17–40, 1976.
  * Ghahramani and Jordan [1994] Z. Ghahramani and M. I. Jordan.  Supervised learning from incomplete data via an em approach.  In _Advances in neural information processing systems_ , pages 120–127, 1994.
  * Gondara and Wang [2018] L. Gondara and K. Wang.  MIDA: Multiple imputation using denoising autoencoders.  In _Pacific-Asia Conference on Knowledge Discovery and Data Mining_ , pages 260–272. Springer, 2018.
  * He and Yuan [2015] B. He and X. Yuan.  On non-ergodic convergence rate of douglas–rachford alternating direction method of multipliers.  _Numerische Mathematik_ , 130(3):567–577, 2015\.
  * Honaker et al. [2011] J. Honaker, G. King, and M. Blackwell.  Amelia ii: A program for missing data.  _Journal of statistical software_ , 45(7):1–47, 2011.
  * Hong et al. [2016] M. Hong, Z.-Q. Luo, and M. Razaviyayn.  Convergence analysis of alternating direction method of multipliers for a family of nonconvex problems.  _SIAM Journal on Optimization_ , 26(1):337–364, 2016.
  * Kim et al. [2005] H. Kim, G. H. Golub, and H. Park.  Missing value estimation for dna microarray gene expression data: local least squares imputation.  _Bioinformatics_ , 21(2):187–198, 2005.
  * Little and Rubin [2019] R. J. Little and D. B. Rubin.  _Statistical analysis with missing data_ , volume 793.  John Wiley & Sons, 2019.
  * Nesterov [1983] Y. E. Nesterov.  A method for solving the convex programming problem with convergence rate o (1/k^ 2).  In _Dokl. akad. nauk Sssr_ , volume 269, pages 543–547, 1983.
  * [19] A. Pedersen, E. Mikkelsen, D. Cronin-Fenton, N. Kristensen, T. Pham, L. Pedersen, et al.  Missing data and multiple imputation in clinical epidemiological research. clin epidemiol. 2017; 9: 157-66.
  * Raghunathan et al. [2001] T. E. Raghunathan, J. M. Lepkowski, J. Van Hoewyk, P. Solenberger, et al.  A multivariate technique for multiply imputing missing values using a sequence of regression models.  _Survey methodology_ , 27(1):85–96, 2001.
  * Schafer [1997] J. L. Schafer.  _Analysis of incomplete multivariate data_.  Chapman and Hall/CRC, 1997.
  * Shivaswamy et al. [2006] P. K. Shivaswamy, C. Bhattacharyya, and A. J. Smola.  Second order cone programming approaches for handling missing and uncertain data.  _Journal of Machine Learning Research_ , 7(Jul):1283–1314, 2006.
  * Sion et al. [1958] M. Sion et al.  On general minimax theorems.  _Pacific Journal of mathematics_ , 8(1):171–176, 1958.
  * Stekhoven and Bühlmann [2012] D. J. Stekhoven and P. Bühlmann.  Missforest—non-parametric missing value imputation for mixed-type data.  _Bioinformatics_ , 28(1):112–118, 2012.
  * Sterne et al. [2009] J. A. Sterne, I. R. White, J. B. Carlin, M. Spratt, P. Royston, M. G. Kenward, A. M. Wood, and J. R. Carpenter.  Multiple imputation for missing data in epidemiological and clinical research: potential and pitfalls.  _Bmj_ , 338:b2393, 2009.
  * Troyanskaya et al. [2001] O. Troyanskaya, M. Cantor, G. Sherlock, P. Brown, T. Hastie, R. Tibshirani, D. Botstein, and R. B. Altman.  Missing value estimation methods for dna microarrays.  _Bioinformatics_ , 17(6):520–525, 2001.
  * Xu et al. [2009] H. Xu, C. Caramanis, and S. Mannor.  Robustness and regularization of support vector machines.  _Journal of machine learning research_ , 10(7), 2009.
  * Yoon et al. [2018] J. Yoon, J. Jordon, and M. Van Der Schaar.  Gain: Missing data imputation using generative adversarial nets.  _arXiv preprint arXiv:1806.02920_ , 2018.
  * Zhang et al. [2018] T. Zhang, S. Ye, K. Zhang, J. Tang, W. Wen, M. Fardad, and Y. Wang.  A systematic dnn weight pruning framework using alternating direction method of multipliers.  In _Proceedings of the European Conference on Computer Vision (ECCV)_ , pages 184–199, 2018.
  * Zhang et al. [2008] X. Zhang, X. Song, H. Wang, and H. Zhang.  Sequential local least squares imputation estimating missing value of microarray data.  _Computers in biology and medicine_ , 38(10):1112–1120, 2008.

## Appendix A.

In this section, we introduce the datasets used in Section 5 to evaluate the
performance of RIFLE on regression and classification tasks. Note that, except
for the NHANES, all other datasets contain no missing values. For those
datasets we generate MCAR and MNAR missing values artificially (for MNAR
patterns we apply Algorithm 10 to the datasets).

### Datasets for Evaluating RIFLE on Regression and Imputation Tasks

  * •

NHANES: The percentage of missing values varies for different features of
NHANES dataset. There are two sources of missing values in NHANES data:
Missing entries during the collection of data, and missing entries resulted by
merging different datasets in the NHANES collection. On average, approximately
$20\%$ of data is missing.

  * •

Super
Conductivity222https://archive.ics.uci.edu/ml/datasets/Superconductivty+Data:
Super Conductivity datasets contains $21263$ samples describing
superconductors and their relevant features ($81$ attributes). All features
are continuous and the assigned task is to predict the critical temperature
based on the given $81$ features. We have used this dataset in experiments
summarized in Figure 3, Figure 4, and Table 2.

  * •

BlogFeedback333https://archive.ics.uci.edu/ml/datasets/BlogFeedback:
BlogFeedback data is a collection of $280$ features extracted from HTML-
documents of the blog posts. The assigned task is to predict the number of
comments in the upcoming $24$ hours based on the features of more than $60K$
data training data points. The test dataset is fixed and is originally
separated from the training data. The dataset is used in experiments described
in Table 2.

  * •

Breast
Cancer(Prognostic)444https://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Wisconsin+(Prognostic):
The dataset consists of $34$ features and $198$ instances. Each record
represents follow-up data for one breast cancer case collected in $1984$. We
have done several experiments to impute the MCAR missing values generated
artificially with different proportions. The results are depicted in Table 1
and Figure 7.

  * •

Parkinson555https://archive.ics.uci.edu/ml/datasets/parkinsons: The dataset
describes a range of biomedical voice recording from $31$ people, $23$ with
Parkinson’s disease (PD). The assigned task is to discriminate healthy people
from those with PD. There are $193$ records and $23$ features in the dataset.
The datasets is processed similar to the Breast Cancer dataset and used in the
same experiments.

  * •

Spam666https://archive.ics.uci.edu/ml/datasets/spambase: The dataset consists
of $4601$ instances and $57$ attributes. The assigned classification task is
to predict whether the email is spam. To evaluate different imputation method,
we randomly mask a proportion of data entries and impute them with different
approaches. The results are depicted in Table 1 and Figure 7.

  * •

Housing777https://www.kaggle.com/c/boston-housing: Boston Housing dataset
contains $506$ instances and $14$ columns. We generate random missing entries
with different proportions and impute it with RIFLE and several state-of-the-
art approaches. The results are demonstrated in Table 1 and Figure 7.

  * •

Clouds888https://archive.ics.uci.edu/ml/datasets/Cloud: The dataset has $1024$
instances and $10$ features extracted from clouds images. We use this dataset
in experiments depicted in Table 1 with $70\%$ artificial MCAR missing values.

  * •

Wave Energy
Converter999https://archive.ics.uci.edu/ml/datasets/Wave+Energy+Converters: We
sample a subset of $3000$ instances with $49$ features from the original Wave
Energy Converter dataset. We have executed several imputation methods on the
dataset and the results are shown in Figure 7.

  * •

Drive
Dataset101010https://archive.ics.uci.edu/ml/datasets/dataset+for+sensorless+drive+diagnosis:
The $49$ continuous features in this dataset are extracted from electric
current drive signals and the associated classification task is to determine
the condition of device’s motor. We choose different random samples with size
$400$ to run experiments (imputation) in Figure 8.

### Datasets for Evaluating Robust LDA on Classification Tasks

  * •

Avila111111https://archive.ics.uci.edu/ml/datasets/Avila: The Avila dataset
consists of $10$ attributes extracting from $800$ images of ”Avila Bible”. The
associated classification task is to match each pattern (an instance of the
dataset) to a copyist. We put $40\%$ of MCAR missing values (both input
features and the target variable) on $10$ different random samples of the
dataset with size $1000$. The average accuracy of the robust LDA method on the
$10$ datasets is demonstrated on Figure 10 for each value of $k$ (the number
of covariance estimations).

  * •

Magic121212https://archive.ics.uci.edu/ml/datasets/MAGIC+Gamma+Telescope: The
dataset consists of $11$ continuous MC-generated features from contributing to
the prediction of the type of event (signal or background). We used the same
procedure as the above dataset for the results depicted in Figure 10 (sampling
random subsets of $1000$ data points out of more than $19000$).

  * •

Glass
Identification131313https://archive.ics.uci.edu/ml/datasets/glass+identification:
This dataset is composed of $10$ continuous features and $214$ instances. The
assigned classification task is to predict the type of the glasses based on
the materials used for making them. We have assigned $40\%$ of MCAR missing
values to the dataset for the experiments reported in Table 4.

  * •

Annealing141414https://archive.ics.uci.edu/ml/datasets/Annealing This dataset
is a mix of categorical and numerical features ($37$ in total), and the
associated task is to predict the class (5 classes) of instances (metals). The
number of instances in this dataset is $798$. We use $500$ data points as the
training data and the rest as the test. We apply $40\%$ of MCAR missing values
to both input features and the target. The accuracy of different models are
reported in Table 4.

  * •

Abalone151515https://archive.ics.uci.edu/ml/datasets/abalone: This dataset
consists of $4177$ instances and $8$ categorical and continuous features. The
goal is to predict the age of abalone based on the physical measurements. The
first $1000$ samples are used as the training data and the rest as the test
data. We applied the same pre-processing procedure as the above dataset for
generating missing values on the training data. he accuracy of different
models are reported in Table 4.

  * •

Lymphography161616https://archive.ics.uci.edu/ml/datasets/Lymphography:
Lymphography is a categorical dataset containing $18$ features and $148$ data
points obtained from the University Medical Centre, Institute of Oncology,
Ljubljana, Yugoslavia. $100$ data points are used as the training data, and
the rest are test data points (with no missing values). We applied the same
pre-processing described for the above dataset to generate MCAR missing
values.

  * •

Adult 171717https://archive.ics.uci.edu/ml/datasets/adult.: Adult dataset
contains the census information of individuals including education, gender,
and capital gain. The assigned classification task is to predict whether a
person earns over 50k annually. The train and test sets are two separated
files consisting of $32,000$ and $16,000$ samples respectively. We consider
gender and race as the sensitive attributes (For the experiments involving one
sensitive attribute, we have chosen gender). Learning a logistic regression
model on the training dataset (without imposing fairness) shows that only 3
features out of 14 have larger weights than the gender attribute.

## Appendix B.

In this section, we prove all lemmas and theorems presented in the main body
of the article in order.

Proof of Lemma 2 Since, the problem is convex in ${\bm{\theta}}$ and concave
in ${\mathbf{C}}$ and ${\mathbf{b}}$, we have:

$\displaystyle\min_{{\bm{\theta}}}\quad\max_{{\mathbf{C}},{\mathbf{b}}}\quad{\bm{\theta}}^{T}{\mathbf{C}}{\bm{\theta}}-2{\mathbf{b}}^{T}{\bm{\theta}}+\lambda\|\theta\|^{2}=-\min_{{\mathbf{C}},{\mathbf{b}}}\quad\max_{{\bm{\theta}}}\quad-{\bm{\theta}}^{T}{\mathbf{C}}{\bm{\theta}}+2{\mathbf{b}}^{T}{\bm{\theta}}-\lambda\|{\bm{\theta}}\|^{2}$

The new objective function is convex in ${\mathbf{C}}$ and ${\mathbf{b}}$, and
strongly concave with respect to ${\bm{\theta}}$. According to Lemma 1 in
[Barazandeh and Razaviyayn, 2020],
$g^{{}^{\prime}}=-g=\max_{{\bm{\theta}}}-{\bm{\theta}}^{T}{\mathbf{C}}{\bm{\theta}}+2{\mathbf{b}}^{T}{\bm{\theta}}-\lambda\|{\bm{\theta}}\|^{2}$
is Lipschitz continuous with the Lipschitz constant equal to:

$L_{g}=L_{g^{{}^{\prime}}}=L_{11}+\frac{L_{12}}{\sigma},$

where $\sigma$ is the strong-concavity modulus of
$-{\bm{\theta}}^{T}{\mathbf{C}}{\bm{\theta}}+2{\mathbf{b}}^{T}{\bm{\theta}}-\lambda\|{\bm{\theta}}\|^{2}$,
and $L_{11}$, $L_{12}$ are defined as follows:

It is easy to see that $L_{11}=0$ and $L_{12}=2\max\|{\bm{\theta}}\|_{2}\leq
2$, which completes the proof.

Proof of Theorem 3 Since the set of feasible solutions for ${\mathbf{b}}$ and
${\mathbf{C}}$ defines a compact set, and function $g$ is a concave function
with respect to ${\mathbf{b}}$ and ${\mathbf{C}}$, the projected gradient
ascent algorithm converges to the global maximizer of $g$ in
$T=\mathcal{O}(\frac{LD^{2}}{\epsilon})$ iterations [Bubeck, 2014], where
$D=\|{\mathbf{C}}^{*}-\hat{{\mathbf{C}}}\|_{F}$ and $L$ is the Lipschitz
constant of function $g$, which is equal to $\frac{2}{\lambda}$ according to
Lemma 2.

Proof of Theorem 4 Algorithm 4 applies the projected Nesterov acceleration
method on the concave function $g$. As it is proved in [Nesterov, 1983], the
rate of convergence of this method conforms to the lower bound of first order
oracles for the general convex minimization (concave maximization) problems,
which is ${\cal O}(\sqrt{\frac{LD^{2}}{\epsilon}})$. The Lipschitz constant
$L$ appeared in the iteration complexity bound is computed by Lemma 2.

Proof of Theorem 5 First, note that if we multiply the objective function by
$-1$, Problem (9) can be equivalently formulated as:

$\begin{array}[]{lll}\displaystyle{\max_{{\bm{\theta}}}\quad\min_{{\mathbf{C}},{\mathbf{b}}}}&\displaystyle{-{\bm{\theta}}^{T}{\mathbf{C}}{\bm{\theta}}+2{\mathbf{b}}^{T}{\bm{\theta}}-\lambda\|{\bm{\theta}}\|_{2}^{2}}\\\
{\rm s.t.}&-{\mathbf{C}}+{\mathbf{C}}_{min}\leq 0,\\\
&{\mathbf{C}}-{\mathbf{C}}_{max}\leq 0,\\\
&-{\mathbf{b}}+{\mathbf{b}}_{min}\leq 0,\\\
&{\mathbf{b}}-{\mathbf{b}}_{max}\leq 0,\\\ &-{\mathbf{C}}\preceq 0\end{array}$
(36)

If we assign
${\mathbf{A}},{\mathbf{B}},{\mathbf{d}},{\mathbf{e}},{\mathbf{H}}$ to the
constraints respectively, then the Lagrangian function can be written as:

$\begin{array}[]{lll}\displaystyle{L({\mathbf{C}},{\mathbf{b}},{\mathbf{A}},{\mathbf{B}},{\mathbf{d}},{\mathbf{e}},{\mathbf{H}})}=&\displaystyle{-{\bm{\theta}}^{T}{\mathbf{C}}{\bm{\theta}}+2{\mathbf{b}}^{T}{\bm{\theta}}+\langle{\mathbf{A}},-{\mathbf{C}}+{\mathbf{C}}_{min}\rangle}\\\
&+\langle{\mathbf{B}},{\mathbf{C}}-{\mathbf{C}}_{max}\rangle+\langle{\mathbf{d}},-{\mathbf{b}}+{\mathbf{b}}_{min}\rangle\\\
&+\langle{\mathbf{e}},{\mathbf{b}}-{\mathbf{b}}_{max}\rangle-\langle{\mathbf{C}},{\mathbf{H}}\rangle-\lambda\|{\bm{\theta}}\|_{2}^{2},\end{array}$
(37)

The dual problem is defined as:

$\begin{array}[]{lll}\displaystyle{\max_{{\mathbf{A}},{\mathbf{B}},{\mathbf{d}},{\mathbf{e}},{\mathbf{H}}}\min_{{\mathbf{C}},{\mathbf{b}}}L({\mathbf{C}},{\mathbf{b}},{\mathbf{A}},{\mathbf{B}},{\mathbf{d}},{\mathbf{e}},{\mathbf{H}})}\end{array}$
(38)

The minimization of $L$ takes the following form:

$\begin{array}[]{lll}\displaystyle{\min_{{\mathbf{C}},{\mathbf{b}}}}&\langle{\mathbf{C}},-{\bm{\theta}}{\bm{\theta}}^{T}-{\mathbf{A}}+{\mathbf{B}}-{\mathbf{H}}\rangle+\langle{\mathbf{b}},2{\bm{\theta}}-{\mathbf{d}}+{\mathbf{e}}\rangle-\lambda\|{\bm{\theta}}\|_{2}^{2}\\\
&-\langle{\mathbf{B}},{\mathbf{C}}_{max}\rangle+\langle{\mathbf{A}},{\mathbf{C}}_{min}\rangle-\langle{\mathbf{e}},{\mathbf{b}}_{max}\rangle+\langle{\mathbf{d}},{\mathbf{b}}_{min}\rangle,\end{array}$
(39)

To avoid $-\infty$ value for the above minimization problem, it is required to
set $-{\bm{\theta}}{\bm{\theta}}^{T}-{\mathbf{A}}+{\mathbf{B}}-{\mathbf{H}}$
and $2{\bm{\theta}}-{\mathbf{d}}+{\mathbf{e}}$ to zero. Thus the dual problem
of (36) is formulated as:

$\begin{array}[]{lll}\displaystyle{\max_{{\mathbf{A}},{\mathbf{B}},{\mathbf{d}},{\mathbf{e}},{\mathbf{H}}}}&\displaystyle{{\mathbf{b}}_{min}^{T}{\mathbf{d}}-{\mathbf{b}}_{max}^{T}{\mathbf{e}}+\langle{\mathbf{C}}_{min},{\mathbf{A}}\rangle-\langle{\mathbf{C}}_{max},{\mathbf{B}}\rangle-\lambda\|{\bm{\theta}}\|_{2}^{2}}\\\
{\rm
s.t.}&-{\bm{\theta}}{\bm{\theta}}^{T}-{\mathbf{A}}+{\mathbf{B}}-{\mathbf{H}}=0,\\\
&2{\bm{\theta}}-{\mathbf{d}}+{\mathbf{e}}=0,\\\
&{\mathbf{A}},{\mathbf{B}},{\mathbf{d}},{\mathbf{e}}\geq 0,\\\
&{\mathbf{H}}\succeq 0\end{array}$ (40)

Since the duality gap is zero, Problem (9) can be equivalently formulated as:

$\begin{array}[]{lll}\displaystyle{\max_{{\bm{\theta}},{\mathbf{A}},{\mathbf{B}},{\mathbf{d}},{\mathbf{e}},{\mathbf{H}}}}&\displaystyle{{\mathbf{b}}_{min}^{T}{\mathbf{d}}-{\mathbf{b}}_{max}^{T}{\mathbf{e}}+\langle{\mathbf{C}}_{min},{\mathbf{A}}\rangle-\langle{\mathbf{C}}_{max},{\mathbf{B}}\rangle-\lambda\|{\bm{\theta}}\|^{2}}\\\
{\rm
s.t.}&-{\bm{\theta}}{\bm{\theta}}^{T}-{\mathbf{A}}+{\mathbf{B}}-{\mathbf{H}}=0,\\\
&2{\bm{\theta}}-{\mathbf{d}}+{\mathbf{e}}=0,\\\
&{\mathbf{A}},{\mathbf{B}},{\mathbf{d}},{\mathbf{e}}\geq 0,\\\
&{\mathbf{H}}\succeq 0.\end{array}$ (41)

We can multiply the objective function by $-1$ and change the minimization to
maximization, which gives the dual problem described in (13).

Proof of Theorem 8

(a) Let $\bm{\Delta}_{n}$ be the estimated confidence matrix obtained from $n$
samples. The first part of the theorem is true, if we show that
$\bm{\Delta}_{n}$ converges to $0$ as $n$, the number of samples goes to
infinity. Assume that $\\{({\mathbf{x}}_{i1},{\mathbf{x}}_{i2})\\}_{i=1}^{n}$
is an i.i.d bootstrap sample over data points that both features
${\mathbf{X}}_{1}$ and ${\mathbf{X}}_{2}$ are available. Since the
distribution of missing values is completely at random (MCAR), we have
$\mathbb{E}[{\mathbf{x}}_{i1}{\mathbf{x}}_{i2}]=\mathbb{E}[{\mathbf{X}}_{1}{\mathbf{X}}_{2}]$.
Therefore,
$\mathbb{E}[\frac{1}{n}\sum_{i=1}^{n}{\mathbf{x}}_{i1}{\mathbf{x}}_{i2}]=\mathbb{E}[{\mathbf{X}}_{1}{\mathbf{X}}_{2}]$.
Moreover, since the samples are drawn independently,
$\textrm{Var}[\frac{1}{n}\sum_{i=1}^{n}{\mathbf{x}}_{i1}{\mathbf{x}}_{i2}]=\frac{1}{n^{2}}\sum_{i=1}^{n}\textrm{Var}[{\mathbf{x}}_{i1}{\mathbf{x}}_{i2}]=\frac{n}{n^{2}}\textrm{Var}[{\mathbf{X}}_{1}{\mathbf{X}}_{2}]=\frac{1}{n}\textrm{Var}[{\mathbf{X}}_{1}{\mathbf{X}}_{2}]$.
Since the variance of the product of every two features is bounded, as $n$
goes to $\infty$, the variance of the bootstrap estimations converges to $0$,
which means $\lim_{n\to\infty}\bm{\Delta}_{n}=0$.

(b) Fix two features $i$ and $j$. Let
$(\hat{{\mathbf{X}}}_{i1},\hat{{\mathbf{X}}}_{j1}),\dots,(\hat{{\mathbf{X}}}_{im},\hat{{\mathbf{X}}}_{jm})$
be $m=n(1-p)$ i.i.d pairs sampled via bootstrap from the entries where both
features $i$ and $j$ are available. Define $Z_{m}=\hat{X}_{im}\hat{X}_{jm}$
(for simplicity we do not consider the dependence of $Z_{m}$ to $i$ and $j$ in
the notation). Assume that we initialize
${\mathbf{C}}_{ij}=\frac{1}{m}\sum_{t=1}^{m}Z_{t}$. Note that,
$\mathbb{E}[Z_{t}]=\mathbb{E}[\hat{X}_{it}\hat{X}_{jt}]={\mathbf{C}}^{*}_{ij}$.
According to Hoeffding’s inequality, we have:

$\textrm{Pr}\Bigg{[}\quad\Bigg{|}\frac{1}{m}\sum_{t=1}^{m}\Big{(}Z_{t}-\mathbb{E}[Z_{t}]\Big{)}\Bigg{|}\geq\bm{\Delta}_{ij}\Bigg{]}\leq\exp(-\frac{2m\bm{\Delta}_{ij}^{2}}{C^{2}})$

Let $\bm{\Delta}=\min\\{\bm{\Delta}_{ij}\\}$. Then the above inequality holds
for all pairs of features:

$\textrm{Pr}\Bigg{[}\quad\Bigg{|}\hat{{\mathbf{C}}}_{ij}-\mathbb{E}[Z_{t}]\Bigg{|}\geq\bm{\Delta}\Bigg{]}\leq\exp\Big{(}-\frac{2n(1-p)\bm{\Delta}^{2}}{C^{2}}\Big{)}$

Using a union bound argument, with the probability of at least
$1-\frac{d^{2}}{2}\exp\Big{(}-\frac{2n(1-p)\bm{\Delta}^{2}}{C^{2}}\Big{)}$, we
have:
$\hat{{\mathbf{C}}}-\bm{\Delta}\leq{\mathbf{C}}^{*}\leq\hat{{\mathbf{C}}}+\bm{\Delta}$,
which means the actual covariance matrix is within the confidence intervals we
have considered. For
$(\tilde{{\bm{\theta}}},\tilde{{\mathbf{b}}},\tilde{{\mathbf{C}}})$, we have:

$\displaystyle\tilde{{\bm{\theta}}}^{T}\tilde{{\mathbf{C}}}\tilde{{\bm{\theta}}}-2\tilde{{\mathbf{b}}}^{T}\tilde{{\bm{\theta}}}=\max_{{\mathbf{C}},{\mathbf{b}}}\quad\tilde{{\bm{\theta}}}^{T}{\mathbf{C}}\tilde{{\bm{\theta}}}-2{\mathbf{b}}^{T}\tilde{{\bm{\theta}}}\geq\tilde{{\bm{\theta}}}^{T}{\mathbf{C}}^{*}\tilde{{\bm{\theta}}}-2{\mathbf{b}}^{*T}\tilde{{\bm{\theta}}}$

## Appendix C.

In this appendix we show how to ADMM schema to Problem (14) to obtain
Algorithm 6.

As we discussed in Section 3, we consider two separate blocks of variables
${\mathbf{w}}=({\bm{\theta}},{\mathbf{d}},{\mathbf{e}},{\mathbf{G}},{\mathbf{B}}^{{}^{\prime}},{\mathbf{A}}^{{}^{\prime}})$
and
${\mathbf{z}}=({\mathbf{d}}^{{}^{\prime}},{\mathbf{e}}^{{}^{\prime}},{\bm{\theta}}^{{}^{\prime}},{\mathbf{B}},{\mathbf{A}})$.
Assigning
$\bm{\Gamma},{\bm{\eta}},{\mathbf{M}}_{A},{\mathbf{M}}_{B},{\bm{\mu}}_{d},{\bm{\mu}}_{e},$
and ${\bm{\mu}}_{\theta}$ to the constraints of Problem (14) in order, we can
write the corresponding augmented Lagrangian function as:

$\begin{array}[]{lll}\displaystyle{\min_{{\bm{\theta}},{\bm{\theta}}^{{}^{\prime}},{\mathbf{A}},{\mathbf{A}}^{{}^{\prime}},{\mathbf{B}},{\mathbf{B}}^{{}^{\prime}},{\mathbf{d}},{\mathbf{d}}^{{}^{\prime}},{\mathbf{e}},{\mathbf{e}}^{{}^{\prime}},{\mathbf{G}}}}&\displaystyle{-\langle{\mathbf{b}}_{min},{\mathbf{d}}\rangle+\langle{\mathbf{b}}_{max},{\mathbf{e}}\rangle-\langle{\mathbf{C}}_{min},{\mathbf{A}}\rangle+\langle{\mathbf{C}}_{max},{\mathbf{B}}\rangle+\lambda\|{\bm{\theta}}\|^{2}}\\\
&+\langle{\mathbf{A}}-{\mathbf{A}}^{{}^{\prime}},{\mathbf{M}}_{A}\rangle+\frac{\rho}{2}\|{\mathbf{A}}-{\mathbf{A}}^{{}^{\prime}}\|_{F}^{2}\\\
&+\langle{\mathbf{B}}-{\mathbf{B}}^{{}^{\prime}},{\mathbf{M}}_{B}\rangle+\frac{\rho}{2}\|{\mathbf{B}}-{\mathbf{B}}^{{}^{\prime}}\|_{F}^{2}\\\
&+\langle{\mathbf{d}}-{\mathbf{d}}^{{}^{\prime}},{\bm{\mu}}_{d}\rangle+\frac{\rho}{2}\|{\mathbf{d}}-{\mathbf{d}}^{{}^{\prime}}\|^{2}\\\
&+\langle{\mathbf{e}}-{\mathbf{e}}^{{}^{\prime}},{\bm{\mu}}_{e}\rangle+\frac{\rho}{2}\|{\mathbf{e}}-{\mathbf{e}}^{{}^{\prime}}\|^{2}\\\
&+\langle{\bm{\theta}}-{\bm{\theta}}^{{}^{\prime}},{\bm{\mu}}_{\theta}\rangle+\frac{\rho}{2}\|{\bm{\theta}}-{\bm{\theta}}^{{}^{\prime}}\|^{2}\\\
&+\langle
2{\bm{\theta}}-{\mathbf{d}}+{\mathbf{e}},{\bm{\eta}}\rangle+\frac{\rho}{2}\|2{\bm{\theta}}-{\mathbf{d}}+{\mathbf{e}}\|^{2}\\\
&+\langle{\mathbf{B}}-{\mathbf{A}}-{\mathbf{G}},\bm{\Gamma}\rangle+\frac{\rho}{2}\|{\mathbf{B}}-{\mathbf{A}}-{\mathbf{G}}\|_{F}^{2}\\\
{\rm
s.t.}&{\mathbf{A}}^{{}^{\prime}},{\mathbf{B}}^{{}^{\prime}},{\mathbf{d}}^{{}^{\prime}},{\mathbf{e}}^{{}^{\prime}}\geq
0,\\\
&\mathbf{G}\succeq{\bm{\theta}}^{{}^{\prime}}{\bm{\theta}}^{{}^{\prime}T},\end{array}$
(42)

At each iteration of ADMM algorithm, the parameters of one block are fixed and
the optimization problem is solved with respect to the parameters of the other
block. For the simplicity of presentation, let
${\mathbf{c}}_{1}^{t}=\rho{\bm{\theta}}^{{}^{\prime}t}-{\bm{\mu}}_{{\bm{\theta}}}^{t}-2{\bm{\eta}}^{t},{\mathbf{c}}_{2}^{t}=\rho{\mathbf{d}}^{{}^{\prime}t}-{\bm{\mu}}_{{\mathbf{d}}}^{t}-{\mathbf{b}}_{min}+{\bm{\eta}}^{t},{\mathbf{c}}_{3}^{t}=\rho{\mathbf{e}}^{{}^{\prime}t}-{\bm{\mu}}_{e}^{t}+{\mathbf{b}}_{max}-{\bm{\eta}}^{t},{\mathbf{D}}_{1}^{t}=\rho{\mathbf{A}}^{{}^{\prime}t}-\rho{\mathbf{G}}^{t}+\bm{\Gamma}^{t}-{\mathbf{M}}_{A}^{t}+{\mathbf{C}}_{min},$
and
${\mathbf{D}}_{2}^{t}=\rho{\mathbf{B}}^{{}^{\prime}t}+\rho{\mathbf{G}}^{t}-\bm{\Gamma}^{t}-{\mathbf{M}}_{B}^{t}-{\mathbf{C}}_{max}$.

We have two non-trivial problems containing positive semi-definite
constraints. The sub-problem with respect to ${\mathbf{G}}$ can be written as:

$\begin{array}[]{lll}\displaystyle{\min_{{\mathbf{G}}}}&\displaystyle{\langle{\mathbf{B}}^{t}-{\mathbf{A}}^{t}-{\mathbf{G}},\bm{\Gamma}^{t}\rangle+\frac{\rho}{2}\|{\mathbf{B}}^{t}-{\mathbf{A}}^{t}-{\mathbf{G}}\|_{F}^{2}}\\\
{\rm
s.t.}&{\mathbf{G}}\succeq{\bm{\theta}}^{{}^{\prime}t}{\bm{\theta}}^{{}^{\prime}tT},\end{array}$
(43)

By completing the square, and changing the variable
${\mathbf{G}}^{{}^{\prime}}={\mathbf{G}}-{\bm{\theta}}^{{}^{\prime}t}{\bm{\theta}}^{{}^{\prime}tT}$,
equivalently we require to solve the following problem:

$\begin{array}[]{lll}\displaystyle{\min_{{\mathbf{G}}^{{}^{\prime}}}}&\displaystyle{\frac{\rho}{2}\|{\mathbf{G}}^{{}^{\prime}}-({\mathbf{B}}^{t}-{\mathbf{A}}^{t}-{\bm{\theta}}^{{}^{\prime}t}{\bm{\theta}}^{{}^{\prime}tT}+\frac{\bm{\Gamma}^{t}}{\rho})\|_{F}^{2}}\\\
{\rm s.t.}&\bm{G}^{{}^{\prime}}\succeq 0,\end{array}$ (44)

Thus,
${\mathbf{G}}^{{}^{\prime}*}=[{\mathbf{B}}^{t}-{\mathbf{A}}^{t}+\frac{\bm{\Gamma}^{t}}{\rho}-{\bm{\theta}}^{{}^{\prime}t}{\bm{\theta}}^{{}^{\prime}tT}]_{+}$,
where $[{\mathbf{A}}]_{+}$ is the projection to the set of PSD matrices, which
can be done by setting the negative eigenvalues of ${\mathbf{A}}$ in its
singular value decomposition to zero.

The other non-trivial sub-problem in Algorithm (6) is the minimization with
respect to ${\bm{\theta}}^{{}^{\prime}}$ (Line 10). By completing the square,
it can be equivalently formulated as:

$\begin{array}[]{lll}\displaystyle{\min_{{\bm{\theta}}^{{}^{\prime}}}}&\displaystyle{\|{\bm{\theta}}^{{}^{\prime}}-({\bm{\theta}}^{t+1}+\frac{{\bm{\mu}}_{{\bm{\theta}}}^{t}}{\rho})\|_{2}^{2}}\\\
{\rm
s.t.}&{\mathbf{G}}^{t+1}\succeq{\bm{\theta}}^{{}^{\prime}}{\bm{\theta}}^{{}^{\prime}T},\end{array}$
(45)

Let ${\mathbf{G}}=\bm{U}\bm{\Lambda}\bm{U}^{T}$ and
$\bm{\alpha}={\bm{\theta}}^{t+1}+\frac{{\bm{\mu}}_{{\bm{\theta}}}^{t}}{2}$.
Since $\bm{U}^{T}\bm{U}=\bm{I}$, we have:

$\displaystyle\|\bm{U}^{T}{\bm{\theta}}-\bm{U}^{T}\bm{\alpha}\|^{2}=\|{\bm{\theta}}-\bm{\alpha}\|_{2}^{2}$

Set $\bm{\beta}=\bm{U}^{T}{\bm{\theta}}^{{}^{\prime}}$, then Problem (45) can
be reformulated as:

$\begin{array}[]{lll}\displaystyle{\min_{\bm{\beta}}}&\displaystyle{\|\bm{\beta}-\bm{U}^{T}\bm{\alpha}\|_{2}^{2}}\\\
{\rm s.t.}&\bm{\beta}\bm{\beta}^{T}\preceq\bm{\Lambda}.\end{array}$ (46)

Note that the constraint of the above optimization problem is equivalent to:

$\displaystyle\bm{\beta}\bm{\beta}^{T}\preceq\bm{\Lambda}\Leftrightarrow\begin{bmatrix}1&\bm{\beta}^{T}\\\
\bm{\beta}&\bm{\Lambda}\end{bmatrix}\succeq
0\Leftrightarrow\bm{\beta}^{T}\bm{\Lambda}^{-1}\bm{\beta}\leq
1\Leftrightarrow\sum_{i=1}^{n}\frac{\bm{\beta}_{i}^{2}}{\lambda_{i}}\leq 1$

Since the block matrix is symmetric, using Schur Complement, it is positive
semi-definite if and only if $\bm{\Lambda}$ is positive semi definite and
$1-\bm{\beta}^{t}\bm{\Lambda}^{-1}\bm{\beta}\geq 0$ (the third inequality
above).

Set $\bm{\gamma}=\bm{U}^{T}\bm{\alpha}$, then we can write Problem (46) as:

$\begin{array}[]{lll}\displaystyle{\min_{\bm{\beta}}}&\displaystyle{\|\bm{\beta}-\bm{\gamma}\|_{2}^{2}}\\\
{\rm s.t.}&\sum_{i=1}^{n}\frac{\bm{\beta}_{i}^{2}}{\lambda_{i}}\leq
1,\end{array}$ (47)

It can be easily shown that, the optimal solution has the form
$\bm{\beta}^{*}_{i}=\frac{\bm{\gamma}_{i}}{1+\frac{{\bm{\mu}}^{*}}{\lambda_{i}}}$,
where ${\bm{\mu}}^{*}$ is the optimal Lagrangian multiplier corresponding to
the constraint of Problem (47). The optimal Lagrangian multiplier can be
obtained by the bisection algorithm similar to Algorithm 7. Having
$\bm{\beta}^{*}$, the optimal ${\bm{\theta}}$ can be computed by solving the
linear equation $\bm{U}^{T}{\bm{\theta}}^{*}=\bm{\beta}^{*}$.