# Building a COVID-19 Vulnerability Index

Dave DeCaprio
ClosedLoop.ai
Corresponding author<EMAIL_ADDRESS>Joseph Gartner
ClosedLoop.ai
Corresponding author<EMAIL_ADDRESS>Carol J. McCall, FSA, MPH
ClosedLoop.ai
Corresponding author<EMAIL_ADDRESS>Thadeus Burgess
ClosedLoop.ai
Kristian Garcia
Healthfirst
Sarthak Kothari
ClosedLoop.ai
Shaayaan Sayed
ClosedLoop.ai

## 1\. Introduction

### 1.1. COVID-19 Illness

Coronaviruses (CoV) are a large family of viruses that cause illnesses ranging
from the common cold to more severe diseases such as Middle East respiratory
syndrome (MERS-CoV) and severe acute respiratory syndrome (SARS-CoV). CoV are
zoonotic, meaning they are transmitted between animals and people. Coronavirus
disease 2019 (COVID-19) is caused by a new strain discovered in 2019, severe
acute respiratory syndrome coronavirus 2 (SARS-CoV-2), that has not been
previously identified in humans[1].

COVID-19 is a highly contagious respiratory infection with common signs that
include respiratory symptoms, fever, cough, shortness of breath, and breathing
difficulties. In more severe cases, infection can cause pneumonia, severe
acute respiratory syndrome, kidney failure, cardiac arrest, and death[17][18].
Experts continue to learn more about COVID-19, including its etiology,
symptoms, complications, and potential treatments.

### 1.2. Flattening the Curve

On March 11, 2020, the World Health Organization (WHO) declared COVID-19 to be
a pandemic[2]. Public health and healthcare experts agree that mitigation is
required in order to slow the spread of COVID-19 and prevent the collapse of
healthcare systems. Health systems in the United States run close to
capacity[3], and so every transmission that can be avoided and every case that
can be prevented has enormous impact.

### 1.3. Identifying Vulnerable People

The risk of severe complications from COVID-19 is higher for certain
vulnerable populations, particularly people who are elderly, frail, or have
multiple chronic conditions. The risk of death has been difficult to
calculate[4], but a small study[5] of people who contracted COVID-19 in Wuhan
along with patterns also seen in early reports from the United States suggest
that the risk of death increases with age, and is also higher for those who
have diabetes, heart disease, blood clotting problems, or have shown signs of
sepsis. With an average death rate of 1%, the death rate rose to 6% for people
with cancer, high blood pressure, or chronic respiratory disease, 7% for
people with diabetes, and 10% for people with heart disease. There was also a
steep age gradient; the death rate among people age 80 and over was 15% [6].

Identifying who is most vulnerable is not necessarily straightforward. More
than 55% of Medicare beneficiaries meet at least one of the risk criteria
listed by the US Centers for Disease Control and Prevention (CDC)[7]. People
with the same chronic condition don’t have the same risk, and many people will
have other comorbidities, which compounds their vulnerability. Simple rules
can’t reflect these differences or capture complex factors like frailty[9]
which makes people more vulnerable to severe infections.

## 2\. Methods

### 2.1. Proxy Outcome

Since real-world data on COVID-19 cases are not readily available, the C-19
Index was developed using close proxy events. A person’s C-19 Index is
measured in terms of their near-term risk of severe complications from
respiratory infections (e.g. pneumonia, influenza). The most direct proxy
event is Acute respiratory distress syndrome (ARDS), identified by ICD-10
diagnosis code J80. ARDS is extremely rare, with an annual occurence of less
than 0.05% in Medicare members. To create a viable machine learning model, the
outcome was broadened to include 4 closely related categories of respiratory
diagnoses from the Clinical Classifications Software Refined (CCSR)[12]
classification system. Patient’s were considered to have the proxy event if
they had any of the following diagnosis codes in any position on a medical
claim associated with a hospital inpatient visit or observation stay:

  * •

ICD-10-CM J80 - Acute respiratory distress syndrome

  * •

RSP002 - Pneumonia (except that caused by tuberculosis)

  * •

RSP003 - Influenza

  * •

RSP005 - Acute bronchitis

  * •

RSP006 - Other specified upper respiratory infections

Machine learning models were created that predict the likelihood that a
patient will have an inpatient hospital stay due to one of the above
conditions in the next 3 months. While claims-based machine learning models
typically focus on longer term predictions such as outcomes with the next 6
months or one year, this model is expected to be used for immediate targeting
decisions around COVID-19, and so a shorter window was deemed appropriate.

### 2.2. Datasets

These models were trained and initially tested using historical medical claims
data. Two different data sets were obtained, representing different segments
of the population. The first was the Center for Medicare & Medicaid Services
(CMS) Limited Data Set (LDS) [14] for 2015 & 2016\. The LDS contains
beneficiary level health information for 5% of the Medicare population, and is
available to the public subject to a Data Use Agreement. The second was
medical claims data for 2.5 million beneficiaries obtained from Healthfirst,
which provides health insurance for New Yorkers. These two data sets
represented different demographics within the US population. The LDS data
contains only Medicare beneficiaries, who are predominantly over the age of 65
or disabled. The Healthfirst data set is primarily a Medicaid population,
which includes more healthy adults.

In order to build a predictive model appropriate for the overall US
population, these two data sets were combined into a single population that
has a demographic profile consistent with the overall US population. Cohorts
were individually created for each data set and then the resulting cohorts
were combined to create the final training and test sets for the models.

####  2.2.1 CMS Data Preparation

The CMS data cohort was created by identifying all living members above the
age of 65 on 9/30/2016. This particular date was chosen because it was 3
months from the end of the data set and was therefore the latest prediction
date possible in the data set. Using the latest date minimized the use of
ICD-9 data in the prediction histories since this data spanned the transition
from International Classification of Diseases version 9 to version 10 (ICD-10)
on October 1, 2015. Members under 65 were excluded because the second
population we had was a better representation of the under 65 population. Only
fee-for-service members were included because medical claims histories for
other members are not reliably complete. We then excluded all members who had
less than 6 months of continuous eligibility prior to 9/30/2016. We also
excluded members who lost coverage within 3 months after 9/30/2016, except for
those members who lost coverage due to death. Because these members lost
coverage, we cannot be confident that a negative label s Table 1 below
summarizes the population selection.

Population Size | Selection Criteria
---|---
3,114,713 | Total members in the CMS dataset
1,867,879 | Fee-for-service members with 6 months of continuous coverage prior to 9/30/2016
1,511,950 | 65 years old or older
1,506,659 | Exclude members who died before 9/30/2016
1,500,700 | Exclude members who lose coverage before 12/31/2016 not due to death.
Table 1: Cohort selection for the CMS population

####  2.2.2 Healthfirst Data Preparation

The Healthfirst cohort contained a longer history, from 2017 to the present.
For this data set members were evaluated at multiple points in time rather
instead of having a single prediction date. Each person was evaluated using
each month of eligibility as a prediction date. This approach was not used
with the CMS data because that data didn’t cover a long enough time period to
allow for multiple months. Members had to have 3 months of eligibility after
and had to be at least 18 years old on the prediction date. Table 2 below
summarizes the population selection.

PopulationSize | Selection Criteria
---|---
3,008,781 | Total members in the Healthfirst dataset
30,898,086 | Total member-months of eligibility
29,388,003 | 3 months of eligibility after the prediction date
19,470,511 | 18 years or older on the prediction date
Table 2: Cohort selection for the Healthfirst population

####  2.2.3 Combined Population

Each data set was split by person 80%/20% into train and test sets. In the
Healthfirst data set, doing the split by person ensures that all prediction
dates for a given individual were in the same set, eliminating the possibility
of data leakage. The data set for training was taken by taking all of the
positive examples from each data set along with a sampling of the negative
examples in the training set. The sampling was designed to meet two criteria.
First, the percentage of the adult population 65 or older was set to match
that of the US population (21% of adults). Second, the difference in
prevalence between those under and over 65 from the Healthfirst data set was
maintained in the overall data set. This difference was 3.9X.

The combined test population was created by unioning the full test set from
the CMS data and a random 20% sample of the Healthfirst test set. The full
test population contained 1,621,149 training examples with a prevalence of
3.86%. The test population had 761,898 examples with a prevalence of 0.36

### 2.3. Models

In order to encourage wide adoption and rapid implementation of the predictive
models, we created three separate models which represent different tradeoffs
between accuracy and ease of implementation. All three models were trained and
tested on the same data set. The models were built using the ClosedLoop
platform on the AWS cloud. The models are summarized in Table 3 below.

Model | Algorithm | Deployment | Features | SLA @ 5%
---|---|---|---|---
Survey | Logistic Regression | Online survey | 14 | 49.8%
Open Source | XGBoost | Open Source | 559 | 53.8%
Full | XGBoost | Hosted | 892 | 54.1%
Table 3: Different models created and their tradeoffs in terms of accuracy
versus ease of implementation. The full model has the highest accuracy while
the survey model has a very constrained feature set that is easy to implement.
Accuracy is measured using Sensitivity and Low Alert Rates (SLA).

The "Survey" model is the simplest and uses only a small number of features
that were designed to be able to be generated from a simple health risk
assessment questionnaire. This model requires no technical implementation, and
we have made it available through a web-based survey[15]. The "Open source"
model uses approximately 600 features derived from medical claims diagnosis
and utilization data. We have made this model available on GitHub[16].
Finally, the "Full" model was created that uses an extensive feature set
derived generated from medical claims data along with linked geographical and
social determinants of health data. This model is being made freely available
to healthcare organization. Information about accessing the platform can be
found at https://cv19index.com.

####  2.3.1 Survey Model

The first model is aimed at using a simple health history survey to enrich the
high-level recommendations from the CDC website[8] for identifying those
individuals who are at risk. The CDC identifies risk factors as:

  * •

Older adults

  * •

Individuals with heart disease

  * •

Individuals with diabetes

  * •

Individuals with lung disease

The purpose of the survey to to let an individual know their risk relative to
the general population with more detail than is available through the CDC. The
way this is achieved is by mapping questions related to an individuals medical
history into diagnosis code categories from the CCSR. We also included age and
gender as well as prior year hospital inpatient or emergency room (ER) visits.
In addition to the conditions coming from the recommendations of the CDC, we
included features that our other modeling efforts surfaced as important. The
mapping between the survey questions and CCSR codes is described in Table 4.
To turn this into a model, we extract ICD-10 diagnosis codes from the claims
in the year before the prediction date and aggregate them using the CCSR
categories. We create indicator features for the presence of any code in the
CCSR category. A logistic regression model is then trained on the available
claims data. A person’s percentile risk score is based on risk relative to the
other values in the training distribution. In addition to the CCSR codes,
Table 4 includes the coefficients associated with these features in the
logistic regression model.

Feature name | Coefficient | CCSR Code
---|---|---
Intercept | -6.74 | n/a
Age | 0.041 | n/a
Gender Male | 0.171 | n/a
Prior Admissions | 0.682 | n/a
Prior ER Visits | 0.413 | n/a
Chronic obstructive pulmonary disease (COPD) or | 1.167 | CCSR:RSP008,CCSR:END012
emphysema, cystic fibrosis, or chronic bronchitis |  |
Asthma | 1.393 | CCSR:RSP009
Obesity | 0.935 | CCSR:END009
Diabetes (other than when you were pregnant) | 0.096 | CCSR:END002,CCSR:END003,
|  | CCSR:END004,CCSR:END005
Hypertension (also called high blood pressure) | 0.832 | CCSR:CIR007,CCSR:CIR008
Congestive Heart Failure | 0.982 | CCSR:CIR019
Heart attack (also called myocardial infarction) | 0.159 | CCSR:CIR009,CCSR:CIR010
Rheumatic heart disease | 0.788 | CCSR:CIR001,CCSR:CIR002,CCSR:CIR011,
|  | CCSR:CIR014,CCSR:CIR015
Stroke | 0.285 | CCSR:CIR020,CCSR:CIR021
Sickle cell anemia/HIV infection/Transplant | 2.582 | CCSR:BLD005,CCSR:INF006,CCSR:FAC023
Chronic kidney disease | 0.966 | CCSR:GEN003
Hemodialysis | 1.369 | CCSR:GEN002
Liver disease | 0.055 | CCSR:DIG019
Pneumonia, acute bronchitis, influenza or other | 0.696 | CCSR:RSP002,CCSR:RSP003,
|  | CCSR:RSP005,CCSR:RSP006
acute respiratory infection |  |
Cancer | 1.091 | CCSR:NEO
Neurocognitive conditions | 0.294 | CCSR:NVS011,CCSR:CIR022,CCSR:CIR025
Pregnancy | 0.789 | CCSR:PRG
COPD x Age | -0.002 | n/a
Asthma x Age | -0.015 | n/a
Obesity x Age | -0.004 | n/a
Diabetes x Age | 0.000 | n/a
Hypertension x Age | 0.005 | n/a
Congestive heart failure x Age | -0.007 | n/a
Myocardial infarction x Age | 0.003 | n/a
Rheumatic heart disease x Age | -0.008 | n/a
Stroke x Age | -0.003 | n/a
Sickle cell/HIV/Transplate x Age | -0.028 | n/a
Chronic kidney disease x Age | -0.008 | n/a
Hemodialysis x Age | -0.018 | n/a
Liver disease x Age | 0.001 | n/a
Pneumonia, acute bronchitis, influenza or other | -0.005 | n/a
acute respiratory infection x Age |  |
Cancer x Age | -0.009 | n/a
Neurocognitive conditions x Age | 0.004 | n/a
Pregnancy, childbirth and the puerperium x Age | -0.003 | n/a
Table 4: Features used associated with risk factors identified by CDC
and their corresponding CCSR codes

####  2.3.2 Open Source Model

The Open Source model uses gradient boosted trees. Gradient boosted trees are
a machine learning method that use an ensemble of simple models to create
highly accurate predictions[10]. The resulting models demonstrate higher
accuracy. A drawback to these models is that they are significantly more
complex; consequently, "by hand" implementations of such models are
impractical. A nice feature of gradient boosted trees is that they are fairly
robust against learning features that are eccentricities of the training data,
but do not extend well to future data. As such, we allow full diagnosis
histories to be leveraged within our simpler XGBoost model. In this approach,
every category in the full CCSR is converted into an indicator feature,
resulting in 559 features. A 3-month delay was imposed on the claims data, so
that claims within the most recent 3 months before the prediction date were
not used to make the predictions. This 3-month delay simulates the delay in
claims processing that usually occurs in practical settings and enables the
model to be used with current claims data. This delay was not imposed in the
survey model since questionnaires would generally use current data. The Github
repository for the open source model contains scripts that automatically
prepare the features for this model from simple CSV files containing claims
data.

####  2.3.3 Full Model

We additionally built a model within the ClosedLoop platform. The ClosedLoop
platform is a software system designed to enable rapid creation of machine
learning models utilizing healthcare data. The full details of the platform
are outside the bounds of this paper; however, using the platform allows us to
leverage engineered features coming from peer-reviewed studies. Examples are
social determinants of health and the Charlson Comorbidity Index[13]. The
computation of these features from claims data is often complex or involves
the linking of additional data. These are operations handled by the CLosedLoop
platform that were not easy to extract into an open source format, but do
provide improved model accuracy.

## 3\. Results and Model Interpretation

We quantify the performance of the C-19 Index using metrics that are relevant
for it’s intended use. The standard receiver operating characteristic curves
show that all 3 models have identical areas under the curve, at 0.87. However,
if we consider the intended use of the model, which is to identify highly
vulnerable populations for additional targeting, a more appropriate metric is
sensitivity at low alerts rates. This metric measures the sensitivity of the
model when looking only at some small percentage of the population. This
sensitivity if plotted for alert rates up to 20% of the population in Figure
1. Additionally, the metrics quantifying the effectiveness of our models are
presented in Table 5. To provide a baseline for comparison, the models are
compared to the Charleson Comorbidity Index (CCI)[13]. The CCI is a typical
risk scoring algorithm based on claims data. In this case, we see that models
explicitly built for the proxy outcome are more accurate than the CCI.

Figure 1: Sensitivity of each model at alert rates up to 20%. Model | ROC AUC | Sensitivity at 5% Alert Rate
---|---|---
Survey | .86 | .498
Open Source | .87 | .538
Full | .87 | .541
Charlson Comorbidity | .79 | .291
Table 5: Measures of effectiveness for the three models

### 3.1. Validation

The model has subsequently been validated by evaluating its results against
approximately 14,000 hospital admissions for known COVID-19 cases in New York
City from 2/1/2020 until mid-May 2020. These admissions were We validated the
Vulnerability Index by comparing the mortality rate for these admissions with
their Vulnerability Index. All cases were insured by Healhfirst, and each
patient’s prior claims data was used to compute their Vulnerability Index.

Members who were admitted and survived had an average score of 2.4% and
members who were admitted and passed away had an average score of 3.3% (a 38%
relative difference). The ROC AUC was 0.68. This is lower than the test ROC,
but that drop is expected since this test set consisted only of patients who
had already been admitted for COVID-19, presumably removing many of the low
vulnerability patients and increasing the difficulty of the prediction
problem.

Decile | Medicaid Mortality % | Medicaid Lift | Medicare Mortality % | Medicare Lift
---|---|---|---|---
Top 10% | 6.48 | 1.30 | 5.56 | 2.99
Top 20% | 6.43 | 1.29 | 3.70 | 2.00
Top 30% | 6.39 | 1.29 | 3.09 | 1.66
Top 40% | 6.39 | 1.28 | 2.78 | 1.50
Top 50% | 6.21 | 1.25 | 2.59 | 1.40
Top 60% | 5.80 | 1.17 | 3.10 | 1.67
Top 70% | 5.51 | 1.11 | 2.65 | 1.43
Top 80% | 5.25 | 1.06 | 2.32 | 1.25
Top 90% | 5.11 | 1.03 | 2.06 | 1.11
Full | 4.97 | - | 1.86 | -
Table 6: Validation of the Vulnerability Index using COVID-19 admissions

The 14,000 cases were divided into a adult Medicaid population and a Medicare
population. Table 6 sorts each population by their vulnerability index and
shows the detah rate of each decile along with the lift relative to the
overall death rate.

### 3.2. Accessing Models

In an effort to make these models as broadly available as possible we have
provided several different avenues for the models to be used, each optimized
for a different user base. The logistic regression model powers a publicly
available web-based survey at http://c19survey.closedloop.ai. The open source
model is available through github at https://github.com/closedloop-
ai/cv19index. This model is written in the Python programming language. We
have included synthetic data for testing an wrapper code that converts tabular
medical claims data to the input format specific for our models. We encourage
the healthcare data science community to fork the repository and adapt it to
their own purposes. We encourage collaboration from the open-source community,
and pull requests will be considered for inclusion in the main branch of the
package. Finally, for those wishing to take advantage of the full model, we
are providing access to the COVID-19 model hosted on the ClosedLoop platform
free of charge. Please visit https://closedloop.ai/cv19index for instructions
on how to gain access.

### 3.3. Limitations

The approach taken in this paper has several limitations. Most notably, no
actual COVID-19 cases were used in the training of the model. The usefulness
of the model in predicting COVID-19 vulnerability is entirely dependent on the
actual occurrence of COVID-19 matching the proxy outcome. While the logic
behind these decisions is defensible, they are not currently backed up with
actual data. As COVID-19 case data becomes available, we expect to validate
the proxy outcome and determine if it is in fact appropriate. Eventually,
enough data will be available to build models on COVID-19 Vulnerability itself
without having to use a proxy.

Anther major limitation of these models is their reliance on claims data,
which does not contain enough clinical detail, such as lab values. For this
reason, we do not recommend using the C-19 Vulnerability Index in inpatient
settings, where more detailed data is likely available. The C-19 Index is most
useful in a population health context where the only data available is claims
data.

Finally, bassed on medical guidance the authors have decided exclude pediatric
populations from the training and test sets for the C-19 Index. At the time of
development, there was so little information available on COVID-19 that we
could not confidently assert that the proxy endpoint we were using was
appropriate for those under 18 years of age.

## 4\. Conclusions

This pandemic has already claimed tens of thousands of lives as of this
writing, and sadly this number is sure to grow. As healthcare resources are
constrained by the same scarcity constraints that affect us all, it is
important to empower intervention policy with the best information possible.
We have provided several implementations of the C-19 Index and means of access
for those individuals with varying levels of technical expertise. It is our
hope that by providing this tool quickly to the healthcare data science
community, widespread adoption will lead to more effective intervention
strategies and, ultimately, help to curtail the worst effects of this
pandemic.

## 5\. Addendum

Due to the early release of the first versions of the C19 Index, several
organizations were able to quickly apply the model to their populations. The
authors have been in contact with several organizations that have used the
vulnerability index to prioritize proactive outreach towards the most
vulnerable members. In many cases, these organizations had existing care
management teams that were able to rapidly shift to COVID-19 by applying a
different prioritization and focus to their interventions. In other cases, new
phone or text messaging campaigns were developed for COVID-19 that used the
C-19 Index as a prioritization mechanism. As the response to the pandemic
continues to develop, we will continue to update the C-19 Index and provide
more information on its usage.

## 6\. Ethics

All relevant ethical guidelines have been followed in this study. For the CMS
LDS data, we submitted a research proposal as part of our data use agreement
that discussed our plans to build predictive models using this data. For the
Healthfirst data, which did contain identifiable PHI, the data was covered
under a Business Associate Agreement between ClosedLoop.ai and Healthfirst and
the data was stored in ClosedLoop’s HIPAA comliant cloud storage. The project
was part of an ongoing Quality Improvement effort at Healthfirst and so did
not require preapproval by an Institutional Review Board. The study protocol
and this manuscript were reviewed and approved by Healthfirst’s compliance
team.

## 7\. Acknowledgements

We would like to thank Healthfirst for their collaboration on this work and
for allowing us to use insights from their dat ain generation of the model. We
would also like to thank Amazon Web Services for sponsoring this work with AWS
platform credits.

## References

  * [1] World Health Organization. “Coronavirus.” Who.Int, 2019, www.who.int/emergencies/diseases/novel-coronavirus-2019. Accessed 15 Mar. 2020.
  * [2] World Health Organization. “WHO Director-General’s Opening Remarks at the Media Briefing on COVID-19 - 11 March 2020.” Who.Int, 11 Mar. 2020, www.who.int/dg/speeches/detail/who-director-general-s-opening-remarks-at-the-media-briefing-on-covid-19—11-march-2020. Accessed 15 Mar. 2020.
  * [3] Specht, Liz. “Simple Math Offers Alarming Answers about Covid-19, Health Care.” STAT, 10 Mar. 2020, www.statnews.com/2020/03/10/simple-math-alarming-answers-covid-19/. Accessed 15 Mar. 2020.
  * [4] Page, Michael Le. “Why Is It so Hard to Calculate How Many People Will Die from Covid-19?” New Scientist, 11 Mar. 2020, www.newscientist.com/article/mg24532733-700-why-is-it-so-hard-to-calculate-how-many-people-will-die-from-covid-19/. Accessed 15 Mar. 2020.
  * [5] Zhou F, Yu T, Du R, et al. Clinical course and risk factors for mortality of adult inpatients with COVID-19 in Wuhan, China: a retrospective cohort study [published online ahead of print, 2020 Mar 11] [published correction appears in Lancet. 2020 Mar 12;:]. Lancet. 2020;S0140-6736(20)30566-3. doi:10.1016/S0140-6736(20)30566-3.
  * [6] Wu Z, McGoogan JM. Characteristics of and Important Lessons From the Coronavirus Disease 2019 (COVID-19) Outbreak in China: Summary of a Report of 72314 Cases From the Chinese Center for Disease Control and Prevention [published online ahead of print, 2020 Feb 24]. JAMA. 2020;10.1001/jama.2020.2648. doi:10.1001/jama.2020.2648
  * [7] Centers for Disease Control and Prevention. “Coronavirus Disease 2019 (COVID-19).” Centers for Disease Control and Prevention, 11 Feb. 2020, www.cdc.gov/coronavirus/2019-ncov/specific-groups/high-risk-complications.html.
  * [8] Centers for Disease Control and Prevention. “Coronavirus Disease 2019 (COVID-19).” Centers for Disease Control and Prevention, 11 Feb. 2020, www.cdc.gov/coronavirus/2019-ncov/specific-groups/high-risk-complications.html.
  * [9] Hubbard, Ruth E, et al. “Frailty Status at Admission to Hospital Predicts Multiple Adverse Outcomes.” Age and Ageing, vol. 46, no. 5, 2017, pp. 801–806, www.ncbi.nlm.nih.gov/pubmed/28531254, 10.1093/ageing/afx081. Accessed 11 Jan. 2020.
  * [10] James, Gareth, et al. An Introduction to Statistical Learning, with Applications in R. Springer, 2013.
  * [11] “Measures and Data Sources | County Health Rankings and Roadmaps.” Countyhealthrankings.Org, University of Wisconsin Population Health Institute, 2019, www.countyhealthrankings.org/explore-health-rankings/measures-data-sources. Accessed 1 Feb. 2020.
  * [12] The Healthcare Cost and Utilization Project. “Clinical Classifications Software Refined (CCSR) for ICD-10-CM Diagnoses.” www.Hcup-Us.Ahrq.Gov, www.hcup-us.ahrq.gov/toolssoftware/ccsr/ccs refined.jsp. Accessed 15 Mar. 2020.
  * [13] Charlson, M E, et al. “A New Method of Classifying Prognostic Comorbidity in Longitudinal Studies: Development and Validation.” Journal of Chronic Diseases, vol. 40, no. 5, 1987, pp. 373–83, 10.1016/0021-9681(87)90171-8.
  * [14] Centers for Medicare & Medicaid Services. “Limited Data Set (LDS) Files.” www.cms.gov, https://www.cms.gov/Research-Statistics-Data-and-Systems/Files-for-Order/LimitedDataSets. Accessed 14 Apr. 2020.
  * [15] ClosedLoop.ai. “COVID-19 Vulnerability Index” closedloop.ai, https://closedloop.ai/cv19index/. Accessed 14 Apr. 2020.
  * [16] ClosedLoop.ai. “COVID-19 Vulnerability Index GitHub” github.com, https://github.com/closedloop-ai/cv19index. Accessed 14 Apr. 2020.
  * [17] Sanche, Steven, et al. “High Contagiousness and Rapid Spread of Severe Acute Respiratory Syndrome Coronavirus 2.” Emerging Infectious Diseases, vol. 26, no. 7, July 2020, wwwnc.cdc.gov/eid/article/26/7/20-0282_article, 10.3201/eid2607.200282. Accessed 21 Apr. 2020.
  * [18] Markian Hawryluk. “Mysterious Heart Damage, Not Just Lung Troubles, Befalling COVID-19 Patients.” Kaiser Health News, 6 Apr. 2020, khn.org/news/mysterious-heart-damage-not-just-lung-troubles-befalling-covid-19-patients/. Accessed 21 Apr. 2020.
  * [19] Chow, Nancy, et al. “Preliminary Estimates of the Prevalence of Selected Underlying Health Conditions Among Patients with Coronavirus Disease 2019 — United States, February 12–March 28, 2020.” MMWR. Morbidity and Mortality Weekly Report, vol. 69, no. 13, 3 Apr. 2020, pp. 382–386, www.cdc.gov/mmwr/volumes/69/wr/mm6913e2.htm?s_cid=mm6913e2_w, 10.15585/mmwr.mm6913e2. Accessed 21 Apr. 2020.

## 8\. Appendix: Full Feature List

We include the list of features used in the Full model. This was a feature set
that could be implemented across the 2 data sets and utilizes just medical
claims data. The ClosedLoop platform is capable of incorporating a wide
variety of other health data, including pharmacy claims, electronic health
records, patient reported data, and linked census data. That additional data
was not used in this model due ot the constraints around our usage of the CMS
data. The majority of features are binary variables indicating if a patient
has had one type of medical event 15 moths prior to the date of prediction,
excluding the 3 most recent months.

Gender Age Admission Count ED Visit Count Office Visit Count Medical
Eligibility Months # Distinct DME Categories (12M) # Days Since Last Hospital
Discharge Use of Preventative Services Procedure History DME # ICU Stays (12M)
Frailty Indicator Discharge Disposition Diagnosis History Monthly Medical Cost
# of Admissions (12M) Inpatient Days # of Observation Stays # of ER Visits
(12M) # of ER Visits (6M) Admission Cause Charlson Comorbidity Index (CCI)
Evidence of tobacco smoking Evidence of smokeless tobacco Evidence of tobacco
non-use Falls risk assessment performed Abnormal BMI Normal BMI LDL LVEF
Systolic Blood Pressure Diastolic Blood Pressure # Distinct CCSR Body Systems
# Distinct CCSR Diagnosis Categories (12M) PCS Procedure History Prior
Respiratory Infections Prior Hospital Acquired Infections Diagnosis History
High-Level
---
Table 7: Full feature list available within platform Evidence of smokeless
tobacco
---
Evidence of tobacco non-use
Falls risk assessment performed
Abnormal BMI
Normal BMI
LDL
LVEF
Systolic Blood Pressure
Diastolic Blood Pressure
Continuity of Care Index (12M)
PQI Any
PQI Diabetes
PQI 5
PQI 7
PQI 8
PQI 11
PQI 12
# Distinct CCSR Body Systems
# Distinct CCSR Diagnosis Categories (12M)
PCS Procedure History
Prior Respiratory Infections
Prior Hospital Acquired Infections
Diagnosis History High-Level
Table 8: Full feature list available within platform (continued)


# Optimizing Codes for Source Separation in Color Image Demosaicing and
Compressive Video Recovery

Alankar Kotwal<EMAIL_ADDRESS>Department of Electrical
Engineering, IIT Bombay  Ajit Rajwade<EMAIL_ADDRESS>Department of
Computer Science and Engineering, IIT Bombay

###### Abstract

There exist several applications in image processing (eg: video compressed
sensing [1] and color image demosaicing [2]) which require separation of
constituent images given measurements in the form of a coded superposition of
those images. Physically practical code patterns in these applications are
non-negative, systematically structured, and do not always obey the nice
incoherence properties of other patterns such as Gaussian codes, which can
adversely affect reconstruction performance. The contribution of this paper is
to design code patterns for video compressed sensing and demosaicing by
minimizing the mutual coherence of the matrix $\boldsymbol{\Phi\Psi}$ where
$\boldsymbol{\Phi}$ represents the sensing matrix created from the code, and
$\boldsymbol{\Psi}$ is the signal representation matrix. Our main contribution
is that we explicitly take into account the special structure of those code
patterns as required by these applications: (1) non-negativity, (2) block-
diagonal nature, and (3) circular shifting. In particular, the last property
enables for accurate and seamless patch-wise reconstruction for some important
compressed sensing architectures.

###### keywords:

video compressed sensing, color image demosaicing, source separation, sensing
matrices, mutual coherence, compressed sensing

††journal: Elsevier Journal of Signal
Processingmymainaddressmymainaddressfootnotetext: AK gratefully acknowledges
support from the Electrical Engineering department at IIT
Bombay.mymainaddress2mymainaddress2footnotetext: AR gratefully acknowledges
support from IIT Bombay Seed Grant number 14IRCCSG012.

##  1 Introduction and Related Work

C

ompressed sensing has been explored as an alternative way of sampling
continuous-time signals. Its success with still images has inspired efforts to
apply it to video. Indeed, [1, 3] achieve compression across time by combining
frames into coded snapshots while sensing and separating them with a pre-
trained over-complete dictionary. [1, 3] make a diagonal choice for the
sensing matrix for each frame. $T$ input frames of size $N_{1}\times N_{2}$
denoted as $\\{\boldsymbol{X_{i}}\\}_{i=1}^{T}$ are sensed so that output
$\boldsymbol{Y}$ of size $N_{1}\times N_{2}$ appears as a pixel-wise coded
superposition (dictated by the sensing matrices $\boldsymbol{\Phi_{i}}$) of
the input frames. The pixel-wise codes represent modulation of incident light
by patterns as per a pre-programmed DMD (digital micromirror device) or LCoS
(liquid crystal on silicon) device. $\boldsymbol{Y}$ is sometimes termed the
coded snapshot image. The sensing framework is

$\textrm{vec}(\boldsymbol{Y})=\sum_{i=1}^{T}\boldsymbol{\Phi_{i}}\textrm{vec}(\boldsymbol{X_{i}})=\boldsymbol{\Phi}\textrm{vec}(\boldsymbol{X}).$
(1)

Here, each $\boldsymbol{\Phi_{i}}$ is a $n\times n$ (where $n=N_{1}N_{2}$ is
the number of pixels per frame) non-negative (possibly, but not necessarily
binary) diagonal sensing matrix with vectorized code elements on the diagonal,
and the overall sensing matrix
$\boldsymbol{\Phi}=(\boldsymbol{\Phi_{1}}|\boldsymbol{\Phi_{2}}|...\boldsymbol{\Phi_{T}})$
has size $n\times nT$. The complete $N_{1}\times N_{2}\times T$ video is
represented as
$\boldsymbol{X}=(\boldsymbol{X_{1}}|\boldsymbol{X_{2}}|...|\boldsymbol{X_{T}})$.
This acquisition model also holds in color image demosaicing with $T=3$.

The sparsifying basis here could be a 3D dictionary learned on video patches
or a universal basis such as 3D-DCT (discrete cosine transform). Given this
dictionary, called $\boldsymbol{D}$, any given signal $\boldsymbol{X}$, and in
particular, its frames $\\{\boldsymbol{X_{i}}\\}_{i=1}^{T}$ can be
approximately reconstructed as a sum of its projections
$\boldsymbol{\alpha_{j}}$ on the $K$ atoms in $\boldsymbol{D}$:

$\textrm{vec}(\boldsymbol{X_{i}})=\sum_{j=1}^{K}\boldsymbol{D_{ji}\alpha_{j}}$
(2)

where $\boldsymbol{D_{ji}}$ is the $i^{\text{th}}$ frame in the
$j^{\text{th}}$ 3-D dictionary atom $\boldsymbol{D_{ji}}$. From the
measurements and the dictionary, the input images are recovered solving the
following optimization problem:

$\min_{\boldsymbol{\alpha}}\|\boldsymbol{\alpha}\|_{0}\text{ subject to
}\left\|\textrm{vec}(\boldsymbol{Y})-\sum_{i=1}^{T}\boldsymbol{\Phi_{i}}\sum_{j=1}^{K}\boldsymbol{D_{ji}\alpha_{j}}\right\|_{2}\leq\epsilon.$
(3)

This problem can be approximately solved with sparse recovery techniques like
orthogonal matching pursuit [4]. The drawback here, though, is that the 3D
dictionary imposes a smoothness assumption on the scene. Since a linear
combination of dictionary atoms cannot speed an atom up, the typical speeds of
objects moving in the video must be roughly the same as the dictionary. Also,
because of the nature of the training data, the dictionary fails to sparsely
represent sudden scene changes caused by, say, lighting or occlusion. Other
techniques like [5] exploit additional structure within the signal, like
periodicity, rigid motion or analytical motion models and cannot be used in
the general video sensing case.

We try relaxing these constraints using a source-separation approach [6],
where precise error bounds on the recovery of the images have been derived,
with possible improvement using the techniques in [7]. Each of the coded
snapshots is treated as a mixture of sources, each sparse in some basis. We
experimented with basis pursuit recovery with Gaussian-random sensing
matrices, getting excellent results with no visible ghosting for both similar
and radically different images. Unfortunately, the actually realizable
positive sensing matrices are more structured, and do not have the nice
incoherence properties of Gaussian-random matrices, which are sufficient
conditions for near-accurate recovery as derived in [6].

Survey of Previous Work: Here, we aim to design such _structured_ sensing
matrices with low mutual coherence (see Section 2.2 for a precise definition),
making them ideal for compressed video. The existing literature on designing
compressed sensing matrices can be broadly divided into two categories: (1)
optimizing mutual coherence and (2) information theoretic approaches.
Techniques in the former category include [8, 9, 10, 11, 12, 13, 14, 15]. Most
approaches [10, 13, 12, 11, 14, 15] involve minimizing the difference between
a Gram matrix and the identity matrix, thereby minimizing some forms of an
average of normalized dot products of effective dictionary columns. Although
minimizing averages doesn’t guarantee minimizing the maximum (which is the
mutual coherence in this case) of the quantities forming this average, the
techniques have been quite effective. Methods such as [8] first estimate
signal support using matrices designed to minimize coherence, followed by
matrices designed to reduce MSE on the known support. [9] specializes to
orthogonal or biorthogonal bases. However these methods are designed for
general compressed sensing matrices and do not account for the _special
structure_ of the sensing matrices used for video compressed sensing or for
demosaicing, a framework which this paper expressly deals with. Techniques in
the second category include [16, 17, 18]. These papers design sensing matrices
$\boldsymbol{\Phi}$ such that the mutual information between a set of small
patches $\\{\boldsymbol{x_{i}}\\}_{i=1}^{N_{p}}$ and their corresponding
projections $\\{\boldsymbol{y_{i}}\\}_{i=1}^{N_{p}}$ where
$\boldsymbol{y_{i}}=\boldsymbol{\Phi x_{i}}$, is maximized. Computing this
mutual information first requires estimation of the probability density
function of the patches and their projections, using Gaussian mixture models
for instance. This can be expensive and is an iterative process. Moreover
these learned GMMs for a class of patches may not be general enough.

Contributions: Quite importantly, the aforementioned techniques in both
categories do not account for the special structure of the sensing matrices
described for video compressed sensing [1, 3] and demosaicing [19, 2], which
is the primary aim of this paper. We design codes to enable the reconstruction
of small patches from $\boldsymbol{Y}$ in the model described in Eq. 1 using
algorithms such as basis pursuit. This is taking into account the fact that
although the acquisition of $\boldsymbol{Y}$ is performed at the level of a
group of $T$ consecutive frames (and specifically, _not_ at the level of
individual patches), the reconstruction can in fact proceed patch-wise as
shown in [1], and as we shall show in Section 2.5. Hence the sensing matrix
can in fact be designed at the level of small patches, followed by a _tiling_
to create the complete sensing matrix. We also demonstrate the benefit of
designing patch-based codes whose _circular shifts_ also have low mutual
coherence. Reconstruction results are presented for video compressive recovery
and color image demosaicing.

Organization: The complete theoretical development is presented in Section 2,
followed by experimental results on video compressed sensing and color image
demosaicing in Section 3, followed by a discussion in Section 4.

##  2 Method and Analysis

###  2.1 Our framework

We propose to use a recovery method different from the one used in [1], within
the same acquisition framework. Thus, our signals are still acquired according
to Eq. 1. However, the choice of the sparsifying basis is different: we use a
2D-DCT basis $\boldsymbol{D}$ to model each frame in the input data, though
our method is general enough to work for any other basis. The dictionary
$\boldsymbol{\Psi}$ sparsifying the entire video sequence, thus, is a block-
diagonal matrix with the $n\times n$ sparsifying basis $\boldsymbol{D}$ on the
diagonal where $n=N_{1}N_{2}$ is the number of pixels per video frame. Thus,

$\displaystyle\textrm{vec}(\boldsymbol{Y})$
$\displaystyle=\begin{pmatrix}\boldsymbol{\Phi_{1}}&\ldots&\boldsymbol{\Phi_{T}}\end{pmatrix}\begin{pmatrix}\boldsymbol{D\alpha_{1}}&\ldots&\boldsymbol{D\alpha_{T}}\end{pmatrix}^{T}$
(4)
$\displaystyle=\begin{pmatrix}\boldsymbol{\Phi_{1}D}&\ldots&\boldsymbol{\Phi_{T}D}\end{pmatrix}\begin{pmatrix}\boldsymbol{\alpha_{1}}&\ldots&\boldsymbol{\alpha_{T}}\end{pmatrix}^{T}$
(5)

Given a measurement $\boldsymbol{Y}$, we recover the input
$\\{\boldsymbol{X_{i}}\\}_{i=1}^{T}$ through the DCT coefficients
$\boldsymbol{\alpha}$ by solving the optimization problem

$\min_{\boldsymbol{\alpha}}\|\boldsymbol{\alpha}\|_{1}\text{ subject to
}\|\textrm{vec}(\boldsymbol{Y})-\boldsymbol{\Phi\Psi\alpha}\|_{2}\leq\epsilon,\
\boldsymbol{\alpha}=\begin{pmatrix}\boldsymbol{\alpha_{1}}&\boldsymbol{\alpha_{2}}&\ldots&\boldsymbol{\alpha_{T}}\end{pmatrix}^{T}.$
(6)

In our implementation we used the CVX [20] solver for solving the convex
optimization problem in Eq. 6.

###  2.2 Calculation of mutual coherence for our matrices

Our aim, then, is to optimize the sensing matrices $\boldsymbol{\Phi_{i}}$
directly for minimum mutual coherence with gradient descent. The criterion of
coherence is motivated by the fact that the error bounds on signal recovery
presented in [6] are monotonically increasing functions of coherence. We now
calculate gradients of the mutual coherence with respect to the (non-zero,
i.e., diagonal) elements of $\boldsymbol{\Phi_{i}}$. As in Eq. 5, with an
$n\times n$ dictionary $\boldsymbol{D}$ for representing a single 2D frame, we
have the effective dictionary

$\displaystyle\boldsymbol{\Phi\Psi}=\begin{pmatrix}\boldsymbol{\Phi_{1}D}&\boldsymbol{\Phi_{2}D}&\ldots&\boldsymbol{\Phi_{T}D}\end{pmatrix}.$
(7)

The mutual coherence of a general matrix $\boldsymbol{H}$, with the
$i^{\text{th}}$ column defined as $\boldsymbol{h_{i}}$, is

$\displaystyle\mu(\boldsymbol{H})=\max_{i\neq
j}\frac{|\left<\boldsymbol{h_{i}},\boldsymbol{h_{j}}\right>|}{\sqrt{\left<\boldsymbol{h_{i}},\boldsymbol{h_{i}}\right>\left<\boldsymbol{h_{j}},\boldsymbol{h_{j}}\right>}}.$
(8)

This expression contains max and abs functions that a gradient-based scheme
cannot handle. Instead, we soften the max and convert the abs to a square by
using the following approximation, which is very accurate for large enough
$\theta$,

$\displaystyle\max_{i}\\{t_{i}^{2}\\}_{i=1}^{n}\approx\frac{1}{\theta}\log\sum_{i=1}^{n}e^{\theta
t_{i}^{2}}.$ (9)

We need to evaluate the mutual coherence of the dictionary
$\boldsymbol{\Phi\Psi}$ as a function of the non-zero elements of
$\boldsymbol{\Phi}$. We will call the time index varying from $1$ to $T$ as
$\mu$ or $\nu$, and the spatial index varying from $1$ to $n$ as $\alpha$,
$\beta$ or $\gamma$. The $\mu^{\text{th}}$ block of $\boldsymbol{\Phi}$ is
thus $\boldsymbol{\Phi_{\mu}}$. Let the $\beta^{\text{th}}$ diagonal element
of $\boldsymbol{\Phi_{\mu}}$ be $\phi_{\mu\beta}$. Define the
$\alpha^{\text{th}}$ column of $\boldsymbol{D}^{T}$ to be
$\boldsymbol{d_{\alpha}}$. Then, it can be shown [A] that the normalized dot
product between the $\beta^{\text{th}}$ column of the $\mu^{\text{th}}$ block
and the $\gamma^{\text{th}}$ column of the $\nu^{\text{th}}$ block is

$\displaystyle M_{\mu\nu}(\beta\gamma)$
$\displaystyle=\frac{\sum_{\alpha=1}^{n}\phi_{\mu\alpha}\phi_{\nu\alpha}d_{\alpha}(\beta)d_{\alpha}(\gamma)}{\sqrt{\left(\sum_{\alpha=1}^{n}\phi_{\mu\alpha}^{2}d^{2}_{\alpha}(\beta)\right)\left(\sum_{\tau=1}^{n}\phi_{\nu\tau}^{2}d^{2}_{\tau}(\gamma)\right)}}.$
(10)

Finally, using the squared soft-max function [Eq. 9] to deal with the max and
the abs in the mutual coherence expression, we get the squared soft mutual
coherence $\mathcal{C}$ to be

$\displaystyle\mathcal{C}=\frac{1}{\theta}\log\left[\sum_{\mu=1}^{T}\sum_{\nu=1}^{\mu-1}\sum_{\beta=1}^{n}\sum_{\gamma=1}^{n}e^{\theta
M_{\mu\nu}^{2}(\beta\gamma)}+\sum_{\mu=1}^{T}\sum_{\beta=1}^{n}\sum_{\gamma=1}^{\beta-1}e^{\theta
M_{\mu\mu}^{2}(\beta\gamma)}\right].$ (11)

In the above, the first term corresponds to all ($\mu>\nu$) blocks that are
below the block diagonal. Here, we consider all terms in the given block for
the maximum. The second term corresponds to ($\mu=\nu$) blocks on the block
diagonal. Here, we consider only consider ($\beta>\gamma$) below-diagonal
elements for the maximum.

Note that lowering the value of $\theta$ in Eq 9 will take us closer toward an
average of the normalized dot-product values (often referred to as ‘average
coherence’) instead of the mutual coherence defined in Eq. 9 and Eq. 8. While
this approach has been adopted in many papers, we observed equally good
results with a large $\theta$ and with a much faster optimization. As such,
our method remains unchanged even with smaller values of $\theta$.

###  2.3 Calculation of coherence derivatives

We note that the $\mathcal{C}$ computed in the section above is a function of
$\boldsymbol{\Phi}$. We differentiate $\mathcal{C}$ with respect to
$\Phi_{\delta\epsilon}$, the elements located on the diagonals of each of the
square components, i.e. on the diagonals of
$\boldsymbol{\Phi_{1}},\boldsymbol{\Phi_{2}},...,\boldsymbol{\Phi_{T}}$. These
are in some sense, the $nT$ ‘degrees of freedom’ of the $n\times nT$ matrix
$\boldsymbol{\Phi}$. For this, we define the numerator of the expression for
$M_{\mu\nu}(\beta\gamma)$ as $\chi_{\mu\nu}(\beta\gamma)$ and the denominator
as $\xi_{\mu\nu}(\beta\gamma)$. The derivative of the objective function can
be found in terms of these quantities. Defining $\uparrow_{\mu\delta}$ to be
the Kronecker delta function that is 1 if and only if $\mu=\delta$, it can be
shown that [B]

$\frac{d\chi_{\mu\nu}(\beta\gamma)}{d\phi_{\delta\epsilon}}=d_{\epsilon}(\beta)d_{\epsilon}(\gamma)\left(\phi_{\mu\epsilon}\uparrow_{\nu\delta}+\uparrow_{\mu\delta}\phi_{\nu\epsilon}\right)$
(12)
$\frac{d\xi_{\mu\nu}(\beta\gamma)}{d\phi_{\delta\epsilon}}=\frac{1}{\xi_{\mu\nu}(\beta\gamma)}\left[\phi_{\mu\epsilon}d_{\epsilon}^{2}(\beta)\uparrow_{\mu\delta}\sum_{\tau=1}^{n}\phi_{\nu\tau}^{2}d^{2}_{\tau}(\gamma)+\phi_{\nu\epsilon}d_{\epsilon}^{2}(\gamma)\uparrow_{\nu\delta}\sum_{\alpha=1}^{n}\phi_{\mu\alpha}^{2}d^{2}_{\alpha}(\beta)\right].$
(13)

Using these, we perform projected gradient descent (to maintain non-negativity
of $\boldsymbol{\Phi}$) with adaptive step-size and use a multi-start strategy
to combat the non-convexity of the problem. We settle for the matrix produced
by the start that yielded the least value of $\mathcal{C}$. In practice, we
needed a smaller number ($\sim 20$) of starts.

For reasons to become clear further in this paper, we term matrices designed
this way as ‘non-circularly designed matrices’.

###  2.4 Time complexity and the need for something more

The calculation of coherence for a matrix requires us to evaluate normalized
dot products between columns of the matrix. In our case, the size of the
matrix is $n\times nT$, and each dot product needs $\mathcal{O}(n)$
operations, warranting the calculation of $\mathcal{O}(n^{3}T^{2})$
quantities. Optimizing this rapidly becomes intractable as $n$ increases. The
performance of gradient descent on this non-convex optimization problem also
worsens as the dimensionality of the search-space ($\mathcal{O}(nT)$)
increases. Empirically, we observe that it is intractable to design codes that
are more than $20\times 20$ in size in any reasonable time, whereas the size
of typical images in video compressed sensing or demosaicing problems is far
higher. This points to the fact that we need something more to make designing
effective codes possible.

###  2.5 Circularly-symmetric coherence minimization

An important aspect of the compressive architecture described in Eq. 1 is as
follows. On one hand, the compressive measurement corresponds to a sequence of
_entire_ video frames producing a coded snapshot image $\boldsymbol{Y}$ of
size $N_{1}\times N_{2}$. On the other hand however, the reconstruction can
actually be performed in a patch-wise manner across overlapping patches in a
sliding window fashion, as described in [1] (or as in [21], for a related
compressive architecture in hyperspectral imaging). This is because for a
$m\times m$ patch $\boldsymbol{y_{i}}$ ($m^{2}\ll n$) at location $i$ in
$\boldsymbol{Y}$, we have

$\textrm{vec}(\boldsymbol{y_{i}})=\sum_{t=1}^{T}\boldsymbol{\phi_{it}}\textrm{vec}(\boldsymbol{x_{it}})=\sum_{t=1}^{T}\boldsymbol{\phi_{it}\tilde{D}_{m^{2}\times
K}\alpha_{it}},$ (14)

where $\boldsymbol{\phi_{it}}$ is a $m^{2}\times m^{2}$ diagonal sub-matrix of
$\boldsymbol{\Phi_{t}}$ representing the portion of the code that modulates
the pixels inside an $m\times m$ patch around location $i$ in frame $t$, and
$\boldsymbol{x_{it}}$ represents a $m\times m$ patch around location $i$ in
frame $\boldsymbol{X_{t}}$. Exploiting the fact that each patch
$\boldsymbol{x_{it}}$ will be sparse in some dictionary
$\boldsymbol{\tilde{D}}_{m^{2}\times K}$, its dictionary codes
$\boldsymbol{\alpha_{it}}$ (and hence the patch itself) can be independently
inferred for each $i$, by solving the following optimization problem:

$\textrm{min}\|\boldsymbol{\alpha_{i}}\|_{1}\textrm{ subject to
}\|\textrm{vec}(\boldsymbol{y_{i}})-\boldsymbol{\phi_{i}\tilde{D}\alpha_{i}}\|_{2}\leq\epsilon$
(15)

where
$\boldsymbol{\phi_{i}}=(\boldsymbol{\phi_{i1}}|\boldsymbol{\phi_{i2}}|...|\boldsymbol{\phi_{iT}})$
has size $m^{2}\times
m^{2}T$,$\boldsymbol{\alpha_{i}}=(\boldsymbol{\alpha_{i1}}|\boldsymbol{\alpha_{i2}}|...|\boldsymbol{\alpha_{iT}})$
has size $TK\times 1$, and $\boldsymbol{\tilde{D}}\in\mathbb{R}^{m^{2}\times
KT}$ is a dictionary with $KT$ columns obtained by concatenating identical
dictionaries $\boldsymbol{\tilde{D}}_{m^{2}\times K}$ with $K$ columns. The
reconstruction of $\boldsymbol{x_{i}}=\boldsymbol{\tilde{D}\alpha_{i}}$ is
then followed by an averaging operation in pixels that are covered by more
than one patch. This avoids patch-seam artifacts that would occur if only non-
overlapping patches were reconstructed.

As a matter of fact, this leads us to designing smaller code-masks (of size
$m\times m$ per frame, where $m^{2}\ll n$) and tiling them to fit the image
size we are dealing with. This helps us get around the issue of computational
tractability. A small coherence for the designed patch guarantees good
reconstruction for patches of size $m\times m$ exactly aligned with the code-
mask; however, other patches see a code that is a circular shift of the
original code. Fig 1 provides a visual explanation. The big outer square
denotes the image. On top of the image we show tiled designed codes. Now, the
patch in red clearly multiplies with the exact designed code; however the
patch in green multiplies with a code that is a version of the original code
but shifted in both the coordinates circularly. Of course, one could seek to
reconstruct only those patches that are aligned exactly with the code mask,
but that will produce patch-seam artifacts and require heuristic post-
processing operations like deblocking.

Figure 1: Motivation behind circularly-shifted optimization. The code-mask for
the green patch is not aligned with the usual patch boundaries unlike the red
patch, and is in fact a circular shift of the original patch-aligned code-mask
(see text for more details).

This points to designing sensing matrices that have small coherence in all
their circular permutations. Note that these permutations happen in two
dimensions and must be handled as such. To this end, we modify the above
objective function to minimize the maximum coherence resulting from all
circularly-shifted vectorized versions of $\boldsymbol{\Phi}$. We thus have

$\displaystyle\mathcal{C}=\frac{1}{\theta}\log\left[\sum_{\zeta\in\text{perm}(\Phi)}\left[\sum_{\mu=1}^{T}\sum_{\nu=1}^{\mu-1}\sum_{\beta=1}^{n}\sum_{\gamma=1}^{n}e^{\theta
M_{\mu\nu}^{(\zeta)2}(\beta\gamma)}+\sum_{\mu=1}^{T}\sum_{\beta=1}^{n}\sum_{\gamma=1}^{\beta-1}e^{\theta
M_{\mu\mu}^{(\zeta)2}(\beta\gamma)}\right]\right]$ (16)

where $M_{\mu\nu}^{(\zeta)}(\beta\gamma)$ represents the normalized dot
product between the $\beta^{\text{th}}$ column of the $\mu^{\text{th}}$ block
and the $\gamma^{\text{th}}$ column of the $\nu^{\text{th}}$ block, resulting
from the instance of the circular permutation $\zeta$ of $\Phi$. Derivatives
of this expression are found exactly like in B, except that the $\mu$, $\nu$,
$\beta$ and $\gamma$ parameters are now subjected to the appropriate circular
permutation. Matrices designed this way are henceforth referred to as
‘circularly designed matrices’

For $m\times m$ patches, the time complexity for determining this maximum
coherence among all circular permutations is $\mathcal{O}(\tilde{m}^{4}T^{2})$
where $\tilde{m}=m\times m$, out of which a $\mathcal{O}(\tilde{m}^{3}T^{2})$
term arises from the calculation of coherence for each circular permutation,
and an extra $\mathcal{O}(\tilde{m})$ factor arises from the fact that there
are $m^{2}$ such permutations. The advantage here, though, is that we don’t
need to optimize masks having a very large size; we can do away with keeping
$m$ a small constant because the scheme works for any $m$ such that $m\times
m$ patches are sparse in some chosen dictionary. This scheme is, thus, more
scalable in terms of the size of the input image. Therefore the effective
dimension of the optimization problem in such a scheme is, in terms of the
variables that matter, $\mathcal{O}(T^{2})$.

It is worth mentioning that this simple idea has been largely ignored in
literature concerning sensing matrix optimization. As mentioned in the
introduction, previous attempts mostly use an average coherence minimization
technique [11, 12, 13] for sensing matrices, and ignore not only the special
structure of the matrices we are concerned with here but also the important
issue of patch overlap. Moreover, none of these techniques are directly
scalable for large images because they involve optimization problems in
variables whose dimensions are at least of the order of image size. Sensing
matrices can be designed at the patch level as well, for instance using
information theoretic techniques as in [16, 17, 18], but the methods therein
are not designed to account for either the special strcuture nor the issue of
overlapping reconstruction. To the best of our knowledge, ours is the first
piece of work to handle this important issue in a principled manner.

##  3 Validation and Results

In this section, we show extensive results on reconstruction in the context of
video compressive sensing and color image demosaicing. However, we start with
testing the proposed framework visually. In each figure caption are values of
the relative root mean square errors between the images and their
reconstructed versions.

###  3.1 Demonstrating a need for optimization

We first show a need for optimizing codes by demonstrating source separation
results on frames from a $301\times 225$ 30 fps video scene. In this
subsection, all reconstructions are performed in a non-overlapping way using
Eq. 15. With positive [0, 1] uniform random codes, the relative root mean
square errors are around 8% for each image. The results are shown in Fig. 2.

Figure 2: Real images, reconstructions with [0,1] uniform-random matrices. Up:
input images, down: reconstructions. RRMSEs = 0.081 and 0.084 for the two
reconstructions respectively. Compare with Fig. 3

Looking at these results, one notices that there is very little to no
ghosting, that is, appearance of features from one image into the other, in
the output images even when the images are very close to each other. This is a
very desirable property in any algorithm that separates images from compressed
video. These results are as good as the ones that Gaussian random codes
generate.

Next, we try separating three images from the same video with uniform
matrices. The third frame contains some significant difference in content as
compared to the earlier frames, and hence the reconstruction from data
involving random codes reveals ghosting artifacts. For instance, where there
is a car in the first two images, there is plain road in the third. See Fig.
3. Here, we notice ghosting artifacts appearing in the third frame (see the
box in the third frame reconstruction). However, with better designed sensing
matrices, one can think of getting rid of this effect, providing a motivation
for matrix design in coded source separation. The relative root mean square
errors here are worse, around 10% for each image.

Figure 3: Separating three images, uniform matrices. Up: input images, down:
reconstructions. RRMSEs = 0.0814, 0.0926 and 0.0859 for the three
reconstructions respectively. Compare with Fig. 2. The ghost of the car from
the first two frames appears in the white box in the reconstruction for the
third frame.

We perform a numerical comparison between our designed codes and random codes
for various values of $s=\|\boldsymbol{X_{i}}\|_{0}/n$ and $T$. We randomly
generate $T$ $s$-sparse (in 2D DCT) $8\times 8$ signals
$\\{\boldsymbol{X_{i}}\\}_{i=1}^{T}$, combine them using either random or
designed codes to get $\boldsymbol{Y}$. Average relative root mean square
errors on recovering the input signals from $\boldsymbol{Y}$ as a function of
$s$ and $T$ are shown in Figs. 11 and 12. Errors are near-zero in the region
where both $T$ and $s$ are small, and one can expect reasonable quality
reconstructions till $T=4$ from random matrices. To increase $T$ further, we
would need to optimize our sensing matrix appropriately, as is shown further
in this paper.

###  3.2 Coherence minimization

The coherence of a uniform random matrix of the type we’re interested in has a
typical value around 0.8 for $8\times 8$ codes. The distribution of these
values is shown in the boxplot in Fig 4.

Figure 4: Distribution of coherences for $8\times 8$ random positive codes as
a function of $T$

The typical profile of descent on coherence from a random initialization is
shown in Fig. 5.

Figure 5: Typical coherence decrease profile

The minimum coherence we have been able to achieve in this scheme has been
around 0.27 (for $T=2$). It is interesting to note that all initialization
instances lead to coherences (for $T=2$) of at the most 0.35, and hence
empirically yield nearly as good matrices.

###  3.3 Circularly-symmetric coherence minimization

Again, we design $8\times 8$ codes for $T=2$. To show coherence improvement
between positive random codes, and codes designed with and without circular
permutations, we plot the distribution of coherences of $\Phi^{(\zeta)}D$ in
Fig. 6 for all circular permutations $\zeta$. Note that even though the
coherences of non-circularly designed matrices are much lower than positive
random matrices, the maximum coherence among all permutations is quite large.
The circularly-designed matrices, however, have permuted coherences clustered
around a low value. We then expect good reconstruction with all circular
permutations, yielding good expected reconstructions for images.

Figure 6: Left to right: Circularly-shifted coherence histograms, in order,
for random matrices, non-circularly optimized matrices and circularly
optimized matrices

###  3.4 Demosaicing

To demonstrate the utility of this scheme, we show results on demosaicing RGB
images. The general demosaicing problem involves addressing the difficulty
that on a camera sensor, a single pixel can sense only one of the three R, G
and B channels. Therefore, raw camera data needs to be interpolated to recover
all the three channels. Traditional approaches to demosaicing involve the use
of the Bayer pattern, which tiles a fixed [B, G; G, R] pattern over the image
and use variants of algorithms like edge-directed interpolation which are
tuned to the Bayer pattern. The Matlab demosaic function, for instance, uses
[22], which takes a gradient-corrected bilinear interpolated approach.
However, recently a case has been made for panchromatic demosaicing [19],
where we sense a linear combination of the three channels. In this paper, we
use techniques from compressive recovery to perform demosaicing, i.e. to
reconstruct the original color. However, it turns out that the Bayer pattern
has very high mutual coherence, so it is unsuitable for compressive recovery.
Here, we propose to design the mosaic patterns by minimizing mutual coherence.

We design $8\times 8$ codes for linearly combining the three channels using
our method and visually compare overlapping reconstructions.

Figure 7: Demosaicing. Left: original images, middle: reconstructions with
random matrices, right: reconstructions with non-circularly designed matrices.
RRMSEs = 0.0189 and 0.0189 for upper image with random matrices and non-
circularly designed matrices respectively. RRMSEs = 0.0517 and 0.0510 for
lower image with random matrices and non-circularly designed matrices
respectively

As Figs. 7 and 8 show, results from the designed case are more faithful to the
ground-truth than the random reconstructions are. The random reconstructions
show (more) color artifacts, especially in areas where the input image varies
a lot (car headlights in the top image, around parrot eyes in the bottom). Our
designed codes do not show as many color artifacts. The relative root mean
square errors don’t differ much for these two cases, but subtle details of
color are better preserved by our matrices. (Indeed, it is well-known in the
image processing community that object quality metrics such as RMSE can differ
from human perception). In Fig. 8, notice in the first case the green
artifacts near car headlights and the leftmost cyclist in the results with
random matrices, while that area is better reconstructed with our designed
matrices. The car headlight area on the car at the right is also better
reconstructed by our matrices. In the bottom, notice less color artifacts in
the densely-varying area near the eye and on the bottom part of the beak.

|
---|---
|
|  |
---|---
|  |
---|---
|  |
---|---
|  |
---|---
Figure 8: Demosaicing closeups: In each group of images - original image (top
left), reconstructions with random matrices (top right), with circularly
optimized matrices (bottom right), and with non-circularly optimized matrices
(bottom left).

###  3.5 Video compressed sensing

We now visually validate that our designed matrices perform better than
positive random matrices. We design $8\times 8$ codes and tile them,
reconstructing patchwise with overlapping patches. An example of this running
on six not necessarily close frames in a video is shown in Fig. 10 (Fig. 9
shows an example for $T=2$). This is a challenging test case involving large
object motion or significant structural changes in the scene. Ghosting
artifacts marked out in Fig. 10 in white boxes in the random matrix
reconstructions are absent or lower in the designed matrix reconstructions.
These outputs show that on the large scale, we do as well as random matrices
for low $T$ and better for high $T$. For a small scale comparison, see Figs.
15 and 16.

Figure 9: Optimized output from combining two close images. Left: inputs,
middle, right: reconstructions, in order, with random matrices and non-
circularly optimized matrices. RRMSEs = 0.0277, 0.0139 for frame 1 and 0.0280,
0.0136 for frame 2 for random matrices and non-circularly optimized matrices
respectively

Figure 10: Optimized output from combining six not necessarily close images.
Left: inputs, middle: reconstructions with random matrices, right:
reconstructions with non-circularly optimized matrices. Respective RRMSES =
0.1053 and 0.0574 for frame 1, 0.1051 and 0.0532 for frame 2, 0.1016 and
0.0558 for frame 3, 0.1090 and 0.607 for frame 4, 0.1012 and 0.0496 for frame
5, 0.1100 and 0.0650 for frame 6 for random matrices and non-circularly
optimized matrices.

Finally, we perform a numerical comparison similar to the one in Fig. 11. The
resulting error map is shown in Fig. 12. On an average, we see that we perform
better than the random case at high $T$ and $s$, which is the region we aim to
optimize for. For a comparison, the difference between the two errors is shown
in Fig. 13. The RRMSE difference does produce significant changes in subtle
texture as seen in Figs. 15 and 16. In the deep blue low $T$ regions in Fig.
13, for instance near $s=0.3$ (see blue arrow), where random matrices seem to
outperform us, the reconstruction errors are small enough for both random and
designed matrices to cause no significant depreciation in our reconstruction
quality in comparison to random.

Figure 11: Average relative root mean square errors in our scheme as a
function of $s$ and $T$ with positive random matrices. Compare with Figures 12
and 13. Figure 12: Error map for optimized codes as a function of $s$ and $T$.
Compare with Figures 11 and 13. Figure 13: Average RRMSE using random codes
minus average RRMSE using optimized codes, as a function of $s$ and $T$.
Notice that the designed matrices outperform random matrices in most places
except the blue regions, where both types of matrices yield low RMSE.

Similar to the above analysis for non-circularly designed matrices, we
validate our circularly designed matrices visually. Following the same
conventions, here is an output for the $T=2$ case [Fig. 14].

Figure 14: Circularly optimized output from combining two close images. Left
to right: reconstructions, in order, with random matrices and circularly
optimized matrices. RRMSEs = 0.0277, 0.0139 for frame 1; and 0.0280, 0.0136
for frame 2 for random matrices and circularly optimized matrices respectively

We now look at reconstructions from random and both classes of our designed
matrices on a small scale. As a first example, we show a close-up from the car
video sequence shown earlier [Fig. 15]. Note, to start off, that the
reconstruction of the numberplate and headlight area is much clearer in our
case than the random matrix case. Further, notice the presence of major
ghosting in the random case, especially near the rear-view mirrors, bonnet
(marked by arrows) and headlights (marked by boxes), while our reconstructions
remain free of these artifacts. Adding circular optimization to the picture
further improves image quality especially in the bonnet area, where the non-
circular reconstruction is slightly splotchy. Next, in Fig. 16, which is a
smaller part of the same image, the superiority of our reconstruction is
clearer, with the circular optimization smoothing out blotchier parts of the
bonnet.

Figure 15: Close-ups showing subtle texture preservation with optimized
matrices, example 1. Left to right: inputs, reconstructions, in order, with
random matrices, non-circularly optimized matrices and circularly optimized
matrices

Figure 16: Close-ups showing subtle texture preservation with optimized
matrices, example 2. Left to right: inputs, reconstructions, in order, with
random matrices, non-circularly optimized matrices and circularly optimized
matrices

###  3.6 Comparison with Other Techniques

We emphasize that a direct comparison with earlier methods [12, 11, 8, 10, 14]
is not possible as they do not tackle the same computational problem that we
deal with in this paper - namely these earlier works do not impose the
structure that we impose on the designed matrices. In fact, the aforementioned
techniques cannot be modified to handle the special structure of the matrices
we deal with. For example, [11, 14] provides a closed form expression for the
sensing matrix in the form of an SVD, thereby making it difficult to impose
additional constraints.

###  3.7 Reproducible Research

Most code used in generating results and optimizing sensing matrices in this
paper lives in the Bitbucket repository at alankarkotwal/coded-sourcesep [23].
Gradient descent lives in the src/descent folder, circularly-symmetric
gradient descent in src/descent-circular and reconstruction code in
src/circular.

##  4 Conclusion

We cast the video compressed sensing or color image demosaicing problem as one
of separation of coded linear combinations of signals sparse in a given basis.
We evaluated this scheme and found it works well for low sparsity levels, and
yields reasonably visually good reconstructions. However, especially at high
$T$, we found that random matrices aren’t good enough, motivating design of
sensing matrices based on optimizing mutual coherence, a well-known measure of
matrix ‘quality’ in compressive recovery.

We then provided an analytical expression for the coherence of the sensing
matrix in the coded source separation scheme and optimized for coherence using
these. Results showed better quality and less ghosting visually, and less
error numerically. However as image size increased, the optimization problem
became rapidly intractable, so we settled for optimizing small masks such that
they have small coherence in all circular permutations, so they can be tiled
for overlapping patch-wise reconstruction. We highlighted the importance of
these circular permutations for improving reconstruction performance.

Future work will involve matrix design using criteria that make explicit use
of specific image priors as well as likelihood functions that reflect specific
noise models occuring in imaging applications. Our method can also be extended
to handle multi-spectral image demosaicing or reconstruction.

##  Appendix A Derivation of coherence expression

Recalling our definitions, we call the index varying from $1$ to $T$ as $\mu$
or $\nu$, and the index varying from $1$ to $n$ as $\alpha$, $\beta$ or
$\gamma$. The $\mu^{\text{th}}$ block of $\boldsymbol{\Phi}$ is thus
$\boldsymbol{\Phi_{\mu}}$. Let the $\beta^{\text{th}}$ diagonal element of
$\boldsymbol{\Phi_{\mu}}$ be $\phi_{\mu\beta}$. Define the
$\alpha^{\text{th}}$ column of $\boldsymbol{D^{T}}$ to be
$\boldsymbol{d_{\alpha}}$. Thus, the Gram matrix
$\boldsymbol{\tilde{M}}=\boldsymbol{\Psi^{T}\Phi^{T}\Phi\Psi}$ has the block
structure

$\displaystyle\boldsymbol{\tilde{M}_{\mu\nu}}$
$\displaystyle=\boldsymbol{D^{T}\Phi_{\mu}^{T}\Phi_{\nu}D}$
$\displaystyle=\boldsymbol{D^{T}\Phi_{\mu}\Phi_{\nu}D}$
$\displaystyle=\begin{pmatrix}\boldsymbol{d_{1}}&\boldsymbol{d_{2}}&\ldots&\boldsymbol{d_{n}}\end{pmatrix}\begin{pmatrix}\phi_{\mu
1}\phi_{\nu 1}&0&\ldots&0\\\ 0&\phi_{\mu 2}\phi_{\nu 2}&\ldots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&\ldots&\phi_{\mu n}\phi_{\nu
n}\end{pmatrix}\begin{pmatrix}\boldsymbol{d_{1}^{T}}\\\
\boldsymbol{d_{2}^{T}}\\\ \vdots\\\ \boldsymbol{d_{n}^{T}}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\boldsymbol{d_{1}}&\boldsymbol{d_{2}}&\ldots&\boldsymbol{d_{n}}\end{pmatrix}\begin{pmatrix}\phi_{\mu
1}\phi_{\nu 1}\boldsymbol{d_{1}^{T}}\\\ \phi_{\mu 2}\phi_{\nu
2}\boldsymbol{d_{2}^{T}}\\\ \vdots\\\ \phi_{\mu n}\phi_{\nu
n}\boldsymbol{d_{n}^{T}}\end{pmatrix}$
$\displaystyle=\sum_{\alpha=1}^{n}\phi_{\mu\alpha}\phi_{\nu\alpha}\boldsymbol{d_{\alpha}d_{\alpha}^{T}}$

The $\beta\gamma^{\text{th}}$ element of $\boldsymbol{\tilde{M}_{\mu\nu}}$,
thus, is (defining $d_{\alpha}(\beta)$ as the $\beta^{\text{th}}$ element of
$\boldsymbol{d_{\alpha}}$)

$\displaystyle\tilde{M}_{\mu\nu}(\beta\gamma)$
$\displaystyle=\sum_{\alpha=1}^{n}\phi_{\mu\alpha}\phi_{\nu\alpha}d_{\alpha}(\beta)d_{\alpha}(\gamma)$
(17)

Now we need to normalize the columns of $\boldsymbol{\Phi\Psi}$. Squared
column norms are diagonal elements of $\boldsymbol{\tilde{M}_{\mu\nu}}$. So
the product of the squared norms of the $\beta^{\text{th}}$ column of the
$\mu^{\text{th}}$ block and the $\gamma^{\text{th}}$ column of the
$\nu^{\text{th}}$ block is (call this $\xi_{\mu\nu}^{2}(\beta\gamma)$)

$\displaystyle\xi_{\mu\nu}^{2}(\beta\gamma)$
$\displaystyle=\left(\sum_{\alpha=1}^{n}\phi_{\mu\alpha}^{2}d^{2}_{\alpha}(\beta)\right)\left(\sum_{\tau=1}^{n}\phi_{\nu\tau}^{2}d^{2}_{\tau}(\gamma)\right)$
(18)

Let the normalized Gram matrix be $\boldsymbol{M}$. Thus, following the same
conventions as above (define the numerator of the expression to be
$\chi_{\mu\nu}(\beta\gamma)$),

$\displaystyle M_{\mu\nu}(\beta\gamma)$
$\displaystyle=\frac{\sum_{\alpha=1}^{n}\phi_{\mu\alpha}\phi_{\nu\alpha}d_{\alpha}(\beta)d_{\alpha}(\gamma)}{\sqrt{\left(\sum_{\alpha=1}^{n}\phi_{\mu\alpha}^{2}d^{2}_{\alpha}(\beta)\right)\left(\sum_{\tau=1}^{n}\phi_{\nu\tau}^{2}d^{2}_{\tau}(\gamma)\right)}}=\frac{\chi_{\mu\nu}(\beta\gamma)}{\xi_{\mu\nu}(\beta\gamma)}$
(19)

Finally, using the square soft-max function to deal with the max in the
coherence expression, we get the squared soft coherence $\mathcal{C}$ to be

$\displaystyle\mathcal{C}=\frac{1}{\theta}\log\left[\sum_{\mu=1}^{T}\sum_{\nu=1}^{\mu-1}\sum_{\beta=1}^{n}\sum_{\gamma=1}^{n}e^{\theta
M_{\mu\nu}^{2}(\beta\gamma)}+\sum_{\mu=1}^{T}\sum_{\beta=1}^{n}\sum_{\gamma=1}^{\beta-1}e^{\theta
M_{\mu\mu}^{2}(\beta\gamma)}\right]$ (20)

In the above, the first term corresponds to all ($\mu>\nu$) blocks that are
‘below’ the block diagonal. Here, we consider all terms in the given block for
the maximum. The second term corresponds to ($\mu=\nu$) blocks on the block
diagonal. Here, we consider only consider ($\beta>\gamma$) below-diagonal
elements for the maximum.

##  Appendix B Derivation of coherence derivatives

Differentiating the expression for the squared soft coherence above, we get

$\begin{split}\frac{d\mathcal{C}(\Phi)}{d\phi_{\delta\epsilon}}=\frac{1}{\theta
e^{\theta\mathcal{C}(\Phi)}}\left[\sum_{\mu=1}^{T}\sum_{\nu=1}^{\mu-1}\sum_{\beta=1}^{n}\sum_{\gamma=1}^{n}2\theta
e^{\theta
M_{\mu\nu}^{2}(\beta\gamma)}M_{\mu\nu}(\beta\gamma)\frac{dM_{\mu\nu}(\beta\gamma)}{d\phi_{\delta\epsilon}}\right.\\\
\left.+\sum_{\mu=1}^{T}\sum_{\beta=1}^{n}\sum_{\gamma=1}^{\beta-1}2\theta
e^{\theta M_{\mu\mu}^{2}(\beta\gamma)}M_{\mu\mu}(\beta\gamma)\frac{\theta
M_{\mu\mu}(\beta\gamma)}{d\phi_{\delta\epsilon}}\right]\end{split}$ (21)

Next, we calculate the derivatives in the above equation,
${dM_{\mu\nu}(\beta\gamma)}/{d\phi_{\delta\epsilon}}$. Define the numerator of
the expression for $M_{\mu\nu}(\beta\gamma)$ as $\chi_{\mu\nu}(\beta\gamma)$,
and thus,
$M_{\mu\nu}(\beta\gamma)=\chi_{\mu\nu}(\beta\gamma)/\xi_{\mu\nu}(\beta\gamma)$.
Clearly,

$\displaystyle\frac{dM_{\mu\nu}(\beta\gamma)}{d\phi_{\delta\epsilon}}=\frac{\xi_{\mu\nu}(\beta\gamma)\frac{d\chi_{\mu\nu}(\beta\gamma)}{d\phi_{\delta\epsilon}}-\chi_{\mu\nu}(\beta\gamma)\frac{d\xi_{\mu\nu}(\beta\gamma)}{d\phi_{\delta\epsilon}}}{\xi_{\mu\nu}(\beta\gamma)^{2}}$
(22)

Next,

$\displaystyle\frac{d\chi_{\mu\nu}(\beta\gamma)}{d\phi_{\delta\epsilon}}$
$\displaystyle=\frac{d}{d\phi_{\delta\epsilon}}\sum_{\alpha=1}^{n}\phi_{\mu\alpha}\phi_{\nu\alpha}d_{\alpha}(\beta)d_{\alpha}(\gamma)$
$\displaystyle=\sum_{\alpha=1}^{n}d_{\alpha}(\beta)d_{\alpha}(\gamma)\frac{d}{d\phi_{\delta\epsilon}}\left(\phi_{\mu\alpha}\phi_{\nu\alpha}\right)$

Notice that a term in the above summation can be non-zero only if
$\alpha=\epsilon$. Thus,

$\displaystyle\frac{d\chi_{\mu\nu}(\beta\gamma)}{d\phi_{\delta\epsilon}}$
$\displaystyle=d_{\epsilon}(\beta)d_{\epsilon}(\gamma)\frac{d}{d\phi_{\delta\epsilon}}\left(\phi_{\mu\epsilon}\phi_{\nu\epsilon}\right)$
$\displaystyle=d_{\epsilon}(\beta)d_{\epsilon}(\gamma)\left(\phi_{\mu\epsilon}\frac{d\phi_{\nu\epsilon}}{d\phi_{\delta\epsilon}}+\frac{d\phi_{\mu\epsilon}}{d\phi_{\delta\epsilon}}\phi_{\nu\epsilon}\right)$

Now, notice that ${d\phi_{\mu\epsilon}}/{d\phi_{\delta\epsilon}}$ is non-zero
only if $\mu=\epsilon$. Denote by $\uparrow_{\mu\epsilon}$ the Kronecker delta
function, which is 1 only if $\mu=\epsilon$, 0 otherwise. Then,

$\displaystyle\frac{d\chi_{\mu\nu}(\beta\gamma)}{d\phi_{\delta\epsilon}}$
$\displaystyle=d_{\epsilon}(\beta)d_{\epsilon}(\gamma)\left(\phi_{\mu\epsilon}\uparrow_{\nu\delta}+\uparrow_{\mu\delta}\phi_{\nu\epsilon}\right)$
(23)

Next,

$\displaystyle\frac{d\xi_{\mu\nu}(\beta\gamma)}{d\phi_{\delta\epsilon}}$
$\displaystyle=\frac{d}{d\phi_{\delta\epsilon}}\sqrt{\left(\sum_{\alpha=1}^{n}\phi_{\mu\alpha}^{2}d^{2}_{\alpha}(\beta)\right)\left(\sum_{\tau=1}^{n}\phi_{\nu\tau}^{2}d^{2}_{\tau}(\gamma)\right)}$
$\displaystyle=\frac{1}{2\xi_{\mu\nu}(\beta\gamma)}\frac{d}{d\phi_{\delta\epsilon}}\left(\sum_{\alpha=1}^{n}\phi_{\mu\alpha}^{2}d^{2}_{\alpha}(\beta)\sum_{\tau=1}^{n}\phi_{\nu\tau}^{2}d^{2}_{\tau}(\gamma)\right)$
$\displaystyle=\frac{1}{2\xi_{\mu\nu}(\beta\gamma)}\left[\sum_{\alpha=1}^{n}\phi_{\mu\alpha}^{2}d^{2}_{\alpha}(\beta)\frac{d}{d\phi_{\delta\epsilon}}\left(\sum_{\tau=1}^{n}\phi_{\nu\tau}^{2}d^{2}_{\tau}(\gamma)\right)\right.$
$\displaystyle\quad\quad\quad\quad\quad\left.+\sum_{\tau=1}^{n}\phi_{\nu\tau}^{2}d^{2}_{\tau}(\gamma)\frac{d}{d\phi_{\delta\epsilon}}\left(\sum_{\alpha=1}^{n}\phi_{\mu\alpha}^{2}d^{2}_{\alpha}(\beta)\right)\right]$

Again, a term in one of the above summations is non-zero only if $\alpha$ or
$\tau$ is the same as $\epsilon$. Thus,

$\displaystyle\frac{d}{d\phi_{\delta\epsilon}}\left(\sum_{\alpha=1}^{n}\phi_{\mu\alpha}^{2}d^{2}_{\alpha}(\beta)\right)$
$\displaystyle=2\phi_{\mu\epsilon}d_{\epsilon}^{2}(\beta)\uparrow_{\mu\delta}$

Thus,

$\displaystyle\frac{d\xi_{\mu\nu}(\beta\gamma)}{d\phi_{\delta\epsilon}}$
$\displaystyle=\frac{1}{\xi_{\mu\nu}(\beta\gamma)}\left[\phi_{\mu\epsilon}d_{\epsilon}^{2}(\beta)\uparrow_{\mu\delta}\sum_{\tau=1}^{n}\phi_{\nu\tau}^{2}d^{2}_{\tau}(\gamma)+\phi_{\nu\epsilon}d_{\epsilon}^{2}(\gamma)\uparrow_{\nu\delta}\sum_{\alpha=1}^{n}\phi_{\mu\alpha}^{2}d^{2}_{\alpha}(\beta)\right]$
(24)

This completes the calculation of derivatives.

## References

## References

  * [1] Y. Hitomi, J. Gu, M. Gupta, T. Mitsunaga, S. K. Nayar, Video from a single coded exposure photograph using a learned over-complete dictionary, in: 2011 International Conference on Computer Vision, 2011, pp. 287–294.
  * [2] A. A. Moghadam, M. Aghagolzadeh, M. Kumar, H. Radha, Compressive framework for demosaicing of natural images, IEEE Trans. Image Processing 22 (6) (2013) 2356–2371.
  * [3] D. Reddy, A. Veeraraghavan, R. Chellappa, P2C2: programmable pixel compressive camera for high speed imaging, in: The 24th IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2011, Colorado Springs, CO, USA, 20-25 June 2011, 2011, pp. 329–336.
  * [4] T. T. Cai, L. Wang, Orthogonal matching pursuit for sparse signal recovery with noise, IEEE Transactions on Information Theory 57 (7) (2011) 4680–4688.
  * [5] A. Veeraraghavan, D. Reddy, R. Raskar, Coded strobing photography: Compressive sensing of high speed periodic videos, IEEE Transactions on Pattern Analysis and Machine Intelligence 33 (4) (2011) 671–686.  doi:10.1109/TPAMI.2010.87.
  * [6] C. Studer, R. Baraniuk, Stable restoration and separation of approximately sparse signals, Applied and Computational Harmonic Analysis 37 (1) (2014) 12 – 35.
  * [7] T. T. Cai, L. Wang, G. Xu, New bounds for restricted isometry constants, IEEE Transactions on Information Theory 56 (9) (2010) 4388–4394.
  * [8] B. Li, L. Zhang, T. Kirubarajan, S. Rajan, A projection matrix design method for {MSE} deduction in adaptive compressive sensing, Signal Processing 141 (2017) 16 – 27.
  * [9] M. Pereira, L. Lovisolo, E. da Silva, M. L. de Campos, On the design of maximally incoherent sensing matrices for compressed sensing using orthogonal bases and its extension for biorthogonal bases case, Digital Signal Processing 27 (2014) 12 – 22.
  * [10] V. Abolghasemi, S. Ferdowsi, S. Sanei, A gradient-based alternating minimization approach for optimization of the measurement matrix in compressive sensing, Signal Processing 92 (4) (2012) 999 – 1009.
  * [11] J. M. Duarte-Carvajalino, G. Sapiro, Learning to sense sparse signals: Simultaneous sensing matrix and sparsifying dictionary optimization, IEEE Transactions on Image Processing 18 (7) (2009) 1395–1408.
  * [12] M. Elad, Optimized projections for compressed sensing, IEEE Transactions on Signal Processing 55 (12) (2007) 5695–5702.
  * [13] M. Mordechay, Y. Y. Schechner, Matrix optimization for poisson compressed sensing, in: Signal and Information Processing (GlobalSIP), 2014 IEEE Global Conference on, 2014, pp. 684–688.
  * [14] T. Hong, H. Bai, S. Li, Z. Zhu, An efficient algorithm for designing projection matrix in compressive sensing based on alternating optimization, Signal Process. 125 (C) (2016) 9–20.
  * [15] X. Li, H. Bai, B. Hou, A gradient-based approach to optimization of compressed sensing systems, Signal Processing 139 (2017) 49 – 61.
  * [16] W. R. Carson, M. Chen, M. R. D. Rodrigues, R. Calderbank, L. Carin, Communications-inspired projection design with application to compressive sensing, SIAM Journal on Imaging Sciences 5 (4) (2012) 1185–1212.
  * [17] F. Renna, M. R. D. Rodrigues, M. Chen, R. Calderbank, L. Carin, Compressive sensing for incoherent imaging systems with optical constraints, in: 2013 IEEE International Conference on Acoustics, Speech and Signal Processing, 2013, pp. 5484–5488.
  * [18] Y. Weiss, H. S. Chang, W. T. Freeman, Learning Compressed Sensing, in: Allerton Conference on Communication, Control, and Computing, 2007.
  * [19] K. Hirakawa, P. J. Wolfe, Spatio-spectral color filter array design for optimal image recovery, IEEE Transactions on Image Processing 17 (10) (2008) 1876–1890.  doi:10.1109/TIP.2008.2002164.
  * [20] M. Grant, S. Boyd, CVX: Matlab software for disciplined convex programming, version 2.1, http://cvxr.com/cvx (Mar. 2014).
  * [21] A. Rajwade, D. Kittle, T.-H. Tsai, D. Brady, L. Carin, Coded hyperspectral imaging and blind compressive sensing, SIAM J. Imaging Sciences 6 (2) (2013) 782–812.
  * [22] R. Malvar, L. wei He, R. Cutler, High-Quality Linear Interpolation for Demosaicing of Bayer-Patterned Color Images, in: International Conference of Acoustic, Speech and Signal Processing, 2004.
  * [23] A. Kotwal, Implementation, http://bitbucket.org/alankarkotwal/coded-sourcesep/.