A Note on the Connection Between Trek Rules and Separable Nonlinear Least Squares in Linear Structural Equation Models

Abstract

We show that separable nonlinear least squares (SNLLS) estimation is applicable to all linear structural equation models (SEMs) that can be specified in RAM notation. SNLLS is an estimation technique that has successfully been applied to a wide range of models, for example neural networks and dynamic systems, often leading to improvements in convergence and computation time. It is applicable to models of a special form, where a subset of parameters enters the objective linearly. Recently, Kreiberg et al. (Struct Equ Model Multidiscip J 28(5):725–739, 2021. https://doi.org/10.1080/10705511.2020.1835484) have shown that this is also the case for factor analysis models. We generalize this result to all linear SEMs. To that end, we show that undirected effects (variances and covariances) and mean parameters enter the objective linearly, and therefore, in the least squares estimation of structural equation models, only the directed effects have to be obtained iteratively. For model classes without unknown directed effects, SNLLS can be used to analytically compute least squares estimates. To provide deeper insight into the nature of this result, we employ trek rules that link graphical representations of structural equation models to their covariance parametrization. We further give an efficient expression for the gradient, which is crucial to make a fast implementation possible. Results from our simulation indicate that SNLLS leads to improved convergence rates and a reduced number of iterations.

In the behavioral and social sciences, structural equation models (SEMs) have become widely accepted as a multivariate statistical tool for modeling the relation between latent and observed variables. Apart from maximum likelihood estimation, least squares (LS) estimation is a common approach for parameter estimation. In LS, parameters are estimated by minimizing a nonlinear function of the parameters and data. In practice, this problem is typically solved by applying generic nonlinear optimization techniques, such as Newton-type gradient descent approaches that iteratively minimize the objective function until convergence is reached. However, for some model classes, generic optimization algorithms can be adapted to make better use of the model structure and thus solve the problem more efficiently. For a particular type of models, the parameters separate, that is, one set of parameters enters the objective in a nonlinear way, while another set of parameters enters the objective linearly. For a vector of observations y and predictors x of size m, the objective is of the form

(1) $\begin{array}{c}\sum _{i=1}^m{\left(y_i,-,\sum _{j=1}^n,{\alpha }_j,{\phi }_j,(\beta ,x_i)\right)}^2\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \sum _{i = 1}^m \left[ y_i - \sum _{j=1}^n \alpha _j \varphi _j(\beta , x_i)\right] ^2 \end{aligned}$$\end{document}

where $\alpha \in R^n,\beta \in R^k$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\alpha \in \mathbb {R}^n, \beta \in \mathbb {R}^k$$\end{document} are parameter vectors and the (nonlinear) functions ${\phi }_j$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varphi _j$$\end{document} are continuously differentiable w.r.t. $\beta$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\beta $$\end{document} . Golub and Pereyra (1973) showed that this kind of objective allows for a reformulation of the optimization problem, such that only the parameters $\beta$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\beta $$\end{document} have to be obtained iteratively, while the parameters $\alpha$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\alpha $$\end{document} can be computed after the optimization in a single step. The procedure has been subsequently called separable nonlinear least squares (SNLLS). It has been successfully applied in many disciplines, and it has been observed that the reduced dimension of the parameter space can lead to reduced computation time, a reduced number of iterations and better convergence properties (Golub & Pereyra, 2003). Inspired by earlier work (Kreiberg et al., 2016, 2021) that showed that this procedure can also be applied to factor analysis models, we generalize their result to the entire class of linear structural equation models and give analytical gradients for the reduced optimization problem, which is central for an efficient implementation.

1. Review of Concepts

In the following, we briefly review the notation for structural equation models, the generalized least squares estimator and the trek rules used to derive the model-implied covariance matrix.

1.1. Linear Structural Equation Models

Linear structural equation models can be defined in RAM notation (reticular action model; McArdle & McDonald, 1984) as follows (we follow the notation from Drton, 2018): Let $x,\epsilon$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$x, \varepsilon $$\end{document} be random vectors with values in $R^m$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\mathbb {R}^m$$\end{document} and

(2) $\begin{array}{c}x=\Lambda x+\epsilon \end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} x = \varvec{\Lambda }x + \varepsilon \end{aligned}$$\end{document}

where $\Lambda \in R^{m\times m}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Lambda }\in \mathbb {R}^{m \times m}$$\end{document} is a matrix of constants or unknown (directed) parameters. Let $\Omega \in R^{m\times m}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Omega }\in \mathbb {R}^{m \times m}$$\end{document} be the covariance matrix of $\epsilon$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varepsilon $$\end{document} and $I$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{I}$$\end{document} the identity matrix. If $I-\Lambda$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{I}- \varvec{\Lambda }$$\end{document} is invertible, Eq. 2 can be solved by $x={(I-\Lambda )}^{-1}\epsilon$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$x = (\textbf{I}-\varvec{\Lambda })^{-1}\varepsilon $$\end{document} with covariance matrix

(3) $\begin{array}{c}V\left[x\right]={\left(I,-,\Lambda \right)}^{-1}\Omega {\left(I,-,\Lambda \right)}^{-T}\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \mathbb {V}[x] = \left( \textbf{I}-\varvec{\Lambda }\right) ^{-1} \varvec{\Omega }\left( \textbf{I}-\varvec{\Lambda }\right) ^{-T} \end{aligned}$$\end{document}

If x is partitioned into a part $x_{\text{obs}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$x_{\text {obs}}$$\end{document} of $m_{\text{obs}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$m_\textrm{obs}$$\end{document} observed variables and $x_{\text{lat}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$x_{\text {lat}}$$\end{document} of $m_{\text{lat}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$m_\textrm{lat}$$\end{document} latent variables, we can reorder x such that $x={\left(x_{\text{obs}}^T,,x_{\text{lat}}^T\right)}^T$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$x = \left( x_{\text {obs}}^T \; x_{\text {lat}}^T\right) ^T$$\end{document} , and the covariance matrix of the observed variables is given by

(4) $\begin{array}{c}\Sigma :=V\left[x_{\text{obs}}\right]=F{(I-\Lambda )}^{-1}\Omega {(I-\Lambda )}^{-T}{F}^T\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \varvec{\Sigma }:=\mathbb {V}[x_{\text {obs}}] = \textbf{F}(\textbf{I}-\varvec{\Lambda })^{-1}\varvec{\Omega }(\textbf{I}-\varvec{\Lambda })^{-T}\textbf{F}^T\end{aligned}$$\end{document}

where $F=\left(I,,|,,0\right)\in R^{m_{\text{obs}}\times (m_{\text{obs}}+m_{\text{lat}})}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{F}= \left[ {\textbf {I}}\,\big |\,{\textbf {0}}\right] \in \mathbb {R}^{m_\textrm{obs} \times (m_\textrm{obs} +m_\textrm{lat})}$$\end{document} is a rectangular filter matrix. We denote the parameters by $\theta ={\left({\theta }_{\Lambda }^T{\Omega }_{\Lambda }^T\right)}^T\in R^q$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta = \left( {\theta _{\Lambda }^{T} \Omega _{\Lambda }^{T}}\right) ^{T}\in \mathbb {R}^{q}$$\end{document} , partitioned into directed parameters from $\Lambda$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Lambda }$$\end{document} and undirected parameters from $\Omega$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Omega }$$\end{document} . (We call them directed or undirected parameters because they correspond to directed or undirected paths in the graph of the model.) If we want to stress that $\Sigma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }$$\end{document} is a function of the parameters, we write $\Sigma \left(\theta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }(\theta )$$\end{document} . If we are also interested in the mean structure, we introduce a vector of (possibly zero) mean parameters $\gamma \in R^m$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma \in \mathbb {R}^m$$\end{document} such that $x=\gamma +\Lambda x+\epsilon$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$x = \gamma + \varvec{\Lambda }x + \varepsilon $$\end{document} and obtain

(5) $\begin{array}{c}\mu :=E\left[x_{\mathrm{obs}}\right]=F{(I-\Lambda )}^{-1}\gamma \end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \mu :=\mathbb {E}[x_{obs}] = \textbf{F}(\textbf{I}- \varvec{\Lambda })^{-1}\gamma \end{aligned}$$\end{document}

1.2. Least Squares Estimation

The least squares objective function for $\theta$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta $$\end{document} is:

(6) $\begin{array}{c}{F}_{\text{LS}}={(s-\sigma )}^{T}V(s-\sigma )\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} F_{\text {LS}} = (s - \sigma )^T \textbf{V}(s - \sigma ) \end{aligned}$$\end{document}

where $\sigma =\text{vech}\left(\Sigma \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\sigma = {{\,\textrm{vech}\,}}(\varvec{\Sigma })$$\end{document} is the half-vectorization of $\Sigma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }$$\end{document} , that is, the vector of non-duplicated elements of the model-implied covariance matrix, $s=\text{vech}\left(S\right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$s = {{\,\textrm{vech}\,}}(\textbf{S})$$\end{document} is the half-vectorization of the observed covariance matrix and $V$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{V}$$\end{document} is a fixed symmetric positive definite weight matrix. Specific forms of the weight matrix $V$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{V}$$\end{document} lead to commonly used special cases of this estimation technique: Generalized least squares estimation uses $V=\frac{1}{2}{D}^T\left(S^{-1},\otimes ,S^{-1}\right)D$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{V}= \frac{1}{2} \, \textbf{D}^T\left( \textbf{S}^{-1} \otimes \textbf{S}^{-1}\right) \textbf{D}$$\end{document} (where $D$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{D}$$\end{document} denotes the duplication matrix from Magnus and Neudecker (2019b)), asymptotic distribution-free estimation uses a consistent estimator of the asymptotic covariance matrix of s, and unweighted least squares estimation uses the identity matrix (Bollen, 1989; Browne, 1982, 1984).

1.3. Trek Rules

To show that in SEM undirected effects enter the least squares objective linearly, we employ trek rules (Drton, 2018), which are path tracing rules used to derive the model-implied covariance between any pair of variables in a SEM (Boker et al., 2002). Various authors have proposed rules to link the graph to the covariance parametrization of the model. Here, we give the rules as put forward by Drton (2018), which are based on treks as basic building blocks (for an overview of alternative formulations see Mulaik, 2009). A trek $\tau$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tau $$\end{document} from a node i to j is a path connecting them, where directed edges can be traveled forwards and backwards, but it is not allowed to walk from one arrowhead into another (without colliding arrowheads). A top node of a trek is a node which has only outgoing edges.

To derive an expression for the model-implied covariance between any two variables i and j based on the postulated SEM, we follow 4 steps:

1. Find all treks from i to j.

2. For each trek, multiply all parameters along it.

3. If a trek does not contain a covariance parameter, factor in the variance of the top node.

4. Add all obtained trek monomials from the different treks together.

Note that a trek is ordered in the sense that two treks containing exactly the same nodes and edges are considered different if they are traveled in a different order. In particular, each trek has a source (i) and a target (j), and a trek from j to i is considered to be a different trek, even if it contains exactly the same nodes and edges. Also note that variances are not considered to be edges in the mixed graph corresponding to the model (i.e., it is not allowed to travel variance edges). Therefore, all graphical representations of SEMs in this article omit variance edges, and it is required to factor them in according to rule 3 after the treks are collected.

1.3.1. Example

To illustrate how the model-implied covariances can be derived using trek rules, we give an example based on the graph shown in the path diagram in Fig. 1. To find the model-implied covariance between nodes $X_2$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$X_2$$\end{document} and $X_6$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$X_6$$\end{document} in the model shown in Fig. 1, we first find all treks from $X_2$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$X_2$$\end{document} to $X_6$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$X_6$$\end{document}

(7)

(8)

We now compute the trek monomials for each trek. The second trek does not contain a covariance parameter, so we need to factor in the variance of the top node. We find the trek’s top node G and denote the variance parameter of G by ${\omega }_G$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\omega _G$$\end{document} . Finally, we add the resulting trek monomials and we find that the model-implied covariance between $X_2$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$X_2$$\end{document} and $X_6$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$X_6$$\end{document} can be expressed as follows:

(9) $\begin{array}{c}\text{cov}(X_2,X_6)={\lambda }_2{\omega }_l{\lambda }_6+{\beta }_2{\omega }_G{\beta }_6\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {{\,\textrm{cov}\,}}(X_2, X_6) = \lambda _2\omega _l\lambda _6 + \beta _2\omega _G\beta _6 \end{aligned}$$\end{document}

As a second example, we derive the model-implied variance of $X_3$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$X_3$$\end{document} . Again, we first find all treks from $X_3$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$X_3$$\end{document} to $X_3$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$X_3$$\end{document} :

(10)

(11)

(12) $\begin{array}{c}{X}_3\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} X_3 \end{aligned}$$\end{document}

All treks do not contain a covariance parameter, so we need to factor in the variance of the respective top nodes ${\zeta }_1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\zeta _1$$\end{document} , G and $X_3$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$X_3$$\end{document} . We denote the variance parameters of ${\zeta }_1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\zeta _1$$\end{document} and $X_3$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$X_3$$\end{document} by ${\omega }_{{\zeta }_1}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\omega _{\zeta _1}$$\end{document} and ${\omega }_3$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\omega _3$$\end{document} and add the resulting trek monomials to obtain

(13) $\begin{array}{c}\text{var}\left(X_3\right)={\lambda }_3^2{\omega }_{{\zeta }_1}+{\beta }_3^2{\omega }_G+{\omega }_3\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {{\,\textrm{var}\,}}(X_3) = \lambda _3^2\omega _{\zeta _1} + \beta _3^2\omega _G + \omega _3 \end{aligned}$$\end{document}

Figure 1. Graph of a bi-factor model with one general factor and two specific factors. Circles represent latent variables, and rectangles represent observed variables. Variances are omitted in this representation.

1.3.2. Formal Definitions

We denote the elements of $\Omega$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Omega }$$\end{document} , the undirected effects between nodes i and j, by ${\omega }_{\mathrm{ij}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\omega _{ij}$$\end{document} and the elements of $\Lambda$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Lambda }$$\end{document} , the directed effects, by ${\lambda }_{\mathrm{ij}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\lambda _{ij}$$\end{document} . Drton (2018) defines a trek monomial of a trek $\tau$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tau $$\end{document} without a covariance parameter as

(14) $\begin{array}{c}\tau (\Lambda ,\Omega )={\omega }_{i_0{i}_0}\prod _{k\to l\in \tau }{\lambda }_{\mathrm{lk}}\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \tau (\varvec{\Lambda }, \varvec{\Omega }) = \omega _{i_0 i_0} \prod _{k \rightarrow l \in \tau } \lambda _{lk} \end{aligned}$$\end{document}

where $i_0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$i_0$$\end{document} is the top node of the trek, and a trek monomial of a trek $\tau$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tau $$\end{document} containing an undirected edge between $i_0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$i_0$$\end{document} and $j_0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$j_0$$\end{document} as

(15) $\begin{array}{c}\tau (\Lambda ,\Omega )={\omega }_{i_0{j}_0}\prod _{k\to l\in \tau }{\lambda }_{\mathrm{lk}}\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \tau (\varvec{\Lambda }, \varvec{\Omega }) = \omega _{i_0 j_0} \prod _{k \rightarrow l \in \tau } \lambda _{lk} \end{aligned}$$\end{document}

(notice the swapped indices of ${\lambda }_{\mathrm{lk}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\lambda _{lk}$$\end{document} compared to the formula in Drton because our $\Lambda$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Lambda }$$\end{document} corresponds to his ${\Lambda }^T$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Lambda }^T$$\end{document} ). With this, the elements of $\Sigma \left(\theta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }(\theta )$$\end{document} are represented as a summation over treks. He proves that

(16) $\begin{array}{c}\Sigma {\left(\theta \right)}_{\mathrm{ij}}=\sum _{\tau \in T(i,j)}\tau (\Lambda ,\Omega )\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \varvec{\Sigma }(\theta )_{ij} = \sum _{\tau \in {\mathcal {T}}(i,j)} \tau (\varvec{\Lambda }, \varvec{\Omega }) \end{aligned}$$\end{document}

where $T(i,j)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {T}}(i,j)$$\end{document} is the set of all treks from i to j. It follows that the model-implied covariance is a sum of monomials of parameters. Because covariances between the error terms are not transitive, exactly one undirected parameter (variance or covariance) is present in each monomial. Therefore, if all the directed parameters were fixed, the model-implied covariance would be a linear function of the undirected parameters. This is what makes the SNLLS procedure applicable to structural equation models.

For later use, we also note that Drton gives the following expression:

(17) $\begin{array}{c}{(I-\Lambda )}_{\mathrm{ij}}^{-1}=\sum _{\tau \in P(j,i)}\prod _{k\to l\in \tau }{\lambda }_{\mathrm{lk}}\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} (\textbf{I}- \varvec{\Lambda })^{-1}_{ij} = \sum _{\tau \in {\mathcal {P}}(j,i)} \prod _{k \rightarrow l \in \tau } \lambda _{lk} \end{aligned}$$\end{document}

where $P(j,i)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {P}}(j,i)$$\end{document} is the set of directed paths from j to i. This is because we can write ${(I-\Lambda )}^{-1}={\sum }_{k=0}^{\infty }{\Lambda }^k$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$(\textbf{I}- \varvec{\Lambda })^{-1} = \sum _{k = 0}^\infty \varvec{\Lambda }^k$$\end{document} , where the geometric series converges iff all eigenvalues of $\Lambda$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Lambda }$$\end{document} lie in $(-1,1)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$(-1, 1)$$\end{document} . (Further explanations about this and an excellent account of the connections between matrix algebra and graphs can be found in Kepner & Gilbert (2011).)

2. Separable Nonlinear Least Squares for SEM

We first outline the proofs for the applicability of SNLLS to CFA as given by Golub and Pereyra (1973) and Kreiberg et al. (2021). Subsequently, we proof that SNLLS is applicable to linear structural equation models. We further extend the existing proofs to subsume models that contain a mean structure. Last, we derive analytic gradients that are central for efficient software implementations.

2.1. Outline of Previous Work

To minimize Eq. 1, Golub and Pereyra (1973) define the matrix function

(18) $\begin{array}{c}{\Phi }_{\mathrm{ij}}:={\phi }_j(\beta ,x_i)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \Phi _{ij} :=\varphi _j(\beta , x_i) \end{aligned}$$\end{document}

such that Eq. 1 can be written as

(19) $\begin{array}{c}{\Vert y-\Phi \left(\beta \right)\alpha \Vert }^2\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \Vert y - \Phi (\beta ) \alpha \Vert ^2 \end{aligned}$$\end{document}

where $\Vert ·\Vert$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\Vert \cdot \Vert $$\end{document} denotes the euclidean norm. For a fixed value of $\beta$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\beta $$\end{document} , a solution for $\alpha$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\alpha $$\end{document} can be obtained as $\alpha ={\Phi }^{+}\left(\beta \right)y$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\alpha = \Phi ^+(\beta ) y$$\end{document} . They further proved that under the assumption that $\Phi \left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\Phi (\beta )$$\end{document} has constant rank near the solution, only the nonlinear parameters $\beta$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\beta $$\end{document} have to be obtained iteratively by replacing $\alpha$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\alpha $$\end{document} and minimizing the modified objective

(20) $\begin{array}{c}{\left(y,-,\Phi ,\left(\beta \right),{\Phi }^{+},\left(\beta \right),y\right)}^2\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \left\Vert y - \Phi (\beta )\Phi ^+(\beta ) y \right\Vert ^2 \end{aligned}$$\end{document}

where ${\Phi }^{+}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\Phi ^+$$\end{document} denotes the Moore–Penrose generalized inverse. Afterward, the least squares solution for the linear parameters $\alpha$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\alpha $$\end{document} can be obtained as the standard least squares estimator ${\arg \min }_{\alpha \in R^n}\Vert \Phi \left(\hat{\beta }\right)\alpha -y\Vert ={\Phi }^{+}\left(\hat{\beta }\right)y$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${{\,\mathrm{arg\,min}\,}}_{\alpha \in \mathbb {R}^n} \Vert \Phi (\hat{\beta })\alpha - y \Vert = \Phi ^+(\hat{\beta }) y$$\end{document} .

Kreiberg et al. (2021) showed that this procedure is applicable for CFA models (we reproduce their main results in our notation), as it is possible to rewrite the model-implied covariances $\sigma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\sigma $$\end{document} as a product of a matrix-valued function $G\left({\theta }_{\Lambda }\right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{G}(\theta _{\varvec{\Lambda }})$$\end{document} (that depends only on the directed parameters) and the undirected parameters ${\theta }_{\Omega }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{\varvec{\Omega }}$$\end{document} , so the LS objective can be written as

(21) $\begin{array}{cc}{F}_{\text{LS}}& ={(s-\sigma )}^{T}V(s-\sigma )\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} F_{\text {LS}}&= (s - \sigma )^T \textbf{V}(s - \sigma ) \end{aligned}$$\end{document}

(22) $\begin{array}{cc}& ={\left(s,-,G,\left({\theta }_{\Lambda }\right),{\theta }_{\Omega }\right)}^{T}V\left(s,-,G,\left({\theta }_{\Lambda }\right),{\theta }_{\Omega }\right)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned}&= \left( s - \textbf{G}(\theta _{\varvec{\Lambda }})\theta _{\varvec{\Omega }}\right) ^T \textbf{V}\left( s - \textbf{G}(\theta _{\varvec{\Lambda }})\theta _{\varvec{\Omega }} \right) \end{aligned}$$\end{document}

(23) $\begin{array}{cc}& ={\left(s,-,G,\left({\theta }_{\Lambda }\right),{\theta }_{\Omega }\right)}_V^2\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned}&= \left\Vert s - \textbf{G}(\theta _{\varvec{\Lambda }})\theta _{\varvec{\Omega }}\right\Vert ^2_\textbf{V}\end{aligned}$$\end{document}

They further stated that if ${\theta }_{\Lambda }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{\varvec{\Lambda }}$$\end{document} is fixed, we know from standard linear least squares estimation that the minimizer for the undirected effects can be obtained as

(24) $\begin{array}{c}{\hat{\theta }}_{\Omega }={\left(G,{\left({\theta }_{\Lambda }\right)}^T,V,G,\left({\theta }_{\Lambda }\right)\right)}^{-1}G{\left({\theta }_{\Lambda }\right)}^{T}Vs\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \hat{\theta }_{\varvec{\Omega }} = \left( \textbf{G}( \theta _{\varvec{\Lambda }})^T \textbf{V}\textbf{G}( \theta _{\varvec{\Lambda }})\right) ^{-1} \textbf{G}( \theta _{\varvec{\Lambda }})^T\textbf{V}s \end{aligned}$$\end{document}

Inserting Eq. 24 into Eq. 22 and simplifying, they obtained a new objective to be minimized:

(25) $\begin{array}{cc}{\hat{\theta }}_{\Lambda }& =\underset{{\theta }_{\Lambda }}{\arg \min }\left(s^T,V,s,-,s^T,V,G,\left({\theta }_{\Lambda }\right),{\left(G,{\left({\theta }_{\Lambda }\right)}^T,V,G,\left({\theta }_{\Lambda }\right)\right)}^{-1},G,{\left({\theta }_{\Lambda }\right)}^T,V,s\right)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \hat{\theta }_{\varvec{\Lambda }}&= \mathop {{{\,\mathrm{arg\,min}\,}}}\limits _{\theta _{\varvec{\Lambda }}} \left[ s^T\textbf{V}s - s^T\textbf{V}\textbf{G}(\theta _{\varvec{\Lambda }}) \left( \textbf{G}(\theta _{\varvec{\Lambda }})^T \textbf{V}\textbf{G}(\theta _{\varvec{\Lambda }})\right) ^{-1} \textbf{G}(\theta _{\varvec{\Lambda }})^T\textbf{V}s \right] \end{aligned}$$\end{document}

(26) $\begin{array}{cc}& =\underset{{\theta }_{\Lambda }}{\arg \min }{F}_{\text{SNLLS}}\left({\theta }_{\Lambda }\right)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned}&= \mathop {{{\,\mathrm{arg\,min}\,}}}\limits _{\theta _{\varvec{\Lambda }}} F_{\text {SNLLS}}(\theta _{\varvec{\Lambda }}) \end{aligned}$$\end{document}

This objective only depends on the directed parameters ${\theta }_{\Lambda }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{\varvec{\Lambda }}$$\end{document} . After minimizing it to obtain a LS estimate ${\hat{\theta }}_{\Lambda }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\hat{\theta }_{\varvec{\Lambda }}$$\end{document} , Eq. 24 can be used to obtain the LS estimate of ${\theta }_{\Omega }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{\varvec{\Omega }}$$\end{document} . We would like to note they assume that $G$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{G}$$\end{document} has full rank, which is not a necessary assumption and can be relaxed using alternative formulations of Eqs. 24 and 25. To extend the method to general structural equation models, we have to derive $G\left({\theta }_{\Lambda }\right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{G}(\theta _{\varvec{\Lambda }})$$\end{document} . We do that in the following for all models formulated in the RAM notation.

2.2. Derivation of $G\left({\theta }_{\Lambda }\right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{G}(\theta _{\varvec{\Lambda }})$$\end{document}

Since $F=\left(I,,|,,0\right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{F}= \left[ {\textbf {I}}\,\big |\,{\textbf {0}}\right] $$\end{document} with $0\in R^{m_{\text{obs}}\times m_{\text{lat}}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{0} \in \mathbb {R}^{m_\textrm{obs} \times m_\textrm{lat}}$$\end{document} , the product $FM{F}^T$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{F}\textbf{M}\textbf{F}^{T}$$\end{document} for any $M\in R^{m\times m}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{M}\in \mathbb {R}^{m \times m}$$\end{document} is equal to just deleting the last $m_{\text{lat}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$m_\textrm{lat}$$\end{document} rows and columns of $M$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{M}$$\end{document} . We also note that for any matrices $M,D\in R^{n\times n}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{M}, \textbf{D}\in \mathbb {R}^{n \times n}$$\end{document} we can write

(27) $\begin{array}{c}{\left(M,D,M^T\right)}_{\mathrm{ij}}=\sum _{k=1}^n\sum _{l=1}^n{m}_{\mathrm{il}}{d}_{\mathrm{lk}}{m}_{\mathrm{jk}}\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \left( \textbf{M}\textbf{D}\textbf{M}^{T}\right) _{ij} = \sum _{k = 1}^n \sum _{l = 1}^n m_{il} \; d_{lk} \; m_{jk} \end{aligned}$$\end{document}

With this in mind, we can rewrite the model-implied covariance matrix $\Sigma \left(\theta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }(\theta )$$\end{document} as

(28) $\begin{array}{cccccc}\Sigma {\left(\theta \right)}_{\mathrm{ij}}& =(F& {(I-\Lambda )}^{-1}& \Omega & {(I-\Lambda )}^{-T}& F^T{)}_{\mathrm{ij}}\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \varvec{\Sigma }(\theta )_{ij}&= \bigg ( \textbf{F}{} & {} (\textbf{I}-\varvec{\Lambda })^{-1}{} & {} \varvec{\Omega }&(\textbf{I}-\varvec{\Lambda })^{-T}{} & {} \textbf{F}^T \bigg )_{ij} \end{aligned}$$\end{document}

(29) $\begin{array}{cccccc}& =(& {(I-\Lambda )}^{-1}& \Omega & {(I-\Lambda )}^{-T}& {)}_{\mathrm{ij}}\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned}&= \bigg ({} & {} (\textbf{I}-\varvec{\Lambda })^{-1}{} & {} \varvec{\Omega }&(\textbf{I}-\varvec{\Lambda })^{-T}{} & {} \bigg )_{ij} \end{aligned}$$\end{document}

(30) $\begin{array}{cccccc}& =\sum _{k=1}^m\sum _{l=1}^m& {(I-\Lambda )}_{\mathrm{il}}^{-1}& {\omega }_{\mathrm{lk}}& {(I-\Lambda )}_{\mathrm{jk}}^{-1}& \end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned}&= \; \; \sum _{k = 1}^m \sum _{l = 1}^m{} & {} (\textbf{I}- \varvec{\Lambda })^{-1}_{il}{} & {} \omega _{lk}&(\textbf{I}- \varvec{\Lambda })^{-1}_{jk}{} & {} \end{aligned}$$\end{document}

(31) $\begin{array}{cccccc}& =\sum _{k=1}^m\sum _{l=1}^m& (\sum _{\tau \in P(l,i)}\prod _{r\to s\in \tau }{\lambda }_{\mathrm{sr}})& {\omega }_{\mathrm{lk}}& (\sum _{\tau \in P(k,j)}\prod _{r\to s\in \tau }{\lambda }_{\mathrm{sr}})& \end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned}&= \; \; \sum _{k = 1}^m \sum _{l = 1}^m{} & {} \Big (\sum _{\tau \in \mathcal {P}(l,i)} \prod _{r \rightarrow s \in \tau } \lambda _{sr}\Big ){} & {} \omega _{lk}&\Big (\sum _{\tau \in \mathcal {P}(k,j)} \prod _{r \rightarrow s \in \tau } \lambda _{sr}\Big ){} & {} \end{aligned}$$\end{document}

with $i,j\in \{1,\dots ,m_{\text{obs}}\}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$i,j \in \{1, \ldots , m_\textrm{obs}\}$$\end{document} . We now immediately see that each entry of $\Sigma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }$$\end{document} is a sum of products of entries of ${(I-\Lambda )}^{-1}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$(\textbf{I}-\varvec{\Lambda })^{-1}$$\end{document} and $\Omega$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Omega }$$\end{document} . More importantly, exactly one entry of $\Omega$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Omega }$$\end{document} enters each term of the sum; if we keep all entries of $\Lambda$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Lambda }$$\end{document} fixed, each element in $\Sigma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }$$\end{document} is a linear function of the entries of $\Omega$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Omega }$$\end{document} and is therefore a linear function of the undirected parameters in $\Omega$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Omega }$$\end{document} (under the assumption that $\Omega$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Omega }$$\end{document} is linearly parameterized). As a result, the parameter vector $\theta$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta $$\end{document} is separable in two parts, ${\theta }_{\Lambda }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{\varvec{\Lambda }}$$\end{document} from $\Lambda$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Lambda }$$\end{document} and ${\theta }_{\Omega }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{\varvec{\Omega }}$$\end{document} from $\Omega$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Omega }$$\end{document} , and ${\theta }_{\Omega }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{\varvec{\Omega }}$$\end{document} enters the computation of the model-implied covariance linearly. As stated before, this is the reason why we will be able to apply separable nonlinear least squares estimation to our problem. Before we proceed, we would like to introduce some notation. If $F$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {F}}$$\end{document} and $G$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {G}}$$\end{document} are tuples of length n and m, and f and g are functions, we define a column vector of length n as

(32) $\begin{array}{c}([f\left(i\right){]}_{i\in F})=\left(\begin{array}{c}f\left(F_1\right)\\ f\left(F_2\right)\\ \dots \\ f\left(F_n\right)\end{array}\right)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \Bigg (\Big [f(i)\Big ]_{i \in {\mathcal {F}}}\Bigg ) = \begin{pmatrix} f({\mathcal {F}}_1) \\ f({\mathcal {F}}_2) \\ \ldots \\ f({\mathcal {F}}_n) \end{pmatrix} \end{aligned}$$\end{document}

and a matrix of size $n\times m$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n \times m$$\end{document} as

(33) $\begin{array}{c}([g(i,j){]}_{i\in F,j\in G})=\left(\begin{array}{cccc}g(F_1,G_1)& g(F_1,G_2)& \dots & g(F_1,G_m)\\ g(F_2,G_1)& g(F_2,G_2)& \dots & g(F_2,G_m)\\ \dots & \dots & \dots & \dots \\ g(F_n,G_1)& \dots & \dots & g(F_n,G_m)\end{array}\right)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \Bigg (\Big [g(i,j)\Big ]_{i \in {\mathcal {F}},\; j \in {\mathcal {G}}}\Bigg ) = \begin{pmatrix} g({\mathcal {F}}_1, {\mathcal {G}}_1) &{}\quad g({\mathcal {F}}_1, {\mathcal {G}}_2) &{}\quad \ldots &{}\quad g({\mathcal {F}}_1, {\mathcal {G}}_m)\\ g({\mathcal {F}}_2, {\mathcal {G}}_1) &{}\quad g({\mathcal {F}}_2, {\mathcal {G}}_2) &{}\quad \ldots &{}\quad g({\mathcal {F}}_2, {\mathcal {G}}_m)\\ \ldots &{}\quad \ldots &{}\quad \ldots &{}\quad \ldots \\ g({\mathcal {F}}_n, {\mathcal {G}}_1) &{}\quad \ldots &{}\quad \ldots &{}\quad g({\mathcal {F}}_n, {\mathcal {G}}_m) \end{pmatrix} \end{aligned}$$\end{document}

To make the subsequent steps easier to follow, we assume that there are no equality constraints between parameters in $\Omega$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Omega }$$\end{document} and no constant terms in $\Omega$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Omega }$$\end{document} different from 0. In Appendices A and B, we show how to lift those assumptions. We now further simplify Eq. 30: Since only nonzero entries of $\Omega$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Omega }$$\end{document} (the parameters ${\theta }_{\Omega }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{\varvec{\Omega }}$$\end{document} ) contribute to the sum, we define $C$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {C}}$$\end{document} as the lower triangular indices of ${\theta }_{\Omega }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{\varvec{\Omega }}$$\end{document} in $\Omega$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Omega }$$\end{document} , i.e., $C_i=(l,k)\in N\times N$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {C}}_i = (l, k) \in \mathbb {N}\times \mathbb {N}$$\end{document} with ${\left({\theta }_{\Omega }\right)}_i={\omega }_{\mathrm{lk}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$({\theta _{\varvec{\Omega }}})_i = \omega _{lk}$$\end{document} and $l\ge k$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$l \ge k$$\end{document} . We now rewrite Eq. 30 by omitting all zero terms:

(34) $\begin{array}{cc}\Sigma {\left(\theta \right)}_{\mathrm{ij}}& =\sum _{(l,k)\in C}\left({(I-\Lambda )}_{\mathrm{il}}^{-1},{\omega }_{\mathrm{lk}},{(I-\Lambda )}_{\mathrm{jk}}^{-1},+,{\delta }_{k\ne l},,{(I-\Lambda )}_{\mathrm{ik}}^{-1},{\omega }_{\mathrm{lk}},{(I-\Lambda )}_{\mathrm{jl}}^{-1}\right)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \varvec{\Sigma }(\theta )_{ij}&= \sum _{(l,k) \in {\mathcal {C}}} \left[ (\textbf{I}- \varvec{\Lambda })^{-1}_{il} \omega _{lk} (\textbf{I}- \varvec{\Lambda })^{-1}_{jk} + \delta _{k \ne l} \, (\textbf{I}- \varvec{\Lambda })^{-1}_{ik} \omega _{lk} (\textbf{I}- \varvec{\Lambda })^{-1}_{jl} \right] \end{aligned}$$\end{document}

(35) $\begin{array}{cc}& ={\left({\left({(I-\Lambda )}_{\mathrm{il}}^{-1},{(I-\Lambda )}_{\mathrm{jk}}^{-1},+,{\delta }_{k\ne l},,{(I-\Lambda )}_{\mathrm{ik}}^{-1},{(I-\Lambda )}_{\mathrm{jl}}^{-1}\right)}_{(l,k)\in C}\right)}^T{\theta }_{\Omega }\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned}&= \left( \left[ (\textbf{I}- \varvec{\Lambda })^{-1}_{il} (\textbf{I}- \varvec{\Lambda })^{-1}_{jk} + \delta _{k \ne l} \, (\textbf{I}- \varvec{\Lambda })^{-1}_{ik} (\textbf{I}- \varvec{\Lambda })^{-1}_{jl}\right] _{(l,k) \in {\mathcal {C}}} \right) ^T \; \theta _{\varvec{\Omega }} \end{aligned}$$\end{document}

where ${\delta }_{k\ne l}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\delta _{k \ne l}$$\end{document} is an indicator function that takes the value 1 if $k\ne l$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$k \ne l$$\end{document} and 0 otherwise. Since we are only interested in the non-duplicated elements $\sigma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\sigma $$\end{document} of $\Sigma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }$$\end{document} , we define another index tuple $D$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {D}}$$\end{document} that denotes the indices of the original position of ${\sigma }_k$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\sigma _k$$\end{document} in $\Sigma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }$$\end{document} , i.e., $D_k=(i,j)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {D}}_k = (i, j)$$\end{document} such that ${\sigma }_k={\Sigma }_{\mathrm{ij}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\sigma _k = \varvec{\Sigma }_{ij}$$\end{document} . This allows us to stack the expression we just found for ${\Sigma }_{\mathrm{ij}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_{ij}$$\end{document} rowwise to get

(36) $\begin{array}{cc}\sigma & =\left([,{\Sigma }_{\mathrm{ij}},{]}_{(i,j)\in D}\right)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \sigma&= \left( \Big [\varvec{\Sigma }_{ij} \Big ]_{(i, j) \in {\mathcal {D}}} \right) \end{aligned}$$\end{document}

(37) $\begin{array}{cc}& =\left({\left({(I-\Lambda )}_{\mathrm{il}}^{-1},{(I-\Lambda )}_{\mathrm{jk}}^{-1},+,{\delta }_{k\ne l},,{(I-\Lambda )}_{\mathrm{ik}}^{-1},{(I-\Lambda )}_{\mathrm{jl}}^{-1}\right)}_{(i,j)\in D,(l,k)\in C},\right){\theta }_{\Omega }\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned}&= \left( \left[ (\textbf{I}- \varvec{\Lambda })^{-1}_{il} (\textbf{I}- \varvec{\Lambda })^{-1}_{jk} + \delta _{k \ne l} \, (\textbf{I}- \varvec{\Lambda })^{-1}_{ik} (\textbf{I}- \varvec{\Lambda })^{-1}_{jl}\right] _{(i, j) \in {\mathcal {D}}, \; (l,k) \in {\mathcal {C}}} \; \right) \theta _{\varvec{\Omega }} \end{aligned}$$\end{document}

(38) $\begin{array}{cc}& =G\left({\theta }_{\Lambda }\right){\theta }_{\Omega }\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned}&= \textbf{G}(\theta _{\varvec{\Lambda }}) \; \theta _{\varvec{\Omega }} \end{aligned}$$\end{document}

where $G\left({\theta }_{\Lambda }\right)\in R^{dim\left(\sigma \right)\times dim\left({\theta }_{\Omega }\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{G}(\theta _{\varvec{\Lambda }}) \in \mathbb {R}^{\dim (\sigma ) \times \dim (\theta _{\varvec{\Omega }})}$$\end{document} . (We let $dim(·)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\dim (\cdot )$$\end{document} of a vector denote its number of elements, i.e., the dimension of the underlying (finite-dimensional) vector space.)

Even though this expression may appear involved, it is in fact easy to compute. Before the optimization procedure starts, we store $C$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {C}}$$\end{document} by looking up the positions of the parameters in $\Omega$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Omega }$$\end{document} and also store $D$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {D}}$$\end{document} . At each step of the optimization procedure, to compute $G\left({\theta }_{\Lambda }\right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{G}(\theta _{\varvec{\Lambda }})$$\end{document} , we now compute ${(I-\Lambda )}^{-1}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$(\textbf{I}- \varvec{\Lambda })^{-1}$$\end{document} first and then loop through the entries $C$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {C}}$$\end{document} and $D$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {D}}$$\end{document} to compute each entry of $G\left({\theta }_{\Lambda }\right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{G}(\theta _{\varvec{\Lambda }})$$\end{document} according to Eq. 37. We note that $G$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{G}$$\end{document} will typically be sparse; therefore, it is advisable to analyze its sparsity pattern previous to the optimization, and only loop through nonzero values.

In Appendix D, we present a different way of obtaining $G\left({\theta }_{\Lambda }\right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{G}(\theta _{\varvec{\Lambda }})$$\end{document} and the gradients, which mimics the approach of Kreiberg et al. (2021). However, the expressions obtained here are computationally more efficient, as the ones in the appendix contain very large Kronecker products.

2.3. Mean Structures

If the model contains mean parameters, we partition the parameter vector $\theta$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta $$\end{document} into three parts: ${\theta }_{\Lambda }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{\varvec{\Lambda }}$$\end{document} and ${\theta }_{\Omega }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{\varvec{\Omega }}$$\end{document} as before, and ${\theta }_{\gamma }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _\gamma $$\end{document} from the mean vector $\gamma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma $$\end{document} . From Eq. 5, we directly see that the model-implied mean vector $\mu \left(\theta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\mu (\theta )$$\end{document} is a linear function of ${\theta }_{\gamma }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _\gamma $$\end{document} . If we let $A$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {A}}$$\end{document} denote the indices of the parameters ${\theta }_{\gamma }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{\gamma }$$\end{document} in $\gamma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma $$\end{document} , i.e., for $i=A_j$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$i = {\mathcal {A}}_{j}$$\end{document} we have ${\left({\theta }_{\gamma }\right)}_j$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$({\theta _{\gamma }})_j$$\end{document} = ${\gamma }_i$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma _i$$\end{document} , we obtain the formula

(39) $\begin{array}{c}\mu =\left({\left({(I-\Lambda )}_{\mathrm{ij}}^{-1}\right)}_{i\in (1,\dots ,m_{\text{obs}}),j\in A}\right){\theta }_{\gamma }\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \mu = \left( \left[ (\textbf{I}- \varvec{\Lambda })^{-1}_{ij}\right] _{i \in (1, \ldots , m_\textrm{obs}), \; j \in {\mathcal {A}}} \right) \; \theta _{\gamma }\end{aligned}$$\end{document}

We now make a slight change in notation: For the previously obtained $G\left(\theta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{G}(\theta )$$\end{document} -matrix, we write $G_{\sigma }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{G}_\sigma $$\end{document} instead and define $G_{\mu }:=\left({\left({(I-\Lambda )}_{\mathrm{ij}}^{-1}\right)}_{i\in (1,\dots ,m_{\text{obs}}),j\in A}\right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{G}_\mu :=\left( \left[ (\textbf{I}- \varvec{\Lambda })^{-1}_{ij}\right] _{i \in (1, \ldots , m_\textrm{obs}), \; j \in {\mathcal {A}}} \right) $$\end{document} . Using a formulation of the least squares objective that also includes a mean structure, we see that

(40)

with

(41) $\begin{array}{c}G:=\left(\begin{array}{cc}{G}_{\sigma }& 0\\ 0& G_{\mu }\end{array}\right)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \textbf{G}:=\left[ \begin{array}{cc} \textbf{G}_\sigma &{}\quad \textbf{0} \\ \textbf{0} &{}\quad \textbf{G}_\mu \end{array}\right] \end{aligned}$$\end{document}

It follows that in addition to the undirected parameters, the mean parameters also do not have to be optimized iteratively but can instead be computed analytically after the iterative optimization is completed.

2.4. Gradient of the SNLLS Objective

There are computationally efficient expression to compute the SNLLS objective and its gradient analytically (Kaufman, 1975; O’Leary & Rust, 2013). Because numerical approximations of the gradient are often slow and may become numerically unstable, we derive an analytical expression for the part of the gradient that is specific to SEMs. We use the notation and methods from Magnus and Neudecker (2019a) and denote the differential by $d$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${{\,\mathrm{\textsf{d}}\,}}$$\end{document} and the Jacobian by $D$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${{\,\mathrm{\textsf{D}}\,}}$$\end{document} . The Jacobian of a matrix function $M$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{M}$$\end{document} with respect to a vector x is defined as $DM=\frac{\partial \text{vec}M}{\partial x^T}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${{\,\mathrm{\textsf{D}}\,}}\textbf{M}= \frac{\partial {{\,\textrm{vec}\,}}\textbf{M}}{\partial x^{T}}$$\end{document} . In the approaches by Kaufman (1975) and O’Leary and Rust (2013), the gradient of the SNLSS objective is expressed in terms of the partial derivatives of the entries of $G$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{G}$$\end{document} w.r.t the nonlinear parameters, i.e., $DG$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${{\,\mathrm{\textsf{D}}\,}}\textbf{G}$$\end{document} . In order to be able to implement such efficient approaches in practice, we derive $DG$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${{\,\mathrm{\textsf{D}}\,}}\textbf{G}$$\end{document} here. We also give the full gradient of Eq. 25 for completeness in Appendix C, although in practice, a more efficient expression from the cited literature can be used (which also does not assume $G$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{G}$$\end{document} to have full rank). For reasons of clarity, we here only consider the case without mean structure, e.g., $G=G_{\sigma }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{G}= \textbf{G}_\sigma $$\end{document} . This is because the derivative of $G_{\mu }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{G}_\mu $$\end{document} is similar to obtain and we do not want to make the derivation unnecessarily technical.

Let $E$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {E}}$$\end{document} denote the indices of ${\theta }_{\Lambda }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{\varvec{\Lambda }}$$\end{document} in $\Lambda$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Lambda }$$\end{document} , i.e., $E_k=(i,j)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {E}}_k = (i,j)$$\end{document} such that ${\Lambda }_{\mathrm{ij}}={\left({\theta }_{\Lambda }\right)}_k$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Lambda }_{ij} = ({\theta _{\varvec{\Lambda }}})_k$$\end{document} . We note that

(42) $\begin{array}{c}\frac{\partial {(I-\Lambda )}_{\mathrm{kl}}^{-1}}{\partial {\Lambda }_{\mathrm{ij}}}={(I-\Lambda )}_{\mathrm{ki}}^{-1}{(I-\Lambda )}_{\mathrm{jl}}^{-1}\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \frac{\partial (\textbf{I}- \varvec{\Lambda })^{-1}_{kl}}{\partial \varvec{\Lambda }_{ij}} = (\textbf{I}- \varvec{\Lambda })^{-1}_{ki}(\textbf{I}- \varvec{\Lambda })^{-1}_{jl} \end{aligned}$$\end{document}

With this, we derive the partial derivatives of each entry of $G$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{G}$$\end{document} in terms of the matrix ${(I-\Lambda )}^{-1}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$(\textbf{I}- \varvec{\Lambda })^{-1}$$\end{document} as

(43) $\begin{array}{cc}\frac{\partial G_{r,s}}{\partial {\left({\theta }_{\Lambda }\right)}_n}& =\frac{\partial }{\partial {\left({\theta }_{\Lambda }\right)}_n}\left({(I-\Lambda )}_{\mathrm{il}}^{-1},{(I-\Lambda )}_{\mathrm{jk}}^{-1},+,{\delta }_{k\ne l},,{(I-\Lambda )}_{\mathrm{ik}}^{-1},{(I-\Lambda )}_{\mathrm{jl}}^{-1}\right)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \frac{\partial \textbf{G}_{r, s}}{\partial ({\theta _{\varvec{\Lambda }}})_n}&= \frac{\partial }{\partial ({\theta _{\varvec{\Lambda }}})_n} \left[ (\textbf{I}- \varvec{\Lambda })^{-1}_{il} (\textbf{I}- \varvec{\Lambda })^{-1}_{jk} +\delta _{k \ne l} \, (\textbf{I}- \varvec{\Lambda })^{-1}_{ik} (\textbf{I}- \varvec{\Lambda })^{-1}_{jl}\right] \end{aligned}$$\end{document}

(44) $\begin{array}{cc}& =\left(\frac{\partial }{\partial {\left({\theta }_{\Lambda }\right)}_n},{(I-\Lambda )}_{\mathrm{il}}^{-1},{(I-\Lambda )}_{\mathrm{jk}}^{-1},+,{(I-\Lambda )}_{\mathrm{il}}^{-1},\frac{\partial }{\partial {\left({\theta }_{\Lambda }\right)}_n},{(I-\Lambda )}_{\mathrm{jk}}^{-1}\right)\\ & +{\delta }_{k\ne l}\left(\frac{\partial }{\partial {\left({\theta }_{\Lambda }\right)}_n},{(I-\Lambda )}_{\mathrm{ik}}^{-1},{(I-\Lambda )}_{\mathrm{jl}}^{-1},+,{(I-\Lambda )}_{\mathrm{ik}}^{-1},\frac{\partial }{\partial {\left({\theta }_{\Lambda }\right)}_n},{(I-\Lambda )}_{\mathrm{jl}}^{-1}\right)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned}&= \left[ \frac{\partial }{\partial ({\theta _{\varvec{\Lambda }}})_n}(\textbf{I}- \varvec{\Lambda })^{-1}_{il} (\textbf{I}- \varvec{\Lambda })^{-1}_{jk} + (\textbf{I}- \varvec{\Lambda })^{-1}_{il} \frac{\partial }{\partial ({\theta _{\varvec{\Lambda }}})_n}(\textbf{I}- \varvec{\Lambda })^{-1}_{jk}\right] \nonumber \\&\quad +\delta _{k \ne l} \left[ \frac{\partial }{\partial ({\theta _{\varvec{\Lambda }}})_n}(\textbf{I}- \varvec{\Lambda })^{-1}_{ik} (\textbf{I}- \varvec{\Lambda })^{-1}_{jl} + (\textbf{I}- \varvec{\Lambda })^{-1}_{ik} \frac{\partial }{\partial ({\theta _{\varvec{\Lambda }}})_n}(\textbf{I}- \varvec{\Lambda })^{-1}_{jl} \right] \end{aligned}$$\end{document}

(45) $\begin{array}{cc}& =\left({(I-\Lambda )}_{\mathrm{iu}}^{-1},{(I-\Lambda )}_{\mathrm{vl}}^{-1},{(I-\Lambda )}_{\mathrm{jk}}^{-1},+,{(I-\Lambda )}_{\mathrm{il}}^{-1},{(I-\Lambda )}_{\mathrm{ju}}^{-1},{(I-\Lambda )}_{\mathrm{vk}}^{-1}\right)\\ & +{\delta }_{k\ne l}\left({(I-\Lambda )}_{\mathrm{iu}}^{-1},{(I-\Lambda )}_{\mathrm{vk}}^{-1},{(I-\Lambda )}_{\mathrm{jl}}^{-1},+,{(I-\Lambda )}_{\mathrm{ik}}^{-1},{(I-\Lambda )}_{\mathrm{ju}}^{-1},{(I-\Lambda )}_{\mathrm{vl}}^{-1}\right)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned}&= \left[ (\textbf{I}- \varvec{\Lambda })^{-1}_{iu} (\textbf{I}- \varvec{\Lambda })^{-1}_{vl} (\textbf{I}- \varvec{\Lambda })^{-1}_{jk} + (\textbf{I}- \varvec{\Lambda })^{-1}_{il} (\textbf{I}- \varvec{\Lambda })^{-1}_{ju} (\textbf{I}- \varvec{\Lambda })^{-1}_{vk} \right] \nonumber \\&\quad +\delta _{k \ne l} \left[ (\textbf{I}- \varvec{\Lambda })^{-1}_{iu} (\textbf{I}- \varvec{\Lambda })^{-1}_{vk} (\textbf{I}- \varvec{\Lambda })^{-1}_{jl} + (\textbf{I}- \varvec{\Lambda })^{-1}_{ik} (\textbf{I}- \varvec{\Lambda })^{-1}_{ju} (\textbf{I}- \varvec{\Lambda })^{-1}_{vl} \right] \end{aligned}$$\end{document}

with $(i,j)=D_r$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$(i, j) = {\mathcal {D}}_r$$\end{document} , $(l,k)=C_s$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$(l, k) = {\mathcal {C}}_s$$\end{document} , and $(u,v)=E_n$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$(u, v) = {\mathcal {E}}_n$$\end{document} . Since $G$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{G}$$\end{document} is of dimension $dim\left(\sigma \right)\times dim\left({\theta }_{\Omega }\right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\dim (\sigma ) \times \dim (\theta _{\varvec{\Omega }})$$\end{document} , with $k=dim\left(\sigma \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$k = \dim (\sigma )$$\end{document} we have

(46) $\begin{array}{c}\text{vec}{\left(G\right)}_t=G_{t-k\lfloor (t-1)/k)\rfloor ,\lceil t/k\rceil }\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {{\,\textrm{vec}\,}}(\textbf{G})_t = \textbf{G}_{t - k\lfloor (t-1)/k) \rfloor , \; \lceil t / k \rceil } \end{aligned}$$\end{document}

and we obtain $DG\in R^{dim\left(\sigma \right)dim\left({\theta }_{\Omega }\right)\times dim\left({\theta }_{\Lambda }\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${{\,\mathrm{\textsf{D}}\,}}\textbf{G}\in \mathbb {R}^{\dim (\sigma )\dim (\theta _{\varvec{\Omega }}) \times \dim (\theta _{\varvec{\Lambda }})}$$\end{document} as

(47) $\begin{array}{c}DG=\frac{\partial \text{vec}G}{\partial {\theta }_{\Lambda }^T}=\left({\left(\frac{\partial G_{t-k\lfloor (t-1)/k)\rfloor ,\lceil t/k\rceil }}{\partial {\left({\theta }_{\Lambda }\right)}_n}\right)}_{t\in (1,\dots ,dim\left(\sigma \right)dim\left({\theta }_{\Omega }\right)),n\in (1,\dots ,dim\left({\theta }_{\Lambda }\right))}\right)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {{\,\mathrm{\textsf{D}}\,}}\textbf{G}= \frac{\partial {{\,\textrm{vec}\,}}\textbf{G}}{\partial \theta _{\varvec{\Lambda }}^T} = \left( \left[ \frac{\partial \textbf{G}_{t - k\lfloor (t-1)/k) \rfloor , \; \lceil t / k \rceil }}{\partial ({\theta _{\varvec{\Lambda }}})_n}\right] _{t \in (1, \ldots , \dim (\sigma )\dim (\theta _{\varvec{\Omega }})), \; n \in (1, \ldots , \dim (\theta _{\varvec{\Lambda }}))}\right) \end{aligned}$$\end{document}

To facilitate software implementation, we give a way to compute $DG$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${{\,\mathrm{\textsf{D}}\,}}\textbf{G}$$\end{document} in pseudocode in Algorithm 1. In practice, $DG$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${{\,\mathrm{\textsf{D}}\,}}\textbf{G}$$\end{document} will typically contain many zero values. Therefore, it is advisable to analyze the sparsity pattern of $DG$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${{\,\mathrm{\textsf{D}}\,}}\textbf{G}$$\end{document} before the optimization procedure begins and to only compute the nonzero values of $DG$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${{\,\mathrm{\textsf{D}}\,}}\textbf{G}$$\end{document} at each iteration. Also note that the entries of $DG$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${{\,\mathrm{\textsf{D}}\,}}\textbf{G}$$\end{document} are continuous w.r.t ${\theta }_{\Lambda }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{\varvec{\Lambda }}$$\end{document} , since they are sums of products of entries of the inverse ${(I-\Lambda )}^{-1}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$(\textbf{I}- \varvec{\Lambda })^{-1}$$\end{document} , which is continuous.

3. Discussion

We have shown that separable nonlinear least squares is applicable to generalized least squares estimation of structural equation models formulated in the RAM notation. We have also shown a connection to path tracing rules in the form of trek rules. Note that when the same weight matrix is used, the point estimates obtained by SNLLS and LS are identical. Therefore, standard errors and test statistics are obtained using the same methods available for regular least squares estimation. In the following, we would like to discuss the two major benefits of using SNLLS for SEM: better convergence properties and a reduction in the computation time for parameter estimation.

3.1. Convergence

An important issue in SEM is convergence problems, especially in small samples (De Jonckere & Rosseel, 2022). If the optimizer fails to converge, no parameter estimates can be obtained. Using the SNLLS objective should lead to fewer convergence problems than LS, since only the directed parameters need to be estimated iteratively. Therefore, only the subset of directed parameters requires starting values. In many models, most of the directed parameters are factor loadings, and we can obtain very good starting values for them with the FABIN 3 estimator (Hägglund, 1982). Also, Ruhe and Wedin (1980) and Golub and Pereyra (2003) give additional proofs and reasons for why the reduced optimization problem of SNLLS should in principle be better behaved than the full LS problem. Additionally, for the class of models without unknown directed parameters, convergence problems should be eliminated altogether, as the estimator of the mean and (co)variance parameters can be computed analytically. Most prominently, this features many types of latent growth curve models.

To investigate the convergence properties of SNLLS in SEM, we ran a small simulation. We used the model in Fig. 2 to draw 1000 random data sets for varying sample sizes (N = 10 to N = 100) under the assumption of multivariate normality with zero expectation and the model-implied covariance induced by the parameters. The sample size and the factor loadings are deliberately chosen to be small to achieve a setting where non-convergence often occurs. We fitted the true model to each sample with generalized least squares (GLS; Bollen, 1989) and SNLLS estimation. All analyses were done in the programming language R (R Core Team, 2021). For GLS estimation, we used lavaan (Rosseel, 2012). The plots were created with ggplot2 (Wickham, 2016), and the data were prepared with dplyr (Wickham et al., 2021). In Fig. 3 we report the number of converged models for each sample size. In Fig. 4, we report the median number of iterations needed until convergence for each sample size. Using SNLLS effectively halved the median number of iterations until convergence for most sample sizes and more than halved the number of non-converged models for most sample sizes. This indicates that SNLLS might be a useful alternative for applied researchers to consider if they encounter convergence problems.

Figure 2. The structural equation model used to compare convergence properties of SNLLS and GLS estimation, with two latent variables, ${\zeta }_1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\zeta _{1}$$\end{document} and ${\zeta }_2$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\zeta _{2}$$\end{document} . Variances are omitted in this representation. The population values are the same as in De Jonckere and Rosseel (2022): ${\lambda }_1={\lambda }_4=1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\lambda _{1} = \lambda _{4} = 1$$\end{document} , ${\lambda }_2={\lambda }_5=0.8$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\lambda _{2} = \lambda _{5} = 0.8$$\end{document} , ${\lambda }_3={\lambda }_6=0.6$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\lambda _{3} = \lambda _{6} = 0.6$$\end{document} , $\beta =0.25$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\beta =0.25$$\end{document} , and all error variances are set to 1

Figure 3. Simulation results—number of converged replications out of 1000. GLS, generalized least squares; SNLLS, separable nonlinear least squares

Figure 4. Simulation results—median number of iterations by sample size. GLS, generalized least squares; SNLLS, separable nonlinear least squares

3.2. Computation Time

The benefits of SNLLS estimation, specifically the reduced dimensionality of the parameter space, better starting values and fewer iterations to convergence, could lead to reduced computation times. However, the computation of the SNLLS objective function and gradient is also more costly, so the cost per iteration can be higher. In sum, the question whether SNLLS estimation is faster in actual time spent in the optimization hinges upon several aspects, such as the actual implementation of the gradient, meta-parameters of the optimizer and model complexity.

Kreiberg et al. (2021) stated that estimation by SNLLS will typically be multiple times faster than LS when the reduced parameter space is much smaller than the original one. They conducted a simulation study, where they fitted a number of CFA models and concluded that the estimation time is bigger for LS than for SNLLS as the number of estimated parameters increases. Even though their simulation is useful to illustrate the potential benefits of SNLLS, it seems unfit to us to make a case for a general reduction in computation time when using SNLLS in modern software. The gradient computation in the simulation was based on a finite difference approximation in both the LS and the SNLLS condition. In existing software (Rosseel, 2012; von Oertzen et al., 2015), analytic gradients are implemented for LS estimation, so the authors compare against a straw man that would not be used in practice if computational efficiency is important. In addition, centered finite differences takes 2q calls to the objective function per computation of the gradient, where q is the number of parameters. Since SNLLS results in a smaller parameter space, their method of differentiation favors the SNLLS procedure.

It remains to implement a competitive version of SNLLS optimization for SEM using the analytic gradients derived in this paper to be able to do a realistic simulation to investigate whether SNLLS outperforms the LS estimator in practice. However, there is a large body of research concerning the efficient implementation of SNLLS (see, for example, Kaufman, 1975; O’Leary and Rust, 2013); writing competitive software for SNLLS in SEMs would be a research topic on its own. Therefore, we only give simulation results concerning the improvement of convergence rates and the number of iterations in this paper. As noted previously, for the class of models without unknown directed parameters, the estimator of the mean and (co)variance parameters can be computed in a single step. As a result, those models should especially benefit from lower computation times.

3.3. An Outlook on Maximum Likelihood Estimation

If the assumption of multivariate normality is tenable, another method of obtaining parameter estimates is maximum likelihood estimation. Here, we briefly discuss to what extent our results may have an impact on maximum likelihood optimization of SEMs. In least squares estimation with a fixed weight matrix, we saw that the undirected parameters ${\theta }_{\Omega }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{\varvec{\Omega }}$$\end{document} and the mean parameters ${\theta }_{\gamma }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _\gamma $$\end{document} enter the objective linearly. For maximum likelihood estimation, we believe it is not possible to factor out the undirected parameters (for most models used in practice). This is because the likelihood of the normal distribution

(48) $\begin{array}{c}\varphi \left(x\right)={\left({\left(2\pi \right)}^{m_{\text{obs}}},det,\Sigma \right)}^{-\frac{1}{2}}exp\left(-,\frac{1}{2},{(x-\mu )}^T,{\Sigma }^{-1},(x-\mu )\right)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \phi (x) = \left( (2\pi )^{m_\textrm{obs}} \det \varvec{\Sigma }\right) ^{-\frac{1}{2}} \exp \left( -\frac{1}{2} (x - \mu )^T \varvec{\Sigma }^{-1} (x - \mu ) \right) \end{aligned}$$\end{document}

depends on the inverse of the model-implied covariance matrix. For the simplistic example model depicted in Fig. 5, we derive the model-implied covariance matrix as

(49) $\begin{array}{c}\Sigma =\left(\begin{array}{cc}{\omega }_l+{\omega }_1& {\omega }_l\\ {\omega }_l& {\omega }_l+{\omega }_2\end{array}\right)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \varvec{\Sigma }= \begin{pmatrix} \omega _l + \omega _1 &{} \omega _l\\ \omega _l &{} \omega _l + \omega _2 \end{pmatrix} \end{aligned}$$\end{document}

and the inverse can be computed as

(50) $\begin{array}{c}{\Sigma }^{-1}={(det\Sigma )}^{-1}\text{adj}\Sigma \end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \varvec{\Sigma }^{-1} = (\det \varvec{\Sigma })^{-1} {{\,\textrm{adj}\,}}\varvec{\Sigma }\end{aligned}$$\end{document}

where $\text{adj}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${{\,\textrm{adj}\,}}$$\end{document} refers to the adjugate matrix, so in our example,

(51) $\begin{array}{c}det\Sigma =({\omega }_l+{\omega }_1)({\omega }_l+{\omega }_2)-{\omega }_l^2={\omega }_1{\omega }_l+{\omega }_2{\omega }_l+{\omega }_1{\omega }_2\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \det \varvec{\Sigma }= (\omega _l + \omega _1)(\omega _l + \omega _2) - \omega _l^2 = \omega _1\omega _l + \omega _2\omega _l + \omega _1\omega _2 \end{aligned}$$\end{document}

and

(52) $\begin{array}{c}{\Sigma }^{-1}={({\omega }_1{\omega }_l+{\omega }_2{\omega }_l+{\omega }_1{\omega }_2)}^{-1}\left(\begin{array}{cc}{\omega }_l+{\omega }_2& -{\omega }_l\\ -{\omega }_l& {\omega }_l+{\omega }_1\end{array}\right)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \varvec{\Sigma }^{-1} = (\omega _1\omega _l + \omega _2\omega _l + \omega _1\omega _2)^{-1} \begin{pmatrix} \omega _l + \omega _2 &{}\quad -\omega _l\\ -\omega _l &{}\quad \omega _l + \omega _1 \end{pmatrix} \end{aligned}$$\end{document}

Figure 5. Graph of a simplistic example model with one latent variable, measured by two indicators. The model contains no unknown directed effects and only two observed variables to allow for an easily traceable computation of the inverse of the model-implied covariance matrix. All variances are treated as unknown parameters

We see that ${\theta }_{\Omega }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{\varvec{\Omega }}$$\end{document} enters the determinant and therefore the inverse of $\Sigma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }$$\end{document} in a nonlinear way. In general, the Leibniz Formula for the determinant gives

(53) $\begin{array}{c}det\Sigma =\sum _{\pi \in S_{m_{\text{obs}}}}\text{sgn}\left(\pi \right)\prod _{i=1}^{m_{\text{obs}}}{\Sigma }_{i,\pi \left(i\right)}\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \det \varvec{\Sigma }= \sum _{\pi \in {\mathcal {S}}_{m_\textrm{obs}}} {{\,\textrm{sgn}\,}}(\pi ) \prod _{i = 1}^{m_\textrm{obs}} \varvec{\Sigma }_{i,\pi (i)} \end{aligned}$$\end{document}

where $S_{m_{\text{obs}}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {S}}_{m_\textrm{obs}}$$\end{document} denotes the symmetric group. Since this formula multiplies entries of $\Sigma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }$$\end{document} , and we saw in Eq. 30 that the entries of $\Sigma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }$$\end{document} depend on the undirected parameters, it is very likely that those form a product and enter the objective in a nonlinear way. However, for the mean parameters, the picture may be different and we leave this for future work. If the model is saturated (e.g., has zero degrees of freedom), the least squares estimates are the same as the maximum likelihood estimates, since $S=\Sigma \left({\hat{\theta }}_{\text{ML}}\right)=\Sigma \left({\hat{\theta }}_{\text{LS}}\right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{S} = \varvec{\Sigma }(\hat{\theta }_\textrm{ML}) = \varvec{\Sigma }(\hat{\theta }_\textrm{LS})$$\end{document} . Also, Lee and Jennrich (1979) showed that maximum likelihood estimation can be obtained as a form of iteratively reweighted least squares if $V$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{V}$$\end{document} is a function of the parameters:

(54) $\begin{array}{c}V=\frac{1}{2}{D}^T\left({\Sigma }^{-1},\otimes ,{\Sigma }^{-1}\right)D\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \textbf{V}= \frac{1}{2}\textbf{D}^T\left( \varvec{\Sigma }^{-1} \otimes \varvec{\Sigma }^{-1}\right) \textbf{D}\end{aligned}$$\end{document}

where $D$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{D}$$\end{document} denotes the duplication matrix from Magnus and Neudecker (2019b). Another way of obtaining ML estimates with SNLLS would therefore be to minimize the SNLLS objective and use the obtained $\Sigma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }$$\end{document} to update the weight matrix $V$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{V}$$\end{document} as given in Eq. 54. SNLLS could then be rerun with the updated weight matrix, and the weight matrix be updated again, until $\Sigma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }$$\end{document} converges to $\Sigma \left({\hat{\theta }}_{\text{ML}}\right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }(\hat{\theta }_\textrm{ML})$$\end{document} . However, we would like to note that this procedure is probably computationally quite inefficient.

3.4. Conclusion

We generalized separable nonlinear least squares estimation to all linear structural equation models that can be specified in the RAM notation, particularly those including a mean structure. We explained this result with the help of trek rules and the non-transitivity of the covariances of the error terms, providing deeper insight into the algebraic relations between the parameters of SEMs. We further derived analytic gradients and explained why they are of central importance to obtain a competitive implementation. Our simulation indicates that SNLLS leads to improvements in convergence rate and number of iterations. It remains for future research to investigate the computational costs empirically. We also showed why it is unlikely that undirected parameters enter the maximum likelihood objective linearly. Thus, another line of research could be concerned with the applicability of SNLLS to the mean parameters in maximum likelihood estimation and the relationship of SNLLS to other decomposition methods for maximum likelihood estimation (Pritikin et al., 2017, 2018). Further research might also examine whether SNLLS is applicable to multilevel models. SNLLS promises better convergence rates for least squares parameter estimation in SEM and, with an efficient implementation, also reduced computation times. This result is important in its own right but may as well serve as a first step for generating starting values for subsequent ML estimation.