extfieldsofstudy:
Computer Science
Medicine

text:
Clinical Decision Support Systems for Comorbidity: Architecture, Algorithms, and Applications

In this paper, we present the design of a clinical decision support system (CDSS) for monitoring comorbid conditions. Specifically, we address the architecture of a CDSS by characterizing it from three layers and discuss the algorithms in each layer. Also we address the applications of CDSSs in a few real scenarios and analyze the accuracy of a CDSS in consideration of the potential conflicts when using multiple clinical practice guidelines concurrently. Finally, we compare the system performance in our design with that in the other design schemes. Our study shows that our proposed design can achieve a clinical decision in a shorter time than the other designs, while ensuring a high level of system accuracy.

Introduction
A clinical decision support system (CDSS) is designed for ubiquitous health monitoring, which provides mobile health service in order to help physicians and doctors monitor patients' conditions anytime and anywhere with the aid of wireless networks [1]. In such an ubiquitous health monitoring system, building on the biosignal data of patients, a CDSS can assess their conditions with a certain algorithm and inform physicians and doctors immediately when an emergency case occurs, being able to reduce the delay for medical treatment of a critical patient [2].
Most work in the design of CDSS has focused on making clinical decisions on one type of disease, either following the clinical guidelines of this specific disease or employing various data mining schemes to extract clinical rules from numerous biosignal data of patients [3][4][5][6][7][8][9][10]. The primary algorithms include -nearest neighbors based algorithms ( -NN), artificial neural networks based algorithms (ANN), decision trees based algorithms (DTs), fuzzy logic based algorithms (FL), and the extension of these algorithms by merging them.
(a) -Nearest Neighbors Based Algorithms. A -nearest neighbors based algorithm is a method for classifying cases by referring to the closest training cases in the feature space. A typical -NN algorithm runs as follows [10]: cases in question are firstly mixed with the other known cases to constitute a case set. Next, all these cases in this set are mapped into a finite-dimension feature space. A dimension in this space characterizes a feature of cases, so a point in this space uniquely represents a case. Finally, all the points in this space are classified into disjoint subsets, and each subset represents one type of cases. The type of cases in question is assigned by the type of known cases which are within the same subset.
(b) Artificial Neural Network Based Algorithms. Quite similar to the structure of a human brain, an ANN based algorithm is composed of a few interconnected layers of neurons and extracts the exact way of connecting layers by training the historical data of patients. Based on the different training schemes, numerous ANN methods are presented in the literature. In the context of wireless health, two main ANN based algorithms are widely employed in a CDSS [9], including clustering, self-organizing map: clustering based training scheme is implemented by extracting a few critical features of cases and then mapping these features into the structure of an ANN [8]; self-organizing map alters the structure of an 2 International Journal of Telemedicine and Applications ANN by adding necessary neurons and deleting unnecessary neurons [7].
(c) Decision Tree Based Algorithms. Generally, a decision tree based algorithm makes a decision through a series of question and answer sessions. Firstly, the tree based algorithm starts with a root node, which represents the most significant feature in determining the condition of a patient. Then, cases with different states of this feature are split into two or more branches. If all the cases fall into one branch, the branch can be denoted as a leaf. The process of splitting is repeated until all the branches are denoted as leaves [5].
(d) Fuzzy Logic Based Systems. The application of fuzzy models into decision support systems is motivated by the uncertainty of thresholds to make a decision on the patient condition [4]. For example, if the normal range of heart rate is 80-90 times/min, a heart rate of 79 or 91 times/min should not be certainly determined as abnormal [3]. Instead, the result should be determined as normal with a certain probability. In general, the fuzzy process is composed of three procedures: (1) Fuzzification, which is characterized as mapping crisp inputs into fuzzy variables by acquiring the membership values of these crisp inputs, the membership functions can take different shapes, and the simplest but quite widely employed is trapezoidal function; (2) inference based on some rules; (3) defuzzification, which is characterized as converting fuzzy variables into crisp outputs with various mapping methods, and one typical defuzzification method is the centroid of area (COA) method.
Also a few CDSSs are designed to implement in the scenario of making clinical decisions on a comorbid patient (e.g., a patient with Atrial Fibrillation in the setting of Wolff Parkinson's White Syndrome), and in such scenarios the decisions on the treatment of multiple diseases might be in conflict with each other. Such a need is especially clear for elderly patients: several studies have shown that about half of people above 65 years old have a comorbid condition [11]. Boyd et al. in [12] discuss the employment of multiple clinical practice guidelines (CPG) to the care of patients in a comorbid condition. By analyzing various hypothetical clinical scenarios of typical clinical diseases and their associated CPGs, the authors conclude that a blind combination of individual CPGs in the presence of comorbid condition may negatively impact a clinical decision due to the potential mutual conflicts between guidelines. They explain this situation by the fact that a CPG is developed by a specialty-dominated committee and thus is not designed to provide help with regard to evaluating and managing comorbid patients. The conclusion in [12] motivates the research on how to develop a computing algorithm in the design of CDSS for making clinical decisions on the patients in a comorbid condition. Research in the stream of computing algorithms in the design of CDSS for patients in a comorbid condition, despite its importance, seems to be in its infancy. Riaño et al. in [13] advocate using an ontological model for guideline personalization. Unfortunately, the development of this model needs to be conducted manually by physicians using a specialized clinical guideline editor instead of automatically by a computer. In [14], the authors propose a methodology as well as its specific implementation for automatically combining multiple clinical guidelines; however, the data of patients need to be manually entered by physicians.
1.1. Motivation. All these above-mentioned papers do not address how to establish a system which can automatically extract the patient data and offer a decision for clinical practice. Also all the aforementioned research focuses on the level of algorithms and do not mention the design of a CDSS at the system level. The importance of designing a CDSS for comorbidity monitoring as well as the lack of efficient design in the system level motivates us to investigate how to design a CDSS for monitoring comorbid patients by merging the clinical guidelines as well as how to assess the system performance.

Main Contributions.
In this paper, we present the design of a CDSS for monitoring comorbid conditions, specifically from the perspectives of the architecture, the algorithms, and the applications of a CDSS. Also we address how to analyze the accuracy of a CDSS in consideration of the potential conflicts when using multiple clinical practice guidelines concurrently.
The primary contributions of this paper rest at the following issues: (i) address the architecture, the algorithms in the design of a CDSS for comorbidity; (ii) present the implementation process of a CDSS for comorbidity; (iii) assess the system accuracy of a CDSS for comorbidity.

Architecture of CDSS for Comorbidity Monitoring
In this section, we address the architecture of a CDSS in our design for comorbidity monitoring, and it is shown in Figure 1. As shown in Figure 1, the architecture of a CDSS is composed of three layers: physical layer, data layer, and application layer.
In the physical layer, the system needs to extract data from various data sources. The primary data source in clinical scenarios is hospital information systems (HIS), which are located in hospitals or medical centers. Also the medical data may be from a few medical sensors, including the medical sensors attached to or even embedded into the body of a patient, for example, Holter and blood pressure sensors. As shown in Figure 2, data from different sources might be in different formats, and we employ a front end processor (FEP) to do the data extraction, transformation, and loading (ETL) process and also convert data into XML or JSON. In consideration of the security issues, we add a virtual private network (VPN) connection into the system, and also we use a firewall and web service server at the end of data center, to ensure that the data center is not directly connected to public networks for security issues.  In the data layer, with the dramatic increase of clinical data, the traditional data processing applications are inadequate to process such a large set of data, so we employ a Hadoop/Hbase based server cluster to store and process data. The data stored in the data center are prepared to offer data APIs to the application layer, and we use a few offthe-shelf tools of Hadoop to achieve the data APIs [15]. The detailed Hadoop ecosystem is shown in Figure 3. Specifically, we employ Pig and Hive to support the application layer by combining them with Hbase [16] and also Sqoop to connect with the traditional relational databases.
In the application layer, the system can support the CDSS for specific clinical diagnosis, and we focus on the diagnosis on the comorbid patients. The CDSS for comorbidity diagnosis needs to merge clinical practice guidelines (CPG) and detect and mitigate the potential conflicts between CPGs. After mitigating all the conflicts, the system can achieve a final decision on the clinical practice. In the design of an application layer, we use a constraint logic programming model to represent all the CPGs and detect the potential conflicts. The detailed algorithms are illustrated in Sections 3 and 4.
Additionally, the issue of security is critical in the application of healthcare, and generally the security issues in Hadoop involve three aspects: authentication, authorization, and network encryption. We can employ a few off-the-shelf (1) Authentication. Kerberos is a secret-key based authentication system, which can provide authentication and Single-Sign-On (SSO) of users and servers in a network. Hadoop uses the Kerberos protocol to authenticate users to Hadoop and to authenticate Hadoop services to each other.
(2) Authorization. Hadoop provides authorization controls to authenticated users with the aid of HDFS file permissions and service-level authorization. Each of the files is associated with a group and has associated permissions of read and write. Based on such file permissions, HDFS controls access to reads and writes to the distributed file system. Service-level authorization is used by a cluster to control the permission of users to submit jobs.
(3) Network Encryption. If a few sensitive data are exchanged in a network, then we need to protect the cluster from eavesdroppers and network sniffers by encrypting these data.

Algorithms of CDSS for Comorbidity Monitoring
As shown in Figure 1, the architecture of a CDSS is composed of three layers: physical layer, data layer, and application layer.
In this section, we present the specific algorithms in each of the layers.

Physical Layer.
In the physical layer, the CDSS extracts data from various data sources and runs the ETL on raw data in order to load the processed data (clean data) into the data center. The detailed algorithm in the physical layer is in Figure 4. The CDSS firstly accesses the data from data sources and extracts data, does the data cleaning, and standardizes the data into the XML or JSON format. Data cleaning includes data transformation, duplication elimination, replacing missing data, and removing outliers. Of course, the process of data cleaning refers to a standard, for example, HL7, and thus the system supports HL7 to transform the data from different sources (in different formats) into the data in a uniform format. Then, the CDSS calls a web service, and once the connection is successfully built up, the formatted clean data would be transmitted to the data center. Usually, the data are transmitted in bulk, so the data transmission would finish International Journal of Telemedicine and Applications 5 PersonExternalInfo * _11112222_ * T9 Name = "Peter" National = "Canada" T8 IDNumber = "11112222" FamilyDisease = "None" T7 Address = "Ottawa" Allergen = "None" * _22223333_ * T6 Name = "Steven" National = "USA" T5 IDNumber = "22223333" FamilyDisease = "Diabetes" T4 Sex = "Male" Allergen = "Sea food" * _33334444_ * T3 Name = "Lee" National = "China" T2 IDNumber = "33334444" FamilyDisease = "None" T1 Sex = "Male" until all of the bulks of data are transmitted. When calling the web service, if the connection is failed, the CDSS would recall the web service 10 times and write the results into the Log if none of 10 tries succeeds. Otherwise, when the connection succeeds, the CDSS would extract more data and repeat the above-mentioned procedures.

Data Layer.
In the data layer, the clinical data are stored in Hbase tables, the design of which is shown in Time stamps are automatically generated by Hbase, and they denote the time of data entry. Column families are usually a set of columns, and one column represents a property of patients. For example, in the first row, the row key "11112222" represents the identification (ID) of such a row. The time stamps "T7," "T8," and "T9" represent the time of data entry. Specifically, the entry of "Name = Peter" in the column "PersonBasicInfo" and the entry of "National = Canada" in the column "PersonExternalInfo" occur at the time stamp of "T9"; the entry of "IDNumber = 11112222" in the column "PersonBasicInfo" and the entry of "FamilyDisease = None" in the column "PersonExternalInfo" occur at the time stamp of "T8"; the entry of "Address = Ottawa" in the column "PersonBasicInfo" and the entry of "Allergen = None" in the column "PersonExternalInfo" occur at the time stamp of "T7."

Application Layer.
In the application layer, the algorithms of our CDSS are developed with CLP-CPG models.
A CLP-CPG model refers to the development of a combined clinical practice guideline for a specific patient with the aid of constraint logic programming (CLP) [14]. The CLP-CPG model is motivated by the fact that the diagnosis and treatment plan for a specific patient may not be made in consideration of only one CPG; instead, the consideration of combining multiple CPGs is necessary to develop a treatment for patients, especially for those with comorbidity [17]. However, the challenge of combining multiple CPGs is the possible conflicts among multiple CPGs, these conflicts constituting a point of contention (POC) [18]. The CLP-CPG model attempts to mitigate all the POCs by modifying the details of a CPG with the approval of physicians, and the modification of CPGs leads to a modified CLP-CPG model. The modified CLP-CPG model can provide a complete diagnosis and treatment plan to physicians if all the POCs are mitigated. Otherwise, physicians can achieve the conflicts between the CPGs, and they will manually make a complete plan based on the information offered by the model. The work flow of a CLP-CPG model can be shown in Figure 5, given that CPGs are considered.
In the following, we show the detailed process of establishing a CDSS for comorbidity monitoring by providing the example of CDSS for both hypertension (HTN) and deep vein thrombosis (DVT), but our design of CDSS for comorbidity is applicable to the monitoring of comorbid patients with other diseases. The work flow is composed of 3 steps (see details in Section 4): (1) For each of CPGs (in the format of text or graph), the system can establish a corresponding CPG-CLP model to represent the specific CPG. (2) Usually, the CDSS for comorbidity requires multiple CPGs, which might be in conflict with each other. Thus, the second step is to combine each of the CPG-CLP models which are relevant to the comorbid patient and investigate the potential conflicts.
(3) When any conflict exists, the CDSS needs to find a solution which can mitigate the conflict. If no conflict exists, the CDSS can output the final plan of making diagnosis decisions on the patient; if all the possible solutions are tried, but cannot mitigate all the conflicts, the CDSS needs to report the conflicts to the physicians. Then, the physicians need to manually make a plan, and the CDSS memorizes the solution by the physicians for future decision making.

Applications of CDSS for Comorbidity Monitoring
In this section we present how to use a CLP-CPG model in a clinical scenario for comorbidity monitoring. Specifically, we take the scenario of a patient suffering from both hypertension (HTN) and deep vein thrombosis (DVT) as an example.
In such example, we use anonymized patient data transcribed from a public data source [19]. Please note that the design of a CDSS for HTN/DVT monitoring can also be applied to the monitoring of the other diseases.

Clinical Guideline.
The clinical guidelines for a patient who has both HTN and DVT are illustrated in Figure 6. Based on the clinical guideline of DVT, the patient would be inquired whether he or she has a history of severe bleeding tendency (sbt). If so, a physician would recommend this patient to use an inferior vena cava filter (ivcf). Otherwise, this patient would be inquired whether he or she has a history of heparin-induced thrombocytopenia (hit). If so, this patient would be required to use anticoagulant agents (ana), and otherwise this patient would be recommended to take warfarin (wa). Finally, after this patient uses ivcf or ana or wa, the physician would transfer this patient to his or her family doctor to follow up (fdfu). Based on the clinical guideline of HTN, a patient would be checked to confirm that the hypertension he or she suffers is controlled (htnc) or uncontrolled (htnuc). If the type of hypertension is in the type of controlled, the patient would be directly transferred to his or her family doctor to follow up. If the type of hypertension is uncontrolled, the patient would be further checked to confirm whether this patient's hypertension is in urgency (htnur) or emergency (htnem). If the patient is in the status of htnur, the patient would be required to take oral antihypertensive agents (oaa) and otherwise to take the IV antihypertensive agents (iaa). Finally, the patient would be transferred to his or her family doctor to follow up. For a patient with both HTN and DVT, he/she must follow the guidelines of both. An adverse interaction between guidelines for HTN and DVT is that a patient who suffers urgency hypertension (htnur) cannot take warfarin (wa) or anticoagulant agents (ana), since wa and ana have negative impact on the treatment of hypertension. In such a case, physicians must replace the wa or ana with the other types of treatment, that is, intra-vena cava filter.

Work Flow of a CDSS for Comorbidity Monitoring.
In this section, we would like to address the detailed design of a CDSS, and usually a CDSS is embedded in a patient device. As the primary component of a patient device, a CDSS is designed by referring to clinical guidelines and can output the diagnosis on patient conditions. The first procedure of designing a CDSS is to combine multiple clinical guidelines by translating each of the clinical guidelines into a logic expression. A clinical guideline is composed of a few clinical rules, and each of the rules can be represented by a logic   (2) Transform Constraints between Guidelines. Taking warfarin (wa) or anticoagulant agents (ana) could result in negative impacts on an urgency hypertension patient. Thus, we need to set a constraint of htnur for the combination of guidelines of HTN and DVT. In this paper, we assume that the constraints are manually analyzed and set by a few physicians.
In the future, when all these constraints, that is, the conflicts between clinical guidelines, can be analyzed and stored in a repository, the CDSS can automatically detect the potential conflicts by searching in the repository. This constraint is not presented in any individual clinical guideline but appears when we need to combine multiple guidelines together for the diagnosis on a comorbid patient. In the form of logic expressions, we can denote this constraint as not (htnur∧wa) and not (htnur ∧ ana). When this constraint is violated, that is, in the setting of htnur ∧ wa, the clinicians must adjust the regular treatment plan for this patient.

Establish a Constraint Logic Programming Model.
Based on Sections 4.1 and 4.2, we can establish a constraint logic programming (CLP) model by transforming rules as well as the constraints between multiple clinical guidelines. Consequently, the problem of designing a CDSS can be characterized as a CLP model, and the detailed model is shown in Table 2. In the CLP model, we can capture the clinical decision by ensuring that at least one of the logic expressions which represent the rules of HTN and DVT is "true." In the stage of implementation, we employ an ECLiPSe platform [20] to represent the above-mentioned CLP model which is established by referring to Table 2. ECLiPSe is a platform designed to support a CLP modelling language-Zinc, which is high-level enough to well characterize most CLP problems and also low-level enough to be easily compatible with many CLP solvers. Also the ECLiPSe offers a special library (called repair) to capture any violations of constraints between CLP expressions. When a constraint in the CLP expressions leads to a conflict, it can be captured by the repair library.
In the following, we present a few sample ECLiPSe codes for capturing the violations of constraints between clinical guidelines: To illustrate the results of the above-mentioned ECLiPSe sample codes running on the platform, we consider the following clinical scenario: a patient takes both htnur and wa; then ECLiPSe will return the results As expected, ECLiPSe captures the violation of constraint Constraint1 (htnur, wa) between clinical guidelines. Then, the CDSS would find a few solutions in order to mitigate the violation of constraints, for example, replacing a few types of medicine by the others. If a CDSS cannot find any solutions, the results would be reported to the physicians who will manually handle the problem.

Analysis on the System Accuracy of a CDSS
The nodes in a CLP-CPG model of a CDSS are composed of decision nodes (i.e., the yellow hexagon box in Figure 6) and process nodes (i.e., the green rectangle box in Figure 6). The errors of a CDSS can only occur at the decision nodes when the physicians make a mistake in diagnosis or recording the testing results. In such a case, the final decision made by a CDSS might also be incorrect. In this section we analyze the system accuracy of a CDSS to investigate the performance of our design of system architecture.
The system accuracy (ACC) refers to the probability of correct decision results of both HTN and DVT, and its value equals the multiplication of HTN and DVT : where HTN represents the probability that the decision for HTN is correct, while DVT represents the probability that the decision for DVT is correct. Given that DVT = (sbt = 1,ŝbt = 1) + (sbt = 0,ŝbt = 0) (hit = 1,ĥit = 1) + (sbt = 0,ŝbt = 0) (hit = 0,ĥit = 0) , where sbt = 1 and sbt = 0 denote that the patient has a history of severe bleeding tendency or not (refer to Figure 6), respectively, whileŝbt = 1 andŝbt = 0 denote that the physicians record that the patient has a history of severe bleeding tendency or not, respectively. Similarly, hit = 1 and hit = 0 denote that the patient has a history of heparininduced thrombocytopenia or not (refer to Figure 6), respectively, whileĥit = 1 andĥit = 0 denote that the physicians record that the patient has a history of heparin-induced thrombocytopenia or not, respectively.
where htnc = 1 denotes that the patient actually suffers the controlled hypertension, andĥtnc = 1 denotes that the patient is diagnosed to suffer the controlled hypertension (refer to Figure 6). Similarly, htnuc = 1 denotes that the patient actually suffers the urgency hypertension, and htnuc = 1 denotes that the patient is diagnosed to suffer the urgency hypertension (refer to Figure 6), respectively. htnem = 1 denotes that the patient actually suffers the emergency hypertension, andĥtnem = 1 denotes that the patient is diagnosed to suffer the emergency hypertension, respectively. The rationale for (2) and (3) is calculating the probability of correct decisions on each rule of Table 2.

Simulation and Results
In consideration of the errors in the clinical diagnosis, we investigate the system accuracy of a CDSS. Also we compare the average running time of achieving a clinical decision by the CDSS in our design with that in the other designs. In the comparison, we employ anonymized patient data transcribed from a public data source [19]. MIT-BIH database is a wellknown standard ECG database, which is used universally for ECG analysis purpose. ECG signals in the database have normal sinus rhythm with a sample rate of 128 Hz and have uniform intervals of 7.8125 ms. We employ all the samples in the above-mentioned database as the dataset to analyze the performance of various systems.  Figure 7 illustrates that both analytical results and simulation results of the average accuracy of a CDSS match well, verifying the accuracy of our analytical results. Also Figure 7 illustrates that our designed system can achieve a 95% of accuracy when the level of SNR reaches 20 dB or above (i.e., the average probability of diagnosis errors is around 1% or below).
As expected, we can find that the system accuracy increases with the values of SNRs. However, when the SNR is increasing, the cost of system may also increase, because it takes cost (both in time and in money) to help physicians or clinical staff reduce their errors in the diagnosis or recording the testing results. Given an acceptable level of system performance, we should search for a necessary SNR to guarantee the performance of a CDSS. For instance, we need to design a CDSS with an acceptable level of system accuracy as 95%; we must guarantee the SNR higher than 20 dB (SNR depends on the communication channel conditions between the transmitter and the receiver. In the simulation, we analyze the accuracy with the change of SNR in order to investigate the accuracy of a system under various channel conditions. However, the level of SNR cannot be manually adjusted.). Such a rule of system design can partially determine whether clinical staff or physicians need to be trained. Thus, the analysis on system accuracy can not only help healthcare staff estimate the accuracy of a system, but enable the administrator of a hospital or medical center to make the plan of training their staff.

Running Time of a CDSS.
In this section, we compare the running time of a CDSS in our design with that in a few current designs. The running time is primarily determined by the design in the physical layer and data layer. More specifically, we investigate the impact of FEP (FEP versus no FEP) in the physical layer as well as Hbase (Hbase versus MySQL) in the data layer on the running time of a CDSS. Figure 8 shows that the running time of a CDSS with Hbase and FEP can dramatically decrease in comparison with that without Hbase and FEP. When the data volume reaches 300000 entries (data volume refers to the number of entries in the data table. The number of entries does not mean the number of users, since a user usually has multiple entries; e.g., a patient may visit a clinic multiple times and therefore have multiple entries of clinical records. Thus, 300000 entries are possible in the real world), the running time of a CDSS with Hbase and FEP is just around a half of that of a CDSS without Hbase and FEP. Also we investigate the proportion of contributions made by the use of Hbase and FEP on the decrease of running time. Figure 8 shows that the design of Hbase in the data layer is the primary contributor, while the design of FEP makes less contributions on reducing running time than the design of Hbase. When the data volume reaches 300000 entries, the use of Hbase can achieve an 8 × 10 4 decrease of running time, while the use of FEP can only attain a 4 × 10 4 decrease of running time.

Discussion and Conclusion
In this paper, we present the detailed design of a CDSS for monitoring comorbid conditions from the perspectives of architecture, algorithms, and applications. Specifically, we establish the architecture of a CDSS with three layers, that is, physical layer, data layer, and application layer, and discuss the detailed algorithms in each of the layers. Also we address the application of our CDSS in a real scenario of monitoring comorbid patients with both DVT and HTN and analyze the system accuracy of a CDSS. Finally, we compare the running time of making a clinical decision in our design with that in the other designs. The results are that our proposed design of CDSS can achieve a clinical decision faster than the other designs, while ensuring a 90%-95% of the system accuracy.
Using typical clinical scenarios, we have shown how our scheme can process two clinical guidelines by developing a computable model to identify the adverse interactions between clinical guidelines. Such a scheme can support the concurrent application of clinical guidelines, and it is an important step in the process of building up an automated CDSS for the management of comorbid patients.
In the following, we present a few open issues we will work on: (1) Expanding our proposed scheme to process more than two CPGs at a time: while the proposed scheme is general enough, the increased complexity in the operations on multiple guidelines requires substantial modifications to the proposed scheme; (2) allowing for dosages of medications for a detailed description of adverse interactions: this requires revisiting definitions of logical models in our proposed scheme; (3) allowing for the automatic extraction of clinical knowledge from diversified sources to update the scheme of processing multiple clinical guidelines.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.


Multi-scale Convolutional Neural Networks for Crowd Counting

Crowd counting on static images is a challenging problem due to scale variations. Recently deep neural networks have been shown to be effective in this task. However, existing neural-networks-based methods often use the multi-column or multi-network model to extract the scale-relevant features, which is more complicated for optimization and computation wasting. To this end, we propose a novel multi-scale convolutional neural network (MSCNN) for single image crowd counting. Based on the multi-scale blobs, the network is able to generate scale-relevant features for higher crowd counting performances in a single-column architecture, which is both accuracy and cost effective for practical applications. Complemental results show that our method outperforms the state-of-the-art methods on both accuracy and robustness with far less number of parameters.

INTRODUCTION
Crowd counting aims to estimate the number of people in the crowded images or videos feed from surveillance cameras. Overcrowding in scenarios such as tourist attractions and public rallies can cause crowd crushes, blockages and even stampedes. It has been much significant to public safety to produce an accurate and robust crowd count estimation using computer vision techniques.
Existing methods of crowd counting can be generally divided into two categories: detection-based methods and regression-based methods.
Detection-based methods generally assume that each person on the crowd images can be detected and located by using the given visual object detector [1, 2,3], and obtain the counting result by accumulating each detected person. However, these methods [4,5,6] need huge computing resource and they are often limited by person occlusions and complex background in practical scenarios, resulting at a relatively low robustness and accuracy. * Xiangmin Xu is the corresponding author.
Regression-based methods regress the crowd count from the image directly. Chan et al. [7] used handcraft features to translate the crowd counting task into a regression problem. Following works [8,9] proposed more kinds of crowd-relevant features including segment-based features, structural-based features and local texture features. Lempitsky et al. [10] proposed a density-based algorithm that obtain the count by integrating the estimated density map.
Recently, deep convolutional neural networks have been shown to be effective in crowd counting. Zhang et al. [11] proposed a convolutional neural network (CNN) to alternatively learn the crowd density and the crowd count. Wang et al. [12] directly used a CNN-based model to map the image patch to its people count value. However, these single-CNNbased algorithms are limited to extract scale-relevant features and hard to address the scale variations on crowd images. Zhang et al. [13] proposed a multi-column CNN to extract multi-scale features by columns with different kernel sizes. Boominathan et al. [14] proposed a multi-network CNN that used a deep and shallow network to improve the spatial resolution. These improved algorithms can relatively suppress the scale variations problem, but they still have two shortages: • Multi-column/network need pre-trained single-network for global optimization, which is more complicated than end-to-end training.
• Multi-column/network introduce more parameters to consume more computing resource, which make it hard for practical application.
In this paper, we propose a multi-scale convolutional neural network (MSCNN) to extract scale-relevant features. Rather than adding more columns or networks, we only introduce a multi-scale blob with different kernel sizes similar to the naive Inception module [15]. Our approach outperforms the state-of-the-art methods on the ShanghaiTech and UCF CC 50 dataset with a small number of parameters.

MULTI-SCALE CNN FOR CROWD COUNTING
Crowd images are usually consisted of various sizes of persons pixels due to perspective distortion. Single-network is hard to counter scale variations with the same sized kernels combination. In [15], a Inception module is proposed to process visual information at various scales and aggregated to the next stage. Motivated by it, we designed a multi-scale convolutional neural network (MSCNN) to learn the scale-relevant density maps from original images.

Multi-scale Network Architecture
An overview of MSCNN is illustrated in Figure. 1, including feature remapping, multi-scale feature extraction, and density map regression. The first convolution layer is a traditional convolutional layer with single-sized kernels to remap the image feature. Multi-Scale Blob (MSB) is a Inceptionlike model (as Figure. 2) to extract the scale-relevant features, which consists of multiple filters with different kernel size (including 9×9, 7×7, 5×5 and 3×3). A multi-layer perceptron (MLP) [16] convolution layer works as a pixel-wise fully connection, which has multiple 1 × 1 convolutional filters to regress the density map. Rectified linear unit (ReLU) [17] is applied after each convolution layer, which works as the activation function of previous convolutional layers except the last one. Since the value in density map is always positive, adding ReLU after last convolutional layer can enhance the density map restoration. Detailed parameter settings are listed in Table 1.

Scale-relevant Density Map
Following Zhang et al. [13], we estimate the crowd density map directly from the input image. To generate a scalerelevant density map with high quality, the scale-adaptive kernel is currently the best choice. For each head annotation of the image, we represent it as a delta function δ (x − x i ) and describe its distribution with a Gaussian kernel G σ so that the density map can be represented as and finally accumulated to the crowd count value. If we assume that the crowd is evenly distributed on the ground plane, the average distance d i between the head x i and its nearest 10 annotations can generally characterize the geometric distortion caused by perspective effect using the Eq. (1), where M is the total number of head annotations in the image and we fix β = 0.3 as [13] empirically.  Figure. 3. The ground truth and estimated density map of test images in ShanghaiTech dataset

Model Optimization
The output from our model is mapped to the density map, Euclidean distance is used to measure the difference between the output feature map and the corresponding ground truth. The loss function that needs to be optimized is defined as Eq.
(2), where Θ represents the parameters of the model while F (X i ; Θ) represents the output of the model. X i and F i are respectively the i th input image and density map ground truth.

EXPERIMENTS
We evaluate our multi-scale convolutional neural network (MSCNN) for crowd counting on two different datasets, which include the ShanghaiTech and UCF CC 50 datasets. The experimental results show that our MSCNN outperforms the state-of-the-art methods on both accuracy and robustness with far less parameter. All of the convolutional neural networks are trained based on Caffe [18].

Evaluation Metric
Following existing state-of-the-art methods [13], we use the mean absolute error (MAE), the mean squared error (MSE) and the number of neural networks parameters (PARAMS) to evaluate the performance on the testing datasets. The MAE and the MSE are defined in Eq. (3) and Eq. (4).
Here N represents the total number of images in the testing datasets, z i andẑ i are the ground truth and the estimated value respectively for the i th image. In general, MAE, MSE and PARAMS can respectively indicate the accuracy, robustness and computation complexity of a method.

The ShanghaiTech Dataset
The ShanghaiTech dataset is a large-scale crowd counting dataset introduced by [13]. It contains 1198 annotated images with a total of 330,165 persons. The dataset consists of 2 parts: Part A has 482 images crawled from the Internet and Part B has 716 images taken from the busy streets. Following [13], both of them are divided into a training set with 300 images and a testing set with the remainder.

Model Training
To ensure a sufficient number of data for model training, we perform data augmentation by cropping 9 patches from each image and flipping them. We simply fix the 9 cropped points as top, center and bottom combining with left, center and right. Each patch is 90% of the original size.
In order to facilitate comparison with MCNN architecture [13], the network was designed similar to the largest column of MCNN but with MSB, of which detailed settings are described in Table 1. All convolutional kernels are initialized with Gaussian weight setting standard deviation to 0.01. As described in Sec. 2.3, we use the SGD optimization with momentum of 0.9 and weight decay as 0.0005.

Results
We compare our method with 4 existing methods on the ShanghaiTech dataset. The LBP+RR method used LBP feature to regress the function between the counting value and the input image. Zhang et al. [11] designed a convolutional network to regress both the density map and the crowd count value from original pixels. A multi-column CNN [13] is proposed to estimate the crowd count value (MCNN-CCR) and crowd density map (MCNN).
In Table 2, the results illustrate that our approach achieves the state-of-the-art performance on the ShanghaiTech dataset. In addition, it should be emphasized that the number of our parameters is far less than other two CNN-based algorithms. MSCNN uses approximately 7× fewer parameters than the state-of-the-art method (MCNN) with higher accuracy and robustness.

The UCF CC 50 Dataset
The UCF CC 50 dataset [19] contains 50 gray scale images with a total 63,974 annotated persons. The number of people range from 94 to 4543 with an average 1280 individuals per image. Following [11,13,14], we divide the dataset into five splits evenly so that each split contains 10 images. Then we use 5-fold cross-validation to evaluate the performance of our proposed method.

Model Training
The most challenging problem of the UCF CC 50 dataset is the limited number of images for training while the people count in the images span too large. To ensure enough number of training data, we perform a data augmentation strategy following [14] by randomly cropping 36 patches with size 225×225 from each image and flipping them as similar in Sec. 3.2.1. We train 5 models using 5 splits of training set. The MAE and the MSE are calculated after all the 5 models obtained the estimated results of the corresponding validation set. During training, the MSCNN model is initialized almost the same as the experiment on the ShanghaiTech dataset except that the learning rate is fixed to be 1e-7 to guarantee the model convergence.

Results
We compared our method on the UCF CC 50 dataset with 6 existing methods. In [20,10,19], handcraft features are used to regress the density map from the input image. Three CNN-based methods [11,14,13] proposed to used multicolumn/network and perform evaluation on the UCF CC 50 dataset. Table 3 illustrates that our approach also achieves the state-of-the-art performance on the UCF CC 50 dataset. Here our parameters number is approximately 5× fewer than the CrowdNet model, demonstrating that our proposed MSCNN can work more accurately and robustly.  [19] 419.5 541.6 -Zhang et al. [11] 467.0 498.5 7.1M CrowdNet [14] 452.5 -14.8M MCNN [13] 377.

CONCLUSION
In this paper, we proposed a multi-scale convolutional neural network (MSCNN) for crowd counting. Compared with the recent CNN-based methods, our algorithm can extract scale-relevant features from crowd images using a single column network based on the multi-scale blob (MSB). It is an end-to-end training method with no requirement for multicolumn/network pre-training works. Our method can achieve more accurate and robust crowd counting performance with far less number of parameters, which make it more likely to extend to the practical application.

REFERENCES
[1] Sheng-Fuu Lin, Jaw-Yeh Chen, and Hung-Xin Chao, "Estimation of number of people in crowded scenes using perspective transformation," IEEE Transactions on


extfieldsofstudy:
Mathematics
Computer Science

s2fieldsofstudy:
Computer Science
Mathematics

text:
Trimmed Density Ratio Estimation

Density ratio estimation is a vital tool in both machine learning and statistical community. However, due to the unbounded nature of density ratio, the estimation procedure can be vulnerable to corrupted data points, which often pushes the estimated ratio toward infinity. In this paper, we present a robust estimator which automatically identifies and trims outliers. The proposed estimator has a convex formulation, and the global optimum can be obtained via subgradient descent. We analyze the parameter estimation error of this estimator under high-dimensional settings. Experiments are conducted to verify the effectiveness of the estimator.

Introduction
estimator. However, it is not clear how to de-bias such an estimator to recover the true density ratio function. [26] took a more direct approach. It estimates a thresholded density ratio by setting up a tolerance t to the density ratio value. All likelihood ratio values bigger than t will be clipped to t. The estimator was derived from Fenchel duality for f -divergence [18]. However, the optimization for the estimator is not convex if one uses log-linear models. The formulation also relies on the non-parametric approximation of the density ratio function (or the log ratio function) making the learned model hard to interpret. Moreover, there is no intuitive way to directly control the proportion of ratios that are thresholded. Nonetheless, the concept studied in our paper is inspired by this pioneering work.
In this paper, we propose a novel method based on a "trimmed Maximum Likelihood Estimator" [17,10]. This idea relies on a specific type of density ratio estimator (called log-linear KLIEP) [30] which can be written as a maximum likelihood formulation. We simply "ignore" samples that make the empirical likelihood take exceedingly large values. The trimmed density ratio estimator can be formulated as a convex optimization and translated into a weighted M-estimator. This helps us develop a simple subgradient-based algorithm that is guaranteed to reach the global optimum.
Moreover, we shall prove that in addition to recovering the correct density ratio under the outlier setting, the estimator can also obtain a "corrected" density ratio function under a truncation setting. It ignores "pathological" samples and recovers density ratio only using "healthy" samples.
Although trimming will usually result a more robust estimate of the density ratio function, we also point out that it should not be abused. For example, in the tasks of two-sample test, a diverging density ratio might indicate interesting structural differences between two distributions.
In Section 2, we explain some preliminaries on trimmed maximum likelihood estimator. In Section 3, we introduce a trimmed DRE. We solve it using a convex formulation whose optimization procedure is explained in Section 4. In Section 5, we prove the estimation error upper-bound with respect to a sparsity inducing regularizer. Finally, experimental results are shown in Section 6 and we conclude our work in Section 7.

Preliminary: Trimmed Maximum Likelihood Estimation
Although our main purpose is to estimate the density ratio, we first introduce the basic concept of trimmed estimator using density functions as examples. Given n samples drawn from a distribution P , i.e., X := x (i) n i=1 i.i.d.
∼ P, x ∈ R d , we want to estimate the density function p(x). Suppose the true density function is a member of exponential family [20], where f (x) is the sufficient statistics, Z(θ) is the normalization function and q(x) is the base measure.
Maximum Likelihood Estimator (MLE) maximizes the empirical likelihood over the entire dataset. In contrast, a trimmed MLE only maximizes the likelihood over a subset of samples according to their likelihood values (see e.g., [10,31]). This paradigm can be used to derive a popular outlier detection method, one-class Support Vector Machine (one-SVM) [24]. The derivation is crucial to the development of our trimmed density ratio estimator in later sections.
Without loss of generality, we can set the log likelihood function as log p(x (i) ; θ) − τ 0 , where τ 0 is a constant. As samples corresponding to high likelihood values are likely to be inliers, we can trim all samples whose likelihood is bigger than τ 0 using a clipping function This optimization has a convex formulation: where is the slack variable measuring the difference between log p x (i) ; θ and τ 0 . However, formulation (2) is not practical since computing the normalization term Z(θ) in (1) is intractable for a general f and it is unclear how to set the trimming level τ 0 . Therefore we ignore the normalization term and introduce other control terms: The 2 regularization term is introduced to avoid θ reaching unbounded values. A new hyper parameter ν ∈ (0, 1] replaces τ 0 to control the number of trimmed samples. It can be proven using KKT conditions that at most 1 − ν fraction of samples are discarded (see e.g., [24], Proposition 1 for details). Now we have reached the standard formulation of one-SVM. This trimmed estimator ignores the large likelihood values and creates a focus only on the low density region. Such a trimming strategy allows us to discover "novel" points or outliers which are usually far away from the high density area.

Trimmed Density Ratio Estimation
In this paper, our main focus is to derive a robust density ratio estimator following a similar trimming strategy. First, we briefly review the a density ratio estimator [27] from the perspective of Kullback-Leibler divergence minimization.

Density Ratio Estimation (DRE)
For two sets of data X p := {x ∼ Q, assume both the densities p(x) and q(x) are in exponential family (1). We know p(x;θp) q(x;θq) ∝ exp [ θ p − θ q , f (x) ] . Observing that the data x only interacts with the parameter θ p − θ q through f , we can keep using f (x) as our sufficient statistic for the density ratio model, and merge two parameters θ p and θ q into one single parameter δ = θ p − θ q . Now we can model our density ratio as where N (δ) is the normalization term that guarantees q(x)r(x; δ)dx = 1 so that q(x)r(x; δ) is a valid density function and is normalized over its domain. Interestingly, despite the parameterization (changing from θ to δ), (4) is exactly the same as (1) where q(x) appeared as a base measure. The difference is, here, q(x) is a density function from which X q are drawn so that N (δ) can be approximated accurately from samples of Q. Let us definê Note this model can be computed for any f even if the integral in N (δ) does not have a closed form .
In order to estimate δ, we minimize the Kullback-Leibler divergence between p and q · r δ : where c is a constant irrelevant to δ. It can be seen that the minimization of KL divergence boils down to maximizing log likelihood ratio over dataset X p . Now we have reached the log-linear Kullback-Leibler Importance Estimation Procedure (log-linear KLIEP) estimator [30,14].

Trimmed Maximum Likelihood Ratio
As stated in Section 1, to rule out the influences of large density ratio, we trim samples with large likelihood ratio values from (6). Similarly to one-SVM in (2), we can consider a trimmed MLEδ = arg max δ where t 0 is a threshold above which the likelihood ratios are ignored. It has a convex formulation: (7) is similar to (2) since we have only replaced p(x; θ) withr(x; δ). However, the ratio modelr(x; δ) in (7) comes with a tractable normalization termN while the normalization term Z in p(x; θ) is in general intractable. Similar to (3), we can directly control the trimming quantile via a hyper-parameter ν: where R(δ) is a convex regularizer. (8) is also convex, but it has n p number of non-linear constraints and the search for the global optimal solution can be time-consuming. To avoid such a problem, one could derive and solve the dual problem of (8). In some applications, we rely on the primal parameter structure (such as sparsity) for model interpretation, and feature engineering. In Section 4, we translate (8) into an equivalent form so that its solution is obtained via a subgradient ascent method which is guaranteed to converge to the global optimum.
One common way to construct a convex robust estimator is using a Huber loss [12]. Although the proposed trimming technique rises from a different setting, it shares the same guiding principle with Huber loss: avoid assigning dominating values to outlier likelihoods in the objective function.
In Section 8.1 in the supplementary material, we show the relationship between trimmed DRE and binary Support Vector Machines [23,4].

Optimization
The key to solving (8) efficiently is reformulating it into an equivalent max min problem. Proposition 1. Assuming ν is chosen such thatt > 0 for all optimal solutions in (8), then δ is an optimal solution of (8) if and only if it is also the optimal solution of the following max min problem: The proof is in Section 8.2 in the supplementary material. We define (δ,ŵ) as a saddle point of (9): where the second ∇ δ means the subgradient if R is sub-differentiable. Now the "trimming" process of our estimator can be clearly seen from (9): The max procedure estimates a density ratio given the currently assigned weights w, and the min procedure trims the large log likelihood ratio values by assigning corresponding w i to 0 (or values smaller than 1 np ). For simplicity, we only consider the cases where ν is a multiple of 1 np . Intuitively, 1 − ν is the proportion of likelihood ratios that are trimmed thus ν should not be greater than 1. Note if we set ν = 1, (9) is equivalent to the standard density ratio Algorithm 1 Gradient Ascent and Trimming Input: X p , X q , ν and step sizes {η it } itmax it=1 ; Initialize δ 0 , w 0 , Iteration counter: it = 0, Maximum number of iterations: it max , Best objective, parameter pair (O best = −∞, δ best , w best ) . while not converged and it ≤ it max do Obtain a sorted set x ) and update (δ best , w best ) accordingly. it = it + 1. end while Output: (δ best , w best ) estimator (6). Downweighting outliers while estimating the model parameter δ is commonly used by robust estimators (See e.g., [3,29]).
The search for (δ,ŵ) is straightforward. It is easy to solve with respect to w or δ while the other is fixed: given a parameter δ, the optimization with respect to w is a linear programming and one of the extreme optimal solutions is attained by assigning weight 1 np to the elements that correspond to the νn p -smallest log-likelihood ratio logr(x (i) , δ). This observation leads to a simple "gradient ascent and trimming" algorithm (see Algorithm 1). In Algorithm 1, In fact, Algorithm 1 is a subgradient method [2,16], since the optimal value function of the inner problem of (9) is not differentiable at some δ where the inner problem has multiple optimal solutions. The subdifferential of the optimal value of the inner problem with respect to δ can be a set but Algorithm 1 only computes a subgradient obtained using the extreme point solution w it+1 of the inner linear programming. Under mild conditions, this subgradient ascent approach will converge to optimal results with diminishing step size rule and it → ∞. See [2] for details. Algorithm 1 is a simple gradient ascent procedure and can be implemented by deep learning softwares such as Tensorflow 1 which benefits from the GPU acceleration. In contrast, the original problem (8), due to its heavily constrained nature, cannot be easily programmed using such a framework.

Estimation Consistency in High-dimensional Settings
In this section, we show how the estimated parameterδ in (10) converges to the "optimal parameters" δ * as both sample size and dimensionality goes to infinity under the "outlier" and "truncation" setting respectively.
In the outlier setting (Figure 1a), we assume X p is contaminated by outliers and all "inlier" samples in X p are i.i.d.. The outliers are injected into our dataset X p after looking at our inliers. For example, hackers can spy on our data and inject fake samples so that our estimator exaggerates the degree of change.
In the truncation setting, there are no outliers. X p and X q are i.i.d. samples from P and Q respectively. However, we have a subset of "volatile" samples in X p (the rightmost mode on histogram in Figure 1b) that are pathological and exhibit large density ratio values.
In the theoretical results in this section, we focus on analyzing the performance of our estimator for high-dimensional data assuming the number of non-zero elements in the optimal δ * is k and use the 1 regularizer, i.e., R(θ) = θ 1 which induces sparsity onδ. The proofs rely on a recent development [35,34] where a "weighted" high-dimensional estimator was studied. We also assume the optimization of δ in (9) was conducted within an 1 ball of width ρ, i.e., Ball(ρ), and ρ is wisely chosen so that the optimal parameter δ * ∈ Ball(ρ). The same technique was used in previous works [15,35].
Notations: We denote w * ∈ R np as the "optimal" weights depending on δ * and our data. To lighten the notation, we shorten the log density ratio model as The proof of Theorem 1, 2 and 3 can be found in Section 8.4, 8.5 and 8.6 in supplementary materials.

A Base Theorem
Now we provide a base theorem giving an upperbound of δ − δ * . We state this theorem only with respect to an arbitrary pair (δ * , w * ) and the pair is set properly later in Section 5.2 and 5.3.
We make a few regularity conditions on samples from Q. Samples of X q should be well behaved in terms of log-likelihood ratio values.
We also assume the Restricted Strong Convexity (RSC) condition on the covariance of X q , i.e., cov( . Note this property has been verified for various different design matrices X q , such as Gaussian or sub-Gaussian (See, e.g., [21,22]).
Theorem 1. In addition to Assumption 1 and 2, there exists coherence between parameter w and δ at a saddle point (δ,ŵ): whereû :=δ − δ * , κ 2 > 0 is a constant and τ 2 (d, n) > 0. It can be shown that if and with probability converging to one.
The condition (11) states that if we swapŵ for w * , the change of the gradient ∇ δ L is limited. Intuitively, it shows that our estimator (9) is not "picky" on w: even if we cannot have the optimal weight assignment w * , we can still use "the next best thing",ŵ to compute the gradient which is close enough. We later show how (11) is satisfied. Note if ∇ δ L(δ * , w * ) ∞ , τ 2 (n, d) converge to zero as n p , n q , d → ∞, by taking λ n as such, Theorem 1 guarantees the consistency ofδ. In Section 5.2 and 5.3, we explore two different settings of (δ * , w * ) that make ||δ − δ * converges to zero.

Consistency under Outlier Setting
Setting: Suppose dataset X p is the union of two disjoint sets G (Good points) and B (Bad points) such that G Figure 1a).
∼ q(x) does not contain any outlier. We set ν = |G| np . The optimal parameter δ * is set such that p(x) = q(x)r(x; δ * ). We set w * i = 1 np , ∀x (i) p ∈G and 0 otherwise.
Remark: Knowing the inlier proportion |G|/n p is a strong assumption. However it is only imposed for theoretical analysis. As we show in Section 6, our method works well even if ν is only a rough guess (like 90%). Loosening this assumption will be an important future work.
This assumption says that the log density ratio model is Lipschitz continuous around its optimal parameter δ * and hence there is a limit how much a log ratio model can deviate from the optimal model under a small perturbation u. As our estimated weightsŵ i depends on the relative ranking ofẑδ(x (i) p ), this assumption implies that the relative ranking between two points will remain unchanged under a small perturbation u if they are far apart. The following theorem shows that if we have enough clearance between "good"and "bad samples", δ converges to the optimal parameter δ * .
Theorem 2. In addition to Assumption 1, 2 and a few mild technical conditions (see Section 8.5 in the supplementary material), Assumptions 3 holds. Suppose min j∈B z δ * (x kλ n with probability converging to 1. It can be seen that δ − δ * = O log d/min(|G|, n q ) if d is reasonably large.

Consistency under Truncation Setting
In this setting, we do not assume there are outliers in the observed data. Instead, we examine the ability of our estimator recovering the density ratio up to a certain quantile of our data. This ability is especially useful when the behavior of the tail quantile is more volatile and makes the standard estimator (6) output unpredictable results.
∼ Q. Truncated densities p δ and q δ are the unbounded densities p and q restricted only on the truncated domain X(δ). Note that the truncated densities are dependent on the parameter δ and ν. We show that under some assumptions, the parameterδ obtained from (9) using a fixed hyperparameter ν will converge to the δ * such that q δ * (x)r(x; δ * ) = p δ * (x). We also define the "optimal" weight assignment p ∈ X(δ * ) and 0 otherwise. Interestingly, the constraint in (9), w * , 1 = ν may not hold, but our analysis in this section suggests we can always find a pair (δ,ŵ) in the feasible region so that δ − δ * converges to 0 under mild conditions. We first assume the log density ratio model and its CDF is Lipschitz continuous.

Assumption 4.
∀u ∈ Ball(ρ), sup Define In this assumption, we define a "zone" T (u, ) near the ν-th quantile t(δ * ) and assume the CDF of our ratio model is upper-bounded over this region. Different from Assumption 3, the RHS of (12) is with respect to 2 norm of u. In the following assumption, we assume regularity on P and Q.
Theorem 3. In addition Assumption 1 and 2 and other mild assumptions (see Section 8.6 in the supplementary material), Assumption 4 and 5 hold.
It can be seen that δ − δ * = O log d/min(|X p (δ * )|, n q ) if d is reasonably large and |X q \X q (δ * )|/n q decays fast.
6 Experiments 6.1 Detecting Sparse Structural Changes between Two Markov Networks (MNs) [14] In the first experiment 2 , we learn changes between two Gaussian MNs under the outlier setting. The ratio between two Gaussian MNs can be parametrized as is the difference between precision matrices. We generate 500 samples as X p and X q using two randomly structured Gaussian MNs. One point [10, . . . , 10] is added as an outlier to X p . To induce sparsity, we set R(∆) = d i,j=1,i≤j |∆ i,j | and fix λ = 0.0938. Then run DRE and TRimmed-DRE to learn the sparse differential precision matrix ∆ and results are plotted on Figure 2a    the ground truth (the position i, j, ∆ * i,j = 0) is marked by red boxes. It can be seen that the outlier completely misleads DRE while TR-DRE performs reasonably well. We also run experiments with two different settings (d = 25, d = 36) and plot True Negative Rate (TNR) -True Positive Rate (TPR) curves. We fix ν in TR-DRE to 90% and compare the performance of DRE and TR-DRE using DRE without any outliers as gold standard (see Figure 2c). It can be seen that the added outlier makes the DRE fail completely while TR-DRE can almost reach the gold standard. It also shows the price we pay: TR-DRE does lose some power for discarding samples. However, the loss of performance is still acceptable.

Relative Novelty Detection from Images
In the second experiment, we collect four images (see Figure 3a) containing three objects with a textured background: a pencil, an earphone and an earphone case. We create data points from these four images using sliding windows of 48 × 48 pixels (the green box on the lower right picture on Figure 3a). We extract 899 features using MATLAB HOG method on each window and construct an 899-dimensional sample. Although our theorems in Section 5 are proved for linear models, here f (x) is an RBF kernel using all samples in X p as kernel basis. We pick the top left image as X p and using all three other images as X q , then run TR-DRE, THresholded-DRE [26], and one-SVM.
In this task, we select high density ratio super pixels on image X p . It can be expected that the super pixels containing the pencil will exhibit high density ratio values as they did not appear in the reference dataset X q while super pixels containing the earphone case, the earphones and the background, repeats similar patches in X q will have lower density ratio values. This is different from a conventional novelty detection, as a density ratio function help us capture only the relative novelty. For TR-DRE, we use the trimming thresholdt as the threshold for selecting high density ratio points.
It can be seen on Figure 3b, 3c and 3d, as we tune ν to allow more and more high density ratio windows to be selected, more relative novelties are detected: First the pen, then the case, and finally the earphones, as the lack of appearance in the reference dataset X q elevates the density ratio value by different degrees. In comparison, we run TH-DRE with top 3% highest density ratio values thresholded, which corresponds to ν = 97% in our method. The pattern of the thresholded windows (shaded in red) in Figure 3e is similar to Figure 3b though some parts of the case are mistakenly shaded. Finally, one-SVM with 3% support vectors (see Figure 3f) does not utilize the knowledge of a reference dataset X q and labels all salient objects in X p as they corresponds to the "outliers" in X p .

Conclusion
We presents a robust density ratio estimator based on the idea of trimmed MLE. It has a convex formulation and the optimization can be easily conducted using a subgradient ascent method. We also investigate its theoretical property through an equivalent weighted M-estimator whose 2 estimation error bound was provable under two high-dimensional, robust settings. Experiments confirm the effectiveness and robustness of the our trimmed estimator.

Appendix
To lighten the notation system, we drop the feature transform f from our equations. The analysis procedure does not change with or without f .

Relationship between Trimmed DRE and Binary SVM [23, 4]
Consider a "symmetrized" extension to the criterion (6): that jointly minimizes the KL divergence from P to Q and from Q to P . Similar to (5), we user 2 to model the ratio q/p: The minus in front of the δ is due to the inversion of the ratio. We can trim the objective function (13) and add a regularization term λR(δ) as we did for the asymmetric one: Proposition 2. If n p = n q , t 0 = 1, R(·) = · 2 2 , the maximizerδ of (14) is the same as the primal solution of a modified SVM using X p and X q as positive and negative class respectively.
It suggests SVM learns an unnormalized and trimmed density ratio function as the decision function.
Proof. By introducing the slack variables as we did in (7). (14) can be rewritten as: After substitutingr andr 2 , (15) can be rewritten as Let n p = n q , t 0 = 1, R(δ) = δ 2 , (16) is an SVM (without a bias term) using X p and X q as positive and negative samples respectively, except the presences of log normalization terms logN (δ) and logN 2 (δ).

Proof of Proposition 1
Proof. To prove the statement, we construct the dual of (8) which has the exactly same form as (9). Denote ∈ R d×nq . The Lagrangian of (8) can be written as where α ∈ R np + , α ∈ R + , α ∈ R + . Now we analyze the KKT condition of the above Lagrangian.
Suppose the optimalt > 0 4 , then α = 0 by the slackness condition that t α = 0. The optimality condition of t in (17) yields: and the optimality condition of yields From (18) and (19), and the slackness condition of optimization (7), we can see x (i) p ∈ X p , if logr(x (i) p ; δ) < t, then i > 0 which leads to α i = 0 (the constraint of i ≥ 0 is ineffective) and thus α i = 1 np . In contrast, if logr(x (i) p ; δ) > t, then we have α i = 0, i = 0 (the constraint of i ≥ 0 is effective). If logr(x (i) p ; δ) falls right on the boundary t, i.e., log r(x (i) p ; δ) = t, α i ∈ [0, 1 np ], since the KKT condition i α i = 0 indicating α i can take non-negative values as long as 1 np · 1 = α + α . We summarize: It can be observed that for (8), (δ = 0, = 0.2 · 1, t = 0.1) is a feasible interior point, and it makes all inequality constraints strict, so the Slater's condition holds for our original primal problem which is also convex. Therefore, the lagrangian dual of the original problem (8) is which is the same as (9) and any points satisfy the KKT condition are both dual (22) and primal (8) optimal.
8.3 Lemma 1 Lemma 1. If Assumptions 1 and 2 hold, then where c is the constant determined by Assumption 2.
Proof. First, we write down -∇ 2 δ L(δ * + u, w * ): , then we have the following inequalities We then invoke Assumption 2 to obtain ν

Proof of Theorem 1
Proof. First, we define the S and S c are the set of indices of non-zero and zero elements of δ * . The cardinlity of S is k. Defineû :=δ − δ * . From the Lemma 1 we can see that, where we set κ 1 := Using Holder's inequality, The introduction of ρ is due to the bounded optimization region. Due to (11), we can convert the above inequality into and because of the setting of λ n , Note that in the first term,δ is obtained at the stationary condition, which implies that there is a subgradient, denoted by ∇ δ 1 , such that (the second ∇ is the subgradient notation) thus we can obtain the upper-bound of ∇ δ L(δ,ŵ),û using the following standard procedure: 1 ) due to convexity of δ 1 and the definition of subgradient.

Proof of Theorem 3
First we verify (11).

Lemma 3. Under Assumptions 4 and 5,
where L is a positive constant. The second term reflects the cost of using the empirical sample to control the ν-th quantile in (28). It can be seen that κ 2 = 2 √ kC CDF C p , τ 2 (n, d) = 2L·Cp √ np . The proof of Lemma 3 uses a fact that only x p in the "zone" T (u, L 1 √ np ) are "dangerous" as they may be mistakenly included or missed out under small perturbation of u. See Section 8.8 in Appendix for the proof.
The reason for this arrangement can be seen from Figure 4. Denote e (j) := exp δ * , x (j) q , s := nq j=1 e (j) ands = i∈G e (i) . Note that Now, as X p (δ * ) and G contains only i.i.d. samples and due to the definition of δ * , we can invoke Lemma 2 again from [14] to bound a(n, d). That is if Assumptions 6 and 7 hold and n q = Ω(n 2 p ), and λ n ≥ where K 1 and c 1 are constants. By taking the extra 2C 2 r |B|Cq nq into account, we obtain Lemma 4.

Proof of Lemma 3
Before we start, we need to define a few empirical counterparts of population quantities used in Section 5.3.
• P n is the empirical distribution of P .

Numerical Analysis
In this section, we present a few numerical experimental results under outlier and truncation setting. In all experiments, we set n p = n q = 5000, λ = 0, and the solution ofδ was obtained using Algorithm 1. We let f (x) = x. Note this is the correct log-ratio model for two Gaussian distributions with different means.
Outlier Setting In this setting, we first generate two "good" datasets G i.i.d.
∼ q(x) = N (−.75, 1). The outlier set B b is generated from a uniform distribution U (−0.4 + b, 0.4 + b), b ∈ [0, 6]. The density ratio estimation is performed using two sets of data: X p,b = {G, B b } and X q , where the cardinality of B is 1000. We repeat the estimation using different choices of b and test its influence on our estimater(x;δ b ). The results can be seen from Figure 6, where the histograms of G and X q are colored red and green respectively. The true density ratio p(x) q(x) is plotted as a dotted line. The histograms of B b with different choices of b was plotted using gradient colors from light blue to purple (we skipped some choices of b for better visualization). For each b, we run the density ratio estimation, and plot learnedr(x;δ b ) using the same gradient color. In the figure, we resaler(x;δ b ) and the true density ratio using a same constant, so they can be plotted alongside with the histogram. Here, we test two methods: the log-Linear KLIEP and the robust estimator proposed in this paper.
It can be easily seen that as b → 6, KLIEP (Figure 6a) tends to significantly overestimate the density ratio and is sensitive to the change of b. The proposed method (Figure 6b), tends to underestimate the density ratio when b is small. However, as b gradually shifts away from the center of X p , leaving the "gap" between inlier and outlier, the robust estimator converges to the true density ratio function.
Truncated Setting In this setting, we generate samples X p i.i.d.
∼ p(x) = N (0, 1) without any contamination. Usually, the ν-th quantile of z(x p ; δ * ) cannot be analytically computed as we do not know the true density ratio. However, it can be seen that for a strictly monotone increasing z(x p , δ * ), samples in the ν-th quantile of z(x p , δ * ) must be in the ν-th quantile of x p since the relative order among x p is preserved after a strictly monotone transform. Thus, we obtain the truncation domain X(δ * ) = {−∞ ≤ x ≤ Φ −1 (ν)}, where Φ −1 is the inverse CDF of N (0, 1). We then generate samples X q ∼ T N (−0.5, 1, −∞, Φ −1 (ν)), where T N is a truncated Gaussian distribution and the last two parameters are the truncation borders. Note we set the mean of Q to be a negative value so that the true density ratiop/q is a monotone increasing function.
The results for ν = 0.5 are plotted on Figure 7 where the true truncated ratio is plotted as a dotted line. It can be seen that the learnedr(x;δ) is fairly close to the true truncated density ratio.