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

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.

** 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?

• I don't follow this at all, because (1) you have already estimated $\Sigma_x$ and (2) the rest is purely a linear algebraic question. Obviously you're not looking at it this way, which suggests there's more to your situation than you have disclosed. Could you clarify? – whuber Aug 18 '16 at 13:58
• Your question is not clear. Please tell us what is known and unknown, and what the certain theoreticla structures are. Does this correspond to r.v. $x = \boldsymbol{\Lambda}Y + Z$ where $Y$ has covariance $\Theta$ and $Z$ has covariance $\boldsymbol{\Psi}$ and $x$ has covariance $\boldsymbol{\Sigma}_x$? – Mark L. Stone Aug 18 '16 at 14:00
• Based on a quick look, this seems so related to your other question stats.stackexchange.com/questions/230102 that you should at least link to it – Juho Kokkala Aug 18 '16 at 18:07
• Non-trivial problem. I don't know whether you have positive degrees of freedom (multiple solutions solve exactly), or negative degrees of freedom (need to find closest fitting, in some sense, solution). Either way, can formulate objective function as a non-convex function which could be some (weighted sums of) norm of differences between LHS and RHS, subject to semidefinite constraints on covariance matrices being solved for. Use PENLAB (free) solver under YALMIP (free) under MATLAB. Ignoring $\tau$, for fixed $\Lambda$, problem may be convex, allowing simpler formulation and solution. – Mark L. Stone Aug 18 '16 at 18:28
• I don't really understand what your structure is, but am guessing you can impose constraints to account for it, or build them into your optimization problem formulation. The solution you obtain may not be globally optimal, and may depend on starting values you can provide for the unknowns. If you have a good basis for doing that, it will help a lot. – Mark L. Stone Aug 18 '16 at 18:31

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