{'fraction_non_alphanumeric': 0.03310820731170166, 'fraction_numerical': 0.010106200753682768, 'mean_word_length': 4.890314211588354, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
VarFA: A Variational Factor Analysis Framework For Efficient Bayesian Learning Analytics

We propose VarFA, a variational inference factor analysis framework that extends existing factor analysis models for educational data mining to efficiently output uncertainty estimation in the model's estimated factors. Such uncertainty information is useful, for example, for an adaptive testing scenario, where additional tests can be administered if the model is not quite certain about a students' skill level estimation. Traditional Bayesian inference methods that produce such uncertainty information are computationally expensive and do not scale to large data sets. VarFA utilizes variational inference which makes it possible to efficiently perform Bayesian inference even on very large data sets. We use the sparse factor analysis model as a case study and demonstrate the efficacy of VarFA on both synthetic and real data sets. VarFA is also very general and can be applied to a wide array of factor analysis models.


INTRODUCTION
A core task for many practical educational systems is student modeling, i.e., estimating students' mastery level on a set of skills or knowledge components (KC) [67,11]. Such estimates allow in-depth understanding of students' learning status and form the foundation for automatic, intelligent learning interventions. For example, intelligent tutoring systems (ITSs) [57,21,3,47,25,52] rely on knowing the students' skill levels in order to effectively recommend individualized learning curriculum and improve students' learning outcomes. In the big data era, student modeling is usually formulated as an educational data mining (EDM) problem [54,4,41] where an underlying machine learning (ML) model estimates students' skill mastery levels from students' learning records, i.e., their answers to assessment questions.
Many student modeling methods have been proposed in prior literature. A fruitful line of research for student modeling follows the factor analysis (FA) approach. FA models usually assume that an unknown, potentially multidimensional student parameter, in which each dimension is associated with a certain skill, explains how a student answers questions and is to be estimated. Popular and successful FA models include item response theory (IRT) [66], multi-dimensional IRT [2], learning factor analysis (LFA) [6], performance factor analysis (PFA) [50], DASH (short for difficulty, ability, and student history) [35,45], DAS3H (short for difficulty, ability, skill and student skill history) [10], knowledge tracing machines (KTM) [68], and so on. Recently, more complex student models based on deep neural networks (DNN) have also been proposed [51,72,43]. Nevertheless, thanks to their simplicity, effectiveness and robustness, FA models remain widely adopted and investigated for practical EDM tasks. Moreover, there is evidence that simple FA models could even outperform DNN models for student modeling in terms of predicting students' answers [69]. Because of FA models' competitive performance and its elegant mathematical form, we focus on FA-based student models in this paper.
Most of the aforementioned FA models compute a single point estimate of skill levels for each student. Often, however, it is not enough to obtain mere point estimates of students' skill levels; knowing the model's uncertainty in its estimation is crucial because it potentially helps improve the model's performance and improve both students' and instructors' experience with educational systems. For example, in ITS, an underlying model can use the uncertainty information in its estimation of students' skill level to automatically decide that its recommendations based on highly uncertain estimations are unreliable and instead notify a human instructor to evaluate the students' performance. This enables collaboration between ITS and instructors to create a more effective learning environment. In adaptive testing systems [7,71], knowing the uncertainty in model's estimation could help the model intelligently pick the next test items to most effectively reduce its uncertainty about estimated students' skill levels. This will help to potentially reduce the number of items needed to have a confident, accurate estimation of the students' skill mastery level, saving time for both students to take the test and instructors to have a good assessment of the student's skills.
All the above applications require the model to "know what it does not know," i.e., to quantify the uncertainty of its estimation. Achieving this does not necessary changes the model. rather, we need a different inference algorithm for inferring not only a point estimate of student's skill level from observed data (i.e., students' answer records to questions) but also uncertainty in the estimations.
Fortunately, there exist methods that both compute point estimates of students' skill levels and quantifies the uncertainty of those estimates. These methods usually follow the Bayesian inference paradigm, where each student's unknown skill levels are treated as random variables drawn from a posterior distribution. Thus, one can use the credible interval to quantify the model's uncertainty on the student's estimated skill level. A classical method for Bayesian inference is Monte Carlo Markov Chain (MCMC) sampling [17] which has been widely used in many other disciplines [19] other than EDM. In the context of student modeling with FA models, existing works such as [33,15,48,42] have applied MCMC methods to obtain credible intervals.
Unfortunately, classic Bayesian inference methods suffer from extensive computational complexity. For example, each computation step in MCMC involves a time-consuming evaluation of the posterior distribution. Making matters worse, MCMC typically takes many more steps to converge than non-Bayesian inference methods. As a concrete illustration, in [33], the Bayesian inference method (10 minutes) is about 100 times slower than the other non-Bayesian inference method based on stochastic gradient descent (6 seconds). The high computational cost prevents Bayesian inference methods from mass-deployment in large-scale, realtime educational systems where timely feedback is critical for learning [14] and tens of thousands of data points need to be processed in seconds or less instead of minutes or hours. It is thus highly desirable to accelerate Bayesian inference so that one can quantify model's uncertainty in its estimation as efficiently as non-Bayesian inference methods.
Contributions. In this work, we propose VarFA, a novel framework based on variational inference (VI) to perform efficient, scalable Bayesian inference for FA models. The key idea is to approximate the true posterior distribution, whose costly computation slows down Bayesian inference, with a variational distribution. We will see in Section 3 that, with this approximation, we turn Bayesian inference into an optimization problem where we can use the same efficient inference algorithms as in non-Bayesian inference methods. Moreover, the variational distribution is very flexible and we have full control specifying it, allowing us to freely use the latest development in machine learning, e.g., deep neural networks (DNNs), to design the variational distribution that closely approximates the true posterior. Thus, we also regard our work as a first step in applying DNNs to FA models for student modeling, achieving efficient Bayesian inference (enabled by DNNs) without losing interpretability (brough by FA models). We demonstrate the efficacy of our framework on both synthetic and real data sets, showcasing that VarFA substantially accelerates classic Bayesian inference for FA models with no compromise on performance.
The remainder of this paper is organized as follows. Section 2 introduces the problem setup in FA and reviews the important FA models and their inference methods, in particular Bayesian inference. Section 3 explains our VarFA framework in detail. Section 4 presents extensive experimental results that substantiate the claimed advantages of our framework. Section 5 concludes this paper and discusses possible extensions to VarFA.

BACKGROUND AND RELATED WORK
We first set up the problem and review related work. Assume we have a data set Y ∈ R N ×Q organized in matrix format where N is the total number of students and Q is the number of questions. This is a binary students' answer record matrix where each entry yij represents whether student i correctly answered question j. Usually, not all students answer all questions. Thus, Y contains missing values. We use {i, j} ∈ Ω obs to denote entries in Y , i.e., the i-th student's answer record to the j-th question, that are observed.
We are interested in models capable of inferring each i-th student's skill mastery level that can accurately predict the student's answers given the above data. These models are often evaluated on the prediction accuracy and whether the inferred student skill mastery levels are easily interpretable and educationally meaningful. We now review factor analysis models (FA), one of the most widely adopted and successful methodologies for the student modeling task.

Factor Analysis For Student Modeling
One of the earliest FA model for student modeling is based on the item response theory (IRT) [66]. It usually has the following form which assumes that each student's answer yij is independently Bernoulli distributed. The above formula says that the students' answers can be explained by an unknown scalar student skill level factor ci for each student i and an unknown scalar question difficulty level factor µj for each question j. σ(·) is a Sigmoid activation function, i.e., σ(x) = 1/(1 + exp(−x)). The multi-dimensional IRT model (MIRT) [2] extends IRT by using a multi-dimensional vector to represent the student skill levels. More recently, [33] proposed sparse factor analysis model (SPARFA) which extends MIRT by imposing additional assumptions, resulting in improved interpretations of the inferred factors.
Other FA models seek to improve student modeling performance by cleverly incorporate additional auxiliary information. For example, the additive factor model (AFM) [6] incorporates students' accumulative correct answers for a question and the skill tags associated with each question: where κ is the total number of skills in the data, b ik is the i-th student's total number of correct answers for skill k and q kj ∈ {0, 1} indicates whether the skill k is associated with question j. In particular, q kj 's form matrix Q ∈ R K×Q which is commonly known as the Q-matrix in literature [63].
The unknown, to-be-inferred factors are β k , a scalar difficulty factor for each skill k, and ρ k , a scalar learning rate of skill k. The performance factor model (PFM) [50] builds upon AFM that additionally incorporate the total number of a student's incorrect answers to questions associated with a skill. The instructor factor model (IFM) [9] further builds on PFM to incorporate the prior knowledge of whether a student has already mastered a skill. More recently, [68] introduces knowledge tracing machines (KTM) that could flexibly incorporate a number of auxiliary information mentioned above, thus generalizing AFM, PFM and IFM.
Another way to improve student modeling performance is to utilize the so-called memory, i.e., using students' historic action data over time. For example, [35,45] proposed DASH (short for difficulty, ability, and student history) that incorporates a student's total number of correct answers and total number of attempts for a question in a given time window where W is the length of the time windows (e.g., number of days), ρijw and γijw are the i-th student's total number of correct answers and total number of attempts for question j at time w, respectively. θ's are parameters that captures the effect of correct answers and total attempts. More recently, DAS3H [10] extends DASH to further include skill tags associated with each question which essentially amounts to adding, within the last summation term in Eq. 3, another summation over the number of skills.

A general formulation for FA models
Although the above FA models differ in their formulae, modeling assumptions and the available auxiliary data used, we argue that the aforementioned FA models can be unified into a canonical formulation below where ci ∈ R K , mj ∈ R K and µj ∈ R are factors whose dimension, interpretations and subscript indices depend on the specific instantiations of the FA model. To illustrate that Eq. 4 subsumes the FA models mentioned above, we demonstrate, as examples, how to turn SPARFA and AFM into the form in Eq. 4 and how to reinterpret the parameters in the reformulation. The equivalence between Eq. 4 and the original SPARFA formula is immediate and we can interpret the parameters in Eq. 4 as follows to recover SPARFA: K is the number of latent skills that group the actual skill tags in the data set into meaningful coarse clusters; ci represents the i-th student's skill level on the latent skills; mj represents the strength of association of the j-th question with the latent skills which is assumed to be nonnegative and sparse for improved interpretation; and µj represents the difficulty of question j. All three factors in SPARFA are assumed to be unknown.
Similarly, for AFM, we can perform the following change of variables and indices to obtain Eq. 4. We first change the indices mj to mij and µj to µi, and remove the index in ci to c. Then, we simply set c = [β1, ..., βκ, ρ1, ..., ρκ] , mij = [q1j, ..., qκj, q1jb ik , ..., qκjbiκ] and µi = αi to obtain Eq. 4. We can interpret the parameters in the reformulation as follows to recover AFM: κ = K/2 is the total number of skill tags in the data; c represents the unknown skill information including skill difficulty and learning rate; mij summarizes known auxiliary information for the i-th student and j-th question; µi represents the unknown i-th student's skill mastery level. Tthe other FA models can be cast into Eq. 4 in a similar fashion. Note that, if the observed data Y is not binary but rather continuous or categorical, we can simply use a Gaussian or categorical distribution for the observed data and change the Sigmoid activation σ(·) to some other activation functions accordingly. In this way, we still retain the general FA model in Eq. 4. We will use this canonical form to facilitate discussions in the rest of this paper.

Inference Methods for FA Models
We first introduce the inference objective and briefly review maximum log likelihood estimation method. Then, we review existing works that uses Bayesian inference for FA and highlight their high computational complexity, paving the way to VarFA, our proposed efficient Bayesian inference framework based on variational inference.
The factors in Eq. 4 may contain known factors that need no inference. Therefore, for convenience of notation, let [ν, θ, ψ] be a partition of the factors ci, mj and µj where ν contains the unknown students' skill level factors to be estimated, θ contains the remaining unknown factors and ψ contains the known factors. For example, for AFM, .., mQ, µ1, ..., µQ} and ψ = ∅. Further, let the subscripts for these partitions indicate the corresponding latent factors in FA; i.e., for SPARFA, θi = [mi, µi] and νi = ci. Since ψ summarizes known factors, we will omit it from the mathematical expositions in the remainder of the paper.
The inference objective in FA is then to obtain a good estimate of the unknown factors, represented by θ and ν, with respect to some loss function L, usually the marginal data log likelihood. There are two ways to infer the unknown factors.

Maximum Likelihood Estimation
The first way is through the maximum likelihood estimate (MLE) which obtains a point estimate of the unknown factors Recall that Ω obs indicates which student i has provided an answer to question j. R(θ, ν) is a regularization term, e.g., 2 regularization on the factors and λ is a hyper-parameter that controls the strength of the regularization. Most of the works in FA use this method of inference [6,50,9,68,10,35,45]. The advantage of MLE is that the above optimization objective allows the use of fast inference algorithms, in particular stochastic gradient descent (SGD) and its many variants, making the inference highly efficient.

Maximum A Posteriori Estimation
The second way is through maximum a posteriori (MAP) estimation through Bayesian inference. This is the focus of our paper. Recall that we wish to obtain not only a point estimate but also credible interval, which MAP allows while MLE does not.
Bayesian inference methods treat the unknown parameters ν and θ as random variables drawn from some distribution.
We are thus interested in inferring the posterior distribution of ν and θ. Using the Bayes rule, we have that We can similarly derive the posterior for the other unknown factor θ, although in this paper we focus on Bayesian inference for the unknown student skill level parameter ν. In Eq. 6 to Eq. 7 we have applied standard Bayes rule. Note that, because ψ contains known factors and thus does not enter the Bayes rule equation, we have omitted it from the above equations. Once we estimate the posterior distribution for ν from data, we can obtain both point estimates and credible intervals by respectively taking the mean and the standard deviation of the posterior distributions. A number of prior works have proposed to use Bayesian inference for factor analysis, mostly developed for the IRT model [16,44]. More recently, Bayesian methods are developed for more complex models. For example, [33] proposed SPARFA-B, a specialized MCMC algorithm that accounts for SPARFA's sparsity and nonnegativity assumptions.
However, Bayesian inference suffers from high computational cost despite its desired capability to quantify uncertainty. The challenge stems from the difficulty to evaluate the denominator in the posterior distribution in Eq. 6 and 7 which involves an integration over a potentially multidimensional variable. This term usually cannot be computed in close form, thus we usually cannot get an analytical formula for the posterior distribution. Monte Carlo Markov Chain (MCMC) method gets around this difficulty by using analytically tractable proposal distributions to approximate the true posterior and sequentially updating the proposal distribution in a manner reminiscent of gradient descent. MCMC methods also enjoy the theoretical advantage that, when left running for long enough, the proposal distribution is guaranteed to converge to the true posterior distribution [17]. However, in practice, it might take very long time for MCMC to converge to the point that running MCMC is no longer practical when data set is very large. The high computational complexity of MCMC is apparently impractical for educational applications where timely feedback is of critical importance [14].

VARFA: A VARIATIONAL INFERENCE FACTOR ANALYSIS FRAMEWORK
We are ready to introduce VarFA, our variational inference (VI) factor analysis framework for efficient Bayesian learning analytics. The core idea follows the variational principle, i.e., we use a parametric variational distribution to approximate the true posterior distribution. VarFA is highly flexible and efficient, making it suitable for large scale Bayesian inference for FA models in the context of educational data mining.
In this current work, we focus on obtaining credible interval for the student skill mastery factor ν as a first step of VarFA because, recall from Section 1, this factor is often of more practical interest than the other unknown factors. Therefore, currently VarFA is a hybrid MAP and MLE inference method: we perform VI on the unknown student factor ν and MLE on all other unknown factors θ. We consider this hybrid nature of VarFA in its current formulation as a novel feature because classic MLE or MAP estimation are not capable of performing Bayesian inference on a subset of the unknown factors. Extension to VarFA to full Bayesian inference for all unknown factors is part of an ongoing research; see 5 for more discussions.

VarFA Details
Now, we explain in detail how to apply variational inference for FA models for efficient Bayesian inference. Because the posterior distribution is intractable to compute (recall Eq. 6 and 7 and the related discussion), we approximate the true posterior distribution for ν with a parametric variational distribution where φ is a collection of learnable parameters that parametrize the variational distribution and yi is all the answer records by student i. Notably, we have removed the dependency of the variational distribution on ψ and θ so that the variational distribution is solely controlled by the variational parameter φ. Thus, the design of the variational distribution is highly flexible. All we need to do is to specify a class of distributions and design a function parametrized by φ to output the parameters of q φ . Common in prior literature is to use a Gaussian with diagonal covariance for q φ : where its mean and variance [u j , v j ] = f φ (yi). We can use arbitrarily complex functions such as a deep neural network for f φ as long as they are differentiable; more details in Section 3.2. With the above approximation, Bayesian inference turns into an optimization problem under the variational principle, where we now optimize a lower bound, known as the evidence lower bound (ELBO) [5], of the marginal data log likelihood. We derive a novel ELBO objective for variational inference applied to FA models in the EDM setting, summarized in the proposition below.
Proposition 1. In VarFA, the ELBO objective for an FA model in the form of Eq. 4 is where DKL is the KullbackâȂŞLeibler (KL) divergence [39] between two distributions.
Proof. We start with the marginal data log likelihood which is the objective we want to maximize and introduce the ran-dom variables νi's: where step (i) follows from Jensen's inequality [28].
The ELBO objective differs from those used in existing literature in that, in our case, because the data matrix is only partially observed, the summation is only over {i, j} ∈ Ω obs , i.e., the observed entries in the data matrix. In contrast, in prior literature, the summation in the ELBO objective is over all all entries in the data matrix.
We form the following optimization objective to estimate φ and θ: where R(θ) is a regularization term. That is, we perform VI on the student factor ν and MLE inference on the remaining factors θ. More explicitly, to perform VI on ν, we simply compute the mean and standard deviation of q φ using f φ with the learnt variational parameters φ.

Why is VarFA efficient?
VarFA supports efficient stochastic optimization. By formulating Bayesian inference as an optimization problem via the ELBO objective, VarFA allows efficient optimization algorithms such as SGD, with a few additional tricks. In particular, we require all terms in L(θ, φ) to be differentiable with respect to θ and φ which enable gradient computation necessary for SGD. This requirement is easily satisfied with a few moderate assumptions. Specifically, we let the variational distribution q φ and the prior distribution p(νi) for each i to be Gaussian (See Eq. 9; we use a standard Gaussian for the prior distribution p(νi), i.e., p(νi) = N (0, I)) following existing literature [31,53]. These two assumptions are not particularly limiting especially given the vast number of successful applications of VI that rely on the same assumptions [75,40,58,12,8,26,29]. Thanks to the Gaussian assumption, for the first term in L(θ, φ), we can easily compute the KL divergence term analytically in closed form. For the second term in L(θ, φ), we can use the socalled reparametrization trick, i.e., νi = ui + vi , where ∼ N (0, I), which allows a low variance estimation of the gradient of the second term in L(θ, φ) with respect to φ. See Section 2.3 in [31], Section 2.3 in [32] and Section 3 in [53] for more details on stochastic gradient computation in VI. As a result, our framework can be efficiently implemented using a number of open source, automatic differentiation packages such as Tensorflow [1] and PyTorch [49].
VarFA supports amortized inference. Classic Bayesian inference methods such as MCMC infer each parameter νi for each student. As the number of students increases, the number of parameters that MCMC needs to infer also increases, which may not be scalable in large-scale data settings. In contrast, unlike classic Bayesian inference methods, VI is amortized. It does not infer each parameter νi. Rather, it estimates a single set of variational parameter φ responsible for inferring all νi's. In this way, once we have trained FA models using VarFA and obtain the variational parameter φ, we can easily infer νi by simply computing q φ , even for new students, without invoking any additional optimization or inference procedures.

Dealing with missing entries
Note from Eq. 10 that the variational distribution q φ for each νi is conditioned on yi ∈ R Q , i.e., an entire row in Y . However, in practice, the data matrix Y is often only partially observed, i.e., some students only answer a subset of questions. Then, yi's will likely contain missing values which computation cannot be performed on. To work around this issue, we use "zero imputation", a simple strategy that transforms the missing values to 0, following prior work on applying VI in the context of recommender systems [34,55,56,46] that demonstrated the effectiveness of this strategy despite its simplicity. We note that there exist more elaborate ways to deal with missing entries, such as designing specialized function f φ for the variational distribution [37,20,36]. We leave the investigation of more effectively dealing with missing entries to future work.

Remarks
Applicability of VarFA. The VarFA framework is general and flexible and can be applied to a wide array of FA models. The recipe for applying VarFA to an FA model of choice is as follows: 1) formulate the FA model into the canonical formulation as in Eq. 4; 2) partition the factors into student factor ν, the remaining unknown factors θ and known factors ψ; 3) perform VI on ν and MLE estimation on θ, following Section. 3.1.
Relation to variational auto-encoders. Our proposed framework can be regarded as a standard variational autoencoder (VAE) but with the decoder implemented not as a neural network but by the FA model. Because the decoder is constraint to a FA model, it is more interpretable than a neural network.
Relation to other efficient Bayesian factor analysis methods. We acknowledge that VI applied to general FA models have been proposed in existing literature [18,74]. However, to our knowledge, little prior work have applied VI to educational FA models nor investigated its effectiveness. One concurrent work applied VI to IRT [70]. Our VarFA framework applies generally to a number of other FA models by following the recipe in the preceding paragraph, including IRT. Thus, our work complements existing literature in providing promising results in applying VI for FA models in the context of educational data mining.

EXPERIMENTS
We demonstrate the efficacy of VarFA variational inference framework using the sparse factor analysis model (SPARFA) as the underlying FA model. This choice of SPARFA as the FA model to investigate is motivated by its mathematical generality: it assumes no auxiliary information is available and all three factors ci's, mj's and µj's need to be estimated, which makes the inference problem more challenging. [33] provides two inference algorithms including SPARFA-M (based on MLE) and SPARFA-B (based on MAP).
We conduct experiments on both synthetic and real data sets. Using synthetic data sets, we compare VarFA to both SPARFA-M and SPARFA-B. We demonstrate that 1) VarFA predicts students' answers as accurately as SPARFA-M and SPARFA-B; 2) VarFA is almost 100× faster than SPARFA-B. Using real data sets, we compare VarFA to SPARFA-M. We demonstrate that 1) VarFA predicts students' answers more accurately than SPARFA-M; 2) VarFA can output the same insights as SPARFA-M, including point estimate of students' skill levels and questions' associations with skill tags; 3) VarFA can additionally output meaningful uncertainty quantification for student skill levels, which SPARFA-M is incapable of, without sacrifice to computational efficiency. Note that SPARFA-B can also compute uncertainty for small data sets but fails for large data sets due to scalability issues and thus we do not compare to SPARFA-B for real data sets.
Specifically for SPARFA, using the notation convention in Section 2.2, the question-skill association factors mj's and the question difficulty factors µj are collected in θ = {m1, ..., mQ, µ1, ..., µQ}. The student skill level factors cj's are collected in ν = {c1, ..., cN }. Because the factors mj's are unknown, the "skills" in SPARFA are latent (referred to as "latent skills" in subsequent discussions) and are not attached to any specific interpretation. However, as we will see in Section 4.2.4, assuming the skill tags are available as auxiliary information, we can associate the estimated latent skills with the provided skills tags using the same approach proposed in [33]. For the regularization term in Eq. 11, because the factors mj's are assumed to be sparse and nonnegative, we use 1 regularization for mj's. When performing MLE for mj's, we use the proximal gradient algorithm for optimization problem with nonnegative and sparsity requirements; see Section 3.2 in [33] for more details. For the other unknown factor µj's in θ, we ap-  ∼ Exponential(1) when s kj = 1 or setting m kj = 0 when s kj = 0 for each j and k. In all synthetic data experiments, we use π = 0.3 and set the true number of latent skills to K = 5.

Synthetic Data Experiments
Experimental settings. We vary the size of the data set and use 5 different data matrix sizes: 100×50, 300×50, 500×50, 700×50 and 900×50. We select a data missing rate of 50%, i.e., we randomly choose 50% of all entries in the data matrix as training set and the rest as test set. Experi- mental results for each of the above data sizes are averaged over 5 runs where we randomize over the train/test data split. We train SPARFA with both SPARFA-M and VarFA using the Adam optimizer [30] with learning rate = 0.05 for 100 epochs. Regularization hyper-parameters are chosen by grid search for each experiment. For VarFA, we use a simple 3-layer neural network for the function f φ in the variational distribution. We use the hyper-parameter settings for SPARFA-B following Section 4.2 in [33].
Evaluation metrics. We evaluate the inference algorithms on their ability to recover (predict) the missing entries in the data matrix given the observed entries. We term this criterion "student answer prediction". Since our data matrix is binary, using prediction accuracy (ACC) alone does not accurately reflect the performance of the inference algorithms under comparison. Thus, we use area under the receiver operating characteristic curve (AUC) and F1 score in addition to ACC, as standard in evaluating binary predictions [23]. In summary, with VarFA, we obtain posteriors that allow uncertainty quantification with roughly the same computation complexity of computing point estimates and without compromising prediction performance.

Setup
Data sets and pre-processing steps. We perform experiments on three large-scale, publicly available, real educational data sets including ASSISTments 2009-2010 (Assistment) [24], Algebra I 2006-2007 (algebra) [60] and Bridge to Algebra 2006-2007 (bridge) [61,62]. The details of the data sets, including data format and data collection procedure can be found in the preceding references. We remove students and questions that have too few answer records from the data sets to reduce the sparsity of the data sets. Specifically, we keep students and questions that have no less than 30, 35 and 40 answer records for Assistment, Algebra and Bridge data sets, respectively. In each data set, each student may provide more than one answer record for each question. Therefore, we also remove student-question answer records except for the first one. Table 1 presents Experimental settings and evaluation metrics. Optimizer, learning rate, number of epochs, neural network architecture and evaluation metrics are the same as in synthetic data experiments. For real data sets, we use 8 latent skills instead of 5. Other choices of the number of latent skills might result in better performance. However, since we are comparing different inference algorithms for the same model, it is a fair to compare the inference algorithms on the same model with the same number of latent dimensions. We use a 80:20 data split, i.e., we randomly sample 80% and 20% of the observed entries in each data set, without replacement, as training and test sets, respectively. Regularization hyper-parameters are selected by grid search and are different for each data set. We only compare VarFA with SPARFA-M because SPARFA-B does not scale to such large data sets. Table 2 shows the average performance on the test set of each data set comparing VarFA and SPARFA-M for all three data sets and additionally run time. We can see that VarFA achieves slightly better student answer prediction on most data sets and on most metrics. Interestingly, recall that, in the preceding synthetic data set experiment, A possible explanation is that, for synthetic data sets, the underlying data generation process matches the SPARFA model, whereas for real data sets, the underlying data generation process is unknown. VarFA's better performance than SPARFA-M for real data sets implies that VarFA is more robust to the unknown underlying data generation process.

Results: Performance Comparison
As a result, VarFA may be more applicable than SPARFA-M in real-world situations.

Results: Bayesian Inference With VarFA
We now illustrate VarFA's capability of outputting credible intervals using the Assistment data set. Fig. 3 presents violin plots that show the sampled student latent skill levels for a random subset of 10 students. Plots 3a, 3b and 3c shows the inferred students ability for the 3rd, 4th and 7th latent skill dimension. In each plot, the bottom axis shows the student ID and the top axis shows the total number of questions answered by the corresponding student. For each student, the horizontal width of the violin represents the density of the samples; the skinnier the violin, the more widespread the samples are, implying the model's less certainty on its estimations. Fig. 3

Results: Post-Processing for Improved Interpretability
SPARFA assumes that each student factor νi identifies a multi-dimensional skill level on a number of "latent" skills (recall that we use 8 latent skills in our experiments). As mentioned earlier, these latent skills are not interpretable without the aid of additional information. To improve interpretability, [33] proposed that, when the skill tags for each question is available in the data set, we can associate each latent skill with skill tags via a simple matrix factorization. Then, we can compute each students' mastery levels on the actual skill tags. We refer readers to Sections 5.1 and 5.2 in [33] for more technical details. Although the above method to improve interpretability has already been proposed, [33] only presented results on private data sets, whereas here we presents results on publicly available data set with code, making our results more transparent and reproducible.
We again use the Assistment data set for illustration. We compute the association of skill tags in the data set with each of the latent skills and show 4 of the latent skills with their top 3 most strongly associated skill tags. We can see that each latent skill roughly identify the same group of skill tags. For example, latent skill 4 clusters skill tags on statistics and probability while latent skill 7 clusters skill tags on geometry. Thus, by simple post-processing, we obtain an interpretation of the latent skills by associating them with known skill tags Table 3: Illustration of the estimated latent skills with the their top 3 most strongly associated skill tags in the Assistment data set. The percentage in the parenthesis shows the association probability (summed to 1 for each latent skill). We see that the tagged skills associated with each estimated latent skill form intuitive and interpretable groups. in the data.
We can similarly obtain VarFA's estimations of the students' mastery levels on each skill tags through the above process. In Fig. 4, we compare the predicted mastery level for each skill tag (only for the questions this student answered) with the percent of correct answers for that skill tag. Blue curve shows the empirical student's mastery level on a skill tag by computing the percentage of correctly answered questions belonging to a particular skill tag. Orange curve shows VarFA's estimated student mastery level on a skill tag, normalized to range [0, 1]. We can see that, when the student gave more correct answers to questions of a particular skill tag, such as skill tag ID 2, 10, 14 and 30, VarFA also predicts a higher ability score for these skill tags. When the student gave more incorrect answers to questions of a particular skill tag, such as skill tag ID 0, 4 and 27, VarFA also predicts a lower ability score for those skill tags. Although the correspondence is not perfect, VarFA's predicted student's mastery levels match our intuition about student's abilities (using the observed correct answer ratio as an empirical estimate) reasonably well. Thus, by post-processing, we can interpret VarFA's estimated students' latent skill mastery levels using the easily understandable skill tags.

CONCLUSIONS AND FUTURE WORK
We have presented VarFA, a variational inference factor analysis framework to perform efficient Bayesian inference for learning analytics. VarFA is general and can be applied to a wide array of FA models. We have demonstrated the effectiveness of our VarFA using the sparse factor analysis (SPARFA) model as a case study. We have shown that VarFA can very efficiently output interpretable, educationally meaningful information, in particular credible intervals, much faster than classic Bayesian inference methods. Thus, VarFA has potential application in many educational data mining scenarios where efficient credible interval computation is desired, i.e., in adaptive testing and adaptive learning systems. We have also provided open-source code to reproduce our results and facilitate further research efforts.
We outline three possible future research directions and extensions. First, VarFA currently performs Bayesian inference on the student skill mastery level factor. We are working on extending the framework to perform full Bayesian inference for all unknown factors. Extensions to some models Figure 4: Comparison between the estimated skill mastery levels using VarFA's predictions and using empirical observations for student with ID 110. Even though the two curves show different numeric values, they nevertheless demonstrate similar trends, showing that the predictions reasonably match our intuition about student's skill mastery levels.
such as IRT is straightforward, as studied in [70], because all factors can be reasonably assumed to be Gaussian. However, some other models involve additional modeling assumptions, making it challenging to design distributions that satisfy these assumptions. For example, in SPARFA, one of the factors is assumed to be nonnegative and sparse [33]. No standard distribution fulfill these requirements. We are exploring mixture distributions, in particular spike and slab models [27] for Bayesian variable selection, and methods to combine them with variational inference, following [64,65].
Second, other methods exist that accelerate classic Bayesian inference. Recent works in large-scale Bayesian inference proposed approximate MCMC methods that scale to large data sets [73,38,59]. Some of these methods apply variational inference to perform the approximation [13,22]. We are investigating extending VarFA to support approximate MCMC methods.
Finally, the goal of VarFA is to improve real-world educational systems and ultimately improve learning. Therefore, it is necessary to evaluate VarFA beyond synthetic and benchmark data sets. We plan to integrate VarFA into an existing educational system and conduct a case study where students interact with the system in real-time to understand the VarFA's educational implications in the wild.