# Solving Low-rank Factor Analysis Problem using `NExOS.jl`
**Shuvomoy Das Gupta**

In this tutorial, we will show how to solve low-rank factor analysis problem using `NExOS`. 

Factor analysis is a very popular method in statistics that reduces linear dimensionality of a particular model. The best way to understand factor analysis is to consider a generative model.

### Basic factor analysis model

A basic factor
analysis model is of the form
$$
\begin{eqnarray}
x & = & Lf+\epsilon,
\end{eqnarray}
$$

where $x\in\mathbf{R}^{n}$ is the observed random vector, $f\in\mathbf{R}^{r}$
(with $r\leq n$) is a random vector of common factor variables or
scores, $L\in\mathbf{R}^{n\times r}$ is a matrix of factor loadings,
and $\epsilon\in\mathbf{R}^{n}$ is a vector of uncorrelated random
variables. 

The observed random vector $x$ may contain series of achievement
tests, psychological evaluation, intellectual performance etc. Without
loss of generality, we assume that $x$ is mean-centered *i.e.*,
$\mathbf{E}(x)=0$, the vectors $f$ and $\epsilon$ are uncorrelated,
and the covariance matrix of $f$ is the identity matrix.

We will denote $\mathbf{cov}(\epsilon)=D=\mathbf{diag}(d_{1},d_{2},\ldots,d_{n}).$
Then the covariance matrix of $x$, denoted by $\Sigma$ can be written
as: 
$$
\Sigma=X+D,
$$
 where $X=LL^{T}$ is the covariance matrix corresponding to the common
factors. The statistical method of factor analysis involves looking
at $N$ samples generated by the model (1), *i.e.*, given $x_{1},\ldots,x_{N}$
generated by (1) we want to estimate the matrices $X$ and $D$.

### Optimization problem in consideration

The optimization problem to determine the matrices $X,D$ can be written
as: 

$$
\begin{equation}
\begin{array}{ll}
\textrm{minimize} & \left\Vert \Sigma-X-D\right\Vert _{F}^{2}\\
\textrm{subject to} & D=\mathbf{diag}(d)\\
 & d\geq0\\
 & X\succeq0\\
 & \mathbf{rank}(X)\leq r\\
 & \Sigma-D\succeq0\\
 & \|X\|_{2}\leq M,
\end{array}
\end{equation}
$$

where $X\in\mathbf{S}^{n}$ and the diagonal matrix $D\in\mathbf{S}^{n}$
with nonnegative diagonal entries are the decision variables, and
$\Sigma\in\mathbf{S}_{+}^{n}$ (a positive semidefinite matrix), $r\in\mathbf{Z}_{+},$
and $M\in\mathbf{R}_{++}$ are the problem data. The term $r$ is called the number of factors for a factor analysis model. 

A proper solution for the optimization problem above requires that
both $X$ and $D$ are positive semidefinite. Furthermore, when $\Sigma-D$,
which is the covariance matrix for the common parts of the variables,
is not positive semidefinite, that would as embarrassing as having
a negative unique variance in $D$, as noted by ten Berge [here, page 326](https://link.springer.com/chapter/10.1007/978-3-642-72253-0_44). To prevent the aforementioned undesriable situation we enforce the constraint $\Sigma-D\succeq0$.



Next, we are going to show how to solve factor analysis problem using `NExOS`, which has a built in function for this purpose. 

First, we load the packages. 


```julia
using NExOS

using LinearAlgebra, Convex, JuMP, MosekTools
```

Next, we load the data. You can change it to any other dataset as desired.


```julia

Σ = [1.0 -0.34019653769952096 -0.263030887801514 -0.14349389289304187 -0.18605860686109255; -0.34019653769952096 1.0 0.4848473200092671 0.3421745595621214 0.38218138592185846; -0.263030887801514 0.4848473200092671 1.0 0.3768343949936584 0.5028863662242727; -0.14349389289304187 0.3421745595621214 0.3768343949936584 1.0 0.3150998750134158; -0.18605860686109255 0.38218138592185846 0.5028863662242727 0.3150998750134158 1.0]

n, _ = size(Σ)

r = convert(Int64, round(rank(Σ)/2))

M = 2*opnorm(Σ ,2)
```




    4.747866256448232



Now, we construct this problem using `NExOS`.


```julia

Z0 = Σ # Initial point in Z

z0 = zeros(n) # Initial point in z

problem =  NExOS.ProblemFactorAnalysisModel(Σ, r, M, Z0, z0)
```




    ProblemFactorAnalysisModel{Array{Float64,2},Int64,Float64,Array{Float64,1}}([1.0 -0.34019653769952096 … -0.14349389289304187 -0.18605860686109255; -0.34019653769952096 1.0 … 0.3421745595621214 0.38218138592185846; … ; -0.14349389289304187 0.3421745595621214 … 1.0 0.3150998750134158; -0.18605860686109255 0.38218138592185846 … 0.3150998750134158 1.0], 2, 4.747866256448232, [1.0 -0.34019653769952096 … -0.14349389289304187 -0.18605860686109255; -0.34019653769952096 1.0 … 0.3421745595621214 0.38218138592185846; … ; -0.14349389289304187 0.3421745595621214 … 1.0 0.3150998750134158; -0.18605860686109255 0.38218138592185846 … 0.3150998750134158 1.0], [0.0, 0.0, 0.0, 0.0, 0.0])



Now, we construct the settings.


```julia
settings = NExOS.Settings(μ_max = 1, μ_min = 1e-8, μ_mult_fact = 0.5, n_iter_min = 100, n_iter_max = 100, verbose = false, freq = 50, tol = 1e-3, γ_updt_rule = :adaptive)
```




    Settings(1.0, 1.0e-8, 0.5, 100, 100, 1.0e99, 1.0e-10, 0.001, false, 50, :adaptive)




```julia
state_final = NExOS.solve!(problem, settings)
```




    StateFactorAnalysisModel{Float64,Array{Float64,2},Array{Float64,1},Int64}([0.9002146797772541 -0.5073098207885949 … -0.01162192244559973 -0.1504992835320832; -0.5073098207885949 0.6031485179871381 … 0.4175939168614036 0.5042520342917494; … ; -0.01162192244559973 0.4175939168614036 … 0.5330983234204246 0.5456366988923106; -0.1504992835320832 0.5042520342917494 … 0.5456366988923106 0.5800701078457016], [0.005250172082997496, 0.021187955387039055, 0.018516601263122615, 0.024938920128344446, 0.02242803582382049], [0.9001997771769747 -0.507340289128668 … -0.011601738359221986 -0.15049153558795636; -0.5073402889938 0.6030783235875922 … 0.41760808612095407 0.5042712190274384; … ; -0.01160173895085777 0.41760808617757567 … 0.533032993581096 0.5456672012480441; -0.1504915350227588 0.5042712188724042 … 0.5456672016930424 0.580001876846992], [0.005252094448759534, 0.021188541148083113, 0.01851724057698888, 0.02493945575366491, 0.022428602243515484], [0.9001657040093316 -0.5073953096772941 … -0.011563299637903743 -0.15047580889563483; -0.5073953095456232 0.6029405753767878 … 0.41763257751863764 0.5043046332629799; … ; -0.011563300228443932 0.4176325775757554 … 0.5329055096773659 0.5457216902191176; -0.15047580833462976 0.5043046331097865 … 0.5457216906628047 0.5798682066817082], [0.005250172082997496, 0.021187955387039055, 0.018516601263122615, 0.024938920128344446, 0.02242803582382049], 1, 7.211676530795817e-5, 1.7948927094613154e-5, 7.450580596923828e-9, 8.6316745750311e-8)



Let us take a look at the quality of the found solution.


```julia
println("fixed point gap = $(state_final.fxd_pnt_gap)")
```

    fixed point gap = 7.211676530795817e-5



```julia
println("feasibility gap = $(state_final.fsblt_gap)")
```

    feasibility gap = 1.7948927094613154e-5



```julia
println("objective value of the found solution = $(norm(vec(Σ-state_final.X-diagm(state_final.x)),2)^2)")
```

    objective value of the found solution = 0.9539784868101482



```julia
U_3, svs_3, Vt_3 = svd(state_final.X)

U_4, svs_4, Vt_4 = svd(Σ-diagm(state_final.x))

explained_variance = sum(svs_3[1:r])/sum(svs_4[1:n])

println("explained variance of the found solution = $explained_variance")
```

    explained variance of the found solution = 0.6661498849182806



```julia

```

353d67c8a02621c182de83ae6369940906735bab
ipynb
Jupyter Notebook
tutorials/Low-rank_factor_analysis_using_NExOS.ipynb
Shuvomoy/NExOS.jl
2058a74915bf56a2baea485ea3bca477bd241209
MIT
2020-11-06T14:10:39.000Z
2021-12-19T16:50:43.000Z
tutorials/Low-rank_factor_analysis_using_NExOS.ipynb
Shuvomoy/NExOS.jl
2058a74915bf56a2baea485ea3bca477bd241209
MIT
2020-11-12T10:37:12.000Z
2021-12-19T14:05:25.000Z
tutorials/Low-rank_factor_analysis_using_NExOS.ipynb
Shuvomoy/NExOS.jl
2058a74915bf56a2baea485ea3bca477bd241209
MIT
Qwen/Qwen-72B
1. YES
2. YES