\section{Introduction}

Reinforcement Learning (RL) has achieved great successes in games and robotics control. For example, the agents trained with RL on StarCraft and Go have surpassed the top level of human~\cite{silver2018a,vinyals2017starcraft}.

Deep RL algorithms such as Deep Q Network (DQN)~\cite{mnih2013playing} usually apply dithering exploration strategy, such as $\varepsilon$-greedy~\cite{osband2016deep}, noise injection to actions~\cite{lillicrap2016continuous}, and action-level regularization~\cite{williams1992simple}. These exploration strategies are all based on local perturbation of actions, hence are not likely to lead to large scale behavior needed for efficient exploration~\cite{osband2019deep}, which is a must for real-world problem such as auto-driving due to enormously large search space.

Recent developments inspired by the duality between control and inference~\cite{todorov2008general} have shown improvement on exploration efficiency, such as VIME~\cite{houthooft2016vime}, BBQ Network~\cite{blundell2015weight}, Bootstrapped DQN~\cite{osband2016deep}, NoisyNet~\cite{Fortunato2017Noisy}, and VDQN~\cite{tang2017variational}. Among these works, either parameters of the value/policy function or the parameters that govern the distribution of the value/policy function are considered as random variables. Variational inference technique is applied to approximate the posterior of these parameters. Experimental results suggested variational inference can indeed result in deep exploration.

In this paper, we propose an amortized variational Deep Q network framework, in which the state-action value function, i.e. Q, is considered to be a random variable. The posterior of Q is approximated by variational inference. Our contributions can be summarized as follows. First, the parameters of the auxiliary posterior distribution of Q are the output of a deep neural network as in the variational autoencoder (VAE)~\cite{kingma2014auto}. Second, to encourage exploration, we add entropy bonus to the DQN loss. The equivalence between the loss function of our new model and the amortized variational inference loss is then established. Third, we propose a two-stage (pre-train and fine-tune) training algorithm. The auxiliary posterior is learned from a family of heavy-tailed Cauchy distribution in the pre-train stage, and Gaussian in the fine-tune stage. Finally, experiments on chain MDP and OpenAI Gym tasks verify the proposed method performs significantly better than DQN, VDQN, and NoisyNet.

\section{Background}\label{Background}

\subsection{Markov Decision Process and Reinforcement learning}

Reinforcement Learning has been playing an important role in the thriving of artificial intelligence. It aims to find a policy for an agent so that it can perform optimally in the environment. RL can be modeled as a Markov Decision Process (MDP). Consider a finite-horizon MDP with state and action space defined by a tuple $(\mathcal{S},\mathcal{A},\mu_0,p,r,\pi,T)$ where $\mathcal{S}\in\mathbb{R}^D$ denotes the state space, $\mathcal{A}\in\mathbb{R}^d$ the action space, $\mu_0$ the initial distribution of the state, $r:\mathcal{S}\rightarrow\mathbb{R}$ the reward, and $T$ the time horizon, respectively. At each time $t$, there are $s_t\in\mathcal{S}$, $a_t\in\mathcal{A}$ and a transition probability $p:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\rightarrow \mathbb{R}$, where $p(s_{t+1}|a_t,s_t)$ denotes the transition probability of $s_{t+1}$ conditionally on $s_t$ and $a_t$. The policy $\pi:\mathcal{S}\times\mathcal{A}\times \{0,1,\cdots\,T\}\rightarrow\mathbb{R}$, where $\pi(a_t|s_t)$ is the probability of choosing action $a_t$ when observing current state $s_t$. Fig.~\ref{RL_introduction} shows the flowchart of the finite-horizon MDP.
\begin{figure}[htbp]
\centerline{\includegraphics[width=0.6\columnwidth]{RL_introduction.png}}
\caption{The flowchart of a finite-horizon Markov chain process.}\label{RL_introduction}
\end{figure}

At each state $s_t$, the agent takes an action $a_t$ according to present policy $\pi(a_t|s_t)$. Given this action, environment responds with a new state $s_{t+1}$ and an immediate reward $r_{t+1}$. The goal is to find a policy $\pi = p(a_t|s_t)$  so as to maximize the expectation of total rewards $\mathbb{E}(\sum_{t=1}^T \alpha_t r_t) $ where $\alpha_t$ denotes time-step dependent decaying factors. In practice, the decaying factor is set to be the exponential power of a constant, i.e. $\alpha_t=\gamma^t$.

There are many RL algorithms, such as Q-learning, sarsa, deep Q network and policy gradient, which are developed to deal with different environments~\cite{Sutton1998Reinforcement}. Among them, the Q-learning is developed for MDP with discrete state and action space. Its critical idea is to use the action value function $Q(s,a)$ to estimate the reward in case $S_t = s$ and $A_t = a$, by using Bellman equation:
\begin{equation}
Q(s,a) = \mathbb{E}\left(R_{t+1}+ \max_{a_{t+1}}\gamma Q(S_{t+1}, a_{t+1})|S_t = s, A_t = a\right)\label{Bellman equation}
\end{equation}

In practice, Monte Carlo sampling is used to approximate the expectation in Eq.~\ref{Bellman equation}. A number of trajectories $\tau^i\doteq\{s_0^i,r_0^i,a_0^i,\cdots,s_T^i,r_T^i\},\ i=1,\cdots,l$ are sampled and used to update $Q$ as follows:
\begin{equation}
Q(s_t^i,a_t^i)\leftarrow(1-\alpha) Q(s_t^i,a_t^i)+\alpha\left(r_{t+1}^i+\max_{a_{t+1}}\gamma Q(s_{t+1}^i, a_{t+1})\right)
\end{equation}

As $\mathcal{A}$ is discrete and finite, the optimal policy can be regarded as $\pi(a|s)=\arg\max_{a\in\mathcal{A}}Q(s,a)$.

In Deep Q Network (DQN) \cite{mnih2013playing},  the action value function is approximated by a deep neural network with parameter $\theta$ (denoted as $Q_{\theta}(s,a)$). Similar to the Q-learning, the loss function of DQN is
\begin{equation}\label{dqn_loss}
L(\theta)=\frac{1}{l}\frac{1}{T}\sum_{i=1}^l\sum_{t=1}^T\left[Q_{\theta}(s_t^i,a_t^i)-\left(r(s_{t+1}^i)+\max_{a_{t+1}\in {\cal A}}\gamma Q_{\theta}(s_{t+1}^i, a_{t+1})\right)\right]^2
\end{equation}

\subsection{Amortized Variational Inference}

Let $p(x,z)$ be a joint distribution for latent variable $z\in \mathcal{Z}$ and observed variable $x\in\mathcal{X}$. Inference aims to compute the posterior $p(z|x)$. However, the posterior is usually intractable to compute, thus approximation is inevitable. Variational inference uses a family of tractable distributions $\mathcal{Q}$ parameterized by $\psi$ over $z$ and finds a member $q_{\psi_x^*}\in\mathcal{Q}$ by minimizing the Kullback-Leibler (KL) divergence: $q_{\psi_x^*}=\arg\min_{q_{\psi}\in\mathcal{Q}}D_{KL}(q_{\psi}(z)||p(z|x))$ where $q_{\psi}(z)$ is called the auxiliary posterior distribution. For a dataset $\mathcal{D}$, we need to find the auxiliary posterior $q_{\psi_x^*}$ for each latent $z$ associated with each $x\in\mathcal{D}$ under i.i.d. assumpution.

Amortized variational inference (AVI) \cite{gershman2014amortized} assumes that the local variational parameter can be predicted by a parameterized function whose parameters are shared by all data points, Namely, to find a parameterized function $f_{\theta}:\mathcal{X}\rightarrow\mathcal{Q}$ to predict $\psi_x^*$ where $\theta$ is the parameter. Neural network can be readily used as the function $f_{\theta}$ due to its universal function approximation capability~\cite{cybenko1989approximation,hornik1989multilayer}. In literature, such a network is often called inference network. AVI has been applied in VAE~\cite{kingma2014auto}. In VAE, during training, the expectation over $q_{\psi_x}(z)$ is needed for the calculation of the variational lower bound:
\begin{equation}
\log p(x) \geq \mathbb{E}_{q_{\psi_x}(z)} [\log p(x|z)] - D_{KL}(q_{\psi_x}(z)||p(z|x))
\end{equation}
The stochastic back-propagation~\cite{rezende2015variational} is used to approximate the expectation by Monte Carlo gradient estimation, including a reparameterization step and a back-propagation with Monte Carlo step.
\section{Related Work}\label{Related Work}

NoisyNet~\cite{Fortunato2017Noisy} proposes to directly add noises to the parameters of the value/policy function. It is shown that doing so can enable consistent exploration. Bootstrapped DQN~\cite{osband2016deep} applies the bootstrap to approximate the posterior of the value function by applying different heads trained with different bootstrapped data, which entails diverse strategies to encourage exploration, but the improvement is limited and may be compute-intensive when training.

Bayes-by-Backprop Q-Network (BBQN)~\cite{blundell2015weight} randomizes the policy space to achieve the balance between exploration and exploitation. It achieves good performance on dialogue tasks when combined with Replay Buffer Spiking (RBS) trick and pipeline of natural language processing system~\cite{DBLP:conf/aaai/LiptonLG00D18}.

Variational Information Maximizing Exploration (VIME)~\cite{houthooft2016vime} proposes an intrinsic reward based on information bonus to encourage exploration, which is a curiosity-driven exploration strategy. It learns a dynamic model of environment represented by a Bayesian neural network and uses a variational family distribution to approximate posterior and compute bonus.

Variational Deep Q Network (VDQN)~\cite{tang2017variational} shares a similar spirit with BBQN. It proposes a surrogate objective to the Bellman error by adding an entropy term which can explicitly encourages exploration. It is shown that the surrogate objective is equivalent to the variational inference loss. Moreover, the algorithm is interpreted as performing approximate Thompson Sampling. Bayesian neural network is used as the posterior and variational inference subroutines are used to minimize the KL divergence.

Recently, distributional RL methods~\cite{tang2018exploration,moerland2018the,mavrin2019distributional,dabney2017distributional} have been proposed, in which the cumulative reward is considered as random variable and the distributional Bellman equation is applied to iteratively update the distribution to minimize the distance between the target distribution and the predicted distribution.  It is more challenging than computing the expectation of the reward but it can offer more information.
\section{Method}\label{Method}

\subsection{Model}\label{avdqn}

In DQN, $\theta$ in Eq.~\ref{dqn_loss} denotes the parameters of the DNN. In VDQN, $\theta$ is regarded as a random variable. The posterior of $\theta$ of a Bayesian neural network is approximated by variational inference.

In our work, we consider $Q_{\theta}(a_t^i,s_t^i)$ to be a random variable and borrow the idea of amortized variational inference for the updating of $q(Q_{\theta}(a_t^i,s_t^i))$. Assume $Q_{\theta}(a_t^i,s_t^i) \sim q(Q_{\theta}(a_t^i,s_t^i);\alpha_{\theta}(a_t^i,s_t^i))$ where $\alpha_{\theta}(a_t^i,s_t^i)$ is the parameter of the posterior distribution $q$. The same as in VAE, $\alpha_{\theta}(a_t^i,s_t^i)$ is the output of a deep neural network (DNN) with parameter $\theta$ and input $a_t^i,s_t^i$. Since $Q_{\theta}(a_t^i,s_t^i)$ is a random variable, Eq.~\ref{dqn_loss} can be re-written as follows:
\begin{eqnarray}
\min_{\theta}\sum_{i,t} \mathbb{E}_{Q_{\theta}(a_t^i,s_t^i)\sim q(Q_{\theta}(a_t^i,s_t^i);\alpha_{\theta}(a_t^i,s_t^i))}\left[\left (Q_{\theta}(a_t^i,s_t^i)-\gamma \cdot
\max_{a_{t+1}\in {\cal A}}Q_{\theta}(a_{t+1},s_{t+1}^i)-r(s_{t+1}^i)\right)^2\right]\label{AVDQN_log}
\end{eqnarray}

If we regard $\gamma \max_{a_{t+1}\in {\cal A}}Q_{\theta}(a_{t+1},s_{t+1}^i)+r(s_{t+1}^i)$ as observed data (denoted as $\mathcal{D}_t^i$), and assume the prior $p(Q_{\theta}(a_t^i,s_t^i))\propto1$, Eq.~\ref{AVDQN_log} is equal to
\begin{eqnarray}
\max_{\theta}\sum_{i,t} \mathbb{E}_{Q_{\theta}(a_t^i,s_t^i)\sim q(Q_{\theta}(a_t^i,s_t^i);\alpha_{\theta}(a_t^i,s_t^i))}\left[ \ln P\left(\mathcal{D}_t^i|Q_{\theta}\left(a_t^i,s_t^i\right)\right)\right]\label{log_posterior}
\end{eqnarray}
where $p\left(\mathcal{D}_t^i | Q_{\theta}(a_t^i,s_t^i)\right)=\mathcal{N}\left(\mathcal{D}_t^i|Q_{\theta}(a_t^i,s_t^i),1\right)$ with $Q_{\theta}(a_t^i,s_t^i)$ as the mean and unit variance.

The same as in \cite{tang2017variational}, adding an entropy term over $q(Q_{\theta}(a_t^i,s_t^i);\alpha_{\theta}(a_t^i,s_t^i))$ on Eq.~\ref{log_posterior} to encourage exploration, we obtain the following loss function
\begin{eqnarray}
\max_{\theta}\sum_{i,t} \mathbb{E}_{Q_{\theta}(a_t^i,s_t^i)\sim q(Q_{\theta}(a_t^i,s_t^i);\alpha_{\theta}(a_t^i,s_t^i))}\left[ \ln p\left(Q_{\theta}\left(a_t^i,s_t^i\right)|\mathcal{D}_t^i\right)\right]+\mathcal{H}[q(Q_{\theta}(a_t^i,s_t^i);\alpha_{\theta}(a_t^i,s_t^i))]\label{AVDQN_LOSS}
\end{eqnarray}where $\mathcal{H}[\cdot]$ denotes the entropy. Eq.~\ref{AVDQN_LOSS} can then be rewritten as
\begin{eqnarray}
\min_{\theta}\sum_{i,t} D_{\text{KL}}\left(q(Q_{\theta}(a_t^i,s_t^i);\alpha_{\theta}(a_t^i,s_t^i))||P(Q_{\theta}(a_t^i,s_t^i)|\mathcal{D}_t^i)\right)\label{KL_divergence}
\end{eqnarray}Eq.~\ref{KL_divergence} indicates that Eq.~\ref{AVDQN_LOSS} actually finds $\theta$ to minimize the sum of all the KL-divergences between $q(Q_{\theta}(a_t^i,s_t^i))$ and the posterior $p\left(Q_{\theta}(a_t^i,s_t^i)|\mathcal{D}_t^i\right)$ at the $i$th episode. That is, if we optimize Eq.~\ref{AVDQN_LOSS}, an optimal $\theta^*$ can be obtained to make $q(Q_{\theta^*}(a_t^i,s_t^i);\alpha_{\theta^*}(a_t^i,s_t^i))$ approximate the true posterior.
\begin{figure}[htbp]
\setlength{\belowcaptionskip}{-2pt}
\setlength{\abovecaptionskip}{-4pt}
\centering\includegraphics[width=0.82\columnwidth]{avdqn1}
\caption{Framework of VDQN and AVDQN. Blue (resp. red) arrow indicates the forward propagation (resp. backpropagation) in the neural network pipeline. Green arrow is the reparameterization trick.}\label{framework}
\end{figure}

Fig.~\ref{framework} shows the framework of VDQN and the proposed amortized variational deep Q network (AVDQN). In VDQN, the weights of the Bayesian neural network are random variables. In the forward propagation, a sampled neural network with parameter $\theta'$ is obtained by sampling the weights of the Bayesian network i.i.d. The action-value $Q_{\theta'}(s,a)$ is then obtained by taking $a$ and $s$ as input to the sampled neural network. In AVDQN, the weights $\theta$ of the network are deterministic values, $Q_{\theta}(s,a)$ is considered to be random. The parameters of the posterior of $Q_{\theta}(s,a)$ (in the figure, $\mu$ and $\sigma$ represent the mean and standard deviation, respectively) are the output of the network, while $Q_{\theta}(a,s)$ is the sampled value from ${\cal N}(\mu,\sigma^2)$. In the backpropagation, reparametrization trick is applied differently. In VDQN, it is applied to obtain a set of Bayesian networks, while in AVDQN, it is to obtain different parameters ($\mu$ and $\sigma$) of the auxiliary posterior distribution.

\subsection{Algorithm}\label{tavdqn}

Applying variational inference aims to enhance the exploration ability. However, balancing the exploration and exploitation is also important. Therefore, the training of AVDQN includes a pre-train stage and a fine-tune stage. In the pre-train stage, the posterior $q$ is assumed to be a Cauchy distribution. Cauchy distribution is heavy-tailed which means it is suitable for exploration. In the fine-tune stage, $q$ is considered to be a Gaussian distribution, which is preferable for exploitation.

The training procedure of AVDQN is summarized in Alg.~\ref{alg:tavdqn}. The same as in DQN~\cite{mnih2013playing}, to stabilize training, we use an evaluation network and a target network with parameter $\theta$ and $\theta^-$, respectively. From line~\ref{a6} to \ref{a11}, we sample a $Q_\theta(s_t,a) $ from distribution $q$ for each $a\in {\cal A}$. From line~\ref{a12} to~\ref{a13}, at each time step $t$, we select action by being greedy w.r.t. $Q_\theta(s_t,a)$ and add experience tuple $ \{ s_t,a_t,r_t,s_{t+1} \} $ to buffer $R$. From line~\ref{a14} to \ref{a21}, when updating parameters, we sample a mini-batch of $M$ tuples from the buffer $R$ and use the target network to obtain the target value. From line~\ref{a22} to~\ref{a23}, the gradient of the proposed loss function is computed and used to update parameter $\theta$. In line~\ref{a25}, the target network parameter $\theta^-$ is updated at every $\tau$ steps. In the first $\omega$ episodes, the Cauchy distribution is applied when sampling Q (line~\ref{sampleC1} and~\ref{sampleC2}) where ${\cal C}(\cdot,\cdot)$ denotes the Cauchy distribution. The Gaussian distribution is used in the fine tune stage (line~\ref{a11} and~\ref{sampleN}).
\begin{algorithm}
\caption{Amortized Variational DQN}\label{alg:tavdqn}
\vspace{-1pt}
\LinesNumbered
\KwIn{the target network update period $\tau$; the learning rate $\alpha$; the number of episodes used for pre-train $\omega$; and the discount factor $\gamma$;}
Initialize: parameter $\theta$ and $\theta ^ -$, replay buffer $R \gets \{ \}$, step counter $counter \gets 0$ and $e\gets 0$

\While{not converge}
{
 $e\gets e+1, t\gets 0$;\\
 \While{\rm{episode not terminated}}
 {
  $counter \gets counter+1$

  \If {$e \le \omega $}
  {
   Compute $\mu_{t}^i,\delta_{t}^i = \text{DNN}(s_t, a_i;\theta)$ for every $a_i \in {\cal A}$; \\\label{a6}
   Sample $Q^i_\theta(s_t,a_i) $ from ${\cal C}\left(\mu_t^i,\delta_t^i\right)$, for every $a_i\in\mathcal{A}$ \label{sampleC1}
  }
  \Else
  {
   Compute $\mu_{t}^i,\sigma_t^i = \text{DNN}(s_t, a_i;\theta)$ for every $a_i \in {\cal A}$; \\
   Sample $Q_\theta(s_t,a_i) $ from ${\cal N}\left(\mu_t,[\sigma_t^i]^2\right)$, for every $a_i\in\mathcal{A}$\label{a11}
   }
  For state $s_t$, set $a_t=\arg\max_{a_i \in {\cal A}} Q^i_\theta(s_t,a_i)$, obtain transition $s_{t+1}$ and reward $r_t$\label{a12}

  Save experience tuple $ \{ s_t,a_t,r_t,s_{t+1} \} $ to buffer $R$\label{a13}

  Sample $M$ tuples $D= \{ s_j,a_j,r_j,s'_j \}$ from $R$\label{a14}

  \If {$e\le\omega$}
  {
    \For{$j = 1\to M$}
  {
   Compute $\mu_{t}^i,\delta^i_{t} = \text{DNN}(s'_j, a_i;\theta^-)$ for all $a_i \in {\cal A}$;\\
   Sample $Q_{\theta^-}(s'_j,a_i)$ from ${\cal C}(\mu^i_{t},\delta^i_{t})$ for $a_i\in\mathcal{A}$\label{sampleC2}
   }
  }
  \Else
  {
  \For{$j = 1\to M$}
  {
   Compute $\mu_{t}^i,{\sigma_t^i} = \text{DNN}(s'_j, a_i; \theta^-)$ for all $a_i \in {\cal A}$; \\
   Sample $Q_{\theta^-}(s'_j,a_i)$ from ${\cal N}\left(\mu_{t}^i,[{\sigma_t^i}]^2\right)$ for all $a_i\in\mathcal{A}$ \label{sampleN}
   }
  }

  Compute the target value for tuples in $D$: $d_j=r_j+\gamma \cdot \max_{a_i}Q_{\theta ^-} (s'_j,a_i), 1\leq j \leq N$\label{a21}

  Take the gradient of Eq.~\ref{AVDQN_LOSS} to obtain $\triangle \theta$\label{a22}

  $\theta = \theta - \alpha \triangle \theta$\label{a23}

  $t \gets t+1$

  \If {$counter$ {\rm mod} $\tau=0$}
  {
   Update the parameters of the target network $\theta ^- \leftarrow \theta$\label{a25}
  }
 }
}

\end{algorithm}
The original DQN samples experiences uniformly from the buffer. In sparse reward scenarios, this may lead to sampling some important but rare experiences with very small probability and results in converging to local optimum~\cite{schaul2015prioritized}. Prioritized Experience Replay (PER)~\cite{schaul2015prioritized} proposes to give each experience a priority. When sampling, experience with higher priority is preferable. The priority is defined either based on proportion or rank. For the proportional priority, the priority of the $i$th experience is $p_i=(\delta_i+\varepsilon)$ where $\delta_i$ is the temporal-difference error. For the rank-based priority, $p_i=\frac{1}{\text{rank}_i}$ where $\text{rank}_i$ is the rank of the transition $i$ when the replay memory is sorted according to $|\delta_i|$. To be more specific, the probability of the $i$th experience to be sampled is $p(i)=\frac{p_i^{\alpha}}{\sum_k {p_k}^\alpha}$. It degenerates to uniform sampling when $\alpha=0$. In our study, the rank-based prioritized replay is used since it is insensitive to outliers and blind to the relative error scales~\cite{schaul2015prioritized} .

\subsection{The Number of Model Parameters}\label{ca}

Assuming the deep neural network has $\ell$ hidden layers, and the $i$th layer has $H_i$ neurons in DQN, VDQN, NoisyNet, and AVDQN. The dimension of the input layer is $I$ and the dimension of the output layer of AVDQN and VDQN is $n_A\cdot n_V$ and $n_V$, respectively, where $n_A$ is the number of parameters of $q(Q)$. In total, the number of parameters of AVDQN is $(I+1)H_1+\sum_{i=1}^{\ell-1}(H_{i}+1)H_{i+1}+n_A(H_{\ell}+1)n_V$. In AVDQN, $n_A=2$ because $q(Q)$ follows Cauchy or Gaussian distribution. In VDQN, each weight has an auxiliary posterior distribution, thus the total number of parameters is $((I+1)H_1+\sum_{i=1}^{\ell-1}(H_{i}+1)H_{i+1}+(H_l+1)n_V)\times n'_A$, where $n'_A=2$ since each weight follows a Gaussian distribution.

In comparison, we see that VDQN has $((I+1)H_1+\sum_{i=1}^{l-1}(H_{i}+1)H_{i+1})\times (n'_A-1)$ more parameters than AVDQN. In case all $H_i$'s are the same ($=H$), VDQN has $(H(\ell-1)+ I + \ell)H$ more parameters than AVDQN.

\section{Experiments}\label{Experiments}

In this section, we test AVDQN on four classic control tasks in OpenAI Gym: CartPole-v0, CartPole-v1, Acrobot-v1, and MountainCar-v0 and four chain MDP tasks with $N=5, 10, 50$ and $100$ (for details about the task please see Appendix). These tasks are very challenging because only through extensive exploration, a proper policy can be learned to control mechanical systems. Particularly, for the CartPole-v0/CartPole-v1 tasks, it is impossible for $\varepsilon$-greedy strategy to learn a policy to balance the pole. On the other hand, learning an optimal policy for chain MDP tasks  becomes much more difficult as $N$ increases~\cite{tang2017variational,osband2016deep}. DQN, VDQN, and NoisyNet are used as baselines in the experiments.

\subsection{Implementation Details}\label{det}

In all the experiments, we use the fully connected network with two hidden layers and ReLU activation as the inference network. Each hidden layer has 100 units. The size of the mini-batch sampled from the buffer $R$ is $128$. The target network is updated every $100$ time steps. All results are averaged over $5$ different random initializations. The maximum size of the replay buffer is $10^6$. DQN uses $\varepsilon$-greedy strategy where $\varepsilon$ linearly decays from $1$ to $0.01$ as carried out in the original reference.

On classic control tasks, we use $1500$ episodes to train the model. The learning rate $\alpha$ for DQN, VDQN, NoisyNet, and AVDQN is $10^{-2},10^{-3},10^{-3},10^{-3}$, respectively. The discount factor $\gamma=0.99$ for all the compared algorithms. For the chain MDP tasks, we use $1000$ episodes to train $5,10$ and $3000$ episodes for $N=50,100$. For AVDQN, the last $200$ episodes are used for fine-tune. The learning rate is set constant in the pre-train stage, and set as $\alpha/(1+0.9 \times (e-\omega))$ in the fine-tune stage where $e$ is the training episode, and $\omega$ is the number of episodes used for the pre-train. The discount factor is $\gamma=1$. For NoisyNet, we use independent Gaussian noise. In AVDQN, the rank-based prioritized replay is applied~\cite{schaul2015prioritized}. We save the experience to a priority queue implemented with an array-based binary heap. The heap array was infrequently sorted every $1000$ time steps to prevent the heap becoming too unbalanced.

\subsection{Results}\label{result}

First we compare the performances of DQN, VDQN, NoisyNet, and AVDQN based on the same running time. Figs.~\ref{fig4} and~\ref{fig5} show the training curves of these compared algorithms. We can see that, AVDQN can reach the highest final reward on all the tasks. Particularly, we find in the chain MDP tasks with $N=5$ and $10$, VDQN is unsteady: there is an obvious declining from $50$ seconds upwards. Possible reason is that VDQN over-emphasizes on exploration so that VDQN cannot converge on some simple tasks. In comparison, on difficult tasks (chain MDP with $N=50$ and $100$) DQN cannot explore deeply while VDQN performs well, but AVDQN performs the best. In addition, NoisyNet reveals unsteady exploration in chain MDP problem when $N$ is large. We may thus conclude that AVDQN can successfully balance the exploration and exploitation.

\begin{figure}[htbp]
\subfigure[CartPole-v0]{\includegraphics[scale = 0.4]{3-a-new.png}}\qquad%
\subfigure[CartPole-v1]{\includegraphics[scale=0.4]{3-b-new.png}}\\
\subfigure[Acrobot-v1]{\includegraphics[scale = 0.4]{3-c-new.png}}
\subfigure[MountainCar-v0]{\includegraphics[scale = 0.4]{3-d-new.png}}
\caption{The training curves of DQN, VDQN, NoisyNet, and AVDQN on the four classic tasks within the same running time (in seconds).}\label{fig4}
\end{figure}

\begin{figure}
\subfigure[$N=5$]{\includegraphics[scale=0.4]{4-a-new.png}}\qquad
\subfigure[$N=10$]{\includegraphics[scale=0.4]{4-b-new.png}}\\
\subfigure[$N=50$]{\includegraphics[scale=0.4]{4-c-new.png}}\qquad
\subfigure[$N=100$]{\includegraphics[scale=0.4]{4-d-new.png}}
\caption{The training curves of DQN, VDQN, NoisyNet, and AVDQN on the four chain MDP tasks averaged over 5 random initializations with same running time.}\label{fig5}
\end{figure}

Table~\ref{tab1} shows the average running time of the compared algorithms within the same number of episodes. It can be observed that AVDQN requires much less time compared with VDQN and NoisyNet. Specifically, training AVDQN is at least two times and at most eight times faster than training VDQN and NoisyNet. Table~\ref{para} shows the number of parameters of DQN, VDQN, NoisyNet, and AVDQN. It is seen that the number of parameters of DQN and AVDQN are similar, and are only half to that of VDQN and NoisyNet. Table~\ref{tab2} shows the final rewards obtained by the compared algorithms within the same running time. From the table we can see that in all the chain MDP tasks, AVDQN has reached the global optimal reward (11) and on the classical control tasks in OpenAI Gym, AVDQN has achieved higher rewards than DQN, VDQN and NoisyNet.
\begin{table}[htbp]
\centering
\caption{The average running time of DQN, VDQN, NoisyNet, and AVDQN on the considered tasks with the same number of episodes (in seconds)}.
    \label{tab1}
    \begin{tabular}{lllllll}
    \toprule
        Tasks & Episodes & DQN & VDQN & NoisyNet & AVDQN \\
    \midrule
        CartPole-v0 & 1500 & 1176 & 7531 & 7955 & 1286 \\
        CartPole-v1 & 1500 & 2350 & 16461 & 16871 & 2196 \\
        Acrobot-v1 & 1500 & 884 & 4241 & 3882 & 1223 \\
        MountainCar-v0 & 1500 & 1241 & 6610 & 5428 & 1790 \\
        MDP $N=5$ & 1000 & 55 & 379 & 389 & 155 \\
        MDP $N=10$ & 1000 & 77 & 544 & 555 & 188 \\
        MDP $N=50$ & 3000 & 796 & 6722 & 7070 & 1448 \\
        MDP $N=100$ & 3000 & 1584 & 16383 & 17193 & 2866 \\
    \bottomrule
    \end{tabular}
\end{table}

\begin{table}[htbp]
\centering
\caption{The number of parameters of DQN, AVDQN, VDQN, and NoisyNet on the considered tasks.}\label{para}
    \begin{tabular}{llllll}
    \toprule
        Tasks & DQN & AVDQN & VDQN & NoisyNet \\
    \midrule
        CartPole-v0 & 10802 & 11004 & 21604 & 21604 \\
        CartPole-v1 & 10802 & 11004 & 21604 & 21604 \\
        Acrobot-v1 & 11103 & 11406 & 22206 & 22206 \\
        MountainCar-v0 & 10703 & 11006 & 21406 & 21406 \\
        MDP $N=5$ & 10902 & 11104 & 21804 & 21804 \\
        MDP $N=10$ & 11402 & 11604 & 22804 & 22804 \\
        MDP $N=50$ & 15402 & 15604 & 30804 & 30804 \\
        MDP $N=100$ & 20402 & 20604 & 40804 & 40804 \\
    \bottomrule
    \end{tabular}
\end{table}

\begin{table}[H]
    \centering
   \caption{The average final rewards in the last 10 episodes of DQN, VDQN, NoisyNet, and AVDQN on the considered tasks under the same running time. }

    \label{tab2}
    \begin{tabular}{llllll}
    \toprule
        Tasks & Time(s) & DQN & VDQN & NoisyNet & AVDQN \\
    \midrule
        CartPole-v0 & 1300 & 198.89 & 196.58 & 185.64 & \textbf{200} \\
        CartPole-v1 & 2350 & 416.71 & 480.33 & 362.64 & \textbf{500} \\
        Acrobot-v1 & 1250 & -104.69 & -104.76 & -92.18 & \textbf{-74.51} \\
        MountainCar-v0 & 1800 & -178.27 & -176.84 & -119.58 & \textbf{-107.09} \\
        MDP $N=5$ & 155 & 9 & 8.62 & 10 & \textbf{11} \\
        MDP $N=10$ & 190 & 8.69 & 6.64 & 8.56 & \textbf{11} \\
        MDP $N=50$ & 1450 & 0.36 & 9.73 & 1.89 & \textbf{11} \\
        MDP $N=100$ & 3000 & 0.04 & 10.67 & 0.24 & \textbf{11} \\
    \bottomrule
    \end{tabular}
\end{table}

\section{Conclusion}\label{Conclusion}

We proposed a framework called Amortized Variational Deep Q Network (AVDQN) to tackle the deep exploration problem in reinforcement learning. In AVDQN, the action value function is considered to be random variables. Under amortized variational inference, an inference network is used to output the parameters of the posterior distribution of $Q$. To train the model, we used the heavy-tailed Cauchy distribution in the pre-train stage to explore the parameter space, and used the Gaussian distribution in the fine-tune stage to exploit for the optimal parameters of the inference network. The experiments on classic control tasks and chain MDP tasks showed that AVDQN performs significantly better than the state-of-the-art RL methods, including DQN, VDQN, and NoisyNet, in terms of the final reward and the training time.
\section*{Acknowledgements}
This work was partly supported by the National Natural Science Foundation of China (grant no.
11991023, 62076197), the Major Project of National Science Foundation of China (grant no. U1811461), and Key
Project of National Science Foundation of China (grant no. 11690011).
\bibliographystyle{unsrt}


\section{Introduction}
	\label{sect:1}
	Over the last several decades, numerous studies have been conducted on the semilinear elliptic equations
	\begin{align}
	\label{mainpro}
	-\Delta u(x)=f(u(x)), ~~x \in \Omega,
	\end{align}
	with appropriate boundary conditions. One such example is the Dirichlet problem
	\begin{align}
	\label{mainpro:d}
	\left\{\begin{array}{l l}
	-\Delta u(x)=f(u(x)), &x \in \Omega,\\
	u(x)=0,  &x \in \partial\Omega,\\
	\end{array}\right.
	\end{align}
	where $\Omega \subset \mathbb{R}^{N}$~$(N=1,2,3,\cdots)$ is a bounded domain,
	$\Delta $ is the Laplacian,
	and $ f: \mathbb{R} \rightarrow \mathbb{R}$ is a given nonlinear map.
	In particular, the investigation of positive solutions of \eqref{mainpro} has attracted significant attention \cite{lions1982existence,gidas1979symmetry,lin1994uniqueness,damascelli1999qualitative,gladiali2011bifurcation,de2019morse}.
	Here are some important examples derived from model problems for many applications in which we are interested.
	Positive solutions of problem \eqref{mainpro:d} with $ f(t) = \lambda t + t|t|^{p-1} $, $ \lambda \in [0,\lambda_1(\Omega)) $, $p\in (1,p^*)$ have been investigated from various points of view --- uniqueness, multiplicity, nondegeneracy, symmetricity, and so on \cite{gidas1979symmetry,lin1994uniqueness,damascelli1999qualitative,gladiali2011bifurcation,de2019morse},
	where
	$p^*=\infty$ when $N=1,2$ and $p^*=(N+2)/(N-2)$ when $ N\geq3 $, and
	$ \lambda_1(\Omega) $ is the first eigenvalue of $ -\Delta $ imposed on the homogeneous Dirichlet boundary value condition;
	the eigenvalue problem is understood in the weak sense at least when $ \Omega $ is not regular.
	Another important nonlinearity is $ f(t)=\varepsilon^{-2}(t-t^3) $, which corresponds to the stationary problem of the Allen-Cahn equation motivated by \cite{allen1979microscopic} and has been investigated by many researchers. The value $\varepsilon>0$ is a small parameter related to the so-called singular perturbation.
	A variational method ensures that problem \eqref{mainpro:d} with this nonlinearity has a positive solution when $ \varepsilon^{-2} \geq \lambda_1(\Omega)$.
	On the other hand, when $ \varepsilon^{-2} < \lambda_1(\Omega)$, no positive solution is admitted as we prove later.
	Despite these results, quantitative information about the positive solutions, such as their shape, has not been clarified analytically.
	Throughout this paper, $ H^k(\Omega) $ denotes the $k$-th order $ L^2 $ Sobolev space.
	We define $ H^1_0(\Omega) := \{ u \in H^1(\Omega) : u = 0~\mbox{on} ~ \partial \Omega\} $, with the inner product $(u, v)_{H^1_0}:=(\nabla u, \nabla v)_{L^2}$ and the norm $\| u \|_{H^1_0} := \sqrt{(u, u)_{\smash{H^1_0}}}$.
	We say that the solution $ u $ is ``positive'' if $ u>0 $ in $ \Omega $, and ``nonnegative'' if $ u\geq 0 $ in $ \Omega $.

	This paper is concerned with numerical verification (also known as verified numerical computation or computer-assisted proof) for positive weak solutions of problem \eqref{mainpro:d}.
 	The target weak form of \eqref{mainpro:d} will be shown explicitly at the beginning of the next section together with further regularity assumptions for nonlinearity $ f $.
	The pioneering research on numerical verification methods for partial differential equations began with \cite{nakao1988numerical,plum1991computer} and has been further developed by many researchers (see, for example, \cite{plum2008,nakao2011numerical,tanaka2017sharp} and the references therein).
	In particular, studies have revealed that methods based on several fixed point theorems for Newton operators (including their simplified versions) are greatly effective.
	This approach is closely related to our method (see Section \ref{sec:nk}).
	These methods enable us to obtain an explicit ball containing exact solutions of \eqref{mainpro:d}.
	For weak solutions, this is typically done in the sense of the norm $\left\|\cdot\right\|_{H_{0}^{1}}$, because $ H^1_0(\Omega)$ is a natural solution-space for \eqref{mainpro:d} in the distributional sense.
	In other words, they allow us to prove, for a numerical approximation $ \hat{u} \in H^1_0(\Omega) $, the existence of an exact weak solution $ u \in H^1_0(\Omega) $ satisfying
	\begin{align}
	\label{h10error}
	\left\|u-\hat{u}\right\|_{H_{0}^{1}} \leq \rho
	\end{align}
	for an explicit error bound $ \rho $.
	Therefore, these methods have the advantage that quantitative information about solutions to a target equation is provided accurately in a strict mathematical sense.
	However, irrespective of how small the error bound is, the positivity of some solutions is not guaranteed without additional considerations.
	In particular, in the homogeneous Dirichlet case \eqref{mainpro:d}, it is possible for a solution that is verified by such methods to be negative near the boundary $\partial\Omega$.
	When $ N \leq 3 $, positivity can be verified if we have an $ H^3 $-norm evaluation $ \|u-\hat{u}\|_{H^3} $ with $\hat{u} \in H^3(\Omega)$ and an explicit bound for the embedding $ H^3(\Omega) \cap H^1_0(\Omega) \hookrightarrow C^1(\overline{\Omega})$.
	The bound can be evaluated by applying the bound for the embedding $ H^2(\Omega) \hookrightarrow C(\overline{\Omega})$ to $ u $ and its derivatives $\partial u/ \partial x_i$ ($i=1,2,\cdots,N$).
	However, for some shapes of $ \Omega $, for example, a nonconvex polygonal domain, the regularity of the solution $ u $ is in general outside $ H^3(\Omega) $.
	Otherwise, even if $ u, \hat{u} \in  H^3(\Omega) $, evaluating $ \|u-\hat{u}\|_{H^3} $ itself is not easy, and troublesome numerical techniques for estimating the slope of $ \hat{u} $ are required to complete the proof of positivity.

	In previous studies, we developed the method of verifying the positivity of solutions of \eqref{mainpro:d} \cite{tanaka2015numerical,tanaka2017numerical,tanaka2017sharp}.
	These methods succeeded in verifying the existence of positive solutions by checking simple conditions, but they required $ L^{\infty} $-error estimation obtained by considering the embedding from a solution set with $ H^2 $-regularity.
	This is primarily for these methods because they need to find a subdomain where $ u $ may be negative and to evaluate the first eigenvalue of $ -\Delta $ on this subdomain.
	For the same reasons as mentioned above, this requirement narrows the range of applications of these methods.

	The main contribution of this paper is proposing methods for verifying the positivity of solutions $ u $ of \eqref{mainpro:d} while assuming only $ H^1_0 $-error estimation \eqref{h10error}, which is more generally applicable (wider class of domains, solutions that lack $ H^2 $-regularity, and so on) than previous methods \cite{tanaka2015numerical,tanaka2017numerical,tanaka2017sharp}.
	Theorems \ref{theo1} and \ref{theo:keeppositive}, and Corollary \ref{theo:linf} enable us to verify the nonnegativity of $ u $ under suitable conditions.
	The positivity of $ u $ follows from its nonnegativity using a maximum principle under appropriate conditions, such as when $ f $ is a subcritical polynomial (see, for example, \cite{drabek2009maximum}).
	Table \ref{table:applicable} summarizes the scope of application of our theorems to the case in which $ f $ is a subcritical polynomial $ f(t) = \lambda t + \sum_{i=2}^{n(<p^*)}a_{i}t|t|^{i-1}$, $ \lambda, a_i \in \mathbb{R}$, $ a_i \neq 0 $ for some $ i $.
	They are also applicable for more general nonlinearities other than polynomials (see Theorems \ref{theo1} and \ref{theo:keeppositive}, and Corollary \ref{theo:linf}).
	The table shows that the coefficient $ \lambda $ of the linear term essentially affects the existence of a positive solution of \eqref{mainpro:d} and the applicability of our theorems.
	In particular, there exists no positive solution in two specific cases.
	The first case is when $ \lambda \geq \lambda_1(\Omega)$ and $a_i \geq 0$ for all $ i $.
	This can be checked by multiplying \eqref{mainpro} with the first eigenfunction of $ -\Delta $ and integrating both sides.
	The second case is when $ \lambda < \lambda_1(\Omega)$ and $a_i \leq 0$ for all $ i $.
	This can be checked in the same way.
	Our theorems are applicable to the other cases where the problem \eqref{mainpro:d} may admit positive solutions (see again Table \ref{table:applicable}).

	\renewcommand{\arraystretch}{1.3}
	\begin{table}[h]
		\label{table:applicable}
		\caption{The applicability of our theorems to a subcritical polynomial $ f(t) = \lambda t + \sum_{i=2}^{n(<p^*)}a_{i}t|t|^{i-1}$, $ \lambda, a_i \in \mathbb{R} $, $ a_i \neq 0 $ for some $ i $. Here $ \lambda_1(\Omega) $ is the first eigenvalue of $ -\Delta $ imposed on the homogeneous Dirichlet boundary value condition.}
		\begin{center}
			\begin{tabular}{l|p{30mm}|p{30mm}}
				\hline
				\backslashbox[30mm]{$a_i$}{$\lambda$}&$\geq \lambda_1(\Omega)$&$<\lambda_1(\Omega)$\\
				\hline
				\hline
				$a_i \geq 0$ for all $ i $ &No positive solution& Theorem \ref{theo1}\\
				\hline
				$a_i \leq 0$ for all $ i $ &Theorem \ref{theo:keeppositive}& No positive solution\\
				\hline
				$a_i a_j < 0$ for some $ i,j $ &Corollary \ref{theo:linf}& Theorem \ref{theo1}\\
				\hline
			\end{tabular}
		\end{center}
	\end{table}
	\renewcommand{\arraystretch}{1}

	Assuming a certain growth condition for $ f $, Theorem \ref{theo1} provides a sufficient condition on the nonnegativity of $ u $ that can be checked only from $ H^1_0 $-evaluation as in \eqref{h10error}.
	Theorem \ref{theo1} is proved by a constructive norm estimation for the negative part $ u_{-}:= \max\left\{-u,0\right\} $ of $ u $ by considering the embedding corresponding to each exponent.
	Theorem \ref{theo:keeppositive} verifies the nonnegativity of $ u $ by imposing another inequality condition on $ f $,
	which is proved using a fundamentally different approach based on the Newton iteration.
	It ensures, under an appropriate condition, that a Newton sequence staring from a nonnegative function $ \omega \not\equiv
	0 $ remains nonnegative at every step, and as a result confirms the nonnegativity of the solution $ u $ to which it converges.
	Known as Newton-Kantorovich's theorem \cite{deuflhard1979affine, kantorovich1982functional},
	the convergence property of Newton's method on function spaces has been investigated in several studies (see, for example, \cite{argyros2008convergence, proinov2010new, argyros2012weaker} for recent results).
	However, little is known about the influence of Newton operators on the sign of functions.
	In this sense, Theorem \ref{theo:keeppositive} increases our understanding of Newton's method for elliptic equations.
	Among the cases listed in Table \ref{table:applicable}, when $ \lambda \geq \lambda_1(\Omega)$ and $a_i a_j < 0$ for some pair $ \{i,j\} $, only Corollary \ref{theo:linf} is applicable.
	In fact, this can be applied to all the cases listed in Table \ref{table:applicable}.
	Although Corollary \ref{theo:linf} is a generalized version of a previous theorem \cite[Theorem 2.2]{tanaka2017numerical}, this still requires $ L^{\infty} $-norm estimation; in this sense, a problem still remains.
	However, Theorems \ref{theo1} and \ref{theo:keeppositive} have wide application including several important problems such as those introduced at the beginning of this section.

	The remainder of this paper is organized as follows.
	In Section \ref{sec:1}, we provide a method for proving the positivity of a solution $ u $ of \eqref{mainpro:d} whose existence is confirmed as in inequality \eqref{h10error}.
	This method does not restrict verification methods, but admits any methods that can prove the existence of a solution with $ H^1_0 $-error estimation.
	Section \ref{sec:nk} provides another positivity-validation method based on the Newton iteration that retains nonnegativity.
	Finally, in Section \ref{sec:ex}, we present numerical examples where our method is applied to elliptic problems with the nonlinearities introduced at the beginning of this section.

	\section{Verification of positivity --- Constructive norm estimation for $ u_{-} $ when $\lambda < \lambda_1(\Omega)$}\label{sec:1}

	We begin by introducing some required notation.
	We denote $ V=H^1_0(\Omega) $ and $ V^*=$ (the topological dual of $ V $).
	For two Banach spaces $X$ and $Y$, the set of bounded linear operators from $X$ to $Y$ is denoted by ${\mathcal L}(X, Y)$ with the usual supremum norm $\| T \|_{{\mathcal L}(X, Y)} := \sup\{ \| T u \|_{Y} / \| u \|_{X} : {u \in X \setminus \{0\}}\}  $ for $T \in {\mathcal L}(X, Y)$.
	The norm bound for the embedding $V\hookrightarrow L^{p+1}\left(\Omega\right)$ is denoted by $C_{p+1}$, that is, $C_{p+1}$ is a positive number that satisfies
	\begin{align}
	\label{embedding}
	\left\|u\right\|_{L^{p+1}(\Omega)}\leq C_{p+1}\left\|u\right\|_{V}~~~{\rm for~all}~u\in V,
	\end{align}
	where $p\in [1,\infty)$ when $N=1,2$ and $p\in [1,p^*]$ when $ N\geq3 $.
	Assuming that $ f $ is a $ C^1 $ function satisfying
	\begin{align*}
		&|f(t)| \leq a_0 |t|^p + b_0 \text{~~~for~all~~} t \in \mathbb{R}\\
		&|f'(t)| \leq a_1 |t|^{p-1} + b_1 \text{~~~for~all~~} t \in \mathbb{R}
	\end{align*}
	for some $ a_0,a_1,b_0,b_1\geq 0 $ and $ p<p^* $,
	we define the operator $ F $ by
	\begin{align*}
		F : \left\{\begin{array}{ccc}{u(\cdot)} & {\mapsto} & {f(u(\cdot))}, \\
		{V} & {\rightarrow} & {V^{*}}.\end{array}\right.
	\end{align*}
	We then define another operator $ \mathcal{F} : V \rightarrow V^{*}$ by $\mathcal{F}(u) :=-\Delta u-F(u) $.
	More precisely, $ \mathcal{F} $ is characterized by
	\begin{align}
		\left<\mathcal F(u),v\right> = \left(\nabla u,\nabla v\right)_{L^2} - \left<F(u),v\right>  \text{~~for~all~~} u,v \in V,
	\end{align}
	where $\left<F(u),v\right> = \int_{\Omega} f(u(x)) v(x) dx$.
	The Fr\'echet derivatives of $ F $ and $ \mathcal{F} $ at $ \varphi \in V $ (respectively denoted by $ {F'_{\varphi}} $ and $ {\mathcal F'_{\varphi}} $) are given by
	\begin{align}
		&\langle F'_{\varphi}u,v\rangle = \int_{\Omega} f'(\varphi(x))u(x) v(x) dx \text{~~for~all~~} u,v \in V
		\shortintertext{and}
		&\langle \mathcal F'_{\varphi}u,v \rangle = \left(\nabla u,\nabla v\right)_{L^2} - \langle F'_{\varphi}u,v \rangle  \text{~~for~all~~} u,v \in V. \label{def:derivativecalf}
	\end{align}
	Under this notation and assumption, we look for positive solutions $ u \in V $ of
	\begin{align}
		\label{main:fpro}
		\mathcal{F}(u)=0,
	\end{align}
	which corresponds to the weak form of \eqref{mainpro:d}.
	We assume that some verification method succeeds in proving the existence of a solution $u \in V$ of \eqref{main:fpro} in $\overline{B}(\hat{u},\rho) := \left\{v\in V : \left\|v-\hat{u}\right\|_{V}\leq\rho\right\}$ for some $ \hat{u} \in V $ and $ \rho >0 $.

	When $\lambda < \lambda_1(\Omega)$, the following theorem is helpful for verifying the nonnegativity of $ u $.

	\begin{theo}\label{theo1}
		Let $ f $ satisfy
		\begin{align}
		\label{f:theo1}
		-f(-t)\leq \lambda t + \displaystyle \sum_{i=1}^{n}a_{i}t^{p_{i}} \text{~~for~all~~}t\geq 0
		\end{align}
		for some $ \lambda < \lambda_1(\Omega) $, nonnegative coefficients $ a_1,a_2,\cdots,a_n $, and subcritical exponents $ p_1,p_2,\cdots,p_n \in (1,p^*)$.
		If
		\begin{align}
		\label{cond:theo1}
		\displaystyle \sum_{i=1}^{n}a_i C_{p_{i}+1}^{2}\left( \left\|\hat{u}_{-}\right\|_{L^{p_{i}+1}}+C_{p_{i}+1}\rho\right)^{p_{i}-1}<1 - \frac{\lambda}{\lambda_1(\Omega)},
		\end{align}
		then the verified solution $u \in V$ of \eqref{main:fpro} in $\overline{B}(\hat{u},\rho)$ is nonnegative.
	\end{theo}

	\begin{rem}
		\label{rem1:theo1}
		The polynomial
		$ f(t) = \lambda t + \sum_{i=2}^{n(<p^*)}a_{i}t|t|^{i-1} $
		with $ \lambda < \lambda_1(\Omega) $ and $a_i \in \mathbb{R} $, which were discussed in the previous section, obviously satisfies the required inequality \eqref{f:theo1}.
		Indeed, for the set of subscripts $ \Lambda_{+} $ for which $ a_i \geq 0~(i \in \Lambda_{+})$ and $ a_i < 0$ $($otherwise$)$,
		we have $-f(-t)\leq \lambda t + \sum_{i \in \Lambda_{+}}a_{i}t^{i} \text{~~for~all~~}t\geq 0 $.
		Moreover, for this polynomial, the positivity of $ u $ follows from the nonnegativity via the maximum principle $ ( $see, for example, {\rm\cite{drabek2009maximum}} for a generalized maximum principle applicable for weak solutions$ ) $.
	\end{rem}

	\begin{rem}
		The formula in parentheses in \eqref{cond:theo1} goes to $0$ as $\hat{u}$ approaches a nonnegative function and $ \rho \downarrow 0 $.
		Therefore, as long as verification succeeds for a nonnegative approximation $\hat{u}$ with sufficient accuracy,
		the nonnegativity of $ u $ can be confirmed using Theorem $\ref{theo1}$.
	\end{rem}

	\begin{rem}
		Even if the approximation $ \hat{u} $ is negative in some parts of $ \Omega $, if it is close enough to a nonnegative function in the sense that $\overline{B}(\hat{u},\rho)$ contains at least one nonnegative function,
		this theorem may work for verifying nonnegativity because it only requires the bound for $ \left\|\hat{u}_{-}\right\|_{L^{p_{i}+1}} $ small enough to satisfy \eqref{cond:theo1}.
		For the same reason, this theorem is applicable for $ \hat{u} $ whose nonnegativity is difficult to prove.
		For example, this theorem is reasonable even when long computation time is required for proving the nonnegativity of $ \hat{u} $ due to its high regularity (see Section {\rm\ref{sec:ex}}).
	\end{rem}

	\subsection*{Proof of Theorem \ref{theo1}}
	First we prove that, for $p\in (1,p^*)$,
	\begin{align}
	\left\|u_{-}\right\|_{L^{p+1}}
	\leq
	\left\|\hat{u}_{-}\right\|_{L^{p+1}}+ C_{p+1}\rho.
	\label{uandu_}
	\end{align}
	We express $u\in V$ as $ u=\hat{u}+ \rho\omega$, where $\omega\in V$ satisfies $\left\|\omega\right\|_{V}\leq 1$.
	Because $b\displaystyle \geq(a-b)_{-}(:=\max\{-(a-b),0\})$ for nonnegative numbers $a,b\in \mathbb{R}$, we have
	\begin{align*}
	0 \leq u_{-}=&\left(\hat{u}+ \rho\omega\right)_{-}
	=\left(\hat{u}_{+}-\hat{u}_{-}+ \rho\omega_{+}- \rho\omega_{-}\right)_{-}
	=\left(\hat{u}_{+}+ \rho\omega_{+}-\left(\hat{u}_{-}+ \rho\omega_{-}\right)\right)_{-}
	\leq \hat{u}_{-}+ \rho\omega_{-},
	\end{align*}
	which implies \eqref{uandu_} because $\left\|\omega_{-}\right\|_{L^{p+1}}\leq C_{p+1}\left\|\omega_{-}\right\|_{V}\leq C_{p+1}$.

	We then prove that the norm of $ u_{-} $ vanishes.
	Because $u$ satisfies
	\begin{align*}
	\left(\nabla u,\nabla v\right)_{L^2}=\left<F(u),v\right>{\rm~~for~all~}v\in V,
	\end{align*}
	by fixing $v=u_{-}$, we have
	\begin{align}
	\left\|u_{-}\right\|_{V}^{2}
	\leq &\displaystyle \int_{\Omega}\left\{\lambda\left(u_{-}(x)\right)^{2}+\sum_{i=1}^{n}a_{i}\left(u_{-}(x)\right)^{p_{i}+1}\right\}dx \nonumber \\
	= &\displaystyle \lambda\left\|u_{-}\right\|_{L^{2}}^{2}+\sum_{i=1}^{n}a_{i}\left\|u_{-}\right\|_{L^{p_{i}+1}}^{p_{i}+1}\nonumber \\
	\leq &\left\{ \displaystyle \frac{\lambda}{\lambda_{1}(\Omega)} +\sum_{i=1}^{n}a_{i}C_{p_{i}+1}^{2} \left\|u_{-}\right\|_{L^{p_{i}+1}}^{p_{i}-1} \right\} \left\|u_{-}\right\|_{V}^{2}.
	\end{align}
	Inequalities \eqref{cond:theo1} and \eqref{uandu_} lead to
	\begin{align}
	\displaystyle \frac{\lambda}{\lambda_{1}(\Omega)} +\sum_{i=1}^{n}a_{i}C_{p_{i}+1}^{2} \left\|u_{-}\right\|_{L^{p_{i}+1}}^{p_{i}-1} < 1,
	\end{align}
	which ensures $ \|u_{-}\|_V = 0 $.
	Therefore, the nonnegativity of $ u $ is proved.
	\qed

	\medskip
	This theorem can be applied to the nonlinearity discussed at the beginning of Section \ref{sect:1} (see again Remark \ref{rem1:theo1}).

	\begin{cor}
		\label{cor:emden}
		Let $ f(t)= \lambda t + t|t|^{p-1} $, with $ \lambda < \lambda_1(\Omega) $ and $ p \in (1,p^*) $.
		If
		\begin{align}
		C_{p+1}^{2}\left(\left\|\hat{u}_{-}\right\|_{L^{p+1}}+C_{p+1}\rho\right)^{p-1}<1 - \frac{\lambda}{\lambda_1(\Omega)},
		\end{align}
		then the verified solution $u \in V$ of \eqref{main:fpro} in $\overline{B}(\hat{u},\rho)$ is positive.
	\end{cor}

	\section{Verification of positivity --- Newton iteration retaining nonnegativity when $\lambda \geq \lambda_1(\Omega)$}
	\label{sec:nk}
	In this section, we discuss another approach to verifying the positivity of a solution $ u $ of \eqref{main:fpro} when $ \lambda \geq \lambda_1(\Omega) $ to which Theorem \ref{theo1} is not applicable.
	For this purpose, we impose another assumption on the nonlinearity $ f $:
	\begin{align}
	\label{asm:f2}
	f(t)\geq f'(t)t ~~~\mbox{for~all}~~~ t\geq 0.
	\end{align}
	The class of functions $ f $ satisfying \eqref{asm:f2} includes polynomials, which are ignored in Theorem \ref{theo1}, of the form
	\begin{align}
	\label{ft:allen}
	f(t)=\lambda t+\displaystyle \sum_{i=2}^{n(<p^*)}a_{i}t|t|^{i-1}~~\mbox{with}~~\lambda \geq \lambda_1(\Omega)~~\mbox{and}~~a_2,\,a_3,\cdots,\,a_n \le 0.
	\end{align}
	This admits some other important cases such as $ f(t) = \varepsilon^{-2}(t-t^3)$ with $ \varepsilon>0 $.
	Recall that there exists no positive solution when $ \lambda < \lambda_1(\Omega) $ and the coefficients $ a_i $ corresponding to super-linearity are nonpositive.

	The method proposed in this section is based on the Newton iteration retaining nonnegativity.
	Before presenting the main theorem of this section, we introduce the affine invariant Newton-Kantorovitch theorem, which ensures the convergence of Newton's method for a ``good'' starting point $ \hat{u} $.
	This theorem has wide applicability to verification methods for nonlinear equations, including differential equations.
	We discuss the positivity of a solution $ u $ whose existence is proved via this theorem.
	\begin{theo}[\cite{deuflhard1979affine}]
		\label{theo:nk}
		Let $ \hat{u} \in V$ be some approximation of a solution $ u $ of $\mathcal F(u)=0$.
		Suppose that there exists some $\alpha>0$ satisfying
		\begin{align}
		\label{alpha}
		||{\mathcal F'_{\hat{u}}}^{-1}\mathcal F(\hat{u})||_V\le\alpha.
		\end{align}
		Moreover, suppose that there exists some $\beta>0$ satisfying
		\begin{align}
		\label{beta}
		||{\mathcal  F'_{\hat{u}}}^{-1}(\mathcal  F'_v-\mathcal  F'_w)||_{{\cal L}(V,V)}\le\beta||v-w||_V~~\mathrm{for~all~} v,w \in D,
		\end{align}
		where $D=B(\hat{u},2\alpha+\delta)$ is an open ball depending on the above value $\alpha>0$ for small $\delta>0$.
		If
		\begin{align}
		\label{alphabeta}
		\alpha\beta \le \frac{1}{2},
		\end{align}
		then there exists a solution $u \in V$ of $\mathcal F(u)=0$ in $\overline{B}(\hat{u},\rho)$ with
		\begin{align}
		\label{nkerror}
		\rho = \frac{1-\sqrt{1-2\alpha\beta}}{\beta}.
		\end{align}
		Furthermore, $ \mathcal  F'_{\varphi} $ is invertible for every $ \varphi \in B(\hat{u},\rho)$, and the solution $u$ is unique in $\overline{B}(\hat{u},2\alpha)$.
	\end{theo}

	The following theorem verifies the nonnegativity of $ u $ the local existence of which is confirmed by Theorem \ref{theo:nk}.
	Here, $ g \in V^* $ is called nonnegative if and only if
	\begin{align}
	\label{def:positive}
	\langle g, \varphi\rangle \geq 0 ~~ \text { for all } \varphi \in V \text { with } \varphi \geq 0.
	\end{align}
	\begin{theo}\label{theo:keeppositive}
		Suppose the following:
		\begin{enumerate}
			\setlength{\parskip}{0cm}
			\setlength{\itemsep}{0cm}
			\item $ f $ satisfies \eqref{asm:f2};
			\item there exist $ \alpha $ and $ \beta $ satisfying \eqref{alpha}, \eqref{beta}, and \eqref{alphabeta};
			\item the minimal eigenvalue $ \mu_1(\hat{u}) $ of ${\cal F}_{\hat{u}}'$ is positive, where $ \mu_1(\hat{u}) $ is given by
			\begin{align}
				\label{def:eigenvalue}
				\mu_1(\hat{u}) = \inf_{v\in V\backslash{\{0\}}} \frac{\langle {\cal F}_{\hat{u}}'v,v \rangle}{\|v\|_{L^2}^2};
			\end{align}
			\item there exists a nonnegative function $ \omega \not\equiv
			0 $  in $ B(\hat{u},\rho) $.
		\end{enumerate}
		Then there exists a nonnegative solution $u \in V$ of $\mathcal F(u)=0$ in $\overline{B}(\hat{u},\rho)$.
	\end{theo}

	\begin{rem}
	Assumption $ 2 $ can be replaced with another condition which proves the convergence of Newton's method in a neighborhood around $ \hat{u} $ and the invertibility $ \mathcal  F' $ in the whole of it.
	For example, in {\rm \cite{nakao2011numerical}} and the references therein, numerical verification methods using several fixed-point theorems were developed.
	By applying an appropriate fixed point theorem such as Banach's fixed point theorem to the Newton operators starting from $ \hat{u} $ on the basis of these methods, another (but probably similar) condition can be used in place of Assumption $ 2 $.
	Note that the conditions proving the convergence of simplified Newton's method, such as in {\rm \cite{plum2008,nakao2011numerical}}, are not directly replaceable with Assumption $ 2 $, because Theorem $ \ref{theo:keeppositive} $ requires the convergence property of the original Newton's method.
	\end{rem}

	\begin{rem}
		Although it is unknown whether Assumption $ 3 $ may hold for all nonlinearities $ f $ satisfying \eqref{asm:f2} for some $ \hat{u} $ near to a positive function,
		Assumption $ 3 $ is expected to be satisfied at least for the nonlinearity \eqref{ft:allen} with $ \hat{u} $ approximating a positive solution of \eqref{main:fpro} with sufficient accuracy because a standard variational method ensures that the desired positive solutions are least-energy solutions.
	\end{rem}

	\begin{rem}
		Assumption $ 4 $ must hold for us to find a nonnegative or positive solution in $\overline{B}(\hat{u},\rho)$.
		In practice, it is useful for us to check $ \max\left\{\hat{u},0\right\} \in B(\hat{u},\rho)$.
	\end{rem}

	\subsection*{Proof of Theorem \ref{theo:keeppositive}}
	Theorem \ref{theo:nk} and Assumption 2 guarantee the existence of a solution $ u $ in $ \overline{B}(\hat{u}, \rho) $.
	Therefore, it remains to prove the nonnegativity of $ u $.

	Assumption 3 ensures that the minimal eigenvalue $ \mu_1(\varphi) $ of $ \mathcal{F}'_{\varphi} $ is positive for all $ \varphi \in B(\hat{u},\rho) $.
	Indeed, if $ \mu_1(\varphi) $ is nonpositive for some $ \varphi \in B(\hat{u},\rho) $, then it follows that there exists $ \varphi_0 \in B(\hat{u},\rho) $ such that $ \mu_1(\varphi_0) = 0 $,
	because $ \mu_1(\varphi) $ is continuous with respect to $ \varphi $.
	This contradicts the result from Theorem \ref{theo:nk}.
	It can be proved from the following discussion that the operator $ \mathcal{F}_{\varphi}^{\prime-1} $ retains nonnegativity for all $ \varphi \in B(\hat{u},\rho) $.
	Let $ w \in V $ satisfy $ \mathcal{F}'_{\varphi} w \geq 0$ in the sense of \eqref{def:positive}.
	By fixing $ v = w_{-}\,(= \max\{-w,0\}) $ in \eqref{def:derivativecalf}, we have
	\begin{align*}
		0 \leq \langle \mathcal{F}'_{\varphi} w, w_{-} \rangle
		= - \left\|w_{-}\right\|_{V}^{2} + \int_{\Omega} f'(\varphi(x))w^2_{-}(x) dx
		= - \langle \mathcal{F}'_{\varphi} w_{-}, w_{-} \rangle.
	\end{align*}
	Therefore, from the definition \eqref{def:eigenvalue} of $ \mu_1(\varphi) $, we have
	\begin{align*}
		\left\|w_{-}\right\|_{L^2}^{2} \mu_1(\varphi) \leq \langle \mathcal{F}'_{\varphi} w_{-}, w_{-} \rangle \leq 0.
	\end{align*}
	The positivity of $ \mu_1(\varphi) $ ensures that $\left\|w_{-}\right\|_{L^2} = 0$.
	Hence, $ w $ is nonnegative.

	We next prove that the Newton iteration in $ B(\hat{u},\rho) $ retains nonnegativity, that is, for every nonnegative $ \varphi \in B(\hat{u},\rho) $, the Newton operator
	\begin{align}
	T(\varphi) := \varphi - \mathcal{F}_{\varphi}^{\prime-1}\mathcal{F}(\varphi)
	\end{align}
	on $ B(\hat{u},\rho) $ maintains nonnegativity.
	For this operator, we have
	\begin{align*}
	T(\varphi)=\varphi-\mathcal{F}_{\varphi}^{\prime-1}\mathcal{F}(\varphi)
	=\varphi-\mathcal{F}_{\varphi}^{\prime-1}\left(\mathcal{F}_{\varphi}'\varphi+F_{\varphi}'\varphi-F(\varphi)\right)
	=\mathcal{F}_{\varphi}^{\prime-1}\left(F(\varphi)-F_{\varphi}'\varphi\right).
	\end{align*}
	Because Assumption 1 ensures that $ F(\varphi)-F_{\varphi}'\varphi \geq 0 $ for $ \varphi \geq 0 $ in the sense of \eqref{def:positive}, it follows that $ T(\varphi) \geq 0$ for nonnegative $ \varphi \in B(\hat{u},\rho) $ (recall that $ \mathcal{F}_{\varphi}^{\prime-1} $ retains nonnegativity).
	Therefore, Assumption 4 ensures the existence of a Newton sequence starting from nonnegative $ \omega \not\equiv 0 $ that converges to a solution $ u \geq0 $.
	\qed

	\medskip
	This theorem can be applied to the special case discussed at the beginning of Section \ref{sect:1}.
	\begin{cor}
		\label{cor:allen}
		Let $f(t)=\varepsilon^{-2}(t-t^3)$, with $ \varepsilon^{-2} \geq \lambda_1(\Omega) $.
		If Assumptions $2,3,$ and $4$ in Theorem $\ref{theo:keeppositive}$ hold, then there exists a positive solution $u \in V$ of $\mathcal F(u)=0$ in $\overline{B}(\hat{u},\rho)$.
	\end{cor}
	\begin{proof}
		The polynomial $ f $ satisfies \eqref{asm:f2}.
		Therefore, Theorem \ref{theo:keeppositive} ensures the nonnegativity of $ u $.
		Its positivity follows from the maximum principle (see, for example, {\rm\cite{drabek2009maximum}}).
	\end{proof}

	\section{Numerical examples}\label{sec:ex}
	In this section, we present examples in which the positivity of solutions to \eqref{main:fpro} are verified via our method.
	All computations were implemented on a computer with 2.90 GHz Intel Core(TM) i9-7920X CPU, 128 GB RAM, and Ubuntu 18.04 using MATLAB 2018a with GCC version 6.3.0.
	All rounding errors were strictly estimated using the toolboxes INTLAB version 10.2 \cite{rump1999book} and kv library version 0.4.47 \cite{kashiwagikv}.
	We constructed approximate solutions of \eqref{main:fpro} for $ \Omega=(0,1)^2 $ from a Legendre polynomial basis.
	More concretely, we constructed $\hat{u}$ as
	\begin{align}
	\displaystyle \hat{u}(x,y)=\sum_{i=1}^{N}\sum_{j=1}^{N}u_{i,j}\phi_{i}(x)\phi_{j}(y),~~u_{i,j}\in \mathbb{R},
	\end{align}
	where each $\phi_{n}$ ($ n=1,2,3,\cdots $) is defined by
	\begin{align}
	\displaystyle \phi_{n}(x)=\frac{1}{n(n+1)}x(1-x)\frac{dQ_{n}}{dx}(x) \text{~~with~~}
	Q_{n}=\displaystyle \frac{(-1)^{n}}{n!}\left(\frac{d}{dx}\right)^{n}x^{n}(1-x)^{n},~~n=1,2,3,\cdots.
	\end{align}
	We define a finite dimensional subspace $ V_N~(\subset V) $ as the tensor product $V_N = \text{span\,}\{\phi_{1},\phi_{2},\cdots,\phi_{N}\} \otimes \text{span\,}\{\phi_{1},\phi_{2},\cdots,\phi_{N}\} $, then defining the orthogonal projection $ P_N $ from $ V $ to $ V_N $ by
	\begin{align*}
		(v- P_N v,v_N)_V = 0 \text{~~~for all~} v \in V \text{~and~} v_N \in V_N.
	\end{align*}
	We used \cite[Theorem 2.3]{kimura1999on} to obtain an explicit interpolation-error constant $ C_N $ satisfying
	\begin{align}
		\label{eq:interpolation}
		\left\|v-P_{N}v\right\|_{V}\leq C_{N}\left\|\Delta v\right\|_{L^{2}} \text{~~for~all~~} v \in V \cap H^2(\Omega).
	\end{align}
	Recall that our method does not limit the basis functions that constitute approximate solutions,
	being applicable to many kinds of bases other than the Legendre polynomial basis, such as the piecewise linear finite element basis or the Fourier basis, etc.

	We proved the existence of solutions $u$ of \eqref{main:fpro} in $ \overline{B}(\hat{u},\rho) $ using Theorem \ref{theo:nk}.
	The key constants $ \alpha $ and $ \beta $ were estimated by
	\begin{align*}
	\alpha \leq \|\mathcal{F}_{\hat{u}}^{\prime-1}\|_{\mathcal{L}\left(V^{*}, V\right)}\|\mathcal{F}(\hat{u})\|_{V^{*}} \text{~~~and~~~}
	\beta \leq \|\mathcal{F}_{\hat{u}}^{\prime-1}\|_{\mathcal{L}\left(V^{*}, V\right)} L,
	\end{align*}
	where $ L $ is a positive number satisfying
	\begin{align*}
	\left\|F_{v}^{\prime}-F_{w}^{\prime}\right\|_{\mathcal{L}\left(V, V^{*}\right)} \leq L\|v-w\|_{V} \text {~~for all~~} v, w \in D.
	\end{align*}
	In the following examples,
	the inverse norm $ \|\mathcal{F}_{\hat{u}}^{\prime-1}\|_{\mathcal{L}\left(V^{*}, V\right)} $ was estimated using the method described in \cite{tanaka2014verified, liu2015framework} in the finite dimensional subspace $ V_N $.
	Moreover, we evaluated the upper bound on $ \|\mathcal{F}(\hat{u})\|_{V^{*}} $ by $C_2 \|\mathcal{F}(\hat{u})\|_{L^2} $,
	where $C_2$ is the constant of embedding $ L^{2}(\Omega) \hookrightarrow V^* $ which in fact coincides with the constant of embedding $ V \hookrightarrow L^{2}(\Omega) $ (see, for example, \cite{plum2008}).
	This $ L^2$-norm was computed using a numerical integration method with strict estimation of rounding errors using \cite{kashiwagikv}.
	The embedding constant $ C_2 $ was calculated as $ C_2=(2 \pi^2)^{-\frac{1}{2}} \approx 0.2251$ with strict estimation of rounding errors.
	Other embedding constants $ C_{p+1}~(p>1) $ were evaluated via the formula in \cite[Corollary A.2]{tanaka2017sharp}.

	For our first example,
	we consider the problem of finding positive solutions to Emden's equation
	\begin{align}
	\left\{\begin{array}{l l}
	-\Delta u=u\left|u\right|^{p-1} &\mathrm{in~} \Omega,\\
	u=0 &\mathrm{on~} \partial\Omega
	\end{array}\right.\label{emden}
	\end{align}
	with $p=3, 5$,
	the approximate solutions of which are displayed in Fig.~\ref{fig1}.
	The Lipschitz constant $ L $ was estimated as
	\begin{align}
	\label{lip:emden}
	L \leq p(p-1) C_{p+1}^{3}\left(\|\hat{u}\|_{L^{p+1}}+C_{p+1} r\right)^{p-2},~~r=2\alpha+\delta~~{\rm for~small}~~\delta>0
	\end{align}
	via a simple calculation from the definition, where we set $ r $ to be the next floating-point number of $ 2\alpha $.
	We verified the positivity of the verified solutions $ u $ using Corollary \ref{cor:emden} (or Theorem \ref{theo1}).
	Table \ref{tab1} shows the verification result.
	In all cases, the positivity of the verified solutions $ u $ were confirmed under the condition $ C_{p+1}^{2}\left(\left\|\hat{u}_{-}\right\|_{L^{p+1}}+C_{p+1}\rho\right)^{p-1} < 1 $.
	It should be noted that the positivity of the approximation $ \hat{u} $ was not proved but, in stead of it, upper bounds for $ \left\|\hat{u}_{-}\right\|_{L^{p+1}} $ were roughly estimated by dividing the domain $ \Omega $ into $ 2^{14} $ smaller congruent squares and implementing interval arithmetic on each of them.
	Whereas one can infer from the shapes of $ \hat{u} $ displayed in Fig.~\ref{fig1} that $ \hat{u} $ ($ p=3,5 $) are positive in $ \Omega $,
	we only used the rough estimation of the negative part to avoid proof of the positivity,
	because it requires much computational cost such as the rigorous calculation for the slope of $ \hat{u} $ near the boundary $ \partial \Omega $.

	\newcommand{0.49\hsize}{0.49\hsize}
	\newcommand{\vspace{0mm}}{\vspace{0mm}}
	\begin{figure}[H]
		\begin{minipage}{0.49\hsize}
			\begin{center}
				\vspace{0mm}\includegraphics[height=45 mm]{p3.pdf}\\
				\footnotesize$p=3$,~~$\displaystyle \max_{x\in\Omega}\hat{u}(x)\approx 6.6232$
			\end{center}
			~
		\end{minipage}
		\begin{minipage}{0.49\hsize}
			\begin{center}
				\vspace{0mm}\includegraphics[height=45 mm]{p5.pdf}\\
				\footnotesize$p=5$,~~$\displaystyle \max_{x\in\Omega}\hat{u}(x)\approx 3.1721$
			\end{center}
			~
		\end{minipage}
		\caption{Approximate solutions to \eqref{emden} on $\Omega=(0,1)^{2}$ for $p=3, 5$.}
		\label{fig1}
	\end{figure}

	\begin{table}[H]
	\caption{Verification results for \eqref{emden} on $\Omega=(0,1)^{2}$ for $p=3, 5$. The values (except for those in row $ N $) represent strict upper bounds in decimal form.}
	\label{tab1}
	\begin{center}
		\renewcommand\arraystretch{1.3}
		\footnotesize
		\begin{tabular}{lll}
			\hline
			$p$& 3& 5\\
			\hline
			\hline
			$N$& 40 & 40 \\
			$ \|\mathcal{F}_{\hat{u}}^{\prime-1}\|_{\mathcal{L}\left(V^{*}, V\right)} $&1.70325176&2.36317681\\
			$ \|\mathcal{F}(\hat{u})\|_{V^{*}} $&$ 2.64173615\ten{-8} $&$ 1.92671579 \ten{-3} $\\
			$ L $&$ 0.67839778 $&$ 6.47198581 $\\
			$\alpha$& $ 4.49954173\ten{-8} $& $ 4.55317005\ten{-3} $\\
			$\beta$& $ 1.15548221 $& $ 15.2944468 $\\
			$\rho$& $ 4.63295216\ten{-8} $&	$ 5.47604979\ten{-3} $\\
			$C_{p+1}$&	$ 0.31830989 $& $ 0.39585400 $\\
			$ \left\|\hat{u}_{-}\right\|_{L^{p+1}} $&$ 4.19109326\ten{-2} $&$ 4.81952900\ten{-2} $\\
			$ C_{p+1}^{2}\left(\left\|\hat{u}_{-}\right\|_{L^{p+1}}+C_{p+1}\rho\right)^{p-1} $& $ 1.77973446\ten{-4} $&$ 1.00813027\ten{-6} $\\
			\hline
		\end{tabular}
	\end{center}
\end{table}

	In our next example, we consider the stationary problem of the Allen-Cahn equation
	\begin{align}
	\left\{\begin{array}{l l}
	-\Delta u=\varepsilon^{-2}(u-u^3) &\mathrm{in~} \Omega,\\
	u=0 &\mathrm{on~} \partial\Omega,
	\end{array}\right.\label{pro:allen}
	\end{align}
	where $\varepsilon>0$.
	We constructed approximate solutions $ \hat{u} $ of this problem using a Legendre polynomial basis in the same way, obtaining the figures displayed in Fig.~\ref{fig3}.
	The Lipschitz constant $ L $ was estimated as
	\begin{align*}
	L \leq 6 \varepsilon^{-2} C_{4}^{3}\left(\|\hat{u}\|_{L^{4}}+C_{4} r\right),~~r=2\alpha+\delta~~{\rm for~small}~~\delta>0
	\end{align*}
	in the same manner as \eqref{lip:emden} with $ p=3 $.
	Using Theorem \ref{theo:nk}, we again obtained $ H^1_0 $-error estimations for solutions of \eqref{pro:allen} centered around these approximations.
	Table \ref{tab3} shows the verification results for $\varepsilon=0.1$, $0.05$, and $0.025$.
	The positivity of the verified solutions was confirmed on the basis of Corollary \ref{cor:allen} (or Theorem \ref{theo:keeppositive}),
	where the required condition $\mu_1(\hat{u}) > 0$ was ensured in all cases.
	The lower bounds on $\mu_1(\hat{u})$ were computed numerically by estimating all rounding errors using the method in \cite{liu2015framework} with the interpolation-error constant $ C_N $ satisfying \eqref{eq:interpolation}.
	Note that proving the positivity of $ \hat{u} $ was also ignored in this example.

	\renewcommand{0.49\hsize}{0.325\hsize}
	\begin{figure}[H]
		\begin{minipage}{0.49\hsize}
			\begin{center}
				\vspace{0mm}\includegraphics[height=36 mm]{allen1.pdf}\\
				\footnotesize{$\varepsilon=0.1$}
			\end{center}
			~
		\end{minipage}
		\begin{minipage}{0.49\hsize}
			\begin{center}
				\vspace{0mm}\includegraphics[height=36 mm]{allen2.pdf}\\
				\footnotesize{$\varepsilon=0.05$}
			\end{center}
			~
		\end{minipage}
		\begin{minipage}{0.49\hsize}
			\begin{center}
				\vspace{0mm}\includegraphics[height=36 mm]{allen3.pdf}\\
				\footnotesize{$\varepsilon=0.025$}
			\end{center}
			~
		\end{minipage}
		\caption{Approximate solutions to \eqref{pro:allen} on $\Omega=(0,1)^{2}$ for $\varepsilon=0.1$, $0.05$, and $0.025$.}
		\label{fig3}
	\end{figure}

	\begin{table}[H]
		\caption{Verification results for \eqref{pro:allen} on $\Omega=(0,1)^{2}$ for $\varepsilon=0.1$, $0.05$, and $0.025$. The values (except for those in rows $ N $ and $\mu_1(\hat{u})$) represent strict upper bounds in decimal form. The values in row $\mu_1(\hat{u})$ represent strict lower bounds in decimal form.}
		\label{tab3}
		\begin{center}
		\renewcommand\arraystretch{1.3}
		\footnotesize
		\begin{tabular}{llll}
			\hline
			$\varepsilon$& 0.1& 0.05&0.025\\
			\hline
			\hline
			$N$& 40 & 40 & 60\\
			$ \|\mathcal{F}_{\hat{u}}^{\prime-1}\|_{\mathcal{L}\left(V^{*}, V\right)} $&$2.85871420$&$4.57367687$&$26.8239159$\\
			$ \|\mathcal{F}(\hat{u})\|_{V^{*}} $&$5.57390453\ten{-10}$&$2.15869521\ten{-6}$&$1.99428443\ten{-6}$\\
			$ L $&$3.00408573$&$5.02704780$&$7.57229904$\\
			$\alpha$&$ 1.59342000\ten{-9}$&$ 9.87317430\ten{-6}$&$5.34945174\ten{-5}$\\
			$\beta$&$ 8.58782250$&$ 22.9920923$&$2.03118712\ten{+2}$\\
			$\rho$&$ 1.59342002\ten{-9}$&$	9.87429519\ten{-6}$&$5.37883476\ten{-5}$\\
			$\mu_1(\hat{u}) \geq$&$ 1.13045870\ten{+2}$&$3.89094819\ten{+2}$&$1.18578735\ten{+3}$\\
			\hline
		\end{tabular}
		\end{center}
	\end{table}

	\section{Conclusion}
	In this paper, we have proposed methods for verifying the positivity of weak solutions $ u $ of the elliptic problem \eqref{mainpro:d} (namely, solutions of \eqref{main:fpro}) using the $ H^1_0 $-error estimation $ \left\|u-\hat{u}\right\|_{H_{0}^{1}} \leq \rho $ for some numerical approximation $ \hat{u} \in V $ and an explicit error bound $ \rho $.
	Theorem \ref{theo1} and \ref{theo:keeppositive} provide sufficient conditions for the solution $ u $ to be nonnegative.
	Using the maximum principle \cite{drabek2009maximum}, the positivity of $ u $ follows from the nonnegativity.
	Our theorems have a wide range of applications, including several important problems such as the elliptic problem \eqref{mainpro:d} with polynomial nonlinearities.
	Numerical examples confirmed that our approach works effectively for several important problems.