\section{Introduction}

In the process of creating a neural network binary classifier, the network is often trained via the binary cross-entropy (BCE) objective function. The network's output is a probability value $p \in [0,1]$ that must be translated into a binary value $\{0,1\}$ indicating set membership to the positive class. The Heaviside step function $H$ is commonly used to determine set membership based on a heuristic threshold $\tau$, where $p \geq \tau$ are considered positive classification outcomes.

It is a common assumption that optimizing a desired evaluation metric is preferable to optimizing a different loss \cite{eban2016scalable, Rezatofighi_2019_CVPR, herschtal2004optimising, song2016training}. However, the Heaviside function, commonly used to compute confusion matrix set membership in terms of true positives, false positives, false negatives, true negatives, has a gradient with properties not conducive to optimization via gradient descent. The Heaviside function’s gradient is not defined at the threshold $\tau$ and is zero everywhere else (see Fig. \ref{fig:heaviside1} (left) for an example with $\tau = 0.7$). Therefore the gradient of the Heaviside function cannot be used to effectively back propagate errors and optimize a neural network according to metrics composed of confusion matrix set values, such as $F_1$-Score. To address this challenge, we propose a novel approximation of the Heaviside function that adheres to important properties of the original Heaviside function while allowing for optimization via gradient descent. Figure \ref{fig:heaviside2} (left) depicts our proposed approximation.

\begin{figure}[tb!p]
    \centering
    \begin{minipage}{.25\textwidth}
        \centering
        \includegraphics[width=0.7\linewidth]{plots/heaviside/heaviside.pdf}
    \end{minipage}%
    \begin{minipage}{0.25\textwidth}
        \centering
        \includegraphics[width=0.7\linewidth]{plots/heaviside/heaviside_sigmoid_approx.pdf}
    \end{minipage}
    \caption{Heaviside function $H$ (left) and \cite{kyurkchiev2015sigmoid} approximation (right) at $\tau=0.7$.}
    \label{fig:heaviside1}
    \vspace{-0.5em}
\end{figure}

\begin{figure}[tb!p]
    \centering
    \begin{minipage}{0.25\textwidth}
        \centering
        \includegraphics[width=0.7\linewidth]{plots/heaviside/heaviside_linear_approx_labeled.pdf}
    \end{minipage}%
    \begin{minipage}{0.25\textwidth}
        \centering
        \includegraphics[width=0.7\linewidth]{plots/heaviside/heaviside_sigmoid_fit_approx.pdf}
    \end{minipage}
    \caption{Proposed Heaviside function approximation $\mathcal{H}$ (left) and a sigmoid fit to our proposed Heaviside approximation (right) at $\tau=0.7$.}
    \label{fig:heaviside2}
    \vspace{-0.5em}
\end{figure}

The proposed approximation is flexible in that it can be used to directly optimize any metrics based on the confusion matrix sets without re-formulating the loss. One may simply substitute the step function with the proposed Heaviside approximation and directly optimize the metric. This flexibility is important because the optimal metric for binary classifiers often changes with the level of class imbalance of the data. For example, while accuracy is popular for evaluating classifiers, accuracy will over-estimate performance \cite{he2009learning} in datasets with high levels of class imbalance. In these cases, better estimators of performance include the harmonic mean of precision and recall ($F_\beta$-Score).

Though metrics based on confusion matrix set values are considered not decomposable, we find that when trained with an adequately large batch size our method performs as well as or better than the commonly used BCE objective function and other methods of directly optimizing the metric. Though this requires use of a large-enough batch size during training, we find these batch sizes to be well within limits of modern deep learning hardware when training on a variety of publicly available datasets. These datasets have varying levels of class balance, reflecting the fact that class imbalance is common in real-world data \cite{reed2001pareto}.

Our approach is scalable with the runtime equivalent to the metric's computational cost. Even for metrics that require computation over multiple thresholds, such as area under the receiver operating characteristic (AUROC) curve, we suggest a constant-time algorithm which ensures training runtime scales linearly with the number of samples in the dataset. Moreover, the proposed method places no limit on the maximum size of the batch used for training, as is the case in the adversarial approach proposed by \cite{fathony2020ap}.

In summary, our main contribution is a novel approximation of the Heaviside function to facilitate batch gradient descent when training classifiers that directly optimize evaluation metrics based on confusion-matrix sets. We explore the practical effects of batch size on non-decomposable metrics in light of our approach. Also, we show the flexibility of the approach by applying our direct optimization method to several confusion-matrix based metrics on various datasets

\section{Preliminaries}
\label{sec:preliminaries}

Once a binary classification network is trained, it is common to choose a threshold value $\tau$ to transform the output probability to a binary value. The threshold is often chosen by searching for an optimal value considering the carnality of the sets in the confusion matrix ($|TP|$, $|FP|$, $|FN|$, $|TN|$). These sets are computed using the Heaviside function $H$:
\begin{equation}
H(p, \tau) =
  \begin{cases}
    1 & p \geq \tau \\
    0 & p < \tau
  \end{cases}
\end{equation}

There are several important properties of the Heaviside function that must be considered when evaluating approximations. First, the complement of $H$,  $1-H$, is the reflection over the line $y=0.5$. Second, at the threshold $p=\tau$, the sample could be a member of either class with equal probability. In practice, to avoid ambiguity either the positive or negative class is chosen -- we chose the positive class in this work. Figure \ref{fig:heaviside1} (left) shows an example Heaviside function.

The confusion matrix set membership can then be computed for a prediction $p$ and for ground truth label $y$ via the set membership functions:
{\small
\begin{center}
    \noindent\begin{minipage}{.44\linewidth}
        $$
            tp =
              \begin{cases}
                H(p, \tau) & y = 1 \\
                0 & \text{else}
              \end{cases}
        $$
    \end{minipage}%
    \begin{minipage}{.56\linewidth}
        $$
            fp =
              \begin{cases}
                H(p, \tau) & y = 0 \\
                0 & \text{else}
              \end{cases}
        $$
    \end{minipage}
\end{center}
\begin{center}
    \noindent\begin{minipage}{.5\linewidth}
        $$
            fn =
              \begin{cases}
                1-H(p, \tau) & y = 1 \\
                0 & \text{else}
              \end{cases}
        $$
    \end{minipage}%
    \begin{minipage}{.5\linewidth}
        $$
            tn =
              \begin{cases}
                1-H(p, \tau) & y = 0 \\
                0 & \text{else}
              \end{cases}
        $$
    \end{minipage}
\end{center}
}

Consider a set of predictions $p \in [0,1]$, ground truth labels $y \in \{0,1\}$ and threshold value $\tau \in (0,1)$, e.g. $\{(p_1, y_1, \tau), (p_2, y_2, \tau), \cdots, (p_n, y_n, \tau)\}$. The cardinality of each confusion matrix set is computed as the sum over the set membership functions:
{\small
\begin{center}
    \noindent\begin{minipage}{.5\linewidth}
        $$
            |TP| = \sum_{i=1}^n tp(p_i, y_i, \tau)
        $$
    \end{minipage}%
    \begin{minipage}{.5\linewidth}
        $$
            |FP| = \sum_{i=1}^n fp(p_i, y_i, \tau)
        $$
    \end{minipage}
\end{center}
\begin{center}
    \noindent\begin{minipage}{.5\linewidth}
        $$
            |FN| = \sum_{i=1}^n fn(p_i, y_i, \tau)
        $$
    \end{minipage}%
    \begin{minipage}{.5\linewidth}
        $$
            |TN| = \sum_{i=1}^n tn(p_i, y_i, \tau)
        $$
    \end{minipage}
\end{center}
}
Common classification metrics based on the confusion matrix then build on the above cardinality sets. Two examples of these metrics are precision and recall. Precision, $|TP|/(|TP|+|FP|)$, is the proportion of positive results that are true positive results. Recall, $|TP|/(|TP|+|FN|)$, indicates the proportion of actual positives examples that are correctly identified.

Precision and recall represent a trade-off between classifier objectives; it is generally undesirable to optimize or evaluate for one while ignoring the other \cite{he2009learning}. This makes summary evaluation metrics that balance between these trade-offs popular. Commonly used objectives include accuracy and $F_1$-Score:
{\small
$$
\text{accuracy} = \frac{|TP|+|TN|}{|TP|+|TN|+|FP|+|FN|}
$$
$$
F_1\text{-Score}= \frac{2}{\text{precision}^{-1} + \text{recall}^{-1}}
$$
}
Accuracy is the rate of correct predictions to all predictions. $F_1$-Score is a specific instance of $F_\beta\text{-Score}=(1+\beta^2)\cdot(\text{precision} \cdot \text{recall})/(\beta^2 \cdot \text{precision}^{-1} + \text{recall}^{-1})$, which is the weighted harmonic mean of precision and recall. The parameter $\beta$ is chosen such that recall is considered $\beta$ times more important than precision.

When a binary classification neural network is trained in parts, independently of $\tau$ and the chosen metric, the network is unaware of the end goal and therefore is unlikely to optimize for the desired evaluation metric. This limitation often becomes evident when the goal is to balance between various confusion matrix set values, e.g.,  $F_\beta\text{-Score}$.

Prior works have approximated the Heaviside function. For instance,  \cite{kyurkchiev2015sigmoid} proposed using a sigmoid function, as shown in Figure \ref{fig:heaviside1} (right). However, this prior work focused on minimizing the error between the approximation and $H$, not on a Heaviside approximation useful for back-propagation. Such sigmoid-based approximations do not adhere to the properties necessary for gradient descent. For example, the Heaviside approximation proposed in \cite{kyurkchiev2015sigmoid}, $s_0(k;p) = (1+e^{-kp})^{-1}$ can be reparameterized to account for $\tau$: $s_{\tau}(k;p) = (1+e^{-k(p-\tau)})^{-1}$. While this approximation approaches $H$ at the limit when $\tau=0.5$, with an appropriate value of $k$, it does not necessarily approach $H$ at the limit when $\tau$ nears the limits, for a given value of $k$. Moreover, with increasing values of $k$, while the error between $s_{\tau}$ and $H$ decreases, the range over which the derivative $s'_{\tau}=0$ increases. A non-zero gradient is important during training so the network can learn from similar samples belonging to different classes. Figure \ref{fig:heaviside1} (right) shows the \cite{kyurkchiev2015sigmoid} approximation. Our proposed approximation ensures, in addition to the above mentioned properties of  $H(p=\tau)=0.5$ and reflection over $y=0.5$, that the  gradient is non-zero and $\mathcal{H}$ approaches $H$ at the limit for any threshold $\tau$.

\section{Heaviside Piecewise Linear Approximation}
\label{sec:method}

To approximate the Heaviside function in a simple manner while adhering to the important properties noted in the prior section, we propose a five-point linearly interpolated function as shown in Figure \ref{fig:heaviside2} (left). Our approximation $\mathcal{H}$ is defined over $[0, 1]$ and is parameterized by a given threshold $\tau$ and a slope parameter $\delta$, which define three linear segments with slopes: $m_1$, $m_2$, and $m_3$. Let $\tau_m = \text{min} \{\tau, 1-\tau\}$ in order to ensure a gradient that adheres to properties desirable for optimization via back propagation. Then, the slope of each line segment is defined as:
\begin{center}
    \noindent\begin{minipage}{.33\linewidth}
    \centering
        $$
        m_1 = \frac{\delta}{\tau-\frac{\tau_m}{2}}
        $$
    \end{minipage}%
    \noindent\begin{minipage}{.33\linewidth}
    \centering
        $$
        m_2 = \frac{1-2\delta}{\tau_m}
        $$
    \end{minipage}%
    \noindent\begin{minipage}{.33\linewidth}
    \centering
        $$
        m_3 = \frac{\delta}{1-\tau-\frac{\tau_m}{2}}
        $$
    \end{minipage}
\end{center}
The Heaviside approximation is thus given by:
{\small
\begin{equation}
\mathcal{H} =
  \begin{cases}
    p \cdot m_1 & \text{if $p < \tau - \frac{\tau_m}{2}$} \\
    p \cdot m_3 +(1-\delta-m_3(\tau+\frac{\tau_m}{2})) & \text{if $p > \tau + \frac{\tau_m}{2}$} \\
    p \cdot m_2 + (0.5 - m_2\tau) & \text{otherwise}
  \end{cases}
\end{equation}
}
The above equation can be trivially derived by computing the slope and intercept for each segment $m_1, m_2, m_3$. Slope-intercept form is shown for clarity.

Including the threshold $\tau$ in our parameterization of the approximation is crucial in that many evaluation metrics compute confusion matrix set cardinality based on a given threshold value or across a range of thresholds to effectively evaluate classifier performance. For example, $F_1$-Score is usually computed at a specific threshold whereas area under the receiver operating characteristic curve (AUROC) is computed over a range of $\tau$ values. AUROC is a commonly used ranking performance measure defined in terms of the true positive rate (TPR) and false positive rate (FPR) as $\text{AUROC}=\int_{\tau=0}^1 \text{TPR}(\text{FPR}^{-1}(\tau)) d\tau$. Our approximation allows for the computation of either type of metric.

We define the slope parameter $\delta$ of $\mathcal{H}$ on the interval $(0,0.5)$. A larger value of $\delta$ providers a smoother derivative but further approximation from the Heaviside function. In practice, we found that this parameter was not crucial in our experiments and that a value of $\delta = 0.1$ works well. In the next sections we abuse notation for simplicity and denote our approximation as $\mathcal{H}(p, \tau)$, although it also depends on $\delta$.

The proposed approximation $\mathcal{H}$ adheres to the important properties of the Heaviside function including reflection over $y=0.5$ and $\mathcal{H}(p=\tau)=0.5$ as well as the properties necessary for effective back-propagation. At the limits, our approximation $\mathcal{H}$ approaches the Heaviside function $H$,
$$
\lim_{p \to \{0, 1\}} \mathcal{H}(p,\tau) = H(p,\tau) \hspace{1em} \forall \, \tau
$$
Our formulation also ensures that the derivative of $\mathcal{H}$ is non-zero: $\mathcal{H}'(p, \tau) \neq 0 \ ,\ \forall \, \tau$.

\subsection{Soft Set Membership with the Approximation}

We define soft set membership functions via the Heaviside approximation. Intuitively, the idea is to output a value proportional to the input prediction value and closer to $1$ when the prediction is a set member, and a value closer to $0$ in all other cases. In the negative cases where the set membership function yields the complement, we use $1-\mathcal{H}(p, \tau)$. Formally, we define the soft set membership functions, where $s$ is for ``soft," $tp_s$, $fp_s$, $fn_s$, and $tn_s$ based on the truth tables for each confusion matrix set:
\begin{align*}
tp_s(p, y, \tau) &=
  \begin{cases}
    \mathcal{H}(p,\tau) & y = 1 \text{ or } p < \tau \\
    1-\mathcal{H}(p,\tau) & \text{else}
  \end{cases}\\
fp_s(p, y, \tau) &=
  \begin{cases}
    \mathcal{H}(p,\tau) & y = 0 \text{ or } p < \tau \\
    1-\mathcal{H}(p,\tau) & \text{else}
  \end{cases}\\
fn_s(p, y, \tau) &=
  \begin{cases}
    1-\mathcal{H}(p,\tau) & y = 1 \text{ or } p \geq \tau \\
    \mathcal{H}(p,\tau) & \text{else}
  \end{cases}\\
tn_s(p, y, \tau) &=
  \begin{cases}
    1-\mathcal{H}(p,\tau) & y = 0 \text{ or } p \geq \tau \\
    \mathcal{H}(p,\tau) & \text{else}
  \end{cases}
\end{align*}
The properties of $\mathcal{H}$ ensure soft set membership weight is equally close to $1$ in the positive case as it is close to $0$ in the negative case. Positive or negative predictions should not disproportionately skew soft set membership magnitude.

Metrics such as accuracy and $F_1$-Score that combine cardinality values from the confusion matrix sets can then be approximated using the relevant soft set membership functions. For instance, we approximate precision using  $|TP_s| = \sum_{i=1}^n tp_s(p_i, y_i, \tau)$ and $|FP_s| = \sum_{i=1}^n fp_s(p_i, y_i, \tau)$. In practice, we sum these values at training time over members of mini-batches while optimizing via gradient descent. As gradient descent and its variants expect a small but representative sample of the broader data \cite{bottou2012stochastic}, the proposed method also expects a representative sample.

\subsection{Metrics and Losses}

Any metric composed of combinations of confusion matrix set values ($TP$, $FP$, $FN$, and $TN$) can be directly optimized using the proposed Heaviside approximation. In our experiments, we train and evaluate on accuracy, $F_1$-Score, AUROC, and $F_\beta$-Score  \cite{Rijsbergen}. Accuracy is a reasonable measure for balanced datasets. In the case of imbalanced data, accuracy tends to over-estimate classifier performance, easily maximizing performance by learning to predict only dominant-class samples achieving accuracy proportional to class imbalance \cite{he2009learning}. Popular alternative metrics for accuracy when dealing with imbalanced data include $F_\beta$-score and AUROC. Our experiments demonstrate that our proposed Heaviside approximation can be used to train a classifier with either objective using the soft membership functions.

\subsection{Computational Efficiency}

At training time, a loss that utilizes our proposed Heaviside approximation has a runtime linear with regard to the number of samples. In practice however, optimizing for a metric with $\mathcal{H}(p, \tau)$ over all samples in each batch could lead to increased run time due to the number of constant-time operations required to compute a metric. To mitigate this and further minimize run time, we observe that $\mathcal{H}(p, \tau)$ can be replaced with an reasonably-sized $O(1)$ lookup table by truncating $p$ to several decimal places and precomputing $\mathcal{H}$ for values of $p$ and $\tau$ over the range $[0,1]$. For example, if an interval between values of $\tau$ is set to $0.1$, as was done in our experiments, and $p$ is truncated at two decimal places, the lookup table has only $1,000$ elements. Using 8-bit storage, this table -- which only needs to be computed once -- consumes $1$kB of memory.

\section{Experiments}

This section presents four experiments to (1) understand the effect of different batch sizes in the proposed approach, (2) compare our proposed Heaviside approximation with a sigmoid approximation in the context of binary classification, (3) evaluate the performance of our approach against optimizing for other relevant classification metrics, and (4) evaluate optimizing directly for $F_\beta$-Score with our Heaviside approximation in order to balance precision and recall.

\subsection{Architecture and Training}

Unless otherwise noted, all experiments use the same network architecture and training scheme. The binary classifier is a feedforward neural network classifier consisting of three fully connected layers of 32 units, 16 units, and 1 unit. The first two layers have a Rectified Linear Unit (ReLU) activation \cite{nair2010rectified} and are followed by dropout \cite{hinton2012improving}. The final layer is a sigmoid-activated single-unit output. We use the ADAM optimizer \cite{kingma2014adam} with a learning rate of $0.001$ and a batch size of 2048. The network is trained using early stopping, monitoring validation loss for improvement over a sliding window of 100 epochs.
The maximum number of training iterations is set to 5000, a number high enough such that early stopping takes effect before finishing all epochs. The models for all losses were trained with TensorFlow, except for \cite{fathony2020ap} and AUROC via the proposed Heaviside approximation which were implemented in PyTorch. Note that we provide open source implementations of our proposed approach for both deep learning libraries.

Training systems used either an NVIDIA Titan X or RTX 2080ti GPU, Intel i7 3.7GHz processor, and 32GB of RAM.

\subsubsection{Datasets}

Experiments were conducted on two synthetic datasets (Section \ref{sig-vs-lin-approx}) and four publicly available datasets with binary classification labels. The latter datasets were chosen for their varying levels of class imbalance in a variety of domains, as described in Section \ref{direct-optimization-experiment}. All of the datasets were minimally pre-processed, the details of which are described in the supplementary material. Features were centered and scaled to unit variance, and the data was split into separate train (64\%), test (20\%) and validation (16\%) sets.
\subsection{Batch Size Effects}

Our first experiment explored the effects of batch size on approximating non-decomposable metrics, such as $F_1$-Score, for classification of microcalcifications in the Mammography dataset \cite{woods1993comparative}. The dataset is imbalanced, having only 2.32\% positive-class examples as detailed in the supplementary.

Throughout training, $F_1$-Score was computed against both the current batch of the train split as well as the entire dataset using the network undergoing training. The difference between the batch $F_1$-Score and the $F_1$-Score of the entire dataset was then aggregated. Finally, we computed the absolute mean and standard deviation for the difference over all batches used during training.
\begin{table}[tb]
    \centering
    {\small
        \begin{tabular}{c|c|c|c}
            $B=128$ & $B=1024$ & $B=2048$ & $B=4096$\\ \hline
            $0.041\pm0.23$ & $0.005\pm0.06$ & $0.009\pm0.05$ & $0.006\pm0.03$
        \end{tabular}
    }
    \caption{Difference between batch $F_1$-Score and the $F_1$-Score of the entire Mammography dataset. Each column corresponds to a batch size $B$. See the text for more details.}
    \label{tbl:batch_results}
\end{table}
We performed this test for batch sizes of $\{128, 1024, 2048, 4096\}$ samples, which corresponds to $\{1\%, 9\%, 18\%, 37\%\}$ of the total samples of the dataset. As can be seen in Table \ref{tbl:batch_results}, the difference between batch $F_1$-Score and total $F_1$-Score decreases as the batch size increases. This trend is observed both in terms of absolute mean difference and standard deviation.

The relative performance of classifiers across batch sizes is shown in Figure \ref{fig:batch-size}. Directly optimizing Accuracy results in zero precision, recall, and $F_1$-Score because all predictions are negative. Training on $F_1$-Score, AUROC, and BCE produce non-zero $F_1$-Scores that are similar across batch sizes with the most noticeable difference being BCE, which is perhaps due to the regularizing effect when using small batches \cite{wilson2003inefficiency}

For the next experiments described in the paper, we chose a batch size of 1024, which worked well with our method and preliminary experiments with BCE.

\begin{figure}[tb!p]
    \centering
    \begin{minipage}{0.35\textwidth}
        \centering
        \includegraphics[width=1\linewidth]{plots/batch/batch_size_comparison.pdf}
    \end{minipage}
    \caption{Comparison of batch sizes across losses (rows) and metrics (columns). The Accuracy loss is computed via the proposed Heaviside approximation $\mathcal{H}$, $F_1$-Score via $\mathcal{H}$, AUROC via Wilcoxon-Mann-Whitney statistic \cite{yan2003optimizing} and BCE is binary cross-entropy loss. The top charts that seem empty indicate zero performance with the Accuracy loss. Bar height corresponds to the score in each case.}
    \label{fig:batch-size}
\end{figure}
\subsection{Sigmoid vs. Piecewise Linear Approximation}
\label{sig-vs-lin-approx}

Prior work proposed a sigmoid approximation for the Heaviside function \cite{kyurkchiev2015sigmoid}. Thus, we compare our proposed piecewise linear Heaviside approximation with a sigmoid parameterization fit to the piecewise function using non-linear least squares (Fig. \ref{fig:heaviside2} right). We fit a sigmoid to our piecewise linear approximation instead of using the original sigmoid approximation from \cite{kyurkchiev2015sigmoid} because the latter approximation poses challenges for gradient descent, as discussed in Sec. \ref{sec:preliminaries}.

Two synthetic datasets were used for this experiment to control for class imbalance. Both synthetic datasets were generated by creating two isotropic gaussian blobs with $\sigma = 10$ in 3D space using a fixed random seed value and resulting in $10000$ samples each. We then removed a randomly-sampled proportion of the positive data points. This led to datasets with a positive to negative ratio of 1:1 (50\%) and 1:3 (33\%). The total number of samples in the final datasets were 7500 and 6250, respectively.

We trained and evaluated models on accuracy and $F_{1}$-Score by directly optimizing for the respective metric over $10$ trials using the sigmoid or piecewise linear Heaviside approximation. The results in Table \ref{tbl:sigmoid-vs-linear} show similarity between approximations on the dataset with balanced classes and on the dataset with $33\%$ positive class balance. We argue in favor of the piecewise linear approximation since it has similar performance to the Sigmoid while avoiding the added complexity of curve fitting.
\begin{table}[b!p]
    \begin{center}
        \begin{tabular}{c|cc}
            \multicolumn{3}{c}{Synthetic 50\%} \\ \hline
            Loss & Accuracy Score & $F_1$-Score \\ \hline

            Sigmoid & $0.839 \pm 0.008$ & $0.845 \pm 0.008$ \\ \hline

            Linear* & $0.838 \pm 0.008$ & $0.846 \pm 0.008$ \\
        \end{tabular}
        \end{center}
    \begin{center}
        \begin{tabular}{c|cc}
            \multicolumn{3}{c}{Synthetic 33\%}  \\ \hline
            Loss & Accuracy Score & $F_1$-Score \\ \hline
            Sigmoid & $0.845 \pm 0.028$ & $0.782 \pm 0.020$ \\ \hline
            Linear* & $0.849 \pm 0.010$ & $0.779 \pm 0.010$ \\
        \end{tabular}
        \end{center}
        \caption{Networks trained to directly optimize for Accuracy score and $F_1$-Scores on synthetic datasets with even class balance (top) and 33\% positive samples (bottom). \textit{* indicates use of the proposed method.}}
        \label{tbl:sigmoid-vs-linear}
\end{table}
\subsection{Direct Optimization of Metrics Based on the Confusion Matrix Sets}
\label{direct-optimization-experiment}

Using four different datasets, we evaluate the performance of neural networks that are trained to optimize directly for accuracy, $F_{1}$-Score, and AUROC using our proposed Heaviside approximation. We compare these models to optimizing $F_1$-score via the adversarial approach by \cite{fathony2020ap}, optimizing AUROC via the Wilcoxon-Mann-Whitney (WMW) Statistic per \cite{yan2003optimizing}, and the classical binary cross entropy (BCE) loss. We evaluate the performance of these networks on accuracy score, $F_{1}$-Score, and AUROC over 10 repeated trials, which allow to measure the variability of the results.

The same batch size and learning rate are used for all methods except for \cite{fathony2020ap}, for which we use the hyperparameters suggested by the first author through personal communication. Direct comparison with \cite{fathony2020ap} is challenging as the maximum batch size for this method is 25 and in practice, set to only 20 samples due to the increased runtime associated with larger batch sizes.

Tables \ref{tbl:results-cocktailparty}, \ref{tbl:results-adult}, \ref{tbl:results-mammography}, and \ref{tbl:results-kaggle} show the results from the different datasets, which are further discussed in the next Sections. We report the mean of results over threshold values $\tau = \{0.1, 0.2, \cdots, 0.9\}$ across all trials.

\subsubsection{CocktailParty Dataset Results}

Results on dyadic, social interaction classification using the CocktailParty dataset are shown in Table \ref{tbl:results-cocktailparty}. This dataset has 4800 samples and a 30.29\% positive class balance making it both the smallest and the most class-balanced dataset among those considered in this experiment. Prior work has evaluated classification models for this dataset via pairwise accuracy \cite{hedayati2019}; we also consider other metrics.

Experimental results show training on accuracy via the proposed Heaviside approximation leads to the best accuracy and $F_1$-Score. When evaluating on the AUROC metric, BCE loss shows the best performance, possibly due to the small dataset size resulting in thresholds for which the true positive rate or false positive rate is ill-defined. In the case of $F_1$-Score evaluation, training on accuracy may outperform training on $F_1$-Score because of the reasonable class-balance in the dataset and the fact that accuracy incorporates the $TN$ set while $F_1$-Score does not take it into account.
\begin{table}[tb!p]
    \begin{center}
        {\small
            \begin{tabular}{c|ccc}
                \multicolumn{4}{c}{\textbf{CocktailParty} Results ($\mu\pm\sigma$)} \\
                \hline
                Loss & Accuracy & $F_1$-Score & AUROC \\ \hline
                Accuracy* & $0.846 \pm 0.01$ & $0.897 \pm 0.01$ & $0.776 \pm 0.02$ \\ \hline
                $F_1$-Score* & $0.782 \pm 0.07$ & $0.865 \pm 0.04$ & $0.650 \pm 0.12$ \\ \hline
                $F_1$-Score\dag & $0.764 \pm 0.01$ & $0.348 \pm 0.04$ & $0.615 \pm 0.01$ \\ \hline
                AUROC* & $0.637 \pm 0.04$ & $0.589 \pm 0.11$ & $0.707 \pm 0.07$ \\ \hline
                AUROC\ddag & $0.522 \pm 0.22$ & $0.424 \pm 0.41$ & $0.555 \pm 0.11$ \\ \hline
                BCE & $0.830 \pm 0.05$ & $0.876 \pm 0.05$ & $0.793 \pm 0.07$ \\ \hline
             \end{tabular}
         }
    \end{center}
    \caption{Losses (rows): Accuracy*, $F_1$-Score*, and AUROC* optimize the proposed linear Heaviside approximation $\mathcal{H}$. $F_1$-Score\dag\ uses the adversarial approach of \cite{fathony2020ap} and AUROC\ddag\ the WMW statistic \cite{yan2003optimizing}. BCE is the typical log loss. \textit{* indicates use of the proposed method.}}
    \label{tbl:results-cocktailparty}
\end{table}
\subsubsection{Adult Dataset Results}

Results on salary classification using the Adult dataset are shown in Table \ref{tbl:results-adult}. This dataset has 23.93\% positive data balance. Prior work has evaluated this dataset on accuracy score \cite{kohavi1996scaling}. In our experiments, when evaluating on accuracy, we obtain best performance by training on accuracy as a loss via the proposed Heaviside approximation. We believe accuracy is a reasonable evaluation metric for this dataset as the class balance is not extreme; however, accuracy may yet overestimate classifier performance \cite{he2009learning}. Evaluating on metrics that better deal with class imbalance, $F_1$-Score and AUROC, we find that optimizing on $F_1$-Score as a loss via $\mathcal{H}$ results in the best performance.
\begin{table}[tb!p]
    \begin{center}
        {\small
            \begin{tabular}{c|ccc}
                \multicolumn{4}{c}{\textbf{Adult} Results ($\mu\pm\sigma$)} \\
                \hline
                Loss & Accuracy & $F_1$-Score & AUROC \\ \hline
                Accuracy* & $0.813 \pm 0.01$ & $0.380 \pm 0.05$ & $0.617 \pm 0.02$ \\ \hline
                $F_1$-Score* & $0.772 \pm 0.04$ & $0.625 \pm 0.02$ & $0.781 \pm 0.02$ \\ \hline
                $F_1$-Score\dag & $0.788 \pm 0.00$ & $0.199 \pm 0.02$ & $0.562 \pm 0.01$ \\ \hline
                AUROC* & $0.470 \pm 0.05$ & $0.424 \pm 0.04$ & $0.586 \pm 0.02$ \\ \hline
                AUROC\ddag & $0.414 \pm 0.25$ & $0.353 \pm 0.15$ & $0.539 \pm 0.08$ \\ \hline
                BCE & $0.748 \pm 0.11$ & $0.315 \pm 0.14$ & $0.599 \pm 0.06$ \\ \hline
            \end{tabular}
        }
    \end{center}
    \caption{Losses (rows): Accuracy*, $F_1$-Score*, and AUROC* optimize the proposed linear Heaviside approximation $\mathcal{H}$. $F_1$-Score\dag\ uses the adversarial approach of \cite{fathony2020ap} and AUROC\ddag\ the WMW statistic \cite{yan2003optimizing}. BCE is the typical log loss. \textit{* indicates use of the proposed method.}}
    \label{tbl:results-adult}
\end{table}
\subsubsection{Mammography Dataset Results}
Results on classification of microcalcifications on the Mammography dataset are presented in Table \ref{tbl:results-mammography}. As the data balance is more skewed in this case (2.32\% positive), we expect accuracy to be less-representative of performance due to the classifier's propensity to learn to predict only dominant-class examples and achieve an accuracy proportional to the class imbalance \cite{he2009learning}. The results confirm this phenomenon with accuracy, highlighting the value of alternative evaluation measures on class-imbalanced data.

Prior work with the Mammography dataset evaluates results on both precision and false positive rate independently \cite{woods1993comparative}. In lieu of a summary metric to directly optimize, we train and evaluate on $F_1$-Score. We find that directly optimizing the metric via the proposed Heaviside approximation results in close to or better than baseline performance. AUROC provides another summary metric, but considers only positive samples.
\begin{table}[tb!p]
    \begin{center}
        {\small

            \begin{tabular}{c|ccc}
                \multicolumn{4}{c}{\textbf{Mammography} Results ($\mu\pm\sigma$)} \\
                \hline
                Loss & Accuracy & $F_1$-Score & AUROC \\ \hline
                Accuracy* & $0.976 \pm 0.00$ & $0.000 \pm 0.00$ & $0.500 \pm 0.00$ \\ \hline
                $F_1$-Score* & $0.986 \pm 0.00$ & $0.699 \pm 0.04$ & $0.838 \pm 0.03$ \\ \hline
                $F_1$-Score\dag & $0.985 \pm 0.00$ & $0.554 \pm 0.15$ & $0.715 \pm 0.06$ \\ \hline
                AUROC* & $0.521 \pm 0.05$ & $0.105 \pm 0.06$ & $0.614 \pm 0.09$ \\ \hline
                AUROC\ddag & $0.636 \pm 0.42$ & $0.143 \pm 0.19$ & $0.606 \pm 0.15$ \\ \hline
                BCE & $0.984 \pm 0.00$ & $0.558 \pm 0.18$ & $0.743 \pm 0.11$ \\ \hline
            \end{tabular}
        }
    \end{center}
    \caption{Losses (rows): Accuracy*, $F_1$-Score*, and AUROC* optimize the proposed linear Heaviside approximation $\mathcal{H}$. $F_1$-Score\dag\ uses the adversarial approach of \cite{fathony2020ap} and AUROC\ddag\ the WMW statistic \cite{yan2003optimizing}. BCE is the typical log loss. \textit{* indicates use of the proposed method.}}
    \label{tbl:results-mammography}
\end{table}
\subsubsection{Kaggle Dataset Results}

The Kaggle Credit Card Fraud Detection dataset has the most extreme class imbalance of our datasets with only 0.17\% positive samples. Again, we see the effects of class imbalance on the accuracy metric and loss: Table \ref{tbl:results-kaggle} shows near perfect accuracy scores for all losses except AUROC. In the case of large class imbalance, AUROC is a challenging metric to directly optimize, as there are thresholds for which there are no positive samples and AUROC is ill-defined. One possible approach to addressing this is to apply ideas from sample mining, a common technique in deep metric learning \cite{Schroff2015facenet,wang2019deepmetriclearning}, selecting samples for a representative batch.

In line with \cite{davis2006relationship}, Kaggle acknowledges the class imbalance issue and recommends evaluation via AUPRC.\footnote{https://www.kaggle.com/mlg-ulb/creditcardfraud} However, the AUPRC summary metric is problematic due to incoherent scale assumptions \cite{flach2015precision}. We find that optimizing on $F_1$-Score directly via the heaviside approximation is effective in comparison with the baselines.
\begin{table}[tb!p]
    \begin{center}
        {\small
            \begin{tabular}{c|ccc}
                \multicolumn{4}{c}{\textbf{Kaggle} Results ($\mu\pm\sigma$)} \\
                \hline
                Loss & Accuracy & $F_1$-Score & AUROC \\ \hline
                Accuracy* & $0.998 \pm 0.00$ & $0.000 \pm 0.00$ & $0.500 \pm 0.00$ \\ \hline
                $F_1$-Score* & $0.999 \pm 0.00$ & $0.820 \pm 0.04$ & $0.901 \pm 0.02$ \\ \hline
                $F_1$-Score\dag & $0.999 \pm 0.00$ & $0.739 \pm 0.15$ & $0.833 \pm 0.07$ \\ \hline
                AUROC* & $0.479 \pm 0.08$ & $0.214 \pm 0.02$ & $0.694 \pm 0.06$ \\ \hline
                AUROC\ddag & $0.882 \pm 0.28$ & $0.470 \pm 0.35$ & $0.814 \pm 0.15$ \\ \hline
                BCE & $0.999 \pm 0.00$ & $0.553 \pm 0.26$ & $0.736 \pm 0.14$ \\ \hline
            \end{tabular}
        }
    \end{center}
    \caption{Losses (rows): Accuracy*, $F_1$-Score*, and AUROC* optimize the proposed linear Heaviside approximation $\mathcal{H}$. $F_1$-Score\dag\ uses the adversarial approach of \cite{fathony2020ap} and AUROC\ddag\ the WMW statistic \cite{yan2003optimizing}. BCE is the typical log loss. \textit{* indicates use of the proposed method.}}
    \label{tbl:results-kaggle}
\end{table}
\subsubsection{Discussion} The choice of evaluation metric for a given dataset depends on the dataset's class balance \cite{he2009learning, davis2006relationship} and the effects of this choice is amplified when optimizing the metric as a loss. Our results show the best classifiers tend to follow from networks trained on accuracy or $F_1$-score, depending class imbalance. Comparing the optimization of a given metric as a loss, our method tends to perform comparably or better than other methods for optimizing the metric as well as the BCE baseline. The AUROC results are unusual in that models trained with BCE tend to outperform both the proposed piecewise linear approximation method and the WMW statistic method for optimizing AUROC. AUROC is computed over a range of thresholds, which is problematic when positive samples for a given threshold do not exist.

\subsection{Balancing between Precision and Recall}

Beyond the ability to directly optimize confusion matrix based metrics, our method also allows training-time optimization that balances between precision and recall using $F_\beta$-Score. This end-to-end approach to training is particularly useful in real-world scenarios where there is a high cost associated with missed detections.

Table \ref{tbl:fb-results} shows the effectiveness of optimizing for $F_{\beta}$-Score on the Mammography dataset.  We find that directly optimizing $F_{\beta}$-Score is an effective way of maintaining maximum classifier performance measured via $F_1$-Score, while balancing between precision and recall at a ratio appropriate for the task at hand. Increasing values of $\beta$ correspond to an increased preference toward recall with small loss of total performance measured by $F_1$-Score.

\begin{table}[tb]
    \begin{center}

            \begin{tabular}{c|ccc}
                \multicolumn{4}{c}{\textbf{Mammography} $F_\beta$-Score Results ($\mu\pm\sigma$)} \\ \hline
                Loss & $F_1$-Score & Precision & Recall \\ \hline
                $F_1$-Score* & $0.679 \pm0.0$ & $0.727 \pm0.1$ & $0.645 \pm0.1$ \\ \hline
                $F_2$-Score* & $0.678 \pm0.0$ & $0.627 \pm0.1$ & $0.746 \pm0.1$ \\ \hline
                $F_3$-Score* & $0.579 \pm0.1$ & $0.475 \pm0.1$ & $0.752 \pm0.1$ \\ \hline
            \end{tabular}
        }
    \end{center}
    \caption{Balancing between precision and recall at training time by optimizing $F_\beta$-Score via the piecewise linear Heaviside approximation. Losses (rows): $F_\beta$-Score over $n=10$ trials, where $\beta = \{1,2,3\}$. \textit{* indicates use of the proposed method.}}
    \label{tbl:fb-results}
\end{table}

\section{Related Work}

Our work is inspired by a rich line of research focused on the direct optimization of evaluation metrics for binary classification. \cite{narasimhan2014statistical} explore plug-in algorithms that learn a classifier by applying an empirically determined threshold to suitable estimates of class probabilities. Their experiments demonstrate the applicability of their approach to optimizing $F_1$-Score with linear models. For metrics based on linear combinations of confusion matrix set cardinalities, \cite{koyejo2014consistent} identify an optimal plug-in classifier with a metric-dependent threshold. Also for linear models, \cite{kar2015surrogate} explore the precision@k metric in the context of ranking.

In the area of online learning, characterized by the sequential availability of data, works such as \cite{kar2014online, busa2015online} optimize specific metrics such as precision@k and $F_1$-Score. Our work does not address online learning, but instead batch methods. Our Heaviside approximation is applicable to any confusion matrix based metrics without reformulation.

Other prior work has focused on optimized a specific metric, like AUROC \cite{yan2003optimizing, herschtal2004optimising}, F-score \cite{dembczynski2011exact, zhao2013beyond}, and AUPRC in the context of ranking \cite{eban2016scalable}. We considered in particular the work by \cite{yan2003optimizing} as a baseline in our experiments. Our results showed that using our proposed Heaviside approximation for computing AUROC results in competitive results to using the WMW statistic for the AUROC loss.

Recently, adversarial approaches have emerged as another active area of research. \cite{wang2015adversarial} uses a structured support vector machine and reports on precision@k as well as $F_1$-Score. \cite{fathony2020ap} improves performance via a marginalization technique to ensure polynomial time convergence and evaluates on accuracy and $F_1$-Score, among other metrics, while reporting performance relative to BCE. One downside of the latter approach is small batch size, which was on the order of 25 samples in their evaluation and in our own experiments.

A parallel line of work has focused on dealing with class imbalance by dropping data in the over-represented class or augmenting the under-represented data by synthesizing new samples. Within this approach of re-sampling to deal with class imbalance, \cite{drummond2003c4} finds under-sampling preferential to over-sampling. However, in both cases it is vital to maintain representative distributions for both the over and under-represented classes \cite{he2009learning}. Therefore, strategic approaches toward handling imbalance via data preprocessing have been proposed. Synthetic Minority Over-sampling Technique (SMOTE) \cite{chawla2002smote} for example, generates synthetic data based on feature space similarity. However, recent work finds re-sampling based approaches inadequate and suggests a larger number of samples is necessary to improve results \cite{juba2019precision}. As an alternative to re-sampling, our work focuses on choosing an appropriate evaluation metric for the given level of class imbalance and directly optimizing this metric as a loss.

Our proposed method is unique in its flexibility. It allows optimization of any metric based on confusion matrix sets, even when computation of a range of thresholds is required, e.g., as in the case of AUROC. This means that if one chooses to optimize AUROC, but later finds that in the given application it would be more appropriate to specify the trade-off between precision and recall, one can easily switch to optimizing $F_\beta$-Score.

\section{Conclusion}

We proposed a Heaviside function approximation that is useful for direct optimization of metrics based on confusion matrix set cardinalities. Our experiments showed the efficacy and flexibility of this approach in training end-to-end binary classifiers compared to classifiers that are trained with the traditional binary cross-entropy loss, unaware of their final evaluation metric. Further, our experiments demonstrated that our proposed method is competitive with other approaches for training binary classifiers that optimize for accuracy, $F_1$-Score, and AUROC. In the future, we are interested in studying the use of our Heaviside approximation on multi-class classification problems.

Our work was motivated by the wide ranging impact that more flexible and robust binary classifiers would have across application domains, especially in the case of class imbalance. In our experiments, we showed applications related to social group dynamics, sociology, economics, and medicine. Improving the tools used in these areas of research has the potential to positively impact human quality of life. It is worth noting that end-to-end classifiers that directly optimize the evaluation metric as a loss remove the disconnect between training and evaluation, but are still susceptible to data bias. In real application scenarios, these tools must be used with care, keeping in mind the goal of building safe and responsible artificial intelligence systems.


\section{Introduction}
The moving finite element method(MFEM) was first developed  in \cite{miller1981movingA,miller1981movingB} about forty years ago.
It  is a typical  $r$-type adaptive method\cite{huang1994moving,cao2002moving,li2001moving,tang2005moving,budd2009adaptivity},
where the mesh vertexes are relocated without changing the mesh topology.
In the MFEM, the mesh relocation
is done by solving a dynamic equation for the vertexes
 coupled with the original partial differential equations.
No interpolation is needed in the method since the mesh is continuous with time.
The MFEM has arisen considerable interest and has been further developed in several directions(c.f \cite{Baines,dupont1982mesh,carlson1998designA,wathen1985structure,baines2005moving,bank2019diagonally} among many others).

However, like  all other $r$-adaptive methods, the theoretical analysis for the MFEM is far from being complete.
The first error analysis for  MFEM was done by Dupont \cite{dupont1982mesh}, where he proved the optimal convergence
of the method
for a linear convection diffusion equation when the solution is smooth. This is not enough since we are more interested in
 non-smooth solutions for adaptive methods.
Later on, Jimack proved the locally optimal approximation for the stationary solution of a linear parabolic equation\cite{jimack1992steady,jimack1996best,jimack1996optimal}. This is remarkable since
the solution is allowed to have lower regularity  in Jimack's proof.
In this study, we aim to do analysis for a nonlinear gradient flow system by using the Onsager variational principle as an approximation tool.

The Onsager variational principle is a fundamental principle for irreversible thermodynamic processes
in statistical physics\cite{Onsager1931,Onsager1931a,DoiSoft}. It has be
used to model many dissipative physical systems\cite{Doi2011,DoiSoft},
 such as  the Stokes equation in hydrodynamics, the Ericksen-Leslie equation in liquid
crystal, and the GNBC boundary condition for moving contact lines\cite{QianTiezheng2006}, etc.
Recently, the Onsager principle has been used as an approximation tool for many problems in two-phase flows
and in material science\cite{Doi2015,XuXianmin2016,DiYana2016,ManXingkun2017,ZhouJiajia2018,Jiang19b}.
In particular, it has  been used to derive an efficient numerical method for wetting dynamics\cite{LuXu2020}.

In this work, we first give a new derivation of the MFEM  for a gradient flow system by using the Onsager principle as an approximation tool.
 The key idea is to approximate the system in a nonlinear approximation space of free-knot piecewise polynomials\cite{devore1998nonlinear}.
Both the mesh vertexes and the nodal values of the finite element function  are regarded as  unknowns.
We  derive a system of ordinary differential equations(ODEs) for them. The ODE system  coincides
with the discrete equation of the MFEM, which has been derived in a totally different way in \cite{miller1981movingA}.
Here we do not need to compute the multiply of a Dirac measure and a discontinuous function,
so that the ``mollification" technique or any other formally calculation is not needed.
Furthermore, our derivation shows that the discretized problem has the same energy dissipation structure of
the continuous system. This makes us to prove the energy decay property of the discrete problem easily.

 Based on the variational formula, we do  error analyse for
the MFEM for a stationary solution of the gradient flow system.
The analysis can be regarded as a generalization of the results in \cite{jimack1996best} to nonlinear equations.
We show that the MFEM gives locally best approximations to the energy.
When a global minimizer is detected, an optimal error estimate is proved using the nonlinear
approximation theory.
Numerical examples show that the optimal convergence
can be obtained for a linear diffusion equation and for  stationary solutions of a nonlinear Allan-Cahn equation.
In this paper, we mainly consider the one dimensional problem. All the results can be generalized
to higher dimensional cases directly.

The structure of the paper is as follows. In section 2, we introduce the Onsager variational principle and show
that it can be used to derive the partial differential equation model for a gradient flow system.
In Section 3, we derive the MFEM by using the Onsager principle as an approximation tool.
In Section 4, we do error analysis for the stationary solution of a nonlinear reaction diffusion equation.
Some numerical examples are illustrated to verify the analytical results in the last section.

\section{The Onsager variational principle for a gradient flow system}
\subsection{The Onsager principle}
Suppose a physical system is described by a time dependent function $u$.
For simplicity, we denote by $\dot{u}=\frac{\partial u}{\partial t}$ the time derivative of $u$.
For a dissipated system, the evolution of $u$ can dissipate energy.
The dissipation function is defined as half of the total energy dissipated with respect to the flux $\dot{u}$(c.f. \cite{DoiSoft}).
For a simple gradient flow system, we assume the dissipation function is
\begin{equation}\label{e:diss}
\Phi(\dot{u})=\frac{\xi}{2}\|\dot{u}\|^2,
\end{equation}
where $\xi$ is a positive friction coefficient and $\|\cdot\|$ is a $L^2$ norm.
Suppose that the free energy of the system is given by a functional $\mathcal E(u)$.
For given $\dot u$, the rate of change of the total energy is calculated by
\begin{equation}\label{e:engRate}
\dot{\mathcal E}(u;\dot{u})=\langle \frac{\delta \mathcal E(u)}{\delta u},\dot{u}\rangle.
\end{equation}
Then  a Rayleighian (functional)  is defined as
\begin{equation}\label{e:Reyl}
\mathcal{R}(u;\dot{u})=\Phi(\dot{u})+\dot{\mathcal E}(u;\dot{u}).
\end{equation}

With these definitions, the Onsager  principle can be stated as follows(\cite{DoiSoft}). For
any given $u$ at the present time, the time derivative $\dot{u}$ is obtained by minimizing the Rayleighian
among its all possible choices. In other words, the  evolution equation of $u$ is determined by
minimizing the Rayleighian $\mathcal{R}(u,\dot{u})$ with respect to $\dot{u}$, i.e.
\begin{equation}\label{e:Onsager}
\min_{\dot{u}\in V} \mathcal{R}(u,\dot{u}).
\end{equation}
Here $V$ is the admissible  space of $\dot{u}$.
Since the Rayleighian is a quadratic form with respect to $\dot{u}$,  the problem~\eqref{e:Onsager} is equivalent to  its
Euler-Lagrange equation
\begin{equation}\label{e:gf1}
\xi (\dot{u},\psi) =- \langle\frac{\delta \mathcal E(u)}{\delta u},{\psi}\rangle,\qquad\qquad \forall \psi\in V,
\end{equation}
or in a simple form,
\begin{equation}\label{e:forcebalance}
\xi\dot{u}=-\frac{\delta \mathcal E(u)}{\delta u}.
\end{equation}
 In physics, the equation indicates a balance between the friction force $-\xi\dot{u}$ and
 the general driven force $\frac{\delta \mathcal E(u)}{\delta u}$.
The equation \eqref{e:forcebalance}  can be simply rewritten as a gradient flow equation,
\begin{equation}\label{e:gf}
\frac{\partial u }{\partial t}=-\xi^{-1} \frac{\delta \mathcal \mathcal E(u)}{\delta u}.
\end{equation}

It is easy to see that the solution $u$ of Equation \eqref{e:gf1}(or \eqref{e:gf}) satisfies  the following energy decay property
\begin{equation}\label{e:engdecay}
\frac{d \mathcal E(u)}{d t}\leq 0.
\end{equation}
Actually, by setting $\psi=\frac{\partial u}{\partial t}$ in \eqref{e:gf1}, we have
\begin{align}\label{e:engdecay1}
\frac{d \mathcal E}{d t} =\langle \frac{\delta \mathcal E(u)}{\delta u},\frac{\partial u }{\partial t} \rangle=-\xi \|\frac{\partial u }{\partial t}\|_{V}^2=-2\Phi(\dot{u})
\leq 0.
\end{align}
The rate of decreasing of the total energy is equal to twice of the dissipation function.

\subsection{The model problem}
For simplicity in presentations,  we mainly consider a specific gradient flow system in this paper,
which corresponds to a (nonlinear) reaction diffusion equation.

 Denote by  $\Omega$ a domain in $R^n$.
Suppose the energy functional is given by
\begin{equation}\label{e:Eng}
\mathcal E(u,\nabla u)=\int_{\Omega} \frac{\alpha}{2}(\nabla u)^2+F(x,u) dx,
\end{equation}
where $\alpha>0$ is the diffusion coefficient and $F(x,u)$ is a function with respect to $x$ and $u$.
We assume that $u\in H^1_0(\Omega)$ and $\dot u\in L^2(\Omega)$.
The dissipation function is given by
\begin{equation}
\Phi(\dot{u})=\frac{\xi}{2}\int_{\Omega} \dot u^2 dx,
\end{equation}
with a positive friction coefficient $\xi$.  Then the Rayleighian is calculated as
\begin{align}
\mathcal{R}(u,\dot u)&=\Phi(\dot u)+\dot{\mathcal E}=\frac{\xi}{2}\int_{\Omega} \dot u^2 dx+ \int_{\Omega} \alpha \nabla u\cdot \nabla \dot u+ \partial_u F(x,u) \dot udx.
\label{e:Rayleig1}
\end{align}
We minimize $\mathcal{R}$ with respect to $\dot u$.

By integral by part, $\mathcal{R}$ is rewritten as
\begin{align}\nonumber
\mathcal{R}(u,\dot u)&=\frac{\xi}{2}\int_{\Omega} \dot u^2 dx+ \int_{\Omega} [-\alpha\Delta u+\partial_u F(x,u) ]\dot udx.
\end{align}
The corresponding Euler-Langrange equation of \eqref{e:Rayleig1}  is
\begin{equation}\label{e:modelpb}
\partial_t u -\frac{\alpha}{\xi}\Delta u+\frac{1}{\xi}  f(x,u) =0,
\end{equation}
where $f(x,u)=\partial_u F(x,u)$.
It is a quasi-linear reaction diffusion equation.
Some typical examples include the linear diffusion equation(when $F=0$) and the Allen-Cahn equation(when $F$ is
a double-well function), etc.

In general,  we assume that $f(x,u)$ is Liptchitz continuous with respect to $u$, i.e.
\begin{equation}\label{e:Lip}
| f(x,v)- f(x,w)|\leq L_0 |v-w|, \qquad \forall v, w,
\end{equation}
for some constant $L_0>0$.  Under this condition and some assumption on the domain $\Omega$,
the equation \eqref{e:modelpb} has a unique  solution (c.f. Theorem 5.1 in \cite{pao2012nonlinear}).

\section{Derivation of the moving finite element method by the Onsager Principle}

The Onsager principle can be used to derive numerical methods for the model problem~\eqref{e:modelpb}.
The  idea is as follows. We  choose a finite dimensional subspace $V_h$ of $H^1_0(\Omega)$.
For a time dependent function $u_h(t,x)\in V_h$,
we can compute  the energy functional $\mathcal E(u_h)$  and the dissipation function
$\Phi(\dot{u}_h)$. By minimizing the discrete Rayleighian
$
\mathcal{R}(u_h;\dot{u}_h):=\dot{E}(u_h;\dot{u}_h)+\Phi(\dot{u}_h)
$
with respect to $\dot{u}_h$, we  obtain a dynamic equation for $u_h(t,x)$,
which is a (semi-)discrete problem for \eqref{e:modelpb}.
For example, if we choose $V_h$ as a standard Lagrangian finite element space,  it is easy to verify that we
can obtain a finite element problem same as that derived by the standard Galerkin approach.
In the following, we  choose $V_h$ as a nonlinear approximation space composed of
free-knot piecewise polynomials. Then we derive the MFEM proposed in \cite{miller1981movingA}.

For simplicity in notations, we consider only the one-dimensional problem hereinafter.
The domain $\Omega$ is simply  a bounded interval $I=(a,b)$.  The model problem~\eqref{e:modelpb} is
reduced to
\begin{equation}\label{e:modelpb1}
\partial_t u -\frac{\alpha}{\xi}\partial_{xx} u+\frac{1}{\xi}  f(x,u) =0.
\end{equation}

\subsection{Nonlinear  approximation space}
We  recall some known results on the nonlinear approximation space of the piecewise linear functions with free-knots \cite{devore1998nonlinear}.
Let $N$ be a positive integer and let
$${X}:=\{a=:x_0< x_1<\cdots<x_N:=b\}$$
be a set   of ordered points in the interval $\bar{I}=[a,b]$. This generates a partition of $I$,
which we denote as
$\mathcal{T}({X}):=\{I_k\}_{k=1}^N$, where $I_{k}=(x_{k-1},x_k)$. Let
\begin{equation}\label{e:finiteElement}
V_h(X):=\{ v_h\in C([a,b]) | v_h \hbox{ is linear in }  I_k, \forall k=1,\cdots, N \},
\end{equation}
be the standard piecewise linear finite element space with respect to the partition $\mathcal T(X)$.
Denote by ${V}_h^N$ the space with $N$ intervals as follows
\begin{equation}\label{e:Nspace}
V_h^N:=\cup_{\#(X)=N+1} V_h(X),
\end{equation}
where $\#(X)$ denotes the cardinality of $X$. It is a function space for piecewise linear
functions with free knots $\{x_k\}$.
Notice that $V_h^N$ is not
a linear  space, since the summation of two functions in $V_h^N$ may
not belongs to the same space when they correspond to different partitions $X$.

For any function $u_h(x)\in V_h^N$, it can be written as
\begin{equation}\label{e:basis}
u_h(x)=\sum_{k=1}^{N-1} u_k\phi_k(x),
\end{equation}
where $\phi_k(x)$ is the standard nodal basis function with respect to a partition
$\mathcal{T}(X)$, namely
\begin{equation*}
\phi_k(x)=\frac{x-x_{k-1}}{x_{k}-x_{k-1}}\chi_{I_k}(x)+\frac{x_{k+1}-x}{x_{k+1}-x_{k}} \chi_{I_{k+1}}(x),
\end{equation*}
where $\chi_{I_k}$ is the characteristic function corresponding to $I_k$,
\begin{equation*}
\chi_{I_k}(x)=\left\{
\begin{array}{ll}
1& \hbox{if } x \in I_{k} \\
0& \hbox{otherwise}.
\end{array}
\right.
\end{equation*}
In Equation~\eqref{e:basis}, both $u_k$ and $x_k$ ($k=1,\cdots,N-1$) can change their values.

 Therefore,
 $V_h^N$ is a $2(N-1)$ manifold.
We remark that $V_h^N$ is not a smooth manifold in general. There exists some function $u_h\in V_h^N$
which corresponds to many different coordinates $\{u_k, x_k\ |\ k=1,\cdots N-1\}$.
 In this sense, the manifold
 degenerates for some functions in $V_h^N$.

The nonlinear approximation theory for free-knot piecewise linear polynomials has been studied extensively (c.f.  \cite{petrushev1988direct,devore1998nonlinear}
and the reference therein).
For a function $u\in H^1_0$, the best approximation of $u$ in $V_h^N$ in energy norm is defined as
\begin{equation}
\sigma_N(u)=\inf_{v_h\in V_h^N}|u-v_h|_{H^1}.
\end{equation}
Here $|u|_{H^1}:=\Big(\int_{I} (\partial_x u)^2 dx\Big)^{\frac{1}{2}}$ is the standard $H^1$ semi-norm.

To characterize the approximation property of functions in $V_h^N$, it is convenient to use the Besov spaces.
We denote by $B_q^r(L_q(I))$ a standard Besov space.
We will not give the details of the definitions of the Besov space here, but refer to \cite{triebel1978interpolation,devore1993besov}.
We  only mention that $B_q^r(L_q(I))$  is a  space consisting of functions with smoothness order $r$ measured in $L_q$.
When $q=2$, $B_2^r(L_2(I))$ is identical to the Sobolev space $H^r$.
The following lemma is known from the literature(c.f.  \cite
 binev2002approximation,binev2004adaptive}).
\begin{lemma} \label{lem:approx}
If $u\in B^{s+1}_q(L_q(I))$ with $0\leq s\leq 1$ and $1/q<s+1/2$, then we have
\begin{equation}
\sigma_N(u)\leq C N^{-s} |u|_{ B^{s+1}_q(L_q(I))}.
\end{equation}
\end{lemma}
The lemma is a one-dimensional version of  Theorem 9.1 in \cite{binev2004adaptive}.
The proof of the lemma can be found in \cite{binev2002approximation}.
 By this lemma, the best approximation in energy norm of a (one-dimensional) function  by a free-knot piecewise linear function
is of order $O(N^{-1})$.
\subsection{The moving finite element method}
Consider a time dependent function $u_h(t,x)=\sum_{k=1}^{N-1}u_k(t)\phi_k(t,x)$ in $V_h^N$,
where
\begin{equation}
\phi_k(t,x)=\frac{x-x_{k-1}(t)}{x_{k}(t)-x_{k-1}(t)}\chi_{I_k}(x)+\frac{x_{k+1}(t)-x}{x_{k+1}(t)-x_{k}(t)} \chi_{I_{k+1}}(x),
\end{equation}
 is the nodal basis function corresponding to a time dependent partition $\mathcal{T}(X(t))$ and the set $X(t)$ is given by
 $$X(t):=\{a=x_0< x_1(t)<\cdots<x_{N-1}(t)<x_N=b\}.$$
In the partition, the interior notes $x_k(t), k=1,\cdots, N-1,$ may change positions when time $t$ evolves.
In the formula of $u_h(t,x)$, there are $2(N-1)$ time dependent parameters
\begin{equation}
\{u_k(t), x_k(t)\ |\ k=1,\cdots N-1\}.
\end{equation}
We aim to approximate the solution $u$ of the model problem~\eqref{e:modelpb1} by a
discrete function $u_h(t,x)$. For that purpose, we will derive a dynamic equation for $u_k(t)$ and $x_k(t)$
by using the Onsager principle. The derivation is similar to that for continuous problems in Section~2.

We first compute the discrete energy functional and the discrete dissipation function as
follows. Notice that the time derivative and space derivative of $u_h(t,x)$ are given by
\begin{align*}
\partial_t u_h&=\sum_{k=1}^{N-1} (\dot{u}_k(t)\phi_k(t,x)+\dot{x}_k(t)\beta_k(t,x)),\\
\partial_x u_h&=\sum_{k=1}^{N-1} {u}_k(t) \partial_x \phi_k(t,x),
\end{align*}
where
$$\beta_k(t,x)=\frac{\partial u_h}{\partial x_k}=-D_h u_{k-1}\frac{x-x_{k-1}(t)}{x_{k}(t)-x_{k-1}(t)}\chi_{I_k}(x)
-D_h  u_{k}\frac{x_{k+1}(t)-x}{x_{k+1}(t)-x_{k}(t)}\chi_{I_{k+1}}(x),$$
with $D_h  u_k=\frac{u_{k+1}(t)-u_{k}(t)}{x_{k+1}(t)-x_{k}(t)}$.
Then the discrete energy functional $\mathcal E$ with respect to  $u_h$ is calculated by
\begin{equation}
\mathcal E_h(u_1,\cdots,u_{N-1},x_1,\cdots,x_{N-1}):=\mathcal E(u_h)=\sum_{k=1}^{N} \int_{I_k}\frac{\alpha}{2}(\partial_x u_h)^2+F(x,u_h) \mathrm{d}x.
\end{equation}
It is a nonlinear function with respect to  $\{u_k(t), x_k(t)\ |\ k=1,\cdots N-1\}$. The discrete dissipation function is given by
\begin{align
&\Phi_h(u_1,\cdots,u_{N-1},x_1,\cdots,x_{N-1};
\dot u_1,\cdots,\dot{u}_{N-1},\dot{x}_1,\cdots,\dot{x}_{N-1})\nonumber \\
&:=\Phi(\partial_t u_h(t,x))= \frac{\xi}{2}\int_{I}(\dot{u}_k(t)\phi_k(t,x)+\dot{x}_k(t)\beta_k(t,x))^2 dx.\label{e:dissip_h}
\end{align}
It is a quadratic function with respect to  $\{\dot{u}_k(t), \dot{x}_k(t)\ |\ k=1,\cdots N-1\}$.

We apply the Onsager
principle for $\{\dot{u}_k(t), \dot{x}_k(t)\ |\ k=1,\cdots N-1\}$. We minimize the discrete Rayleighian $\mathcal{R}_{h}$
with respect to $\dot{u}_k$ and $\dot{x}_k$:
\begin{equation}
\min_{\dot{u}_k,\dot{x}_k} \mathcal{R}_{h}:=\Phi_h+\dot{\mathcal E}_h.
\end{equation}
The problem is equivalent to  its Euler-Lagrange equation
\begin{equation}\label{e:Onsager1}
\left\{
\begin{array}{ll}
\frac{\partial \Phi_h}{\partial \dot{u}_k}+\frac{\partial \mathcal E_h}{\partial u_k}=0, &k=1,\cdots, N-1;\\
\frac{\partial \Phi_h}{\partial \dot{x}_k}+\frac{\partial \mathcal E_h}{\partial x_k}=0, &k=1,\cdots, N-1.
\end{array}
\right.
\end{equation}
This gives an  ordinary differential system for $\{u_k(t), x_k(t)\ |\ k=1,\cdots N-1\}$.

We now derive the explicit formula for the system \eqref{e:Onsager1}.
Denote  $\mathbf{u}=(u_1,\cdots, u_{N-1})^T$ and
$\mathbf{x}=(x_1,\cdots, x_{N-1})^T$.
Notice that $\Phi_h$ is a quadratic function with respect to $\{\dot{u}_k(t),\dot{x}_k(t)\ |\ k=1,\cdots N-1\}$.
Then the equation~\eqref{e:Onsager1} can be rewritten as
\begin{equation}\label{e:ODE0}
\left(\begin{array}{cc}
\mathbf{A} &\mathbf{B} \\
\mathbf{B} &\mathbf{C}
\end{array}
\right)
\left(
\begin{array}{l}
\dot{\mathbf{u}}\\
\dot{\mathbf{x}}
\end{array}
\right)=
\left(
\begin{array}{l}
{\mathbf{f}}\\
{\mathbf{g}}
\end{array}
\right).
\end{equation}
The right hand side terms are given by
$$\mathbf{f}=(f_1,f_2,\cdots,f_{N-1})^T, \qquad \hbox{with }$$
$$ f_k=-\frac{1}{\xi}\frac{\partial \mathcal E_h}{\partial u_k}=-\frac{1}{\xi}\int_{I}\alpha\partial_x u_h \partial_x \phi_k+f(x,u_h)\phi_k dx.
$$
and
$$\mathbf{g}=(g_1,g_2,\cdots,g_{N-1})^T, \qquad \hbox{with }$$
$$
g_k=-\frac{1}{\xi}\frac{\partial \mathcal E_h}{\partial x_k}=-\frac{1}{\xi}\Big(\int_{I}\alpha\partial_x u_h \partial_x \beta_k+f(x,u_h)\beta_k dx-\frac{\alpha}{2}[(\partial_x u_h)^2]|_{x_{k}}\Big),
$$
where the jump $[(\partial_x u_h)^2]|_{x_{k}}:=(\partial_x u_h)^2|_{I_{k}}-(\partial_x u_h)^2|_{I_{k-1}}$.
The  blocks $\mathbf{A}$, $\mathbf{B}$ and $\mathbf{C}$ in the coefficient matrix are $(N-1)\times(N-1)$ tridiagonal matrices, whose
nonzero elements  are functions of $\{u_k\}$ and $\{x_k\}$.
Direct computations for the nonzero elements of $\mathbf{A}$ give
\begin{align*}
&a_{k,k}=\int_{I_k\cup I_{k+1}} \phi_k^2 dx,\qquad \qquad\qquad\quad \ \  k=1,\cdots,N-1;\\
&a_{k,k+1}=a_{k+1,k}=\int_{I_{k+1}} \phi_k\phi_{k+1} dx, \qquad k=1,\cdots,N-2.
\end{align*}
This implies  $\mathbf{A}$ is the mass matrix in the standard finite element method on
a given triangulation $\mathcal{T}(X)$ with $X=\{a=x_0<x_1<\cdots<x_N=b\}$.
The nonzero elements of the matrix $\mathbf{B}$ are computed as
\begin{align*}
&b_{k,k}=\int_{I_k\cup I_{k+1}} \phi_k\beta_k dx,\qquad \qquad\qquad\qquad\quad \  k=1,\cdots,N-1;\\
&b_{k,k+1}=b_{k+1,k}=\int_{I_{k+1}}\phi_k\beta_{k+1}dx,\qquad\qquad \ \  k=1,\cdots,N-2.
\end{align*}
We can also compute the nonzero elements of the matrix $\mathbf{C}$,
\begin{align*}
&c_{k,k}=\int_{I_k\cup I_{k+1}}\beta_k^2 dx, \qquad \qquad\qquad\qquad\quad\  k=1,\cdots,N-1;\\
&c_{k,k+1}=c_{k+1,k}=\int_{I_{k+1}}\beta_k\beta_{k+1}dx,\qquad \qquad   \!  k=1,\cdots,N-2.
\end{align*}

The ordinary differential system~\eqref{e:ODE0} is exactly the same as the moving finite element scheme  in~\cite{miller1981movingA}.
The  scheme was originally derived through a $L^2$ projection of a partial differential equation in the tangential space of $u_h$.
There one needs to compute the inner product of a Dirac function and a discontinuous function, which is not well-defined even in
a distribution sense.
Our derivation is much simpler than that in \cite{miller1981movingA}.  In this sense, the Onsager
principle gives a natural variational framework for the moving finite element method.

\subsection{The stabilized scheme}
The coefficient matrix in \eqref{e:ODE0} may degenerate for some function $u_h\in V_h^N$. This is because
the dissipation function $\Phi_h$ is a semi-positive definite quadratic form with respect to $\dot{u}_k$ and $\dot{x}_k$.

The semi-positive definiteness is related to the degeneracy of  the manifold $V_h^N$.
For example,  if we set
$u_0(t)=u_1(t)=\cdots=u_{N}(t)=0,$
then $u_h(x,t)\equiv 0$ for all possible choice of $X$.
That implies $\partial_t u_h\equiv 0$ and  $\Phi_h\equiv 0$  for some nonzero $\dot{x}_{k}$.

 To overcome the degeneracy of the system \eqref{e:ODE0}, we
add a stabilized term as follows,
\begin{equation}
\Phi_h^{\delta}=\Phi_h+ \frac{\delta\xi}{2} \sum_{k=1}^{N-1}\dot{x}_k^2,
\end{equation}
where $\delta>0$ is a small stabilization parameter.
Other stabilizations can also apply, c.f. \cite{miller1981movingB}.
With the modified dissipation function $\Phi^{\delta}_h$,  by the Onsager principle,
$\dot{u}_k$ and $\dot{x}_k$ are obtained by
\begin{equation}
\min_{\dot{u}_k,\dot{x}_k} \Phi_h^{\delta} +\dot{\mathcal E}_h.
\end{equation}
This leads to a modified  system
\begin{equation}\label{e:Onsager2}
\left\{
\begin{array}{ll}
\frac{\partial \Phi_h^\delta}{\partial \dot{u}_k}+\frac{\partial \mathcal E_h}{\partial u_k}=0,&k=1,\cdots, N-1;\\
\frac{\partial \Phi_h^\delta}{\partial \dot{x}_k}+\frac{\partial \mathcal E_h}{\partial x_k}=0,&k=1,\cdots, N-1.
\end{array}
\right.
\end{equation}
Once again it is an ordinary differential system for $\{u_k(t), x_k(t)\ |\ k=1,\cdots N-1\}$.
The explicit form of the system \eqref{e:Onsager2} is
\begin{equation}\label{e:ODE}
\left(\begin{array}{cc}
\mathbf{A} &\mathbf{B} \\
\mathbf{B} &\mathbf{C}+\delta \mathbf{I}
\end{array}
\right)
\left(
\begin{array}{l}
\dot{\mathbf{u}}\\
\dot{\mathbf{x}}
\end{array}
\right)=
\left(
\begin{array}{l}
{\mathbf{f}}\\
{\mathbf{g}}
\end{array}
\right).
\end{equation}
The positive definiteness of the coefficient matrix
$\mathbf{M}_{\delta}:=\left(\begin{array}{cc}
\mathbf{A} &\mathbf{B} \\
\mathbf{B} &\mathbf{C}+\delta \mathbf{I}
\end{array}\right)$
is given by the following lemma.
\begin{lemma}\label{lem:posdef}
For any given $\delta>0$, $X(t)=\{a=x_0<x_1(t)<\cdots<x_{N-1}(t)<x_N=b\}$, and $\{u_k(t) | k=1,\cdots, N-1\}$,
the coefficient matrix $\mathbf{M}_{\delta}$ is positive definite.
\end{lemma}
\begin{proof}
For any $\mathbf{y}\in R^{2(N-1)}$, we can denote it as $\mathbf{y}=\left(
\begin{array}{l}
{\mathbf{y}_1}\\
{\mathbf{y}_2}
\end{array}
\right)$, where $\mathbf{y}_i\in R^{N-1}$, $i=1,2$. We suppose $\mathbf{y}\neq 0$.
If $\mathbf{y}_2\neq 0$, we can have
\begin{equation}
\mathbf{y}^T \mathbf{M}_{\delta} \mathbf{y}=  \mathbf{y}^T \mathbf{M}_{0} \mathbf{y}+\delta \mathbf y_2^T\mathbf{y}_2
\geq \delta \mathbf y_2^T\mathbf{y}_2>0.
\end{equation}
where $ \mathbf{M}_{0} = \left(\begin{array}{cc}
\mathbf{A} &\mathbf{B} \\
\mathbf{B} &\mathbf{C}
\end{array}\right)$ and we have used the fact that
$$\mathbf{y}^T \mathbf{M}_{0} \mathbf{y}=\Phi_h(u_1,\cdots,u_{N-1},x_1,\cdots,x_{N-1};\mathbf{y}_1,\mathbf{y}_2)$$
 is semi-positive definite by its definition~\eqref{e:dissip_h}.
Otherwise, if $\mathbf{y}_2=0$, then $\mathbf{y}_1\neq 0$. We have
\begin{equation}
\mathbf{y}^T \mathbf{M}_{\delta} \mathbf{y}=  \mathbf{y}_1^T \mathbf{A} \mathbf{y}_1>0,
\end{equation}
for any given  partition $\mathcal{T}(X(t))$ of the interval $I$, since the matrix $\mathbf{A}$ is
the standard mass matrix of the linear finite element method with respect to the partition $\mathcal{T}(X(t))$
and thus positive definite.
Combine the analysis together, we show that  $\mathbf{M}_{\delta}$ is positive definite.
\end{proof}

By this lemma,  the ODE system \eqref{e:ODE} can be rewritten as
\begin{equation}\label{e:ODE1}
\left(
\begin{array}{l}
\dot{\mathbf{u}}\\
\dot{\mathbf{x}}
\end{array}
\right)=
\mathbf{M}_{\delta}^{-1}
\left(
\begin{array}{l}
{\mathbf{f}}\\
{\mathbf{g}}
\end{array}
\right).
\end{equation}
The  equation has a unique solution  for any initial value $\mathbf{u}(0)$ and $\mathbf{x}(0)$ when the right hand side
function $\left(
\begin{array}{l}
{\mathbf{f}(\mathbf{u},\mathbf{x})}\\
{\mathbf{g}(\mathbf{u},\mathbf{x})}
\end{array}
\right)$ is Lipschitz continuous with respect to $\mathbf{u}$ and $\mathbf{x}$.
The stiffness of the
  system and its numerical solution has been analyzed in \cite{miller1981movingB,wathen1985structure}.
When the numerical solution is non-degenerate in $V_h^N$, the stabilization parameter $\delta$ can be chosen as zero.

\section{Numerical analysis}
One important advantage of nonlinear approximations is that they have better  accuracy
than linear approximations. This is illustrated in Lemma~\ref{lem:approx} for the interpolation error.
One would expect that the MFEM also has better accuracy than the standard FEM.
However, the theoretical analysis is very difficult.
In this section, we consider only the stationary solution when $t$ goes to infinity.

\subsection{The energy decay property}
We first prove the following discrete energy decay property of the ordinary differential system~\eqref{e:ODE}.
\begin{theorem}\label{th:decay}
Let $(\mathbf{u}(t),\mathbf{x}(t))$ be
the solution of the equation \eqref{e:ODE} and $u_h(t,x)\in V_h^N$ be the corresponding piecewise linear function.
Then we have
\begin{equation}
\frac{d \mathcal E(u_h)}{d t}
\leq 0,
\end{equation}
where the equality holds only when $\dot{x}_k=0$ and $\dot{u}_k=0$ for $k=1,\cdots,N-1$.
\end{theorem}
\begin{proof}
Notice that
$$\frac{d \mathcal E(u_h)}{d t}=\sum_{k=1}^{N-1}(\frac{\partial \mathcal E_h }{ \partial u_k}\dot{u}_k+\frac{\partial \mathcal E_h }{\partial x_k}\dot{x}_k)=
-\xi \sum_{k=1}^{N-1}(f_k\dot{u}_k+g_k\dot{x}_k)=-\xi (\dot{\mathbf{u}}^T,\dot{\mathbf{x}}^T) \mathbf M_{\delta} \left(
\begin{array}{l}
\dot{\mathbf{u}}\\
\dot{\mathbf{x}}
\end{array}
\right)
$$
where we have used the equation \eqref{e:ODE}.
By Lemma~\ref{lem:posdef}, we finish the proof of the theorem.
\end{proof}

From the proof, we easily see that
\begin{equation}\label{e:disStru}
\frac{d \mathcal E(u_h)}{d t}=-2\Phi_h^{\delta}.
\end{equation}
This implies the semi-discrete problem~\eqref{e:ODE} preserves the
energy dissipation structure \eqref{e:engdecay1} of the original gradient flow system.

\subsection{Error analysis for the equilibrium state}
We now consider the stationary solution of \eqref{e:modelpb1}.
When $t$ goes to infinity, we suppose the solution $u(t,x)$ of \eqref{e:modelpb1} converge to a function $u^\infty(x)$
satisfying
\begin{equation}\label{e:EL_S}
\int_I \alpha \partial_x u^\infty\partial_x v+ f(x,u^\infty) v \mathrm{d}x=0,\qquad \forall v\in H^1_0(I).
\end{equation}
 To do analyse for this problem, we need one more assumption for the Lipschitz constant $L_0$ in \eqref{e:Lip} that
\begin{itemize}
\item[(H1)]{\qquad\qquad\qquad\qquad ${\alpha}- L_0 c_0^2 >c_1$, \qquad\ \  for some constant $c_1>0$,}
\end{itemize}
where $c_0$ is a constant in the Poincare inequality
\begin{equation}\label{e:Sob}
\|u\|_{L^2(I)}\leq c_0|u|_{H^1(I)}.
\end{equation}
Under this assumption, it is easy to see that the equation~\eqref{e:EL_S} is elliptic and has a unique solution.
Furthermore, $u^{\infty}$ is the unique minimizer for the energy minimization problem
\begin{equation}\label{e:conEng}
\inf_{v\in H^1_0} \mathcal E(v).
\end{equation}

In the following, we consider the approximation of the problem~\eqref{e:EL_S} by the MFEM.
For the stationary solution for  the discrete problem~$\eqref{e:ODE}$, we have the following theorem.
\begin{theorem}\label{th:engRel}
Let $(\mathbf{u}(t),\mathbf{x}(t))$ be
the solution of the equation \eqref{e:ODE}  which corresponds to a nondegenerate function $u_h(t,x)\in V_h^N$.
 When $t$ goes to infinity,  $(\mathbf{u}(t),\mathbf{x}(t))$ will converge to a vector $(\mathbf{u}^{\infty},\mathbf{x}^{\infty})$.
 Suppose that $\mathbf{x}^{\infty}$  corresponds to a partition of $I$, then
$u_h(t,x)$ converges to a function $u_h^{\infty}(x)\in V_h^N$ which satisfies
\begin{equation}\label{e:discreteEL}
\left\{
\begin{array}{l}
\int_{I}\alpha\partial_x u_h^{\infty} \partial_x \phi_k+f(x,u_h^{\infty})\phi_k dx=0, \qquad\qquad\qquad\qquad\quad\  k=1,\cdots, N-1,\\
\int_{I}\alpha\partial_x u_h^{\infty} \partial_x \beta_k+f(x,u_h^{\infty})\beta_k dx-\frac{\alpha}{2}[(\partial_x u_h^{\infty} )^2]|_{x_{k}^{\infty}}=0, \qquad k=1,\cdots, N-1.
\end{array}
\right.
\end{equation}
Furthermore, if $u_h^{\infty}$ is nondegenerate and asymptotic stable, it is a local minimizer of  $\mathcal E(u_h)$ in $V_h^N$.
\end{theorem}

\begin{proof}
By Equation~\eqref{e:disStru}, we have the energy decay property
$$
\frac{d \mathcal E(u_h(t))}{dt}=-2\Phi(\partial_t u_h)-\delta\xi \sum_{k=1}^N\dot{x}_k^2\leq 0, \qquad\qquad \forall t>0.
$$
Therefore, we have
\begin{align}\label{e:temp1}
\mathcal E(u_h(T))=\mathcal E(u_h(0))-\xi\int_0^T\int_I(\partial_t u_h)^2dx dt-\delta\xi \sum_{k=1}^N\int_0^T\dot{x}_k^2dt, \qquad \hbox{for } T>0.
\end{align}
Notice that
$$
\mathcal E(u_h(T))\geq \inf_{v_h\in V_h^N}\mathcal E(v_h)\geq\inf_{v\in H^1_0(I)}\mathcal E(v)=\mathcal E(u^{\infty}).
$$
Thus $\mathcal E(u_h(T))$ is bounded from below when $T$ goes to infinity.
By Equation~\eqref{e:temp1}, we have $\lim_{t\rightarrow \infty}\dot{x}_k(t)=0$ and $\lim_{t\rightarrow \infty}\int_I(\partial_t u_h)^2dx=0$.
Since $x_k$ is in the bounded interval $I$, $x_k(t)$ converges when $t$ goes to infinity.
Similarly, $u_k(t)$ also converges when $t$ goes to infinity.
When the limit $\mathbf{x}^{\infty}$ generates a partition of $I$,
then $(\mathbf{u}^{\infty},\mathbf{x}^{\infty})$ corresponds to a function $u_h^{\infty}$ in $V_h^N$.
From Equation~\eqref{e:ODE}, we have
\begin{equation*}
\lim_{t\rightarrow\infty} f_k(\mathbf{u},\mathbf{x}) =0, \quad \lim_{t\rightarrow\infty} g_k(\mathbf{u},\mathbf{x})=0, \qquad \hbox{for } k=1,\cdots,N-1.
\end{equation*}
This leads to the equation~\eqref{e:discreteEL}, i.e.
$u_h^\infty$ is a critical point of the  energy $\mathcal E(\cdot)$ in the manifold $V_h^N$.

When $u_h^{\infty}$ is nondegenerate and asymptotic stable,  it is a local minimizer of $\mathcal E$ in $V_h^N$.
Otherwise, if it is a critical point but not a local minimizer of $\mathcal E$, there exists a non-trivial trajectory near $u_h^{\infty}$ which makes
the energy nondecreasing. This contradicts with the energy decay property in Theorem~\ref{th:decay}.
\end{proof}

\begin{remark}
If $\mathcal E(u)$ has a quadratic form, Equation~\eqref{e:EL_S} is a linear equation.  Then the theorem implies that $u_h^{\infty}$ is
a locally best approximation in energy norm. For example, suppose
$$\mathcal E(u)=\int_{I}\frac{\alpha}{2}(\partial_x u)^2+f(x) u \mathrm{d}x,$$
then by  Equation~\eqref{e:EL_S}, we have
\begin{align*}
\mathcal E(v_h)&= \int_{I}\frac{\alpha}{2}(\partial_x v_h)^2+f(x) v_h \mathrm{d}x\\
&=\int_{I}\frac{\alpha}{2}(\partial_x v_h)^2- \alpha \partial_x u^\infty \partial_x v_h \mathrm{d}x\\
&=\int_{I}\frac{\alpha}{2}(\partial_x v_h-\partial_x u^\infty )^2dx -\int_{I}\frac{\alpha}{2} (\partial_x u^\infty)^2 \mathrm{d}x.
\end{align*}
Therefore, Theorem~\ref{th:engRel} implies that $u_h^\infty$ is a locally best approximation of $u^\infty$ in the energy norm
when $u_h^{\infty}$ is nondegenerate and asymptotic stable.
This is one main result obtained in \cite{jimack1996best}.
\end{remark}

Notice the admissible set $V_h^N$ for the discrete function is a $2(N-1)$ dimensional manifold in $H^1_0(I)$.
In general the energy minimization problem $\mathcal E(\cdot)$ in the manifold $V_h^N$ is not convex, even when $\mathcal E(\cdot)$ is
a convex functional. This is very different from the standard FEM, noticing that the minimizer of
a convex functional in a linear space is unique.
Therefore, the local best approximations given in Theorem~\ref{th:engRel} is optimal.
However, when the initial value
is in the domain of attraction of a global minimizer of $\mathcal E(\cdot)$ in $V_h^N$, one can approximate the global minimizer
by the MFEM.
The following theorem shows that the global minimizer has an optimal convergence order.

\begin{theorem}\label{theo:main}
Let $u_h^{\infty}$  be a global minimizer of $\mathcal E(\cdot)$ in $V_h^N$. Assume that $f(x,u)$ is differentiable with respect to $u$.
 If $u^\infty \in B^{s+1}_q(L_q(I))$ with $0\leq s\leq 1$ and $1/q<s+1/2$,
then we have
\begin{equation}
 |u^{\infty}-u_h^{\infty}|_{H^1}\lesssim N^{-s} |u^\infty|_{ B^{s+1}_q(L_q(I))}.
 \end{equation}
\end{theorem}
\begin{proof}
Denote $w= u_h^{\infty}-u^{\infty}$. Since $u_h^\infty\in V_h^N\subset H^1_0(I)$ and $u^{\infty}$ is the minimizer of $E(u)$ in $H^1_0(I)$,  then we have
$w \in H^1_0(I)$ and
\begin{align}
&\mathcal E(u_h^{\infty})-\mathcal E(u^{\infty}) \nonumber\\
&=\int_{I}\Big(\frac{\alpha}{2}(\partial_x u_h^\infty )^2-\frac{\alpha}{2}(\partial_x u^\infty)^2\Big)+\Big(F(x,u_h^{\infty})
-  F(x,u^{\infty})\Big)\mathrm{d}x\nonumber\\
&=\int_{I}\frac{\alpha}{2}(\partial_x w)^2+\alpha\partial_x u^\infty \partial_x w+
 f(x,u^{\infty}) w +\frac{1}{2}\partial_{u}f(x,u^\infty+s w)w^2 \mathrm{d}x\nonumber\\
&=\int_{I}[\alpha\partial_x u^\infty \partial_x w+
 f(x,u^\infty) w ]+\frac{\alpha}{2}(\partial_x w)^2+\frac{1}{2}\partial_{u}f(x,u^\infty+s w)w^2 \mathrm{d}x\nonumber\\
&=\int_{I}\frac{\alpha}{2}(\partial_x w)^2+\frac{1}{2}\partial_{u}f(x,u^\infty+s w)w^2 \mathrm{d}x,\label{e:temp2}
\end{align}
where $0\leq s\leq 1$ and we have used the equation~\eqref{e:EL_S} in the last equation. From the Lipschitz condition~\eqref{e:Lip} of $f$,
we know that $|\partial_{u}f(x, u)|\leq L_0$. This further leads to
\begin{align}
\mathcal E(u_h^{\infty})-\mathcal E(u^{\infty})&\geq \int_{I}\frac{\alpha}{2}(\partial_x w)^2\mathrm{d}x -\frac{L_0}{2}\int_I w^2 \mathrm{d}x \nonumber \\
&\geq \frac{\alpha-L_0 c_0^2}{2} \int_I (\partial_x w)^2\mathrm{d}x\nonumber \\
&\geq \frac{c_1}{2} \int_I (\partial_x w)^2\mathrm{d}x=\frac{c_1}{2}|\partial_x  u^{\infty}-\partial_x u^{\infty}_h |^2_{H^1}.\label{e:temp}
\end{align}
Here we have used  assumption (H1) and the Poincare inequality~\eqref{e:Sob}.

For any $v_h\in V_h^N$,  we denote $\tilde{w}=v_h-u^\infty$. Noticing that $u_h^{\infty}$ is a global minimizer
of $\mathcal E$, similar calculations as \eqref{e:temp2} gives
\begin{align*}
\mathcal E(u_h^{\infty})-\mathcal E(u^{\infty})&\leq  \mathcal E(v_h)-\mathcal E(u^{\infty}) \\
&=\int_{I}\frac{\alpha}{2}(\partial_x \tilde{w})^2+\frac{1}{2}\partial_{u}f(x,u^\infty+s\tilde{w})\tilde{w}^2 \mathrm{d}x\\
&\lesssim  \|u^{\infty}-v_h\|_{H^1}^2,\qquad \forall v_h\in V_h^N.
\end{align*}
Combine  the inequality with \eqref{e:temp} and using the Poincare inequality again, we obtain
\begin{equation}
 |u^{\infty}-u_h^{\infty}|_{H^1}^2\lesssim \mathcal E(u_h^{\infty})-\mathcal E(u^{\infty})\lesssim \inf_{v_h\in V_h^N}\|u^{\infty}-v_h\|_{H^1}^2 \lesssim (\sigma_N(u))^2.
\end{equation}
By Lemma~\ref{lem:approx}, we have finished the proof.

\end{proof}
In general it is not clear how to determine the domain of  attraction of a global minimizer.
There is no guarantee that  all initial values can tend to a  globally stationary solution.
As suggested by the numerical results in \cite{jimack1996best}, it is possible that the
solution of a MFEM does not go to the global minimizer when some elements in the partition
tend to degenerate.

\section{Numerical Examples}
As we mentioned in Section 3, the equation~\eqref{e:ODE} is a stiff ODE system.
One reason is that some mesh size in the MFEM can be very small.
In \cite{miller1981movingB}, a technique by adding ``internodal spring force" was
 developed to keep the grid points at least slightly separated.
Here we use a  similar technique by adding a penalty term to the total energy.
We replace $\mathcal E(v_h)$ by
\begin{equation}\label{e:modifyEng}
\mathcal E_{\tilde{\delta}} (v_h)=\mathcal E(v_h)+\frac{\tilde{\delta}}{N}\sum_{k=0}^N(\ln(N h_k/L))^2.
\end{equation}
where $h_k=x_{k+1}-x_k$ and $\tilde{\delta}$ is a small parameter.
The penalty term goes to infinity when $h_k$ goes to zero.

In the following, we show some numerical examples. In our simulations, we solve the problem \eqref{e:ODE} by a simple
forward Euler scheme with very small time step. More efficient
solvers(e.g. some implicit stiff ODE solver) have been studied in \cite{miller1981movingB}.
\subsection{A linear diffusion equation}
We first consider a linear diffusion problem:
\begin{equation}\label{e:lineardiffusion}
\left\{
\begin{array}{ll}
 \partial_t u-\partial_{xx} u=f(x), & x\in (a,b),t>0;\\
 u(x,0)=u_0(x),& x\in (a,b);\\
 u(a,t)=u(b,t)=0 & t>0.
\end{array}
\right.
\end{equation}
This corresponding energy density is given by
\begin{equation}\label{e:energydiffusion}
\mathcal E(u)=\int_{I}\frac{1}{2}|\partial_x u|^2 -f(x)u(x) \mathrm{d}x,
\end{equation}
where $I=(a,b)$.
We choose $\delta=0.0001$ and $\tilde{\delta}=0.01$. We test for  several different $u_0$ and $f$.

{\bf Example 1.} In the first example, we set $I=(0,1)$ and $f(x)=\delta_{a_0}(x)$
which is a Dirac measure at the middle point $a_0=0.5$.
We have $\int_{0}^1 f(x)v(x)\mathrm{d}x=v(x_0)$, for any $v\in C_0^1([0,1])$.
The initial solution is given by $u_0= \sin(\pi x)+x\chi_{\{x<x_0\}}+(1-x)\chi_{\{x>x_0\}}$.
In this case the equation~\eqref{e:lineardiffusion} has an analytic solution
$$u(x,t)=\sin(\pi x)e^{-\pi^2 t}+x\chi_{\{x<x_0\}}+(1-x)\chi_{\{x>x_0\}}.$$
The spacial derivative of the  solution has a jump  at the middle point $a_0=0.5$.

The initial partition  of $I$ is uniform and given by
$0=x_0<\frac{1}{N}<\cdots<\frac{N-1}{N}<x_N=1$.
Here we choose $N=5,9,19,39,79$ so that the non-smooth point $a_0$ is
not a node in the initial partition.
 We solve the problem \eqref{e:ODE} until $T=0.04$.
 The $H^1$ error and the $L^2$ error between the discrete solution $u_h$ and the exact solution $u$ at $T$
 are computed by
\begin{align*}
&err_{H^1}:=\Big(\int_I(\partial_x u(x,T)-\partial_x u_h(x,T))^2 \mathrm{d}x\Big)^{1/2},\\
 &err_{L^2}:=\Big(\int_I( u(x,T)- u_h(x,T))^2 \mathrm{d}x\Big)^{1/2}.
\end{align*}

\begin{table}[h]\small
\caption{\small The $H^1$-norm and $L^2$-norm of the error in Experiment~1.}\label{tab:ErrorEx1}
\vspace{-0.2cm}
\begin{center}
\begin{tabular}{l|lc|lc}
\hline
{\bf Adaptive} & $err_{H^1}$  & order &  $err_{L^2}$   & order \\
\hline
$N=5$&0.3326 & --  &0.0338 &-- \\
$N=9$ &0.1571& 1.28 & 0.0098&2.11\\
$N=19$ &0.0842 &0.83& 0.0022 &2.00 \\
$N=39$&0.0470& 0.81&0.000523&1.99 \\
$N=79$&0.0265& 0.81&0.000141& 1.86 \\
\hline
\end{tabular}\\
\vspace{0.2cm}
\begin{tabular}{l|lc|lc}
\hline
{\bf Uniform} &  $err_{H^1}$  & order &  $err_{L^2}$   & order  \\
\hline
$N=5$&0.6363 & --  &0.0435 &-- \\
$N=9$ &0.1571& 0.72 & 0.0153&1.77 \\
$N=19$ &0.2587 &0.64& 0.0043 &1.70 \\
$N=39$&0.1703& 0.58&0.0013&1.66 \\
$N=79$&0.1159& 0.54&0.000438& 1.54 \\
\hline
\end{tabular}
\end{center}
\end{table}
In Table~\ref{tab:ErrorEx1}, we show the numerical errors for the different choice of $N$. The convergence order is computed
by $s_i:=\ln(err_i/err_{i+1})/\ln(N_{i+1}/N_{i})$, which implies that
the errors decrease with order $O(N^{-s})$. We can see that the $H^1$-error is almost of order $O(N^{-0.8})$, close to the optimal
convergence rate $O(N^{-1})$. Meanwhile, the $L^2$-error is of order $O(N^{-2})$, which is  optimal in
approximation theory. In comparison, we also show the errors computed by the standard finite element method
 on a uniform partitions of $I$. We could see that the $H^1$ norm is almost of order $O(N^{-0.5})$ and
the $L^2$-error is of order $O(N^{-1.5})$. This is reasonable since the exact solution is not smooth.
In this case, the MFEM has better accuracy than the standard FEM.

In Figure~\ref{fig:Exampl1} we illustrate the numerical solutions for the case $N=9$ at different time $t$.
We can see that initially the interval $I$ is uniformly divided into $N-1$ cells and $u_k$ are simply the interpolation of $u(0,x)$.
There is a large error at $a_0=0.5$ since the initial function is not smooth there. With time increasing,
we can see that the  notes $x_k$ move with time. In particular, the two notes near the non-smooth point $a_0$ approach
to the point gradually. The distance between them become smaller and smaller. At about $t=0.0048$, the right point almost arrives at
$a_0$. After that, the left point moves away from the point to further decrease the computational error.
 \begin{figure}[ht!]
 \centering
  \subfigure[$t=0$]{

   \includegraphics[width=0.31\textwidth]{diffusion0.eps}
   }
  \subfigure[$t=0.0012$]{

   \includegraphics[width=0.31\textwidth]{diffusion0_0012.eps}
}
  \subfigure[$t=0.0024$]{

   \includegraphics[width=0.31\textwidth]{diffusion0_0024.eps}
}
  \subfigure[$t=0.0036$]{

   \includegraphics[width=0.31\textwidth]{diffusion0_0036.eps}
}
  \subfigure[$t=0.0048$]{

   \includegraphics[width=0.31\textwidth]{diffusion0_0048.eps}
}
  \subfigure[$t=0.006$]{

   \includegraphics[width=0.31\textwidth]{diffusion0_0060.eps}
}
  \subfigure[$t=0.009$]{

   \includegraphics[width=0.31\textwidth]{diffusion0_0090.eps}
}
  \subfigure[$t=0.014$]{

   \includegraphics[width=0.31\textwidth]{diffusion0_015.eps}
}
  \subfigure[$t=0.04$]{

   \includegraphics[width=0.31\textwidth]{diffusion0_040.eps}
}
 \caption{The profile of the numerical solution and the exact solution at different time $t$ when $N=9$ for Example 1.
 The small circles in the figures show the computational solutions $u_k$ corresponding to $x_k$, for $k=0,\cdots,9$.
 }
 \label{fig:Exampl1}
\end{figure}

{\bf Example 2.} To further illustrate the behaviour of the MFEM, we consider a diffusion problem with
an almost singular initial condition. We set $f=0$ and $I=(-3,3)$.
 We choose $u_0=\frac{1}{\sqrt{0.004\pi}}(e^{-x^2/(0.004)}-e^{-9/(0.004)})$, which is a very sharp
Gaussian function as shown in the first sub-figure in Figure~\ref{fig:Example2}.  $u_0$ is almost equal to
zero outside of a narrow interval. For this problem, we do not have an explicit formula for the analytic solution.

Initially we choose a grid whose points concentrated in a small interval $(-0.2,0.2)$ to resolve the profile of $u_0$:
\begin{equation*}
X(0)=\{-3,-0.2,\cdots,-0.2+i\frac{0.4}{N-2},\cdots,0.2,3\}.
\end{equation*}

In Figure~\ref{fig:Example2}, we show the numerical solution when $N=19$.
We can see that the numerical solutions spread outwards gradually while
the height of the profile decreases.
In particular, the grid points gradually adjust their positions
so that  the Gaussian profile of the numerical solution is always well resolved.
 \begin{figure}[ht!]
 \centering
  \subfigure[$t=0$]{

   \includegraphics[width=0.45\textwidth]{kernal0_0.eps}
   }
  \subfigure[$t=0.001$]{

   \includegraphics[width=0.45\textwidth]{kernal0_001.eps}
}
  \subfigure[$t=0.005$]{

   \includegraphics[width=0.45\textwidth]{kernal0_005.eps}
}
  \subfigure[$t=0.02$]{

   \includegraphics[width=0.45\textwidth]{kernal0_02.eps}
}
  \subfigure[$t=0.1$]{

   \includegraphics[width=0.45\textwidth]{kernal0_1.eps}
}
  \subfigure[$t=0.2$]{

   \includegraphics[width=0.45\textwidth]{kernal0_2.eps}
}

 \caption{The profile of the numerical solution at different time $t$ when $N=19$ for Example 2.
 The small circles in the figures show the computational solutions $u_k$ corresponding to $x_k$, for $k=0,\cdots,19$.}
 \label{fig:Example2}
\end{figure}
\subsection{The Allen-Cahn equation}
We  then consider an Allen-Cahn equation in a interval $I=(0,1)$,
\begin{equation}\label{e:AC}
\left\{
\begin{array}{ll}
\partial_t \phi -\varepsilon\partial_{xx} \phi +\frac{\phi^3-\phi}{\varepsilon}=0,& x\in (0,1),\\
\phi(x,0)=x, & x\in (0,1),\\
\phi(0,t)=0, \phi(1,t)=1,
\end{array}
\right.
\end{equation}
 which is a nonlinear example for the model problem~\eqref{e:modelpb1}.
The corresponding free energy is
\begin{equation}
\mathcal E(\phi)=\int_0^1\frac{\varepsilon}{2}(\partial_x \phi)^2 +\frac{(1-\phi^2)^2}{4\varepsilon} \mathrm{d}x.
\end{equation}
We are interested in the stationary profile of the problem. When $\varepsilon$ is small enough, the  profile is approximated nicely by
\begin{equation}
\phi^{\infty}(x)=\tanh(\frac{x-0.5}{\sqrt{2}\varepsilon}).
\end{equation}
It has a very thin inner layer at the middle point $x=0.5$.
In simulations, we choose $\delta=0.0001$, $\tilde{\delta}=0.0001$ and test for  different values of $\varepsilon$.
The initial condition is given by
$$\phi(0,x)=2(x-0.5).$$

{\bf Example 3.} In this example, we set $\varepsilon=0.05$ and test the method for different choice of $N$.
Initially, we give a uniform partition of $I$ that
 $X(0)=\{0=x_0<\frac{1}{N}<\cdots<\frac{N-1}{N}<x_N=1\}$,

We solve the problem for various choices of $N=5,10,20,40,80$.

In the computation, the profile of $\phi_h$ changes gradually to a stationary state $\phi_h^\infty$.
Meanwhile the grid points move
into a highly non-uniform distribution to resolve the inner layers of the solution.
We illustrate some plots of the final solution $\phi_h^{\infty}$ as well as the distributions of $x_k$ in Figure~\ref{fig:Example3}.
They are obtained by solving the problem \eqref{e:ODE} until the decrease of the discrete energy
 in one step is smaller than a tolerance $TOL=1e^{-10}$.
The profile of the ``exact" solution $\phi^{\infty}$ is also shown.
 We can see that the numerical solution fits the exact solution very well
 even on a very coarse mesh($N=5$).
Almost all the grid points are distributed in the inner layer of the solution.
Interestingly,  it seems that the nodes are dense where the second order derivative
of the solution is  large.  This indicates a uniform distribution of the errors in $H^1$ norm.

 \begin{figure}[ht!]\label{fig:Example3}
 \centering
  \subfigure[$N=5$]{

   \includegraphics[width=0.22\textwidth]{AC05_5.eps}
}
  \subfigure[$N=10$]{

   \includegraphics[width=0.22\textwidth]{AC05_10.eps}
}
  \subfigure[$N=20$]{

   \includegraphics[width=0.22\textwidth]{AC05_20.eps}
}
  \subfigure[$N=40$]{

   \includegraphics[width=0.22\textwidth]{AC05_40.eps}
}
 \caption{The numerical solutions of the Allen-Chan equation in stationary sate with $\varepsilon=0.05$
 for various $N=5,10,20,40$.   The small circles in the figures mark the points $(x_k,\phi_{h,k})$.
 }
\end{figure}

In Table~\ref{tab:ErrorEx3}, we illustrate the $H^1$~error between $\phi_h^\infty$ and $\phi^{\infty}$:
\begin{equation*}
err_{H^1}:=\Big(\int_I(\partial_x \phi^\infty(x)-\partial_x \phi_h^{\infty}(x))^2 \mathrm{d}x\Big)^{1/2}.
\end{equation*}
It seems that the $H^1$~errors decrease with order $O(N^{-1})$.
This is optimal for free-knot piecewise linear functions.
This is consistent with the analytical results in the Section 4.
In addition, we also compute the errors of the minimal energies by computing
\begin{equation*}
err_{eng}=|\mathcal E(\phi_h^\infty)-\sigma|,
\end{equation*}
where $\sigma=\frac{2\sqrt{2}}{3}$ is the total energy corresponding to the profile $\phi^{\infty}$.
We  see  that the errors of the energy decrease with order $O(N^{-2})$.

\begin{table}[h]\small
\caption{\small The $H^1$ errors of the stationary profile and the errors of the minimial energies in Example~3.}\label{tab:ErrorEx3}
\vspace{-0.2cm}
\begin{center}
\begin{tabular}{l|lc|lc}
\hline
{\bf Adaptive} & $err_{H^1}$  & order &  $err_{eng}$   & order \\
\hline
$N=5$&1.0004 & --  &0.0264 &-- \\
$N=10$ &0.5025& 0.99  & 0.00646&2.03\\
$N=20$ &0.2641 & 0.93 & 0.00175 &1.88 \\
$N=40$&0.1268& 1.06 &0.000402&2.12 \\
$N=80$&0.0691& 0.88&0.000120& 1.87 \\
\hline
\end{tabular}
\end{center}
\end{table}

{\bf Example 4.} In the last example, we test for a smaller $\varepsilon=0.01$, which corresponds to a sharper inner layer.
All other setups are the same as in the previous example.

In Figure~\ref{fig:Example4}, we show the the stationary solutions of the Allen-Cahn equation for various different
choice of $N$. We see that the very sharp inner layer is resolved nicely by the moving finite element method even when
$N=5$.
In Table~\ref{tab:ErrorEx4}, we show the $H^1$~errors between $\phi_h^{\infty}$ and $\phi^\infty(x)$
and also the errors of the minimal energy. We  see that the $H^1$ errors decrease with order $O(N^{-1})$
and the errors of the energies are of order $O(N^{-2})$. Both of them are optimal.
This again verifies the theoretical results in the previous section.

 \begin{figure}[ht!]\label{fig:Example4}
 \centering
  \subfigure[$N=5$]{

   \includegraphics[width=0.22\textwidth]{AC01_5.eps}
}
  \subfigure[$N=10$]{

   \includegraphics[width=0.22\textwidth]{AC01_10.eps}
}
  \subfigure[$N=20$]{

   \includegraphics[width=0.22\textwidth]{AC01_20.eps}
}
  \subfigure[$N=40$]{

   \includegraphics[width=0.22\textwidth]{AC01_40.eps}
}

 \caption{The numerical solutions of the Allen-Chan equation in stationary sate with $\varepsilon=0.01$
 for various $N=5,10,20,40$.   The small circles in the figures mark the points $(x_k,\phi_{h,k})$.
 }
\end{figure}

\begin{table}[h]\small
\caption{\small The $H^1$ errors of the stationary profile and the errors of the minimial energies in Example~4.}\label{tab:ErrorEx4}
\vspace{-0.2cm}
\begin{center}
\begin{tabular}{l|lc|lc}
\hline
{\bf Adaptive} & $err_{H^1}$ & order &  $err_{eng}$   & order \\
\hline
$N=5$&2.2779 & --  &0.0281 &-- \\
$N=10$ &1.7365& 0.39 & 0.0067&2.07\\
$N=20$ &0.8043 &1.11 & 0.0017 &1.98 \\
$N=40$&0.3728& 1.11&0.000467&1.86 \\
$N=80$&0.1641& 1.18&0.000135& 1.79 \\
\hline
\end{tabular}
\end{center}
\end{table}

\section{Conclusions}
We show that the Onsager principle
can act as a  variational framework for the MFEM for a gradient flow system.
The derivation of the method using the principle is much easier than the original approach in \cite{miller1981movingA}.
The discrete problem has the same energy dissipation structure as the
continuous one.
This helps us to do numerical analysis for
a long time(stationary) solution of the gradient flow system. Under some conditions,
we prove that the MFEM gives locally best approximation for the energy.
The optimal convergence rate can be obtained if a global minimizer is detected in
the free-knot piecewise linear function space. Although we consider only the one dimensional
case, the analysis can be generated to high dimensional
problems directly. In addition, it is also possible to consider more complicated gradient flows.

 In this paper, we restrict our analysis on  the stationary solution, although
 numerical experiments show that the method has optimal convergence rate
  for time dependent solutions
 (at least) for the linear diffusion problem.
Numerical analysis for the dynamic problems  will be left for future study.

\section*{Acknowledgment}
This work was supported in part by NSFC grants DMS-11971469
 and  the National Key R\&D Program of China under Grant 2018YFB0704304 and Grant 2018YFB0704300.
\bibliographystyle{abbrv