extfieldsofstudy:
Computer Science
Mathematics

text:
Deep Batch Active Learning by Diverse, Uncertain Gradient Lower Bounds

We design a new algorithm for batch active learning with deep neural network models. Our algorithm, Batch Active learning by Diverse Gradient Embeddings (BADGE), samples groups of points that are disparate and high-magnitude when represented in a hallucinated gradient space, a strategy designed to incorporate both predictive uncertainty and sample diversity into every selected batch. Crucially, BADGE trades off between diversity and uncertainty without requiring any hand-tuned hyperparameters. We show that while other approaches sometimes succeed for particular batch sizes or architectures, BADGE consistently performs as well or better, making it a versatile option for practical active learning problems.

INTRODUCTION
In recent years, deep neural networks have produced state-of-the-art results on a variety of important supervised learning tasks. However, many of these successes have been limited to domains where large amounts of labeled data are available. A promising approach for minimizing labeling effort is active learning, a learning protocol where labels can be requested by the algorithm in a sequential, feedback-driven fashion. Active learning algorithms aim to identify and label only maximally-informative samples, so that a high-performing classifier can be trained with minimal labeling effort. As such, a robust active learning algorithm for deep neural networks may considerably expand the domains in which these models are applicable.
How should we design a practical, general-purpose, label-efficient active learning algorithm for deep neural networks? Theory for active learning suggests a version-space-based approach (Cohn et al., 1994;Balcan et al., 2006), which explicitly or implicitly maintains a set of plausible models, and queries examples for which these models make different predictions. But when using highly expressive models like neural networks, these algorithms degenerate to querying every example. Further, the computational overhead of training deep neural networks precludes approaches that update the model to best fit data after each label query, as is often done (exactly or approximately) for linear methods (Beygelzimer et al., 2010;Cesa-Bianchi et al., 2009). Unfortunately, the theory provides little guidance for these models.
One option is to use the network's uncertainty to inform a query strategy, for example by labeling samples for which the model is least confident. In a batch setting, however, this creates a pathological scenario where data in the batch are nearly identical, a clear inefficiency. Remedying this issue, we could select samples to maximize batch diversity, but this might choose points that provide little new information to the model. For these reasons, methods that exploit just uncertainty or diversity do not consistently work well across model architectures, batch sizes, or datasets. An algorithm that performs well when using a ResNet, for For all examples x in U \ S: 1. Compute its hypothetical labelŷ(x) = h θt (x). 2. Compute gradient embedding g x = ∂ ∂θout CE (f (x; θ),ŷ(x))| θ=θt , where θ out refers to parameters of the final (output) layer.

5:
Compute S t , a random subset of U \ S, using the k-MEANS++ seeding algorithm on {g x : x ∈ U \ S} and query for their labels. 3 ALGORITHM BADGE, described in Algorithm 1, starts by drawing an initial set of M examples uniformly at random from U and asking for their labels. It then proceeds iteratively, performing two main computations at each step t: a gradient embedding computation and a sampling computation. Specifically, at each step t, for every x in the pool U , we compute the labelŷ(x) preferred by the current model, and the gradient g x of the loss on (x,ŷ(x)) with respect to the parameters of the last layer of the network. Given these gradient embedding vectors {g x : x ∈ U }, BADGE selects a set of points by sampling via the k-MEANS++ initialization scheme (Arthur and Vassilvitskii, 2007). The algorithm queries the labels of these examples, retrains the model, and repeats.
We now describe the main computations -the embedding and sampling steps -in more detail.
The gradient embedding. Since deep neural networks are optimized using gradient-based methods, we capture uncertainty about an example through the lens of gradients. In particular, we consider the model uncertain about an example if knowing the label induces a large gradient of the loss with respect to the model parameters and hence a large update to the model. A difficulty with this reasoning is that we need to know the label to compute the gradient. As a proxy, we compute the gradient as if the model's current prediction on the example is the true label. We show in Proposition 1 that, assuming a common structure satisfied by most natural neural networks, the gradient norm with respect to the last layer using this label provides a lower bound on the gradient norm induced by any other label. In addition, under that assumption, the length of this hypothetical gradient vector captures the uncertainty of the model on the example: if the model is highly certain about the example's label, then the example's gradient embedding will have a small norm, and vice versa for samples where the model is uncertain (see example below). Thus, the gradient embedding conveys information both about the model's uncertainty and potential update direction upon receiving a label at an example.
The sampling step. We want the newly-acquired labeled samples to induce large and diverse changes to the model. To this end, we want the selection procedure to favor both sample magnitude and batch diversity. Specifically, we want to avoid the pathology of, for example, selecting a batch of k similar samples where even just a single label could alleviate our uncertainty on all remaining (k − 1) samples.
A natural way of making this selection without introducing additional hyperparameters is to sample from a k-Determinantal Point Process (k-DPP; (Kulesza and Taskar, 2011) (2018) showed that in experimental design for least square linear regression settings, learning from samples drawn from a k-DPP can have much smaller mean square prediction error than learning from iid samples. In this process, when the batch size is very low, the selection will naturally favor points with a large length, which corresponds to uncertainty in our space. When the batch size is large, the sampler focuses more on diversity because linear independence, which is more difficult to achieve for large k, is required to make the Gram determinant non-zero.
Unfortunately, sampling from a k-DPP is not trivial. Many sampling algorithms (Kang, 2013;Anari et al., 2016) rely on MCMC, where mixing time poses a significant computational hurdle. The state-of-the-art algorithm of Dereziński (2018) has a high-order polynomial running time in the batch size and the embedding dimension. To overcome this computational hurdle, we suggest instead sampling using the k-MEANS++ seeding algorithm (Arthur and Vassilvitskii, 2007), originally made to produce a good initialization for k-means clustering. k-MEANS++ seeding selects centroids by iteratively sampling points in proportion to their squared distances from the nearest centroid that has already been chosen, which, like a k-DPP, tends to select a diverse batch of high-magnitude samples. For completeness, we give a formal description of the k-MEANS++ seeding algorithm in Appendix A.
Example: multiclass classification with softmax activations. Consider a neural network f where the last nonlinearity is a softmax, i.e. σ(z) i = e z i / K j=1 e z j . Specifically, f is parametrized by θ = (W, V ), where θ out = W = (W 1 , . . . , W K ) ∈ R K×d are the weights of the last layer, and V consists of weights of all previous layers. This means that f (x; θ) = σ(W · z(x; V )), where z is the nonlinear function that maps an input x to the output of the network's penultimate layer. Let us fix an unlabeled sample x and define p i = f (x; θ) i . With this notation, we have , y) for a label y and g x = gŷ x as the gradient embedding in our algorithm, wherê y = argmax i∈[K] p i . Then the i-th block of g x (i.e. the gradients corresponding to label i) is Based on this expression, we can make the following observations: 1. Each block of g x is a scaling of z(x; V ), which is the output of the penultimate layer of the network. In this respect, g x captures x's representation information similar to that of Sener and Savarese (2018). 2. Proposition 1 below shows that the norm of g x is a lower bound on the norm of the loss gradient induced by the example with true label y with respect to the weights in the last layer, that is g x ≤ g y x . This suggests that the norm of g x conservatively estimates the example's influence on the current model. 3. If the current model θ is highly confident about x, i.e. vector p is skewed towards a standard basis vector e j , thenŷ = j, and vector (p i − I(ŷ = i)) K i=1 has a small length. Therefore, g x has a small length as well. Such high-confidence examples tend to have gradient embeddings of small magnitude, which are unlikely to be repeatedly selected by k-MEANS++ at iteration t. Proposition 1. For all y ∈ {1, . . . , K}, let g y x = ∂ ∂W CE (f (x; θ), y). Then Consequently,ŷ = argmin y∈[K] g y x .
Proof. Observe that by Equation (1), The second claim follows from the fact thatŷ = argmax y∈[K] p y .
This simple sampler tends to produce diverse batches similar to a k-DPP. As shown in Figure 1, switching between the two samplers does not affect the active learner's statistical performance but greatly improves its computational performance. Appendix G compares run time and test accuracy for both k-MEANS++ and k-DPP based sampling based on the gradient embeddings of the unlabeled examples. Figure 2 illustrates the batch diversity and average gradient magnitude per selected batch for a variety of sampling strategies. As expected, both k-DPPs and k-MEANS++ tend to select samples that are diverse (as measured by the magnitude of their Gram determinant) and high magnitude. Other samplers, such as furthest-first traversal for k-Center clustering (FF-k-CENTER), do not seem to have this property. The FF-k-CENTER algorithm is the sampling choice of the CORESET approach to active learning, which we describe in the proceeding section (Sener and Savarese, 2018). Appendix F discusses diversity with respect to uncertainty-based approaches.
Appendix B provides further justification for why BADGE yields better updates than vanilla uncertainty sampling in the special case of binary logistic regression (K = 2 and z(x; V ) = x).

EXPERIMENTS
We evaluate the performance of BADGE against several algorithms from the literature. In our experiments, we seek to answer the following question: How robust are the learning algorithms to choices of neural network architecture, batch size, and dataset?
To ensure a comprehensive comparison among all algorithms, we evaluate them in a batch-mode active learning setup with M = 100 being the number of initial random labeled examples and batch size B varying from {100, 1000, 10000}. The following is a list of the baseline algorithms evaluated; the first performs representative sampling, the next three are uncertainty based, the fifth is a hybrid of representative and uncertainty-based approaches, and the last is traditional supervised learning. We consider three neural network architectures: a two-layer Perceptron with ReLU activations (MLP), an 18-layer convolutional ResNet (He et al., 2016), and an 11-layer VGG network (Simonyan and Zisserman, 2014). We evaluate our algorithms using three image datasets, SVHN (Netzer et al., 2011), CIFAR10 (Krizhevsky, 2009 and MNIST (LeCun et al., 1998) 1 , and four non-image datasets from the OpenML repository (#6, #155, #156, and #184). 2 We study each situation with 7 active learning algorithms, including BADGE, making for 231 total experiments.

k-DPP k-means++
For the image datasets, the embedding dimensionality in the MLP is 256. For the OpenML datasets, the embedding dimensionality of the MLP is 1024, as more capacity helps the model fit training data. We fit Learning curves. Here we show examples of learning curves that highlight some of the phenomena we observe related to the fragility of active learning algorithms with respect to batch size, architecture, and dataset.
Often, we see that in early rounds of training, it is better to do diversity sampling, and later in training, it is better to do uncertainty sampling. This kind of event is demonstrated in Figure 3a, which shows CORESET outperforming confidence-based methods at first, but then doing worse than these methods later on. Element Pi,j corresponds roughly to the number of times algorithm i outperforms algorithm j. Column-wise averages at the bottom show overall performance (lower is better).

BADGE ALBL Coreset Conf
In this figure, BADGE performs as well as diversity sampling when that strategy does best, and as well as uncertainty sampling once those methods start outpacing CORESET. This suggests that BADGE is a good choice regardless of labeling budget.
Separately, we notice that diversity sampling only seems to work well when either the model has good architectural priors (inductive biases) built in, or when the data are easy to learn. Otherwise, penultimate layer representations are not meaningful, and diverse sampling can be deleterious. For this reason, CORESET often performs worse than random on sufficiently complex data when not using a convolutional network ( Figure 3b). That is, the diversity induced by unconditional random sampling can often yield a batch that better represents the data. Even when batch size is large and the model has helpful inductive biases, the uncertainty information in BADGE can give it an advantage over pure diversity approaches ( Figure 3c). Comprehensive plots of this kind, spanning architecture, dataset, and batch size are in Appendix C.
Pairwise comparisons. We next show a comprehensive pairwise comparison of algorithms over all datasets (D), batch sizes (B), model architectures (A), and label budgets (L). From the learning curves, it can be observed that when label budgets are large enough, all algorithms eventually reach similar performance, making the comparison between them uninteresting in the large sample limit. For this reason, for each combination of (D, B, A), we select a set of labeling budgets L where learning is still progressing. We experimented with three different batch sizes and eleven dataset-architecture pairs, making the total number of (D, B, A) combinations 3 × 11 = 33. Specifically, we compute n 0 , the smallest number of labels where RAND's accuracy reaches 99% of its final accuracy, and choose label budget The calculation of scores in the penalty matrix P follows the following protocol: For each (D, B, A, L) combination and each pair of algorithms (i, j), we have 5 test errors (one for each repeated run), e 1 i , . . . , e 5 i and e 1 j , . . . , e 5 j respectively. We compute the t-score as t = We use the two-sided t-test to compare pairs of algorithms: algorithm i is said to beat algorithm j in this setting if t > 2.776 (the critical point of p-value being 0.05), and similarly algorithm j beats algorithm i if t < −2.776. For each (D, B, A) combination, suppose there are n D,B,A different values of L. Then, for each L, if algorithm i beats algorithm j, we accumulate a penalty of 1/n D,B,A to P i,j ; otherwise, if algorithm j beats algorithm i, we accumulate a penalty of 1/n D,B,A to P j,i . The choice of the penalty value 1/n D,B,A is to ensure that every (D, B, A) combination is assigned equal influence in the aggregated matrix. Therefore, the largest entry of P is at most 33, the total number of (D, B, A) combinations. Intuitively, each row i indicates the number of settings in which algorithm i beats other algorithms and each column j indicates the number of settings in which algorithm j is beaten by another algorithm.
The penalty matrix in Figure 4 summarizes all experiments, showing that BADGE generally outperforms baselines. Matrices grouped by batch size and architecture in Appendix D show a similar trend.
Cumulative distribution functions of normalized errors. For each (D, B, A, L) combination, we compute the average error for each algorithm i asē i = 1 5 5 l=1 e l i . To ensure that the errors of these algorithms are on the same scale in all settings, we compute the normalized error of every algorithm i, defined as ne i =ē i /ē r , where r is the index of the RAND algorithm. By definition, the normalized errors of the RAND algorithm are identically 1 in all settings. Like with penalty matrices, for each (D, B, A) combination, we only consider a subset of L values from the set M + 2 m−1 B : m ∈ [ log((n 0 − M )/B) ] . We assign a weight proportional to 1/n D,B,A to each (D, B, A, L) combination, where there are n D,B,A different L values for this combination of (D, B, A). We then plot the cumulative distribution functions (CDFs) of the normalized errors of all algorithms: for a value of x, the y value is the total weight of settings where the algorithm has normalized error at most x; in general, an algorithm that has a higher CDF value has better performance.
We plot the generated CDFs in Figures 5, 22 and 23. We can see from Figure 5 that BADGE has the best overall performance. In addition, from Figures 22 and 23 in Appendix E, we can conclude that when batch size is small (100 or 1000) or when an MLP is used, both BADGE and MARG perform best. However, in the regime when the batch size is large (10000), MARG's performance degrades, while BADGE, ALBL and CORESET are the best performing approaches.

RELATED WORK
Active learning is a been well-studied problem (Settles, 2010;Dasgupta, 2011;Hanneke, 2014). There are two major strategies for active learning-representative sampling and uncertainty sampling.
Representative sampling algorithms select batches of unlabeled examples that are representative of the unlabeled set to ask for labels. It is based on the intuition that the sets of representative examples chosen, once labeled, can act as a surrogate for the full dataset. Consequently, performing loss minimization on the surrogate suffices to ensure a low error with respect to the full dataset. In the context of deep learning, Sener and Savarese (2018); Geifman and El-Yaniv (2017)  There is also a large body of literature on batch mode active learning, where the learner is asked to select a batch of samples within each round (Guo and Schuurmans, 2008;Wang and Ye, 2015;Chen and Krause;Wei et al., 2015;Kirsch et al., 2019). In these works, batch selection is often formulated as an optimization problem with objectives based on (upper bounds of) average log-likelihood, average squared loss, etc.
A different query criterion based on expected gradient length (EGL) has been proposed in the as well (Settles et al., 2008). In recent work, Huang et al. (2016) show that the EGL criterion is related to the T -optimality criterion in experimental design. They further demonstrate that the samples selected by EGL are very different from those by entropy-based uncertainty criterion. Zhang et al. (2017a) use the EGL criterion in active sentence and document classification with CNNs. These approaches differ most substantially from BADGE in that they do not take into account the diversity of the examples queried within each batch.
There is a wide array of theoretical articles that focus on the related problem of adaptive subsampling for fully-labeled datasets in regression settings (Han et al., 2016;Wang et al., 2018;Ting and Brochu, 2018). Empirical studies of batch stochastic gradient descent also employ adaptive sampling to "emphasize" hard or representative examples (Zhang et al., 2017b;Chang et al., 2017). These works aim at reducing computation costs or finding a better local optimal solution, as opposed to reducing label costs. Nevertheless, our work is inspired by their sampling criteria, which also emphasize samples that induce large updates to the model.
As mentioned earlier, our sampling criterion has resemblance to sampling from k-determinantal point processes (Kulesza and Taskar, 2011). Note that in multiclass classification settings, our gradient-based embedding of an example can be viewed as the outer product of the original embedding in the penultimate layer and a probability score vector that encodes the uncertainty information on this example (see Section 3). In this view, the penultimate layer embedding characterizes the diversity of each example, whereas the probability score vector characterizes the quality of each example. The k-DPP is also a natural probabilistic tool for sampling that trades off between quality and diversity (See Kulesza et al., 2012, Section 3.1). We remark that concurrent to our work, Bıyık et al. (2019) develops k-DPP based active learning algorithms based on this principle by explicitly designing diversity and uncertainty measures.

DISCUSSION
We have established that BADGE is empirically an effective deep active learning algorithm across different architectures and batch sizes, performing similar to or better than other active learning algorithms. A fundamental remaining question is: "Why?" While deep learning is notoriously difficult to analyze theoretically, there are several intuitively appealing properties of BADGE: Stephen Mussmann and Percy S Liang. Uncertainty sampling is preconditioned stochastic gradient descent on zero-one loss. In Neural Information Processing Systems, 2018.
A THE k-MEANS++ SEEDING ALGORITHM Here we briefly review the k-MEANS++ seeding algorithm by (Arthur and Vassilvitskii, 2007). Its basic idea is to perform sequential sampling of k centers, where each new center is sampled from the ground set with probability proportional to the squared distance to its nearest center. It is shown in (Arthur and Vassilvitskii, 2007) that the set of centers returned is guaranteed to approximate the k-means objective function in expectation, thus ensuring diversity.
Algorithm 2 The k-MEANS++ seeding algorithm (Arthur and Vassilvitskii, 2007) end for return C k .

B BADGE FOR BINARY LOGISTIC REGRESSION
We consider instantiating BADGE for binary logistic regression, where Y = {−1, +1}. Given a linear classifier w, we define the predictive probability of w on x as p w (y|x, θ) = σ(yw · x), where σ(z) = 1 1+e −z is the sigmoid funciton.
Recall thatŷ =ŷ(x) is the hallucinated label: The binary logistic loss of classifier w on example (x, y) is defined as: (w, (x, y)) = ln(1 + exp(−yw · x)). Now, given model w and example x, we defineĝ x = ∂ ∂w (w, (x,ŷ)) = (1 − p w (ŷ|x, θ)) · (−ŷ · x) as the loss gradient induced by the example with hallucinated label, andg x = ∂ ∂w (w, (x, y)) = (1−p w (y|x, θ))·(−y·x) as the loss gradient induced by the example with true label. Suppose that BADGE only selects examples from region S w = {x : w · x = 0}, then as p w (+1|x, θ) = p w (−1|x, θ) = 1 2 , we have that for all x in S w ,ĝ x = s x · g x for some s x ∈ {±1}. This implies that, sampling from a DPP induced byĝ x 's is equivalent to sampling from a DPP induced by g x 's. It is noted in Mussmann and Liang (2018) that uncertainty sampling (i.e. sampling from D |Sw ) implicitly performs preconditioned stochastic gradient descent on the expected 0-1 loss. In addition, it has been shown that DPP sampling over gradients may reduce the variance of the mini-batch stochastic gradient updates (Zhang et al., 2017b); this suggests that BADGE, when restricted its sampling over low-margin regions (S w ), improves over uncertainty sampling by collecting examples that together induce lower-variance updates on the gradient direction of expected 0-1 loss.

C ALL LEARNING CURVES
We plot all learning curves (test accuracy as a function of the number of labeled example queried) in Figures 6 to 12. In addition, we zoom into regions of the learning curves that discriminates the performance of all algorithms in Figures 13 to 19.

D PAIRWISE COMPARISONS OF ALGORITHMS
In addition to Figure 4 in the main text, we also provide penalty matrices (Figures 20 and 21), where the results are aggregated by conditioning on a fixed batch size (100, 1000 and 10000) or on a fixed neural network model (MLP, ResNet and VGG). For each penalty matrix, the parenthesized number in its title is the   total number of (D, B, A) combinations aggregated; as discussed in Section 4, this is also an upper bound on all its entries. It can be seen that uncertainty-based methods (e.g. MARG) perform well only in small batch size regimes (100) or when using MLP models; representative sampling based methods (e.g. CORESET) only perform well in large batch size regimes (10000) or when using ResNet or VGG models. In contrast, BADGE's performance is competitive across all batch sizes and neural network models. The k-centers sampler finds points that are not as diverse or high-magnitude as other samplers. Notice also that k-MEANS++ tends to actually select samples that are both more diverse and higher-magnitude than a k-DPP, a potential pathology of the k-DPP's degree of stochastisity. Among all algorithms, CONF has the largest average norm of gradient embeddings within a batch; however, in OpenML #6, and the first few interations of SVHN, some batches have a log Gram determinant of −∞ (shown as gaps in the curve), which shows that CONF sometimes selects batches that are inferior in diversity.
many gaps in the CONF plot seem to correspond to situations in which CONF performs poorly in terms of accuracy (see Figure 13 for the corresponding learning curve). Both k-DPP and FF-k-CENTER (an algorithm that approximately minimizes k-center objective) select batches that have lower diversity than k-MEANS++  and space complexity of O(k 2 + kd), where τ is the number of sampling steps. We set τ as 5k ln k in our experiment. The comparisons for batch size 10000 are not shown here as the implementation of k-DPP sampling runs out of memory.
It can be seen from the figures that, although k-DPP and k-MEANS++ are based on different sampling criteria, the classification accuracies of their induced active learning algorithm are similar. In addition, when large batch sizes are required (e.g. k = 1000), the running times of k-DPP sampling are generally much higher than those of k-MEANS++.


extfieldsofstudy:
Computer Science
Medicine

text:
Analyzing and Visualizing Knowledge Structures of Health Informatics from 1974 to 2018: A Bibliometric and Social Network Analysis

Objectives This paper aims to provide a theoretical clarification of the health informatics field by conducting a quantitative review analysis of the health informatics literature. And this paper aims to map scientific networks; to uncover the explicit and hidden patterns, knowledge structures, and sub-structures in scientific networks; to track the flow and burst of scientific topics; and to discover what effects they have on the scientific growth of health informatics. Methods This study was a quantitative literature review of the health informatics field, employing text mining and bibliometric research methods. This paper reviews 30,115 articles with health informatics as their topic, which are indexed in the Web of Science Core Collection Database from 1974 to 2018. This study analyzed and mapped four networks: author co-citation network, co-occurring author keywords and keywords plus, co-occurring subject categories, and country co-citation network. We used CiteSpace 5.3 and VOSviewer to analyze data, and we used Gephi 0.9.2 and VOSviewer to visualize the networks. Results This study found that the three major themes of the literature from 1974 to 2018 were the utilization of computer science in healthcare, the impact of health informatics on patient safety and the quality of healthcare, and decision support systems. The study found that, since 2016, health informatics has entered a new era to provide predictive, preventative, personalized, and participatory healthcare systems. Conclusions This study found that the future strands of research may be patient-generated health data, deep learning algorithms, quantified self and self-tracking tools, and Internet of Things based decision support systems.

implementation was started in the 1990s, but they become truly clinically viable after 2000 [3]. Health informatics is also known as healthcare informatics, medical informatics, nursing informatics, or biohealth informatics [4]. Recent advances in healthcare IT, health data standards, Electronic Health Records (EHRs), and health information exchange (HIE) have increased the growth of this scientific field [5], which has attracted interests in both academic and professional contexts.
The health informatics field has grown over the past quarter century, and many efforts have been made to define it in a scientific and formal language [6]. As health informatics has been practically implemented across medical settings [7], it has also become a major focus of scientific research. All of the technological innovations within the domain of health have been made possible through scientific knowledge. However, the healthcare sector is a field in which the development of new scientific knowledge is 'hectic' and technological expansion is profoundly 'rapid' [8]. Health informatics is also a new discipline [9] whose development is linked with technological trends.
Despite its increasing growth and interest among scholars and practitioners, currently, there is no comprehensive characterization of the knowledge structure of this field, and studies on its evolution are scarce. This understanding is necessary to facilitate the relevant technology growth and academic endeavors. Most existing works have evaluated health informatics research by conducting qualitative research methods, such as systematic literature review [10,11], and most have focused on a period of time limited to the last 10 to 20 years. These studies do not offer a complete and objective overview of the current state of research. A systematic literature review answers a specific question and is more focused and narrow in its approach with having a hypothesis to support or reject [12]. Qualitative reviews can be accompanied by several biases, such as publication bias, search bias, and selection bias [13]; such biases threaten objectivity since qualitative analysis requires personal judgment and the expertise of researchers [14]. On the other hand, the stream of publications on health informatics history suggests that this scientific field is multidisciplinary, and "one of the common challenges of multidisciplinary research is a lack of common language" [15]. This paper is intended to address these challenges and provide a theoretical clarification of the health informatics field. A comprehensive quantitative and more objective review of scientific articles can provide academia with valuable information without the intervention of a researcher's bias about the knowledge structures, hidden trends, information flow, and future research orientation [16]. Moreover, it can help the multidisciplinary scientific fields to find their 'common language' . The findings of this research will supplement previous studies that have attempted to portray the thematic evolution of the field. In this study, we reviewed the literature on health informatics by conducting text mining methods, scientometric analysis, and social network visualization to find "the communities embedded in the social network datasets, and moreover, (to analyze) the evolutions of the communities in dynamic networks" [17]. Text mining in social networks enables the discovery of new patterns as well as existing relations and trends among various unstructured documents [18] by methods, such as keyword mapping or clustering of networks with similar content [19]. Keyword co-occurrence networks as part of bibliometric networks based on the context of citations [20] also enable the identification of differences and similarities of knowledge structures and sub-structures in health informatics. These findings will enable researchers to better understand the current state of health informatics, prevailing subjects, and future lines of research. This work is ultimately intended to extend the theoretical development and clarify the conceptual background of health informatics.
This research mapped the scientific networks, uncovered the prominent and hidden patterns in scientific networks, tracked the flow and burst of scientific subjects, and discovered what effects they have had on the scientific growth of health informatics. The originality of this study is related to its methodology and the timeframe used. We studied the evolution of health informatics during the past 44 years from 1974 to 2018 by applying three research methods: text mining, scientometric analysis, and social network analysis. These methods have not been previously used in studying the knowledge evolution of health informatics field. The scope of our timeframe (44 years) helped us analyze bursts and interactions between keywords, between countries, between authors, and between scientific subject categories since the first works were published in 1974 to provide a broader mapping than those used in previous qualitative studies.
This study aimed to illuminate the knowledge structure of health informatics by (1) reviewing a large number of publications (more than 30 thousand documents); (2) identifying hidden patterns during the last 44 years and visualizing them; (3) identifying the emerging scientific and technological trends since 1974, identifying the relations of keywords, authors, countries and scientific subject categories and their bursts; (4) identifying key studies and visualizing their relations; and (5) suggesting future lines of research.

Data Set Extraction and Filtration
This work was a quantitative study of health informatics science based on text mining and scientometric analysis of articles. This study analyzed and mapped four networks: co-occurrence of keywords, country co-citation network, co-occurrence of subject categories, and author co-citation network.
As the first step, we collected data from the Web of Science database by searching papers that included 'health informatics' in their subject. We limited our search to this keyword only and did not include other keywords, such as medical informatics or nursing informatics. The result of this search in August 2018 was 30,115 articles published during the period from 1974 to 2018. Regarding our inclusion and exclusion criteria, we included all papers from all disciplines and subjects and did not apply any specific exclusion criteria regarding time and discipline. Our inclusion criteria were to include all papers from all of the Web of Science subject categories and all document types. This enabled us to collect all relevant information of all documents on health informatics. However, we excluded non-English papers. Indexes of the search were SCI-Expanded, SSCI, A&HCI, CPCI-S, CPCI-SSH, ESCI, CCR-Expanded, and IC. Contents of records that were saved were full record and cited references. In sum, the population of this research was all scientific documents that included 'health informatics' as their subject and were indexed from 1974 until the end of August 2018 in the Web of Science Core Collection. To increase the quality of the data, we applied some pre-processing. Some examples of pre-processing steps were the removal of duplicates (126 duplicates were deleted), and stop words, tokening, and stemming. Stop words included the most common words like 'and' , 'if ' , and so forth. Tokenization included converting a sequence of characters into a sequence of tokens, and stemming was conducted to reduce inflected words to their root form. For instance, the words 'treatment' , 'treats' , and 'treated' were reduced to 'treat' .

Analytical Tools
We used multiple types of software for analysis and visualization as shown in Table 1. We used CiteSpace 5.3 and VOSviewer as analysis tools, and we used Gephi 0.9.2 and VOSviewer to visualize the networks.

1) Overlay visualization of keywords by VOSviewer
To conduct this analysis, 60,926 keywords were identified. We set the minimum number of occurrences to 40. Around 698 keywords met the threshold. 2) Word co-occurrence analysis and burst analysis by CiteSpace and visualization by Gephi To analyze word co-occurrence in CiteSpace, we set the number of years per slice to 5. We then selected the top 30% of most frequently occurring items from each slice. The common practice among previous studies was to select between the top 50% to 20% of the items. As a result, the CiteSpace software chose the 30 most cited or most frequently occurring items from each slice to construct the networks. We used the same criteria for all our analysis on CiteSpace (i.e., co-occurring subject categories, co-author citation network, word co-occurrence analysis). We used the Gephi software to visualize the network. The Gephi software identified 131 nodes and 466 edges. Our partitioning parameter was modularity. Modularity is used to identify clusters. Modularity results in grouping of nodes that are far more strongly connected and it gives insights into the strength of networks [21]. To calculate the modularity score, we chose the randomization option to produce a better decomposition, and we marked the 'use edge weights' option as well. We set the resolution to 0.7. The goal was to get more but smaller communities. The modularity result was 0.526, which is relatively average, and indicates reasonable relationships within the same clusters and reasonable relationships across the clusters. The number of detected clusters was 24. The graph layout was Force Atlas.

3) Country co-citation network by VOSviewer
To analyze and map the country co-citation network, we set the minimum number of documents of a country to 5. Out of the 157 countries, 103 countries met the threshold.

4) Co-occurring subject category by CiteSpapce and Gephi
After analyzing the network on the CiteSpace based on the criteria mentioned above, we used Gephi to visualize this network. We used the modularity parameter; the score was 0.622, and 11 clusters were identified. We used the Force Atlas as our layout.

5) Author co-citation network by CiteSpace and Gephi
We followed the same procedure on the CiteSpace software to analyze the network. The modularity score of the network was 0.754, and 81 clusters were identified. We used the Force Atlas layout for this as well.

Word Co-occurrence Analysis
From the analysis of co-occurring keywords, six clusters and nine research themes were derived, which are depicted in Table 2 and Figure 1. The research themes were identified by the authors. Cluster 1 (red color) is associated with studies on the utilization of computer science in the healthcare industry as well as the impact of health informatics on patient safety and the quality of healthcare. These studies claim that health information technology improves patient safety by reducing medication errors, reducing adverse drug reactions, and improving compliance with practice guideline (e.g., [22,23]). HIE also improves patient safety by measures such as improving medication information processing or improving laboratory information processing (e.g., [24]).
Cluster 2 (pink color) is associated with decision support, knowledge representation, and management in medicine. Studies claim that using these systems will improve clinical practice (e.g., [25]), improve the practice of evidence-based medicine (e.g., [26]), and reduce errors in medicine (e.g., [27]).
Cluster 3 (light blue color) is associated with the professional behavior change (e.g., [28]) and computer-based guideline implementation systems (e.g., [29]). Cluster 4 (dark blue color) is associated with the quality of health information on the internet (e.g., [30,31]). Cluster 5 (green color) is associated with health informatics education, nursing informatics [32,33], and clinical history taking [34]. The studies on health informatics education deal with the impact of health informatics on curriculum, education, and training of health care professionals as well as healthcare information systems research and development (e.g., [32,33]).
Cluster 6 (yellow color) is associated with the telehealth innovations in health education and healthcare.

Burst Analysis
The analysis of subject categories with the strongest citation bursts (Table 3) by using Kleinberg's burst detection algorithm shows the emergent research front concepts. Before 1991, no burst terms were identified. This analysis shows that in 1991, hospital information system and health informatics were burst terms, while in 1992, health information system, telematics, and primary healthcare were burst terms. Three years after Berners-Lee posted a short summary of the World Wide Web (WWW) project on the alt.hypertext newsgroup in 1991, the WWW became one of the burst terms of health informatics in 1994. In 1994, other terms, such as health informatics, computer-based patient record, and expert systems, became burst subjects. In 1995, some of the burst subjects were patient education and public health. In 1997, electronic patient record became a popular keyword. In 1999, subjects such as information retrieval, recommendation, medical information, health information, medical records system, security, confidentiality, patient record, and preventive care were burst terms. In 2002, bioinformatics and information system became burst subjects. In 2004, biohealth informatics received a burst. In 2007, the concept of e-health and in 2009, clinical decision support system became burst terms. No burst terms were identified between 2010 and 2015. In 2015, mobile health, big data, telehealth, prediction, machine learning, algorithm, social media, and mobile health became popular.
Overlay visualization of keywords ( Figure 2) with a minimum occurrence of 40 shows that since 2016, scholarly subjects that are of great interest in the health informatics field are precision medicine, big data, deep learning, machine learning, patient engagement, patient portals, engagement, m-health, social media, mobile applications, and the Internet of Things.

Countries with the Highest Numbers of Citations and Citation Link Strength
Visualization of countries with the highest numbers of cita-

Co-occurring Subject Category
The analysis of co-occurring subject categories (based on the Web of Science category) resulted in the clusters depicted in Figure 4. The figure shows that, in general, some of the categories of the blue cluster, including health informatics, healthcare science services, computer science, interdisciplinary applications, and information systems, have the highest contribution and interaction in the field of health informatics. Some of the categories of the purple cluster that have more interaction with the blue cluster are engineering, computer science, electrical engineering, biomedical engineering, and artificial intelligence. At the top of the network, there is a green cluster that has some interaction with the category of interdisciplinary applications in the blue cluster. Some of the categories of the green clusters are neurosciences, cell biology, biochemical research methods, biochemistry, and genetics.

Author Co-citation Network
The author co-citation analysis is visualized in Figure 5. Some of the key research strands of the highly cited authors are the following:

IV. Discussion
Based on the papers included in the WoS database service, it was clarified that health informatics is an information engineering field that is applied to healthcare [35]. The study showed that health informatics is applied to various subject categories, such as nursing, public health, biomedical research, and occupational therapy. The co-occurrence of the keywords showed that the overall goal is to improve the effectiveness of care delivery to patients [36]. This study shows that research on health informatics is not only concerned with engineering aspects but also with non-engineering sides of the health informatics. Issues such as the adoption of medical professionals of health informatics and behavioral changes are also key research themes [37]. Moreover, on the human side, as the burst analysis table  (Table 3) shows, concerns such as safety, security, surveillance, and privacy have been of great importance since the early stages of the development of health informatics [38]. The study showed that, since 2016, health informatics have entered a new era, which is predictive, preventative, personalized, and participatory. Health informatics has entered an era in which greater patient engagement with the support of information technologies is incorporated to improve health outcomes [39]. The connection of providers with patients is facilitated by the emerging health technologies, such as patient portals, social media, social health communities, wearables, self-tracking sensors and so forth [40]. On the other hand, patients will have access to consolidated medication management. With the emergence of 'quantified self ' and patients' awareness of their genetic profile, one possible strand of research can be on the changes on electronic health records and new forms of patient's engagement [41]. One other possible future research strand is related to patient-generated health data, increasing the literacy of patients on social and self-tracking tools, and on top of that the ethical issues of biometric and patient generated data. Another important strand of research is precision or personalized medicine to understand how a person's genetics, environment, and lifestyle can assist physicians to best treat and prevent diseases. One future strand of research could be the role of deep learning, new machine learning algorithms, and advanced big data analytics on precision medicine. On the other hand, with the advent of cutting-edge technologies, it is also necessary to conduct new studies on technology adoption and behavioral changes to improve healthcare management. While e-learning is studied highly in the literature, it is also necessary to study the effectiveness of mobile learning and peer-to-peer learning on patient outcomes.
This study also showed that the research strands of highly cited authors are medicine and information science; and the United States has the highest number of citations and the highest citation link strength compared to the rest of the world. The study also showed that diagnosis systems and preventive care were early scholarly subjects. However, most of the diagnosis systems have been for recognition; future studies could focus on early and preventive diagnosis systems with the aid of big data and machine learning methods. Other important subjects that have not been studied enough and could be of future research interest are open-source soft-


Interactive Deep Colorization With Simultaneous Global and Local Inputs

Colorization methods using deep neural networks have become a recent trend. However, most of them do not allow user inputs, or only allow limited user inputs (only global inputs or only local inputs), to control the output colorful images. The possible reason is that it's difficult to differentiate the influence of different kind of user inputs in network training. To solve this problem, we present a novel deep colorization method, which allows simultaneous global and local inputs to better control the output colorized images. The key step is to design an appropriate loss function that can differentiate the influence of input data, global inputs and local inputs. With this design, our method accepts no inputs, or global inputs, or local inputs, or both global and local inputs, which is not supported in previous deep colorization methods. In addition, we propose a global color theme recommendation system to help users determine global inputs. Experimental results shows that our methods can better control the colorized images and generate state-of-art results.

INTRODUCTION
Image colorization refers to the technique that adds colors to monochrome images or videos. Generally speaking, colorization is a ill-posed problem, it does not have a unique solution. To get satisfactory colorized results, two categories of method have been proposed: user-guided edit propagation and datadriven automatic colorization.
The user-guided edit propagation methods [1,2,3,4,5,6,7,8] require the user to draw colored strokes and propagate the colors across the image by solving a global optimization problem. These methods can achieve impressive colorized images, but often require a very large number of scribbles for images with complex textures. This is because each different color region must be explicitly marked by a different colored stroke, even regions with obvious semantic hints, such as a blue sky or green trees, need to be specified by the user.
To address this problem, early data-driven colorization methods proposed to automatically colorize a grayscale image by learning the color hints from one or several exemplar color images with similar semantics [9,10,11,12,13]. Unfortunately, it may be time-consuming or hard to find a suitable exemplar image sometimes. With the popularity of deep learning, recent data-driven colorization methods using deep neural networks have become a recent trend [14,15,16]. Using a large number of grayscale and color image pairs, the deep colorization methods learn a parametric mappings for fully automatic colorization. These methods can generate plausible colorful images in most of time. However, since an semantic region can have multiple choices of colors, the results can contain colors or styles which users do not expected. For example, users may want a green mountain in spring, but get a yellow mountain in autumn.
The excellent recent work, user-guide deep colorization, by Zhang et al [17] combines the advantages of user-guided and data-driven methods. It provide better user controls by taking either global inputs or local inputs in the deep network training. The color of the colorized image can be controlled by a global color histogram, or a few local color points. The user-guide deep colorization [17] can generate plausible colorized images according to the users expectation with only a few inputs. However, The user-guide deep colorization method does not support simultaneous global inputs and local inputs, it can only allow one kind of inputs at one time.
We argue that supporting simultaneous global inputs and local inputs can provide better control on the output images. The ideal case is that a user can control the overall color style of the image with global color inputs, and meanwhile assign local colors to certain regions with local inputs.
Supporting multiple kinds of inputs simultaneously in deep networks is not straightforward. It's difficult to differentiate the influence of different kind of user inputs in network training. An example is shown in Figure 6, the influence of local inputs is lost if the loss function is not designed appropriately. To solve this problem, we present a novel deep colorization method, which allows simultaneous global and local inputs to better control the output colorized images. Our basic model neural network is a U-net network including three parts, the feature extraction part, the fusion part, and the reconstruction part. The feature extraction part, consisting of 3 convolution layers and 3 pooling layers, extracts features from the gray-scale images and the local input images. The fusion part fuse the global inputs (color themes) and the extracted feature maps. The reconstruction part, consisting 3 deconvolution layers, then reconstructs the two chrominance channels (Lab space). The key step is to design an appropriate loss function that can differentiate the influence of input data, global inputs and local inputs. With this design, our method accepts no inputs, or global inputs, or local inputs, or both global and local inputs, which is not supported in previous deep colorization methods.
Different with previous works, which use exemplar images or histograms as global inputs, we choose to use the color theme, a template of colors possibly associated verbal description [18], as the global input. It is more easier for users to assign a color theme than to choose exemplar images or histograms. In our model, we accept color themes consist of 3 -7 colors. More colors are possible but difficult to use. To further save work time for users, we propose a global color theme recommendation system to help users determine the global input. Inspired by [18], we build a mapping between grayscale texture features and color histograms from a large image data base, and use the mapping to predict the color themes. Also, the user can define its own color theme.
Experimental results show that the color images generated by our method look real and natural, and the detail of the image are preserved well. Using different color themes and local inputs, we can output different style images (See Figure 1). To summarize, our methods can better control the colorized images and generate state-of-art results.
Our contributions in this paper are as follows: • We propose a novel deep colorization methods that supports no inputs, or global inputs, or local inputs, or both global and local inputs.
• Our methods enables global inputs using color themes of variable color numbers.
• We propose a color theme recommendation system which can suggest the color themes for users.

RELATED WORK
Image colorization is always a hot topic in the filed of image process. After the advent of deep convolution neural network, it is devoted to train images on network to achieve the goal of great colorization. At present, there are two main ways of colorization. One is user-guided edit propagation colorization, and the other is automatic colorization.
About user-guided coloring method, most of them are devoted to add scribbles on gray images. In the early paper, like [1], it is necessary for users to draw some desired color graffitis on certain regions, through the color propagating among pixels which are similar in intension to achieve the effect of colorization. [3] and [19] worked on the issue of edge detection and enhancement, which can improve the quality of color transmission. In [5], users select the area of interest and then specify some color points within the area. [20] and [21] are optimized for edit propagation and shading. Its also popular to use color theme, such as using data driven color theme enhancements in [18], which is used to change the color style of color images through different color themes. We take this method and using color theme as global input to color gray images. [22] is a popular method to use the palette to change the corresponding color on color image as same as [18]. In [23], it is proposed to automatically learn the degree of similarity between user strokes and input images to propagate graffiti color. [24] uses generative adversarial networks to achieve transforming sketch to real images, which also supports to add user strokes. In [17], by adding global features and local features on deep convolution neural network to train a colorization model, separately, users can specify the color at any location to render the color, or use a reference image to transfer global color. Our method can achieve colorization using both of them simultaneously. There are online applications developed mainly for line-drawing colorization, like [25] and [26].
For more convenient to color grayscale images globally, there are some method using reference images as global inputs early. [27] transfers colors between images with a simple algorithm. [9] and [10] transmit the color emotion of reference images to grayscale images by matching the brightness and texture information between images. [11] needs users to provide the semantic information of the foreground image. The system searches for and downloads pictures with the same foreground semantics on the network, and use them to color the foreground and the background. [12] can quickly transfer the color of a reference image to the target image by using a fast cascade feature matching scheme to exploit multiple image features.
But using reference images which is similar with target images to render the color needs a large image database. Users also have to spend a lot of time to choose satisfactory reference image. So, the research of automatic colorization gets more and more attention. [28] uses algorithm to reserve grayscale images or videos. In this method, some representative color points are reserved and used to restore these medias. [29] uses machine learning tools to extract as much information as possible from the color sample dataset and then estimates the probability color distribution for every pixel. [30] matches grayscale image with reference images downloaded automatically from web to transfer color. [31] train the objective function based on image features on the LEARCH framework and achieve the coloring effect by minimizing the objective function. [14] introduces the concept of feature descriptor, and take the feature descriptors extracted from the grayscale image as input, and finally output the corresponding color values of the UV channel. In [15], a convolutional neural network (CNN) which accepts black and white images as input is designed and constructed, and a statistical learning driven method is used to solve the problem of grayscale colorization. The end-to-end network in [32] adds global feature and uses classification label to optimize output results. In addition, [33], [34], [16], [35], [36] also achieve the effect in coloring gray image without users' intervention. [37] uses conditional adversarial networks and [25] proposes a setup utilizing two networks in tandem to achieve colorization, which both are worth thinking about.
Neural network also has shown us some surprising results in other fields of image processing. [38] reaches the target of automatic style transferring between images. As well as large-scale image recognition [39], automatic photo adjustment [40], sketch simplification [41], context encoders [42] and cartoon colorization [43]. Accuracy and practicability of neural network is the reason that we chose neural network to research and construct satisfactory colorization method.

OUR METHOD
Our deep network model takes in a gray-scale image (L channel in CIE Lab space), a global input and a local input and output the corresponding ab channels. The structure of the model is shown in Figure 2. The global input includes a color theme of several colors and a mask indicating the number of colors (all black means no global inputs). The local input includes a image consisting of the assigned local color points and a mask indicating the position of color points (all black means no local inputs). Therefore, by setting the masks of global inputs and local inputs, our model accepts four combination of user inputs: no inputs, global inputs, local inputs, global inputs and local inputs.

Problem Formulation
The inputs of our model are a grayscale image X ∈ R H×W ×1 , a global input tensor U g ∈ R 1×K×3 , and a local input tensor U l ∈ R H×W ×3 , where H, W is the width and height of the input image, and K is the number of colors in the color theme. The output is a tensor O ∈ R H×W ×2 . We aims to train a convolution neural network (CNN), denoted by F(X, U g , U l ; θ), to approximate the mapping between gray and color, under the constraint of user inputs. Therefore, the colorization problem can be formulated as where D denotes the training data set, L denotes the loss function, and Y ∈ R H×W ×2 is ground truth image color. More details of the loss function will be presented in Section 3.4. It should be noticed that we train a single model to handle simultaneous global input and local input. Moreover, our loss function is explicitly related to the global input and the local input. These are the two key points which differ our method from [17].

User inputs
Global inputs: To control the overall color style of the output, we design the global input as a color theme with K colors plus a one-channel mask indicating the number of colors.   Figure 2. To prepare training data for the global input, we use the K-mean clustering algorithm to find the K representative colors for each color image in the data set. The K representative colors form the global color theme. To enable color theme with variable number of colors, K randomly varies in the interval [3,7]. Color theme with more colors are also supported by modifying the interval, but may be harder to use. The ab channels of each color theme U c g ∈ R 1×K×2 and its mask M g ∈ R 1×K×1 forms the global input U g = {U c g , M g } ∈ R 1×K×3 For each color image, we generate a K-color map by decoding the color image with its representative colors. Figure 3 shows an example. The ab channels of K-color map, I ∈ R H×W ×2 is used to calculate the loss function which measures the similarity between the color theme and the output image. More details of loss functions will be described in Section 3.4.1.
Local inputs: Similar to [17], we prepare the training data for the local inputs colors U c l ∈ R H×W ×2 by randomly selecting some points from the ab channels of every color image. If there are no samples in certain position, the ab values are set to zeros. A one-channel mask M l ∈ R H×W ×1 is also generated to indicate the positions of the local inputs. "1" in certain position means there is an input there, while "0" means there are no inputs in the position. Finally, the local

NetWork model
As shown in Figure 2, our network model uses a U-Net structure [44], which is already shown to works well in colorization [37,17]. It is mainly composed of four parts: the feature extraction module, the global input module, the fusion module and the reconstruction module.

Feature extraction module
The feature extraction module corresponds to the blue part in Figure 2. The feature extraction module takes in a grayscale image X ∈ R H×W ×1 and a local input U l ∈ R H×W ×3 . Initially, X and U l are respectively convolved to a tensor of size H × W × 32. These two tensors are merged to a single tensor of size H × W × 32 by linear interpolation, whose weight is treated as a parameter of the network to be trained. The merged tensor is then feeded to the following convolution blocks.
In convolution block 2 to 7, all the convolution kernels are 3 × 3. Through every layer, only the tensor sizes of input are halved spatially with the stride 2 × 2, or only the tensor dimensions are doubled with the stride 1 × 1. After conv7, the tensor is H 8 × W 8 × 512. We use convolution layers with stride 2 × 2 instead of using max-pooling layers to reduce the tensor sizes, which is helpful in increasing the spatial support of each layer [32]. Conv8 and conv9 process the tensor further with convolution layers to reduce the 512-channel tensors to 256-channel tensor. The output of feature extraction module is a tensor of size H 8 × W 8 × 256, which will be fused with the global input.

Global input module
The global input module, the pink part in Figure 2, is one of the characteristics of our method. It takes in the global inputs U g ∈ R 1×K×3 , which consists of the ab channels of the color theme U c g and the corresponding mask M g . To to unify its size with the size of the fusion module, we reshape it to a tensor of size 1 × 1 × 3K, which is then processed by three fully connected layers. This operation is similar to the global feature network in [32], but we use it to handle user global inputs instead of the features exacted from the input gray-scale image.

Fusion module
The fusion layer, shown with yellow color in Figure 2, is an important step to fuse the global input and the extracted features. Similar operations are also used in [32,17]. The global input is fused with the output of the feature extraction module by linear interpolation, whose weight is also treated as a parameter in the network to be trained. Finally the fusion module outputs a feature tensor of size H 8 × W 8 × 256 .

Reconstruction module
After the fusion module, the feature tensor is processed by a set of convolution layers and upsampling layers for ab channel reconstruction, shown with purple color in Figure 2. Convolution layers reduce tensor dimensions by half and upsampling layers doubles the wide and height of the tensor. Conv17 is the last convolution layer with a Sigmoid transfer function, followed by an upsampling layer. The advantage of Sigmoid is that the output range is limited. The output tensor is of size H × W × 2, which combines the input gray-scale image X ∈ R H×W ×1 to generate the final output O ∈ R H×W ×3 .

Loss function
Designing an appropriate loss function is the most important part of our work. It is not straightforward to design a loss function which can differentiate the impact of different inputs. We interpret the four combinations of inputs as follows: 1 No inputs means the user wants the colors to be assigned with the "experiences" of training data.
2 Only the global input means the user wants the colors to be conditionally assigned with the "experience" of training data.
3 Only the local input means the user wants to assign his/her preferred color to certain region, but let the colors of other regions be assigned with the "experiences" of training data. 4 Both inputs mean that the user wants to assign his/her preferred color to certain region, but let the colors of other regions be conditionally assigned with the "experiences" of training data.

Loss function for no input or global input
For no inputs or only the global input, the straightforward choice is to measure the differences between ab channels of the output image O ∈ R H×W ×2 and the ground truth color Y ∈ R H×W ×2 . There are many kind of loss functions as discussed in [17]. We choose the Huber loss as it produces relative high saturation effect. The Huber loss is given by where δ is the parameter of the Huber loss. The value of δ will slightly affect the results as shown in Figure 4. We set δ = 0.5, as this value works well in our experiments. However, in some cases, the impact of the global input (color theme ) is not so obvious. Therefore, we add an additional part which calculates the Huber loss of the output image and the K-color map I, which is defined by decoding the ground truth with the color theme (see Section 3.2 and Figure 3). Therefore, the loss function for this part is given by where α 1 and α 2 are two parameters to balance the influence of two parts. The influence of the two parameters is shown in Figure 5. As we expect that output images not only keep natural color but also reflect colors of the color theme, we set α 1 = 0.7 and α 2 = 0.3, with which the output can preserve the characteristics of the color theme with almost no color overflows, as showed in Figure 5. When there is no color theme input, the K-color map image will be replaced by the ground truth. In this case, L g = L y .

Loss function for local input
The local input in our method is U l = {U c l , M l } ∈ R H×W ×3 , which is combined by the ab channels of the user input and the corresponding mask, as introduced in Section 3.2. If a user provides a local input in a certain position, we think the user wants the preferred colors instead of colors assigned by the "experience" of training data. In there is only the local input, the Hubber loss L H (O, Y ) defined in Equation (2) is already shown to works well [17], since the differences of the output image and the local inputs at positions with local inputs are already contained in Equation (3).
However, if there are simultaneous global input and local input, the case becomes more complex. The Hubber loss may seems to work sometimes, when the local inputs are consistent or close with the global color theme, as shown in Figure 6(f). For example, the input brown color in the back and swing is close to one color in the color theme (Figure 6(a)). In this case, the local input can enhance the colors of local regions compared to results with only global inputs (Figure 6(b)). But there exist some colors which overflow from the bird body to the background (TODO, indicate in the figure).
The color overflow can be removed by adding the gradient of the output image and the ground truth in the loss function, given by where MSE denotes the mean squared error, Sobel denotes the Sobel gradient operator. As shown in Figure 6(g), the color overflows are removed at the cost of removing the local color impacts on the head and abdomen. When the local input colors are quite different with the colors in the global color theme, the phenomenon is more obvious, as shown in Figure 6(j)(k). Using only the loss L g cannot guarantee the local color impact, and results in more obvious color overflows (Figure 6(j)). The gradient loss L s almost remove the color overflows, but also remove the impact of local color inputs (Figure 6(k)). Since we expect the colorized image preserves not only the global input color theme, but also the local input colors, without few color overflows. We emphasize the impact the local inputs by adding the loss of the output image and the local inputs at positions with local inputs, given by where MSE denotes the mean squared error, M l is the mask of the local input. With L p , we can see from Figure 6(h)(i) that the aforementioned problems are solved. To summarize, our final loss function L is L(F(X, U g , U l ; θ), U g , U l , Y )

Color theme recommendation system
To help user determine the possible color theme of a grayscale image, we propose a color theme recommending system. Our system is inspired by the data driven method in [18], which consists of an off-line process to builds a texturecolor distribution mapping and an on-line process to lookup the color distribution for the segmented input color image. Different with the off-line process in [18], which built its texture-color distribution mapping on color images to constrain the color distribution of textures, we build the mapping on the gray-color image pairs to predict the colors themes. The framework of color theme commended system is illustrated in Figure 7. It consists of an off-line phase(lower row) and an online phase(upper row). The mapping of texture-color distribution is obtained in the off-line stage. The color themes of gray-scale image that input by users can be predicted by computing the main texture of gray-scale image in online phase.
During the off-line phase, we start with building an image database, by randomly choosing thousands of colorful images from Imagenet, and then convert the colorful images to grayscale images to form a grayscale image database. Next we adapt the graph-based method [11] to segment each grayscale image into segmentations. We use the mean and standard deviation of the pixelwise Gabor filter (scale=4, rota-tion=6) responses within every segmentation, resulting a 24 dimensional texture descriptor for every segmentation. We then build a texture library by clustering all the texture features into 120 textures with k-means algorithm. For each cluster, we accumulate a 2D color distribution with 100 bins by quantizing the corresponding ab channels of the segmentations. Hence, we build a texture-color distribution mapping.
In the online-phase, we segment the input gray-scale image, extract the texture feature vectors, find the closet clusters in the texture library, and lookup the corresponding color distributions. We choose the most frequent color of the K largest segmentations (with more number of pixels) as our recommended color theme. We display some images colored by color themes which are generated by this system in Figure 8.

EXPERIMENT
We implemented the network model on the NIVIDIA GTX1080Ti GPU, using Google's TensorFlow architecture. We use the Place dataset to train our model. Figure 9 shows some colored results of this dataset. The K-color map images, color themes, local inputs and their corresponding masks are generated on Matlab. Then, we randomly adjust the order of training sets as input of our model. We train the model using a batch size of 50 for 40,000 iterations, which takes about two days on the GPU. To enable the four combinations of input: no inputs, only global input, only local input, and both global and local input, we randomly select one kind of combination in the training. The ratio to select each combination is equal to 25%.
In the previous works, the normalization of input has been shown to speed up the learning speed. In Lab space, the range of L is [0 − 100], and the range of ab channel is [127, −128]. In order to accelerate the learning speed, we normalize the range of L and ab to [0,1] as L = L/100 ab = (ab + 128)/255 (7)

Automatic Colorization
In these section, we compare our automatic colorization method with several state-of-art methods, including [17], which train coloring with no inputs or local inputs , [32] and [33], which focus on automatic colorization. The colorized results of these four methods are shown in Figure 10.
As shown in the figure, all the automatic methods can colorize the gray-scale images. But due to the different training sets, the colorized images may have different color styles. In terms of numerical evaluation, our method and [32] get better PSNRs without compared with the other two methods. All PSNRs of our automatic results are higher than [17]. Note that, although automatic colorization is only our by-product, it still has good performances.

Interactive Colorization
In this section, we compare our interactive colorization method with the state-of-art user guide method [17], and the classical optimization method [1]. The results are shown in Figure 11. The Butterfly images in the first three columns show results colorized by different local input points. With only local inputs, the results of our method are similar to those of the user-guide method [17]. In contrast, the results of [1] fail to diffuse the color to the whole butterfly. For the next five images, our method uses global input (color theme) and local input together. The other two methods use local inputs only. For images with less clear boundaries, like the fourth column, our method can faithfully assign green to the trees, where the other methods cause color overflows to the background. Using our method, users can control the colorized image by choosing a color theme, and add a small number of local input points to certain regions, while other methods maybe need to add more local input points to produce the same results. For example, in the fourth column, users maybe need to add several extra color points on the house and sky to avoid color overflow and assign certain colors for the user guide method [17]. [1] needs to assign more color points for everywhere in every image. Otherwise, the whole image is going to be one color as shown in the figure. Besides, our method also shows good performance on other type images, like human and natural pictures showed in the last three columns of the Figure 11.  7. The overall pipeline of color theme commended system framework. The (c) is a texture library which stores the texture center of 120 textures. The image in (d) is a mapping between each pixel and texture, so the value of the pixel is the texture type of the pixel belongs. In the same way, (f) can be converted into (g) by calculating the nearest texture center (from (c)) of each segmentation. So the main texture in (f) can be counted. Combine with the main texture, the desired color theme into consideration to generate final recolored image in (h) can be found from (e).  In general, our method provides a more convenient and efficient tool for the users by supporting four combinations of input in a single network: no inputs, only global inputs, only local inputs, and simultaneous global input and local input.

Numerical Comparisons
Besides visual comparisons, we also evaluate the PSNRs of the related methods. We randomly choose another 20 images and colorize them using the automatic methods and interactive methods mentioned above. We use original images as global inputs for [17], and use color themes extracted from original images as global inputs for our method. Local inputs are 3 to 20 colorful points randomly extracted from original images. Then we calculate the average PSNR of them as showed in Table 1. Through Table 1, we can see that results of automatic methods have similar performances, which have almost identical PSNRs. When we add global inputs, our method works much better and gets a higher PSNR.
For method with global inputs, the user guide method [17] shows good performance and get the highest PSNR. Our method also works well. But the PSNR of method is a bit lower than the user guide method [17]. This is possibly because we train our network to support four combination of inputs , the effect of local inputs is not optimized in our method.
For method with local inputs, the user guide method [17] shows good performance and get the highest PSNR. Our method also works well. But the PSNR of method is a bit lower than the user guide method [17]. This is possibly also because we train our network to support four combination of inputs , the effect of local inputs is not optimized in our method. The optimization method [1] gets a lower PSNR, probably because it was designed for strokes not point inputs. When both global inputs and local inputs are used, our method get a higher PSNR, which is very close to the highest in Table 1.

Colorization for past
We test our method on some historic images, which are whiteand-black images. Although these images are different from our dataset, as that they may have rough edges and their picture quality may be unclear, they also can be colorized by our method. Compared with our method, other methods, which use local input only, maybe need a lot of local input points to get the same results. Some examples of coloring historic images are shown in Figure 12.

Computation time
We test computation time of two image sizes on the NIVIDIA GTX1080Ti GPU. We take the average time of 100 computations on 50 images to get a reliable testing value, as shown in Table 2. The data shows that our method is comparable to real-time colorization.

Limitations and discussion
Our model is trained by the Place dataset, which can not include everything in this world. If an gray-scale image contains certain contents, which is not learned from the data, the colors may be not assigned appropriately. In this case, the user may need to give more inputs to the system. For some images with unclear edges, the network may produce unexpected results, like color overflow on background. Users can add another local input points in the background region. In addition, When users add local input points with unusual colors, they may have to add more points to achieve the desired effect or directly change to use an appropriate color theme.  [1] Fig. 11. Comparisons of interactive colorization methods. We only add local inputs on the first grayscale image and change the color of some points on the Butterfly wings, without changing their position. For the last five images, we add color theme and local inputs simultaneously, while other methods add local inputs only.

CONCLUSION AND FUTURE WORK
In this paper, we propose a novel interactive deep colorization method. By designing a suitable loss function, our method allows four combinations of input in a single network model, including no inputs, only global input, only local input, and both global input and local input. Furthermore, our global input is a color theme with variable number of colors, which is more easy and straightforward for user to use. Together with recommendation system, our method can reduce time and cost of coloring images for users. In addition, we evaluate our method on images of outdoor, human picture, past pictures, and show that it can produce satisfactory results for all them. In the future, we would like to extend our method to gray-scale video sequence.