DRL-based Distributed Resource Allocation for Edge Computing in Cell-Free Massive MIMO Network

In this paper, with the aim of addressing the stringent computing and quality-of-service (QoS) requirements of recently introduced advanced multimedia services, we consider a cell-free massive MIMO-enabled mobile edge network. In particular, benefited from the reliable cell-free links to offload intensive computation to the edge server, resource-constrained end-users can augment on-board (local) processing with edge computing. To this end, we formulate a joint communication and computing resource allocation (JCCRA) problem to minimize the total energy consumption of the users, while meeting the respective user-specific deadlines. To tackle the problem, we propose a fully distributed solution approach based on cooperative multi-agent reinforcement learning framework, wherein each user is implemented as a learning agent to make joint resource allocation relying on local information only. The simulation results demonstrate that the performance of the proposed distributed approach outperforms the heuristic baselines, converging to a centralized target benchmark, without resorting to large overhead. Moreover, we showed that the proposed algorithm has performed significantly better in cell-free system as compared with the cellular MEC systems, e.g., a small cell-based MEC system.

I. INTRODUCTION
Over the past few years, there has been a surge in computation hungry, and yet highly delay-intolerant multimedia applications. To support these applications in resource-limited end-user devices, it is essential to augment on-board (local) processing by offloading part of the intensive computation to more powerful computing platforms at the network edge. To this end, the ubiquitous radio and computing resources at the mobile edge network should be jointly optimized in a way to improve certain system objectives such as execution delay and energy consumption. In this regard, a plethora of joint communication and computing resource allocation (JCCRA) schemes based on conventional optimization methods have been investigated. However, most of these approaches either assume accurate knowledge of network-wide information, which is difficult to obtain in practice, or are limited to quasi-static systems, such as deterministic task size and timeinvariant channel conditions. Thus, their application to support time-sensitive applications in dynamic mobile edge network might be critically limited.
Driven by the recent advances in machine learning and deep reinforcement learning (DRL), several learning-based schemes have been proposed to provide flexible joint resource allocation. For instance, deep Q-network (DQN) is employed in [1,2] for designing joint offloading and resource allocation decisions, while [3] presents a deep deterministic policy gradient (DDPG)-based scheme in heterogeneous networks with multiple users. However, these schemes rely on central entities to make joint resource allocation decision, suffering from considerable signaling overhead and scalability issues. [4] investigates distributed JCCRA problem in a single-cell mobile edge computing (MEC) system, however, the adopted computation model is not applicable for tasks with extreme reliability and ultra-low hard deadline requirements. We note that the large amount of effort on optimizing joint resource allocation in MEC networks focuses on cellular system. However, as we move towards the next decade, the current cellular MEC systems may not keep up with the diverse and more stringent requirements of the envisaged advanced applications, such as holographic displays, multi-sensory extended reality, and others.
In light of the recently proposed cell-free massive MIMO architecture, envisioned for beyond-5G and 6G networks, however, there are only limited works. Notably, cell-free architecture can provide sufficiently fast and reliable access links without cell-edge by serving a relatively small number of users from several geographically distributed access points (APs) [5]. It, therefore, opens up a new opportunity to support energy-efficient and consistently-low latency computational task offloading, in contrast to cellular MEC systems. [6] presents several performance analyses in a cell-free MEC system considering coverage radius of the APs, while the issues of active user detection and channel estimations are discussed in [7]. In this paper, we consider a JCCRA problem that minimizes total energy consumption of the users while meeting the user-specific application deadlines by jointly optimizing the allocation of local processor clock speed and uplink transmission power for each user. We then present a fully distributed solution approach based on cooperative multi-agent reinforcement learning, alleviating the signaling and communication overheads associated with centralized implementations. In sharp contrast to our work, both [6] and [7] Fig. 1: User-centric cell-free massive MIMO-enabled mobile edge network: Illustrative system model with the JCCRA problem. To the best of our knowledge, this is the very first attempt to solve JCCRA problem in a fully distributed fashion for cell-free network-enabled mobile edge network. The fully distributed and intelligent JCCRA in our framework combined with the reliable access links of the cellfree massive MIMO architecture can be a promising means to handle the stringent requirements of the newly introduced multimedia services. The rest of the paper is organized as follows: In Section II, we discuss the system model and problem formulation for the JCCRA framework. Section III presents the proposed distributed JCCRA scheme based on cooperative multi-agent reinforcement learning (MARL) framework. Simulation results are presented in Section IV, followed by some concluding remarks in Section V.

II. SYSTEM MODEL AND PROBLEM FORMULATION
A. User-centric Cell-free Massive MIMO: Overview We consider a cell-free massive MIMO-enabled mobile edge network, depicted in Fig. 1, where a central processing unit (CPU) with an edge server of finite capacity f CP U (in cycles per second) is connected to M single-antenna access points (APs) via ideal fronthaul links. Furthermore, M APs serve K single-antenna users, such that M K, entailing to the widely known system model for cell-free massive MIMO [5]. Let M = {1, 2, ..., M } and K = {1, 2, ..., K} denote sets of AP and user indices, respectively. Besides, each user k ∈ K is also equipped with limited local processor with maximum capability f max k (in cycle per second). Let the channel between m-th AP and k-th user is given as g mk = β 1 /2 mk h mk where β mk is a large-scale channel gain coefficient that models the effects of path loss and shadowing while h mk ∼ CN (0, 1) represents small-scale channel fading. Let τ c denote a channel coherence period in which h mk remains the same. The coherence period τ c is divided into pilot transmission interval of τ p samples and uplink data transmission interval of (τ c − τ p ) samples for offloading computation to the edge server in the CPU.
At the beginning of each coherence time, each user k ∈ K simultaneously transmit a pilot sequence ψ k ∈ C τp× 1 , such that ψ k 2 = 1, which is used to estimate the respective channelsĝ mk at each AP m ∈ M. Assuming pairwise orthogonal sequences {ψ 1 , ψ 2 , ....., ψ K }, i.e., τ p = K, we ignore pilot contamination problem in the current model. Then, the received pilot vector at the m-th AP, y p m ∈ C τp× 1 , can be represented as where q p k is the pilot transmit power, and ω p m is a τ pdimensional additive noise vector with independent and identically distributed entries of ω p m ∼ CN 0, σ 2 m . Based on the received vector y p m , the least-square (LS) channel estimatê g mk can be expressed as: The channel estimates are used to decode the uplink transmitted data of the users. After pilot transmission, the users transmit offloaded data to the APs. Let x k denote the uplink data of user k, with transmit power of p k , which is determined as p k = η k p max k , where η k , and p max k represent power control coefficient and maximum uplink power of the k-th user, respectively. Then, the received signal y u m at the m-th AP is represented as It is important to ensure the scalability of cell-free massive MIMO system, in terms of complexity for pilot detection and data processing. To this end, we limit the number of APs that serve the k-th user to N k < M , and form a user-centric cluster of APs C k ⊂ {AP 1 , AP 2 , · · · , AP N k }, by including all APs with the largest β mk to C k until the maximum allowable cluster size is reached, i.e., N k = C max . Thus, the transmitted data by user k is only decoded by the APs in C k . Then, after performing maximum ratio combining (MRC) at each AP m ∈ C k the quantityĝ * mk y u m is transmitted to the CPU via a fronthaul link. The received soft estimates are combined at the CPU to decode the data transmitted by user k as follows: Then, the uplink signal-to interference and noise ratio (SINR) γ k for user k can be expressed as: Given the system bandwidth W , the uplink rate of user k is given as R k = W log 2 (1 + γ k ) by Shannon's theory.

B. Parallel Computation Model
We assume each user k has a computationally intensive task of T k (t) bits at the beginning of every discrete time step t, whose duration ∆t = τ c . The task sizes in different time steps are independent and uniformly distributed over [T min, T max, ], for every user k ∈ K. Let t d k denote the maximum tolerable delay of user k to complete execution of the given task. Considering the delay constraint and energy consumption, the k-th user processes T local k (t) bits locally while offloading the remaining T of f load k (t) bits to the edge server at the CPU. 1], denote the k-th user processor speed (in cycle per second) allocated for local task execution. The entire task is offloaded to the edge server if α k = 0, while α k = 1 implies that the whole local processing capability is fully utilized and the remaining task bits are offloaded to the edge server. Accordingly, the proportion of locally computed task size at time step t is expressed as T local where N cpb corresponds to the number of cycles required to process one-bit task. Then, the time taken for local execution at time step t is expressed as The corresponding energy consumption is given by where ς corresponds to the effective switched capacitance.
2) Computation Offloading: The remaining T of f load k (t) = max 0, T k (t) − T local k (t) bits are offloaded to the edge server at the CPU for parallel computation. While offloading computation to the edge server, the experienced latency can be broken down into transmission delay, computing delay at the edge server, and transmission delay for retrieving the processed result. Assuming the retrieved data size after computation is much smaller compared to the offloaded data size, we ignore the retrieving delay. Then, the transmission delay to offload T of f load , where R k (t) is the uplink rate of user k at time step t. The computing time at the edge server is given by t comp is computation resource at the CPU allocated for user k, which is in proportion to T of f load k (t). The total offloading delay for edge computation is then given as The corresponding energy consumption for offloading computation is then expressed as where p k (t) = η k (t) p max k , ∀η k ∈ [0, 1]. Then, the overall experienced latency t k (t) by user k to execute T k bits locally and at the edge is given as: Similarly, the energy consumption E k (t) of user k at time step t can be expressed as C. JCCRA Problem Formulation According to the aforementioned communication and parallel computation models, we intend to jointly optimize the local processor speed f local k (t) and uplink transmission power p k (t) for every user k ∈ K at each time step t, in order to minimize the total energy consumption of all users subject to the respective user-specific delay requirements. The JCCRA problem can mathematically be formulated to jointly determine (α k (t) , η k (t)) , ∀k ∈ K as follows: (12) is a stochastic optimization problem in which the objective function and the first condition involve constantly changing random variables, i.e., task size and channel conditions. Consequently, the problem must be solved frequently, i.e., at every time step t, and the solution must converge rapidly within the ultra-low delay constraint. Therefore, it is challenging to solve the problem using iterative optimization algorithms with reasonable complexity. Furthermore, relying on a one-shot optimization, it is difficult to guarantee a steady long-term performance with conventional algorithms since they cannot adapt to the changes in the mobile edge network.

III. THE PROPOSED MULTI-AGENT REINFORCEMENT
LEARNING-BASED DISTRIBUTED JCCRA To cope up with the dynamics in the MEC while supporting efficient and adaptive joint resource allocation for every user, we propose a fully distributed JCCRA based on cooperative multi-agent reinforcement learning (MARL). We note that solving the JCCRA problem centrally at the CPU might allow for efficient management of the available radio and computing resources due to globally processing of network-wide information. However, the associated signaling and communication overheads are forbiddingly significant, especially given the ultra-low delay constraints. Accordingly, we transform the JCCRA problem into a partially observable cooperative Markov game that consists of K learning agents. One of the main challenges in multi-agent domain is the nonstationarity of the environment since the agents update their policies concurrently, which might lead to learning instability. To deal with this challenge, the agents are trained with the state-of-the-art multi-agent deep deterministic policy gradient (MADDPG) algorithm under the framework of centralized training and decentralized execution [8].

A. DRL Formulation for JCCRA
Each user is implemented as a DRL agent that learn to map local observation of the environment to efficient JCCRA actions by sequentially interacting with the mobile edge network in discrete time steps. Let O k , and S denote the local observation space of agent k ∈ K, and the complete environment state space, respectively. At the beginning of each step t, agent k perceives the local observation of the environment state o k (t) : S → O k , and takes an action a k (t) ∈ A k from its action space A k according to the current JCCRA policy µ k . The interaction with the environment produces the next observation of the agent o k (t + 1) ∈ O k and a realvalued scalar reward r k (t) : S × A k → R that guides the agent to constantly improve the policy until it converges. The goal of the agent is, therefore, to derive an optimal JCCRA policy µ * k that maximizes the expected long-term discounted cumulative reward, defined as where ε ∈ [0, 1] is a discount factor and T is the total number of time steps (horizon). In the sequel, we define the local observation, action, and reward of agent k for the JCRRA at a given time t.
1) Local observation: At the beginning of each time step t, the local observation of agent k is defined as o , and R k (t−1) are the incoming task size, user-specific application deadline, and the uplink rate from the previous time step, respectively. Note that the local observation is only a partial view of the environment, i.e., it does not include information about the other agents.
2) Action: Based on the local observation of the environment, the agent takes action a k (t) , and uplink transmit power p k (t) = η k (t) p max k .

3) Reward:
To reflect the design objective of the JCCRA problem and enforce a cooperative behavior among the agents, we define joint reward as where ξ k = 1 if t k (t) ≤ t d k , otherwise ξ k = 10 to punish potentially selfish behavior of the agents that leads to failure in meeting the delay constraint.

B. MADDPG Algorithm for Distributed JCCRA
During the offline training phase, the agents can be trained in a central unit, which presents an opportunity to share extra information among the agents for easing and accelerating the training process. Specifically, agent k has now access for the joint action a (t) = (a 1 (t) , ..., a K (t)), and full observation of the environment state s is the local observation of other agents at time step t. The extra information endowed to the agents, however, is discarded during the execution phase, meaning the agents rely on their local observation entirely to make JCCRA decisions in a fully distributed manner.
As an extension of the single-agent deep deterministic policy gradient (DDPG) algorithm [9], the MADDPG is an actorcritic policy gradient algorithm. Specifically, the MADDPG agent k employs two main neural networks: actor network with parameters θ µ k to approximate a JCCRA policy µ k (o k |θ µ k ) and a critic network, parametrized by θ Q k , to approximate a state-value function Q k s k , a|θ Q k , along with their respective time-delayed copies, θ µ k and θ µ k , which serve as targets. At time step t, the actor deterministically maps local observation o k (t) to a specific continuous action µ k (o k (t) |θ µ k ) and then, a random noise process N k is added to generate an exploratory policy such that a k (t) = µ k (o k (t) |θ µ k ) + N k (t). The environment, which is shared among the agents, collects the joint action a (t) = {a k (t) , ∀k ∈ K} and returns the immediate reward r k (t) and the next observation o k (t + 1) to the respective agents. To make use of the experience in later decision-making steps efficiently and improve the stability of the training, the agent's transition along with the extra information e k (t) = (s k (t) , a (t) , r k (t) , s k (t + 1)) is saved in the replay buffer D k of the agent.
To train the main networks, a mini batch of B samples , is randomly drawn from the replay buffer D k , where i denotes a sample index. The critic network is updated to minimize the following loss function: where y i k is the target value expressed as y i On the other hand, the parameters of the actor network are updated along the deterministic policy gradient given as The target parameters in both actor and critic networks are updated as θ µ k ← τ θ µ k + (1 − τ ) θ µ k , and θ Q k ← τ θ Q k + (1 − τ ) θ Q k , respectively, where τ is a constant close to zero. It is worth mentioning that the proposed JCCRA scheme is trained offline. Once it converges, real-time decision making can be conducted in a fully distributed fashion by simply performing inference with the trained actor network of each agent, relying on the respective local observation only.

IV. PERFORMANCE ANALYSIS
We consider M = 100 APs and K = 10 users uniformly distributed at random over an area of 1 km 2 . The APs are connected to the CPU via ideal fronthaul links. We assume only 30% of the entire APs are clustered to serve each user, i.e. N k = 0.3M, ∀k ∈ K. The system bandwidth is set to 5MHz and all active users share this bandwidth without channelization. Denoting the standard deviation of the shadow fading by σ sh , the large-scale channel gain is given as β mk = 10 P L mk 10 10 σ sh z mk 10 , where z mk ∼ N (0, 1). The path loss P L mk is given according to the three-slope model [10] as follows: Here, f is the carrier frequency (in MHz), h AP is the AP antenna height (in m), and h u denotes the user antenna height (in m). Further, no shadowing is considered unless d mk > d 1 .
The edge server at the CPU has computation capacity f CP U = 100 GHz, while the users are equipped with local processor with f max k = 1 GHz. The length of a time step ∆t is set to 1ms, and each user generates a task of a random size that is uniformly distributed in the range of 2.5 to 7.5 kbits at the beginning of every time step t. Furthermore, Furthermore, additional simulation parameters are set to N cpb = 500, ς = 10 −27 , p max k = 0.1W, σ sh = 10dB, and t d k = ∆t = τ c = 1ms. For the proposed MADDPG-based JCCRA implementation, the actor and critic networks of each agent are implemented with fully connected neural networks, which consist of 128, 64, and 64 nodes. All the hidden layers are activated by ReLu function, while the outputs of the actor are activated by sigmoid. The parameters of critic and actor networks are updated with adaptive moment (Adam) optimizer with the learning rates of 0.001 and 0.0001. Further, we set the discount factor γ k = 0.99, target update parameter τ = 0.005, and mini batch size B = 128. In our simulations, the agents are trained for several episodes, each consisting of 100 steps.
To evaluate the performance of the proposed MADDPGbased JCCRA, we compare it with three benchmarks: 1) Centralized DDPG-based JCCRA scheme: This approach refers to a DDPG-based centralized resource allocation scheme implemented at the CPU. We adopt the same neural network structure and other hyperparameters as the MADDPG-based scheme to train the actor and critic in the single DDPG agent. Since global information of the entire network is processed centrally at the CPU to make JCCRA decisions for all users, we can potentially obtain the most efficient resource allocation. Consequently, this baseline serves as a target performance benchmark. However, as it involves significant signaling and communication overheads, this scheme might be infeasible to support time-sensitive applications.
2) Offloading-first with an uplink fractional power control (FPC) scheme: This approach preferably offloads the computation to the edge server with the aim of aggressively exploiting the reliable access link provided by cell-free massive MIMO while saving the local processing energy consumption. The uplink transmit power for the k-th user is determined by the standard fractional power control (FPC) [11] as follows: where p 0 = −35 dBm, λ k = N k m=1 β mk , and ν = 0.5.
3) Local execution-first with an uplink fractional power control (FPC) scheme: The entire local processing capability is fully utilized, i.e., f local k = f max k , and the remaining task bits are offloaded to the edge server with uplink transmit power given according to (18).
In Fig. 2, we investigate the convergence and performance of the proposed MADDPG-based JCCRA scheme by evaluating the total reward accumulated during the training process. It can be observed that the performance of the proposed distributed scheme has converged to the target performance benchmark, i.e. the centralized DDPG-based JCCRA, while relieving the large overhead of the latter. This implies that by relying on local observation and additional information provided during the offline centralized training phase, the agents manage to learn efficient joint resource allocation. Moreover, the proposed scheme has significantly outperformed the heuristic baselines in terms of total reward accumulated.
In Fig. 3, we compare the average success rate of the users in attaining the delay constraint of the time-sensitive applications. One can observe that the proposed MADDPGbased JCCRA has achieved more than 99% success rate in accomplishing the tasks within the deadline as the training episodes progress. Even though local execution-first policy achieves 100% success rate, the performance comes at the cost of high energy consumption, as reflected in Fig. 2. Moreover, the fact that the policy of the learning agents outperforms the offloading-first baseline implies even if cell-free massive MIMO can provide reliable access links to the users, aggressive computation offloading can result in performance degradation due to ineffective use of the limited resources. We then implemented the proposed MADDPG-based JC-CRA scheme in two different cellular MEC architectures, including a small-cell system and a single-cell system with colocated massive MIMO, to further investigate the gain obtained from cell-free architecture. For the small-cell MEC system, each user is served by a single AP with the largest largescale fading coefficient according to the methodology in [5], whereas in co-located massive MIMO case, each user connects to a base station with total number of N k antennas. From Fig.'s 4(a) and 4(b), we can observe that the proposed algorithm over a cell-free system has significantly outperformed the cellular MEC systems, in terms of average success rate and total reward accumulated. This is attributed to the channel quality degradation due to spatial co-channel interference from other near-by cells in the small-cell system case, while in the singlecell MEC system, users are subjected to large interference from all co-located antennas under MRC, hence limiting the uplink rate for computation offloading.

V. CONCLUSION AND FUTURE WORKS
In this paper, we presented a novel distributed joint communication and computing resource allocation (JCCRA) scheme based on cooperative multi-agent reinforcement learning framework for the cell-free massive MIMO-enabled mobile edge network. More specifically, each trained agent can make intelligent and adaptive joint resource allocation in real-time, based on local observation only, in order to minimize the total energy consumption while meeting the respective delay constraints of the dynamic MEC system. The performance of the proposed distributed JCCRA scheme is shown to converge to the centralized target baseline, without resorting to large overhead. Finally, we have shown that the proposed algorithm in cell-free system has outperformed two cellular MEC systems by wide margins. In the future, it would be interesting to extend the current formulation to determine a set of APs that are adaptive to dynamic user mobility.


CoLI-Machine Learning Approaches for Code-mixed Language Identification at the Word Level in Kannada-English Texts

The task of automatically identifying a language used in a given text is called Language Identification (LI). India is a multilingual country and many Indians especially youths are comfortable with Hindi and English, in addition to their local languages. Hence, they often use more than one language to post their comments on social media. Texts containing more than one language are called"code-mixed texts"and are a good source of input for LI. Languages in these texts may be mixed at sentence level, word level or even at sub-word level. LI at word level is a sequence labeling problem where each and every word in a sentence is tagged with one of the languages in the predefined set of languages. In order to address word level LI in code-mixed Kannada-English (Kn-En) texts, this work presents i) the construction of code-mixed Kn-En dataset called CoLI-Kenglish dataset, ii) code-mixed Kn-En embedding and iii) learning models using Machine Learning (ML), Deep Learning (DL) and Transfer Learning (TL) approaches. Code-mixed Kn-En texts are extracted from Kannada YouTube video comments to construct CoLI-Kenglish dataset and code-mixed Kn-En embedding. The words in CoLI-Kenglish dataset are grouped into six major categories, namely,"Kannada","English","Mixed-language","Name","Location"and"Other". The learning models, namely, CoLI-vectors and CoLI-ngrams based on ML, CoLI-BiLSTM based on DL and CoLI-ULMFiT based on TL approaches are built and evaluated using CoLI-Kenglish dataset. The performances of the learning models illustrated, the superiority of CoLI-ngrams model, compared to other models with a macro average F1-score of 0.64. However, the results of all the learning models were quite competitive with each other.

INTRODUCTION
The measure of mineable information is increasing quickly with the rapid growth of social media. In a country like India where multilingualism is popular, people are comfortable in using more than one language and hence usually use a combination of two or more languages to post their comments or messages on social media. However, these comments may be using single script or multiple scripts. The combination of two or more languages in any text is called code-mixing and is gaining popularity among younger generations mainly to use on social media. English is considered as one of the languages for communication in many countries and the keyboard layout of computers and smartphones by default is of Roman script. Even though there are many apps which can be used to write the text in local languages, however, due to technological glitches most of the users prefer Roman script to write the comments in local or code-mixing languages. Analysis of code-mixed text defines a new research trend due to many challenges.
As social media content is not governed by the syntax of any of the languages, short sentences are quite common in addition to incomplete sentences and even words. Words may have a high level of typographical errors intentionally holding creative spellings (gr8 for 'great'), phonetic typescript, word play (goooood for 'good'), and abbreviations (OMG for 'Oh my God!'). Generally, the non-English speakers use English words/sentences (through code-mixing and Anglicism) instead of composing online media text using unicode in their languages. They frequently mix multiple languages in comments/messages to express their thoughts on social media making the analysis of code-mixed text an extremely challenging task. The preliminary step in analyzing code-mixed texts for various applications is identifying the languages used in these texts efficiently as accuracy of the applications depend on the proper identification of languages. Languages may be mixed at paragraph level, sentence level, word level, or even within a word. Despite a lot of work being done in LI, the problem of LI in code-mixed scenario is still a long way from being illuminated [1]. A code-mixed scenario where words of one language are transcribed with words of other languages as prefix or suffix has lot more troubles, particularly due to conflicting phonetics. In such case, proper context can help in tackling issues like ambiguity. However, capturing context in such data is extremely hard. Furthermore, LI faces the problem of accessible code-mixed dataset to build and evaluate the learning models. The bottleneck of data crisis affects the performance of systems quite a lot, generally because of the issue of over-fitting. India being a multilingual country has a rich heritage of languages and Kannada is one of the Dravidian languages as well as the official language of Karnataka state. People of Karnataka read, write and speak Kannada but many find it difficult to use Kannada script to post messages or comments on social media. While, technological limitations like keyboards of computers and smartphones is one reason, another reason may be the complexity of framing words with consonant conjuncts (vattakshara in Kannada). Hence, most of them use only Roman script or a combination of both Kannada and Roman script to post comments on social media. Kn-En code-mixed text on social media is increasing rapidly. Identifying the language of the words in code-mixed social media text is not only interesting but also challenging. LI at word level, is a sequence labeling problem where each and every word in a sentence is tagged with one of the languages in the predefined set of languages. Sequence labeling problem is a special case of Text Classification (TC). Based on ML, DL and TL, this paper explores Learning Approaches for Code-mixed LI (LA-CoLI) at word level for code-mixed Kn-En text. This study includes: • Developing learning models, namely, CoLI-vectors and CoLI-ngrams based on ML, CoLI-BiLSTM based on DL and CoLI-ULMFiT based on TL approaches, • Developing a code-mixed Kn-En annotated dataset for LI task at word level called as CoLI-Kenglish, • Creating code-mixed Kn-En embeddings for each word by merging word, sub-words and char vectors to build a Skipgram model which will be used as features in learning models to determine the efficiency of combination of vectors in ML and DL approaches, • Training a general domain Language Model (LM) using raw code-mixed Kn-En texts for ULMFiT model.
Comments in Kannada YouTube videos are used to create code-mixed Kn-En annotated dataset, code-mixed Kn-En word embeddings and train the LM. Kn-En annotated dataset and Kn-En word embeddings which will be released on request for research purpose. Overall results illustrate the competitive performance among the learning approaches.

RELATED WORK
In the ongoing history, a lot of works have been explored on code-mixed data of various language pairs for various applications such as LI, Part-of-Speech (POS) tagging etc. Soumil et al. [1] introduced a novel design for LI of code-mixed Bengali-English (Bn-En) and Hindi-English (Hi-En) data using context information. Their dataset consists of 6000 instances each selected from the datasets prepared by Mandal et al. [2] and Patra et al. [3] for Bn-En and Hi-En language pairs respectively. They performed multichannel neural associations merging CNN and LSTM coupled with BiLSTM-CRF for word-level LI of code-mixed data to achieve 93.28% and 93.32% accuracies on the test sets of two language pairs. A novel strategy for incremental POS tagging of code-mixed Spanish/English corpus is proposed by Paul et al. [4]. Utilizing dynamic model switching to get an indicator function which emits term-by-term LI tags, their baseline framework obtained an overall accuracy of 77.27%. The indicator function also regulates the output and picks the most reasonable tagging model to use for a given term. Nguyen et al. [5] introduced experiments on LI of individual words in multilingual conversational data crawled from one of the biggest online networks in Netherlands for Turkish-Dutch speakers during May 2006 to October 2012. Albeit Dutch and Turkish language words rule the discussion, English fixed phrases (e.g. 'no comment', 'come on') are incidentally observed. They evaluated strategies from different points of view on how language recognizable proof at word level can be utilized to analyze multilingual data. The highly informal spelling in online conversations and the events of named substances was used as test set. For their experiments with multilingual online conversations, they first tag the language of individual words utilizing language models and dictionaries and then incorporate context to improve the performance and achieved an accuracy of 98%. Results uncover that language models are more robust than dictionaries and adding context improves the performance. Sarkar et al. [6] proposed a Hidden Markov Model (HMM) dependent POS tagger for code-mixed Bengali-English (Bn-En), Hindi-English (Hn-En) and Tamil-English (Ta-En) shared task datasets of ICON 2015 1 . They used information from dictionary based methodologies and some word level features to additionally improve the observation probabilities for prediction. Their framework obtained an average overall accuracy (averaged over all three language sets) of 75.60% in constrained mode and 70.65% in unconstrained mode. Yashvardhan et al. [7] presents the methodologies to classify Dravidian code-mixed comments according to their polarity in the evaluation of the track 'Sentiment Analysis for Dravidian Languages in Code-Mixed Text' organized by the Forum of Information Retrieval Evaluation (FIRE) 2020 2 . They trained, validated, and tested the model using the Tamil [8] and Malayalam [9] code-mixed datasets provided by the organizers. Tamil code-mixed dataset consists of 11335 comments for the train set, 1260 for the validation set and 3149 comments for testing the model. Malayalam code-mix dataset consists of 4851 comments for training, 541 for validating, and 1348 for testing the model. Using Long Short-Term Memory (LSTM) network alongside language-explicit pre-processing and sub-word level portrayal to catch the assumption of the content, they obtained F1-scores of 0.61 and 0.60 and overall ranks of 5 and 12 for Tamil and Malayalam datasets respectively.

Construction of Dataset and Tools
This section describes the functionality used for data collection, preprocessing, training code-mixed word embeddings and building the first ever code-mixed LM for Kn-En language pairs.

Data Collection
Data is the most important part of any study and data for NLP tasks are in form of text and speech. As code-mixed text in Kn-En language pair is required for the proposed work, an efficient module that can scrap data from various sources such as social media platforms, online shopping website, etc. is required. youtube-comment-downloader 3 is modified to download 100000 comments from 373 Kannada YouTube videos which amounts to 72815 sentences after preprocessing. The comments were written only in Kannada or only in English or a combination of Kannada and English and in few cases in other languages namely, Hindi, Telugu and Tamil in addition to Kannada or English or both. However, the script of these comments is either Roman or Kannada or a combination of Roman and Kannada. The workflow of data collection module is shown in Figure 1. Data collection module accepts a list of Kannada YouTube video ids as input, downloads the comments, preprocesses them and provides as output a list of sentences extracted from the comments posted on each video.

Preprocessing
Comments in social media are unstructured, messy, contain incomplete sentences and words in short forms in addition to code-mixing of two or more languages. All these features increase the complexity of analyzing code-mixed text.
Hence, the first step in analyzing these texts is preprocessing, which includes removing duplicate comments, comments in Kannada script, short comments (less than 3 words) and comments consisting of only English words, emojis and unprintable characters. After preprocessing, roughly 90% of the data is used as raw data to train Kn-En tokenizer, code-mixed Kn-En word embeddings and code-mixed LM for Kn-En language pairs. Remaining 10% of the data is processed further to create annotated dataset for LI, at the word level. The major problem faced in analyzing code-mixed text is lack of normalization of words.

Creation of CoLI-Kenglish Dataset
A small portion (10%) of the preprocessed code-mixed texts are selected randomly and tokenized into words. These words are tagged manually by two native Kannada speakers (these people are trained about concepts of code-mixed texts and LI task) to generate CoLI-Kenglish dataset. 19432 unique words extracted from nearly 7000 sentences are categorized into 6 classes namely, 'Kannada', 'English', 'Mixed-language', 'Name', 'Location' and 'Other'. While the first two classes represent Kannada and English words respectively, 'Mixed-language' class represents word created using a combination of Kannada and English in any order. 'Name' class represents the names of persons and 'Location' class the names of locations or places. Any other words are represented as 'Other' class. The words described by 'Mixed-language' pose a real challenge to LI task as these words are framed by various combinations of English/Kannada words and Kannada/English affixes (prefixes and suffices). Beauty and also the complexity of these mixed-language words is that the word pattern depends on an individual and users posting comments on social media is increasing day-by-day. Description and samples of tokens are given in Table 1.

Word Embeddings
Word embeddings are seen as the key ingredient for many NLP tasks and has been proved as an efficient representation for characterizing the statistical properties of natural languages [10]. In addition to providing text to numeric vector conversion that is understandable to Neural Networks (NN), they model the complex characteristics of words, such as syntax and semantics which vary across linguistic contexts. Word embeddings consisting of word, sub-words, and char vectors are trained on 90% of the preprocessed Kn-En code-mixed raw data which is in the form of sentences. The steps to train the vectors as follows: • The sizes of the vector's dimensions selected for the proposed word embeddings are set based on the average unique tokens of each type in the dataset (words, sub-words, and characters). A sentence is made up of several words and each word can be decomposed into several sub-words and several characters. Hence, a word vector is extended by sub-words vectors and character vectors. In order to have a fixed length vector representation for words, the number of sub-words is fixed as the maximum of the number of sub-words of all the words in the vocabulary and similarly the number of characters is fixed as the maximum of the number of characters of all the words. For each word, word2vec, sub-word2vec, and char2vec Skipgram based models are trained as mentioned above. As the number of sub-words is not the same for all words, sub-word2vec of a word is padded with zeros depending on the difference between the maximum of the number of sub-words of all the words and the number of sub-words in a word. Similarly, char2vec is padded with zeros depending on the difference between the maximum of the number of characters of the words and number of characters in a word. Finally, word2vec, sub-word2vec, and char2vec vectors are merged together to obtain one vector for each word as shown in Figure 2. Table 2 gives a glimpse of the size of the all vectors used to obtain a vector for a word.

Kn-En Tokenizer
Tokenization is an initial but very crucial step in many token level classification tasks such as POS [11], Named Entity Recognition (NER), and token level LI [12]. Many pre-trained tokenizers are available in NLTK 6 and iNLTK 7 libraries for tokenizing Indian languages but tokenizers for code-mixed text are rarely found. SentencePiece 8 is an unsupervised text tokenizer that utilizes sub-words units e.g., Byte-Pair-Encoding (BPE) [13] and unigrams [14] with the extension of directly training from raw sentences. A Kn-En code-mixed tokenizer is trained on 90% of the preprocessed Kn-En code-mixed raw texts with a vocabulary size of 10000 using SentencePiece tools. Figure 3 illustrates the procedure of training Kn-En tokenizer and generating vocabulary.  Total size of word2vec = size of word2vec + 8 x size of sub-word2vec + 10 size of char2vec = 200 + 8x100 + 10x30 = 1300 If a word 'w' has 5 sub-words and 6 characters, then word vector will be a combination of word2vec + 5 sub-word2vec + (8-5) sub-word2vec zero paddings + 6 char2vec + (10-6) char2vec zero paddings Sub-word2vec zero paddings will be of the size of sub-word2vec and char2vec zero paddings will be of the size of char2vec.

Language Model
LM is a probability distribution over the sequence of words, in other words, LM is able to predict next word(s) in a given sequence of words and window [15]. It has applications in NLP tasks such as "Smart Compose" feature in Gmail that suggests next words in sequence. Voice to text conversion, speech recognition, sentiment analysis, text summarization, and spell correction are other NLP tasks where LMs can be used. Further, a LM can be seen as a statistical tool that can learn and analyze the natural languages' patterns. LM has got more attention with TL where the knowledge of one (source) model is transferred to another (target) model. Raw text collected from YouTube video comments (as mentioned in section 3.1.1) have been used to train a tokenizer using SentencePiece library as explained

N-grams Model
One of the challenges of LI is the structure of words in natural language. For example, it is very common in English to see the letter "q" to be followed by letter 'u' in words such as question, quarrel, qualifications, quietness, etc. However, this rule is not followed in many code-mixed texts. Since one of the primary advantages of character n-grams is language independence [16] it can be utilized for any language including code-mixed texts to capture the structure of words that has been written in a different script. In this study, a feature engineering module that generates a feature set for a given text is implemented. The feature set consists of prefixes and suffixes of length 1, 2 and 3 along with char ngrams (n = 2, 3, 5) from words, and char ngrams (n = 1, 2, 3) from sub-words.

Learning Models
Four learning models, namely, CoLI-ngrams, CoLI-vectors, CoLI-BiLSTM, and CoLI-ULMFiT are proposed for the Kn-En code-mixed LI task at word level. The learning models based on ML, DL, and TL approaches are constructed and evaluated using CoLI-Kenglish dataset and the tools constructed as mentioned above. All the four learning models are explained below:

CoLI-ngrams
This model is an ensemble of three ML classifiers namely, Linear SVC (LSVC), Multi-Layer Perceptron (MLP) and Logistic Regression (LR) with 'soft' voting. Values of the parameters used in these classifiers are given in Table  3. Figure 4 presents the structure of CoLI-ngrams model which is fed with count vectors of ngrams obtained from a feature engineering module described in section 3.1.7. Char ngrams from sub-words are extracted in two steps: i) extracting sub-words from words using BPEmb and ii) generating char ngrams for extracted sub-words. BPEmb provides pre-trained sub-words embeddings for 275 languages that are trained on texts from Wikipedia [17]. An embedding with a vocabulary size of 10000 is downloaded for English language to encode and extract sub-words from code-mixing text which helps to extract exact English words from code-mixed words. In code-mixed words, one part of the word may be an English word and rest can be Kannada suffix or prefix or with some characters which do not have any meaning in any language. In other words, sub-Words help in the generation of words that are rarely been seen in training set. Table 4 illustrates the samples of words and corresponding features generated for CoLI-ngrams.

CoLI-vectors
This model uses estimators as in CoLI-ngrams model but trained on vectors for words in the training set by utilizing embedding module that generates word embed-dings for words, sub-words and characters from raw text as discussed earlier. The purpose of developing CoLI-vectors model is to compare the performances of voting classifiers with different features and also to compare the efficiency of proposed word embedding architecture using ML and DL approaches. Figure 5 gives the structure of CoLI-vectors model.

CoLI-BiLSTM
Learning models based on DL approach have excelled conventional models based on ML approach in various NLP tasks, such as Sentiment Analysis, NER etc. [18]. CoLI-BiLSTM model is a sequence processing model based on Bidirectional Long Short Term Memory (BiLSTM) architecture. It utilizes the feature vectors obtained from proposed word embedding model. A BiLSTM comprises of two LSTMs that take the input in forward as well as in backward direction. In other words, at every time step BiLSTM networks have both backward and forward information about the sequence [19,20]. CoLI-BiLSTM model consists of layers summarized in Table 5. It includes input and embedding layers to load training data and weights from word embedding model and a BiLSTM layer followed by time_distributed layer. The purpose of using time_distributed layer is to keep one-to-one relations on input and output on RNNs including LSTM and BiLSTM. This scenario is commonly used in NNs in sequence classification tasks such as POS, NER, etc. The structure of CoLI-BiLSTM model is shown in Figure 6.

CoLI-ULMFiT
The approach of transferring knowledge of one model called source model to improve the performance of the other model called target model is called TL. Universal Language Model Fine-Tuning (ULMFiT) is one of architectures that utilize the concept of TL [21]. It consists of training a LM and then transferring the obtained knowledge and fine-tuning the target model with the dataset provided for the given task. Usually in NLP tasks the data used for training LM will be a large corpus with same or different domain from the dataset used for target task. The benefit of training a LM is that once a pre-trained LM is ready, its knowledge can be utilized in different NLP tasks including token level or text level classification, summarization, etc. A pre-trained LM understands the general features of language and then fine-tuning the LM using target task dataset helps in obtaining more properties of specific task. Following the ULMFiT architecture adopted from [15], CoLI-ULMFiT model includes training a LM from preprocessed code-mixed texts (section 3.1.2), transferring and then fine-tuning the weights using training set (section 3.1.3-CoLI-Kenglish) and finally using weights and knowledge obtained from LM in target LI model. Figure 7 presents the overview of CoLI-ULMFiT model. Fast.ai library provides necessary modules for the implementation of ULMFiT model. text.models tools from Fast.ai library is used to construct both LM and target LI models. An encoder for Average-Stochastic Gradient Descent (SGD) Weight-Dropped LSTM (AWD-LSTM) implemented using text.models tools consists of a word embedding of size 400, 3 hidden layers and 1150 hidden activations per layer-plugged in with a decoder plus classification layers to create a TC [22].

Datasets
Inspired by [23,24,25] in utilization of YouTube code-mixed comments, CoLI-Kenglish dataset has been developed. The construction of the CoLI-Kenglish dataset for LI at word level is mentioned in Section 3.1.3 and the distribution of labels in CoLI-Kenglish dataset is shown in Figure 8. Statistics of raw data and CoLI-Kenglish dataset is summarized in Table 6. Since texts in social media generally do not follow any rules the tagged dataset is highly imbalanced which may result in less F1 score. The dataset also illustrates that nearly 44.8% words are Kannada words, about 7.5% words are Kn-En mixed language words like "Dogsgalige" (meaning 'for dogs', dogs is an English word and 'galige' is a suffix in Kannada) and about 32.32% words are English words. Approximately, 70% of the tagged dataset is used for training and remaining 30% for testing.

Results
This study provides a comparison of the performances of the proposed models for word level LI task in Kn-En codemixed texts and the results are shown in terms of macro average metrics. Kannada and English are two completely different languages in various terms such as grammar, script, structure, etc., but still performances of models are promising considering the noisiness of the data. However, it is expected in closely related languages, e.g. English-German or Spanish-Italian mixed texts, LI task will be more challenging but availability of more tools for such languages enable models to have more efficient performances. The performances of the proposed learning models in addition to the performances of individual estimators in case ML models are shown in Table 7. CoLI-Kenglish dataset for word level LI consists of 6 categories and category-wise results in terms of Precision, Recall and F1-score of all the proposed models are shown in Table 8. Further, category-wise comparison of macro average F1-scores of the proposed models is illustrated in Table 9. Results of ML models illustrate that ML classifiers (both individual and ensembled) with character ngrams and affixes outperformed the ML classifiers with proposed word embeddings.
CoLI-BiLSTM model has been trained for 200 epochs (100 epochs with batch size of 128 and 100 epochs with batch size of 64. The results illustrate that that CoLI-ngrams model based on ML approach trained on a subset of morphological features including char ngrams from words and sub-words along with affixes beats the other models. CoLI-ULMFiT model has obtained relatively good performance except for words belonging to "location" class which is due to lack of sufficient samples in tagged dataset. Training an LM for ULMFiT architecture efficiently requires very huge dataset.

Conclusions
This study explores four learning models, the CoLI-ngrams, the CoLI-vectors, the CoLI-BiLSTM and the CoLI-ULMFiT for Kn-En code-mixed LI at the word level. While CoLI-ngrams and CoLI-vectors are based on ML approaches, CoLI-BiLSTM and CoLI-ULMFiT are based on DL and TL approaches respectively. Due to lack of Kn-En code-mixed dataset at word level for LI and tools to process Kn-En code-mixed data, comments in Kannada YouTube videos were scrapped and processed to construct CoLI-Kenglish tagged dataset, Kn-En code-mixed word embeddings and Kn-En code-mixed LM. CoLI-Kenglish dataset was manually tagged by Kannada speakers and grouped into six categories. Kn-En code-mixed word embeddings was constructed by merging word, sub-words, and characters vectors which were built using Skipgram model. This embedding is used as features in CoLI-vectors and CoLI-BiLSTM models and a subset of morphological characteristics are used as features in CoLI-ngrams model. CoLI-ULMFiT utilizes ULMFiT architecture to transfer the knowledge of a pre-trained Kn-En LM to a LI model. The results of the proposed models illustrate that CoLI-ngrams utilizing morphological features outperformed all other models with an average macro F1-score of 0.64. Further, CoLI-ULMFiT model also obtained similar overall performance, except for the "location" category. The results obtained by CoLI-vectors and CoLI-BiLSTM models illustrate the superiority of DL approach over ML approach in using proposed embedding. As the generated Kn-En code-mixed LM and Kn-En code-mixed word embedding can be used for other Ka-En code-mixed NLP tasks they will be released publicly along with CoLI-Kenglish dataset. In the future, it is planned to enrich both unannotated and annotated dataset, construct a balanced label distribution and to explore different feature sets and models based on different learning approaches. It is also planned, to bring morphological features into the vector space, to determine the enhancement potential of the proposed models.


extfieldsofstudy:
Computer Science
Mathematics

text:
An Optimal Mass Transport Method for Random Genetic Drift

We propose and analyze an optimal mass transport method for a random genetic drift problem driven by a Moran process under weak-selection. The continuum limit, formulated as a reaction-advection-diffusion equation known as the Kimura equation, inherits degenerate diffusion from the discrete stochastic process that conveys to the blow-up into Dirac-delta singularities hence brings great challenges to both the analytical and numerical studies. The proposed numerical method can quantitatively capture to the fullest possible extent the development of Dirac-delta singularities for genetic segregation on one hand, and preserves several sets of biologically relevant and computationally favored properties of the random genetic drift on the other. Moreover, the numerical scheme exponentially converges to the unique numerical stationary state in time at a rate independent of the mesh size up to a mesh error. Numerical evidence is given to illustrate and support these properties, and to demonstrate the spatio-temporal dynamics of random generic drift.

1. Introduction. In population genetics, genetic drift describes random fluctuations in the numbers of gene variants (alleles) over time. Allele frequency, expressed as a percentage, measures the relative fraction of an allele at a particular locus in the population, and its change quantifies the intensity of the random genetic drift [17]. Typically, when genetic drift begins, it will continue until either i) the involved allele completely disappears from the population or ii) the allele establishes permanently at 100% frequency (called fixed). In either case, genetic drift causes gene variants to disappear because infrequently occurring alleles face a greater chance of being lost in a small population, or causes a new population genetically distinct from its original population such that initially rare alleles become much more frequent and even fixed [17,32,33,37]. Both events indicate that genetic drift can decrease the population's genetic diversity, and it plays a role in the evolution of new species.
Mathematical modeling of genetic drift dates back to the pioneering works of Ronald Fisher [15] and Sewall Wright [39,40,41]. The Wright-Fisher model employs a discrete stochastic process to model dynamics of finite populations at the individual level under the restrictions that the generations do not overlap and that each copy of the gene of the new generation is selected independently and randomly from the whole gene pool of the previous generation. It is later modified and extended by Patrick Moran [28,27] (allowing generations overlap) and Komoo Kimura [24,21,22,23] (allowing the genetic mutation to spread across the population). In particular, they show that, in the limit of a large population and weak selection, these processes can be approximated by the same diffusion approximation, namely the Kimura equation [22] which describes the probability of fixation of a mutant with frequency-independent fitness. The continuum framework makes a systematic qualitative and quantitative analysis of the new model possible thanks to the tools from modern analytical and numerical analysis.
In this paper, we consider the Moran process as a paradigm and introduce its large population limits with different drift-diffusion scalings assumption [11]. Consider the dynamics of a population with N individuals that are distinguished by two neutral alleles labeled A and B such that they do not affect the survival and reproduction ability of the individual. In light of the balance between selection and drift, Traulsen et al. [35] summarized the process into three simple steps: (a) selection-an individual is randomly selected for reproduction with a probability proportional to its fitness; (b) reproduction-the selected individual produces one identical offspring; (c) replacement-the offspring replaces a randomly selected individual in the population. The process is then repeated after each time step ∆t. The fitness or reproduction rate of an individual depends on environmental conditions as well as the type and frequency of its competitors. For each type, we associate a fitness function depending on the type frequency: In particular, when Ψ is constant, i.e., the fitness of an individual is genetically determined and not affected by interactions, one recovers the classical frequency-independent Moran process.
Under the assumption N −1 ∝ (∆t) µ and a suitable time rescaling, as N → ∞, one collects from [11,23,35] the following thermodynamical limit for the density f (x, t) of allele A when ν = µ = 1 as two special cases such that genetic drift is the only evolutionary force in (1.1a), and the evolutionary force includes genetic mutation, migration and selection in (1.1b). They correspond to the limits of discrete process with ν > µ = 1 2 and µ ∈ ( 1 2 , 1], ν = 1 − µ, respectively. The biological significance urges us to impose the non-flux boundary condition to each of (1.1) such that the following conservation holds where θ(x) is the fixation probability function that satisfies and it can be explicitly given by One notes that (1.4) recovers the conservation of total population for (1.1) with θ(x) ≡ 1, and conservation of mass center for (1.1a) with θ(x) = x.
In this work, we will take advantage of the fact that the following free energy admits (1.1) as a gradient flow with respect to a variation of optimal transport distances (see e.g. [36,6,4,10]) and that (1.1) can be rewritten as denotes the first variation of the energy E with a fixed mass constraint. We now recall several theoretically relevant results on the well-posedness of (1.1)-(1.2) by [12,11,34]. Let us denote by BM + ([0, 1]) the space of functions with positive Radon measure in [0,1].
According to Theorem 1.1, we always expect (1.1) to collapse into the Dirac-delta singularity regardless of the initial distribution. The long-time dynamics along with this spiky spatial profile well demonstrate and capture the formation of gene segregation, i.e., an asymptotic gene fixation with allele A (resp. allele B) occupying the whole population when L 0 f 0 (x)θ(x)dx = 0 (resp. 1).
More recently, the gradient flow structure (1.7) has been discussed in [10]. The analysis of the model is built upon the classical steepest descent variational schemes for nonlinear Fokker-Planck equations which are introduced in [20] and generalized in [1,2]. The transport distance between probability measures has to be adapted to the degenerate diffusion coefficient of (1.1). They introduced a suitable distance called the Shahshahani distance showing the convergence of the variational scheme for the specific case of replicator dynamics, i.e. (1.1) without the degenerate diffusion term. It is an open problem to show the convergence of this variational scheme including the degenerate diffusion term.
The spatial-temporal dynamics of the Kimura equation are well understood in the purely diffusive case. However, the singularity induced by the convergence towards Dirac-delta concentrations, which well model the biologically realistic gene segregation or fixation through the drift, imposes substantial challenges to their numerical approximation, in particular for long-time computations. Therefore, one of the crucial features of any numerical scheme to solve (1.7), while keeping its biological relevance, is to accurately capture the concentration phenomena at the discrete level. Lagrangian numerical schemes for these gradient flow interpretations of one-dimensional Fokker-Planck, thin-film and quantum drift equations have been devised by various authors, see e.g. [19,18,2,38,9,8,25,31,26,30]. We here build upon the strategy of solving an equivalent equation satisfied by the diffeomorphisms mapping the initial data to the solution at later times as performed in [2,7,8,3], see next section for details. We refer to [5] for a very recent review on the state of the art of this kind of numerical schemes. The main advantage of the present approach is being able to deal with Dirac-delta concentrations easily and rigorously now that we work in mass variables. This merit is also utilized by [2,8,3] to study blow-ups in Keller-Segel models for chemotactic cellular aggregation.
Other numerical schemes have recently been proposed in the literature to solve the Kimura equation (1.1). In [13] the authors introduced a numerical method based on Lagrangian maps that preserves the free energy decay of the system. They also analyze a convex-concave splitting approach which leads to an implicit method, and prove the unique solvability of this method. However, the equation is solved in the original variables therein, so they have to devise heuristic criteria to capture the concentration of mass towards the endpoints. Moreover, they do not take full advantage of the gradient flow structure (1.7) and the variational schemes at the discrete level. As an alternative, our approach is to construct numerical schemes directly on the optimal maps and can describe the Dirac-delta to the fullest extent. See details in Section 2. [42] performs an interesting "horse race" comparison of a serial of finite volume and finite element schemes for (1.1a). Their critical comparison of the long-time asymptotic performance urges carefulness in choosing a numerical method for this type of problem, especially when the main properties of the model are not kept by the scheme. We also want to mention that the genetic drift problem with multi-alleles, cast as a multi-dimensional PDE, is investigated in [43] through finite-difference methods, where the authors propose a numerical scheme with absolute stability and conserve several biologically/physically relevant quantities such as positivity, total probability, and conservation.
The rest of this paper is organized as follows. Section 2 introduces the evolution problem for the map in mass variables as in [2,8,5] together with the needed prerequisites in optimal transport theory. Since the free-energy functional shows again in the new variables a non-convex structure, we present our numerical method and then apply it for the genetic-drift problem by introducing the convex-splitting technique as the preprocessing step of the Euler implicit scheme. Section 3 is devoted to the analysis of the properties of the numerical scheme. We show that the discrete problem converges exponentially fast to a unique stationary state with a monotone drift V for any fixed discretization parameter modulo error terms in spatial discretization, and we characterize all possible stationary states and limit points otherwise. See Theorem 3.9 and Theorem 3.12 for our main results. Finally, Section 4 conducts several sets of numerical experiments to validate the shown properties and to accurately illustrate the long-time dynamics of random genetic drift.
2. Numerical Methods. In this section, we recast the random generic drift models into an evolution problem for map of the mass variable, and then propose a numerical scheme to solve the new equation. Since the associated new free-energy functional is non-convex, we introduce the convex-splitting technique as the preprocessing step of the Euler implicit scheme in its numerical solver.
2.1. Optimal transport and the Wasserstein distance. We first introduce a suitable Wasserstein distance in the probability space P([0, 1]) such that equation (1.1) can be interpreted as a gradient flow of the free energy (1.6). Due to the presence of the variable coefficient x(1 − x) in (1.1), the quadratic Wasserstein distance will not be based on the usual Euclidean distance but on the induced generalized Shahshahani distance dt, for x, y ∈ Ω = (0, 1). Chalub et al. [10,Lemma 10] prove that the infimum in (2.1) is achieved at a unique constant-speed geodesic, and d can be uniformly extended toΩ ×Ω as moreover, d defines a distance inΩ and the metric space (Ω, d) is Polish. This Shahshahani distance is locally equivalent to the Euclidean one in the interior, but behaves differently close to the boundary. This difference is reflected by the dynamics of replicator-diffusion equation (1.1), which is locally uniformly parabolic in the interior, but degenerate at the boundaries. We now present several concepts from transport theory essential for this paper. Let µ and ν be two absolutely continuous measures with respect to the Lebesgue in P([0, 1]), and T be a measurable map from [0, 1] → [0, 1]. We say that T transports µ onto ν and write ν = T µ if ν(B) = µ • T −1 (B) for any measurable set B ⊂ [0, 1]. We also say ν is the push-forward or the image measure of µ by T . To introduce the corresponding Wasserstein distance between µ and ν, one defines as soon as the source measure µ has no atoms. In fact, by Brenier's theorem, if µ is absolutely continuous with respect to the Lebesgue measure, then there exists a measurable nondecreasing map T such that ν = T µ. The proper definition of the associated Wasserstein distance needs a relaxed variational scheme given by where Π runs over the set of transference plans Γ between µ and ν. One important simplification in 1D is that the optimal plan can be characterized fully in terms of the inverse of cumulative distribution functions, see [36]. To be specific, let F and G be the cumulative of the 1D functions f and g, and define the pseudo-inverse By Brenier's theorem and the definition of the image measure, we have Then it is straightforward to obtain ϕ = G −1 • F , and the Wasserstein distance in (2.2) becomes The proof of this fact without relying to Brenier's theorem can be found in [36, Section 2.2]. By the definition of pseudo-inverse function, we can derive the following evolution equation satisfied by Φ(η, t) We would like to remark that, the Dirichlet boundary condition of Φ(1, t) = 1 applies here since one have from the strong maximum principle that f (x, t) in strictly positive in [0, 1] for all time t if f 0 (x) ≥, ≡ 0. This, on the other hand, indicates that Φ(η, t) must be strictly increasing in η, which is an important property to preserve for numerical schemes. Actually, according to the non-flux boundary condition of f (x, t), we have This implies the Dirichlet boundary condition holds Φ(0, t) = 0 and Φ(1, t) = 1 for all time t. However, a free boundary condition should be adopted if f (x, t) remains compactly supported when studying problems with diffusion degenerate inside the domain. Note that (2.4) has the following free energy Then we will connect this evolution problem for the map Φ(η, t) pushing forward the initial data f 0 to the solution f (x, t) at time t with the continuum limit of implicit Euler steps obtained as Euler-Lagrange conditions for suitable variational problems. Finally, let us restate Theorem 1.1 in terms of the map Φ(η, t).
In the sequel, we will design a numerical scheme capable of accurately capturing the long-time behavior described in Theorem 2.1.

Discretization for Euler-Lagrange of Replicator-Diffusion (1.1).
We now consider the spatio-temporal discretization for the Euler-Lagrange problem (2.4) of the full Replicator-Diffusion equation (1.1). Throughout this paper, we assume that the time step τ and space step h are constant in the discretization. In terms of the Wasserstein distance (2.3), the Jordan-Kinderlehrer-Otto (JKO) steepest descent scheme implies that finding the inverse distribution function in (2.4) corresponds to solve the following for a fixed time step τ > 0 In light of (2.2), we find that Hence an approximated Euler-Lagrange equation (2.4) associated to this minimization problem is If we denote Φ k i = Φ(ih, kτ ) for i = 0, 1, · · · , N and N h = 1, k ∈ N, our full finite difference discretization of (2.6) is the following implicit scheme: with the Dirichlet boundary condition Φ k+1 The solution at each time step is computed by an iterative Newton's procedure with initial diffeomorphism obtained by a preprocessing step, and one does not need a CFL condition for this implicit in time discretization.

Preprocessing
Step via a Convex Splitting Technique. System (2.4) can be viewed as a gradient flow associated with the energy functional (2.5). This energy has a non-convex structure, and it prohibits the direct application of the proposed implicit scheme due to the singularity of the numerical scheme on the boundary and local convergence of Newton's method. However, the merits brought by the convex-splitting technique enable us to construct a numerical scheme in the preprocessing step and then calculate the initial diffeomorphism for Newton's method in the Euler implicit scheme.
To this end, we write V = V c − V e for V c and V e being smooth convex functions and then apply the convex splitting method of Eyre [14] to obtain the following semi-discrete scheme Here, the map from Φ k to Φ k+1 is an optimal transport in the sense that Φ k+1 minimizes the functional where W is a convex functional explicitly given by Let us introduce the discrete space domain and its closureQ := Q ∂Q with boundary The full finite difference discretization of (2.8) is formulated as follows: with the boundary condition Φ k 0 = 0 and Φ k N = 1 for each k. At each time step, we solve a system of nonlinear equations by a damped Newton's iteration. For notational simplicity, let us denote the nonlinear functional F in (2.11) as then we will calculate the Jacobian matrix DF of (2.12) and determine the Newton update γ k+1,n+1 by where index k corresponds to the temporal discretization and n to the Newton iteration.
3. Numerical Analysis. This section analyzes the implicit Euler scheme (2.7) and the convex splitting scheme (2.11). We first collect several important properties of the splitting technique applied in [13] that includes the unique solvability, convergence of the Damped Newton method (2.13) and dissipation of energy. Then we prove that there exists a lower bound for the discrete energy, which implies the existence of numerical solutions. Finally, in light of the fact that any steady state must be a Heaviside-type step function, we show the convergence of of the implicit scheme (2.7).

Several Properties of the Convex Splitting Scheme.
Lemma 3.1. The numerical scheme (2.11) is uniquely solvable in space Q given by (2.10). Proof. We first introduce the discretization of (2.9) as is a convex function on the closed convex setQ and J N (y) = +∞ on boundary ∂Q, it is straightforward to obtain that there exists a unique minimizer x ∈ Q.
To show the unique solvability of the scheme, it suffices to prove that x ∈ Q is the minimizer of J N (y) if and only if it is a solution of (2.11). Suppose that x ∈ Q minimizes J N (y). Define j 1 ( ) := J N (x + (y − x)). Then there exists 0 > 0 small enough such that x + (y − x) ∈ Q for any ( , y) ∈ (− 0 , 0 ) ×Q. Now that j 1 ( ) achieves its minimizer at zero, one has j 1 (0) = 0 hence h Therefore, x is a solution of (2.11).
To prove the "only if" part, we assume that x ∈ Q solves the scheme (2.11) for any y ∈ Q. Then where the last inequality holds since m − ln m − 1 > 0 for any positive m, and V c is convex function with positive second derivative. This completes the proof.
We now introduce the following concept to prove the convergence of the Newton's method of the scheme (2.11).
Definition 3.2. ( [29]) Let G be a finite-dimensional real vector space and Q be an open nonempty convex subset of G. Then a convex function Λ ∈ C 3 : Q → R is called self-concordant on Q if there exists a constant a 0 > 0 such that the following inequality holds for all x ∈ Q and all u ∈ G: where (D k Λ(x)[u 1 , · · · , u k ] is its k-th differential taken at x alone the collection of direction (u 1 , · · · , u k ).
Then we have the following theorem.
Since both the linear and quadratic functions have zero third-order derivatives, one can easily find that J 1 N (y) is selfconcordant for any a 0 in Q.
We proceed to prove that j 2 (ξ) := J 2 N (y + ξu) is a self-concordant function of ξ along every line u in Q. To this end, we have from direct calculations that Since V c (x) is convex and smooth, for any y i ∈ Q there exists a constant M v > 0 such that Let us define then we proceed to find that where the fourth inequality follows from the following Now, let us choose with M v satisfies (3.2), then J 2 N (y) is self-concordant for y ∈ Q with parameter a 0 . This implies that J N (y) is selfconcordant and the Newton's iteration is convergent thanks to Theorem 2.2.3 in [29].
Let {Φ i } i=0,1,...,N be any strictly increasing sequence with Φ 0 = 0 and Φ N = 1. Then the discrete free energy functional is defined as Then the following theorem states that the energy dissipation is preserved through the discretization.
Lemma 3.4. The discrete energy dissipation for the evolution of the discrete energy (3.4) holds for the scheme (2.11).
Proof. We rewrite According to the inequality ln( x y ) ≤ (x − y) 1 y and convexity of V c , we can obtain Then it follows from the summation-by-parts that Similarly, in light of the convexity of E 2 N we find that On the other hand, let us rewrite scheme (2.11) into the following then one infers from (3.6)-(3.7) that which is claimed.

3.2.
Convergence of the Euler Implicit Scheme. We proceed to study convergence of the discrete solution of (2.7) to a unique steady state in the long time. For this purpose we first prove the existence of this steady state. and our argument begins with the following observation. Proof. The trick of the proof is to find the "median" of the sequence. Let k 0 be the integer such that Therefore, one finds that These inequalities, together with the uniform boundedness |V (x)| ≤ M in [0, 1], enable us to estimate where we apply the fact 0 < Φ k 0 +1 − Φ k 0 < 1 for the second inequality and the Jensen inequality for the third one.
As a consequence we establish the existence of a stationary state realized as a minimizer of the free energy functional. Notice that lower bound of the discrete energy functional (3.5) diverges to −∞ as h → 0, and this is consistent with the fact that the continuous energy functional (2.5) is not bounded from below. Now, we are ready to show that the steady state of (2.4) must be a Heaviside step function.
Lemma 3.6. Denote {Φ i } i=0,1,...,N as the discrete solution of (2.7). Define {F k i } k∈N i=0,1,··· ,N as Then for any fix time step size τ > 0, the following statements hold: (i) the conservation law holds for the scheme (2.7) as follows If V (·) does not change sign in [0, 1], the steady state of the numerical scheme (2.7) must be {Φ * i } i=0,1,...,N . Proof. Summing over i, one can find from (2.7) that which implies that the conservation law in (i) holds for the scheme. Suppose that V (·) is of one sign in [0, 1]. Let {Φ * i } i=0,1,...,N be an equilibrium of scheme (2.7) such that Let m 1 be the largest spatial index for Φ * m1 = 0 and m 2 be the smallest for Φ * m2+2 = 1. Then (3.12) implies ) ≤ 0 leads a contradiction to V (·) > 0 since we then must have Φ * i ≡ 1 for i = m 1 + 1, ..., N . One also gets a contradiction if V (·) < 0. Therefore, the definitions of m 1 and m 2 imply that in either case the steady state is a nondecreasing step function from zero to one, whereas the jump location is determined by the conservation law (3.9). Now that {F k i } k∈N i=0,1,··· ,N is bounded in time thanks to (3.8), the Bolzano-Weierstrass theorem finds a subsequence convergence to {F * i } i=0,1,...,N and then the steady state can be written as (3.10) with a = 1 and m determined by (3.11). Furthermore, for any > 0, we can choose τ small enough such that for all n 2 > n 1 > 1 which implies that {F k i } k∈N i=0,1,··· ,N is Cauchy sequence with respect to k. Therefore, as τ goes to zero, we have the unique existence of {F * i } i=0,1,...,N , which implies the existence and uniqueness of the steady state as expected. Remark 3.7. For diffusion equation (1.1a), (3.13) and (3.14) imply that m 1 = m 2 as V ≡ 0, hence there exists a ∈ (0, 1) such that the steady state is given by (3.10). Moreover, the conservation of mass center within the discrete scheme implies that it the integer-valued function. The assumption V (·) being one sign is technical and our numerics suggest that the conclusion still holds if otherwise.
..,N be the discrete solution of (2.7). Define {F k i } k∈N i=0,1,··· ,N as follows Then for any smooth function V (·) the following statements are true: (i) the conservation law (3.9) holds for the scheme (2.7); Proof. The proof of the first statement is quite similar to that for Lemma 3.6 and it is omitted. To show that (3.10) under (3.11) is a steady state of (2.7), we first recall that the discrete function {Φ * i } i=0,1,...,N is an equilibrium if and In order to show that (3.10) supplemented by (3.11) is a steady state of the scheme (2.7), one only needs to verify the condition of Φ * m+1 . If V (·) is strictly positive or negative, (3.11) implies that a = 0 or a = 1. If V (·) = 0, the steady state is uniquely determined by (3.15). If V (·) changes sign in [0, 1], then (2.7) implies Substituting Φ * m = 0, Φ * m+1 = a and Φ * m+2 = 1 into (3.16) gives that V (a) = 0, which is expected. If {F k i } k∈N i=0,1,··· ,N has the limit {F * i } i=0,1,...,N , the parameter m can be determined by (3.11). In either case, we prove the statement in (ii). We next analyze the long time properties of the numerical scheme concerning its convergence to the steady state. We will first show the exponential decay for the diffusion equation (1.1a) for any fixed time step size modulo error terms. As for the replicator-diffusion equation (1.1) we will be able to obtain a similar result in the case of monotone drift fitness potentials.
for h small enough and C 0 both depending only on the distance of the initial data to the steady state Φ 0 − Φ * , that is, the fully discrete numerical scheme converges exponentially to the steady state in time modulo O(h) terms.
Proof. To show the convergence to the steady state, we recall from Lemma 3.6 that the discrete solution takes the form 0, for i = 0, 1, · · · , m, a, for i = m + 1, 1, for i = m + 2, · · · , N, with a and m to be determined. The time evolution of the L 2 -distance between Φ k and the stationary state Φ * reads To estimate A n , B n and C n , we apply (2.7) with V (·) = 0 and κ = 2 and deduce from the summation-by-part that where we skip the index k + 1 for simplicity. On the other hand, the conservation of mass center implies that Therefore, for any fixed τ > 0, there exists some k 0 ∈ N such that b k+1 ≤ 2 for any k > k 0 , and Finally, let us choose the space step small such that h ≤ Φ 0 − Φ * 2 /5. Then one concludes that 1 1+2τ k ≈ exp(−2kτ ) ≈ exp(−2t) for τ small, and this completes the proof.
Remark 3.10. When h goes to zero, we can deduce a posteriori uniqueness of the equilibrium from the convergence property. Indeed, letΦ * be another equilibrium state and set Φ k = Φ k+1 =Φ * in (3.18) with h = 0. Then we have from the evolution that Φ * − Φ * 2 ≤ 0 which proves the uniqueness of the stationary solution.
Remark 3.11. In the continuum case one can linearize the equation around the steady states (i.e., the Heaviside step function) and collects the following eigenvalue problem (3.19) −ϕ = λ x(1−x) ϕ, x ∈ (0, 1), ϕ(0) = ϕ(1) = 0, the principle eigenvalue of which determines the exponential convergence rate of (1.1) to the unique steady state (e.g., Theorem 3 in [11]). One notes that λ in (3.19) has a Rayleigh's variational quotient as and Theorem 262 in [16] implies that λ ≥ 2, with the principal eigenvalue λ = 2 achieved at ϕ(x) = Cx(1 − x) for some C ∈ R. This indicates that the exponential convergence rate in time of our discrete scheme in Theorem 3.9 is optimal since the decay rate cannot be better for generic initial data. Our numerical simulations in the coming section support this conclusion.
The second main result of this paper goes as follows.
Theorem 3.12. The solution of the numerical scheme (2.7) for the Replicator-Diffusion equation (1.1) with a monotone fitness potential (V either non-positive or non-negative) converges exponentially to the steady state of (1.1) with exponential rate κ (the diffusion rate) as the time step size τ goes to zero modulo O(h) terms. More precisely, the numerical scheme (2.7) for the general replicator-dynamics equation (1.1) with a monotone fitness potential satisfies for h small enough and C 0 both depending only on the distance of the initial data to the steady state Φ 0 − Φ * and the bound of V .
Proof. First of all, Lemma 3.8 implies that the sequence given by (3.10)-(3.11) is an equilibrium of (2.7). With |V (·)| ≤ M 1 in [0, 1] for some constant M 1 > 0, the same calculations in (3.17) lead to Therefore, one finds In light of the concavity of x(1 − x), we can obtain from Taylor's expansion that We next show the convergences of , then summing equation (2.7) with respect to i yields The fact that 0 ≤ c 1 k ≤ m + 1 for each k ∈ N implies that {c 1 k } is Cauchy with the limit a. Similarly, since 0 ≤ d k ≤ N − 1 and d k+1 − d k < 0, one can show that is Cauchy and its limit is (N + a − m − 2). Case ii). If V (·) < 0. Then we denote c 2 k+1 := , then summing equation (2.7) with respect to i yields The fact that 0 ≤ c 2 k ≤ N − m − 2 for each k ∈ N implies that {c 2 k } is Cauchy with the limit (N − m − 2). Similarly, since 0 ≤ d k ≤ N − 1 and d k+1 − d k > 0, one can show that is Cauchy and its limit is (N + a − m − 2), then {c 1 k } is Cauchy with the limit a.
Since c 1 k → a and d k → (N + a − m − 2) with 0 ≤ a ≤ 1, there exists some k 1 ∈ N such that for any k > k 1 Therefore from (3.21) we deduce that 1

Numerical Studies.
We now present several sets of numerical simulations to illustrate and verify the analytical properties of the proposed scheme and to demonstrate further spatio-temporal dynamics within the generic drift which are not captured by the analysis. Note that the time-dependent solution converges to the Dirac-delta function in the long-time limit. In our map variables, this corresponds to the convergence to a stationary map Φ ∞ (η) given by a step function. To deal with the issue that discrete values near the boundary are too close to be distinguished from each other under machine precision, we set the following additional boundary criteria throughout our simulations

Diffusion Equation.
We first numerically study the purely diffusive model (1.1a). The initial data are chosen to be either the symmetric f 0 (x) = π π−2 (1 − sin(πx)) or asymmetric f 0 (x) = 2x so we test the robustness of the proposed scheme with respect to initial data regarding the long-time dynamics. We remind the reader that the conserved quantity (1.4) determines asymptotically the proportion of the mass allocated to each end of the interval, that is the jump location in mass variables, see Theorem 2.1. To study the convergence of the numerical scheme, we consider the discrete L 1 error for Φ k,τ i,h , i = 0, 1, ..., N as follows is obtained with mesh size h = 1/1600 and time increment τ = 0.0005. Notice that this quantity is equivalent to an approximation of the Monge-Kantorovich or 1-Wasserstein transport distance between the associated densities. Moreover, the experimental order of convergence (EOC) in spatial-discretization h is defined as and the experimental order of convergence (EOC) in time-discretization τ is defined as EOC τ := log 2 (E τ ) − log 2 (E τ /2 ). Table 4.1 presents the errors and EOCs in spatial-discretization computed at a final time T = 0.992, a multiple of τ = 0.016, with constant time increment τ = 0.0005. We can see that difference errors decrease as the mesh is refined and a finer spatial-discretization results in higher accuracy. Besides, the orders of convergence in spatial-discretization indicate a convergence order of greater than one in both cases, and EOCs increase for a finer mesh. The results from both symmetric and asymmetric initial data indicate that the proposed scheme is robust to initial data. Table 4.2 shows the errors and EOCs in time-discretization with a fixed space mesh h = 1/1600. Similar robust convergence witnesses that the method converges as the time mesh shrinks.     1a) subject to symmetric initial data on the left and asymmetric initial data on the right. Both experiments capture the dominance of dynamics by the Dirac-delta as stated in Theorem 2.1. In particular, one finds that the degeneracy of diffusion at the endpoints makes the spatio-temporal dynamics of (1.1a) relatively simple in light of the proposed mass transport method now that the singularities can be quantified.
Next we present the evolution of (1.1a) out of the initial data f 0 (x) = π−2 π (1 − sin(πx)) and f 0 (x) = 2x with the step size h = 1/999 and time step τ = 0.001. Figure 4.1 captures the convergence to the steady-state of the numerical scheme given by a Heaviside step function, with at most one intermediate value, as stated in Lemma 3.6 or Lemma 2.1. We note that the jump locations η 0 = 1 2 and 1 3 for these two initial data, respectively, and this, together with the conservation of center of mass, are well preserved by our method. One finds that the (exponential) decay rate of the discrete scheme approaches 1 as space or time mesh shrinks. These findings agree well with Theorem 3.9.
transport distance between the associated densities. We can find that the convergence is faster than the rate obtained in Theorem 3.9 -exponential convergence in time modulo O(h) terms. This confirms the derivation of Theorem 3.9 and the existence of time modulo O(h) terms. We can also see that the convergence rates depend on τ and h, and as either one shrinks, the convergence of the L 2 -norm fits better e −t initially, which is consistent with the theoretical analysis, however, the decay rate is dramatically influenced by the O(h)-error in Theorem 3.9 as time evolves.

Replicator-Diffusion Equation.
We proceed to study the full Replicator-Diffusion equation (1.1) that includes the forward Kimura equation with frequency selection. In particular, we restrict our attention to the choice that V (x) = αx + β for some constant α and β. Throughout the simulations we fix f 0 (x) = 6x(1 − x) with the space mesh size h = 1/999 and time step size τ = 10 −3 , and then study the variation of the diffusion rate and fitness potential, starting with κ = 2, V (x) = 4x + 2 and κ = 4, V (x) = −4x + 2 as examples.
One observes from Figure 4.3-(a) that the numerical solution converges to the step function with the discrete conservation laws holding. As shown in Figure 4.3-(b), there exist some jumps in the dissipation of the entropy, which can be attributed to the additional boundary criteria (4.1). One can observe the smearing of these jumps by decreasing the tolerance in (4.1). However, this tolerance value cannot be chosen too small either for Newton's iterations to converge. Finally, Figure 4.3-(c) shows the decay in the L 2 -distance or 2-Wassertein transport distance towards equilibration. One again finds that the (exponential) decay rate of the discrete scheme approaches 1 as the space or time mesh shrinks for short times. This finding agrees well with Theorem 3.12. We also observe in Figure 4.3-(Top) that the asymmetric initial data leads to an interesting phenomena of different equilibration time scales, a first transient slower time scale in which all the mass from the left is pushed towards 0 followed by a faster decay towards the equilibrium for large times. We remark that the proofs of our main results rely on the assumption that V (x) is of one sign for technique purpose, and the numerics in the bottom of Figure 4.3 favor the application of our scheme otherwise. However, we have to exclude the relaxation of this assumption in this work, unfortunately.
4.3. Spatial-Temporal Dynamics. We conclude this section by presenting several additional sets of numerical experiments on the full model (1.1) to illustrate and verify our main results, and to capture the spatio-temporal dynamics of the random genetic drift when analytical tools lack. For these purposes, we choose featured initial data to showcase the dynamics of generic shrift to the variation of the potential function. For simplicity, we assume κ = 2 and consider two different fitness potential V (x) = 2 and V (x) = −3x + 1. The step size h = 1/999 and time step τ = 0.001 are fixed for all the numerical experiments. Each row in Figure 4.4 shows the evolution of the inverse cumulative distribution function Φ of (1.1) corresponding to the initial data depicted in each column. The dynamics are soon dominated by the Dirac-delta concentrations at the ends of the interval in each case, though one observes relatively different evolutionary transient phenomena such as merging of bumps and different equilibration speeds. All these findings point to the computational/mathematical evidence that the evolutionary changes at the molecular level are caused by random genetic drift when described by a Moran process. Finally, we report in Table 4.3 the error of the computed jump location against the theoretical one given in Theorem 2.1. They add another layer in the robustness of our proposed scheme as the location in Φ determines the portion of the alleles in the genetic aggregation.