\section{Introduction}

Student satisfaction is regarded as a vital component for measuring the success of education service providers such as universities and post-secondary institutions. Proving that higher education (HE) provides a good experience for students is important for promoting and advertising the University both in the home country and internationally. It is considered important to measure the efficiency, quality, impact of the teaching provided by evaluating its effect on students. Student satisfaction ratings play an important role in identifying main areas of strength and weakness for further improvement by the HE service providers. Eventually, it can help them to better meet student expectations by identifying areas of student interests and needs. The HE providers are required to understand the key elements or factors that constitute student satisfaction,  not only around the overall program experience but student satisfaction about their individual course units as well.

Since measuring the overall satisfaction cannot provide an accurate estimation of individual student satisfaction, all instructors are encouraged to survey their students at the end of every course. The collected data can be utilized to enhance the quality of academic institutions \cite{Klemen2015}. The two surveys,  the Student Experience Survey (SES) and the Course Experience Survey (CES) are increasingly being used to measure teaching and university performance in the provision of HE. However, the data collected from student experience and engagement surveys raise questions in terms of quality, reliability, and validity, and may not be used as a piece of concrete evidence in decision-making for higher education institutions \cite{Klem2015}. A recent research investigation has identified the incorporation of many irrelevant questions in the course evaluation surveys which influence the responses and reduce the quality of collected data to be used as an indicator of desirable performance measure \cite{7469785}. Therefore, it is important to identify and analyze relevant factors for informed decision making.

The use of various data-driven solutions are becoming increasingly popular for understanding educational and administrative problems in higher education \cite{5524021,7469785}.
The aim of educational data mining (EDM) is to discover valuable hidden insights from data related to educational services which may be difficult and very time consuming if done manually \cite{MOHAMAD2013320}. The use of EDM can lead towards more benefits and impacts for students, educators and education service providers by upgrading the levels of academic achievement and success for students in a more effective and efficient manner \cite{SHAHIRI2015414}.

In this paper, we utilize state-of-the-art machine learning techniques to investigate five aspects of student satisfaction around course evaluation, class participation, course expectation, course relevancy and professional development (as it relates to real-world employment). Specifically, we build models using student evaluation data to predict each of these satisfaction indicators. We further investigate the crucial attributes associated with each aspect of student satisfaction. Since the satisfaction is very challenging to define and there might be various other factors that influence student satisfaction, we limit our investigation within the course and instructor related predictors. The main research question we address in this paper is as below: \\

\noindent
\textit{``What are the predictors related to instructor and course that influence different aspects of student satisfaction?"}\\

The rest of this paper is organized as follows: Section 2 provides a review of related literature. Section 3 provides a formal definition to the research problem. Section 4 describes the analytic approach used for our experiment and analysis. In Section 5, we show experimental results. Finally, Section 6 concludes the paper.

\section{Review of Related Literature}

The wide use of EDM is evident from the increasing development of methods for exploring large-scale data generated within the higher education sector. A comprehensive survey of the key components of EDM has been presented in \cite{PENAAYALA20141432}. The research identified three categories: tasks, methods, and algorithms as three core building blocks of EDM related research. The use of data mining and machine learning is well studied in educational technology research as well where the aim is to understand how students interact with technology \cite{ANGELI2017226}. The researchers can highlight the existence of several distinct patterns in students' learning experiences.

To analyze the performance of undergraduate students, different data mining techniques have been used in \cite{ASIF2017177}. The study predicted the final academic achievements of the students at the end of their four-year bachelor degree. These authors also analyzed the progression of students throughout their academic years. Another study used predictive analysis to infer the academic performance of students \cite{FERNANDES2018}. Yet another research paper applied data mining techniques on historical student course grade data to model student performance \cite{BURGOS2018541}, demonstrating the determination of course drop out of a student. A tutoring action plan is proposed based on the modeling for the prevention of academic dropouts. Different classification techniques are evaluated to predict the slow learners in an educational setup \cite{KAUR2015500}. The effectiveness of using machine learning techniques to predict the performance of students has been discussed in \cite{7684167}. A atudent evaluation dataset has been leveraged to estimate the performance of the course instructor in \cite{AHMED2016137}. This research uses various machine learning techniques including J48, Multilayer Perception, Na\"{i}ve Bayes, and Sequential Minimal Optimization. Another research paper that investigates the performance of instructors is presented by leveraging the data from course evaluation questionnaire and applying data mining algorithms in \cite{7469785}. This research mainly considers the perception of students and found that there are many irrelevant questions in the course evaluation survey which cannot be related to the instructors' performance. Another research study utilizes clustering and sequential pattern mining techniques to analyze the careers of university graduates by \cite{CAMPAGNI20155508}. This research is based on a real case study which introduces the concept an ideal career to represent the career of a model student. Another review to predict student performance using several machine learning techniques is presented in \cite{SHAHIRI2015414}.

Several research projects have been conducted on various factors to measure student satisfaction after the course has finished. Student satisfaction in a typical computer science course has been researched in
\cite{DBLP:journals/corr/AfrinRRR15}. This research identifies the proper course planning, student encouragement, satisfactory coverage and relevance of the course to the real-world software development scenario as the four key factors of student satisfaction. To understand general measures of student satisfaction from student-opinion data, the regression and decision tree models are utilized in \cite{Thomas2004}. This research identifies distinct predictors for each measure which cannot be used interchangeably. The social integration is found to have an immense effect on student satisfaction for those students who are not highly engaged academically. Another comprehensive use of data mining techniques for gaining insight into student satisfaction is presented in \cite{DEJAEGER2012548}. This research highlights the fact that data mining techniques can be useful to identify and act on a small number of predictors that require attention in order to increase and manage student satisfaction.

The student satisfaction research discussed above considers mainly one aspect of satisfaction. This paper infers and provides an analysis of five different aspects of student satisfaction by applying educational data mining techniques. We aim to investigate course and instructor related factors to identify the ones with the largest impact on student satisfaction.

\section{Research Problem Formulation}
We have two main goals in this research. First, we investigate whether different aspects of student satisfaction can be predicted using various combinations of features related to course and instructor. To achieve our goal, we train a set of machine learning algorithms using a portion of the historical student evaluation data and log the prediction performance for unknown instances. Second, we infer and analyze influential predictors for student satisfaction. We formalize our research problems as below:

\begin{itemize}
    \item Let, $S_{i} = \{1, 2, 3, 4, 5\}$ be the set of five satisfaction levels of satisfaction aspect $S_{i}$, where 1 indicates lowest satisfaction and 5 indicates highest satisfaction. In the student evaluation dataset, each instance $x$ is described by a $d$-dimensional vector of attributes $R^{d}$ related to course and instructor, and one satisfaction level. If $f(.)$ is the satisfaction level prediction function for an unknown instance of $d$-features, $f(.)$ predicts $\hat{S}(x_{q})$ such as $f(x_{q}):R^{d} \rightarrow \hat{S}(x_{q})$  where $\hat{S}(x_{q})$ is the predicted satisfaction level of the query instance $x_{q}$.

    \item Given a set of predictors related to course and instructor, we can identify the factors from the student evaluation dataset that influence the prediction outcomes. We further investigate the extent to which these factors, which are either course or instructor related, drive the student satisfaction.
\end{itemize}

\section{Analytic Approach}

The analytic approach used in this paper has three core components that are interrelated with each other: 1) Investigation of a student evaluation dataset, 2) Definition and prediction of student satisfaction utilizing the student evaluation dataset, and 3) Inference of influential predictors for student satisfaction.

\subsection{The Student Evaluation Data}

\begin{table*}[h!]
\centering
\caption{Description of Fields in the Student Evaluation Dataset}
\label{data-table}
\renewcommand*{\arraystretch}{1.1}
\small
\begin{tabular}{p{1.5cm}p{13cm}p{1.5cm}}
\toprule
Fields  & Description Description of Fields                               &  Values            \\ \midrule
instr      & Instructor identifier                                       & \{1,2,3\}          \\
class      & Course code descriptor                                      & \{1-13\}           \\
repeat     & Number of times the student has taken this course             & \{0,1,2,3\}    \\
attendance & Code of the level of attendance                               & \{0,1,2,3,4\}  \\
difficulty & Level of difficulty of the course as perceived by the student & \{1,2,3,4,5\}      \\
Q1 & The semester course content, teaching method and evaluation system were provided at the start & \{1,2,3,4,5\}\\
Q2 & The course aims and objectives were clearly stated at the beginning.& \{1,2,3,4,5\}\\
Q3 & The course was worth the amount of credit assigned to it. & \{1,2,3,4,5\}\\
Q4 & The course was taught according to the announced syllabus. & \{1,2,3,4,5\}\\
Q5 &	The class discussions, homework assignments, applications and studies were satisfactory. & \{1,2,3,4,5\}\\
Q6 & The textbook and other courses resources were sufficient and up to date.	& \{1,2,3,4,5\}\\
Q7 & The course allowed field work, applications, laboratory, discussion and other studies. & \{1,2,3,4,5\}\\
Q8 & The quizzes, assignments, projects and exams contributed to the learning.	& \{1,2,3,4,5\}\\
Q9 & I greatly enjoyed the class and was eager to actively participate. & \{1,2,3,4,5\}\\
Q10 & Initial expectations about the course were met at the end. & \{1,2,3,4,5\}\\
Q11 & The course was relevant and beneficial to my professional development. & \{1,2,3,4,5\}\\
Q12 & The course helped me look at life and the world with a new perspective. & \{1,2,3,4,5\}\\
Q13 & Instructor's knowledge was relevant and up to date. & \{1,2,3,4,5\}\\
Q14 & Instructor came prepared for classes. & \{1,2,3,4,5\}\\
Q15 & Instructor taught in accordance with the announced lesson plan. & \{1,2,3,4,5\}\\
Q16 & Instructor was committed to the course and was understandable. & \{1,2,3,4,5\}\\
Q17 & Instructor arrived on time for classes. & \{1,2,3,4,5\}\\
Q18 & Instructor has a smooth and easy to follow delivery/speech. & \{1,2,3,4,5\}\\
Q19 & Instructor made effective use of class hours. & \{1,2,3,4,5\}\\
Q20 & Instructor explained the course and was eager to be helpful to students. & \{1,2,3,4,5\}\\
Q21 & Instructor demonstrated a positive approach to students. & \{1,2,3,4,5\}\\
Q22 & Instructor was open and respectful of the views of students. & \{1,2,3,4,5\}\\
Q23 & Instructor encouraged participation in the course. & \{1,2,3,4,5\}\\
Q24 & Instructor gave relevant homework assignments/projects, and helped/guided students. & \{1,2,3,4,5\}\\
Q25 & Instructor responded to questions about the course inside and outside. & \{1,2,3,4,5\}\\
Q26 & Instructor's evaluation (exams, projects, assignments) effectively measured the course objectives. & \{1,2,3,4,5\}\\
Q27 & Instructor provided solutions to exams and discussed them with students. & \{1,2,3,4,5\}\\
Q28 & Instructor treated all students in a right and objective manner. & \{1,2,3,4,5\}\\
\bottomrule
\end{tabular}
\end{table*}

For our experiments and evaluation, we utilize a student evaluation dataset which is publicly available from the UCI repository and was contributed by the students of Gazi University in Ankara, Turkey \cite{GunduzFokoue:2013}. This dataset contains a total of 5820 evaluation scores provided by the students in response to 28 specific questions (Q1 to Q28) related to course and instructor. All the responses to these 28 questions are logged in a 5-point Likert's scale where 1 indicates lowest and 5 indicates highest. The dataset also contains student responses in regards to an additional 5 attributes: instructor code, course code, number of repeats (of this course), level of attendance and perceived difficulty of the course. The response to the difficulty of the course was given between 1 and 5. All these responses were captured for a total of 13 courses. The response to the number of repeats is a number between 0 and 4. The description of different fields of the dataset is given in Table \ref{data-table}. A comprehensive analysis of this student evaluation dataset is provided at the beginning of the experiment and discussion section.

\subsection{Define and Predict Student Satisfaction}
We consider five aspects of student satisfaction. After careful investigation, we decode the meanings of different questions in the student evaluation dataset. We consider that Q8, Q9, Q10, Q11 and Q12 (see Table \ref{data-table}) are the representative questions that capture five aspects of student satisfaction in terms of \textit{course content}, \textit{class participation}, \textit{achievement of initial expectations about the course}, \textit{relevancy towards professional development}, and \textit{the extent to which the course connects and helps the student explore real-world situations}. To predict these five aspects of student satisfaction, we train a set of machine learning techniques including Support Vector Machine (SVM), Multilayer Perceptron (MLP), Decision Tree, Random Forest, Decision Table, $k$-Nearest Neighbour ($k$-NN). Further details will be explained in the Experiment and Discussion section.

\subsection{Calculating Influential Predictors}
\label{influential-predictors}

There are several techniques applied to various domains to compute important predictors for informed decision making \cite{al2012construction,abdullah2012stock,rahaman2017queue,al2012knowledge, rahaman2018wait}. We identify the most influential predictors of student satisfaction using an information theoretic measure known as `Mutual Information' which measures the mutual dependence between two variables. In our case, one of the two variables is a feature while the other variable is the target (i.e. satisfaction level). We consider one aspect of satisfaction at a time during the calculation of mutual information between a prediction and the target. The higher the value of computed mutual information the better the influence of the predictor considered for computation. The mutual information between a predictor and the satisfaction level is calculated as below:

\begin{equation}
I(a_{j};S_i)=-\sum_{j,i} p(a_{j},S_i) \log \frac{p(a_{j},S_i)}{p(a_{j})p(S_i)}
\end{equation}

Here, $a_{j}$ is $j^{th}$ predictor of $i^{th}$ satisfaction aspect and $Si$ is any satisfaction aspect $i$; $p(a_{j},S_i)$ is the joint probability.

\section{Experiment and Discussion}

\subsection{Preliminary Analysis}

\begin{figure*}[t]
\centering
\includegraphics[width=6in, height=3in]{Responses}
 \caption{Distribution of response frequencies for five satisfaction aspects (i.e. Q8, Q9, Q10, Q11, Q12) along with attendance levels, number of repeats and perceived difficulty levels associated with different courses.}
\label{fig:distribution:all-responses}
\end{figure*}

At the beginning of our experiment with student evaluation data, we analyze the distribution of responses logged in the dataset against the five different aspects of student satisfaction. Figure \ref{fig:distribution:all-responses} illustrates the distribution of response frequencies for five satisfaction aspects (i.e. Q8, Q9, Q10, Q11, Q12) along with attendance levels, number of repeats and perceived difficulty levels associated with different courses. It should be noted that a response `5' in a satisfaction aspect indicates the frequency of respondents who are `very satisfied' while a response `1' indicates very `unsatisfied'. A response of `3' indicates the neutral views of the students. We can see an imbalance in the response frequencies for five satisfaction aspects i.e. response `3' is the most dominant response selected by the participant students followed by the response `4'. This issue needs to be considered while training machine learning algorithms since the training performance can be compromised by the label imbalanced problem if present in the dataset. The following subsection explains how we address this challenge.

We also can see from Figure \ref{fig:distribution:all-responses} that a large number of students had different attendance ratings as indicated by the attendance codes between 0 and 4. The attendance code 0 contributes most of the frequencies compared to others. Figure \ref{fig:distribution:all-responses} also highlights that the students mostly undertook the course only once. There are some students who repeated the course for two or three times. However, there are no student participants in this dataset who did not complete this course.

Students also provided their perceptions of the difficulty level of their specific course. As can be seen from Figure \ref{fig:distribution:all-responses}, a large number of students are neutral in their responses in terms of the difficulty level of the course which is represented by a difficulty rating of `3', followed by a large proportion who found it easy represented by `1' (i.e. least difficult). Also, there are a large number of students who found the course 'difficult' or `very difficult' as can be seen from the frequencies of responses `4' and `5'.

\subsection{Student Satisfaction Prediction}

In our experiments, we use the Weka \cite{hall09:_weka_data_minin_softw} toolkit for a set of classifiers including Support Vector Machine (SVM), Multilayer Perceptron (MLP), Decision Tree, Random Forest, Decision Table, and $k$-Nearest Neighbor ($k$-NN). In our Decision Tree implementation, we used information gain measure for top-down induction. An ensemble of 100 trees was used for the Random Forest. Default parameters were considered for the implementation of SVM, MLP and Decision Table. We used $k=5, 10$ for our implementation of the $k$-NN classifier.

To validate the performance of our satisfaction prediction, we use a stratified cross-validation. The reason of using stratified cross-validation is to ensure the presence of an equal number of target (i.e. satisfaction) labels in the training phase of each cross-validation. In other words, this is to reduce the effect of label imbalance in our data. For student satisfaction prediction, we conduct three sets of experiments as follows:

\begin{itemize}
    \item \textit{Predicting satisfaction using predictors related to GTS (Good Teaching Score) only}: GTS is one of the widely used metrics in many countries such as Australia to measure student satisfaction \cite{gts}. A GTS score consists of 6 main features that represent 6 key teaching areas. The responses for calculating the GTS score are collected by asking students to rate the following 6 statements: GTS-1) The teaching staff are extremely good at explaining things; GTS-2) The teaching staff normally give me helpful feedback on how I am going; GTS-3) The teaching staff motivate me to do my best work; GTS-4) The teaching staff work hard to make the course interesting; GTS-5) The teaching staff make a real effort to understand difficulties I might be having with my work; and GTS-6) The staff put a lot of time into commenting on my work. The GTS scores have been shown to correlate well with student satisfaction \cite{gts1}.

    Since the student evaluation dataset used in this paper does not contain all of these 6 questions directly, we consider a total of 9 features that together can be used to describe the actual 6 GTS features. We consider that GTS-1 is equivalent to Q16 and Q18; GTS-2 is equivalent Q25; GTS-3 is equivalent Q28; GTS-4 is equivalent Q22 and Q23; GTS-5 is equivalent Q21; and GTS-6 is equivalent Q24 and Q27.

    The outcome of the student satisfaction prediction using GTS features only is given in Table \ref{table:gtsaccuracy}. Here the maximum prediction accuracies produced by the SVM algorithm for the satisfaction aspect features of Q8, Q10, Q11, and Q12 are 76.01\%, 77.85\%, 76.10\% and 75.34\% respectively. For the overall student satisfaction aspect represented by Q9, $k$-NN ($k$ =10) produces the highest prediction accuracy of 75.94\%.

    \item \textit{Predicting satisfaction using predictors related to instructor only}: In this set of experiments, we consider only those predictors related to instructors. These features are from Q13 to Q28 in the student evaluation dataset. The outcome of the student satisfaction prediction using instructor features only is given in Table \ref{table:q13-q28}. We can see that the maximum prediction accuracy is produced by the SVM algorithm for satisfaction aspect Q8 and is 76.98\%. For the overall student satisfaction aspect Q9, and Q11 $k$-NN ($k$ =10) produces highest prediction accuracies of 75.96\% and 76.89\% respectively. The highest accuracy for Q10 is 78.40\% which is given by the Decision Table algorithm while $k$-NN ($k$ =5) produces the highest prediction accuracy for Q12 (76.03\%).

    \item \textit{Predicting satisfaction using predictors related to both course and instructor}: In these experiments, we consider predictors related to the course (Q1 to Q7) and instructors (Q13 to Q28) together. We can see from table \ref{table:accuracy} that the maximum prediction accuracies produced by the SVM algorithm for the satisfaction aspects Q8, Q10, Q11, and Q12 are 84.16\%, 84.02\%, 80.82\% and 79.86\% respectively. For Q9, the Random Forest algorithm produces the highest accuracy of 80.86\%.
\end{itemize}

\begin{table}[h]
\centering
\caption{Accuracy using Good Teaching Score (GTS) features only}
\label{table:gtsaccuracy}
\renewcommand*{\arraystretch}{1}
\begin{tabular}{llllll}
\toprule
               & Q8    & Q9     & Q10   & Q11   & Q12   \\\midrule
SVM            & \textbf{76.01} & 75.58  & \textbf{77.85} & \textbf{76.10} & \textbf{75.34} \\
MLP            & 73.37 & 72.81 & 76.0 & 73.12 & 73.88 \\
Decision Tree  & 75.24  & 74.52  & 76.82 & 74.97 & 74.78 \\
Random Forest  & 74.91 & 73.42  & 76.77 & 74.28 &  74.07\\
Decision Table & 75.77 & 75.26  & 77.34 & 75.09 & 74.74\\
k-NN (k = 5)   & 75.48 & 75.15  & 77.37 & 75.33 &  74.62\\
k-NN (k = 10)  & 75.96 & \textbf{75.94}  & 77.66 & 76.03 &  75.0\\\bottomrule
\end{tabular}
\end{table}

\begin{table}[h]
\centering
\caption{Accuracy using predictors related to instructors only (Q13-Q28)}
\label{table:q13-q28}
\renewcommand*{\arraystretch}{1}
\begin{tabular}{llllll}
\toprule
               & Q8    & Q9     & Q10   & Q11   & Q12   \\\midrule
SVM            & \textbf{76.98} & 75.89  & 78.37 & 76.54 & 75.96 \\
MLP            & 73.52 & 71.85  & 76.15 & 73.77 & 73.17 \\
Decision Tree  & 75.26 & 74.12  & 77.64 & 75.38 & 74.91 \\
Random Forest  & 75.67 & 74.47  & 77.73 & 76.17 & 75.17 \\
Decision Table & 75.74 & 74.71  & \textbf{78.40} & 75.96 & 75.93\\
k-NN (k = 5)   & 76.27 & 75.57  & 77.92 & 76.36 & \textbf{76.03} \\
k-NN (k = 10)  & 76.39 & \textbf{75.96}  & 78.20 & \textbf{76.89} & 75.84 \\\bottomrule
\end{tabular}
\end{table}

\begin{table}[h]
\centering
\caption{Accuracy of student satisfaction prediction using both course and instructor predictors. Note: GTS features are a subset of instructor related features.}
\label{table:accuracy}
\renewcommand*{\arraystretch}{1}
\begin{tabular}{llllll}
\toprule
               & Q8    & Q9     & Q10   & Q11   & Q12   \\\midrule
SVM            & \textbf{84.16} & 80.69  & \textbf{84.02} & \textbf{80.82} & \textbf{79.86} \\
MLP            & 81.51 & 78.281 & 81.46 & 77.75 & 77.23 \\
Decision Tree  & 82.0  & 78.83  & 81.05 & 78.61 & 77.97 \\
Random Forest  & 84.11 & \textbf{80.86}  & 83.76 & 80.67 & 79.07 \\
Decision Table & 82.32 & 79.85  & 82.65 & 79.79 & 78.45\\
k-NN (k = 5)   & 82.73 & 80.36  & 83.33 & 80.67 & 79.19 \\
k-NN (k = 10)  & 82.44 & 80.45  & 83.13 & 80.27 & 79.38 \\\bottomrule
\end{tabular}
\end{table}

In summary, we see from Tables \ref{table:gtsaccuracy}, \ref{table:q13-q28} and \ref{table:accuracy} that the combination of course and instructor related predictors produces the highest prediction accuracies for student satisfaction prediction. In other words, an improvement of prediction accuracies can be seen by considering the course and instructor related predictors in comparison to results using only GTS features or only instructor related features. This improvement is approximately 5\%-7\% across all the aspects of student satisfaction. This indicates that the course related features have more influence in predicting student satisfaction. We further investigate the most influential predictors of student satisfaction in the following subsection.

\subsection{Inferring Influential Predictors}

\begin{figure*}[t]
\centering
\includegraphics[width=6.2in]{important-predictors-all}
 \caption{Important features for predicting satisfaction aspects (Q8 - Q12)}
\label{fig:important:all}
\end{figure*}

As explained in Section \ref{influential-predictors}, we calculate the mutual information between two variables (i.e. a predictor variable and the satisfaction label) to estimate the influence of any predictor on different aspects of student satisfaction. In this section, we list the top-10 influential predictors related to course (i.e. Q1 to Q7) and instructors (i.e. Q13 to Q28). It should be noted that the corresponding importance scores (i.e. mutual information scores) for all predictors were calculated by building a Random Forest classifier each time. To identify the top-10 predictors, we sorted the importance scores in descending order. We find that the top-10 predictors are a subset of a total of 14 predictors which includes Q1 to Q7, Q13 to Q16, and Q22 to 24. We plot the results in Figures \ref{fig:important:all} (a)-(e) to show the top-10 predictors of Q8, Q9, Q10, Q11, and Q12. There is only one occasion when one predictor related to instructor was found in the top-5 while the other four relate to the course.

 Figure \ref{fig:important:all} (a) shows that there are only four predictors related to instructor in the top-10 list that influence student satisfaction in terms of course content (Q8). All of the top-5 predictors are related to the course itself. Figure \ref{fig:important:all} (b) shows that student satisfaction in terms of class participation (Q9) is influenced by four predictors related to instructor as found in the list of top-10 predictors. Once again, as for Q8, all of the top-5 influential predictors are related to the course itself. There is only one instructor related predictor in the top-5 for predicting student satisfaction in terms of the achievement of initial expectations about the course (Q10) while the other four relate to the course.  This is seen in Figure \ref{fig:important:all} (c). In total there are only four predictors related to instructor in the  top-10 list of influential predictors. None of the instructor related predictors were identified in the top-5 predictors of student satisfaction in terms of the course relevancy for professional development (Q11).  This is illustrated in Figure \ref{fig:important:all} (d) where it can also be seen that only a total of four predictors related to instructor are in the list of top-10 influential predictors. Similarly, no instructor related predictor was found to influence the student satisfaction in terms of the course connecting them to real-world situations (Q12), see Figure \ref{fig:important:all}(e).

We further analyze the inter-correlation between the different aspects of student satisfaction using a Pearson's inter-correlation matrix. We found that there are high positive correlations between the different aspects of student satisfaction (i.e. Q8, Q9, Q10, Q11, Q12) as can be seen from the high correlation coefficients (> 0.8) illustrated in Figure \ref{fig:corr:matrix}. We also found a positive correlation (0.44) between the level of attendance in a course and its perceived level of difficulty.

\begin{figure}[h]
\centering
\includegraphics[width=3.3in]{Capture}
 \caption{Correlations among different aspects of student satisfaction.}
\label{fig:corr:matrix}
\end{figure}

\section{Conclusion}
In this paper we used machine learning techniques to infer student satisfaction from five different perspectives. We showed that the machine learning techniques deployed in our study can predict different perspectives of student satisfaction with accuracies between 80\% and 85\% when trained with both course and instructor related factors. We also investigated the top-10 influential predictors of student satisfaction. We found that the course related predictors dominate the list of top-10 as they fill the top-5 positions. This also implies that the course related factors mostly drive student satisfaction.

This means that for this dataset, the students were satisfied in all five aspects of their course if it allowed fieldwork, applications, laboratory, discussion.  They were also happy if their class discussions, homework assignments, applications and studies were satisfactory. There is a notable influence of up to date course resources and textbooks along with the quality of class discussions, homework, and assignments. The students also prefer the course to be taught according to the syllabus announced on the first day of class with the course aims and objectives stated clearly. Only after these course considerations came class participation under committed understandable instructors who gave relevant assignments and helped, encouraged and guided the students.

\balance
\bibliographystyle{IEEEtran}


\section*{Broader Impact}

{
\small

\section{Method}
\label{sect3:method}

\subsection{Preliminaries}

Reinforcement learning focuses on the setting of an agent exploring in an environment with the aim of learning optimal behavior that can maximize the expected return. At each timestep $t$, the agent selects actions $a_t\in \mathcal A$ w.r.t. its policy $\pi_\theta: \mathcal S\to \mathcal A$ which is parameterized by $\theta$ and receives a reward $r(s_t, a_t)$ and the new state of the environments $s_{t+1}$. The return of a trajectory $\chi=(s_0, a_0,...)$ is defined as the discounted sum of rewards: $R_\chi=\sum_{i=0} \gamma^{i}r(s_i, a_i)$ wherein $\gamma$ is the discount factor.

There are two streams of policy optimization methods that can achieve the goal. The first stream includes the on-policy algorithms, where the behavior policy that samples data is identical to the target policy that needs to be updated. The other stream contains the off-policy algorithms which allow different behavior policy and target policy in sampling and learning. Though these two streams of algorithms have difference properties and performance,
we can still abstract the same procedure in policy optimization:
sampling the data, computing the gradients of the objective w.r.t. the policy's parameters, e.g. the gradient of the expected return $\nabla_\theta J(\theta) = \nabla_\theta \mathop{\mathbb E}_{\chi} [R_\chi]$, and updating the target policy based on the gradients.

In the on-policy setting, taking the standard policy gradient algorithm~\cite{Richard1998Reinforcement} as an example, the gradients is computed based on the advantage: the exceeding return of an action $a_t$ compared to the expected return $v^{\pi}(s_t)= \mathop{\mathbb E} \sum_{t'=t} \gamma^{(t'-t)}r(s_{t'}, a_{t'})$ of all possible actions following the policy $\pi(\cdot|s_t)$ at state $s_t$.
$v^\pi$ is called the state values.
The advantage and the gradient are as follows:
\begin{eqnarray}
  A^{\pi}(s_t, a_t) & = & r(s_t, a_t) + \gamma v^{\pi}(s_{t+1}) - v^{\pi}(s_t), \\
  \nabla_\theta J_{\text{on}}(\theta) & = & \mathop{\mathbb E}_{i} [\nabla_\theta\log \pi_\theta(a_i|s_i)A^{\pi_\theta}(s_i, a_i)].
  \label{eq:on-policy-loss}
\end{eqnarray}
The expectation is conducted over the samples collected by $\pi_\theta$ with index $i$. In modern RL, the advantage and objective are computed in more sophisticated manner via techniques like GAE~\cite{gae_schulman2015high} and clipped surrogate loss~\cite{ppo_algorithm}.
On the other hand, in the off-policy algorithm like Deep Deterministic Policy Gradient (DDPG)~\cite{ddpg_algorithm}, the policy is updated based on the following gradient:
\begin{equation}
    \nabla_\theta J_{\text{off}}(\theta) = \mathop{\mathbb E}_{i} [\nabla_\theta \mu_\theta(a|s_i)\nabla_a Q^{\mu}(s_i, a)|_{a=\mu_\theta(s_i)}],
\end{equation}
wherein $\mu$ is the deterministic policy and $Q^\mu$ denotes the state-action values. SAC~\cite{haarnoja2018soft} follows a similar pipeline that the policy optimization is guided by the approximation of state-action values, but comes with a stochastic policy and an entropy maximization mechanism.

\begin{figure}[!t]
\centering
\centerline{\includegraphics[width=\linewidth]{fig/arch.pdf}}
  \caption{The procedure of (a) conventional single-agent learning and (b) the proposed DiCE. The Collaborative exploration in DiCE leverages multiple agents to sample a shared training batch, while the conventional RL, e.g., PPO~\cite{ppo_algorithm} or SAC~\cite{haarnoja2018soft}, uses a single agent to collect data. Apart from the original task reward, the Diversity Regularization creates another stream of incentives to coordinate the exploration of agents, preventing them from converging to identical behaviors. Note that DR is executed for each agent since the diversity reward is different per agent.
  }
  \label{fig:arch-overall}
\end{figure}

\subsection{Collaborative Exploration}
\label{sect:collaborative-exploration}

As shown in Fig.~\ref{fig:arch-overall}, the Collaborative Exploration (CE) employs of a team of $K$ agents and shares the experiences of them. Different from the multi-agent RL setting~\cite{lowe2017multi} where agents exist in a shared environment and their actions can affect others' states, in our setting each agent lives in an independent copy of the same environment and has no interaction with others when sampling data.
In each training iteration, each of $K$ agents collects a batch of $n$ transitions in the environment. CE merges the experience of each agent and forms a shared training batch with $N = nK$ samples, which is then used to compute the objective as:
\begin{equation}
 J_{\text{CE}}^{(k)}(\theta_k) = \cfrac{1}{K}\sum_{j=1}^K  J^{(j)}(\theta_k),\ \ \forall k=1, ..., K,
\end{equation}
wherein the superscript in $J^{(j)}(\theta_k)$ denotes the data used to compute the objective is collected from agent $j$.
This is the intuition behind the term ``Collaborative Exploration'' since each agent explores in its own part and the aggregation of their efforts breaks the local exploration.

CE introduces mismatching between the behavior policy and the target policy, since the action sampled by an agent $k$ may not subjected to another agent $k'$.
In the off-policy setting, the algorithms already consider the mismatching and thus the mixing of data would not harm the learning.
In on-policy setting, since we use PPO~\cite{ppo_algorithm} as the base of implementation,
the mismatching is not a problem because it has been addressed by the clipped surrogate loss.
It should be reminded that, in the on-policy scenario, the advantage $A^{(j)}$ is computed based on the state values of agent $a^{(j)}$, not the agent $a^{(k)}$ we are training. Indeed, the bias is introduced since the advantage is not subject to $a^{(k)}$. However, it serves as a regularization that can reduce the variance of policy gradient and relieve the estimation error of $A^{(j)}$. Intuitively, it's natural that agent $a^{(k)}$ is not as familiar with the exploration area as agent $a^{(j)}$, thus the advantage calculated by $a^{(k)}$ has high extrapolation error~\cite{fujimoto2018off}.
Accepting the updates that can improve the behavior policy is a good choice that has low variance.

\subsection{Diversity Regularization}

Fujimoto et al.~\cite{fujimoto2018off} point out that the learned policy is highly correlated with the data that trains the policy. Therefore it is inevitable that all the constituent agents gradually become identical during the training with CE, since all the agents are trained with the same shared training data.
If a team of agents has all identical behavior, it becomes equivalent to using parallel actors to sample a large training batch as distributed RL~\cite{espeholt2018impala,dqn_mnih2015human,a3c_algorithm,adamski2018distributed,apex_algorithm}, and CE reduces to the local exploration.
To ensure the effectiveness of CE, we propose the Diversity Regularization (DR), which regularizes the exploration and preserves the diversity of agents.

We first introduce the measure of diversity.
Considering the continuous control tasks we focus on, we use the Mean Square Error (MSE) between the means of Gaussian action distributions produced by two agents as the diversity reward~\cite{hong2018diversity}. Concretely, let $\mu^{(k)}(s_t)$ denote the mean of the action of agent $k$, the diversity reward of agent $k$ against a team with $K$ agents is calculated as:
\begin{equation}
  \label{eq:diversity-reward}
  r^{(k)}_{d}(s_t) =  \cfrac{1}{K-1}\sum_{j=1, j\ne k}^{K}
  ||\mu^{(k)}(s_t) - \mu^{(j)}(s_t)||_2^2  .
\end{equation}

\begin{figure}[!t]
\centering
\centering
\hfill
\subfigure[]{
\centering
\includegraphics[width=0.42\linewidth]{fig/gradient-v2.pdf}
\label{fig:gradient-bisector}
}
\hfill
\subfigure[]{
\centering
\includegraphics[width=0.33\linewidth]{fig/loss}
\label{fig:loss}
}\hfill
  \caption{Illustrations of (a) the aggregation of task gradient and diversity gradient and (b) the proposed Two-side Clip loss for on-policy DiCE.}
  \label{fig:two-side-clip-loss}
\vskip -1em
\end{figure}

DR uses the diversity reward to compute a diversity objective by simply replacing the reward with the diversity in original objective equation. The gradients of such objective forms the \textit{diversity gradient}. DR fuses two streams of gradients using the Feasible Direction Method (FDM) to get the final gradient, which is then used to update the policy network.

We flatten the gradients w.r.t. all parameters into a vector for both objectives and get the task gradient ${\mathbf{g}_t} \in \mathbb R^{|\theta|}$ and the diversity gradient ${\mathbf{g}_d}\in \mathbb R^{|\theta|}$.
As illustrated in Fig.~\ref{fig:gradient-bisector}, we compute the angular bisector of two flattened gradients in the parameter space:
$\mathbf{d}  = \mathcal Z( \mathcal Z({\mathbf{g}_{t}}) + \mathcal Z({{\mathbf{g}_b}}) )$,
wherein $\mathcal Z(\mathbf{x}) =  \mathbf{x} / || \mathbf{x}||_2$ normalizes the input vector to a unit vector.
We use this bisector as the direction of the final gradient, as suggested by \cite{tnb_zhang2019learning} that optimizing toward such bisector can improve both objectives. Please refer to Appendix A for a brief proof.
The magnitude of the final gradient is the average length of two projections of gradients when projected onto $\mathbf{d}$.
We clip the projection of the diversity gradient if it is longer than that of the task gradient to avoid the case when diversity gradient overwhelms the task gradient.
The final gradient after fusing is computed as follows:
\begin{equation}
\mathbf{g}_{\text{final}} = \cfrac{{\mathbf{g}_{t}} \cdot \mathbf{d} + \min({\mathbf{g}_{d}} \cdot \mathbf{d}, {\mathbf{g}_{t}} \cdot \mathbf{d}) }{2} \cdot \mathbf{d}.
\end{equation}
The final gradient is used to update the policy following the conventional process, for example, a Stochastic Gradient Ascent step as $\theta \gets \theta + \alpha \cdot \mathbf{g}_{\text{final}}$ , where $\alpha$ is the learning rate.

\subsection{Implementation Details}
\label{sect:impl}
We implement the framework in both on-policy and off-policy settings.
In on-policy setting, we choose PPO~\cite{ppo_algorithm} as the base and realize CE by mixing the data of $K$ agents after the sampling period in each training iteration.
In off-policy setting, we choose SAC~\cite{haarnoja2018soft} as base and propose two ways to equip SAC with CE. The first way is the \textit{share batch}: in each iteration, we sample $n=N/K$ transitions from the buffer of each of $K$ agents and then merge to a shared training batch with size $N$. The second way is the \textit{share buffer}, where each of $K$ agents can access others' buffers and samples $n$ transitions from each of $K$ buffers to form a personal training batch with size $N$.

The implementation of DR is straightforward for both PPO and SAC. In PPO, we use the discounted diversity return $R^{(k)}_{d,t} = \sum_{t'=t}^{\infty} \gamma^{t'-t}r^{(k)}_{d}(s_{t'})$ to replace the advantage term on Eq.~\ref{eq:on-policy-loss}.
We can also leverage an diversity value network (DVN) to compute the diversity advantage which is used to compute the diversity objective.
However, ablation study presented in Sect.~\ref{sect:ablation-study} shows that using DVN is not the optimal choice since the diversity values are hard to learn due to the instability of diversity.
In SAC, we leverage an extra ``diversity critic'' along with the original critics and use its predicted values to form the diversity gradient.
We also introduce some training designs that make DiCE perform better:

\textbf{Delayed Update Target}. To stabilize the diversity metric, we use a common trick called Delayed Update Targets~\cite{ddpg_algorithm}: a set of target policies is maintained by Polyak averaging the parameters of the learning policy of each agent over the course of training:
\begin{equation}
  \theta^{(k)}_{\text{target}} \gets (1 - \tau) \theta^{(k)}_{\text{target}} + \tau \theta^{(k)}_{\text{latest}},\ \ \forall k=1,...,K,
\end{equation}
wherein $\tau$ is a hyper-parameter that close to 0. We then compute the diversity of a given agent $k$ against such target policies, including the delayed update target of $k$-th agent itself.

\textbf{Two-side Clip Loss}. We notice that the original PPO loss is not bounded when the advantage is less than zero and can goes to negative infinity if the action probabilities ratio is large enough, as shown in Fig.~\ref{fig:loss}.
The unbounded loss is vulnerable and leads to high variance of the policy gradient, which is recently noticed~\cite{ye2019mastering}.
To tackle this issue, we utilize a Two-side Clip (TSC) loss that mitigates the variance of advantage by clipping the loss if the ratio is too large when the advantage is negative, as the TSC loss in Fig.~\ref{fig:loss}. Let $\rho^{(k, j)}$ denotes the action probabilities ratio between agent $k$ and agent $j$, the TSC loss equipped by CE is computed as follows:
\begin{equation}
 L^{(k)}_{TSC}(\theta) = \cfrac{1}{K}\sum_{j=1}^{K}  \mathop{\mathbb E} [\text{clip}(\rho^{(k,j)}, 0, 1+\epsilon)A^{(j)} ].
 \label{eq:ce-loss-2-tsc}
\end{equation}

\section{Experiments}
\label{sect:experiment}

\subsection{Setup}

We implement DiCE framework in both on-policy and off-policy settings using RLLib~\cite{liang2017rllib}.
In our preliminary experiments, we find A3C is unstable and sometimes fails the training, so we use A2C, the synchronized version of A3C instead~\cite{a3c_algorithm}. Therefore, we compare the on-policy DiCE with A2C and PPO~\cite{ppo_algorithm}, as well as a diversity-seeking method TNB~\cite{tnb_zhang2019learning}, though it trains agents sequentially so the time consumption is much larger than ours. For off-policy setting, we compare the off-policy DiCE (DiCE-SAC) with SAC~\cite{haarnoja2018soft}, a powerful and sample efficient method with strong exploration capacity.
A DiCE team will collect the same amount of samples in one training iteration as single agent does. Concretely, supposing the agent in single-agent algorithm samples $N$ transitions in each training iteration, then each agent in DiCE team with $K$ agents is asked to collect $n = N/K$ samples to form a shared training batch with totally $N$ data points. In on-policy setting, the samples are from interactions with environment, while in off-policy setting, they are collected by the share batch (denoted as ``bat.'') or the share buffer (``buf.'') approaches as described in Sect.~\ref{sect:impl} .
We train our agents in five locomotion tasks HalfCheetah-v3, Ant-v3, Walker2d-v3, Hopper-v3, and Humanoid-v3 in MuJoCo simulator~\cite{mujoco_todorov2012mujoco}. All the experiments are repeated three times with different random seeds. The shadow of each curve denotes the standard deviation of the values.
In the default setting, DiCE contains five agents in the team, and we present the average episode reward of the best agent in the team at each timestep. The hyper-parameters can be found in Appendix B.

\subsection{Comparisons}

We first compare DiCE framework in on-policy setting.
As shown in Fig.~\ref{fig:main-performance}, in all the five tasks, our method achieves better results compared to the baselines PPO and A2C. In Table~\ref{table:performance}, we see that in four environments DiCE achieves a substantial improvement over the baselines,
while in the Hopper-v3 PPO and TNB achieve higher score than DiCE.
In Hopper-v3, PPO collapses after a long time of training, while DiCE maintains its performance until the end of the training, which shows that DiCE is robust in training stability.

As shown in Fig.~\ref{fig:main-performance-sac} and Table~\ref{table:performance}, in off-policy setting, DiCE-SAC outperforms the SAC baseline in Hopper-v3 and Humanoid-v3 with faster convergence while achieves comparable performance in HalfCheetah-v3 and Walker2d-v3. In Ant-v3, the DiCE-SAC fails to progress compared to SAC. This might because that Ant-v3 environment has loose constraints on action and has larger action space, thus the structure of diversity is more complex than other environments, making the learning of diversity critic harder. We have the similar observation for on-policy DiCE when utilizing a diversity value network as discussed in Sect.~\ref{sect:ablation-study}.
The performance improvements brought by DiCE in on-policy and off-policy settings shows its generalization ability.

We also conduct experiments on four sparse reward environments in a set of mini-grid worlds~\cite{gym_minigrid}. Experiment results show that our method achieves comparable performance with PPO algorithm and presents better robustness than baselines.
We point out that DiCE is not exclusive with intrinsic reward methods~\cite{rnd_burda2018exploration,eysenbach2018diversity} that can boost performance in sparse reward tasks, and thus a hybrid of two approaches can take advantages from both worlds. Please refer to Appendix C for more information.

\begin{figure}[!t]
\centering
\begin{minipage}{0.22\linewidth}
\centering
\includegraphics[height=1in]{fig/diceppo/ant}
\end{minipage}\hfill
\begin{minipage}{0.19\linewidth}
\centering
\includegraphics[height=1in]{fig/diceppo/halfcheetah}
\end{minipage}\hfill
\begin{minipage}{0.19\linewidth}
\centering
\includegraphics[height=1in]{fig/diceppo/hopper}
\end{minipage}\hfill
\begin{minipage}{0.19\linewidth}
\centering
\includegraphics[height=1in]{fig/diceppo/humanoid}
\end{minipage}\hfill
\begin{minipage}{0.19\linewidth}
\centering
\includegraphics[height=1in]{fig/diceppo/walker}
\end{minipage}
  \caption{The learning curves of A2C, PPO, and on-policy DiCE in 5 MuJoCo environments.}
  \label{fig:main-performance}
\end{figure}
\begin{figure}[!t]
\centering
\begin{minipage}{0.22\linewidth}
\centering
\includegraphics[height=1in]{fig/dicesac/dicesac_result_Ant-v3.pdf}
\end{minipage}\hfill
\begin{minipage}{0.19\linewidth}
\centering
\includegraphics[height=1in]{fig/dicesac/dicesac_result_HalfCheetah-v3.pdf}
\end{minipage}\hfill
\begin{minipage}{0.19\linewidth}
\centering
\includegraphics[height=1in]{fig/dicesac/dicesac_result_Hopper-v3.pdf}
\end{minipage}\hfill
\begin{minipage}{0.19\linewidth}
\centering
\includegraphics[height=1in]{fig/dicesac/dicesac_result_Humanoid-v3.pdf}
\end{minipage}\hfill
\begin{minipage}{0.19\linewidth}
\centering
\includegraphics[height=1in]{fig/dicesac/dicesac_result_Walker2d-v3.pdf}
\end{minipage}
  \caption{The learning curves of SAC and off-policy DiCE in 5 MuJoCo environments.}
  \label{fig:main-performance-sac}
\end{figure}
\begin{table}[!t]
\centering
\caption{The episode reward of the best agents occurring in the training at five MuJoCo tasks. We report the averaged results of three repeated experiments.}
\label{table:performance}
\begin{center}\begin{small}
\begin{tabular}{@{}cccccc@{}}
\toprule
                        & Ant-v3  & HalfCheetah-v3 & Hopper-v3 & Humanoid-v3 & Walker2d-v3 \\ \midrule
TNB                     & 2211.02 & 1623.68        & \textbf{2916.85}   & 493.20      & 2906.36     \\
A2C                     & 951.15  & 600.26         & 280.89    & 352.26      & 477.80      \\
PPO                     & 1230.82 & 1216.41        & 2632.40   & 453.09      & 2670.99     \\
DiCE                    & \textbf{2517.65} & \textbf{2860.45}        & 2512.59   & \textbf{566.88}      & \textbf{3537.80}     \\
\midrule
SAC                     & \textbf{4825.35} & \textbf{7847.73}        & 3536.95   & 5046.93     & 4293.26     \\
DiCE-SAC (bat.)  & 2729.52 & 7807.74        & 3671.16   & 5598.92     & \textbf{4784.59}     \\
DiCE-SAC (buf.) & 1906.16 & 6205.99        & \textbf{3683.04}   & \textbf{5627.55}     & 4640.48     \\ \bottomrule
\end{tabular}%
\end{small}\end{center}
\end{table}

\subsection{Impact of the team size}
To verify the effectiveness of collaborative exploration, we vary the agent number $K$ in DiCE team to see the impact of the team size.
Note that when $K=1$, the agent computes diversity against the delayed update target of itself.
The performance is improved when the number of agents is increased, as shown in Fig.~\ref{fig:varing-the-number-of-agents}. When the team size exceeds a certain number, then the performance decreases.
The phenomenon supports our intuition: the joint exploration area expands as the team size increases, leading to more efficient learning. However, when the number of agents is too large, encouraging diversity among them jeopardizes the performance because one or more agents would learn diverse but weak behaviors, dragging down the whole training process.

\begin{figure}[!t]
\centering
\subfigure[\hskip -2em]{
\centering
  \includegraphics[width=0.45\linewidth]{fig/num-agents-walker}
\label{fig:varing-the-number-of-agents}
}\hfill
\subfigure[\hskip -2em]{
\centering
  \includegraphics[width=0.45\linewidth]{fig/ablation-1}
  \label{fig:main-ablation}
}\\%
\vskip -0.5em%
\subfigure[\hskip -2em]{
\centering
\includegraphics[width=0.45\linewidth]{fig/ablate-other-walker-formal-0601.pdf}
\label{fig:other-ablation-walker}
} \hfill
\subfigure[\hskip -2em]{
\centering
\includegraphics[width=0.45\linewidth]{fig/ablation-2-ant.pdf}
\label{fig:other-ablation-ant}
}%
\vskip -0.5em%
\caption{(a) The impact of agent number $K$ in DiCE team. The performance reaches the highest point when $K=7$, but then decreases when $K=10$. (b) Ablation studies on the major components of DiCE. (c) and (d) Ablation studies on some minor designs in two environments.
}
\label{fig:ablation-2}
\end{figure}

\subsection{Ablation Studies}
\label{sect:ablation-study}

To understand which component catalyzes DiCE, we conduct ablation study on some components in our framework. Results in Fig.~\ref{fig:main-ablation} show that
purely using CE already boosts the exploration and leads to considerable improvement. Disabling CE (``w/o CE'') drops the performance drastically and soon collapses the learning as infinity values happen.
On the other hand, employing DR achieves relatively small improvement.
This result confirms the previous statement that diversity introduced by different random initialization can promote performance~\cite{dqn_bootstrapped_osband2016deep}.

We also conduct ablation studies on other minor components as shown in Fig.~\ref{fig:other-ablation-walker} and Fig.~\ref{fig:other-ablation-ant}.
We find that the Two-side Clip Loss (TSC Loss) and the Delayed Update (DU) of target policies have similar and relatively small impact to the results. The Normalization of Advantage (NA) greatly damages the performance and downgrades DiCE to baseline.
The influence of NA indicates that the information contained in the relative magnitudes of the task gradient and the diversity gradient is crucial.
The Diversity Value Network (DVN) harm the policy seriously in Ant-v3, though this is not notable in Walker2d-v3. This is because the Ant-v3 has loose action constraints that causes the diversity to become more unstable and hard to be learned by the diversity value network. The poor performance of DiCE-SAC in this environment shown in Fig.~\ref{fig:main-performance-sac} supports the statement, since DiCE-SAC leverages an extra diversity critic to approximate the diversity values. Finally, we conduct a case study analyzing the reason why DiCE works in Appendix D.

\section{Conclusion}
\label{sect5:conclusion}
In this work we propose a novel non-local policy optimization framework called Diversity-regularized Collaborative Exploration (DiCE). DiCE implements the concept of teamwork in reinforcement learning and utilizes a team of heterogeneous agents to collaboratively explore the environment while maintaining the diversity of the agents. We implement DiCE in both on-policy and off-policy settings, where the experimental results show that DiCE improves the exploration and outperforms the PPO, A2C, SAC and diversity-seeking baselines in most of the cases in MuJoCo locomotion tasks.

\section{Related Work}
\label{sect:related-work}

\textbf{Exploration and exploitation.} It is a long-standing problem to balance the exploration and the exploitation in RL. The exploration enables the agent to visit under-discovered space of the environment by taking action following a non-optimal policy. The exploitation, on the contrary, lets the agent behave greedily to maximize the expected return.
To enable more thorough exploration, many works add noise to several parts of the policy, such as the output action~\cite{td3_fujimoto2018addressing,ppo_algorithm,dqn_mnih2015human,ddpg_algorithm,apex_algorithm}, the parameters of the policy~\cite{parameter_noise_plappert2017parameter,parameter_noise2_fortunato2017noisy,parameter_noise3_ruckstiess2010exploring}, or the value function~\cite{dqn_bootstrapped_osband2016deep,random_value_function_osband2017deep,xu2019vase}.
The other stream of works encourages exploration through explicit optimization, such as maximizing the entropy of actions to promote new action \cite{soft_q_haarnoja2017reinforcement,haarnoja2018soft}, maintaining a visitation map, which records the frequency of states being visited, to motivate to rarely-visited states~\cite{tang2017exploration,ostrovski2017count}, and maximizing information gain about the agent's learned environmental dynamics~\cite{houthooft2016vime}.
Furthermore, many works introduce concepts like intrinsic reward~\cite{rnd_burda2018exploration,stadie2015incentivizing,bellemare2016unifying}, curiosity~\cite{burda2018large,pathak2017curiosity} or surprise~\cite{achiam2017surprise} to encourage exploration.
In all of these practices, though the agent can access the under-explored area due to the randomness or additional incentives, the agent's \textit{exploration area} is still based on the agent's previous experience because in one way or another, the stochasticity is applied on top of the agent's current optimal policy.
On the contrary, the proposed DiCE is a non-local exploration framework where the diversity of the team serves as a source of extra information that an independent agent is impossible to access.

\textbf{Exploration with shared experience.} The straightforward idea to expand the exploration area is to leverage multiple agents to explore the environment collaboratively and share the experience with each other, as if how people team up and solve the problem together~\cite{cohen1994restructuring,smith1992collaborative}.
Schimitt et al.~\cite{schmitt2019off} present an algorithm that shares the experience of previous agents to the learning agent, which can accelerate the hyper-parameter tuning. We don't maintain a shared replay buffer inherited across a sequence of agents, instead we share the experiences of concurrently training agents.
Mnih et al.~\cite{a3c_algorithm} propose A3C and A2C that use multiple independent actors to explore the environments. Both methods maintain a global agent who receives the gradients from parallel actors and episodically broadcasts the latest parameters to them. A3C receives the gradient asynchronously while A2C synchronously.
IMPALA~\cite{espeholt2018impala} is similar to A3C but is different in that it mixes samples instead of gradients from different actors and updates the global agent based on the shared samples.
Compared to A3C, A2C, and IMPALA, our work draws a decentralized framework which maintains a group of heterogeneous agents that neither contribute to nor reset parameters from a global agent.
Besides, our algorithm introduces a strong regularization that maintains the the diversity of the agents and keeps the shared data informative and diverse for efficient learning.

Bootstrapped DQN~\cite{dqn_bootstrapped_osband2016deep} (BDQN) shares a similar idea with ours. It uses multiple heads building upon a shared network to fit multiple Q functions, each Q function corresponds to a policy.
A shared replay buffer built by the samples collected by all policies is used to train all Q functions.
DiCE follows a similar pattern: it samples data by all agents and trains all agents by the shared data.
Different from BDQN, we do not replay the value network to correct the advantages when training agent $k$ using the samples belonging to another agent $k'$. Besides, we use an entirely independent neural network for each agent.
Additionally, we explicitly encourage diversity among different agents, while BDQN achieves diversity from the different initialization of Q functions.

\textbf{Learning with diversity.}
It has a long history in the evolution algorithms community to address the exploration problem in the perspective of population-based learning~\cite{conti2018improving,doncieux2013behavioral,pugh2016quality,lehman2011abandoning}.
Instead of generating new agents via mutation and gene crossover~\cite{beyer2002evolution}, Conti et al.~\cite{conti2018improving} introduce a diversity measure and use it to add perturbation on the policy's parameters. We also introduce diversity rewards among agents, but use the diversity gradient to update policy, which is more efficient than evolution. Besides, our method does not require hand-craft behavioral representation to compute diversity.
In deep RL community, many works encourage diversity explicitly through adding extra loss to make an agent behave differently to a population~\cite{hong2018diversity,masood2019diversity}; adding noise in parameters to create conjugate policies that used for improving the main policy~\cite{cohen2019diverse}; or adding diversity as explicit reward~\cite{eysenbach2018diversity}.
Zhang et al.~\cite{tnb_zhang2019learning} propose a method called task novelty bisector (TNB) which boosts the diversity (called novelty in the original paper) at the gradient level and improves the policy in both performance and novelty.
 The diversity regularization in DiCE is different in the following aspects: (1) It trains multiple agents concurrently. (2) DiCE uses a simplified and computationally efficient formulation of diversity. (3) Our method does not seek to encourage diversity but only uses it as a regularizer to maintain a certain level of diversity and boost the learning of the shared training batch.

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

\begin{wrapfigure}{R}{0.45\textwidth}
    \centering
\includegraphics[width=\linewidth]{fig/teaser.pdf}
\caption{The proposed DiCE encourages the agents to explore state-action space collaboratively and share the experience while regularizes their behaviors to be sufficiently diverse.
}
\label{fig:teaser}
\vskip -1em
\end{wrapfigure}

Working together in a team towards a common goal brings synergy for higher efficiency and better performance~\cite{cohen1994restructuring,smith1992collaborative}.
However, in most of the existing Reinforcement Learning (RL) algorithms, usually only one agent~\cite{Richard1998Reinforcement} or a global agent with several replicas~\cite{a3c_algorithm,espeholt2018impala} explore the environment by trial and error for optimizing the policy.
The agent usually limits its exploration within a small region of the state-action space due to the initialization and previous experience, as illustrated by the light area in Fig.~\ref{fig:teaser}.

Many works have been proposed to encourage sufficient exploration in the single-agent setting. For example, some works direct the agent to access unknown areas through taking action based on a stochastic distribution like Gaussian distribution~\cite{td3_fujimoto2018addressing,ppo_algorithm,dqn_mnih2015human} or following the Ornstein-Uhlenbeck process~\cite{ddpg_algorithm}.
Some other works use maximum entropy~\cite{soft_q_haarnoja2017reinforcement,haarnoja2018soft}, intrinsic reward~\cite{rnd_burda2018exploration}, extra loss~\cite{hong2018diversity}, parameter noise~\cite{parameter_noise_plappert2017parameter,parameter_noise2_fortunato2017noisy,parameter_noise3_ruckstiess2010exploring}, surprise~\cite{achiam2017surprise}, and curiosity~\cite{burda2018large} to encourage the agent to visit rarely visited states. However, for the above algorithms, the possible \textit{exploration area} covered is bounded in a limited subspace determined by the agent's behavior, which is learned from its previous experience.  We call this as the \textit{local exploration} since the exploration area inevitably misses those rarely visited states.

Another way to expand exploration space is through distributed RL~\cite{dqn_mnih2015human,adamski2018distributed,apex_algorithm}. In typical algorithms like A3C~\cite{a3c_algorithm} and IMPALA~\cite{espeholt2018impala}, multiple homogeneous actors generate a large amount of trajectory data independently by interacting with the parallel environments.
However, inflating the training batch with similarly-behaving actors is still the local exploration since only the improvements at the nearby area of current policy can be realized~\cite{tessler2019distributional,ciosek2019better}.
Motivated by the distributed RL, we address the non-local exploration problem with a new method called the \textit{Collaborative Exploration} (CE) which employs a team of multiple heterogeneous agents and shares the experience of each agent among the team.
CE is distinguished from existing distributed RL algorithms in the following aspects:
(1) CE is \textit{decentralized} so that it neither maintains a global agent that receives updates from other agents nor broadcasts the latest policy parameters to others~\cite{a3c_algorithm,espeholt2018impala}. (2) CE leverages multiple \textit{heterogeneous} agents that can update individually based on the shared experiences and the diversity among them is preserved via a regularization mechanism.
(3) CE enables the \textit{non-local exploration} thus breaks the limits of the previous single-agent setting, which is verified by the performance improvement in the extensive experiments.

Since all the agents are trained with the same shared data, one issue is that the diversity of agents gradually vanishes as the training goes.
To resolve this, we further introduce the \textit{Diversity Regularization} (DR) based on Feasible Direction Method that can improve the policies while preserves the diversity of agents. Compared to previous works that train a set of diverse agents in a sequential manner~\cite{hong2018diversity,masood2019diversity,tnb_zhang2019learning,conti2018improving}, DR modulates the diversity of agents that are concurrently trained, which takes much less time than the sequential training. Besides, we design a compact diversity reward, avoiding maintaining a set of auto-encoders~\cite{tnb_zhang2019learning} or hand-crafted representation~\cite{conti2018improving}.

To summarize, we formulate a novel policy optimization framework called Diversity-regularized Collaborative Exploration (DiCE). DiCE combines the Collaborative Exploration that shares knowledge across multiple agents as well as the Diversity Regularization that directs the exploration of each agent. DiCE is implemented in both on-policy and off-policy settings and is compared with baselines e.g. PPO~\cite{ppo_algorithm} and SAC~\cite{haarnoja2018soft}. The experimental results show that DiCE outperforms both on-policy and off-policy baselines in most cases in the MuJoCo locomotion tasks

\section*{Appendix}
\renewcommand{\thesubsection}{\Alph{subsection}}

 In this document, we present the appendix of paper \textit{Non-local Policy Optimization via Diversity-regularized Collaborative Exploration}.

\subsection{Feasible Direction Method}
We resemble the Feasible Direction Method (FDM) in constrained optimization to address the multi-objective optimization issue in our work.
The key idea of FDM is to select a direction that can optimize both objective functions, i.e., the reward objective $J_r(\theta)$, and the diversity objective $J_d(\theta)$. Given a sufficiently small $\lambda>0$, the FDM considers the Taylor series of the two objectives at point $\bar\theta$ as:

\begin{equation}
    J_r(\bar\theta + \lambda \vec{p}) = J_r(\bar\theta) + \nabla_\theta J_r(\bar\theta)^{\rm{T}} \lambda \vec{p} + O(||\lambda \vec{p}||),
\end{equation}

\begin{equation}
    J_d(\bar\theta + \lambda \vec{p}) = J_d(\bar\theta) + \nabla_\theta J_d(\bar\theta)^{\rm{T}} \lambda \vec{p} + O(||\lambda \vec{p}||).
\end{equation}

A gradient ascent direction $\vec{p}$ is called a feasible directioin if it satisfies $\nabla_\theta J_r(\bar\theta)^{\rm{T}}  \vec{p} >0 $ and $\nabla_\theta J_d(\bar\theta)^{\rm{T}}  \vec{p} >0$. As $O(||\lambda \vec{p}||)$ can be omitted, we have
$
    J_r(\bar\theta + \lambda \vec{p}) > J_r(\bar\theta)
$
and
$
    J_d(\bar\theta + \lambda \vec{p}) > J_d(\bar\theta)
$.

In the RL scenario, we wish to maximize the objectives.
It's obvious that the final gradient $\vec{p}$ we used in Eq.~(6), whose direction is the bisector of two gradients, satisfies $\nabla_{\theta}J_r(\bar\theta)^{\rm{T}} \vec{p} >0 $ and $\nabla_{\theta}J_d(\bar\theta)^{\rm{T}} \vec{p} >0$, and therefore the proposed DR module can increase the reward objective as well as diversity objective.

\subsection{Hyper-parameters}

For all experiments, we use fully-connected neural networks with two layers, 256 hidden units per layer, for both policy networks and value networks. We use ReLU as the activation function. Other hyper-parameters in both on-policy and off-policy setting are listed as follows:

\begin{table}[H]
\centering
\caption{Environment-agnostic hyper-parameters of on-policy DiCE}
\label{tab:env-agn-hyperpara}
\begin{tabular}{@{}ll@{}}
\toprule
Hyper-parameter             & Value  \\ \midrule
Number of Agents      &  5 \\
KL Coefficient              & 1.0    \\
$\lambda$                   & 1.0 \\
$\gamma$                    & 0.99 \\
Number of SGD iterations    & 10     \\
Learning Rate               & 0.0001 \\
Delayed Update Coefficient ($\tau$)   & 0.005  \\
Max Norm of Gradient        & 10.0   \\
Use Diversity Value Network & False  \\
Use Normalize Advantage     & False  \\
Use Delayed Update          & True   \\
Use Two-side Clip Loss      & True   \\
Random Seeds                & 0, 100, 200 \\ \bottomrule
\end{tabular}
\end{table}

\begin{table}[H]
\centering
\caption{Environment-related hyper-parameters of on-policy DiCE}
\label{tab:env-rel-hyperpara}
\begin{tabular}{@{}llllll@{}}
\toprule
Hyper-parameter             & HalfCheetah-v3 & Ant-v3 & Walker2d-v3 & Hopper-v3 & Humanoid-v3 \\ \midrule
Sample Batch Size           & 50             & 50     & 200         & 200       & 50          \\
Training Batch Size         & 2048           & 2048   & 10000       & 10000     & 2048        \\
SGD Mini-batch Size         & 64             & 64     & 100         & 100       & 64          \\
Maximum Steps & 5$\times 10^6$          & 5$\times 10^6$  & 5$\times 10^7$       & 5$\times 10^7$     & 2$\times 10^7$       \\ \bottomrule
\end{tabular}
\end{table}

\begin{table}[H]
\centering
\caption{Environment-agnostic hyper-parameters of off-policy DiCE}
\label{tab:env-agn-hyperpara-off}
\begin{tabular}{@{}ll@{}}
\toprule
Hyper-parameter             & Value  \\ \midrule
Number of Agents      &  5 \\
$\lambda$                   & 1.0 \\
$\gamma$                    & 0.99 \\
Training Batch Size    & 256     \\
Maximum Steps & 1$\times 10^6$ \\
Steps that Learning Starts & 10000 \\
Learning Rate               & 0.0003 \\
Delayed Update Coefficient ($\tau$)   & 0.005  \\
Random Seeds                & 0, 100, 200 \\ \bottomrule
\end{tabular}
\end{table}

\subsection{Experiments on sparse reward environments}

\begin{figure}[H]
\centering
\begin{minipage}{0.45\linewidth}
\centering
\includegraphics[height=0.6\linewidth]{fig/sparse-reward/spare-reward-MiniGrid-DoorKey-5x5-v0.pdf}
\end{minipage}
\begin{minipage}{0.45\linewidth}
\centering
\includegraphics[height=0.6\linewidth]{fig/sparse-reward/spare-reward-MiniGrid-Empty-8x8-v0.pdf}
\end{minipage}\\\vskip 1em%
\begin{minipage}{0.45\linewidth}
\centering
\includegraphics[height=0.6\linewidth]{fig/sparse-reward/spare-reward-MiniGrid-FourRooms-v0.pdf}
\end{minipage}
\begin{minipage}{0.45\linewidth}
\centering
\includegraphics[height=0.6\linewidth]{fig/sparse-reward/spare-reward-MiniGrid-UnlockPickup-v0.pdf}
\end{minipage}
  \caption{The learning curves of PPO and on-policy DiCE agents trained in 4 sparse reward environments. Note that only when reward close to 1 can we say a policy solves the task. PPO and DiCE solve DoorKey and Empty environments but fail in FourRooms and UnlockPickup when trained with one million steps from the environments.}
  \label{fig:sparse-reward-performance}
\end{figure}

We conduct experiments on few MiniGrid environments~\cite{gym_minigrid} which the agent navigates in a grid world to achieve the goal. The environments provide sparse rewards: the agent would only receive +1 reward when it reaches the goal, otherwise 0. Different environments have different obstacles and mechanisms to create difficulty for agents to reach the goal. For example, in MiniGrid-UnlockPickup-v0 environment, the agent has to pick up a box which is placed in another room, behind a locked door. The sparse reward setting makes the tasks extremely hard to be solved.

The Fig.~\ref{fig:sparse-reward-performance} shows the performance of PPO and on-policy DiCE in 4 MiniGrid environments. Experimental results show that our method achieves comparable performance to PPO. In MiniGrid-Empty-8x8-v0 and MiniGrid-DoorKey-5x5-v0 environment, the agent achieves 100\% success rate as what PPO does. We observe two phenomenons on the results: (1) DiCE trains a little slower than PPO. It achieves 0.95 reward at 100k steps while PPO does this at 60k steps. (2) PPO shows strong oscillation even after it first achieves 0.95 reward but DiCE maintains its performance. This shows that DiCE is more robust than PPO in these cases, which is also found in Hopper-v3 environment.

We point out that DiCE is not exclusive with intrinsic reward methods~\cite{rnd_burda2018exploration,eysenbach2018diversity} that can boost performance in sparse reward tasks and thus a hybrid of two ideas can take advantage from both methods. We leave such investigation to future works.

\subsection{Empirical Analysis on DiCE}

\begin{figure}[H]
\centering
\hfill
\subfigure[\hskip -2em]{
\centering
\includegraphics[height=1.25in]{fig/walker-ratio}
\label{fig:case-study-ratio}
}\hfill
\subfigure[\hskip -2em]{
\centering
\includegraphics[height=1.25in]{fig/walker-entropy}
\label{fig:case-study-entropy}
}\hfill
\subfigure[\hskip -2em]{
\centering
\includegraphics[height=1.25in]{fig/walker-cosine}
\label{fig:case-study-cosine}
}
\hfill
\caption{(a) The ratio of the action probabilities in Eq.~9. (b) The average entropy of action distribution. (c) The cosine of two gradients.
}
\vskip -1em
\end{figure}

We conduct case studies to investigate how DR and CE help DiCE achieve better performance in on-policy setting. Fig.~\ref{fig:case-study-ratio} presents the ratio between the probabilities of the same action when sampling from the target policy and the behavior policy, as the one in original PPO loss~\cite{ppo_algorithm}.
As a reminder, the behavior policy in the denominator of the ratio is not only the non-delayed-update version of the learning policy: it may also be other agent's policy due to CE.
The ratio represents the discrepancy between the target policy and the behavior policy. If the ratio is closed to 1, then the information in the sample can be efficiently learned since PPO clipping mechanism is not active.
Fig.~\ref{fig:case-study-ratio} shows that our method, though part of the learning samples is from other agents, can still efficiently utilize the data.

We observe that training agents using DiCE without CE can break the learning, which is also shown in the ablation studies. Without the balance of CE, simply encouraging agents to purse diversity would drive them far away from each other, which further increases the magnitude of the diversity. This forms a self-reinforced loop and leads to infinite diversity. Concretely, the action distributions, which are Gaussian in our setting, become extremely narrow and peak at different means. This is supported by the descending entropy in Fig.~\ref{fig:case-study-entropy} and the drastic increasing ratio in the orange line in Fig.~\ref{fig:case-study-ratio}, since a narrow distributions indicates a low entropy and concentrated actions.

DiCE without DR is another story. We observe the entropy is not descent and even ascends in the green ``w/o DR'' curve in Fig.~\ref{fig:case-study-entropy}. Intuitively, suppose there are two policies, say policy $\pi_A$ and $\pi_B$, that are distinct from each other. Consider a state that two good actions $a_A$, $a_B$ from two policies exist, which have high advantages and also high probabilities $\pi_A(a_A|s_t)$, $\pi_B(a_B|s_t)$. An optimizer of $\pi_A$ would try to flatten the action distribution and makes the probabilities of taking both actions high. For a Gaussian distribution, such optimization would lead to high standard deviation and thus high entropy as shown in the green line in Fig.~\ref{fig:case-study-entropy}. Let's consider the impact of DR. The diversity would be large since two policies are far away, and at this time, the diversity gradient of $\pi_A$ is directing to the opposite of $\pi_B$ since it needs to consistently increase the diversity. Therefore the diversity gradient balances the task gradient that seeks to move $\pi_A$ toward $\pi_B$. This hypothesis is verified by Fig.~\ref{fig:case-study-cosine}, which shows that the cosine of two gradients in parameter space are not closed to 1 and that, in most cases, the diversity gradient is directing to different directions as task gradient does. That's the reason why we call DR a ``regularizer'' instead of a ``booster'' --- we use DR to maintain a certain degree of diversity instead of improving it without limit, in which case it may be harmful because it distracts agents from the prime task and leads to the same consequence just like the case without CE.


\section{Introduction}
Ehrenfest's celebrated theorem showed that the classical Newton
concept of force applies equally well to the acceleration of
isolated quantum systems, but what about forces between parts
of extended quantum systems?  The question is relevant to
the description of the dynamics of nuclear fission. At some
point the nucleus is elongated into
a shape resembling the nascent fission fragments joined by
a small neck.  The Coulomb force between the two nascent
fragments is counterbalanced by an attractive nuclear force
transmitted through the neck region; scission only occurs
when neck force is sufficiently weakened.

In this work, we show that in an extended system the quantum force
of one part on the rest can be consistently defined.
In general, the interaction has long-range components
such as the Coulomb and short-range ones such as the exchange-correlation
energy in density-functional theory.  Except for the
long-range components of the interaction, the  force
reduces to an  operator acting at
the interface between the two subsystems. The operator
contains derivatives of the wave function at the interface
due to the kinetic term in the Hamiltonian.

We first consider the Schr\"odinger equation for a  particle in one
dimension.  We then generalize the operator scope  to many-particle
systems in three dimensions.

\section{Particle in one dimension}
\subsection{Kinetic force}
Consider a particle in one dimension and governed by the Schr\"odinger
Hamiltonian $H = K + V$. Here $K$ is the kinetic Hamiltonian
\begin{equation}
K = -{\hbar^2\over 2M} {\partial^2\over\partial x^2}.
\label{K-def}
\end{equation}
where $M$ is the mass of the particle.
$V$ in the Hamiltonian is a potential depending on $x$. We divide
the space
into two subsystems at some point $x_0$.  The particle number
$n_{R,L}$ on each side and the positions of their centers of
mass $x_{cm(R,L)}$ are evaluated by restricting the integrations
over the wave function
to one side or
the other of the point $x_0$.  For the right hand side, this is
achieved by the operators

\begin{equation}
\hat {\mathbb{1}}_{R} = \Theta(x-x_0)
\end{equation}
\begin{equation}
\hat x_R = x \Theta(x - x_0)
\end{equation}
where $\Theta$ is the Heaviside step function.  We write the expectation values
in the wave function $\phi(x)$ as
\begin{equation}
\langle \phi | \hat {\mathbb{1}}_R | \phi\rangle = n_R
\end{equation}
\begin{equation}
\label{xR}
\langle \phi | \hat x_R | \phi \rangle = n_R x_{cm(R)}.
\end{equation}

Depending on the wave function, both of these quantities can
vary with time under the Schr\"odinger dynamics.  However, it is
not useful to call the effect on a subsystem a force if the
particle number is changing. We therefore restrict the
definition to wave functions for which $ d n_R / d t$ and $d^2
n_R / d t^2$ are zero at the time when the force is computed.
Then the time dependence of Eq. (\ref{xR}) will provide the
acceleration $a_R$ of the right-hand side computed as
\begin{equation}
a_R = {1\over n_R} {d^2 \over d t^2} \langle \phi | \hat x_R | \phi
\rangle.
\end{equation}
Given the acceleration, he Newtonian force can be defined as
\begin{equation}
\label{F-def}
F_R = M {d^2 \over d t^2} \langle \phi | \hat x_R | \phi \rangle.
\end{equation}
We can now carry the
derivatives in the Heisenberg representation by twice applying the
usual commutator formula
\begin{equation}
{d \over dt} \langle \hat {\cal O}\rangle   = {1\over i\hbar} \langle
[{\cal O},H]\rangle.
\end{equation}
The resulting operator for the kinetic energy term in $H$  is
\begin{equation}
\label{FK-op}
\hat F_k(x_0)  =  -M[K,[K,\hat x_R]] =
 -{\hbar^2 \over 4 M}
\left(2{\overleftarrow\partial \over \partial
x}\delta(x-x_0){\overrightarrow\partial \over \partial x}
- {\overleftarrow {\partial^2}  \over \partial x^2}\delta(x-x_0)
-\delta(x-x_0){\overrightarrow {\partial^2}  \over
\partial x^2}\right)
\end{equation}

To show how Eq. (\ref{FK-op}) works, we apply it to some very
simple Hamiltonians.  The first is a particle in a box.  The
n-th stationary state in the box has the wave function
\begin{equation}
\phi = \sqrt{2\over L} \sin (k_n x )
\end{equation}
where $L$ is the length of the box and $k_n = n\pi/L$.  Its
energy is
\begin{equation}
E_n = {\hbar^2 k_n^2 \over 2 M}.
\end{equation}
The force exerted by the particle on the boundary can be
calculated by the usual relation between force, energy,
and displacement
\begin{equation}
\label{F-box}
F = - {d E_n(L)\over d L} = {\hbar^2 \pi^2 n^2\over M L^3}.
\end{equation}
This force must be sustained throughout the interior of the box.  If we
arbitrarily divide the space in two the same force must act between
the two sides.  Applying Eq. (\ref{FK-op}) we find
\begin{equation}
F_R = {\hbar^2\over M L} \left(\left( \left. {d\sin (k_n x)\over dx}
\right|_{x_0}\right)^2
 - \sin(x_0) \left. {d^2 \sin(k_n x)\over d x^2}\right|_{x_0}\right)
=  {\hbar^2 \pi^2 n^2\over M L^3}
\end{equation}
in agreement with Eq. (\ref{F-box}).

Another simple example is the expansion of a Gaussian wave packet.
The initial wave function may be written
\begin{equation}
\label{phi-ho}
\phi(x) = {1\over b^{1/2}\pi^{1/4}} e^{- x^2/2 b^2}.
\end{equation}
Here $b$ is a parameter controlling the width of
the Gaussian.  The time-dependent Schr\"odinger equation has the
exact solution
\begin{equation}
\label{phit}
\phi(x,t) = {1\over b^{1/2}\pi^{1/4}}(1 + i \hbar t/M b^2)^{-1/2} \exp(
-{\textstyle\frac{1}{2}}  x^2/
 b^2(1+i \nu t/M))).
\end{equation}
There is no particle transfer at $x_0=0$ so we can calculate the
force at that point. The center-of-mass position of the right-hand
side particle density is
\begin{equation}
\label{xRt}
\langle x_{R}\rangle = \int_{x_0}^\infty dx x |\phi(x,t)|^2 = {b \over 2 }
\pi^{1/2} (1 + (\hbar t/Mb^2)^2)^{1/2}
\end{equation}
By explicit differentiation of Eq. (\ref{xRt}), the acceleration of
right-hand density distribution is
\begin{equation}
\label{ddotxR}
{d^2 \over d t^2}\langle \hat x_R \rangle = {1 \over 2 \pi^{1/2} M^2b^3} (1 + (\hbar
t/M b^2)^2)^{-3/2}.
\end{equation}
The force calculated from this acceleration together with  Eq. (\ref{F-def})
agrees with that obtained by the expectation value of
Eq. (\ref{FK-op}) in the wave function $\phi(x,0)$.
The important point is that Eq. (\ref{FK-op}) only
requires information about the
wave function around the point $x_0=0$, and yet it
perfectly describes the cm acceleration
of all the matter to the right of that point.

\subsubsection{Potential contribution}
Now we add an external potential field $V$ to the Hamiltonian and treat the
associated force $F_V$ acting on a region in the same way.
The operator requires the double commutator
\begin{equation}
\hat a_R = -\frac{1}{\hbar^2} [V,[K,\hat x_R]]
\end{equation}
and is evaluated as
\begin{equation}
\label{FV-def}
\langle \hat a_R\rangle =\int_{x_0}^\infty dx\, \rho(x) {\partial V \over \partial
x};
\end{equation}
Not surprisingly, this formula
is very similar to
Ehrenfest's second equation, $d \langle p \rangle /dt =
-\langle d V/dx \rangle$. To illustrate
the application of
Eq. (\ref{FV-def}), consider a particle in the ground state of the
harmonic oscillator Hamiltonian.  It is convenient to express the
potential in the form
\begin{equation}
V(x) = \frac{\hbar^2}{2M b^4}   x^2.
\end{equation}
Its ground state is the wave function of Eq. (\ref{phi-ho}).
Since it is stationary, the acceleration is zero everywhere and there
should be an exact force balance
at all points:
\begin{equation}
F_K + F_V = 0 \,\,\,\,\,\  {\rm for~ all} \,\,x_0.
\end{equation}
It is an easy exercise to verify that
\begin{equation}
\label{FK-ho}
F_K =   {\hbar^2 \nu^{3/2} \over 2 M \pi^{1/2}} e^{-\frac{1}{2}\nu x_0^2}.
\end{equation}
Carrying out the integration
in Eq. (\ref{FV-def}),  one obtains the negative of Eq.
(\ref{FK-ho}).  This  verifies the
force balance for the ground state, but in fact
it must be true for excited states as well.

\subsection{Beyond the one-particle Hamiltonian}
The formulas of the last subsections are  easily generalized to
multiparticle systems when particles interact through  mean-field
potentials.  For the kinetic quantum force, all the physics is governed
by one-body operators and so all of the forces are additive.  The summation
can be carried out at the level of the wave function to obtain
the the single-particle density matrix $\rho(x,x')$.  The corresponding kinetic force is then
given by the compact expression
\begin{equation}
\label{density-matrix}
F_K = -{\hbar^2\over 2 M} \left.{ d^2 \over d s^2} \rho(x_0 +s/2,x_0-s/2)\right|_{s=0}.
\end{equation}
An even more compact expression makes use of the Wigner representation
of the density matrix,
$f(x,p) = \frac{1}{2 \pi}\int_\infty^\infty d s
\exp(i p s) \rho(x +s/2,x-s/2)$.  Then Eq. (\ref{density-matrix})
becomes
\begin{equation}
F_K = {\hbar^2\over M}\int_\infty^\infty d p p^2 f(x_0,p)
\end{equation}

To introduce particle-particle interactions,
we first note that the instantaneous force can be derived by
expanding
the Hamiltonian evolution operator to second order in time. In
mean-field theory the single-particle potential can be expanded
as well. However, the time-dependent corrections to the
instantaneous potential start with terms beyond second order in time
and so can be dropped.
All that is required to treat interaction effects is determine
their contribution with an instantaneous single-particle potential.

The simplest case to deal with is a two-body finite-range interaction
$v(x -x')$.  By Eq. (\ref{FV-def}) the force $F_V$ may be expressed
\begin{equation}
F_V = M\int_{x_0}^\infty dx \,\rho(x) \int_{-\infty}^\infty dx'
\rho(x')
{d \over d x} v(x - x').
\end{equation}
Next split the integral over $x'$ into two at the integration point $x$.
The contribution with $x'$ in the range $x<x'<\infty)$ can be shown to
vanish;
physically the interactions between particles in
the same region do not affect the center-of-mass motion in that region.
Thus $F_V$ reduces to
\begin{equation}
\label{long-range}
F_V = M \int_{x_0}^\infty dx \,\rho(x) \int_{-\infty}^{x} dx' \rho(x')
{d \over d x'} v(x'' - x')
\end{equation}
This form is appropriate as it stands for long-range forces such as
the Coulomb interaction.   It needs the entire density distribution
to calculate it, but at least it doesn't require numerical differentiation
of global energies.

As the range of the interaction decreases, it is clear that only the
the density near $x_0$ contributes to the integral.  An easy way to
derive the force is to follow the spirit of
density function theory, where the interaction energy density
${\cal V}$ is treated as a local function of position.
We start with a simple form for the interaction energy functional,
${\cal V} = \frac{1}{2} v_0 \rho^2(x)$.
The corresponding single-particle potential is
\begin{equation}
V(x) = {d \,{\cal V}(x) \over d \rho} = v_0 \rho(x).
\end{equation}
Inserting this in Eq (\ref{long-range}), we have
\begin{equation}
\label{FV-v0}
F_V = v_0\int_{x_0}^\infty  dx  \rho(x)  {d\rho \over dx} =
{\textstyle\frac{1}{2}} v_0  \int_{x_0}^\infty  dx  {d \over dx} \rho^2(x)
= -{\textstyle\frac{1}{2}} v_0  \rho^2(x_0).
\end{equation}
To obtain the far right-hand equality,
we have assumed that $\rho \rightarrow 0$ at large
$x$.

Eqs.  (\ref{density-matrix}) and (\ref{FV-v0}) can be easily tested in the
one-dimensional Fermi gas model.   We first derive the force from the
total energy in a box. The total force on the box wall can
be calculated as before by taking the derivative of the total energy
with respect to box size.  The result is
\begin{equation}
E =  {\hbar^2 k_F^2 \over 6 M}N + {v_0 \over 2 L} N^2
\end{equation}
where $N$ is the number of particles, $k_F \approx  \pi N/L$ and
$\rho = N/L$.
 Taking the derivative, the force on the wall is
\begin{equation}
F_b = -{d E \over d L} =  {\hbar^2 k_F^2 \over 3 M } {N\over L} + {\textstyle\frac{1}{2}} v_0
{N^2 \over L^2}.
\end{equation}
The first term is the kinetic contribution;  it may be calculated
from Eq. (\ref{FK-op}) taking $(\rho(x+s/2,x-s/2)= \sin(k_F s)/(\pi s)$.
The second is identical to Eq. (\ref{FV-v0}).

The same method can be used to derive the force associated with
any energy density functional $\cal V$
that can be expanded in powers of the local density $\rho(x)$.
We write
\begin{equation}
{\cal V}(\rho) = \sum_m {1\over m!} v_m \rho^m.
\end{equation}
Then $F_V$ can be evaluated similarly to Eq. (\ref{FV-v0}) as
\begin{equation}
F_V = \sum_m {1\over m!} (m-1) v_m \rho(x_0)^m.
\end{equation}
This can also be expressed as \cite{be77}
\begin{equation}
F_V = V(x_0)\rho(x_0) - {\cal V}(x_0).
\end{equation}
In the last form, the force is seen to depend only on the density
at the interface.

\section{Three dimensions}

The generalization to three dimensions is trivial if the interface
between the two subsystems is a plane.  Then the kinetic
force operator acts perpendicular to the plane and just requires an
integration over the transverse coordinates.  Formally, one can
define a stress tensor $\Pi$ that transmits momentum from one part
of the system to another.  In the co-moving frame of the medium,
the stress tensor $\Pi$ associated with density-functional dynamics
is given by an expression very similar to the one-dimensional formula,
\begin{equation}
\Pi_{ij}(\vec r)  =  {\hbar^2\over M}\nabla_{s_i} \nabla_{s_j} \rho(\vec r - \vec s/2,
\vec r + \vec s/2) + \left(V(\vec r) \rho(\vec r) -{\cal V}(\vec
r)\right)\delta_{ij}.
\end{equation}
The two terms represent the kinetic and interaction contributions,
respectively.
The interaction term is isotropic, but that need not be the case for
the kinetic term.  A perpendicular Newtonian force can be calculated across any
plane by the integral
\begin{equation}
\vec F\cdot \vec u  = \int \Pi \cdot d \vec A
\label{3D-force}
\end{equation}
where $\vec u$ is a unit vector perpendicular to the plane $\vec A$.
In practice, one would choose a plane going through
the neck region that joins the two nascent fragments.  Note that
there is no mechanism here to generate a transverse force between
subsystems.

In practice in calculating dynamics in nuclear physics,
condensed matter physics and
quantum chemistry, one defines the configurations
by minimizing an energy functional in the presence of
a fixed external field.
In such situations the wave function has no currents
so the conditions for calculating the
force across a plane are satisfied.  However, the constrained
minimization procedure also permits a third way to calculate
force.  This is to use the Feynman-Hellman theorem and  calculate
the energy derivative as an expectation value of the derivative
of the constraining field.  This is a much easier task than to
explicitly calculate total energy derivatives numerically.  Still,
the quantum force operator might still be useful in some situations
and as an independent check on computations carried out by other
methods.

\section{Appendix}

\subsection{Pairing}

Pairing is very important in low-frequency nuclear dynamics.
At present, the
most well-justified models incorporating pairing are based on the
Hartree-Fock-Bogoliubov (HFB) extension of mean-field theory.  It is now possible to carry out the
integration of the HFB equations of motion without uncontrolled
approximation \cite{bu16}.  It is also clear from simplified implementations
that lifetimes are strongly dependent on the pairing field
\cite{ne78,ro14}.  The force associated with pairing can be derived in the
same way as we treated the ordinary interactions in the Hamiltonian.
The dynamics is governed by the HFB equation of motion, also known as
the Bogoliubov-de Gennes equation.  This is written
\begin{equation}
i \hbar {d \over d t} \left(\begin{array}{c}
\vec v_\alpha \\
\vec u_\alpha \\
\end{array}
\right)=
\left(\begin{array}{cc}
H & \Delta\\
-\Delta^* &  -H\\
\end{array}\right)
\left(\begin{array}{c}
\vec v_{\alpha} \\
\vec u_{\alpha} \\
\end{array}
\right),
\end{equation}
in the usual notation.
  The pairing energy associated with a two-body
interaction $\bar v$
can be expressed $ E_p  = \frac{1}{4} Tr_2 \kappa^* \bar v \kappa $
 \cite{RS} where
$\kappa = \sum_\alpha v_\alpha u^T_\alpha$ ; the
$\Delta$ field is given by $\Delta = \frac{1}{2} \bar v \kappa$.  The expressions for
the energy and field are very
similar to those for the ordinary interaction with the replacement of
$\rho$ by $\kappa$.  We expect that the derived force will come out
in a similar way if the interaction $\bar v$ is short-ranged.  Then
the force would be equal to the pairing energy density at the
division point. For most physical systems, the pairing energy density is
small compared to  other interaction terms,
so the
pairing force can be neglected in practice.

\subsubsection{Adiabaticity}

For the complete dynamics of 3-D media
it is essential to understand the time scale of the collective
motion with respect to the time scale to establish local
equilibrium.  For slow collective motion, the presence of
interactions beyond mean field keeps the Fermi surface nearly spherical
and the resulting stress tensor is nearly isotropic.  In the opposite limit,
the stress tensor remembers the strain history of the system,
and the Fermi surface can have quadrupolar distortion.
In terms of the Lam\'e parameterization of the stress tensor,
in the adiabatic case the compressibility is governed by a pressure
field
\begin{equation}
P = \lambda + \textstyle\frac{2}{3} \mu
\end{equation}
whereas in the diabatic case the longitudinal stress for a
strain field in the $z$ direction is given by
\begin{equation}
\Pi_{zz} = \lambda + 2 \mu.
\end{equation}
The kinetic stress tensor for a diabatically deformed Fermi surface
has Lam\'e coefficients
\begin{equation}
\lambda = \mu = {k_f^2 \over 5 M}.
\end{equation}
A number of Fermionic systems exhibit diabatic
dynamics in high-frequency oscillations.  We mention zero
sound in liquid $^3$He \cite{wi67}, the wave-length dependence
of plasmons in conductors \cite{ja67}, and the giant quadrupole resonance
in large nuclei \cite{lh98}.

These considerations only indirectly affect the forces we have
calculated here.  As stated earlier, the force or stress tensor
depends only on the instantaneous state of the system.  If that
wave function is obtained by a constrained mean field solution with
density constraint operators, it will have no local
currents although  it may have a deformed Fermi surface as a result
of the constraints.  If one releases the constraints and allows
the system to evolve, the proto-fragments will at first be
accelerated away from each other by the force dynamics treated
here.  But later as the state of system changes, the compressibility
will play a role.  A stronger restoring force will be present
in the diabatic dynamics than in the adiabatic.  The considerations
discussed here cannot tell us whether the system eventually
come apart in fragments.

\section{Acknowledgment}
The author would like to thank T. Kawano, W. Younes, and A.
Bulgac for discussions on nuclear fission dynamics.