How best to fit factor loadings, mean vectors, and covariance matrices to a nonlinear system of matrix equations?

Question

I will note before starting that this question is related to my own work, in which i am aiming to extend and replicate the methods used in Attanasio et. al. (2015), "Human Capital Development and Parental Investment in India"**.

I am carrying out a large confirmatory factor analysis where the relationship between observed measurements and latent variables can be written as follows:

$$ y_{i,j} = \lambda_{i,j} \theta_j + \varepsilon_{i,j} $$

where $y_{i,j}$ is measurement $i$ on latent variable $j$, $\lambda_{i,j}$ is a factor loading, and $\varepsilon_{i,j}$ a measurement error. My aim is to estimate the joint distribution of the latent factors - the $\theta_j$. For my application i cannot assume normality of this distribution, so instead i assume their joint distribution to be a mixture of two normals. I also assume that the latent variables are independent of measurement errors and that the latter are normally distributed. First, writing the measurement system in matrix form, we have

$$ \textbf{y} = \Lambda \boldsymbol{\theta} + \boldsymbol{\varepsilon} $$

where $\textbf{y}$ is a vector of all measurements on each latent factor, $\boldsymbol{\Lambda}$ is a matrix containing all the factor loadings, $\boldsymbol{\theta}$ a vector containing all latent variables, and $\boldsymbol{\varepsilon}$ is a vector containing all the measurement errors. Using basic rules about the distribution of sums of independent random variables and the convolution operator, the theoretical joint distribution of the of the measurements, $p(\textbf{y})$ is given by:

$$p(\textbf{y}) = \tau \underbrace{\int g(\textbf{y} -\boldsymbol{\Lambda }\boldsymbol{\theta})f_A(\boldsymbol{\Lambda}\boldsymbol{\theta})d\boldsymbol{\theta}}_{p_A({\cdot})} + (1-\tau) \underbrace{ \int g(\textbf{y} -\boldsymbol{\Lambda }\boldsymbol{\theta})f_B(\boldsymbol{\Lambda}\boldsymbol{\theta})d\boldsymbol{\theta}}_{p_B({\cdot})}$$

with $f_A(\cdot)$ and $f_B(\cdot)$ being multivariate normal probability density functions representing the two mixture components of the joint distribution of the latent factors and $g(\cdot)$ is the joint density of the measurement errors. With the assumption made $p(\textbf{y})$ is also a mixture of two normals.

I then use a ML approximation of the joint distribution of measurements (assuming it is in fact a mixture of two normals) and obtain estimates of the moments of $p_A(\cdot)$ and $p_B({\cdot})$, denoted $\boldsymbol{\mu}^y_A$, $\boldsymbol{\mu}^y_B$, $\boldsymbol{\Sigma}^A_y$, $\boldsymbol{\Sigma}^B_y$, and $\tau$ where A,and B represent the first a second mixture components..

Once these are estimated, i would like to set them equal to the implied structure of the distribution given in the above equation. As such i am left with the following equalities:

$$ \tau \boldsymbol{\mu}_A^{\boldsymbol{\theta}} + (1-\tau)\boldsymbol{\mu}_b^{\boldsymbol{\theta}} = 0 \nonumber \\[5pt] \boldsymbol{\Lambda}\boldsymbol{\mu}^{\boldsymbol{\theta}}_A = \hat{\boldsymbol{\mu}}_A^{\tilde{\textbf{y}}} \nonumber \\[5pt] \label{momentconditions} \boldsymbol{\Lambda}\boldsymbol{\mu}^{\boldsymbol{\theta}}_B = \hat{\boldsymbol{\mu}}_B^{\tilde{\textbf{y}}} \\[5pt] \boldsymbol{\Lambda} \boldsymbol{\Theta}_A \boldsymbol{\Lambda}' + \boldsymbol{\Psi} = \hat{\boldsymbol{\Sigma}}_A \nonumber \\[5pt] \boldsymbol{\Lambda} \boldsymbol{\Theta}_B \boldsymbol{\Lambda}' + \boldsymbol{\Psi} = \hat{\boldsymbol{\Sigma}}_B \nonumber $$

Where $\boldsymbol{\Lambda}$ is still the matrix of factor loadings, and $\boldsymbol{\mu}_j^{\boldsymbol{\theta}}$ $\boldsymbol{\Theta}_j$ for $j \in \{ A,B \}$ are the vector of means and the covariance matrix from the corresponding component of the mixed joint distribution of the latent factors. All of the quantities on the right hand side of the above equations and $tau$ have been estimated.

The matrices $\boldsymbol{\Lambda}$, $\boldsymbol{\Theta}_j$ and $\boldsymbol{\Psi}$ have a structure implied by the measurement equation, and so my question is: can their parameters be estimated via minimum distance or any other suitable method in Matlab or STATA? If anyone knows of existing functions that can be used to do so, or how to go about this manually, then any comments or advice would be much appreciated.

Thanks in advance.

** URL to the original paper: http://egcenter.economics.yale.edu/sites/default/files/files/cdp1052.pdf

This question is linked to a previous question of mine that inquired about the first stage of this procedure - deriving and estimating the mixture of normal distribution for the observed measurements. Link:

How to calculate the sum of a "function" of a Gaussian mixture and a Gaussian variable?

Answer 1

This can be formulated as a Nonlinear SDP (Semidefinite Program) using YALMIP, and solved by YALMIP calling PENLAB as a solver.

It is a non-trivial problem. It is non-convex due to multiplication of optimization variables, specifically involving $\Lambda$ and $\tau$. If $\Lambda$ were assumed known, and $\tau$ held at a fixed value, it could be formulated and easily solved in CVX or YALMIP as a convex Linear SDP (then an overarching 1-D search on $\tau$ could be conducted, with a convex Linear SDP solved at each value of $\tau$). For the remainder of this answer, I will presume $\Lambda$ is not known, and show the formulation and solution as a nonlinear SDP.

Regardless of whether an exact solution to all equations and constraints can be found, the problem can be formulated as an optimization problem having an objective function which is a non-convex function, which could be some (weighted sums of) norm of differences between LHS and RHS of the variosu equations, subject to semidefinite constraints on covariance matrices being solved for, and any other constraints. The covariance matrix optimization variables are declared so as to be symmetric, and the semidefinite constraint imposed on them is that the matrix is positive semi-definite (PSD), which is equivalent to it being a covariance matrix, hence the semidefinite constraint constrains the matrix to be a valid covariance matrix (i.e., symmetric PSD).

Because the problem is non-convex, the solution you obtain may not be globally optimal, and may depend on the initial (starting values) you can provide for the optimization variables (unknowns). The initial values which you said in a comment you can obtain via basic estimation of the factor model, could help a lot.

Below is a sample YALMIP formulation. To enter it, you will need to have installed YALMIP http://users.isy.liu.se/johanl/yalmip/, and to run it (the optimize command), you will need to have obtained and installed PENLAB http://web.mat.bham.ac.uk/kocvara/penlab/ . To get started on YALMIP, read http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Tutorials.Basics and follow links. You can get help from YALMIP's developer at https://groups.google.com/forum/#!forum/yalmip .

The comments provide some explanations. You may wish to change the formulation of the objective function, for instance how to combine the norms of the various errors, as well as adding or deleting constraints as appropriate. But it should get you started. I chose the frobenius norm for matrices and the 2-norm for vectors. At some point you may wish to set PENLAB options to non-default values, which can be done by adding fields to sdpsettings. No guarantees that I didn't made any typos in the code.

% n is the relevant dimension, assumed already set to a numerical value
% Optimization variable declarations in YALMIP.
% Because matrices are square and not declared 'full', they are by default symmetric.
% sdpvar is declaring an optimization variable
tau = sdpvar
mu_a = sdpvar(n,1)
mu_b = sdpvar(n,1)
theta_a = sdpvar(n,n)
theta_b = sdpvar(n,n)
psi = sdpvar(n,n)
lambda = sdpvar(n,n,'diagonal')
% f is the relative weighting factor for matrix vs. vector norms, assumed already set
% mu_a_hat, mu_b_hat, sigma_a_hat, sigma_b_hat assumed already set
% >= 0 constraint on the square (symmetric) matrices are interpreted in LMI sense, i.e., constrains matrix to be positive semidefinite
% I put in constraint tau >= 0, but omit it if it doesn't belong
$ I put  in constraint diag(lambda) >= 0, but omit it if it doesn't belong
Constraints = [theta_a >= 0,theta_b >= 0,psi >= 0,diag(lambda) >= 0, tau>= 0]
% Put the unnecessary transpose on post-multiplying lambda in case lambda is ever made non-symmetric
Objective = f*norm(tau*mu_a+(1-tau)*mu_b) + f*norm(mu_a_hat-lambda*mu_a) + f*norm(mu_b_hat-lambda*mu_b) + norm(lambda*theta_a*lambda'+psi-sigma_a_hat,'fro') + norm(lambda*theta_b*lambda'+psi-sigma_b_hat,'fro')
% Specify PENLAB as the solver, and 'usex0' set to 1 means use assigned initial values
ops = sdpsettings('solver','penlab','usex0',1)
% All variables ending in _initial are the initial values to use in the optimizer and assumed already set
assign(tau,tau_initial)
assign(mu_a,mu_a_initial)
assign(mu_b,mu_b_initial)
assign(lambda,diag(vector_of_lambda_diagonal_initial))
assign(psi,psi_initial)
assign(sigma_a,sigma_a_initial)
assign(sigma_b,sigma_b_initial)
sol = optimize(Constraints,Objective,ops) % invokes the optimizer to solve the problem