\section{Introduction}

Repeated measurements data consists of observations from several experimental units, subjects, or groups under different conditions.
This grouping or clustering of the individual responses into experimental units
typically introduces dependencies:
the different units are assumed to be independent,
but there may be heterogeneity between units
and correlation within units.

Mixed-effects models
provide a powerful and flexible tool to analyze grouped data by incorporating fixed and random effects. Fixed effects are associated with the entire population, and random effects are associated with individual
groups
and model the heterogeneity across them and the dependence structure within them~\citep{Pinheiro2000}.
Linear mixed-effects models~\citep{Laird-Ware1982, Pinheiro2000, Verbeke-Molenberghs2009, Demidenko2004}
 impose a linear relationship between all covariates and the response.
Partially linear mixed-effects models~\citep{Zeger-Diggle1994}  extend the linear ones.

We consider the partially linear mixed-effects model
\begin{equation}\label{eq:initPLMM}
	\mathbf{Y}_i = \Xi\beta_0 + g(\mathbf{W}_i) + \mathbf{Z}_i\mathbf{b}_i+ \boldsymbol{\varepsilon}_i
\end{equation}
for groups $i\in\{1,\ldots,N\}$. There are $n_i$ observations per group $i$.
The unobserved random variable $\mathbf{b}_i$, called random effect, introduces
correlation
within its group $i$ because all $n_i$ observations within this group share $\mathbf{b}_i$.
We make the assumption generally made that both the random effect $\mathbf{b}_i$ and the error term $\boldsymbol{\varepsilon}_i$
follow a Gaussian distribution~\citep{Pinheiro2000}.
The matrices $\mathbf{Z}_i$ assigning the random effects to group-level observations are fixed.
The linear covariables $\Xi$ and the nonparametric and potentially high-dimensional covariables $\mathbf{W}_i$ are observed and random, and they may be dependent. Furthermore, the nonparametric covariables may contain nonlinear transformations and interaction terms of the linear ones.
Please see Assumption~\ref{assumpt:Distribution} in Section~\ref{sect:ModelAndDML} for further details.

Our aim is to estimate and make inference for the so-called fixed effect $\beta_0$ in~\eqref{eq:initPLMM} in the presence of a highly complex $g$ using general machine learning
algorithms.
The parametric component $\beta_0$
provides a simple summary of the covariate effects that are of main scientific interest. The nonparametric component $g$ enhances model flexibility  because time trends and further covariates with possibly nonlinear and interaction effects can be modeled nonparametrically.
\\

Repeated measurements, or longitudinal, data is omnipresent
in empirical research.
For example, assume we want to study the effect of a treatment over time. Observing the same subjects repeatedly presents three main advantages over having cross-sectional data.
First, subjects can serve as their own controls.
Second, the between-subject variability is explicitly modeled and can be excluded from the experimental error. This yields more efficient estimators of the relevant model parameters.
Third, data can be collected more reliably~\citep{Davies2002, Fitzmaurice2011}.
\\

Various approaches have been considered in the literature to estimate the nonparametric component
$g$ in~\eqref{eq:initPLMM}:
kernel methods~\citep{Hart-Wehrly1986, Zeger-Diggle1994, Taavoni2019b, Chen2017},
backfitting~\citep{Zeger-Diggle1994, Taavoni2019b},
spline methods~\citep{Rice-Silverman1991, Zhang2004, Qin-Zhu2007, Qin-Zhu2009, Li-Zhu2010, Kim2017, Aniley2019},
and local linear regression~\citep{Taavoni2019b, Liang2009}.

Our aim is to make inference for $\beta_0$ in the presence of potentially highly complex effects of $\mathbf{W}_i$ on $\Xi$ and $\mathbf{Y}_i$. First, we adjust $\Xi$ and $\mathbf{Y}_i$ for $\mathbf{W}_i$ by regressing $\mathbf{W}_i$ out of them using  machine learning algorithms. These machine learning algorithms may yield biased results, especially if regularization methods are used, like for instance with the lasso~\citep{Tibshirani1996}.
Second, we fit a linear mixed-effects model to these regression residuals to estimate $\beta_0$. Our estimator of $\beta_0$ converges at the parametric rate, follows a Gaussian distribution asymptotically, and is semiparametrically efficient.

We rely on the double machine learning
framework of~\citet{Chernozhukov2018} to estimate $\beta_0$ using general machine learning algorithms.
To the best of our knowledge, this is the first work to allow the nonparametric nuisance components of a
partially linear mixed-effects model to be estimated with arbitrary machine learners like random forests~\citep{Breiman2001} or the lasso~\citep{Tibshirani1996, Buehlmann2011}.
In contrast to the setting and proofs of~\citet{Chernozhukov2018}, we have dependent data and need to incorporate this accordingly.

\citet{Chernozhukov2018} introduce double machine learning to estimate a low-dimensional parameter in the presence of nonparametric nuisance components using machine learning
algorithms.
This estimator converges at the parametric rate and is asymptotically Gaussian due to Neyman orthogonality and sample splitting with cross-fitting.
We would like to remark that nonparametric components can be estimated without sample splitting and cross-fitting if the underlying function class satisfies some entropy conditions;
see for instance~\citet{Geer-Mammen1997}. However, these regularity conditions limit the complexity of the function class, and
machine learning algorithms usually do not satisfy them. Particularly, these conditions fail to hold if the dimension of the nonparametric variables increases with the sample size~\citep{Chernozhukov2018}.

\subsection{Additional Literature}

Expositions and overviews of mixed-effects modeling techniques can be found in~\citet{Pinheiro1994, Davidian-Giltinan1995, Vonesh-Chinchilli1997, Pinheiro2000, Davidian-Giltinan2003}.
\\

\citet{Zhang1998} consider partially linear
mixed-effects models and estimate the nonparametric component with natural cubic splines. They treat the smoothing parameter as an extra variance component that is jointly estimated with the other variance components of the model.
\citet{Masci2019} consider partially linear
mixed-effects models for unsupervised classification with discrete random effects.
\citet{Schelldorfer2011} consider  high-dimensional linear mixed-effects models where the number of fixed effects coefficients may be much larger than the overall sample size.
\citet{Taavoni2019a} employ a regularization approach in generalized
partially linear mixed-effects models using regression splines to approximate the nonparametric component.
\citet{gamm4} use penalized regression splines where the penalized components are treated as random effects.
\\

The unobserved random variables in the partially linear mixed-effects model~\eqref{eq:initPLMM}
are assumed to follow a Gaussian distribution.
\citet{Taavoni2021} introduce multivariate $t$
partially linear
mixed-effects models for longitudinal data. They consider $t$-distributed random effects to account for outliers in the data.
\citet[Chapter 4]{Fahrmeir2011} relax the assumption of Gaussian random effects in generalized linear mixed models. They consider nonparametric Dirichlet processes and Dirichlet process mixture priors for the random effects.
\citet[Chapter 3]{Ohinata2012} consider partially linear mixed-effects models and make no distributional assumptions for the random terms,
and the nonparametric component is estimated with kernel methods.
\citet{Lu2016} consider a partially linear mixed-effects model that is nonparametric in time and that features asymmetrically distributed errors  and missing data.

Furthermore, methods have been developed to analyze repeated measurements data that are robust to outliers.
\citet{Guoyou2008} consider robust estimating equations and estimate the nonparametric component with a regression spline.
\citet{Tang2015} consider median-based regression methods in a partially linear model with longitudinal data to account for highly skewed responses.
\citet{Lin2018} present an estimation technique in partially linear models for longitudinal data that is doubly robust in the sense that  it simultaneously accounts for missing responses and mismeasured covariates.
\\

It is prespecified in
the  partially linear mixed-effects model~\eqref{eq:initPLMM} which covariates are modeled with random effects.
Simultaneous variable selection for fixed effects variables and random effects has been developed by~\citet{Bondell2010, Ibrahim2011}. They use penalized likelihood approaches.
\citet{Li-Zhu2010} use a nonparametric test to test the existence of random effects in partially linear mixed-effects models.
\citet{Zhang2020} propose a variable selection procedure for the linear covariates of a generalized partially linear model with longitudinal data.
\\

\textit{Outline of the Paper.}
Section~\ref{sect:ModelAndDML} presents our double machine learning estimator of the linear coefficient in a partially linear mixed-effects model. Section~\ref{sect:numerical-experiments} presents our numerical results.
\\

\textit{Notation.}
We denote by $\indset{N}$ the set $\{1,2,\ldots,N\}$. We add the probability law as a subscript to the probability operator $\Prob$ and the expectation operator $\E$ whenever we want to emphasize the corresponding dependence.
We denote the $L^p(P)$ norm by $\normP{\cdot}{p}$ and the Euclidean or  operator
norm  by $\norm{\cdot}$, depending on the context.
We implicitly assume that given expectations and conditional expectations exist. We denote by $\stackrel{d}{\rightarrow}$ convergence in distribution.
The symbol $\independent$ denotes independence of random variables.
We denote by $\mathds{1}_{n}$ the $n\times n$ identity matrix and omit the subscript $n$ if we do not want to emphasize the dimension.
We denote the $d$-variate Gaussian distribution by $\mathcal{N}_d$.

\section{Model Formulation and the Double Machine Learning Estimator}\label{sect:ModelAndDML}

We consider repeated measurements data that is grouped according to experimental units or subjects.
This grouping structure introduces dependency in the data. The individual experimental units or groups are assumed to be independent, but there may be some
between-group heterogeneity and within-group correlation.
We consider the partially linear mixed-effects model
\begin{equation}\label{eq:PLMM}
	\mathbf{Y}_i = \Xi\beta_0 + g(\mathbf{W}_i) + \mathbf{Z}_i\mathbf{b}_i+ \boldsymbol{\varepsilon}_i, \quad i\in\indset{N}
\end{equation}
for groups $i$ as in~\eqref{eq:initPLMM} to model the
between-group heterogeneity and within-group correlation with random effects.
We have $n_i$ observations per group that are concatenated row-wise into $\mathbf{Y}_i\in\R^{n_i}$, $\Xi\in\R^{n_i\times d}$, and $\mathbf{W}_i\in\R^{n_i\times v}$.
The nonparametric  variables may be high-dimensional, but $d$ is fixed.
Both $\Xi$ and $\mathbf{W}_i$ are random.
The  $\Xi$ and  $\mathbf{W}_i$ belonging to the same group $i$ may be dependent.
For groups $i\neq j$, we assume $\Xi\independent\mathbf{X}_j$, $\mathbf{W}_i\independent\mathbf{W}_j$, and $\Xi\independent\mathbf{W}_j$.
We assume that $\mathbf{Z}_i\in\R^{n_i\times q}$ is fixed.
The random variable $\mathbf{b}_i\in\R^q$ denotes a group-specific vector of random regression coefficients that is assumed to follow a Gaussian distribution.
The dimension $q$ of the random effects model is fixed.
Also the error terms are assumed to follow a Gaussian distribution as is commonly used
in a mixed-effects models framework~\citep{Pinheiro2000}.
All groups $i$ share the common linear coefficient $\beta_0$ and the potentially complex function $g\colon\R^{v}\rightarrow\R$. The function $g$ is applied row-wise to $\mathbf{W}_i$, denoted by $g(\mathbf{W}_i)$.

We denote the total number of observations by $\NN_T:=\sum_{i=1}^{N}n_i$.
We assume that the numbers $n_i$ of within-group observations are uniformly upper bounded by $n_{\mathrm{max}}<\infty$.
Asymptotically, the number of groups, $N$,  goes to infinity.

Our distributional and independency assumptions are summarized as follows:

\begin{assumptions}\label{assumpt:Distribution} Consider the partially linear mixed-effects model~\eqref{eq:PLMM}. We assume that there is some $\sigma_0>0$ and some symmetric positive definite matrix
$\Gamma_{0}\in\R^{q\times q}$
such that the following conditions hold.
\begin{enumerate}[label={\theassumptions.\arabic*}]
	\item \label{assumpt:D1}
		The random effects $\mathbf{b}_1,\ldots,\mathbf{b}_{\NN}$ are independent and identically distributed\ $\mathcal{N}_q(\boldsymbol{0}, \Gamma_{0})$.
	\item \label{assumpt:D2}
		The error terms $\boldsymbol{\varepsilon}_1,\ldots,\boldsymbol{\varepsilon}_{\NN}$ are independent
		and follow a Gaussian distribution, $\boldsymbol{\varepsilon}_i\sim\mathcal{N}_{n_i}(\boldsymbol{0}, \sigma_0^2\mathds{1}_{n_i})$ for $i\in\indset{N}$, with the common variance component $\sigma_0^2$.
	\item \label{assumpt:D3}
		The variables $\mathbf{b}_1,\ldots,\mathbf{b}_{\NN}, \boldsymbol{\varepsilon}_1,\ldots,\boldsymbol{\varepsilon}_{\NN}$ are independent.
	\item \label{assumpt:D4-2}
		For all $i, j \in\indset{N}$, $i\neq j$, we have $(\mathbf{b}_i,\boldsymbol{\varepsilon}_i)\independent(\mathbf{W}_i, \Xi)$ and $(\mathbf{b}_i,\boldsymbol{\varepsilon}_i)\independent(\mathbf{W}_j, \mathbf{X}_j)$.
	\item \label{assumpt:D6}
	         For all $i, j \in\indset{N}$, $i\neq j$, we have $\Xi\independent\mathbf{X}_j$, $\mathbf{W}_i\independent\mathbf{W}_j$, and $\Xi\independent\mathbf{W}_j$.
\end{enumerate}
\end{assumptions}

We would like to remark that the distribution of the error terms $\boldsymbol{\varepsilon}_i$ in Assumption~\ref{assumpt:D2} can be generalized to $\boldsymbol{\varepsilon}_i\sim\mathcal{N}_{n_i}(\boldsymbol{0},\sigma_0^2\Lambda_i(\boldsymbol{\lambda}))$, where $\Lambda_i(\boldsymbol{\lambda})\in\R^{n_i\timesn_i}$ is a symmetric positive definite matrix parametrized by some finite-dimensional parameter vector $\boldsymbol{\lambda}$ that all groups have in common.
For the sake of notational simplicity, we restrict ourselves to Assumption~\ref{assumpt:D2}.

Moreover, we may consider stochastic random effects matrices $\mathbf{Z}_i$. Alternatively, the nonparametric variables $\mathbf{W}_i$ may be part of the random effects matrix. In this case,
we consider
the random effects matrix $\widetilde{\mathbf{Z}}_i = \zeta(\mathbf{Z}_i,\mathbf{W}_i)$ for some known function $\zeta$ in~\eqref{eq:PLMM} instead of $\mathbf{Z}_i$.
Please see Section~\ref{sec:nonfixedZi} in the appendix for further details.
For simplicity, we restrict ourselves to fixed random effects matrices $\mathbf{Z}_i$ that are disjoint from $\mathbf{W}_i$.
\\

The unknown parameters in our model are $\beta_0$, $\Gamma_{0}$, and $\sigma_0$.
Our aim is to estimate
$\beta_0$ and
make inference for it.
Although the variance parameters $\Gamma_{0}$ and $\sigma_0$ need to be estimated consistently to construct an estimator of $\beta_0$, it is not our goal to perform inference for them.

\subsection{The Double Machine Learning Fixed-Effects Estimator}

Subsequently, we describe our estimator of $\beta_0$ in~\eqref{eq:PLMM}.
To motivate our procedure, we first consider the population version with the residual terms
\begin{displaymath}
	\mathbf{R}_{\Xi} := \Xi - \E[\Xi | \mathbf{W}_i] \quad\textrm{and}\quad
	\mathbf{R}_{\Yi} := \mathbf{Y}_i - \E[\mathbf{Y}_i | \mathbf{W}_i] \quad\textrm{for}\quad
	i\in\indset{N}
\end{displaymath}
that adjust $\Xi$ and $\mathbf{Y}_i$ for $\mathbf{W}_i$. On this adjusted level, we have the linear mixed-effects model
\begin{equation}\label{eq:LMM}
	\mathbf{R}_{\Yi} = \mathbf{R}_{\Xi} \beta_0 + \mathbf{Z}_i\mathbf{b}_i + \boldsymbol{\varepsilon}_i, \quad i\in\indset{N}
\end{equation}
due to~\eqref{eq:PLMM} and Assumption~\ref{assumpt:D4-2}.
In particular, the adjusted and grouped responses in this model are independent in the sense that we have $\mathbf{R}_{\Yi} \independent \mathbf{R}_{\Yj}$ for $i\neq j$.
The strategy now is to first estimate the residuals with machine learning algorithms and then use
linear mixed model techniques to infer $\beta_0$. This is done with sample splitting and cross-fitting, and the details are described next.
\\

Let us define $\Sigma_{0} := \sigma_0^{-2}\Gamma_{0}$
and $\mathbf{V}_{0, i}:=(\mathbf{Z}_i\Sigma_{0}\mathbf{Z}_i^T + \mathds{1}_{n_i})$ so that we have
\begin{equation}\label{eq:resNormal}
	(\mathbf{R}_{\Yi} | \mathbf{W}_i,\Xi) \sim \mathcal{N}_{n_i}\big(\mathbf{R}_{\Xi}\beta_0, \sigma_0^2\mathbf{V}_{0, i}\big).
\end{equation}
We assume that there exist functions $m_X^0\colon\R^{v}\rightarrow\R^d$ and $m_Y^0\colon\R^{v}\rightarrow\R$ that we can apply row-wise to $\mathbf{W}_i$ to have
$\E[\Xi | \mathbf{W}_i] = m_X^0(\mathbf{W}_i)$ and $\E[\mathbf{Y}_i | \mathbf{W}_i] = m_Y^0(\mathbf{W}_i)$.
In particular, $m_X^0$ and $m_Y^0$ do not depend on the grouping index $i$.
Let $\eta^0:= (m_X^0,m_Y^0)$ denote the true unknown nuisance parameter.
Let us denote by $\theta_0 := (\beta_0, \sigma_0^2,\Sigma_{0})$ the complete true unknown parameter vector and by $\theta := (\beta,\sigma^2,\Sigma)$ and $\mathbf{V}_{i}:=\mathbf{Z}_i\Sigma\mathbf{Z}_i^T+\mathds{1}_{n_i}$ respective general parameters.
The log-likelihood of group $i$ is given by
\begin{equation}\label{eq:log-likelihood}
\begin{array}{rcl}
	\libig{\theta, \eta^0} &=& -\frac{n_i}{2} \log(2\pi) - \frac{n_i}{2}\log(\sigma^2) - \frac{1}{2}\log\big( \det(\mathbf{V}_{i})\big) \\
	&&\quad- \frac{1}{2\sigma^2} (\mathbf{R}_{\Yi}-\mathbf{R}_{\Xi}\beta)^T\mathbf{V}_{i}^{-1}(\mathbf{R}_{\Yi}-\mathbf{R}_{\Xi}\beta)  - \log \big(p(\mathbf{W}_i,\Xi)\big),
	\end{array}
\end{equation}
where $p(\mathbf{W}_i,\Xi)$ denotes the joint density of $\mathbf{W}_i$ and $\Xi$. We assume that $p(\mathbf{W}_i,\Xi)$ does not depend on $\theta$.
The true nuisance parameter $\eta^0$ in the log-likelihood~\eqref{eq:log-likelihood} is unknown and estimated with machine learning algorithms (see below).
Denote by $\eta := (m_X,m_Y)$ some general nuisance parameter. The terms that adjust $\Xi$ and $\mathbf{Y}_i$ for $\mathbf{W}_i$ with this general nuisance parameter are given by $\Xi-m_X(\mathbf{W}_i)$ and $\mathbf{Y}_i-m_Y(\mathbf{W}_i)$. Up to additive constants that do not depend on $\theta$ and $\eta$, we thus consider maximum likelihood estimation with the likelihood
\begin{displaymath}
	\begin{array}{rl}
	&\ell_i(\theta,\eta)
	= -\frac{n_i}{2}\log(\sigma^2) - \frac{1}{2}\log\big(\det(\mathbf{V}_{i}) \big)\\
	&\quad- \frac{1}{2\sigma^2} \Big(\mathbf{Y}_i-m_Y(\mathbf{W}_i)-\big(\Xi-m_X(\mathbf{W}_i)\big)\beta\Big)^T\mathbf{V}_{i}^{-1}\Big(\mathbf{Y}_i-m_Y(\mathbf{W}_i)-\big(\Xi-m_X(\mathbf{W}_i)\big)\beta\Big),
	\end{array}
\end{displaymath}
which is a function of both the finite-dimensional parameter $\theta$ and the infinite-dimensional
nuisance parameter $\eta$.
\\

Our estimator of $\beta_0$ is constructed as follows using double machine learning.
First, we estimate $\eta^0$ with machine learning
algorithms and plug these estimators into the estimating equations for $\theta_0$, equation~\eqref{eq:setToZero} below,
to obtain an estimator for $\beta_0$. This procedure is done with sample splitting and cross-fitting as explained next.

Consider repeated measurements from $N$ experimental units, subjects, or groups as in~\eqref{eq:PLMM}. Denote by $\mathbf{S}_i :=(\mathbf{W}_i,\Xi,\mathbf{Z}_i,\mathbf{Y}_i)$ the observations of group $i$.
First, we split the group indices $\indset{N}$ into $K\ge 2$ disjoint sets $I_1, \ldots, I_{K}$ of approximately equal size;
please see Section~\ref{sect:AssumptionsDefinitions} in the appendix for further details.

For each $k\in\indset{K}$, we estimate the conditional expectations $m_X^0(W)$ and $m_Y^0(W)$ with data from $I_{\kk}^c$. We call the resulting estimators  $\hat m_X^{I_{\kk}^c}$ and $\hat m_Y^{I_{\kk}^c}$, respectively.
Then, the adjustments $\widehat{\mathbf{R}}_{\Xi}^{I_{\kk}}:=\Xi-\hat m_X^{I_{\kk}^c}(\mathbf{W}_i)$, and $\widehat{\mathbf{R}}_{\mathbf{Y}_i}^{I_{\kk}}:=\mathbf{Y}_i-\hat m_Y^{I_{\kk}^c}(\mathbf{W}_i)$ for $i\inI_{\kk}$ are evaluated on $I_{\kk}$, the complement of $I_{\kk}^c$.
Let $\hat\eta^{I_{\kk}^c} := (\hat m_X^{I_{\kk}^c}, \hat m_Y^{I_{\kk}^c})$ denote the estimated nuisance parameter.
Consider the score function
	$\psi(\mathbf{S}_i; \theta,\eta) := \nabla_{\theta}\ell_i(\theta,\eta)$,
where $\nabla_{\theta}$ denotes the gradient with respect to $\theta$ interpreted as a vector.
On each set $I_{\kk}$, we consider an estimator $\hat\theta_{\kk} = (\hat\beta_{\kk},\hat\sigma_{\kk}^2, \hat{\Sigma}_{\kk})$ of $\theta_0$ that, approximately, in the sense of Assumption~\ref{assumpt:Theta2} in the appendix, solves
\begin{equation}\label{eq:setToZero}
	\frac{1}{\nn_{T, k}}\sum_{i\inI_{\kk}} \psibig{\mathbf{S}_i;\hat\theta_{\kk},\hat\eta^{I_{\kk}^c}}
	= \frac{1}{\nn_{T, k}}\sum_{i\inI_{\kk}} \nabla_{\theta}\ell_i(\theta,\eta)
	\stackrel{!}{=}  \boldsymbol{0},
\end{equation}
where $\nn_{T, k} := \sum_{i\inI_{\kk}}n_i$ denotes the total number of observations from experimental units that belong to the set $I_{\kk}$.
These $K$ estimators $\hat\theta_{\kk}$ for $k\in\indset{K}$ are assembled to form the final cross-fitting estimator
\begin{equation}\label{eq:betahat}
	\hat\beta := \frac{1}{K} \sum_{k=1}^{K} \hat\beta_{\kk}
\end{equation}
of $\beta_0$.
We remark that one can simply use linear mixed model computation and software to compute $\hat\beta_{\kk}$ based on the estimated residuals $\widehat{\mathbf{R}}^{I_{\kk}}$.
The estimator $\hat\beta$ fundamentally depends on the particular sample split. To alleviate this effect, the overall procedure may be repeated $\mathcal{S}$ times~\citep{Chernozhukov2018}. The $\mathcal{S}$ point estimators are aggregated by the median, and an additional term accounting for the random splits is added to the variance estimator of $\hat\beta$; please see
Algorithm~\ref{algo:Summary} that presents the complete procedure.

\begin{algorithm}[h!]
	\SetKwInOut{Input}{Input}
   	\SetKwInOut{Output}{Output}

 \Input{$N$ independent and identically distributed\  observations $\{\mathbf{S}_i=(\mathbf{W}_i,\Xi, \mathbf{Z}_i, \mathbf{Y}_i)\}_{i\in\indset{N}}$ from the model~\eqref{eq:PLMM} satisfying Assumption~\ref{assumpt:Distribution}, a natural number $K$, a natural number $\mathcal{S}$.}
 \Output{An estimator of $\beta_0$ in~\eqref{eq:PLMM} together with its estimated asymptotic variance.}

 \For{$s\in\indset{\mathcal{S}}$}
 {
 Split the grouped observation index set $\indset{N}$ into $K$ sets $I_1, \ldots, I_{K}$ of approximately equal size.

 \For{$k\inK$}
 {
 Compute the conditional expectation estimators $\hat m_X^{I_{\kk}^c}$ and $\hat m_Y^{I_{\kk}^c}$ with some
 machine learning algorithm and data from $I_{\kk}^c$.

 Evaluate the adjustments $\widehat{\mathbf{R}}_{\Xi}^{I_{\kk}}=\Xi-\hat m_X^{I_{\kk}^c}(\mathbf{W}_i)$ and $\widehat{\mathbf{R}}_{\mathbf{Y}_i}^{I_{\kk}}=\mathbf{Y}_i-\hat m_Y^{I_{\kk}^c}(\mathbf{W}_i)$ for $i\inI_{\kk}$.

Compute $\hat\theta_{\kk, s} = (\hat\beta_{\kk, s}, \hat\sigma_{\kk, s}^2, \hat\Sigma_{\kk, s})$ using, for instance, linear mixed model techniques.
 }

Compute $\hat\beta_{s} = \frac{1}{K}\sum_{k=1}^{K}\hat\beta_{\kk, s}$ as an approximate solution to~\eqref{eq:setToZero}.

Compute an estimate $\hat T_{0,s}$ of the asymptotic variance-covariance matrix $T_0$ in Theorem~\ref{thm:asymptoticGaussian}.
 }

 Compute $\hat\beta = \mathrm{median}_{s\in\indset{\mathcal{S}}}(\hat\beta_{s})$.

 Estimate $T_0$ by $\hat T_0 = \mathrm{median}_{s\in\indset{\mathcal{S}}}(\hat T_{0,s} + (\hat\beta-\hat\beta_{s})(\hat\beta-\hat\beta_{s})^T)$.\label{algo:varianceCorrection}

 \caption{Double machine learning in a partially linear mixed-effects model with repeated measurements.}\label{algo:Summary}
\end{algorithm}

\subsection{Theoretical Properties of the Fixed-Effects Estimator}

The estimator $\hat\beta$ as in~\eqref{eq:betahat} converges at the parametric rate, $N^{-1/2}$, and is asymptotically Gaussian distributed and semiparametrically efficient.

\begin{theorem}\label{thm:asymptoticGaussian}
Consider grouped observations $\{\mathbf{S}_i=(\mathbf{W}_i,\Xi,\mathbf{Y}_i)\}_{i\in\indset{N}}$ from the partially linear mixed-effects model~\eqref{eq:PLMM} that satisfy Assumption~\ref{assumpt:Distribution} such that $p(\mathbf{W}_i,\Xi)$ does not depend on $\theta$.
Let $\NN_T:=\sum_{i=1}^{N}n_i$ denote the total number of unit-level observations.
Furthermore, suppose the assumptions in Section~\ref{sect:AssumptionsDefinitions} in the appendix hold,
and consider the symmetric positive-definite matrix $T_0$ given in Assumption~\ref{assumpt:regularity7} in the appendix.
Then, $\hat\beta$ as in~\eqref{eq:betahat} concentrates in a $1/\surd{\NN_{T}}$ neighborhood of $\beta_0$ and is centered Gaussian, namely
\begin{equation}\label{eq:thmEquation}
	\surd{\NN_{T}}T_0^{\frac{1}{2}}(\hat\beta-\beta_0)
	\stackrel{d}{\rightarrow}\mathcal{N}_d(\boldsymbol{0}, \mathds{1}_{d}) \quad (N\rightarrow\infty),
\end{equation}
and semiparametrically efficient.
The convergence in~\eqref{eq:thmEquation} is in fact uniformly over the law $P$ of $\{\mathbf{S}_i=(\mathbf{W}_i,\Xi,\mathbf{Y}_i)\}_{i\in\indset{N}}$.
\end{theorem}

Please see Section~\ref{sect:asymptoticDistribution} in the appendix for a proof of Theorem~\ref{thm:asymptoticGaussian}.
Our proof builds on~\citet{Chernozhukov2018}, but we have to take into account the correlation within units that is introduced by the random effects.

The inverse asymptotic variance-covariance matrix $T_0$ can be consistently estimated;
see Lemma~\ref{lem:multiplierMatrix} in the appendix.
The estimator $\hat\beta$ is semiparametrically efficient because the score function comes from the log-likelihood of our data and because $\eta^0$ solves a concentrating-out equation for fixed $\theta$; see~\citet{Chernozhukov2018, Newey1994b}.

The assumptions in Section~\ref{sect:AssumptionsDefinitions} of the appendix specify regularity conditions and  required
convergence rates of the machine learning estimators. The
machine learning errors need to satisfy the product relationship
\begin{displaymath}
	\normP{m_X^0(W)-\hat m_X^{I_{\kk}^c}(W)}{2}\big(\normP{m_Y^0(W)-\hat m_Y^{I_{\kk}^c}(W)}{2}+ \normP{m_X^0(W)-\hat m_X^{I_{\kk}^c}(W)}{2} \big)\llN^{-\frac{1}{2}}.
\end{displaymath}
This bound requires that only the products of the
machine learning estimation errors
$\normP{m_X^0(W)-\hat m_X^{I_{\kk}^c}(W)}{2}$ and $\normP{m_Y^0(W)-\hat m_Y^{I_{\kk}^c}(W)}{2}$
but not the individual ones need to vanish at a rate smaller than
$N^{-1/2}$. In particular, the individual estimation errors  may vanish at the  rate
smaller than $N^{-1/4}$. This is achieved by many machine learning
methods (cf.
\citet{Chernozhukov2018}):
$\ell_1$-penalized and related methods in a variety of sparse models
\citep{Bickel2009, Buehlmann2011, Belloni2011, Belloni-Chernozhukov2011, Belloni2012, Belloni-Chernozhukov2013}, forward selection in sparse models
\citep{Kozbur2020}, $L_2$-boosting in sparse linear models
\citep{Luo2016}, a class of regression trees and random forests
\citep{Wager2016}, and neural networks \citep{Chen1999}.

We note that so-called Neyman orthogonality makes score functions insensitive to inserting potentially biased machine learning
estimators of the nuisance parameters. A score function is Neyman orthogonal if its Gateaux derivative vanishes at the true $\theta_0$ and the true $\eta^0$. In particular, Neyman orthogonality is a first-order property. The product relationship of the machine learning
estimating errors described above is used to bound second-order terms. We refer to Section~\ref{sect:asymptoticDistribution} in the appendix for more details.

\section{Numerical Experiments}\label{sect:numerical-experiments}

We apply our method to an empirical and a pseudorandom dataset and in a simulation study.
Our implementation is available in the \textsf{R}-package \texttt{dmlalg}~\citep{dmlalg}.

\subsection{Empirical Analysis: CD4 Cell Count Data}

Subsequently, we apply our method to longitudinal CD4 cell counts data collected from human immunodeficiency virus (HIV) seroconverters. This data has previously been analyzed by~\citet{Zeger-Diggle1994} and is available in the \textsf{R}-package \texttt{jmcm}~\citep{jmcm} as \texttt{aids}.
It contains $2376$ observations of CD4 cell counts measured on $369$ subjects. The data was collected during a period ranging from $3$ years before to $6$ years after seroconversion. The number of observations per subject ranges from $1$ to $12$, but for most subjects,  $4$ to $10$ observations are available.
Please see~\citet{Zeger-Diggle1994} for more details on this dataset.

Apart from time, five other covariates are measured: the age at seroconversion in years (age), the smoking status measured by the number of cigarette packs consumed per day (smoking), a binary variable indicating drug use (drugs), the number of sex partners (sex), and the depression status measured on the Center for Epidemiologic Studies Depression (CESD) scale (cesd), where higher CESD values indicate the presence of more depression symptoms.

We incorporate a random intercept per person.
Furthermore, we consider a square-root transformation of the CD4 cell counts to reduce the skewness of this variable as proposed by~\citet{Zeger-Diggle1994}.
The CD4 counts are our response. The covariates that are of scientific interest are considered as $X$'s,
and the remaining covariates are considered as $W$'s in the partially linear mixed-effects model~\eqref{eq:PLMM}.
The effect of time is modeled nonparametrically, but there are several options to model the other covariates.
Other models than partially linear mixed-effects model
have also been considered in the literature to analyze this dataset.
For instance, \citet{Fan2000} consider a functional linear model where the linear coefficients are a function of the time.
\\

We consider two partially linear mixed-effects models for this dataset.
First, we incorporate all covariates except time linearly.
Most approaches in the literature considering a partially linear mixed-effects model for this data that model time nonparametrically
report that sex and cesd are significant and that either smoking or drugs is significant as well; see for instance~\citet{Zeger-Diggle1994, Taavoni2019b, Wang2011}.
\citet{Guoyou2008} develop a robust estimation method
for longitudinal data and estimate
nonlinear effects from time with regression splines.
With the CD4 dataset,
They find that smoking and cesd are significant.

We apply our method with $K=2$ sample splits, $\mathcal{S}=100$ repetitions of splitting the data, and learn the conditional expectations with random forests that consist of $500$ trees whose minimal node size is $5$.

Like~\citet{Guoyou2008},
we conclude that smoking and cesd are significant; please see the first row of Table~\ref{tab:empirical} for a more precise account of our findings.
Therefore,
we can expect
that our method implicitly performs robust estimation. Apart from sex, our point estimators are larger or of about the same size in absolute value as what~\citet{Guoyou2008} obtain. This suggests that our method incorporates potentially less bias. However, apart from age, the standard deviations are slightly larger with our method.
This can be expected because random forests are more complex than the regression splines~\citet{Guoyou2008} employ.
\\

\begingroup
\begin{table}
\footnotesize
\centering
\begin{tabular}{| L{2.5cm} | C{1.9cm}  | C{1.9cm} | C{1.9cm} | C{1.9cm} | C{2.2cm} |}
\hline
&  \textbf{age} & \textbf{smoking} & \textbf{drugs} & \textbf{sex} & \textbf{cesd}\\
\hline
\hline
$W = (\mathrm{time})$ &  $0.004$ ($0.027$) & $0.752$ ($0.123$) & $0.704$ ($0.360$) & $0.001$ ($0.043$)    &   $-0.042$ ($0.015$)\\
\hline
$W = (\mathrm{time}, \mathrm{age}, \mathrm{sex})$ & - & $0.620$ ($0.126$) & $0.602$ ($0.335$) &- & $-0.047$ ($0.015$) \\
\hline
\hline
\citet{Zeger-Diggle1994}
& $0.037$ ($0.18$) & $0.27$ ($0.15$) & $0.37$ ($0.31$) & $0.10$ ($0.038$) &  $-0.058$ ($0.015$) \\
\hline
\citet{Taavoni2019b}
& $1.5\cdot 10^{-17}$ ($3.5\cdot 10^{-17}$)& $0.152$ ($0.208$) & $0.130$ ($0.071$) &  $0.0184$ ($0.0039$) & $-0.0141$ ($0.0061$)  \\
\hline
\citet{Wang2011}
& $0.010$ ($0.033$) & $0.549$ ($0.144$) & $0.584$ ($0.331$) & $0.080$ ($0.038$) & $-0.045$ ($0.013$)\\
\hline
\citet{Guoyou2008}
& $0.006$ ($0.038$) &  $0.538$ ($0.136$) & $0.637$ ($0.350$)& $0.066$ ($0.040$) & $-0.042$ ($0.015$)\\
\hline
\end{tabular}
\caption{\label{tab:empirical}Estimates of the linear coefficient and its standard deviation in parentheses with our method
for nonparametrically adjusting for time (first row) and for time, age, and sex (second row).
The remaining rows display the results from~\citet[Section 5]{Zeger-Diggle1994}, \citet[Table 1, ``Kernel'']{Taavoni2019b}, \citet[Table 2, ``Semiparametric efficient scenario I'']{Wang2011}, and \citet[Table 5, ``Robust'']{Guoyou2008}, respectively.
}
\end{table}
\endgroup

We consider a second estimation approach where we model the variables time, age, and sex nonparametrically and allow them to interact.
It is conceivable that these variables are not (causally) influenced by smoking, drugs, and cesd and that they are therefore exogenous.
The variables smoking, drugs, and cesd are modeled linearly,
and they are considered as treatment variables. Some direct causal effect interpretations are possible if one is willing to assume, for instance, that the nonparametric adjustment variables are causal parents of the linear variables or the response.
However, we do not pursue this line of thought further.

We estimate the conditional expectations given the three nonparametric variables time, age, and sex again with random forests that consist of $500$ trees whose minimal node size is $5$ and use $K=2$ and $\mathcal{S}=100$ in Algorithm~\ref{algo:Summary}.
We again find that smoking and cesd are significant; please see the second row of Table~\ref{tab:empirical}.
This cannot be expected a priori because this second model incorporates more complex adjustments, which can lead to less significant variables.

\subsection{Pseudorandom Simulation Study: CD4 Cell Count Data}\label{sect:pseudorandom}

Second, we consider the CD4 cell count data from the previous subsection and perform a pseudorandom simulation study. The variables smoking, drugs, and cesd are modeled linearly and the variables time, age, and sex nonparametrically.
We condition on these six variables in our simulation. That is, they are the same in all repetitions.
The function $g$ in~\eqref{eq:PLMM} is chosen as a regression tree that we built beforehand.
We let $\beta_0=(0.62, 0.6, -0.05)^T$, where the first component corresponds to smoking, the second one to drugs, and the last one to cesd, consider a  standard deviation of the random intercept per subject of $4.36$, and a standard deviation of the error term of $4.35$. These are the point estimates of the respective quantities obtained in the previous subsection.

Our fitting procedure uses random forests consisting of $500$ trees whose minimal node size is $5$ to estimate the conditional expectations, and we use $K=2$ and $\mathcal{S}=10$ in Algorithm~\ref{algo:Summary}. We perform $5000$ simulation runs.
We compare the performance of our method with that of the spline-based function \texttt{gamm4} from the package~\texttt{gamm4}~\citep{gamm4} for the statistical software \textsf{R}~\citep{R}.
This method represents the nonlinear part of the model by smooth additive functions and estimates them by penalized regression splines. The penalized components are treated as random effects and the unpenalized components
as fixed.

The results are displayed in Figure~\ref{fig:pseudorandom}.
With our method,  \texttt{mmdml}, the two-sided confidence intervals for $\beta_0$ are of about the same length  but achieve a coverage that is closer to the nominal $95\%$ level than with \texttt{gamm4}. The \texttt{gamm4} method largely undercovers the packs component of $\beta_0$, which can be explained by the incorporated bias.

\begin{figure}[]
	\centering
	\caption[]{\label{fig:pseudorandom}
	Coverage and  length of two-sided confidence intervals at significance level 5\% and  bias for our method, \texttt{mmdml}, and
	\texttt{gamm4}.
	In the coverage plot, solid dots represent point estimators, and circles represent $95\%$ confidence bands with respect to the $5000$ simulation runs.
	The confidence interval length and bias are displayed with boxplots without outliers.}
	\includegraphics[width=0.75\textwidth]{pics/pseudoRand_paper_pics__25_January_2022__17_35_42__pseudorandom__Xforest_Yforest_S10_reps5000_FINAL.pdf}
\end{figure}

\subsection{Simulation Study}\label{sect:simul-study}

We consider a partially linear mixed-effects model with $q=3$
random effects and where $\beta_0$ is one-dimensional.
Every subject has their own random intercept term and a nested random effect with two levels.
Thus, the random effects structure is more complex than in the previous two subsections because these models only used a random intercept.
We compare three data generating mechanisms:
One where the function $g$ is nonsmooth and the number of observations per group is balanced,
one where the function $g$ is smooth and the number of observations per group is balanced,
and one where the function $g$ is nonsmooth and the number of observations per group is unbalanced;
please see Section~\ref{sect:dataSimulation} in the appendix for more details.

We estimate the nonparametric nuisance components, that is, the conditional expectations, with
random forests consisting of $500$ trees whose minimal node size is $5$.
Furthermore, we use $K=2$ and $\mathcal{S}=10$ in Algorithm~\ref{algo:Summary}.

We perform $1000$ simulation runs and consider different numbers of groups $N$. As in the previous subsection, we compare the performance of our method with
\texttt{gamm4}.

The results are displayed in Figure~\ref{fig:simulation}. Our method, \texttt{mmdml}, highly outperforms \texttt{gamm4} in terms of coverage for nonsmooth $g$ because the coverage of \texttt{gamm4} equals $0$ due to its substantial bias. Our method overcovers slightly due to the correction factor that results from the $\mathcal{S}$ repetitions. However, this correction factor is highly recommended in practice.
With smooth $g$, \texttt{gamm4} is closer to the nominal coverage and has shorter confidence intervals than our method. Because the underlying model is smooth and additive, a spline-based estimator is better suited. In all scenarios,  our method outputs longer confidence intervals than \texttt{gamm4} because we use random forests; consistent with theory, the difference in absolute value decreases though when $N$ increases.

\begin{figure}[]
	\centering
	\caption[]{\label{fig:simulation}
	Coverage and median length of two-sided confidence intervals for $\beta_0$ at significance level 5\% (true $\beta_0 = 0.5$) and median bias
	for three data generating scenarios for our method, \texttt{mmdml}, and
	\texttt{gamm4}.
	The shaded regions in the coverage plot represent $95\%$ confidence bands with respect to the $1000$ simulation runs.
	The dots in the coverage and bias plot are jittered, but neither are their interconnecting lines nor their confidence bands.}
	\includegraphics[width=\textwidth]{pics/coverage_paper_pics__26_January_2022__99_99_99__simulation__Xforest_Yforest_S10_reps1000_FINAL_merged_1.pdf}
\end{figure}

\section{Conclusion}\label{sect:conclusion}

Our aim was to develop inference for the linear coefficient $\beta_0$ of a partially linear mixed-effects model
that includes a linear term and potentially complex nonparametric terms.
Such models can be used to describe heterogenous and correlated data that feature
some grouping structure, which may result from taking repeated measurements.
Traditionally, spline or kernel approaches are used to cope with the nonparametric part of such a model. We presented a scheme that uses the double machine learning
framework of~\citet{Chernozhukov2018} to estimate any
nonparametric components with arbitrary machine learning algorithms. This allowed us to consider complex nonparametric components with interaction structures and high-dimensional variables.

Our proposed method is as follows. First, the nonparametric variables are regressed out from the response and the linear variables. This step adjusts the response and the linear variables  for the nonparametric variables and may be performed with any machine learning algorithm.
The adjusted variables satisfy a linear mixed-effects model, where the linear coefficient $\beta_0$ can be estimated with standard linear
mixed-effects
techniques. We showed that the estimator of $\beta_0$
asymptotically follows a Gaussian distribution, converges at the parametric rate, and is semiparametrically efficient.
This asymptotic result
allows us to perform inference for $\beta_0$.

Empirical experiments demonstrated the performance of our proposed method.
We conducted an empirical and pseudorandom data analysis and a  simulation study.
The simulation study and the pseudorandom experiment confirmed the effectiveness of our method in terms of coverage, length of confidence intervals, and estimation bias
compared to a penalized regression spline approach relying on additive models.
In the empirical experiment, we analyzed longitudinal CD4 cell counts data collected from
HIV-infected individuals.
In the literature, most methods only incorporate the time component nonparametrically to analyze this dataset.
Because we estimate nonparametric components with machine learning algorithms, we can allow several variables to enter the model nonlinearly, and we can allow these variables to interact.
A comparison of our results with the literature suggests that our method
may perform robust estimation.

Implementations of our method are available in the \textsf{R}-package \texttt{dmlalg}~\citep{dmlalg}.

\section*{Acknowledgements}

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 786461).

\phantomsection
\addcontentsline{toc}{section}{References}


\section{Introduction}
\label{sect:intro}
Grayscale channels, which reflect image luminance, are used for various applications such as printing, tone mapping, data compression, and feature extraction. Thus, obtaining luminance along with human perception has a key role for decolorization, which converts RGB channels to high-quality gray ones. For example, High Dynamic Range (HDR) compression is ideally performed by tone mapping the luminance channel for the lower computational and memory cost. However, applying well-known luminance channels such as Y of YCbCr or V of HSV does not guarantee appropriate tone mapping, as these channels do not reflect human perceptions. Therefore, decolorization has gathered considerable attention and various sophisticated methods to achieve perceptual decolorization have recently been proposed. These methods can be classified into global and local methods. Global methods can define only one conversion function for all pixels, and most of these methods use all pixels in the image to determine the function. On the other hand, local ones process the target from neighboring pixels in the same way as a spatial filter, the function is different for each pixel. However, both types of methods face the issue of calculation cost, which comes from optimization iterations or spatial filter processing.

We have developed a fast decolorization method that reflects the perception of warm and cool colors which is well known in psychophysics studies\cite{morton1998color}. Colors are arranged according to their wavelengths on the color wheel, the ones with longest wavelengths are on the right side of the wheel and are known as warm colors, as they evoke warmth. These hues include shades of red, yellow, and orange. On the other hand, green, blue and violet which have shorter wavelengths are placed on the left side of the color wheel, and are perceived as cool colors. The color of an object in a scene affects our perception of its apparent depth and this phenomenon has been exploited by many artists. This optical illusion has been studied by psychologists, and early researchers explored color-depth relationship. One of the widely accepted theory explains this phenomenon is due to fact that shorter wavelengths of visible light are refracted more than longer wavelengths \cite{sundet1978effects}. In other words, an equidistant source of different wavelengths cannot be focused simultaneously onto our retina. This phenomena is called as chromatic aberration and we discuss it in detail in section \ref{ca}. In our decolorization method, we combine the idea of warm/cool colors with the Helmholtz-Kohlrausch effect \cite{Nayatani1}. On that account, we make two assumptions, which are: $(i)$ warm colors (mainly including $R$) are lighter than $Y$ of YCbCr and $(ii)$ mixed colors are darker than the $Y$ or $L$ of CIE with the same luminance. To satisfy these assumptions, we use a weighted blending of RGB channels and remap them to warm/cool colors on the luminance channel. Following are our main contributions:

\begin{itemize}
	\item We achieve high-speed color mapping by exact pixel-by-pixel processing.
\end{itemize}

\begin{itemize}
	\item We obtain luminance comparable to that of optimization-based methods.
\end{itemize}

\begin{itemize}
	\item We demonstrate effective luminance distribution by performing objective evaluations and a subjective user study.
\end{itemize}

\begin{figure*}[t]
	\centering
	\includegraphics[width = \textwidth]{teaser1.pdf}
	\caption{Main concept of our decolorization method: Human perception of warm and cool colors. Warm colors ``advance'' toward the eye, while cool colors ``recede''. In this work we are able to accurately reflect the human perception of color temperatures, whereas this phenomenon is non-existent in conventional YCbCr color space.}
	\label{fig:1}
\end{figure*}

There are many well defined methods to convert any color image to a grayscale image. An effortless procedure is to assign different weights to color channels, in order to have the same luminance in the grayscale image as the original color image. For example, in the MATLAB function \texttt{rgb2gray}, it converts any RGB values to grayscale ($Gray$) values by forming a weighted sum of the $R$, $G$, and $B$ components as $Gray = 0.2989\times R + 0.5870\times G + 0.1140\times B$. This function operates under an assumption that human visual system is more sensitive to green color. When operating with CIELab and YUV color spaces, one could directly obtain luminance channel as the grayscale version of the color image as they consider the luminance and color channel to be independent. But, such crude approaches will fail to preserve image contrast as shown in these examples (see Fig. \ref{fig:1} and Fig. \ref{fig:2}).

\begin{figure}[!t]
	\centering
	\includegraphics[width = \linewidth]{color.pdf}
	\caption{Comparison of luminance components obtained using YCbCr, CIELAB and our proposed method.}
	\label{fig:2}
\end{figure}

In several real-world image/video processing applications like detail enhancement, image matching, and segmentation under different illumintion a 1-D grayscale image has to be obtained from its corresponding 3-D color image. However, mapping the 3-D color information onto a 1-D grayscale image while retaining the original contrast and fine details is a challenging problem. Additionally, implementing decolorization algorithms with a reasonable computational efficiency is pivotal for realising their real-time applications. Many studies have been carried out to develop novel decolorization methods. These mapping methods can be categorized into global
\cite{gooch2005color2gray,grundland2007decolorize,smith2008apparent,kim2009robust,lu2012contrast} and local methods \cite{ancuti2010image,ancuti2016laplacian,liu2015gcsdecolor,song20161}. In local mapping methods, the same color pixel within an image could be mapped into different grayscale values depending on its spatial location. Ideally this undesirable as such output images may be perceived as unnatural. On the other hand, in global
mapping methods same color pixels within an image irrespective of its spatial location are mapped to same grayscale values. Thus, global methods are more likely to produce grayscale images that are perceived to appear natural.

In the global methods category, Gooch et al. \cite{gooch2005color2gray} proposed a global decolorization algorithm that can be implemented by solving the optimization problem for all image pixels. Then, Kim et al. \cite{kim2009robust} aimed at high-speed processing by simplifying Gooch's method. Smith et al. \cite{smith2008apparent} used unsharp masking and the Helmholtz-Kohlrausch effect (H-K effect) model of Nayatani et al. \cite{Nayatani1}. Nayatani's model \cite{Nayatani1} is merely an experimental model for the effect of the CIELUV chrominance component for human perception. On the other hand, the method of Lu et al. \cite{lu2012contrast} is focused on converting a color image into a gray image with high contrast. The main advantage of these methods is transformation consistency, i.e., the same color is converted to the same grayscale. However, speed remains a problem for  these methods. Most of the local methods are aimed at speeding up the method of Lu et al \cite{lu2012contrast}. To enhance image contrast, Ancuti et al. \cite{ancuti2016laplacian} adopted the strategy of using a Laplacian filter and Song et al. \cite{song20161} used a Sobel filter. Although local methods are effective in terms of contrast emphasis, they are disadvantageous in terms of tone mapping because
conversion consistency is not maintained and it differs from human perception.

Recently, some machine learning-based techniques have also been proposed for image decolorization  \cite{liu2019image,cai2018perception,hou2017deep,lin2008learning,zhang2018contrast}. Cai et al. \cite{cai2018perception} proposed a method, that used the perceptual loss function to pretrain VGG-Net \cite{simonyan2018very}. However, it is difficult to control and many of their output images are far from human perception additionally, the computational cost are high. Processing an image of $256 \times 256$ size it requires roughly 30 seconds on a single Nvidia GeForceGTX 1080 GPU. Zhang et al., proposed a CNN framework that combines local and global image features \cite{zhang2018contrast}. However, their network framework do not account for exposure features \cite{liu2019image}. Lin et al.'s method \cite{lin2008learning} by utilizing a database of 50 images from the Corel dataset produced 50 grayscale images using the Color2Gray algorithm \cite{gooch2005color2gray}. With these 50 input/output image pairs as training examples for their partial differential equations-based (PDE)learning system, the learn color2gray mapping. The PDE system produced images of comparable quality to \cite{gooch2005color2gray}. However, for an input image of size $n \times n$ their PDE color mapping algorithm's computational complexity is $\textit{O}(n^{2})$. Under these circumstances, a high-speed method for
generating grayscale images that accurately captures human perception has not yet been developed. To make tone maps valid, it is necessary to develop such a method.

\section{Proposed Method}

\subsection{Problem Definition}
Luminance components such as $Y$ of YCbCr and CIE $L$ have been used in various image processing applications; however, they do not accurately reflect human perception (Fig.\ref{fig:1} and Fig. \ref{fig:2}). Figure \ref{fig:1}  shows how the warm colored flower (red) advances towards the eye of the observer, while the background mainly green recedes. Using the luminance channel of the conventional YCbCr color space this phenomenon is absent. But, in our method we are able to capture the warmth $R$ (red) component as perceived in human perception. Figure \ref{fig:2} shows that the $R$ (red) and $B$ (blue) components do not come even close to the perception in the luminance component of CIELAB. Moreover, mixed color components such as mud yellow tend to appear dark for people. In this study, we conceived the idea that RGB weight functions for alpha blending can reproduce this phenomenon.

\subsection{Luminance mapping using red and blue weighting function}

\begin{figure}[!t]
	\centering
	\includegraphics[width = \linewidth]{figure3PDF.pdf}
	\caption{(a) Remodeling of luminance space based on our two weighting functions one each for warm and cool colors. Here, luminance is defined as the Euclidean distance of warm/cool color. (b) Decolorization of a color chart illustrating our method's effectiveness in generating brighter shades of gray for warm colors.}
	\label{fig:3}
\end{figure}

We considered the idea of color mapping by combining the warm/cool color and the Helmholtz-Kohlrausch \cite{Nayatani1}. Psychophysical studies find that, warm and cool colors impact our visual perception of the objects that we see. For example, the red color associated with fire/sun advances toward the eye, creates an illusion of heat and therefore perceived as warmth and comforting. On the other hand, cool colors have reverse effects of warm colors. Receding from the eye of the observer, cool colors reminds of the earthy objects, like meadows and oceans. These hues often are perceived as cool and refreshing \cite{morton1998color}\cite{ou2004study}. In our decolorization method, we developed two weighting functions as shown in Fig \ref{fig:3}. One function for remapping warm colors and the other for remapping cool colors. In our method, actual luminance is defined as the Euclidean distance of weighted warm/cool luminance including the $W$ (white) channel. Essential luminance is given by

\begin{footnotesize}
	\begin{equation}
		L_{WHITE}= \sqrt{\frac{R^{2}+G^{2}+B^{2}}{3}}
		\label{eq:Lw1}
	\end{equation}
	\begin{equation}
		L_{B} = B
		\label{eq:Lb}
	\end{equation}
	\begin{equation}
		L_{G,R}= \sqrt{\beta_{R}R^{2}+(1-\beta_{R})G^{2}}, \; (0.5 < \beta_{R} < 1)
		\label{eq:Lgr}
	\end{equation}
\end{footnotesize}
As we know, $L_{WHITE}:W$ is the Euclidean distance of $RGB$ channels; $L_{B}$ is the $B$ channel as it is. The warm color function is also defined as the Euclidean distance by $\beta_{R}$ weighted $R$ and $G$(green) vectors.
The $\beta_{R}$ determines the bias of $R$ components in the replacement; when the $\beta_{R}$ becomes higher, more $R$ components than $G$ components are rated. The remappings made by the blending function using the $R$ and $B$
ratio of $RGB$ values are given by

\begin{footnotesize}
	\begin{equation}
		L_{WARM} = \frac{R}{R+G+B}\cdot L_{G,R} + \left(1-\frac{R}{R+G+B}\right)\cdot L_{WHITE}
		\label{eq:Lwarm}
	\end{equation}
	\begin{equation}
		L_{COOL} = \left(1-\frac{B}{R+G+B}\right)\cdot L_{B} + \frac{B}{R+G+B}\cdot L_{WHITE}
		\label{eq:Lcool}
	\end{equation}
\end{footnotesize}

\begin{figure}[!t]
	\centering
	\includegraphics[width = \linewidth]{betaRK1.pdf}
	\caption{Impact of $\beta_{R}$ and $\beta_{k}$ on output luminance value. (top) $\beta_{k}$ set to 0.8 and for changing $\beta_{R}$ we move from $G_{max}$ to $R_{max}$. (middle) $\beta_{R}$ set to 0.55 and for changing $\beta_{k}$ we move from $B_{max}$ to $R_{max}$. (bottom) $\beta_{R}$ set to 0.55 and for changing $\beta_{k}$ we move from $B_{max}$ to $G_{max}$}
	\label{fig:4}
\end{figure}

\begin{figure*}[!t]
	\centering
	\includegraphics[width = \textwidth]{Compimage.pdf}
	\caption{Color to gray conversion comparison with other decolorization methods.}
	\label{fig:5}
\end{figure*}

In the $L_{WARM}$, using $R$ weighting, we can map $\frac{G}{R}$ to $low/high$ luminance separately. In Inverse $B$ weighting, pure $B$ components are assigned to low luminance in the $L_{COOL}$ . Since both functions are
blended with $L_{WHITE}$, bright orange/yellow and sky blue, which include high white components, are mapped to higher luminance. Finally, we obtain the luminance channel $L$, which is given by

\begin{footnotesize}
	\begin{equation}
		L= \sqrt{\beta_{k}L^{2}_{WARM}+(1-\beta_{k})L^{2}_{COOL}} \; (0.5 < \beta_{k} < 1)
		\label{eq:Lw2}
	\end{equation}
\end{footnotesize}
In this study, we mainly used warm-color weighting luminance in experiments and set the $\beta_{k}$ higher. We set two parameters relating to color component emphasis as follows: $(\beta_{R} = 0.55; \beta_{k} = 0.8)$, and from Fig \ref{fig:4} it can be easily understood how $\beta_{R}$ and $\beta_{k}$ can bias the resulting luminance values.

\subsection{Evaluation of our algorithm implementation}

When evaluating our method, we focused on the following points:
 \begin{itemize}
 	\item Visual comparison with other decolorization methods
 	\item Compare their processing speeds.
\end{itemize}
 	The evaluations demonstrated that our method does reflect human perception better than or equal to other optimized methods (see perceptual test). It delivers high-speed processing, and is a useful tool for many image processing applications. Figure \ref{fig:5} shows a comparison between our proposed method and other decolorization methods, and the subjective evaluation is shown in Fig. \ref{fig:vt}. For this comparison we have used images from Cadik’s dataset \cite{cadik2008perceptual}. In Fig. \ref{fig:5}  first and second columns from the left are original images and images obtained with our method respectively. The images in the third column from the left were obtained with the $L$ component of the CIELAB color space, which is a reversible model reflecting human visual characteristics. The images in the fourth to sixth columns were obtained from global decolorization methods, which means only one conversion function is applied to each pixel. For example, in these methods, optimization techniques referring to whole pixels are applied to luminance conversion without regard to the brightness perceived by human perception. The images in the far right column were obtained with a local method that refers to pixel values in the local patch of the image for contrast enhancement.

\begin{figure}[!t]
	\centering
	\includegraphics[width = \linewidth]{subjective.pdf}
	\caption{Subjective visual test survey and mean opinion score}
	\label{fig:vt}
\end{figure}

\begin{figure*}[!t]
	\centering
	\includegraphics[width = \textwidth]{VisTestTMO2.pdf}
	\caption{Effects of decolorization techniques on tone mapped images that are generated using Shan et al.'s tone mapping operator \cite{shan2009globally}.}
	\label{fig:tmo_1}
\end{figure*}

In the $L$ component of the CIELAB color space, the base luminance is $Y$ in the YCbCr color space; the luminance is also mapped along with the color order of $Y$. In the image (third column, fourth row), since the luminance of all colors are the same in $Y$, output values are also the same in the CIELAB color space. However, we can perceive the contrast in this image. Thus, an appropriate conversion is required; our method can generate the perceived contrast in this image. Gooch et al. \cite{gooch2005color2gray} obtained the highest average C2G-SSIM score but clearly a step artifact occurs in the gradation image (fourth column, second row). Kim et al. \cite{kim2009robust} proposed an improved version of Gooch's method that achieves high-speed optimization and reflects the Helmholtz-Kohlrausch \cite{Nayatani1}. In this method, since a weighting function is applied for expanding luminance distribution of whole pixel colors in the image along with chrominance, its conversion becomes different in each image. Thus, over-enhancement is observed in the images (fifth column, first row) and (fifth column, fourth row). The output images
of our method are similar to those of Smith et al. \cite{smith2008apparent}. Their method also uses the Helmholtz-Kohlrausch, but it also requires a lot of processing time for post-unsharp-mask filtering. In Lu et al. \cite{lu2012contrast}, their method does not reflect human perception; they try generating high contrast images for mask images that are input to an edge-preserving filter such as a guided filter.

\begin{table*}[!t]\scriptsize
	\renewcommand{\arraystretch}{1.2}
	\caption{Run-time comparison table with other decolorization algorithms.}
		\centering
	\label{table:1}
	\begin{tabular}{|p{6em}|p{6.5em}|p{7em}|p{7.5em}|p{5em}|p{6.5em}|p{8em}|}
		\hline
		Algorithm &  Processing time  & Image size W$\times$H & CPU Clock Speed &  Optimization &  Process & Normalized Time$^\star$ \\ \hline
		Gooch et al. \cite{gooch2005color2gray}   &  25.7s & 200$\times$200    & -GPU- & \checkmark &Global & N/A\\ \hline
		Kim et al.  \cite{kim2009robust} &  102ms & 320$\times$240    &  2.66GHz & \checkmark & Global & 1.30$\mu$s \\ \hline
		Smith et al. \cite{smith2008apparent}   &  6.7s & 570$\times$593    & 3.0GHz &$\times$ &Global+Local & 22.02$\mu$s \\ \hline
		Lu et al.  \cite{lu2012contrast} &  800ms & 600$\times$600    &  3.80GHz & \checkmark &Global Contrast & 3.12$\mu$s \\ \hline
		Song et al.  \cite{song20161} &  40ms & 320$\times$240    & N/A & $\times$ &Local Contrast&N/A\\ \hline
		Ancuti et al. \cite{ancuti2016laplacian}  &  100ms & 800$\times$600   & 2.5GHz & $\times$ &Local Contrast & 0.19$\mu$s \\ \hline %
		L of CIE Lab    &  25.57ms & 800$\times$600    &2.7GHz &$\times$ & Global & 0.053$\mu$s \\ \hline
		Ours low res &   \textbf{16.71ms}    & 800$\times$600    & 2.7GHz &$\times$ &Global& \textbf{0.034$\mu$s}\\
		Ours high res &  \textbf{202.05ms}    & 3008$\times$2008   & 2.7GHz &$\times$ &Global& \textbf{0.033$\mu$s}\\
		\hline
	\end{tabular}
	\flushleft
	$^\star$ Normalized time is the processing time normalized by frequency (2.7 GHz) and divided by the number of pixels. It indicates the effective processing time per pixel.
\end{table*}

For our perceptual evaluation we conducted a user study with two objectives:
\begin{itemize}
	\item Evaluate the decolorization process by measuring accuracy and preference \cite{cadik2008perceptual,lu2014contrast}.
\end{itemize}
\begin{itemize}
	\item Assess the effectiveness of our decolorization method as an image pre-processing tool for tone mapping application.
\end{itemize}
Our study group of 15 students (9 males, 6 females, average age = 23) were shown twelve sets of images (refer Figs. \ref{fig:5} and \ref{fig:tmo_1} ). They were asked to evaluate Fig. \ref{fig:5} and assign points to these images, on a scale of 1 (low) to 5 (high) for accuracy and preference. For our first evaluation, to measure the accuracy they compared the original color images to their decolorized output image. For the preference measurements, the user group compared the decolorized images only, i.e., without the corresponding color image. In our next user study we measured the naturalness and overall quality of the tone mapped images. For tone mapping we chose a MATLAB algorithm of Shan et al. \cite{shan2009globally} as they can tonemap standard images and it can be performed by simple parameter settings. Color tone mapped images for the same set shown in Fig. \ref{fig:tmo_1} were evaluated for naturalness and perceptual quality, i.e., the users rated 5 for images that were perceived as natural and observed minimal artifacts. In our test setup we used a Microsoft Surface Pro 7 set to its native resolutions and the lighting of the test room was slightly dim. The compiled response of the volunteers are presented in Fig. \ref{fig:vt}.

Table \ref{table:1} lists the processing speed of various methods, our proposed implementation and the CIELAB were implemented in C++; these codes were executed on an Intel Core i5-5257U (2.70GHz) CPU without any multicore, multithread or SIMD operations. The results confirmed that our method had the fastest runtime among the methods compared. It is worth noticing that it exceeded CIELAB in runtime; this indicates it also has advantages in total calculation cost including post-processing. Its computational complexity is only $\textit{O}(1)$ because it performs exact pixel by pixel processing, referring only to the RGB value at each pixel. Among other methods, the one developed by Gooch et al. was reported to have $\textit{O}(n^{4})$ computational complexity. Using $\textit{O}(1)$ spatial filtering it is possible to develop $\textit{O}(1)$ local methods \cite{ancuti2016laplacian,song20161}, but the filter calculations required would degrade their run-time in comparison to our decolorization method. The proposed color to gray technique has demonstrated faster run-time than local methods by maintaining global coherence, which means the conversions were the same in all pixels.

\begin{figure}[t]
	\centering
	\includegraphics[width =\linewidth]{epsFig.pdf}
	\caption{Comparison of six color to gray methods based on the E-score.}
	\label{fig:plotE}
\end{figure}

\subsection{Experimental results}
In our experiments, we utilized the color250 dataset which comprises of 250 natural and synthetic color images \cite{lu2014contrast}. To quantitatively evaluate our decolorization algorithm we choose two objective metrics: E-score and C2G-SSIM by Ma et al \cite{ma2015objective}.
E-score is a joint measure proposed by Lu et al., a harmonic mean which is computed by combining two metrics: Color Contrast Preserving Ratio (CCPR), and Color Content Fidelity Ratio (CCFR) \cite{lu2014contrast}. The CCPR is useful in maintaining the color contrast in decolorization images which is perceivable to humans. Specifically, when the color difference is smaller than a certain threshold value, it becomes undetectable to humans. Furthermore, CCFR estimates if the decolorization image is accurate in terms of structures when compared to the original color image. C2G-SSIM is new color to gray objective evaluation metric based on the Structural Similarity (SSIM) index quality metric \cite{wang2004image}. The C2G-SSIM generates quality map and has good correlation with HVS subjective preference. Table \ref{tab_om} presents the average E-score and C2G-SSIM for the color 250 dataset in comparison with other decolorization methods.

\begin{figure}[t]
	\centering
	\includegraphics[width = \linewidth]{SSIM2.pdf}
	\caption{Comparison of six color to gray methods based on the C2G-SSIM score.}
	\label{fig:ssim}
\end{figure}

\begin{figure*}[t]
	\centering
	\includegraphics[width = \textwidth]{teaser2.pdf}
	\caption{Chromatic aberration which results from differential refraction of light depending on its wavelength, it causes some rays (green) to converge before other (red). This results in a perception of red ``advancing'' toward the eye, while green to be ``receding''.}
	\label{fig:t2}
\end{figure*}

\begin{figure*}[!t]
	\centering
	\includegraphics[width = \textwidth]{Histogram.pdf}
	\caption{Global tone mapping performed in three different luminance spaces: Our's, Y of YCbCr and V of HSV}
	\label{fig:GTMO}
\end{figure*}

\begin{table}[!t]
\scriptsize
	\renewcommand{\arraystretch}{1.2}
	\caption{Comparison of six color to gray methods based on average objective metrics for 250 images.}
		\centering
	\begin{tabular}{| l | l | l | l | l | l |}
		\hline
		& Decolorization & CCFR & CCPR & E-score  & C2G-SSIM \\ \hline
		(a) & Lu et al. \cite{lu2014contrast} & 0.9922 & 0.9645 & 0.9781 & 0.8900 \\ \hline
		(b) & Nafchi et al. \cite{nafchi2017corrc2g} & 0.9880 & 0.9555 & 0.9710 & 0.8900 \\ \hline
		(c) & Grundland et al \cite{grundland2007decolorize} & 0.9811 & 0.9584 & 0.9747 & 0.8749 \\ \hline
		(d) & Smith et al. \cite{smith2008apparent} & 0.9880 & 0.9555 & 0.9710 & 0.9062 \\ \hline
		(e) & Gooch et al. \cite{gooch2005color2gray} & 0.9839 & 0.9545 & 0.9714 & 0.8935 \\ \hline
		(f) & Ours & 0.9890 & 0.9576 & 0.9728 & 0.9018 \\ \hline
	\end{tabular}
	\label{tab_om}
\end{table}

In our experiment we computed the average CCPR for the 250 images in the dataset by varying $\tau$ from 1 to 15 \cite{lu2014contrast}. As can be seen from Fig. \ref{fig:plotE}, our algorithm's performance is reasonable and practicable when compared to other color to gray algorithms. Figure \ref{fig:ssim} shows the average C2G-SSIM score for the six decolorization methods. According to the plot in Fig. \ref{fig:plotE}, Lu et al.'s method shows best performance based on the E-score, however, our method delivers high average C2G-SSIM measure which is better correlated to human perception \cite{nafchi2017corrc2g}.

\section{Characteristics of visualization of our algorithm}
\subsection{Color temperature and chromatic aberration}\label{ca}
From Snell's law we know that the refraction of light is dependent on its wavelength. As the frequency of light increases, its refractive index becomes larger, causing more refraction of the shorter wavelengths. Therefore, when an image is captured through a lens, all colors do not focus at the same distance, and these imperfections are known as chromatic aberration. In cameras this imperfection is removed by using a combination of second achromatic lens which is made of different material(glass) than the first lens. This second lens would reverse the color dispersion caused by the first lens. The human eye, also employs a lens and does exhibit this phenomena as shown in Fig. \ref{fig:t2}. From this figure we can observe that red light forms the image farthest from the lens as it has the smallest refractive index. Colors with higher index of refraction would ideally bend more thereby forming images closer to the lens. Therefore, it would be impossible to focus on all colors simultaneously, resulting in ``somewhat fuzzy'' images that are not in focus. Colors that are closer to red end of the electromagnetic spectrum are said to be warm colors and are perceived as closer to the observer \cite{niemeyer1988ec88,eastlake2020goethe}. Colors that are around the blue end are said to be cooler colors and are perceived to be receding away from the observer. This phenomenon has been exploited by traditional artists to add depth information in artwork \cite{nguyen2005detection}, display devices \cite{hodges1988chromostereoscopic}, and 3D imagery \cite{hong2012depth}.

\begin{figure*}[!t]
	\centering
	\includegraphics[width = \textwidth]{TMOApplication}
	\caption{Application: (a) Ideal image processing pipeline (tone mapping) (b) Color boosted by local tone mapping and corresponding grayscale image. (c) Demonstration example for easy control features using our decolorization algorithm. TM1: background boost TM2: background suppression}
	\label{fig:SLHE_TMO}
\end{figure*}

\subsection{Tone mapping application}
The color to grayscale conversion is a dimensionality reduction problem whose significance has been underestimated. Usually the tone mapping is performed on the grayscale image because of the lower computational and memory requirements, when compared to tone mapping on the RGB channels. In this section we will discuss global and local tone mapping application using the proposed decolorization algorithm and demonstrate its effectiveness.
\subsubsection{Global tone mapping}

We discuss our decolorization method's applicability for post-processing tasks such as tone mapping. Conventional methods have not shown such applicability because of their long run-time; thus, the use of $Y$ of YCbCr or $V$ of HSV for tone mapping is popular today. We evaluated our method's output of global tone mapping because it has the ability to adequately remap colors to $1-D$ luminance. Thus, even using global tone mapping confirms the effectiveness of our method. Figure \ref{fig:GTMO} shows tone mapped results using $Y$, $V$, or our luminance channel. In the fish image, using $Y$, the fish color became the same as the background and could not be controlled separately, because the weight of $R$ is small in the YCbCr color channel. Our method was able to separate fish color from background in the same way as human perception and directly mapped the color of the background to dark. In the sky and snow scene, the global tone curve for contrast enhancement is a centered sigmoid curve, in which highlight/shadow clippings occur naturally. Since HSV color space treats primary colors and whites the same way, colors of sky and snow were mapped to luminance that was too light. Thus, the sigmoid function degraded the contrast of the outputted image. YCbCr treats $B$ as too low luminance; the sky became too dark. We confirmed that our method generates well-balanced images maintaining the contrast and the colors.

\subsubsection{Local tone mapping}
Local histogram equalization based local tone mapping converts target pixels by using tone curves constructed from local cumulative histograms. Smoothed local histogram equalization (LHE) are also used as a smoothed LH filter \cite{kass2010smoothed}. The Apical's (ARM\textsuperscript{\textregistered}) Iridix algorithm \cite{chesnokov2007image}, which is based on smoothed LHE are used by a range of camera makers, including Nikon, Olympus and Sony. For our local tonemap application we selected smoothed LHE-based function as they are very suitable for practical real-time applications, in Fig. \ref{fig:SLHE_TMO}(a) we show an ideal pipeline for such a system. We implemented the local tonemap application similar to one presented by Ambalathankandy et al. \cite{ambalathankandy2019adaptive}, as their implementation has a linear $\textit{O}(1)$ computational complexity and produces output images with fine quality as shown in Fig. \ref{fig:SLHE_TMO}(b, c). The total computational time including our proposed decolorization and the local tonemap operation was only $14.7\mu s$ per pixel (normalized time $@2.7$GHz CPU). This time utilization is $20\times$ less than Lu et al.'s work \cite{lu2012contrast}. Additionally, using our decolorization method has an advantage, which is to individually control the back/foreground as shown in Fig. \ref{fig:SLHE_TMO}(c).

\subsection{Limitation of our method}
In our decolorization method, pure green is likely to be mapped to a dark luminance value. Light green will also be mapped to a lesser dark luminance part as shown in Fig. \ref{fig:4}. Therefore, certain scenes are likely to be perceived as unnatural. For example:
(a) Vegetables (e.g. leafy greens like cabbage) (b) Green meadows under bright sky. (c) Bright green neon lights. However, in Fig. \ref{fig:Limit} we perceive them as natural. We postulate the following as the possible reasons: (i) There are rarely any pure bright green
(like G255) scene in nature. (ii) In our color space, green color with white components follow Eq.(\ref{eq:Lw1}) by weighting function. Therefore, the color keeps a balance among other color channels. (iii) Vegetation scenery with dark green are well perceived as healthy plants.

\begin{figure*}[!t]
	\centering
	\includegraphics[width = \textwidth]{limitC2G.pdf}
	\caption{Examples to demonstrate limitation of our decolorization method when mapping bright green/neon color.}
	\label{fig:Limit}
\end{figure*}

\section{Conclusion}

In this paper, we present a color temperature-based RGB to gray conversion model by taking in to account the chromatic aberration phenomena. This anomaly results from differential refraction of light depending on its wavelength, it causes some of the rays (cool colors) to converge before others (warm colors). This results in a perception of warmer colors ``advancing'' towards the eye, while the cooler ones to be ``receding''. Essentially, since decolorization is expected to have a key role in the pre-processing of tone mapping or edge preserving filters, low calculation cost and fast operation for the processing are required. To address this requirement, we have developed a high-speed $\textit{O}(1)$ decolorization method that is based on color temperature-based perception. It refers to RGB values in one pixel and performs weighted blending of the Euclidean distances of warm/cool color vectors.
This simple conversion outputs a gray channel that is comparable to the conventional optimization methods using iterations. When our method is applied to tone mapping, it achieves better results than one could be obtained with YCbCr/HSV color space.

\section{Acknowledgments}
We would like to sincerely thank Cadik et al., Lu et al. for making available their datasets, and other authors for sharing their source codes.

\bibliographystyle{plain}


\section{Introduction}
\label{intro}

The presence of heavy flavors in cold quark matter might have non-trivial effects on matter subject to extreme conditions of density. For instance, the case of charm quarks immediately leads to the possibility of charm stars, i.e. self-bound strange stars with some amount of charm quarks in their cores. However, in spite of some constraints on their existence \cite{Lattimer:2010uk,Lattimer:2019eez}, they might still occur in ultradense hybrid neutron stars.
Perturbative quantum chromodynamics (pQCD) can be used much more confidently at these very high densities \cite{Kapusta:2006pm} to build the equation of state (EoS) for charm matter\footnote{This complements advancements in pQCD \cite{Laine:2006cp} and lattice QCD at high temperatures \cite{Borsanyi:2016ksw}.} in electroweak equilibrium.

We investigate cold quark matter with heavy quarks using in-medium pQCD. Specifically, we extend the framework developed in Ref. \cite{Kurkela:2009gj} for $N_{f}=N_{l}+1$ flavors, i.e. $N_{l}$ massless quarks and $1$ massive flavor, to include any number of massive flavors and determine the associated EoS. As an application, we discuss charm matter that is in $\beta$ equilibrium and electrically neutral. This allows us to consider the existence charm quark stars from the pQCD viewpoint. In the past, these exotic stars were investigated only within the bag model, being unstable against radial oscillations \cite{Kettner:1994zs}. Using the method developed in Ref. \cite{Gondek:1997fd} for a couple of first-order differential equations for the Lagrangian radial displacements, our calculations indicate that such self-bound stars are unstable.

\section{Cold and dense pQCD with heavy quarks}
\label{sec1:framework}

The calculation of the EoS for cold quark matter within pQCD has a long history \cite{Kapusta:2006pm}. The three-loop result of Ref. \cite{Kurkela:2009gj} includes renormalization-group effects on the mass and strong coupling parameters. As mentioned above, this framework deals with $N_{f}=N_{l}+1$ quark flavors. The corresponding  thermodynamic potential has a massless contribution plus a massive term, together with the mixed vacuum-matter (VM) diagrams and the corresponding ring terms, i.e. $\Omega=\Omega^{m=0}(\vec{\mu})+\Omega^{m}(\tilde{\mu},m)+\Omega^{x}_{{\rm VM}}(\vec{\mu}, m)+\Omega_{{\rm ring}}(\vec{\mu}, \tilde{\mu}, m)$. Here $m$ is the mass, $\tilde{\mu}$ corresponds to the massive flavor chemical potential and $\vec{\mu}~{\equiv}~(\mu_{1}, ...,\mu_{N_{l}})$ is the vector chemical potential for the massless quarks. $\Omega$ also depends on the renormalization scale parameter $\bar{\Lambda}$ and  $\alpha_{s}(\bar{\Lambda})=({4\pi}/{\beta_{0}L})\left(1-({2\beta_{1}}/{\beta^{2}_{0}}){\ln{L}}/{L}\right)$, where $\beta_{0}=11-2N_{f}/3$, $\beta_{1}=51-19N_{f}/3$, $L=2\ln\left(\bar{\Lambda}/\Lambda_{\rm \overline{MS}}\right)$ and $\Lambda_{\rm \overline{MS}}$ is the $\overline{\rm MS}$ point scale. Usually the renormalization scale is defined as $\bar{\Lambda}=X\sum_i{\mu_{i}}/N_{f}$, where the sum runs over the active quark flavors and the dimensionless parameter $X$ sits between $1$ and $4$ \cite{Kurkela:2009gj}.

On the other hand, at high densities one expects that light quarks could be accompanied by heavy flavors. In this situation, the total number of flavors is $N_{f}=\sum^{N_{m}}_{i=1}(N_{l}+1)^{(i)}$, where $N_{m}$ is the number of massive quarks. This restricts us to add at least one massless quark for each massive flavor. Notice that further physical conditions are required in order to control when a heavy partner is active. Such conditions affect the renormalization scale $\bar{\Lambda}$, depending on the chosen heavy quark to be added. Thus, $\Omega$ at $\mathcal{O}(\alpha^{2}_{s})$ for $N_{l}$ massless and $N_{m}$ massive quarks is $\bar{\Omega}[{N_{f}}]= \sum^{N_{m}}_{i=1}\left\lbrace\Omega[{N^{(i)}_{l}}]+\Omega[{1^{(i)}}]\right\rbrace$, choosing first the massless flavors when adding a massive one, so that $\Omega[{N^{(i)}_{l}}]~{\equiv}~(\Omega^{m=0}(\vec{\mu}))^{(i)}$ is the massless contribution and $\Omega[{1^{(i)}}]~{\equiv}~(\Omega^{m}+\Omega^{x}_{{\rm VM}}+\Omega_{{\rm ring}})^{(i)}$ is the mixed massive contribution. Here $\vec{\mu}^{(i)}=(\mu_{1}, ...,\mu_{i})$ is the massless vector chemical potential, $\tilde{\mu}^{(i)}$ the massive (heavy) quark chemical potentials and $m^{(i)}$ their corresponding masses.  $\Omega[...]$ indicates just the implicit parameter dependence (e.g. on $N_{f}$), whereas $\Omega(...)$ represents an explicit function dependence.

\section{Charm matter thermodynamics}
\label{sec2:charm}

We now apply this formalism to charm quark matter. Below the charm threshold, the usual strange matter conditions apply \cite{Glendenning:2000}. Above the threshold, the reactions ${u+d}\leftrightarrow{c+d}$ yield $\mu_{c}=\mu_{u}$. Furthermore, electric neutrality implies $({2}/{3})n_{u}+({2}/{3})n_{c}-({1}/{3})n_{d}-({1}/{3})n_{s}-n_{e}-n_{\mu}=0$, with $\mu_{\mu}=\mu_{e}$, so that all chemical potentials are functions of $\mu_{s}$. The strange and charm quark mass parameters (up to $\mathcal{O}(\alpha^{2}_{s})$) are \cite{Vermaseren:1997fq} $m_{s}(\bar{\Lambda})=\hat{m}_{s}\left({\alpha_{s}}/{\pi}\right)^{4/9}\left(1+0.895062\left({\alpha_{s}}/{\pi}\right) +1.37143\left({\alpha_{s}}/{\pi}\right)^{2}\right)$ and $m_{c}(\bar{\Lambda})=\hat{m}_{c}\left({\alpha_{s}}/{\pi}\right)^{12/25}\left(1+1.01413\left({\alpha_{s}}/{\pi}\right)+1.38921\left({\alpha_{s}}/{\pi}\right)^{2}\right)$, respectively, $\lbrace\hat{m}_{q}\rbrace$ being the invariant quark masses. By noting that $\alpha_{s}$ depends on $N_{f}$, this gives different values of $\Lambda_{\overline{\rm MS}}$. Choosing $\alpha_{s}(\bar{\Lambda}=1.5~{\rm GeV},~N_{f}=3,4)=0.336^{+0.012}_{-0.008}$ \cite{Bazavov:2014soa} gives $\Lambda^{2+1}_{\overline{\rm MS}}=343^{+18}_{-12}~$MeV and $\Lambda^{2+1+1}_{\overline{\rm MS}}=290^{+18}_{-12}~$MeV, thus defining $\alpha^{2+1}_{s}(\bar{\Lambda})$ and $\alpha^{2+1+1}_{s}(\bar{\Lambda})$, respectively. Fixing the strange quark mass at $m_{s}(2~{\rm GeV}, N_{f}=3,4)=92.4(1.5)~$MeV \cite{Chakraborty:2014aca} gives $\hat{m}^{2+1}_{s}~{\approx}~246.2~$MeV when using $\alpha^{2+1}_{s}$, and $\hat{m}^{2+1+1}_{s}~{\approx}~243.7~$MeV with $\alpha^{2+1+1}_{s}$. Additionally, fixing the charm quark mass at $m_{c}(3~ {\rm GeV}, N_{f}=4)=0.9851(63)~{\rm GeV}{~\equiv~}m^{0}_{c}$ \cite{Chakraborty:2014aca}, gives $\hat{m}^{2+1+1}_{c}~{\approx}~3.0895~$GeV when using $\alpha^{2+1+1}_{s}$, being $m^{0}_{c}$ the vacuum charm mass. Furthermore, charm quarks must satisfy $\mu_{c}=\mu_{s}-\mu_{e}>m^{\rm medium}_{c}>m^{0}_{c}$ being $m^{\rm medium}_{c}$ the unknown in-medium charm mass. The renormalization scale parameter below and above the charm threshold is ${\bar{\Lambda}=}X{(\mu_{u}+\mu_{d}+\mu_{s}+0)}/{3}$ for $\mu_{s}~{\lesssim}~m^{0}_{c}$ and $X^{*}{(\mu_{u}+\mu_{d}+\mu_{s}+\mu_{c})}/{3}$ for $\mu_{s}~{\gtrsim}~m^{0}_{c}$,
where the approximations in the inequalities come from noticing that before the threshold point the electron chemical potential takes its lowest value compared to the strange one, thus allowing us to make the approximation $\mu_{c}~{\approx}~\mu_{s}$.
\begin{figure*}[t!]
\hbox{
\includegraphics[width=0.5\textwidth]{particledensitiesX3.pdf} 	\includegraphics[width=0.5\textwidth]{NormCharmPressureNf4x2.pdf}}
\caption{Left panel: Quark and lepton populations, $n_{i}/n$, present in charm matter for $X=3$. Above the charm threshold, the lepton fractions increase to ensure electric charge neutrality at high densities. Right panel: Total pressure of $N_{f}=2+1+1$ quarks plus leptons in $\beta$ equilibrium and electrically charge neutral normalized by the Stefan-Boltzmann gas with $N_{f}=4$. The case $X=5$ is shown only to verify how the EoS depends on $X$ when including one additional heavy flavor. Taken from Ref. \cite{Jimenez:2019kji}.}
\label{fig:thermo}
\end{figure*}
To build the EoS, we choose $N^{(1)}_{l}=1$ for the up, $N^{(2)}_{l}=1$ for the down, and $N_{m}=2$ for the strange and charm quarks in the thermodynamic potential $\bar{\Omega}[N_{f}=2+1+1]=\left\lbrace\Omega[N^{(1)}_{l}=1]+\Omega[1^{(1)}]\right\rbrace+\left\lbrace\Omega[N^{(2)}_{l}=1]+\Omega[1^{(2)}]\right\rbrace$, so that the flavors are counted as $N_{f}=(1+1)^{(1)}+(1+1)^{(2)}=(u+c)^{(1)}+(d+s)^{(2)}$. We define the total quark number density for charm matter, for a given $X$, as $n_{q}~{\equiv}~n_{u}+n_{d}+n_{s}+n_{c}$ and the total particle density as $n=n_{q}+n_{L}$, where $n_{L}=n_{e}+n_{\mu}$. In Fig. \ref{fig:thermo}, we show the behavior of the relative particle populations in $\beta$ equilibrium and electrically neutral charm quark matter in the case of $X=3$. The total pressure and energy density are obtained by keeping thermodynamic consistency systematically. This allows us to build numerically the EoS, $P=P(\epsilon, X)$. In Fig. \ref{fig:thermo}, the total pressure is displayed for charm matter normalized by a Stefan-Boltzmann gas of quarks with $N_{f}=4$ as a function of the strange quark chemical potential. In the same figure, we show that the charm EoS is largely softened, generating an apparent instability which could have astrophysical effects, e.g. the possibility of charm stars.
\begin{figure*}[t!]
\hbox{
\includegraphics[width=0.5\textwidth]{MassCentralDensNf4x.pdf} 	\includegraphics[width=0.5\textwidth]{MassRadiiNf4.pdf}}
\caption{Left panel: Gravitational masses vs central energy densities for a system with $N_{f}=2+1+1$ quarks plus leptons. Right panel: Mass-radius diagram for the associated quark stars. The black dots indicate the appearance of the charm quarks. Stars not satisfying the general-relativistic (GR, magenta), causality (green), and finite pressure ($P<\infty$, light blue) limits are excluded \cite{Lattimer:2006xb}. Taken from Ref. \cite{Jimenez:2019kji}.}
\label{fig:diagrams}
\end{figure*}
\section{Results and Conclusion}
\label{sec3:results}

The possibility of bare charm quark stars was studied in Ref. \cite{Kettner:1994zs} using the simplest version of the bag model. Their conclusion was that charm stars would be unstable against radial perturbations. We revisit this question using our first-principle pQCD results for the EoS for charm quark matter. We also choose the range $X\geq3$, which satisfies the Bodmer-Witten hypothesis of strange quark matter. Although the structure equations ensure the relativistic hydrostatic equilibrium of the stars, they must also satisfy the condition $\partial{M}/\partial{\epsilon_{c}}~{\geq}~0$ \cite{Glendenning:2000}. The maximal masses occur at $\partial{M}/\partial{\epsilon_{c}}=0$ (see Fig. \ref{fig:diagrams}). For the case $X=2$, this condition is not satisfied when charm quarks appear. On the other hand, for $X>3$, this condition is satisfied when charm quarks are present, which would correspond to stable charm stars. In Fig. \ref{fig:diagrams}, we show the mass-radius diagram for quark stars made of $N_{f}=2+1+1$ quarks plus electrons and muons with varying $X$. Performing the dynamical stability analysis against radial pulsations using the method of Ref. \cite{Gondek:1997fd}, amplitudes oscillate harmonically when the frequencies are such that ${\rm Re(\omega_{n})>0}$ and ${\rm Im(\omega_{n})=0}$, or increase exponentially if ${\rm Re(\omega_{n})~{\geq}~0}$ and ${\rm Im(\omega_{n})>0}$.  Besides, if ${\rm Im(\omega_{0})>0}$ [and ${\rm Re(\omega_{0})=0}$] from some value of central energy density $\epsilon_{c}$, then all the higher modes will become complex too, representing the onset of the instability. Our calculations show that for $X=3$ the amplitudes of the stellar configurations increase even in the region where charm stars would be expected, making them dynamically unstable. Since the same behavior was obtained for larger values of $X$, one can conclude from a perturbative QCD analysis that charm stars are unstable. In the future, it would be interesting to study the presence of charm quarks in neutron-star mergers.

\section*{Acknowledgement}

This work was partially supported by INCT-FNA (Process No. 464898/2014-5). J. C. J. acknowledges support from FAPESP (Processes No. 2020/07791-7 and No. 2018/24720-6). E. S. F. is partially supported by CAPES (Finance Code 001), CNPq, and FAPERJ.