3.1 Item parameter updates via the truncated average stochastic Newton algorithm

Taking the M2PL model as an example, we first introduce the recursive SNA for estimating item parameters. To align with commonly accepted conventions, we assume that the latent vector $\boldsymbol {\theta }_n$ for examinee $n~(n=1,\ldots ,N)$ follows a multivariate normal distribution with a mean vector of $\boldsymbol {0}$ and an identity matrix as its variance–covariance matrix. When the response data $\boldsymbol {y}_{n}=(y_{n1},\ldots ,y_{nJ})^{\top }$ of the nth examinee are received, the marginal probability for examinee n is defined as follows:

(4) $$ \begin{align} f(\boldsymbol{y}_{n}|\boldsymbol{\eta})=\int_{\boldsymbol{\theta}_n} P(\boldsymbol{y}_{n}|\Theta_n,\boldsymbol{\eta})\phi(\boldsymbol{\theta}_n)d\boldsymbol{\theta}_n, \end{align} $$

where $P(\boldsymbol {y}_{n}|\Theta _n,\boldsymbol {\eta })=\prod _{j=1}^J\pi (\Theta _n^{\top }\boldsymbol {\eta }_j)^{y_{nj}}(1-\pi (\Theta _n^{\top }\boldsymbol {\eta }_j))^{1-y_{nj}}$ . The integral of $f(\boldsymbol {y}_{n}|\boldsymbol {\eta })$ is typically evaluated using an integral approximation procedure. When the dimensionality is five or fewer, the fixed-point Gauss–Hermite quadrature method (Bock & Aitkin, Reference Bock and Aitkin1981) is commonly employed. For higher-dimensional cases, adaptive quadrature techniques are used (Schilling & Bock, Reference Schilling and Bock2005). To simplify the explanation without losing generality, the expression of $f(\boldsymbol {y}_{n}|\boldsymbol {\eta })$ after approximation using the fixed-point Gauss–Hermite quadrature technique is given below:

(5) $$ \begin{align} f(\boldsymbol{y}_{n}|\boldsymbol{\eta})\approx\sum_{k_Q=1}^{K}\dots\sum_{k_1=1}^{K} P(\boldsymbol{y}_{n}|\boldsymbol{X}_k,\boldsymbol{\eta})w(x_{k_1})\dots w(x_{k_Q}). \end{align} $$

Here, $\boldsymbol {X}_k = (1, \boldsymbol {x_k}^{\top })^{\top }$ , $k=1,\ldots ,K$ , where $\boldsymbol {x_k} = (x_{k_1}, \ldots , x_{k_Q})^{\top }$ represents the vector of values at the specified quadrature nodes and K denotes the number of quadrature nodes in a single dimension. Each $x_{k_q}$ denotes the value of a quadrature node in dimension q, and $w(x_{k_q})$ is the weight associated with the quadrature node in dimension q. The weight $w(x_{k_q})$ depends on the height of the normal density function at that node and the spacing between adjacent quadrature nodes. The index $k_q$ represents the position of the K quadrature nodes in dimension q. Traditional offline methods typically rely on batch data to compute the global objective function. In contrast, the SNA updates the parameters using only the current single sample data rather than the entire dataset. Consequently, the objective function (defined in this article as the negative log-likelihood function) for the nth examinee is

(6) $$ \begin{align} \ell(\boldsymbol{y}_{n},\boldsymbol{\eta})=-\log f(\boldsymbol{y}_{n}|\boldsymbol{\eta})\approx-\log\sum_{k=1}^{K^Q} P(\boldsymbol{y}_{n}|\boldsymbol{X}_k,\boldsymbol{\eta})\boldsymbol{w}(\boldsymbol{x}_{k}), \end{align} $$

where $P(\boldsymbol {y}_n|\boldsymbol {X}_k,\boldsymbol {\eta }) = \prod _{j=1}^J\pi (\boldsymbol {X}_k^{\top }\boldsymbol {\eta }_j)^{y_{nj}}(1-\pi (\boldsymbol {X}_k^{\top }\boldsymbol {\eta }_j))^{1-y_{nj}}$ with $\pi (\boldsymbol {X}_k^{\top }\boldsymbol {\eta }_j) = \frac {\exp (\boldsymbol {X}_k^{\top }\boldsymbol {\eta }_j)}{1+\exp (\boldsymbol {X}_k^{\top }\boldsymbol {\eta }_j)}$ . The summation $\sum _{k=1}^{K^Q}=\sum _{k_1=1}^K\dots \sum _{k_Q=1}^K$ , and $\boldsymbol {w}(\boldsymbol {x}_{k}) = w(x_{k_1})\dots w(x_{k_Q})$ is the product of the weights for the respective quadrature nodes.

The SNA updates the parameters incrementally as data arrive. For each newly arrived response data, the item parameters are updated using the current response data and the parameter estimates from the previous iteration. When the response data of the nth examinee arrive, the empirical gradient vector and Hessian matrix (i.e., the first and second derivatives computed with respect to the current sample) for the jth item parameter $(j=1, \ldots , J)$ are derived based on the parameters $\widehat {\boldsymbol {\eta }}_{n-1}$ from the previous iteration, as follows (refer to Chapter 6 of Baker and Kim (Reference Baker and Kim2004) for detailed computations):

(7) $$ \begin{align} \nabla_{j}\ell(\boldsymbol{y}_{n},\widehat{\boldsymbol{\eta}}_{n-1})&=\frac{\partial\ell(\boldsymbol{y}_{n},\boldsymbol{\eta})}{\partial\boldsymbol{\eta}_j }\Big|_{\boldsymbol{\eta}=\widehat{\boldsymbol{\eta}}_{n-1}}=\sum_{k=1}^{K^Q}f_{k}(\widehat{\boldsymbol{\eta}}_{n-1})\left[\pi(\boldsymbol{X}_k^{\top}\widehat{\boldsymbol{\eta}}_{n-1,j})-y_{nj}\right]\boldsymbol{X}_k, \end{align} $$

(8) $$ \begin{align} \nabla_{j}^2\ell(\boldsymbol{y}_{n},\widehat{\boldsymbol{\eta}}_{n-1})&=\frac{\partial^2\ell(\boldsymbol{y}_{n},\boldsymbol{\eta})}{\partial\boldsymbol{\eta}_j\partial\boldsymbol{\eta}_j^{\top} }\Big|_{\boldsymbol{\eta}=\widehat{\boldsymbol{\eta}}_{n-1}} =\sum_{k=1}^{K^Q}f_{k}(\widehat{\boldsymbol{\eta}}_{n-1})\pi(\boldsymbol{X}_k^{\top}\widehat{\boldsymbol{\eta}}_{n-1,j})(1-\pi(\boldsymbol{X}_k^{\top}\widehat{\boldsymbol{\eta}}_{n-1,j}))\boldsymbol{X}_k\boldsymbol{X}_k^{\top}, \end{align} $$

where

(9) $$ \begin{align} f_{k}(\widehat{\boldsymbol{\eta}}_{n-1})=\frac{\prod_{j=1}^J\pi(\boldsymbol{X}_k^{\top}\widehat{\boldsymbol{\eta}}_{n-1,j})^{y_{nj}}\left(1-\pi(\boldsymbol{X}_k^{\top}\widehat{\boldsymbol{\eta}}_{n-1,j})\right)^{1-y_{nj}}\boldsymbol{w}(\boldsymbol{x}_{k})}{\sum_{k=1}^{K^Q}\prod_{j=1}^J\pi(\boldsymbol{X}_k^{\top}\widehat{\boldsymbol{\eta}}_{n-1,j})^{y_{nj}}\left(1-\pi(\boldsymbol{X}_k^{\top}\widehat{\boldsymbol{\eta}}_{n-1,j})\right)^{1-y_{nj}}\boldsymbol{w}(\boldsymbol{x}_{k})}. \end{align} $$

Note that Equation (8) represents the second-order derivative matrix for the newly arrived data and does not approximate the true Hessian matrix. However, for large-scale datasets or real-time data streams, directly computing the true Hessian matrix is computationally expensive. To address this, historical Hessian information is accumulated using a recursive formula, allowing the approximation of the Hessian matrix to be updated incrementally as new data arrive. This approach enables the algorithm to update model parameters efficiently while incorporating new data, without the need to re-process the entire dataset. Leveraging historical information accelerates convergence and enhances the algorithm’s robustness and computational efficiency. The recursive formula for updating the Hessian matrix is provided below for $j = 1, \ldots , J$ :

(10) $$ \begin{align} \boldsymbol{S}_{n,j}&=\boldsymbol{S}_{n-1,j}+\nabla_{j}^2\ell(\boldsymbol{y}_{n},\widehat{\boldsymbol{\eta}}_{n-1})\notag\\ &=\boldsymbol{S}_{n-1,j}+\sum_{k=1}^{K^Q}f_{k}(\widehat{\boldsymbol{\eta}}_{n-1})\pi(\boldsymbol{X}_k^{\top}\widehat{\boldsymbol{\eta}}_{n-1,j})(1-\pi(\boldsymbol{X}_k^{\top}\widehat{\boldsymbol{\eta}}_{n-1,j}))\boldsymbol{X}_k\boldsymbol{X}_k^{\top}\notag\\ &=\boldsymbol{S}_{0,j}+\sum_{i=1}^{n}\sum_{k=1}^{K^Q}f_{k}(\widehat{\boldsymbol{\eta}}_{i-1})\pi(\boldsymbol{X}_k^{\top}\widehat{\boldsymbol{\eta}}_{i-1,j})(1-\pi(\boldsymbol{X}_k^{\top}\widehat{\boldsymbol{\eta}}_{i-1,j}))\boldsymbol{X}_k\boldsymbol{X}_k^{\top}, \end{align} $$

where $\boldsymbol {S}_{0,j}$ represents the initial value of the Hessian matrix, which is symmetric and positive definite, and $\widehat {\boldsymbol {\eta }}_{0} = (\widehat {\boldsymbol {\eta }}_{0,1}^\top , \dots , \widehat {\boldsymbol {\eta }}_{0,J}^\top )^\top $ denotes the initial item parameter estimates. By recursively updating the Hessian matrix estimate, the SNA significantly reduces the computational cost of each iteration, making it more suitable for large-scale datasets and online learning context. The updated formula for the SNA within the IRT framework is provided below:

(11) $$ \begin{align} \widehat{\boldsymbol{\eta}}_{n,j}&=\widehat{\boldsymbol{\eta}}_{n-1,j}-\frac{1}{n}\overline{\boldsymbol{S}}_{n,j}^{-1}\nabla_{j}\ell(\boldsymbol{y}_{n},\widehat{\boldsymbol{\eta}}_{n-1}). \end{align} $$

Here, $\overline {\boldsymbol {S}}_{n,j} = \frac {1}{n+1}\boldsymbol {S}_{n,j}$ . Boyer and Godichon-Baggioni (Reference Boyer and Godichon-Baggioni2023) and Cénac et al. (Reference Cénac, Godichon-Baggioni and Portier2025) argued that the original SNA decreases at a step rate of $\frac {1}{n}$ , which may hinder the dynamics of the algorithm and lead to suboptimal results in cases of improper initialization. To address this issue, we introduce a new step size sequence $\gamma _n = \frac {1}{n^{\gamma }}$ , where $\gamma \in (\frac {1}{2},1)$ , and incorporate an averaging step to ensure the asymptotic efficiency of the algorithm. However, the convergence properties of the parameters for the original SNA cannot be established. Therefore, following the approach of Bercu et al. (Reference Bercu, Godichon and Portier2020), we present a truncated averaged SNA (TASNA). The update process for the jth item parameter ( $j=1,\ldots ,J$ ) is given below:

(12) $$ \begin{align} &f_{k}(\widehat{\boldsymbol{\eta}}_{n-1})=\frac{\prod_{j=1}^J\pi(\boldsymbol{X}_k^{\top}\widehat{\boldsymbol{\eta}}_{n-1,j})^{y_{nj}}\left(1-\pi(\boldsymbol{X}_k^{\top}\widehat{\boldsymbol{\eta}}_{n-1,j})\right)^{1-y_{nj}}\boldsymbol{w}(\boldsymbol{x}_{k})}{\sum_{k=1}^{K^Q}\prod_{j=1}^J\pi(\boldsymbol{X}_k^{\top}\widehat{\boldsymbol{\eta}}_{n-1,j})^{y_{nj}}\left(1-\pi(\boldsymbol{X}_k^{\top}\widehat{\boldsymbol{\eta}}_{n-1,j})\right)^{1-y_{nj}}\boldsymbol{w}(\boldsymbol{x}_{k})}, \end{align} $$

(13) $$ \begin{align} &\boldsymbol{Z}_{n,j} =\sum_{k}f_{k}(\widetilde{\boldsymbol{\eta}}_{n-1})\left[\pi(X_k^{\top}\widetilde{\boldsymbol{\eta}}_{n-1,j})-y_{nj}\right]\boldsymbol{X}_k, \end{align} $$

(14) $$ \begin{align} &\widetilde{\boldsymbol{\eta}}_{n,j}=\widetilde{\boldsymbol{\eta}}_{n-1,j}-\gamma_{n}\overline{\boldsymbol{S}}_{n-1,j}^{-1}\boldsymbol{Z}_{n,j} \end{align} $$

(15) $$ \begin{align} &\widehat{\boldsymbol{\eta}}_{n,j}=(1-\kappa_{n})\widehat{\boldsymbol{\eta}}_{n-1,j}+\kappa_{n}\widetilde{\boldsymbol{\eta}}_{n,j}, \end{align} $$

(16) $$ \begin{align} &\boldsymbol{S}_{n,j}=\boldsymbol{S}_{n-1,j}+\sum_{k=1}^{K^Q}f_{k}(\widehat{\boldsymbol{\eta}}_{n-1})\alpha_{n,k}\boldsymbol{X}_k\boldsymbol{X}_k^{\top}, \end{align} $$

where $\widetilde {\boldsymbol {\eta }}_{n,j}$ and $\widehat {\boldsymbol {\eta }}_{n,j}$ represent the jth item parameter estimates before and after applying the averaging step, respectively. The step sequence is defined as $\gamma _{n} = \frac {1}{(n + c_{\eta })^{\gamma }}$ , where the step parameter $\gamma \in (\frac {1}{2}, 1)$ . The averaging step is defined as

$$ \begin{align*}\widehat{\boldsymbol{\eta}}_{n,j} = (1-\kappa_{n})\widehat{\boldsymbol{\eta}}_{n-1,j} + \kappa_{n}\widetilde{\boldsymbol{\eta}}_{n,j} = \frac{1}{n}\sum_{i=1}^{n}\widetilde{\boldsymbol{\eta}}_{i,j}, \end{align*} $$

with the averaging sequence $\kappa _{n} = \frac {1}{n}$ . The term $\boldsymbol {Z}_{n,j}$ denotes the updated gradient $\nabla _{j}\ell (\boldsymbol {y}_{n},\widetilde {\boldsymbol {\eta }}_{n-1})$ , and $\overline {\boldsymbol {S}}_{n-1,j} = \frac {1}{n}\boldsymbol {S}_{n-1,j}$ represents the scaled Hessian matrix from the previous iteration. It is important to note that the update of $\boldsymbol {S}_{n,j}$ depends on the averaged parameter $\widehat {\boldsymbol {\eta }}_{n-1,j}$ rather than $\widetilde {\boldsymbol {\eta }}_{n-1,j}$ , enabling faster convergence. The value of $\alpha _{n,k}$ is given by

$$ \begin{align*}\alpha_{n,k} = \max\left\{\pi(\boldsymbol{X}_k^{\top}\widehat{\boldsymbol{\eta}}_{n-1,j})(1-\pi(\boldsymbol{X}_k^{\top}\widehat{\boldsymbol{\eta}}_{n-1,j})), \frac{c_{\beta}}{n^{\beta}}\right\}, \end{align*} $$

which is used for truncating the Hessian matrix.

Remark 3. For some positive constant $c_{\beta }$ , the sequence $\alpha _{n,k}$ is defined as

(17) $$ \begin{align} \alpha_{n,k}&=\max\{\pi(\boldsymbol{X}_k^{\top}\widehat{\boldsymbol{\eta}}_{n-1,j})(1-\pi(\boldsymbol{X}_k^{\top}\widehat{\boldsymbol{\eta}}_{n-1,j})),\frac{c_{\beta}}{n^{\beta}}\}\notag\\ &=\max\left\{\frac{1}{4(\cosh (\boldsymbol{X}_k^{\top}\widehat{\boldsymbol{\eta}}_{n-1,j}/2))^2},\frac{c_{\beta}}{n^{\beta}}\right\}, \end{align} $$

where $\beta < \gamma - \frac {1}{2}$ . Since the hyperbolic cosine function $\cosh (x)$ takes values in the range $[1, +\infty )$ , we assume that $c_{\beta } < \frac {1}{4}$ , which ensures $\alpha _{n,k} < \frac {1}{4}$ for $n \geq 1$ . Under this truncation operation, the convergence of the Hessian matrix estimation can be proven. Please see Section 4 for details.

3.2 Updating ability parameters by expected a posteriori estimation

For estimating the examinees’ ability parameters, researchers typically rely on methods, such as maximum likelihood estimation (MLE), maximum a posteriori estimation (MAP), or expected a posteriori estimation (EAP) after the item parameters have been estimated (Bock & Aitkin, Reference Bock and Aitkin1981; Wang, Reference Wang2015). It is well known that MLE relies on batch data to compute the global likelihood, which indicates the entire likelihood function must be recalculated whenever new data are added. This makes MLE computationally expensive in online environments where data continuously stream in, rendering it unsuitable for efficient online updating. In contrast, MAP incorporates prior information by maximizing the posterior probability. Although it is more flexible than MLE due to its ability to include prior information, it typically requires a large amount of historical data for global optimization. Furthermore, because of the structure of IRT models, MAP estimation often relies on iterative optimization methods, such as Newton’s method or gradient descent, to maximize the posterior probability. This iterative process further increases computational complexity and time. Detailed implementations of these two methods can be found in Wang (Reference Wang2015) and Bock and Aitkin (Reference Bock and Aitkin1981).

Unlike the two methods mentioned above, EAP estimation computes the posterior expectation of the ability parameter $\boldsymbol {\theta }$ as a point estimate. The process for EAP estimation of $\boldsymbol {\theta }$ has been outlined by Bock and Aitkin (Reference Bock and Aitkin1981), Muraki and Engelhard Jr. (Reference Muraki and Engelhard1985), and Wang (Reference Wang2015). EAP estimation has the advantage of ensuring that the estimated ability parameter derived from item scores, even with short test items, does not diverge toward infinity. This makes it particularly suitable for cases where data are limited or strong prior information is available (Reckase, Reference Reckase2009). The mathematical expression for this estimation is as follows:

(18) $$ \begin{align} \widehat{\theta}_{nq}=\mathbb{E}(\theta_{nq}|\boldsymbol{y}_n)=\frac{\int\theta_{nq}P(\boldsymbol{y}_n|\Theta_n,\boldsymbol{\eta})\phi(\boldsymbol{\theta}_n)d\boldsymbol{\theta}_n}{\int P(\boldsymbol{y}_n|\Theta_n,\boldsymbol{\eta})\phi(\boldsymbol{\theta}_n)d\boldsymbol{\theta}_n}=\int\theta_{nq}\frac{P(\boldsymbol{y}_n|\Theta_n,\boldsymbol{\eta})\phi(\boldsymbol{\theta}_n)}{f(\boldsymbol{y}_n|\boldsymbol{\eta})}d\boldsymbol{\theta}_n. \end{align} $$

Here, $\theta _{nq}$ represents the qth dimension of the latent trait for the nth examinee, where $q = 1, \ldots , Q$ . Given known item parameters $\boldsymbol {\eta }$ , the integral in Equation (18) can be approximated using the Gauss–Hermite quadrature method. The approximated expression for $\widehat {\theta }_{nq}$ is

(19) $$ \begin{align} \widehat{\theta}_{nq} \approx \frac{\sum_{k=1}^{K^Q}x_{k_q}P(\boldsymbol{y}_{n}|\boldsymbol{X}_k,\boldsymbol{\eta})\boldsymbol{w}(\boldsymbol{x}_{k})}{\sum_{k=1}^{K^Q}P(\boldsymbol{y}_{n}|\boldsymbol{X}_k,\boldsymbol{\eta}) \boldsymbol{w}(\boldsymbol{x}_{k})}. \end{align} $$

Therefore, $\widehat {\theta }_{nq}$ depends not only on the quadrature node vector and its weights but also on the data $\boldsymbol {y}_{n}$ and the item parameters $\boldsymbol {\eta }$ . In other words, when the response data for the nth examinee arrive, the examinee’s latent traits can be quickly updated using EAP estimation, provided that the item parameters have already been updated. Given the real-time estimates of the item parameters $\widehat {\boldsymbol {\eta }}_{n-1}$ , the formula for updating the ability parameters is given by

(20) $$ \begin{align} \widehat{\boldsymbol{\theta}}_{n} = \sum_{k=1}^{K^Q}\boldsymbol{x}_{k}f_{k}(\widehat{\boldsymbol{\eta}}_{n-1}). \end{align} $$

3.3 Implementation of stochastic Newton online estimation algorithms

The specific algorithmic flow of the online parameter estimation method in the IRT framework is outlined in Algorithm 1.

5.1 Simulation designs

Five factors were manipulated to vary different simulation conditions: (1) The number of examinees, i.e., $n = (5,000, 10,000, 20,000)$ . (2) The number of items, i.e., $J = (20, 40)$ . (3) The dimensions of latent traits, i.e., $Q = (1, 2, 3,4)$ . $Q = 1$ represents the unidimensional 2PL model. (4) The number of Gauss–Hermite quadrature nodes, i.e., $K = (10, 20)$ . The number of Gauss–Hermite quadrature nodes significantly affects the accuracy of numerical integration and the convergence of the algorithm. By simulating and comparing model performance at different values of K, we aim to identify an appropriate number of nodes that balances accuracy and computational efficiency. (5) Different step sizes, i.e., $\gamma = (0.65, 0.75, 0.9, 1)$ . An appropriate step size can accelerate the algorithm’s convergence while ensuring its stability and enhancing the accuracy of the final solution. By comparing different step size settings, the optimal step size strategy can be determined. Each simulation condition was replicated 50 times.

5.2 Data generation and model identifiability conditions

We generated item response data for the 2PL model and M2PL model based on Equation (1) and Equation (2), respectively. To simplify the computation and enhance the interpretability of the model, we assumed that the latent ability dimensions are independent. This assumption also eliminates the indeterminacy of scale, which is a common practice in other research (e.g., Cui et al., Reference Cui, Wang and Xu2024). Therefore, the ability parameters $\theta _{iq}$ for each examinee i were sampled from a standard normal distribution $N(0, 1)$ , where $q = 1, \ldots , Q$ . Following Shang et al. (Reference Shang, Xu, Shan, Tang and Ho2023), for each item j and dimension q, the discrimination parameters $a_{jq}$ were generated from the uniform distribution $U(0.5, 2)$ , while the difficulty parameters $b_j$ (for the 2PL model) and the intercept parameters $d_j$ (for the M2PL model) were generated from the standard normal distribution $N(0, 1)$ , for $j = 1, \ldots , J$ and $q = 1, \ldots , Q$ .

To ensure the identifiability of the 2PL model, we imposed constraints on the ability parameters by assuming that the population distribution of abilities follows a standard normal distribution. For the M2PL model, we applied the second constraint method proposed by Béguin and Glas (Reference Béguin and Glas2001). This method ensures identifiability by constraining the item parameters as follows: the first Q items’ intercept parameters $d_1, \ldots , d_Q$ are set to zero, and the discrimination parameter $a_{jq}$ is set to 1 when $j = q$ and to 0 when $j \neq q$ , for each $j = 1, \ldots , Q$ and $q = 1, \ldots , Q$ . Meanwhile, the ability parameters are freely estimated.

5.3 Parameter settings and initialization values

In the SNA, several parameters require adjustment, including $c_{\eta }$ , $c_{\beta }$ , $\beta $ , and the step size parameter $\gamma $ . Among these, the step size parameter is the most critical, as it directly controls the magnitude of each update and significantly influences the algorithm’s convergence speed and stability. A step size that is too large may lead to unstable convergence or even divergence, while a step size that is too small can result in slow convergence or becoming trapped in a local optimum. Clearly, the choice of step size impacts both the update direction and magnitude, thereby affecting the accuracy of the final solution and the speed of convergence. To simplify the parameter tuning process and reduce complexity, most parameters were fixed, allowing the primary focus to remain on step size adjustment, which enhances the efficiency and practicality of simulations. Specifically, following Boyer and Godichon-Baggioni (Reference Boyer and Godichon-Baggioni2023), we set $c_{\eta } = 100$ . Based on Bercu et al. (Reference Bercu, Godichon and Portier2020), we set $c_{\beta } = 10^{-10}$ and $\beta = \gamma - \frac {1}{2} - 0.01$ . Additionally, we conducted simulation studies with different values of the tuning parameters ( $c_{\eta }, c_{\beta }, \beta $ ) to systematically assess their impact on the algorithm’s performance. The results indicate that the current configuration of these parameters performs well in maintaining stable and satisfactory estimation accuracy and convergence speed, thereby validating their rationality and robustness. The detailed simulation results are provided in Section A.2 of the Supplementary Material.

In addition, we initialized all estimates of the discrimination parameter to 1 and the difficulty/intercept parameter to 0. This approach simplifies the algorithm’s initialization process and helps address convergence difficulties caused by excessively high or low initial estimates during iterations. To take advantage of the numerical efficiency of the Gauss–Hermite quadrature method for approximate integration, we set the initial Hessian matrix value as

$$ \begin{align*}\boldsymbol{S}_{0,j} = \frac{1}{K^Q}\sum_{k=1}^{K^Q}\boldsymbol{X}_k \boldsymbol{X}_k^{\top}=\text{diag}\left(1,\frac{1}{K}\sum_{k_1=1}^{K}x_{k_1}^2,\dots,\frac{1}{K}\sum_{k_Q=1}^{K}x_{k_Q}^2\right), \quad \text{for } j=1,\ldots,J. \end{align*} $$

This approach addresses common numerical instability issues and provides a robust foundation for subsequent iterations of the algorithm. While this initialization was chosen, it is equally feasible to initialize the Hessian matrix as the identity matrix, provided it ensures numerical stability and avoids singular matrices, which may lead to computational failures. It is worth noting that, in both the simulation and the subsequent empirical study, the batch size used in our algorithm consistently corresponds to the sample size of a single examinee, with parameter updates occurring after observing each individual’s response data.

5.4 Evaluation criteria

To evaluate parameter recovery, the bias and RMSE of the item parameter estimates and ability parameter estimates are calculated. Computational efficiency is assessed based on the average running time of each replication. The following formulas were calculated:

(34) $$ \begin{align} &\text{Bias}(\boldsymbol{\eta})=\frac{1}{R}\sum_{r=1}^R(\widehat{\boldsymbol{\eta}}_r-\boldsymbol{\eta}^*),~~~~\text{RMSE}(\boldsymbol{\eta})=\sqrt{\frac{1}{R}\sum_{r=1}^R(\widehat{\boldsymbol{\eta}}_r-\boldsymbol{\eta}^*)^2},~~~~T=\frac{1}{R}\sum_{r=1}^RT_{r}, \end{align} $$

where R represents the number of replications, $\widehat {\boldsymbol {\eta }}_r$ denotes the parameter estimate from the rth replication, $\boldsymbol {\eta }^*$ indicates the true value of the parameter, and $T_r$ represents the running time of the rth replication. The estimation of ability parameters is carried out as follows: first, the item parameters estimated using the mirt package are treated as known values. Then, the EAP method is applied to estimate $\boldsymbol {\theta }$ . These estimates are subsequently compared with the real-time ability parameter estimates obtained from the two online algorithms.

5.5 Result analysis

5.5.1 Analysis of parameter recovery results

Tables 1–3 present the RMSEs of the model parameter estimates for three different latent dimensions when the number of items is $J=20$ . The results for $J=40$ and $Q=4$ can be found in Tables A1–A4 in the Supplementary Material. From Tables 1–3 and Tables A1–A4 in the Supplementary Material, we can observe the impact of different factors on parameter estimation.

Table 1 Average RMSEs for the 2PL model parameters with $J=20$ across TSNA, TASNA, and EM algorithm implemented in the mirt package

Note: Boldface values indicate estimation results with relatively smaller RMSEs under the same simulation condition. In some settings, more than one tuning parameter value shows comparable performance and is therefore highlighted.

Table 2 Average RMSEs for the M2PL model parameters with $J=20$ , $Q=2$ across TSNA, TASNA, and EM algorithm implemented in the mirt package

Note: Boldface values indicate estimation results with relatively smaller RMSEs under the same simulation condition. In some settings, more than one tuning parameter value shows comparable performance and is therefore highlighted.

Table 3 Average RMSEs for the M2PL model parameters with $J=20$ , $Q=3$ across TSNA, TASNA, and EM algorithm implemented in the mirt package

Note: Boldface values indicate estimation results with relatively smaller RMSEs under the same simulation condition. In some settings, more than one tuning parameter value shows comparable performance and is therefore highlighted.

First, as the sample size increases, the two proposed online algorithms and the EM algorithm show improved parameter estimation accuracy across different latent dimensions. However, as the latent dimensionality increases, the performance of the EM algorithm deteriorates significantly, while the two online algorithms (TASNA and TSNA) exhibit relatively stable performance under various settings. Specifically, for latent dimensions $Q = 1, 2, 3$ , the EM algorithm with fixed quadrature constraints achieves the best overall accuracy, and both TASNA with $\gamma = 0.65, 0.75$ and TSNA with $\gamma = 0.75, 0.9$ yield comparable estimation results. However, as shown in Tables A3 and A4 in the Supplementary Material, when the latent dimension increases to $Q = 3$ and $Q = 4$ , the EM algorithm’s RMSE for parameter estimation exceeds 0.1 in the case of $J = 40$ items, while TASNA with $\gamma = 0.65, 0.75$ and TSNA with $\gamma = 0.75, 0.9$ maintain lower RMSE values, significantly outperforming the EM algorithm.

Secondly, for fixed sample sizes and item numbers, the optimal value of $\gamma $ largely depends on the specific online algorithm used. In the TASNA algorithm, $\gamma = 0.75$ and $\gamma = 0.65$ yield better outcomes, with $\gamma = 0.75$ generally outperforming $\gamma = 0.65$ across various conditions. Under these configurations, TASNA not only achieves estimation accuracy comparable to the EM algorithm, but also converges faster and more stably, even with relatively small sample sizes. In contrast, for the TSNA algorithm, $\gamma = 0.9$ and $\gamma = 0.75$ typically lead to better estimation results, with $\gamma = 0.9$ performing best in most scenarios. Overall, TASNA tends to achieve higher estimation accuracy than TSNA. As the sample size increases, TSNA with $\gamma = 0.75,0.9$ gradually reaches estimation performance comparable to TASNA with $\gamma = 0.65,0.75$ and the EM algorithm. The poorest performance is observed for TASNA when $\gamma = 1$ , with noticeably higher estimation errors and slower convergence. This result aligns with the theoretical analysis, which suggests that TASNA performs better in terms of stability and convergence when $\gamma < 1$ , while setting $\gamma = 1$ may lead to excessive variance accumulation and degraded estimation accuracy.

Moreover, given a fixed sample size, the effectiveness of the quadrature node setting (K) is primarily influenced by the number of items. Simulation results show that when the number of items is small (e.g., $J = 20$ ), the estimation accuracy achieved with $K = 10$ and $K = 20$ is nearly identical. This suggests that, with a relatively small number of items, increasing the number of quadrature nodes does not significantly improve estimation performance and instead results in additional computational cost. As the number of items increases, the number of parameters to be estimated also grows, and in such cases, increasing K can enhance estimation accuracy. However, the improvement from $K = 10$ to $K = 20$ is considerably smaller than the improvement achieved by increasing the sample size. Furthermore, we observed that in most settings, the accuracy obtained with $K = 10$ remains within an acceptable range. Therefore, from the perspective of balancing accuracy and computational efficiency, $K = 10$ is generally the most cost-effective choice.

In addition, regarding the recovery of the ability parameter $\boldsymbol {\theta }$ , both stochastic algorithms show very similar estimation results across different $\gamma $ settings, with the estimates close to the EAP estimation of $\boldsymbol {\theta }$ , as presented in Tables 1 and 2 and Tables A1, A2, A6, and A7 in the Supplementary Material. The recovery results for each method improve as the number of items J increases. Finally, it is worth noting that, under the four-dimensional ( $Q = 4$ ) setting, due to increased model complexity and the “curse of dimensionality,” the computational burden of numerical integration rises sharply when the number of Gauss–Hermite quadrature nodes is set to $K = 20$ . This results in a substantial increase in total running time, even surpassing that of the mirt package under the same conditions. Considering both simulation efficiency and computational resource constraints, we ultimately decided not to run simulations for $K = 20$ in this setting to avoid unnecessary computational costs. However, from Tables A4 and A5 in the Supplementary Material, we found that with $K = 10$ , the proposed algorithms not only demonstrate desirable computational efficiency but also achieve significantly better estimation accuracy compared to mirt. Therefore, even in a four-dimensional setting, the proposed methods remain highly applicable and competitive.

5.5.2 Real-time evaluation of parameter estimation accuracy

Figures 1 and 2, respectively, show the real-time estimation results of item parameters for the 2PL model and the M2PL model with $Q=2$ under the conditions of $N = 20{,}000$ , $J = 20$ , and $K = 10$ . Considering that the EM algorithm is not an online method, to facilitate comparison between the proposed method and the traditional EM algorithm across various sample sizes, the figure presents the estimation results of the EM algorithm at multiple discrete sample points ( $N =$ 500, 1,000, 2,000, 2,500, 5,000, 7,500, 10,000, 15,000, and 20,000). The top two subplots in Figures 1 and 2 display the RMSEs of the real-time estimates for the item parameters. The plots reveal that as the sample size increases, the RMSEs for all $\gamma $ settings decrease and converge rapidly. For example, when the sample size is 2,500, the RMSE of the TASNA with $\gamma = 0.65$ and $\gamma = 0.75$ has already decreased to approximately 0.1. When the sample size reaches 5,000, the RMSE further reduces to a range of around 0.050 to 0.075. Additionally, the TASNA demonstrates smoother and more stable performance compared to the TSNA, which lacks an averaging step and exhibits overall higher fluctuations, especially with smaller $\gamma $ . Furthermore, the TASNA with $\gamma = 0.65$ and $\gamma = 0.75$ performs the best, closely aligning with those from the EM algorithm, which is consistent with the findings presented in Tables 1–3. This is followed by the TSNA with $\gamma = 0.9$ and $\gamma = 0.75$ . It is worth noting that, in small sample scenarios (e.g., $N =$ 500, 1,000, and 2,000), the TASNA method with $\gamma = 0.65$ and $\gamma = 0.75$ achieves accuracy nearly equivalent to that of the EM algorithm. In contrast, the TASNA with $\gamma = 1$ yields the worst performance for the discrimination parameters, while the TSNA with $\gamma = 0.65$ shows the worst performance for the difficulty/intercept parameters, exhibiting substantial instability. The lower subplots in Figures 1 and 2 present the bias of the item parameter estimates, which align with the RMSE results. The bias of TASNA with $\gamma = 0.65$ and $\gamma = 0.75$ remains very stable and close to zero, whereas the bias of TSNA under the same $\gamma $ settings fluctuates more significantly. These findings highlight the advantages of TASNA, particularly with $\gamma = 0.65$ and $\gamma = 0.75$ , in achieving accurate and stable parameter estimates under the given conditions.

Figure 1 Bias and RMSE of real-time estimates of item parameters in the 2PL model across TSNA, TASNA, and EM algorithm for different step sizes with $N=20{,}000$ , $J=20$ , and $K=10$ .

Figure 2 Bias and RMSE of real-time estimates of item parameters in the M2PL model with $Q=2$ across TSNA, TASNA, and EM algorithm for different step sizes with $N=20{,}000$ , $J=20$ , and $K=10$ .

Additionally, the real-time estimation results for the item parameters under other conditions are provided in Figures A1–A6 in the Supplementary Material. These results are consistent with those in Figures 1 and 2, further confirming the robustness and generalizability of the findings across different conditions.

5.5.3 Comparative analysis of running time

Figure 3 presents the average running time of the three estimation methods for all conditions with a fixed $\gamma = 0.75$ . The results clearly indicate that the computation time of the proposed online algorithms is significantly faster than the EM algorithm implemented in the mirt package. Among the online algorithms, TSNA is slightly faster than TASNA, as it does not include the averaging step, while the EM algorithm is the slowest overall. Furthermore, the number of quadrature nodes (K) affects the computation time of the online algorithm: increasing the number of nodes leads to slower computation time, which is particularly noticeable in the M2PL model ( $Q = 2$ and $Q=3$ ). However, even with a higher number of nodes, the online algorithm remains faster than the mirt package. Therefore, our method consistently maintains short computation time with acceptable estimation error for small and medium sample sizes, and demonstrates a smoother, more predictable increase in running time. As sample sizes and item numbers grow, and as the latent trait dimension increases to $Q = 2$ and $Q = 3$ , its advantages in computational efficiency and stability become even more pronounced. In contrast, the computational cost of the mirt method rises significantly with larger sample sizes and higher dimensionality, particularly in high-dimensional contexts.

Figure 3 Running time for the TSNA, TASNA, and EM algorithms across 50 simulation replications, displayed as median values with error bars representing the 25th and 75th percentiles.

For instance, under the three-dimensional setting, with $N = 20{,}000$ and $J = 20$ , the computation time for both online algorithms is approximately 10 seconds when $K = 10$ , whereas the mirt package takes nearly 200 seconds under the same conditions. Although running time increases when using $K = 20$ , it remains shorter than that of mirt. Given that estimation accuracy is nearly identical for $K = 10$ and $K = 20$ in this setting, using $K = 10$ offers a clear advantage in both accuracy and efficiency when $J = 20$ . This trend becomes even more pronounced as the number of items increases to $J = 40$ : with $K = 10$ , the online algorithms take about 20 seconds, while the mirt package requires around 1,000 seconds under the same conditions. Even with $K = 20$ , the online methods take only about 300 seconds—still much shorter than mirt. Therefore, TSNA and TASNA demonstrate superior computational stability and scalability as sample size and model complexity increase. These results highlight the remarkable computational efficiency of the proposed online algorithms, particularly in large-scale applications. Please refer to Table A5 in the Supplementary Material for detailed average running time.

6.1 Data description

In this section, real data from the TIMSS were analyzed to evaluate the performance of the proposed TASNA. TIMSS (Martin et al., Reference Martin, von Davier and Mullis2020; Mullis & Martin, Reference Mullis and Martin2017) is a well-established international assessment of mathematics and science achievement for fourth- and eighth-grade students conducted every four years. The math and science item pools at each grade level are divided into item blocks, which are assembled into 14 distinct student achievement booklets. Each student completes one booklet. We selected items identical to those used by Cui et al. (Reference Cui, Wang and Xu2024), which include one mathematics block and one science block from the eighth-grade Booklets 5 and 6. These blocks contain a total of 28 items, comprising 13 mathematics items and 15 science items. The specific item codes are provided in Appendix C of Cui et al. (Reference Cui, Wang and Xu2024). After excluding students with missing responses, the final sample size was 9,874. The mathematics block contained two polytomous items, while the science block included three polytomous items, all with a maximum score of 2. To simplify the analysis, the polytomous items were treated as dichotomous by recoding scores of 2 as correct (i.e., 1) and scores below 2 as incorrect (i.e., 0). This recoding ensures consistency with the scoring framework used for dichotomous items in the analysis.

6.2 Analysis of empirical data results fitted with the unidimensional 2PL model

In the operational analysis of TIMSS data, the unidimensional 2PL IRT model was applied separately to the mathematics and science domains. Cui et al. (Reference Cui, Wang and Xu2024) used their proposed modified information criterion and demonstrated that both domains are essentially unidimensional when analyzed independently. Thus, we first present the results under a unidimensional framework. Consistent with the simulation studies, the estimation methods include the EM algorithm with fixed quadrature (implemented in the mirt package), the TASNA, and the TSNA. Previous simulation results indicate that when the number of items is fewer than 40, TASNA with step size parameters $\gamma = 0.65$ and $0.75$ , as well as TSNA with $\gamma = 0.75$ and $0.9$ , provides better parameter recovery. While the performance of these methods is comparable when the number of quadrature nodes is $K=10$ or $K=20$ , the computation time is considerably shorter for $K=10$ . Therefore, we fixed $K=10$ and selected TASNA with $\gamma = 0.65, 0.75$ and TSNA with $\gamma = 0.75, 0.9$ . All other parameter settings are consistent with those used in the simulation studies. Furthermore, to address scale indeterminacy, we assumed that the covariance matrix of the latent traits is an identity matrix, following the approach used by Cui et al. (Reference Cui, Wang and Xu2024). This assumption standardizes the scale of the latent variables, ensuring consistency across analyses.

Figure 4 shows the parameter estimates for discrimination ( $\boldsymbol {a}$ ) and difficulty ( $\boldsymbol {b}$ ) obtained using the EM algorithm and the SNAs under the 2PL model. The results indicate that the estimates from both methods are quite similar. However, for some items, the EM algorithm tends to produce slightly higher estimates than the SNAs. Notably, the TSNA method with $\gamma = 0.75$ is more sensitive to extreme parameter values. This happens because, as parameter values become very large, the gradient of the objective function can also grow disproportionately. With $\gamma = 0.75$ , the step size also becomes relatively large, causing each update to be large magnitude. This, in turn, can lead to “jumps” in the parameter estimates that may overshoot the true optimal point of the objective function, resulting in inflated parameter estimates. In contrast, TASNA reduces the influence of extreme outliers and noise by averaging information across all iterations. This results in more stable updates, helping the algorithm converge more consistently toward the optimal solution, and ultimately providing more accurate and reliable parameter estimates across different conditions.

Figure 4 Estimates of item parameters for the 2PL model across TSNA, TASNA, and EM algorithm with corresponding step sizes.

To demonstrate the real-time estimation performance of the algorithm, we present the real-time trajectory plots of discrimination parameter estimation in Figure D1 in the Supplementary Material and difficulty parameter estimation in Figure D2 in the Supplementary Material. These plots show that the parameter estimates obtained using the TASNA demonstrate more stable convergence, with no extreme values, even in the initial stages of estimation. In contrast, Figure D2 in the Supplementary Material shows that the TSNA method exhibits more noticeable fluctuations, particularly for some items in the difficulty parameter $\boldsymbol {b}$ (e.g., ME62002, SE62284, and SE62173B), where unusually large estimates are observed. These findings highlight the necessity of including the averaging step in the algorithm. For large-scale educational testing scenarios that require real-time estimation, the TASNA is notably more stable and accurate compared to the TSNA. Moreover, TASNA effectively avoids extreme parameter estimates, making it a more practical and reliable choice for real-time applications.

6.3 Analysis of empirical data results fitted with the multidimensional 2PL model

To investigate the relationship between latent abilities in different domains, we adopted the M2PL model to estimate all 28 items in the TIMSS dataset, assuming two latent dimensions ( $Q=2$ ). The mathematics proficiency and science proficiency are conceptually distinct, each measuring different cognitive skills. Mathematics proficiency primarily involves numerical reasoning, algebraic thinking, and problem-solving skills, while science proficiency focuses on understanding scientific principles, experimental design, and analytical reasoning. These differences suggest that mathematics and science abilities should be treated as separate latent traits. The modified information criterion introduced by Cui et al. (Reference Cui, Wang and Xu2024) further supports this distinction, showing that a two-dimensional latent structure provides a significantly better fit compared to a unidimensional model. This finding aligns with the inherent design characteristics of the TIMSS data, where the items are specifically developed to assess the distinct domains of mathematics and science proficiency. To ensure model identifiability, we adopt the constraints suggested by Béguin and Glas (Reference Béguin and Glas2001). For the first item in the math domain, we constrain the discrimination parameters as $a_{11} = 1$ and $a_{12} = 0$ , with the corresponding intercept parameter set to $d_1 = 0$ . Similarly, for the first item in the science domain, we constrain the discrimination parameters as $a_{14,1} = 0$ and $a_{14,2} = 1$ , with the corresponding intercept parameter set to $d_{14} = 0$ .

Figure 5 presents the parameter estimation results for both the SNA and the EM algorithm in the M2PL model. As observed in Cui et al. (Reference Cui, Wang and Xu2024), the estimated parameters reveal a clear bi-factor loading structure. The first plot in Figure 5 illustrates the first dimension of the loadings parameter (i.e., the discrimination parameters $\boldsymbol {a}$ ), which primarily reflects students’ mathematical problem-solving ability. The second plot in Figure 5 depicts the second dimension of the loadings, which predominantly measures students’ science ability. These two dimensions represent distinct but correlated aspects of overall cognitive ability, capturing the unique contributions of mathematics and science proficiencies. The third plot in Figure 5 demonstrates that the intercept parameters $\boldsymbol {d}$ estimated using both the SNA and the EM algorithm are very similar. Overall, the estimated factor structure aligns with existing literature, supporting the validity of the proposed method. These results suggest that the SNA is a reliable and efficient tool for analyzing large-scale assessment data, making it well-suited for real-time or near-real-time educational assessment applications.

Figure 5 Estimates of item parameters for the M2PL model across TSNA, TASNA, and EM algorithm with corresponding step sizes.