\section{Introduction}

Programming robots is very difficult, in part because the real world is inherently rich and---to some degree---unpredictable. In addition, our expectations for physical agents are quite high and often difficult to articulate. Nevertheless, for robots to have a significant impact on the lives of individuals, even non-programmers need to be able to specify and customize behavior.
Because of these complexities, relying on end-users to provide instructions to robots programmatically seems destined to fail.

Reinforcement learning (RL) from human trainer feedback provides a compelling alternative to programming because agents can learn complex behavior from very simple positive and negative signals. Furthermore, real-world animal training is an existence proof that people can train complex behavior using these simple signals. Indeed, animals have been successfully trained to guide the blind, locate mines in the ocean, detect cancer or explosives, and even solve complex, multi-stage puzzles.

However, traditional reinforcement-learning algorithms have yielded limited success when the reward signal is provided by humans and have largely failed to benefit from the sophisticated training strategies that expert animal trainers use with animals. This failure has led to the development of new reinforcement-learning algorithms that are designed to learn from human-generated rewards and investigations into how people give interactive feedback~\citep{knox2009interactively,thomaz2006reinforcement,thomaz2007robot,thomaz08,griffith2013policy,loftin15}. In general, many of these human-centered learning algorithms are built on the insight that people tend to give feedback that reflects the policy the agent should be following, rather than as a numeric value that is meant to be maximized by the agent. While this insight seems accurate, existing approaches assume models of feedback that are independent of the policy the agent is currently following. We present empirical results that demonstrate that this assumption is incorrect and further argue that policy-dependent feedback enables effective training strategies, such as
differential feedback and policy shaping, from which we would like a learning agent to benefit. Following this result, we present {\em Convergent Actor-Critic by Humans} (COACH), an algorithm for learning from policy-dependent human feedback that is capable of benefiting from these strategies. COACH is based on the insight that the TD-error used by actor-critic algorithms is an unbiased estimate of the {\em advantage function}, which is a policy-dependent value roughly corresponding to how much better or worse an action is compared to the current policy and which captures these previously mentioned training strategies. To validate that COACH scales to complex problems, we train five different behaviors on a TurtleBot robot that makes decisions every 33ms from noisy image features that are not visible to the trainer.

\section{Background}

For modeling the underlying decision-making problem of an agent being taught by a human, we adopt the Markov Decision Process (MDP) formalism. An MDP is a 5-tuple: $\langle S, A ,T, R, \gamma \rangle$, where $S$ is the set of possible states of the environment; $A$ is the set of actions available to the agent; $T(s' | s, a)$ is the transition function, which defines the probability of the environment transitioning to state $s'$ when the agent takes action $a$ in environment state $s$; $R(s, a, s')$ is the reward function specifying the numeric reward the agent receives for taking action $a$ in state $s$ and transitioning to state $s'$; and $\gamma \in [0, 1]$ is a discount factor specifying how much immediate rewards are preferred to more distant rewards.

A {\em stochastic policy} $\pi$ for an MDP is a per-state action probability distribution that defines a mode of behavior; $\pi : S \times A \rightarrow [0, 1]$, where $\sum_{a \in A} \pi(s,a) = 1, \forall s \in S$. In the MDP setting, the goal is to find the optimal policy $\pi^*$, which maximizes the expected future discounted reward when the agent selects actions in each state according to $\pi^*$; $\pi^* = \argmax_\pi E[\sum_{t=0}^\infty \gamma^t r_t | \pi]$, where $r_t$ is the reward received at time $t$. Two important concepts in MDPs are the value function ($V^\pi$) and action--value function ($Q^\pi$). The value function defines the expected future discounted reward from each state when following some policy
and the action--value function defines the expected future discounted reward when an agent takes some action in some state and then follows some policy $\pi$ thereafter.
These equations can be recursively defined via the Bellman equation: $V^\pi(s) = \sum_a \pi(s,a) Q^\pi(s, a)$ and $Q^\pi(s, a) = \sum_{s'} T(s' | s, a) \left[ R(s, a, s') + \gamma V^\pi(s')\right]$. For shorthand, the value functions for the optimal policies are usually denoted $V^*$ and $Q^*$.

In reinforcement learning (RL), an agent interacts with an environment modeled as an MDP, but does not have direct access to the transition function or reward function and instead must learn a policy from environment observations. A common class of RL algorithms are {\em actor-critic} algorithms. Bhatnagar et al.~\citep{bhatnagar2009natural} provides a general template for these algorithms. Actor-critic algorithms are named for the two main components of the algorithms: The actor is a parameterized policy that dictates how the agent selects actions; the critic estimates the value function for the actor and provides critiques at each time step that are used to update the policy parameters. Typically, the critique
is the temporal difference (TD) error: $\delta_t = r_t + \gamma V(s_t) - V(s_{t-1})$, which describes how much better or worse a transition went than expected.

\section{Human-centered Reinforcement Learning}
In this work, a {\em human-centered reinforcement-learning} (HCRL) problem is a learning problem in which an agent is situated in an environment described by an MDP but in which rewards are generated by a human trainer instead of from a stationary MDP reward function that the agent is meant to maximize.  The trainer has a target policy $\pi^*$ they are trying to teach the agent. The trainer communicates this policy by giving numeric feedback as the agent acts in the environment.
The goal of the agent is to learn the target policy $\pi^*$ from the feedback.

To define a learning algorithm for this problem, we first characterize how human trainers typically use numeric feedback to teach target policies. If feedback is stationary and intended to be maximized, it can be treated as a reward function and standard RL algorithms used.
Although this approach has had some success~\citep{pilarski2011online,isbell2001social}, there are complications that limit its applicability. In particular, a trainer must take care that the feedback they give contains no unanticipated exploits, constraining the feedback strategies they can use. Indeed, prior research has shown that interpreting human feedback like a reward function often induces {\em positive reward cycles} that leads to unintended behaviors~\citep{knoxthesis,ho15}.

The issues with interpreting feedback as reward have led to the insight that human feedback is better interpreted as a comment on the agent's behavior; for example, positive feedback roughly corresponds to ``that was good'' and negative feedback roughly corresponds to ``that was bad.'' Existing HCRL work adopting this perspective includes TAMER~\citep{knox2009interactively}, SABL~\citep{loftin15}, and Policy Shaping~\citep{griffith2013policy}, discussed in more detail in the Related Work section. We note, though, that all assume that human feedback is independent of the agent's current policy. We provide empirical results that show this assumption to be incorrect

\section{Policy-dependent Feedback}

Evidence that human feedback is influenced by the agent's current policy can be seen in previous work. For example, it was observed that trainers taper their feedback over the course of learning~\citep{ho15,knox2012humans,isbell2001social}.
One explanation for decreasing feedback rates is policy-dependent feedback, but trainer fatigue is another.
We provide a stronger result showing that
trainers---for the same state--action pair---choose positive or negative feedback depending on their perception of the learner's behavior. This finding serves as a warning for algorithms that rely on an assumption
policy-independent feedback.

\subsection{Empirical Results}
\begin{figure}
\centering
\includegraphics[width=0.8\columnwidth]{interface.png}
\caption{The training interface shown to AMT users.}
\label{fig:interface}
\end{figure}

We had Amazon Mechanical Turk (AMT) participants teach an agent in a simple sequential task, illustrated
in Figure~\ref{fig:interface}. Participants were instructed to train a virtual dog to walk to the yellow goal location in a grid world as fast as possible but without going through the green cells. They were additionally told that, as a result of prior training, their dog was already either ``bad'', ``alright'', or ``good'' at the task and were shown examples of each behavior before training. In all cases, the dog would start in the location shown in Figure~\ref{fig:interface}. ``Bad'' dogs walked straight through the green cells to the yellow cell. ``Alright'' dogs first moved left, then up, and then to the goal, avoiding green but not taking the shortest route. ``Good'' dogs took the shortest path to yellow without going through green.

During training, participants saw the dog take an action from one tile to another and then gave feedback after every action using a continuous labeled slider as shown. The slider always started in the middle of the scale on each trial, and several points were labeled with different levels of reward (praise and treats) and punishment (scolding and a mild electric shock). Participants went through a brief tutorial using this interface. Responses were coded as a numeric value from $-50$ to $50$, with ``Do Nothing'' as the zero-point.

During the training phase, participants trained a dog for three {\em episodes} that all started in the same position and ended at the goal. The dog's behavior was pre-programmed for all episodes in such a way that the first step of the final episode would reveal if feedback was policy-dependent. The dog always performed the same behavior in the first two episodes, and then performed the ``alright'' behavior in the third episode. Each user was placed into one of three different conditions that determined how the dog would behave in the first two episodes: either ``bad,'' ``alright,'' or ``good.'' If  feedback is policy dependent, we expect more positive feedback in the ``bad'' condition than in the ``alright'' or ``good'' condition; ``alright'' behavior is an improvement over the previous ``bad'' behavior whereas it is either no improvement or a deterioration compared to ``alright'' or ``good'' behavior.

Figure~\ref{fig:feedDist} shows boxplots and individual responses for the first step of the final episode under each of the three conditions. These results indicate that the sign of feedback is sensitive to the learner's policy, as predicted. The mean and median feedback under the ``bad'' condition is slightly positive (Mean = $9.8$, Median = $24$, S.D. = $22.2$; planned Wilcoxon one-sided signed-rank test: $Z = 1.71, p < 0.05$), whereas it is negative for the ``alright'' condition (Mean = $-18.3$, Median = $-23.5$, S.D. = $24.6$; planned Wilcoxon two-sided signed-rank test: $Z = -3.15, p < 0.01$) and ``good'' condition (Mean = $-10.8$, Median = $-18.0$, S.D. = $20.7$; planned Wilcoxon one-sided signed-rank test: $Z = -2.33, p < 0.05$). There was a main effect across the three conditions ($p < 0.01$, Kruskal-Wallace Test), and pairwise comparisons indicated that only the ``bad'' condition differed from ``alright'' and ``good'' conditions ($p < 0.01$ for both, Bonferroni-corrected, Mann-Whitney Pairwise test).
\begin{figure}
\centering
\includegraphics[width=0.8\columnwidth]{feedback_dist.png}
\caption{Box plots with median, interquartile range, and minimum/maximum values of human responses to first step of the final episode for each condition. Feedback tended to be positive when the prior behavior was \textit{bad} and negative otherwise.}
\label{fig:feedDist}
\end{figure}

\subsection{Training Strategies}
Beyond the fact that our previous evidence suggests that people give policy-dependent feedback, we argue that policy-dependent feedback affords desirable training strategies. Specifically, we consider three different feedback schemes that can be viewed as
operationalizations of well-studied behavior analysis reinforcement schedules~\citep{miltenberger2011behavior}: {\em Diminishing Returns}: gradually decrease positive feedback for good actions as the agent adopts those actions. {\em Differential Feedback}: vary the magnitude of feedbacks w.r.t.\ the degree of improvement or deterioration in behavior. {\em Policy Shaping}: provide positive feedback for suboptimal actions that improve behavior and then negative feedback after the improvement has been made.

Diminishing returns
is a useful strategy because it decreases the burden of how actively a trainer must supply feedback and removes the need for explicit training and execution phases.
Differential feedback is useful because it can serve to highlight the most important behaviors in the state space and communicate a kind of urgency for learning them.
Finally, policy shaping concerns feedback that signals an improvement relative to the current baseline---as in the AMT study above. It is useful for providing a direction for the learner to follow at all times, even when the space of policies is continuous or otherwise impractical to search.

\section{Convergent Actor-Critic by Humans}

In this section, we introduce {\em Convergent Actor-Critic by Humans} (COACH), an actor-critic-based algorithm capable of learning from policy-dependent feedback. COACH is based on the insight that the {\em advantage function} is a good model of human feedback and that actor-critic algorithms update a policy using the critic's TD error, which is an unbiased estimate of the advantage function. Consequently, an agent's policy can be directly modified by human feedback without a critic component. We first define the advantage function and describe how it relates to the three previously mentioned training strategies. Then, we present the general update rule for COACH and its convergence. Finally, we present {\em Real-time COACH}, which includes mechanisms for providing variable magnitude feedback and learning in problems with a high-frequency decision cycle.

\subsection{The Advantage Function and Training Strategies}
The advantage function~\citep{baird1995residual} $A^\pi$ is defined as
\begin{equation}
A^\pi(s, a) = Q^\pi(s, a) - V^\pi(s).
\end{equation}
Roughly speaking, the advantage function describes how much better or worse an action selection is compared to the agent's performance following policy $\pi$. We now show that feedback assigned according to the advantage function follows the patterns of all three training strategies. Showing that the advantage function captures the differential feedback strategy is trivial, because the advantage function is defined by how much better taking an action over its current policy is expected to be.

To show that the advantage function induces diminishing returns, consider what happens when the learner improves its behavior by shifting to a higher scoring action $a$ in some state $s$. As its probability of selecting $a$ goes to 1,
$A^\pi(s, a) = Q^\pi(s, a) - Q^\pi(s,a) = 0$, because the value function $V^\pi(s) = \sum_a \pi(s,a) Q^\pi(s,a)$.
Since the expected value is a smooth linear combination of the Q-values, as the agent adopts action $a$, $A^\pi(s, a) \rightarrow 0$, resulting in gradually decreasing feedback.

To show that the advantage function induces policy shaping, let us first assume w.l.o.g.\ that there are three possible actions, where action $a_3$ is exclusively optimal, and action $a_2$ is better than action $a_1$. That is, $Q^*(s, a_3) > Q^*(s, a_2) > Q^*(s, a_1)$. For illustrative purposes, let us also assume that the agent has learned the optimal policy in all states except state $s$, wherein the agent with near-certain probability selects the worst action $a_1$ ($\pi(s, a_1) \approx 1$) and that all actions lead to some other state. In this scenario, $A^\pi(s, a) = Q^*(s, a) - \sum_{a'} \pi(s, a') Q^*(s, a') \approx Q^*(s, a) - Q^*(s, a_1)$.
Since $Q^*(s, a_2) > Q^*(s, a_1)$, it follows that $A^\pi(s, a_2)$ is positive. Consequently, in this condition, suboptimal action $a_2$ would receive positive feedback.

Now consider the case when the agent with near-certainty selects optimal action $a_3$. Under this condition, $A^\pi(s, a) \approx Q^*(s, a) - Q^*(s, a_3)$ and since $Q^*(s, a_3) > Q^*(s, a_2)$, $A^\pi(s, a_2)$ is negative. Once again, since $V^\pi$ is a smooth linear function of $Q^\pi$, as the agent adopts optimal action $a_3$,
$A^\pi(s, a_2)$ becomes negative. Therefore, suboptimal actions can be rewarded and then punished as the agent improves at the task, producing policy shaping.

\subsection{Convergence and Update Rule}

Given a performance metric $\rho$, Sutton et al.~\citep{sutton1999policy} derive a policy gradient algorithm of the form:
$
\Delta \bm{\theta} = \alpha \nabla_{\bm \theta} \rho
$.
Here, ${\bm \theta}$ represents the parameters that control the agent's behavior and $\alpha$ is a learning rate.
Under the assumption that $\rho$ is the discounted expected reward from a fixed start state distribution, they show that
$$
\nabla_{\bm \theta} \rho = \sum_s d^\pi(s) \sum_a \nabla_{\bm \theta} \pi(s,a) Q^\pi(s,a),
$$
where $d^\pi(s)$ is the component of the (discounted) stationary distribution at $s$.
A benefit of this form of the gradient is that, given that states are visited according to $d^\pi(s)$ and actions are taken according to $\pi(s,a)$, the update at time $t$ can be made as:
\begin{equation}
\Delta {\bm \theta}_t = \alpha_t \nabla_{\bm \theta} \pi(s_t,a_t) \frac{f_{t+1}}{\pi(s_t,a_t)},
\label{e:update}
\end{equation}
where $E[f_{t+1}] = Q^\pi(s_t,a_t) - v(s)$ for any action-independent function $v(s)$.

In the context of the present paper, $f_{t+1}$ represents the feedback provided by the trainer. It follows trivially that if the trainer chooses the policy-dependent feedback $f_t = Q^\pi(s_t,a_t)$, we obtain a convergent learning algorithm that (locally) maximizes discounted expected reward. In addition, feedback of the form $f_t = Q^\pi(s_t,a_t)-V^\pi(s_t) = A^\pi(s, a)$ also results in convergence.
Note that for the trainer to provide feedback in the form of $Q^\pi$ or $A^\pi$, they would need to ``peer inside'' the learner and observe its policy. In practice, the trainer estimates $\pi$ by observing the agent's actions.

\subsection{Real-time COACH}

There are challenges in implementing Equation~\ref{e:update} for real-time use in practice. Specifically, the interface for providing variable magnitude feedback needs to be addressed, and the question of how to handle sparseness and the timing of feedback needs to be answered. Here, we introduce {\em Real-time COACH}, shown in Algorithm~\ref{alg:coach}, to address these issues.

For providing variable magnitude reward, we use reward aggregation~\citep{knox2009interactively}. In reward aggregation, a trainer selects from a discrete set of feedback values and further raises or lowers the numeric value by giving multiple feedbacks in succession that are summed together.

While sparse feedback is not especially problematic (because no feedback results in no change in policy), it may slow down learning unless the trainer is provided with a mechanism to allow feedback to affect a history of actions. We use {\em eligibility traces}~\citep{6313077} to help apply feedback to the relevant transitions. An eligibility trace is a vector that keeps track of the policy gradient and decays exponentially with a parameter $\lambda$. Policy parameters are then updated in the direction of the trace, allowing feedback to affect earlier decisions.
However, a trainer may not always want to influence a long history of actions. Consequently, real-time COACH maintains multiple eligibility traces with different temporal decay rates and the trainer chooses which eligibility trace to use. This trace choice may be handled implicitly with the feedback value selection or explicitly.

Due to reaction time, human feedback is typically delayed by about $0.2$ to $0.8$ seconds from the event to which they meant to give feedback~\citep{knoxthesis}. To handle this delay, feedback in Real-time COACH is associated with events from $d$ steps ago to cover the gap. Eligibility traces further smooth the feedback to older events.

Finally, we note that just as there are numerous variants of actor-critic update rules, similar variations can be used in the context of COACH.

\begin{algorithm}[t]
\begin{algorithmic}
\Require{policy $\pi_{\bm{\theta}_0}$, trace set $\bm{\lambda}$, delay $d$, learning rate $\alpha$}
\State Initialize traces $\bm{e}_\lambda \gets \bm{0} \mbox{ } \forall \lambda \in \bm{\lambda}$
\State observe initial state $s_0$
\For{$t=0 \text{ to } \infty$}
\State select and execute action $a_t \sim \pi_{\bm{\theta}_t}(s_t, \cdot)$
\State observe next state $s_{t+1}$, sum feedback $f_{t+1}$, and $\lambda$
\For{$\lambda' \in \bm{\lambda}$}
\State $\bm{e}_{\lambda'} \gets \lambda' \bm{e}_{\lambda'} + \frac{1}{\pi_{\bm{\theta}_t}(s_{t-d}, a_{t-d})}\nabla_{\bm{\theta}_t}\pi_{\bm{\theta}_t}(s_{t-d}, a_{t-d})$
\EndFor
\State $\bm{\theta}_{t+1} \gets \bm{\theta}_t + \alpha f_{t+1} \bm{e}_\lambda$
\EndFor
\end{algorithmic}
\caption{Real-time COACH}
\label{alg:coach}
\end{algorithm}

\section{Related Work}

An inspiration for our work is the TAMER framework~\citep{knoxthesis}.
In TAMER, trainers provide interactive numeric feedback as the learner takes actions. The feedback is interpreted as an exemplar of the reward function for the previous state--action pair and is used to learn the reward function. When the agent makes rapid decisions, TAMER divides the feedback among the recent state--action pairs according to a probability distribution. TAMER makes decisions by myopically choosing the action with the highest reward estimate. Because the agent myopically maximizes reward, the feedback may be interpreted as exemplars of $Q^*$. Later work investigated non-myopically maximizing the learned reward function with a planning algorithm~\citep{knox2013learning}, but this approach requires a model of the environment
and special treatment of termination conditions. Because TAMER expects feedback to be policy independent, it does not support the diminishing returns or policy-shaping strategies, and handles differential feedback only in so far as it uses numeric feedback.
In our robotics case study, we provide an explicit example where TAMER's failure to support diminishing returns can result in unlearning.

Two other closely related approaches are SABL~\citep{loftin15} and Policy Shaping~\citep{griffith2013policy} (unrelated to the policy shaping feedback {\em strategy} defined above). Both of these approaches treat feedback as discrete probabilistic evidence of a parametrized policy. SABL's probabilistic model additionally includes (learnable) parameters for describing how often a trainer is expected to give explicit positive or negative feedback.
Both of these approach assume policy-independent feedback and do not support the three training strategies described.

There have also been some domains in which treating human feedback as reward signals to maximize has had some success, such as in shaping the control for a prosthetic arm~\citep{pilarski2011online} and learning how to interact in an online chat room from multiple users' feedback~\citep{isbell2001social}. An interesting area of future work is to test whether performance in these domains can be improved with COACH given our insights into the nature of human feedback---work on the chat-room learning domain did report challenges due to diminishing feedback

Some research has examined combining human feedback with more traditional environmental rewards~\citep{AAMAS10-knox,tenorio2010dynamic,clouse1992teaching,maclin2005giving}. A challenge in this context in practice is that rewards do not naturally come from the environment and must be programmatically defined. However, it is appealing because the agent can learn in the absence of an active trainer. We believe that COACH could be straightforwardly modified to learn in this setting as well.

Finally, a related research area is {\em learning from demonstration} (LfD), in which a human provides examples of the desired behavior. There are a number of different approaches to solving this problem surveyed by Argall et al.~\citep{argall2009survey}. We see these approaches as complementary to HCRL because it is not always possible, or convenient, to provide demonstrations. LfD approaches that learn a parametrized policy could also operate with COACH, allowing the agent to have their policy seeded by demonstrations, and then fine tuned with interactive feedback.

\section{Robotics Case Study}

In this section, we present qualitative results applying Real-time COACH to a TurtleBot robot. The goal of this study was to test that COACH can scale to a complex domain involving multiple challenges, including training an agent that operates on a fast decision cycle (33ms), noisy non-Markov observations from a camera, and agent perception that is hidden from the trainer. To demonstrate the flexibility of COACH, we trained it to perform five different behaviors involving a pink ball and cylinder with an orange top using the same parameter selections. We discuss these behaviors below.  We also contrast the results to training with TAMER. We chose TAMER as a comparison because, to our knowledge, it is the only HCRL algorithm with success on a similar platform~\citep{knox2013training}.

The TurtleBot is a mobile base with two degrees of freedom that senses the world from a Kinect camera. We discretized the action space to five actions: forward, backward, rotate clockwise, rotate counterclockwise, and do nothing. The agent selects one of these actions every 33ms. To deliver feedback, we used a Nintendo Wii controller to give $+1$, $+4$, or $-1$ numeric feedback, and pause and continue training. For perception, we used only the RGB image channels from the Kinect. Because our behaviors were based around a relocatable pink ball and a fixed cylinder with an orange top, we hand constructed relevant image features to be used by the learning algorithms. These features were generated using techniques similar to those used in neural network architectures. In the future, we will investigate learning these features along with the policy. The features were constructed by first transforming the image into two color channels associated with the color of the ball and cylinder. Sum pooling to form a lower-dimensional $8\times 8$ grid was applied to each color channel. Each sum-pooling unit was then passed through three different normalized threshold units defined by $T_i(x) = \min(\frac{x}{\phi_i},1)$, where $\phi_i$ specifies how quickly the $i$th threshold unit saturates. Using multiple saturation parameters differentiates the distance of objects, resulting in three ``depth'' scales per color channel.
Finally, we passed these results through a $2\times 8$ max-pooling layer with stride 1.

The five behaviors we trained were push--pull, hide, ball following, alternate, and cylinder navigation. In push--pull, the TurtleBot is trained to navigate toward the ball when it is far, and back away from it when it is near. The hide behavior has the TurtleBot back away from the ball when it is near and turn away from it when it is far. In ball following, the TurtleBot is trained to navigate to the ball. In the alternate task, the TurtleBot is trained to go back and forth between the cylinder and ball.
Finally, cylinder navigation involves the agent navigating to the cylinder. We further classify training methods for each of these behaviors as {\em flat}, involving the push--pull, hide, and ball following behaviors; and {\em compositional}, involving the alternate and cylinder navigation behaviors.

In all cases, our human trainer (one of the co-authors) used differential feedback and diminishing returns to quickly reinforce behaviors and restrict focus to the areas needing tuning. However,
in alternate and cylinder navigation, they attempted more advanced compositional training methods.
For alternate, the agent was first trained to navigate to the ball when it sees it, and then turn away when it is near. Then, the same was independently done for the cylinder. After training, introducing both objects would cause the agent to move back and forth between them. For cylinder navigation, they attempted to make use of an animal-training method called {\em lure training} in which an animal is first conditioned to follow a lure object that is then used to guide it through more complex behaviors. In cylinder navigation, they first trained the ball to be a lure, used it to guide the TurtleBot to the cylinder, and finally gave a $+4$ reward to reinforce the behaviors it took when following the ball (turning to face the cylinder, moving toward it, and stopping upon reaching it). The agent would then navigate to the cylinder without requiring the ball to be present.

For COACH parameters, we used a softmax parameterized policy, where each action preference value was a linear function of the image features, plus $\tanh(\theta_a)$, where $\theta_a$ is a learnable parameter for action $a$, providing a preference in the absence of any stimulus. We used two eligibility traces with $\lambda=0.95$ for feedback $+1$ and $-1$, and $\lambda=0.9999$ for feedback $+4$. The feedback-action delay $d$ was set to 6, which is $0.198$ seconds. Additionally, we used an actor-critic parameter-update rule variant in which action preference values are directly modified (along its gradient), rather than by the gradient of the policy~\citep{sutton1998reinforcement}. This variant more rapidly communicates stimulus--response preferences. For TAMER, we used typical parameter values for fast decision cycle problems: delay-weighted aggregate TAMER with uniform distribution credit assignment over $0.2$ to $0.8$ seconds, $\epsilon_p=0$, and $c_{min}=1$~\citep{knoxthesis}. (See prior work for parameter meaning.) For TAMER's reward function approximation, we used the same parameters as the action preferences in COACH.

\subsection{Results and Discussion}

COACH was able to successfully learn all five behaviors and a video showing its learning is available online at \url{https://vid.me/3h2s}. Each of these behaviors were trained in less than two minutes, including the time spent verifying that a behavior worked.  Differential feedback and diminishing returns allowed only the behaviors in need of tuning to be quickly reinforced or extinguished without any explicit division between training and testing. Moreover, the agent successfully benefited from the compositional training methods, correctly combining subbehaviors for alternate, and quickly learning cylinder navigation with the lure.

TAMER only successfully learned the behaviors using the flat training methodology and failed to learn the compositionally trained behaviors. In all cases, TAMER tended to forget behavior, requiring feedback for previous decisions it learned to be resupplied after it learned a new decision. For the alternate behavior, this forgetting led to failure: after training the behavior for the cylinder, the agent forgot some of the ball-related behavior and ended up drifting off course when it was time to go to the ball. TAMER also failed to learn from lure training, which was expected since TAMER does not allow reinforcing a long history of behaviors.

We believe TAMER's forgetting is  a
result of interpreting feedback as reward-function exemplars in which new feedback in similar contexts can change the target.
To illustrate this problem, we constructed a well-defined scenario in which TAMER consistently unlearns behavior. In this scenario, the goal was for the TurtleBot to always stay whenever the ball was present, and move forward if just the cylinder was present. We first trained TAMER to stay when the ball alone was present using many rapid rewards (yielding a large aggregated signal). Next, we trained it to move forward when the cylinder alone was present. We then introduced both objects, and the TurtleBot correctly stayed. After rewarding it for staying with a single reward (weaker than the previously-used many rapid rewards), the TurtleBot moved forward. This counter-intuitive unlearning is a consequence of the small reward decreasing its reward-function target for the stay action to a point lower than the value for moving forward. COACH does not exhibit this problem---any reward for staying will strengthen the behavior.

\section{Conclusion}

In this work, we presented empirical results that show that the numeric feedback people give agents in an interactive training paradigm is influenced by the agent's current policy and argued why such policy-dependent feedback enables useful training strategies. We then introduced COACH, an algorithm that, unlike existing human-centered reinforcement-learning algorithms, converges to a local optimum when trained with policy-dependent feedback. We finally we showed that COACH scales up in the context of a robotics case study in which a TurtleBot was successfully taught multiple behaviors with advanced training methods.

There are a number of exciting future directions to extend this work. In particular, because COACH is built on the actor-critic paradigm, it should be possible to combine it straightforwardly with learning from demonstration and environmental rewards, allowing an agent to be trained in a variety of ways. Second, because people give policy-dependent feedback, greater gains may be possible by investigating how people model the current policy of the agent and how their model differs from the agent's actual policy.

\bibliographystyle{IEEEtranN}


\section{Introduction}
\label{sec:intro}

These proceedings contain a status update of an ongoing calculation of $D^0$-$\bar{D}^0$ mixing matrix
elements~\cite{Chang:2014cea}, similar to our published work on $B^0$-$\bar{B}^0$
mixing~\cite{Bazavov:2016nty}.
We present nearly final results for all five matrix elements, sufficient to describe $D^0$-$\bar{D}^0$
mixing not only in the Standard Model, but also in any high-energy extension that modifies only the local
$\Delta C=2$ interaction.

In the Standard Model, neutral-meson mixing is mediated by one-loop, GIM-suppressed processes, shown in
Fig.~\ref{fig:box}.
\begin{figure}[b]
    \centering
    \vspace*{-1em}
    \includegraphics[width=0.7\textwidth]{figs/box}
    \vspace*{-1em}
    \caption{Box diagrams mediating  $D^0$-$\bar{D}^0$ mixing in the Standard Model.}
    \label{fig:box}
\end{figure}
In extensions of the Standard Model, other particles could appear in the boxes; there could even be
tree-level flavor-changing neutral currents.
Mixing has been observed in all four neutral-meson systems---$K^0$, $D^0$, $B^0$, and $B^0_s$---but the
pattern of internal quark masses and CKM\ factors explains why the phenomenology
differs so greatly from one system to another.

Because the $W$ bosons and $b$ quarks have masses well above the QCD scale, mixing can be
re-expressed as stemming both from a local $\Delta C=2$ interaction and two $\Delta C=1$ interactions
separated by a distance of order $1/\Lambda_\text{QCD}$.
From degenerate perturbation theory, the off-diagonal term in the mass-width matrix is~\cite{Artuso:2008vf}
\begin{equation}
    M_{12} - \frac{i}{2}\Gamma_{12} \propto \langle D^0|\mathcal{L}^{\Delta C=2}|\bar{D}^0\rangle +
        \sum_n \frac{\langle D^0|\mathcal{L}^{\Delta C=1}|n\rangle\langle n|\mathcal{L}^{\Delta C=1}|\bar{D}^0\rangle}{M_D-E_n+i0^+}.
    \label{eq:mix}
\end{equation}
The second term is very difficult to estimate.
For $D^0$ mesons it is also not negligible, unlike for $B^0$ and $B^0_s$, where $t$, $c$, and~$u$ quarks
appear in the box.
(For kaons, the second term is important but not dominant.) \@ %
One can relate the measured mass and width differences, $\Delta M$ and $\Delta\Gamma$, to $|M_{12}|$,
$|\Gamma_{12}|$, and the relative phase $\arg(\Gamma_{12}/M_{12})$~\cite{Buras:1984pq}.
In some extensions of the Standard Model, only the first term and, thus, $M_{12}$
is altered~\cite{Golowich:2009ii}.

The effective Lagrangian $\mathcal{L}^{\Delta C=2}$ (at energies below the $b$-quark mass) is built out of
the following operators (and their Wilson coefficients)~\cite{Hagelin:1992tc,Gabbiani:1996hi,Bagger:1997gg}:
\begin{align}
    \ensuremath{\mathcal{O}}_1 &= \bar{c}\gamma^\mu Lu\,\bar{c}\gamma_\mu Lu, &
        \tilde{\ensuremath{\mathcal{O}}}_1 &= \bar{c}\gamma^\mu Ru\,\bar{c}\gamma_\mu Ru,
    \label{eq:O1} \\
    \ensuremath{\mathcal{O}}_2 &= \bar{c} Lu\,\bar{c} Lu, &
        \tilde{\ensuremath{\mathcal{O}}}_2 &= \bar{c} Ru\,\bar{c} Ru,
    \label{eq:O2}  \\
    \ensuremath{\mathcal{O}}_3 &= \bar{c}^\alpha Lu^\beta\,\bar{c}^\beta Lu^\alpha, &
        \tilde{\ensuremath{\mathcal{O}}}_3 &= \bar{c}^\alpha Ru^\beta\,\bar{c}^\beta Ru^\alpha,
    \label{eq:O3}  \\
    \ensuremath{\mathcal{O}}_4 &= \bar{c} Lu\,\bar{c} Ru,
    \label{eq:O4}  \\
    \ensuremath{\mathcal{O}}_5 &= \bar{c}^\alpha Lu^\beta\,\bar{c}^\beta Ru^\alpha,
    \label{eq:O5}
\end{align}
where $L$ ($R$) denotes a left-(right-)handed projector on the Dirac indices, and $\alpha$ and $\beta$ are
color indices.
By parity conservation in QCD, %
$\langle D^0|\tilde{\ensuremath{\mathcal{O}}}_i|\bar{D}^0\rangle=\langle D^0|\ensuremath{\mathcal{O}}_i|\bar{D}^0\rangle$, $i=1$, 2,~3.
Thus, the five matrix elements $\langle D^0|\ensuremath{\mathcal{O}}_i|\bar{D}^0\rangle$, $i=1,\ldots5$, suffice to
describe the short-distance part of all $\Delta C=2$ processes, whether their origin is $W$-$b$ box or
something else.
In these proceedings, we report on a calculation of all five matrix elements using lattice QCD with 2+1
flavors of sea quarks.

\section{Lattice-QCD calculation}
\label{sec:qcd}

Our $D$-meson calculations have much in common with our published $B$-meson work~\cite{Bazavov:2016nty}.
We use the same ensembles (generated by the MILC collaboration) with 2+1 flavors of sea
quark~\cite{Bazavov:2009bb}.
The light quarks (valence and sea) are based on the staggered asqtad action; the heavy $c$ (or $b$) quark on
the Fermilab interpretation of the clover action.
The lattice spacings for the ensembles satisfy $a\approx0.045$~fm, $\approx0.06$~fm, $\approx0.09$~fm, and
$\approx0.12$~fm.
The sea-quark masses yield pions with
\begin{eqnarray}
    177~\text{MeV} \lesssim & M_\pi            & \lesssim 555~\text{MeV}, \label{eq:MpiG} \\
    257~\text{MeV} \lesssim & M_\pi^\text{rms} & \lesssim 670~\text{MeV}, \label{eq:Mpi-rms}
\end{eqnarray}
The ensembles contain 600--2200 gauge-field configurations, and we use 4 or 8 sources/config.

To carry out the chiral-continuum extrapolation, we take into account the subtle way in which spin emerges
for staggered fermions with staggered-Wilson four-fermion lattice operators.
The three-point correlation function, it turns out, contains contributions not only from the continuum-limit
operator of desired spin, but also some of the wrong spin~\cite{Bernard:2013dfa}.
Because only the five operators in Eqs.~(\ref{eq:O1})--(\ref{eq:O5}) can arise, we automatically have the
information needed to disentangle this effect.
We use the one-loop chiral-perturbation-theory formulas of Ref.~\cite{Bernard:2013dfa} to remove the
wrong-spin contribution in the course of our chiral-continuum fit.

The operators in Eqs.~(\ref{eq:O1})--(\ref{eq:O5}) require renormalization for any ultraviolet regulator.
We carry out the renormalization of the lattice operators corresponding to Eqs.~(\ref{eq:O1})--(\ref{eq:O5})
together with matching to \ensuremath{\overline{\text{MS}}}\ schemes in continuum~QCD.
We use a mostly nonperturbative method to handle the largest lattice-to-continuum matching
corrections~\cite{ElKhadra:2001rv,Harada:2001fi}, supplemented with a one-loop calculation of the remaining,
small renormalization parts~\cite{Evans:2009du,Bazavov:2016nty}.
We choose the renormalization scale for $D$-meson matrix elements to be $3~\text{GeV}$, while we chose $m_b$
for $B_{(s)}$ mesons.

\begin{figure}[b]
    \centering
    \vspace*{-0.125em}
    \includegraphics[width=0.8\textwidth,trim=0.75in 2.65in 0.75in 2.65in,clip]{figs/3pt}
    \vspace*{-0.125em}
    \caption{Fitting ranges for three-point correlators: triangular (green) and-or fan-shaped (magenta)
        regions for $B$ mixing (left); two-strip diagonal region (green) for $D$ mixing (right).
        Background color shows the signal-to-noise ratio from good (blue) to bad (red).}
    \label{fig:3pt}
\end{figure}
The main difference between our work on $D$ vs.\ $B_{(s)}$ mesons is the analysis of the correlation
functions.
The signal-to-noise ratio is much better for $D$-meson correlators.
For the two-point correlators, the optimal time range $t_\text{min}\lesssim t\lesssim t_\text{max}$ differs:
$t_\text{min}\approx0.7\,(0.2)$~fm, $t_\text{max}\approx3.0\,(2.4)$~fm for $D$ ($B_{(s)}$) mesons.
The difference for the three-point correlators is more striking.
We fix the four-quark operators at $t=0$ and the meson creation (annihilation) operator at time $t_x<0$
\pagebreak ($t_y>0$).
As shown in Fig.~\ref{fig:3pt}, we use a triangular and-or fan-shaped region in the $|t_x|$-$t_y$ plane for
$B_{(s)}$ mesons~\cite{Bazavov:2016nty}, while we use a long diagonal of width~2 for $D$ mixing,
$\{|t_x|=t_y\}\cup\{|t_x|=t_y+1\}$.
The long diagonal makes it easier to disentangle the lowest-lying state, if the signal persists that far.
A simultaneous fit to two- and three-point functions is used to extract the matrix elements
$\mx{\ensuremath{O}_i}\equiv\langle D^0|\ensuremath{O}_i|\bar{D}^0\rangle$.

\section{Chiral-continuum extrapolation}
\label{sec:ch-cont}

\begin{figure}[b]
    \centering
    \includegraphics[height=0.475\textheight]{figs/chipt_stability}
    \caption[fig:stab]{Stability of the chiral-continuum extrapolation for several variants of the fit
        function~$F_i$: \linebreak
        $\mx{\ensuremath{\mathcal{O}}_1}$~(left), minimized $\chi^2_\text{aug}/\text{dof}$ (right).
        Stability plots for the other $\mx{\ensuremath{\mathcal{O}}_i}$ look similar.}
    \label{fig:stab}
\end{figure}

To carry out the chiral-continuum extrapolation, we develop a fit function based on chiral perturbation
theory ($\chi$PT), Symanzik effective field theory, and heavy-quark effective theory (HQET).
It takes the form
\begin{equation}
    F_i = F_i^\text{logs} + F_i^\text{analytic} + F_i^\text{HQ disc} + F_i^{\alpha_sa^2~\text{gen}} +
        F_i^\text{renorm} + F_i^\kappa ,
    \label{eq:Fi}
\end{equation}
where $F_i^\text{logs}$ denotes the next-to-leading order description from heavy-meson rooted staggered
$\chi$PT, with nonanalytic terms including those that disentangle the wrong-spin
contributions~\cite{Bernard:2013dfa}; %
$F_i^\text{analytic}$ is a polynomial of various terms that arise in $\chi$PT at next-to-leading or higher
order;
$F_i^\text{HQ disc}$ describes heavy-quark discretization effects using HQET as a theory of cutoff
effects~\cite{Harada:2001fi}; %
$F_i^{\alpha_sa^2~\text{gen}}$ parametrizes generic cutoff effects of light quarks and gluons,
\`a~la~Symanzik; %
and $F_i^\text{renorm}$ allows the fit to be sensitive to higher orders in $\alpha_s$ for matching and
renormalization.
Finally, $F_i^\kappa$ incorporates a correction for tuning the charm-quark hopping parameter~$\kappa$,
based on extra runs at $a\approx0.12$~fm.

Both the renormalization and wrong-spin effects mix operators 1, 2, and~3 with each other, and also 4 and~5
with each other.
It is thus natural to fit the matrix elements in each sector simultaneously.
Some ingredients in $F_i^\text{logs}$ are common for all~$i$, such as masses, $f_\pi$, light-meson $\chi$PT
constants~\cite{Aubin:2004fs}, and the $D^*$-$D$-$\pi$ coupling.
We introduce these external inputs with Gaussian priors, for example $g_{D^*D\pi}=0.53\pm0.8$.
Because of these common ingredients, we choose to fit all five matrix elements simultaneously.
We form a $\chi^2$ function from $F_i-\mx{\ensuremath{\bar{O}}_i}$ and the sample covariance matrix of the
$\mx{\ensuremath{\bar{O}}_i}$, where $\ensuremath{\bar{O}}_i$ denotes the renormalized lattice operators (which differ from the
continuum $\ensuremath{\mathcal{O}}_i$ by discretization effects and higher-order matching effects).
We then augment this $\chi^2$ with Gaussian priors for the fit parameters implied in Eq.~(\ref{eq:Fi}),
choosing a central value of 0 and width of $\pm1$ in natural units for $\chi$PT and
HQET~\cite{Bazavov:2011aa} and minimize the resulting~$\chi^2_\text{aug}$.
We reconstitute the fit function at zero lattice spacing and physical quark masses to obtain our estimate of
the $\mx{\ensuremath{\mathcal{O}}_i}$ and their uncertainty.

We have 510 data points for $\mx{\ensuremath{\bar{O}}_i}$, ranging over the ensembles, valence-quark masses, and five
operators.
In our base version of $F_i$, there are 127 parameters.
To check whether the final results are robust, we repeat the procedure with several variants of~$F_i$, as
illustrated in Fig.~\ref{fig:stab}.
We express the $\chi$PT with $f_K$ instead of $f_\pi$; we choose different orders of $\alpha_s$ in
$F_i^\text{renorm}$ and even replace the mostly nonperturbative (mNPR) matching with a fully perturbative
(PT) one; we check various alternatives for the polynomial~$F_i^\text{analytic}$ (NLO, NNLO, N$^3$LO); we
check what happens when the $\chi$PT prior widths in~$F_i^\text{analytic}$ are doubled; we check
alternatives for the heavy-quark discretization errors; we substitute infinite-volume one-loop integrals for
the finite-volume sums in one-loop $\chi$PT; we omit the data from the coarsest or finest lattice spacing;
we fit each matrix element separately, thereby ignoring data constraints on wrong-spin contributions.
As one can see from Fig.~\ref{fig:stab}, the results for the $\mx{\ensuremath{\mathcal{O}}_1}$ are very stable, so we take these
variations in the fit as cross checks.
The same applies to the other $\mx{\ensuremath{\mathcal{O}}_i}$.
The largest deviations are $\sim1\sigma$ and come from fits that omit important information.
Our nearly final results for $D$ mixing are given in Table~\ref{tab:results}, together with published
results for $B_{(s)}$ mixing from Ref.~\cite{Bazavov:2016nty}.
\begin{table}[t]
    \centering
    \newcommand{\hphantom{-}}{\hphantom{-}}
    \newcommand{\hspace*{1.25em}\hphantom{-}}{\hspace*{1.25em}\hphantom{-}}
    \newcommand{\hspace*{1.25em}-}{\hspace*{1.25em}-}
    \begin{tabular}{clll}
        \hline\hline
    BBGLN~\cite{Beneke:1998sy} & \multicolumn{1}{r}{$\mx{\ensuremath{\mathcal{O}}_i}/M_D~\left(\text{GeV}^{3}\right)$} &
    \multicolumn{2}{c}{$f_{B_q}^2B_{B_q}~\left(\text{GeV}^2\right)$} \\
                    &                           & \multicolumn{1}{c}{$q=d$} & \multicolumn{1}{c}{$q=s$} \\
        \hline
    $\mathcal{O}_1$ &     $\p0.0432(29)(9)$     &       $0.0342(29)(7)$     &       $0.0498(30)(10)$    \\
    $\mathcal{O}_2$ &     $\m0.0833(38)(17)$    &       $0.0303(27)(6)$     &       $0.0449(29)(9)$     \\
    $\mathcal{O}_3$ &     $\p0.0248(16)(5)$     &       $0.0399(77)(8)$     &       $0.0571(77)(11)$    \\
    $\mathcal{O}_4$ &     $\p0.1469(69)(30)$    &       $0.0390(28)(8)$     &       $0.0534(30)(11)$    \\
    $\mathcal{O}_5$ &     $\p0.0554(38)(11)$    &       $0.0361(35)(7)$     &       $0.0493(36)(10)$    \\
    $\mu$           & \multicolumn{1}{c}{3~GeV} & \multicolumn{1}{c}{$m_b$} & \multicolumn{1}{c}{$m_b$} \\
        \hline\hline
    \end{tabular}
    \caption[tab:results]{Results for $D$ [this work] and $B$~\cite{Bazavov:2016nty} mixing in the
        renormalization scheme of Ref.~\cite{Beneke:1998sy}.}
    \label{tab:results}
\end{table}
These matrix elements (as noted above) depend on the renormalization scheme; the tabulated results are in
the \ensuremath{\overline{\text{MS}}}\ scheme with naive (fully commuting)~$\gamma^5$ and the evanescent-operator basis used by Beneke,
Buchalla, Greub, Lenz, and Nierste (BBGLN)~\cite{Beneke:1998sy}.

The MILC asqtad ensembles omit the charmed-quark sea.
As in Ref.~\cite{Bazavov:2016nty}, we assign an additional 2\% uncertainty to account for this omission.
This uncertainty is given separately, in the second set of parentheses, in Table~\ref{tab:results}.

\section{Outlook}
\label{sec:outro}

Our results agree well with and have similar uncertainty as previous lattice-QCD results from the ETM
collaboration, with 2~\cite{Carrasco:2014uya} or 2+1+1~\cite{Carrasco:2015pra} flavors in the sea.
The comparison of these results tests not only the flavor-dependence of the matrix elements but also the
sensitivity to lattice fermion formulation: ETM employs twisted-mass Wilson fermions, while we employ
staggered fermions.
All these calculations use several lattice spacings and take the continuum limit.
References~\cite{Carrasco:2014uya,Carrasco:2015pra} report the so-called ``bag factors'' often used in
phenomenology~\cite{Gabbiani:1996hi}; a detailed comparison would require choices of quark masses and decay
constants (and their uncertainties) that would obscure the error budget of one or the other set of results.
We have a set of calculations underway~\cite{Kronfeld:2015xka} to compute the $D$- and $B_{(s)}$-meson decay
constants on the same ensembles and will report the bag factors then.

Estimates of the contribution to $M_{12}$ of the second term in Eq.~(\ref{eq:mix}) range over
$(10^{-3}\text{--}10^{-2})\Gamma$ \cite{Falk:2004wg}, where $\Gamma$ is the total width of the neutral $D$
meson.
It turns out, however, that all Standard-Model phases appearing in Eq.~(\ref{eq:mix}) are small.
Thus, in a TeV-scale model that might produce a large phase in $M_{12}$, the results for the $\mx{\ensuremath{\mathcal{O}}_i}$
can be used to constrain the model's parameters.
Furthermore, until a method is developed to tame the second term in Eq.~(\ref{eq:mix}), the accuracy
achieved in this work and Refs.~\cite{Carrasco:2014uya,Carrasco:2015pra} should suffice for this purpose.

\acknowledgments

We thank our collaborators in the Fermilab Lattice and MILC Collaborations.
Computations for this work were carried out with resources provided by the USQCD Collaboration, the National
Energy Research Scientific Computing Center, and the Argonne Leadership Computing Facility, which are funded
by the Office of Science of the U.S.\ Department of Energy; and with resources provided by the National
Institute for Computational Science and the Texas Advanced Computing Center, which are funded through the
National Science Foundation's Teragrid/XSEDE Program.
This work was supported in part by the U.S.\ Department of Energy under grants
No.~DE-FG02-13ER42001 (C.C.C., A.X.K.), No.~DE-SC0015655 (A.X.K.),
No.~DE-SC0010005 (E.T.N.);
by the U.S.\ National Science Foundation under grant
PHY14-17805~(J.L.);
by the Fermilab Fellowship in Theoretical Physics (C.M.B., C.C.C.);
by the URA Visiting Scholars' program (C.M.B., C.C.C., A.X.K.);
by the MICINN (Spain) under grant FPA2010-16696 (E.G.); by the Junta de Andaluc\'ia (Spain) under Grants
No.~FQM-101 and No.~FQM-6552 (E.G.);
by the European Commission (EC) under Grant No.~PCIG10-GA-2011-303781 (E.G.);
by the German Excellence Initiative and the European Union Seventh Framework Program under grant agreement
No.~291763 as well as the European Union's Marie Curie COFUND program (A.S.K.).
Fermilab is operated by Fermi Research Alliance, LLC, under Contract No.\ DE-AC02-07CH11359 with the United
States Department of Energy.
Brookhaven National Laboratory is supported by the United States Department of Energy under contract No.\
DE-SC0012704.


\section{Introduction}
{\em Introduction.---} Diffusion phenomena are ubiquitous in nature, e.g.  in
the Brownian motion \cite{Einstein:1905}, biological membranes \cite{Siwy:2002,Arcizet:2008}, fluid systems \cite{Shraiman:2000}, engineering systems \cite{Collins:1993}, material science \cite{Fujita:2014},  and so on.
These diffusion phenomena are sources of noise and are unavoidable in experimental systems \cite{Tabard:2007,Butt:1995}.
Noise can prevent the precise manipulation of small scale systems.
However, it can also be helpful in engineering a system by using
physical phenomena such as
Logical Stochastic Resonance, where a certain amount of noise is
required for logical operation \cite{Murali:2009}.
In both cases, it is important to control noise in terms of variance of stochastic processes,
because noise is normally adjusted according to its strength.
If the variance of noise is controlled, it would be useful for various kinds of applications.

When we consider systems with noise from a theoretical viewpoint, mathematical models
are convenient, because they can numerically simulate an experiment.
A mathematical model of the Brownian motion is a typical example of a noisy system,
which can be modeled by a simple stochastic description, like the Langevin equation.
For practical purposes, it is important to consider how diffusion can be  suppressed
in the Brownian motion model.
The most commonly used approach for suppressing diffusion is to  decrease the temperature of the whole system, resulting in a suppression of thermal fluctuations,
namely the relation $D\sim~T$, where $D$ is the diffusion coefficient and $T$ is temperature.
However, the temperature is decreased for all elements in the system and is not selective for specific targets. Therefore, such a control is expensive.
We introduce an alternative method that can control the diffusion processes of one particle in the Brownian motion.

In terms of cooling nano- microscopic systems, a mathematical model of feedback cooling in electromechanical oscillator \cite{Bonaldi:2009} and a generalized model of Langevin dynamics with non-Markovian feedback \cite{Munakata:2014} have been explored. In these systems, time delay
in feedback loop is inevitable due to a time lag between the sensor and the manipulator.
Furthermore, in nonlinear systems, a time-delayed feedback control (DFC) method
has been proposed to stabilize regular
motion in fluctuated dynamics \cite{Pyragas:1992}.
In general, time delay leads to system instability. However, DFC proves that time delay can stabilize deterministic systems. On the other hand, in stochastic systems, e.g.
a random walk model with time delay, destabilization has been observed as a result of delay  \cite{Ohira:1995}.
Thus, we consider whether DFC can stabilize non-stationary stochastic systems.

In this Letter, we focus on stochastic diffusion of random walks,
corresponding to the Brownian motion in discrete time,
and apply the DFC method to control the extent of  this diffusion.
Firstly, we observe diffusion processes in a one-dimensional (1D) random walk model,
controlled by the DFC method and analyze the control in terms of stochastic delay differential equations. Secondly, we apply the method to a molecular dynamics model of two flocculated particles with thermal fluctuations
and that interact via the  Lennard-Jones (LJ) potential.
We numerically confirm that it is possible to suppress diffusion in the molecular dynamics model by the DFC method.
It is noticed that, in conventional control theory, a single path is controlled to a certain state. In contrast, the purpose of this study was to statistically control  stochastic processes.

{\em Delayed feedback control of random walk.---} In general, nonlinear functions with neutral fixed points show noise-induced diffusion phenomena.
If these functions are linearized around the neutral fixed points, the following random walk model
can be derived.
\begin{equation}
x_{n+1}=x_n+D\xi_n, \label{random},
 \end{equation}
 where $x_n$ is a state variable for discrete time $n$,  $\xi_n$ is a $m$-dimensional random variable following the normal distribution $N(0,\ 1)$ for each entry, and
 $D$ is an amplitude vector of the noise.  Without loss of generality, we consider a
 1D  case, where $m=1$.

 In order to control the random walk in $x_n$, we introduce
 a time-delayed feedback control method proposed by Pyragas \cite{Pyragas:1992} that is defined by the expression:
 \begin{eqnarray}
x_{n+1}&=&x_n+D\xi_n - K(x_n-x_{n-\tau}), \label{DFCR}
 \end{eqnarray}
 In eq. (\ref{DFCR}), the DFC is applied to the model (\ref{random}), where $K$ is
 feedback gain with a value of $0\le K<1$ and $\tau$ is the delay time.
 Any nonlinear function with noise can be linearized around its stable fixed point and
 described in the form of (\ref{DFCR}) as follows.

If we consider the system:
 \begin{equation}
 x_{n+1}=f(x_n,a),
 \end{equation}
 where the map $f$ is nonlinear and has a stable or neutrally stable fixed point. $a$ is a parameter.
We assume that the fixed point is denoted by $\tilde{x}$ and linearize
$f$ around $\tilde{x}$. This results in:
 \begin{equation}
 x_{n+1}=f(\tilde{x}) + f'(\tilde{x}) (x_n-\tilde{x}),
 \end{equation}
 which can be described as:
 \begin{equation}
 x_{n+1}=f'(\tilde{x}) x_n +(1-f'(\tilde{x}))\tilde{x},   \label{line}
 \end{equation}
 where we use the relation $\tilde{x}=f(\tilde{x})$.
 For this equation, if we denote $K=1-f'(\tilde{x})$:
 \begin{equation}
 x_{n+1}=(1-K)x_n+K \tilde{x}. \label{modify}
 \end{equation}
Regarding this linear map, we add the noise term and
 rewrite (\ref{modify}) as:
   \begin{equation}
 x_{n+1}=x_n+D \xi_n - K(x_n-\tilde{x}).  \label{no_delay}
 \end{equation}
  The dynamics of the system (\ref{no_delay}) is simple such that a trajectory
 fluctuates around the fixed point $\tilde{x}$.
It is possible to suppress diffusion of  a random walk by the form of (\ref{no_delay}), in which
 the diffusion coefficient is zero.
 On the other hand, in order to control a non-zero diffusion coefficient, we add  the control term $K(x_{n-\tau}-\tilde{x})$. Consequently, the system
 is in the form of DFC, namely eq.  (\ref{DFCR}).

 It should be noted that this discussion is relevant for systems with only stable fixed points and thus, with no diffusion observed.
 By the procedure above, such stable systems are converted to systems showing diffusion because of the
 neutral fixed points.

Figure \ref{sup_RW} shows the time series of the  (\ref{DFCR}) system dynamics.
It is observed that a longer $\tau$ suppresses the diffusion of a random walk.
This result is counterintuitive, since long time delay should  make a system
unstable. However, in the system (\ref{DFCR}), the opposite is observed.
Thus, we quantify the diffusion decay  with respect to $\tau$.
In Figure \ref{DC}, the diffusion coefficients ${\cal D}$, defined by $\lim_{T \to T_{\infty}} \langle (X_T - X_0)^2 \rangle/T$, is shown as a function of $\tau$. The result is averaged over $1000$ realizations.
We takes $T_{\infty}=100000$. It is evident that ${\cal D}$ decays by $\tau^{-2}$. The decay order is illustrated in the following.

\begin{figure}[htbp]
\begin{center}
\rotatebox{-90}{\includegraphics[width=5.5cm]{fig1.eps}}
\end{center}
\caption{Suppression of diffusion in random walks of the system (\ref{DFCR}). The red solid and green dashed lines represent when $\tau=2$ and
$\tau=20$, respectively. $D=0.001$.}\label{sup_RW}
\end{figure}

\begin{figure}[htbp]
\begin{center}
\rotatebox{-90}{\includegraphics[width=5.5cm]{fig2.eps}}
\end{center}
\caption{${\cal D}$ with respect to $\tau$, averaged over $1000$ realizations. The dashed line is proportional to $1/(1+K\tau)^2$. $K=0.5$. }\label{DC}
\end{figure}

We consider a continuous time (\ref{DFCR}) system:
\begin{equation}
\dot{y}_t=D\xi_t - K(y_t - y_{t-l}), \label{conti}
\end{equation}
where $y$ is a real variable and $t$ and $t-l$ are continuous time points.
For this system, let $\Delta=l/\tau$, $\tau \in \{2,3,\ldots\}$. Then, we define $t_n$ as $n\Delta$,  that is $n\frac{l}{\tau}$, where
$n\in \{-\tau,-\tau+1,\ldots,0,1,\ldots \}$.   If we define $Y_n=y_{t_n}$ and $\xi_n=\xi_{t_n}$,  eq. (\ref{conti}) can
be rewritten as $Y_{n+1}=Y_n+D\Delta \xi_n - \Delta K (Y_n - Y_{n-\tau})$.  This form is
equivalent to eq. (\ref{DFCR}).

We consider the continuous time system (\ref{conti}) using stochastic delay differential equations.
In the standard setting of stochastic calculus, eq. (\ref{conti}) is now redefined as
\begin{equation}
dy_t = K(y_{t-\tau}-y_t) dt + D dw_t;\ t \geq 0 \label{SDDE0_main}
\end{equation}
where $\{w_t\}_{0\leq t \leq T}$ is a standard Wiener process on a filtered probability space $(\Omega, {\cal F}, \{{\cal F}_t\}, P)$,
and $y_s$ for $s\in [\tau,0]$ is assumed to be
a given ${\cal F}_0$-measurable random variable $Z_s$. It should be noted that
${\cal F}_t$ contains all information for the realization of $w_s$ and $y_s$ for $s\leq t$.
Thus, given ${\cal F}_t$, all of these variables can be all treated as constants.
In the following, the method for determining the diffusion coefficient ${\cal D}$ is briefly described.
For details, see Supplemental Materials.

Equation (\ref{SDDE0_main}) has the unique solution of
\begin{eqnarray}
y_t &=& y^Z_t + D \int_0^t y^0_{t-s }dw_s; 0 \leq t, \label{Y_main} \\
y_t &=& Z_t; -\tau \leq t \leq 0 \nonumber
\end{eqnarray}
where
\begin{eqnarray}
y^Z_t &:=& y^0_t Z_0 + K \int_{-\tau}^{0} y^0_{t-s-\tau} Z_s ds \nonumber\\
y^0_t &:=& \sum_{n=0}^{\lceil t/\tau \rceil} \left( K^n \over n!\right)
(t- n\tau)^n \exp(-K(t-n\tau)). \label{y0_main}
\end{eqnarray}
(See \cite{Kuchler:1992}).
It should be noted that  $y^0_t$ is a deterministic, rather than a  stochastic process.

Then, we can express ${\cal D}$ in terms of $y^0_t$ as:
\begin{equation}
{\cal D} =
\lim_{T\rightarrow \infty} {D^2 \over T}  \int_0^T (y^0_{s})^2 ds. \label{D_main}
\end{equation}
Moreover, we examined how the function $y^0_t$ behaves in our numerical computations of the equation (\ref{y0_main}),
and found that it seemed converging to a constant of $\bar{y}\propto {1\over 1+ K\tau}$,
as described in the Supplemental Materials and in \cite{Kobayashi:prep}.
As a consequence, ${\cal D}$ is given by:
\begin{eqnarray*}
{\cal D} = (D\bar{y})^2.\nonumber
\end{eqnarray*}
In addition, it should be noted that this $\bar{y}$ and ${\cal D}$ decreases with increasing $\tau$,
which is consistent with Figure \ref{DC}.

Here, the characteristics of our model should be stressed.
In \cite{Kuchler:1992},  the equivalent conditions to stationary solutions of $y=\{y_t\}$ are given,
which are not met in our system (\ref{Y_main}).
In the case of the stationary $y$, its diffusion coefficient ${\cal D}$ goes to zero as $T \rightarrow \infty$.
Conversely, ${\cal D}$ seems to explode when the solution $y$ is non-stationary.
In this sense, our model with the non-zero ${\cal D}$ lies at the boundary between a stationary and non-stationary state.

This can be seen in another intuitive way.
If a limiting process of our model (\ref{SDDE0_main}) is considered with $\tau \rightarrow 0$:
\begin{eqnarray*}
dy^{\downarrow}_t = D dw_t.\nonumber
\end{eqnarray*}
It should be noted that this is just a scaled random walk, which still includes  ${\cal D}=D$.
On the other hand,
${\cal D} = 0$,
under another limiting process with $\tau \rightarrow \infty$:
\begin{eqnarray*}
dy^{\uparrow}_t = K(y_{-\infty} - y^{\uparrow}_t )dt + D dw_t\nonumber
\end{eqnarray*}
where $y_{-\infty}:=\lim_{t\rightarrow -\infty} y_t$ is given as an initial condition
and as a constant target control level.

As derived before , in our model,  the diffusion coefficient is given by ${\cal D}=D^2\times \left({1\over 1+K\tau} \right)^2$.
Thus, any level of ${\cal D}$ between $(0,D]$ can be achieved with an appropriate choice of $K$ and $\tau$.

{\em Application to Molecular Dynamics. ---} In the following, we apply the control scheme which is shown in the above to a molecular dynamics model.

We consider the dynamics of two particles with noise, and attractive and dissipative forces in a 1D space:
\begin{eqnarray}
  {dx_i \over dt}&=&v_i, \\
  {dv_i \over dt}&=&P_i(\{ x_j\}^2_{j=1})+D_i(\{ v_j\}^2_{j=1})+ R \xi, \label{second}
\end{eqnarray}
where $x_i$ and $v_i$ are the position and velocity of the $i$-th particle ($i=1, 2$),
and $\xi$ is noise from the Gaussian distribution with intensity, $R$.
For simplicity, we assume that the mass of the $i$-th particles, $m_i$=1.
The first term in eq. (\ref{second}) is a conservative force derived as a derivative of the potential $U$ written as:
\begin{eqnarray}
  P_i &=&-{d U(\{ x_j\}^2_{j=1}) \over dx_i}, \\
  U(\{ x_j\}^2_{j=1}) &=& \sum_{i=1}^{2} \sum_{j=i+1}^2 \tilde{V}(x_{ij}).
\end{eqnarray}
where $\tilde{V}$ is a central force potential, and $x_{ij}=x_i-x_j$.
In this paper, we use the Lennard-Jones potential as the central potential as given by:
\begin{eqnarray}
  \tilde{V}(x)=\epsilon(({\sigma \over x})^{12} - ({\sigma \over x})^6).
\end{eqnarray}
The second term in eq. (\ref{second}) represents a dissipative force:
\begin{eqnarray}
  D_i=-\gamma v_i,
\end{eqnarray}
where $\gamma$ is a coefficient of drag.
We set $\sigma=2^{-1/6}$ so as that the LJ potential $\tilde{V}(x)$ has a minimum at $x=1$ and $\epsilon=1$.

\begin{figure}
\begin{center}
\includegraphics[width=8.4cm]{fig3.eps}
\caption{Time evolutions of the positions of two particles without control inputs.
Particles flocculate because of the interacting force, and move together diffusively.
\label{fig:diffusion}}
\end{center}
\end{figure}

We calculate dynamics of two particles with initial conditions $x_1=0.8, v_1=0.0, x_2=5.8$, and $v_2=0.0$.
In this case, the particles flocculate becase of the attracting force from the LJ potential, and move together diffusively (see Fig. \ref{fig:diffusion}).

In this application, we control fluctuations in the system by using the time-delayed feedback control method.
Velocities fluctuate randomly around $0$ because of the Gaussian noise and are stable at 0 without the noise.
Thus, the dynamics of the velocities is the same as the Gaussian noise.
This Gaussian noise-like behavior  makes the dynamics of the positions analogous to a simple random walk, $ {dx_i \over dt} \sim \xi$.
Thus, control inputs are given only according to the dynamics of the positions, which are driven by noise originally from the velocities.
The control scheme for the molecular dynamics model is written as:
\begin{eqnarray}
  {dx_i \over dt}&=&v_i-K(x_i-x_i(t-\tau)), \\
  {dv_i \over dt}&=&P_i(\{ x_j\}^2_{j=1})+D_i(\{ v_j\}^2_{j=1})+ R \xi.
 \end{eqnarray}
Control signals are injected  when the  particles are flocculated. Figure \ref{fig:particle-control} shows the dynamics of the center positions of two particles without control (green dotted line) and with control with $\tau=10$ (red solid line).
The movement of the flocculated particles is diffusive because of noise.
The figure shows that the diffusion of the flocculated particles is suppressed by our control method.

Figure \ref{fig:tau-vs-diff} shows the dependence of the ${\cal D}$ with  $\tau$.
As is the case with the simple random walk,
the ${\cal D}$ decreases with a change in $\tau$.

\begin{figure}
\vspace{-15mm}
\begin{center}
\includegraphics[width=8.4cm]{fig4.eps}
\caption{Time evolutions of the center position of two particles. The green and the red lines are the time evolutions without control inputs and with control input ($\tau=50$) , respectively.}
\label{fig:particle-control}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\rotatebox{-90}{\includegraphics[width=6cm]{fig5.eps}}
\caption{The diffusion coefficient ${\cal D}$ of the molecular dynamics model as a function of  $\tau$. }
\label{fig:tau-vs-diff}
\end{center}
\end{figure}

{\em Summary.---} In this work, we have applied  time-delayed feedback control
to a 1D random walk model.
The interaction between noise and time delay is non-trivial.
We have counterintuitively observed that the diffusion coefficient decreases with increasing delay time. Thus, we have analytically explained the decay of  the diffusion coefficient by solving stochastic delay differential equations derived from the controlled system. In principle, it is possible to control the amount of diffusion from zero to a specific value.

We have applied the proposed method to a molecular dynamics model described by
continuous time differential equations. We have numerically demonstrated the control of
two flocculated particles with the LJ potential.
Consequently, diffusion generated by noise was successfully controlled and
decay of the diffusion coefficient was observed in this system.

The proposed method does not change physical properties of the system in comparison with the conventional technique of reducing diffusion, i.e.  decreasing
temperature, which can change the thermophysical properties.
Moreover, the previous analysis of delayed random walks is based on the assumption of a stationarity  \cite{Ohira:1995,Ohira:2000}. In contrast,
our method is non-stationary and  a conventional analysis by stochastic delay differential equations in \cite{Kuchler:1992} cannot be applied.
In addition to such theoretical findings, we remark that
the method can be applied to
nano- microscopic resonators for reducing the thermal noise.
In those experimental systems, time delay cannot be avoided or shortened less than some length \cite{Bonaldi:2009,Munakata:2014}. However, making longer time delay can be easy so that our method is effective in the systems.

Other types of time-delayed feedback control are also feasible. For example, we have proposed the adaptive time delay system as follows \cite{Ando:2015}:
$x_{n+1}=x_n+D\xi_n-1/x_{n-\tau}(x_n-x_{n-\tau})$. This form of adaptive time-delayed feedback control shows almost similar performance to the original system (\ref{DFCR})
in terms of control of the diffusion coefficient. The difference between the original and
the adaptive control is in determining $K$ or determining the initial values for
the system. We can apply the adaptive control scheme to real world systems depending on
the availability of the control parameters. If we cannot access $K$, adaptive feedback gain
$1/x_{n-\tau}$ can be alternatively used.