Constrained Least Squares for Extended Complex
Factor Analysis

arXiv:1804.00430v1 [stat.CO] 2 Apr 2018

Ahmad Mouri Sardarabadi 1 , Alle-Jan van der Veen2 and L.V.E. Koopmans1
Abstract—For subspace estimation with an unknown colored
noise, Factor Analysis (FA) is a good candidate for replacing the
popular eigenvalue decomposition (EVD). Finding the unknowns
in factor analysis can be done by solving a non-linear least
square problem. For this type of optimization problems, the
Gauss-Newton (GN) algorithm is a powerful and simple method.
The most expensive part of the GN algorithm is finding the
direction of descent by solving a system of equations at each
iteration. In this paper we show that for FA, the matrices
involved in solving these systems of equations can be diagonalized
in a closed form fashion and the solution can be found in
a computationally efficient way. We show how the unknown
parameters can be updated without actually constructing these
matrices. The convergence performance of the algorithm is
studied via numerical simulations.
Index Terms—Factor Analysis, Non-Linear Optimization, Covariance Matching

subspace while D can model the noise covariance matrix.
Given a sample covariance matrix R̂, the objective of FA is
to estimate A and D.
In this paper we use the EFA which was introduced by [9],
[2] and present a new algorithm based on the Gauss-Newton
(GN) method. By showing that a closed-form Singular Valude
Decomposition (SVD) of the Jacobian matrix can be found, we
develop an efficient algorithm that has competitive complexity
and performance.
The structure of the paper is as follows: in Sec. II we discuss
the EFA model, in Sec. III we introduce the least squares
(LS) formulation of the problem and discuss our method
and in Sec. IV we use numerical simulations to evaluate the
performance of the algorithm.

I. I NTRODUCTION
The eigenvalue decomposition (EVD) of the data covariance
matrix (a.k.a. principal component analysis) is a powerful tool
for subspace based and low-rank approximation techniques.
Without noise, the data covariance matrix is considered to
be rank-deficient, and its column span is called the signal
subspace. In the presence of additive noise an implicit assumption is that this noise is white. If this is not the case but the
noise covariance matrix is known from calibration, whitening
techniques can be used as a pre-processing step. However,
in many array processing applications this knowledge is not
available. A preferable approach is to replace the EVD by
techniques that jointly estimate the signal subspace and the
noise covariance matrix. Using FA as a substitute for EVD
in these cases was suggested by [1] and [2]. The latter also
includes several extension of the classical FA model which
were denoted as Extended Factor Analysis (EFA).
The FA model was introduced by Spearman [3] in 1904
and became an established multivariate technique mainly due
to the work of Lawley, Anderson, Jöreskog and others between
1940 and 1970 [4], [5], [6], [7], [8]. The FA model assumes
a covariance matrix R of the data under study (e.g., samples
acquired from an array of sensors) can be modeled as
R = AAH + D,

(1)

H

where
is the Hermitian transpose, A is a “tall” matrix
(AAH has low rank), and D is a positive diagonal matrix.
In terms of subspace-based techniques, A captures the signal
1 A. Mouri Sardarabadi (ammsa@astro.rug.nl) and L.V.E. Koopmans are
affiliated with Kapteyn Astronomical Institute, University of Groningen, The
Netherlands. 2 Alle-Jan van der Veen is affiliated with Delft University of
Technology, Delft, The Netherlands

II. DATA M ODEL
We assume to have access to N samples from the outputs
of P receivers such that a sample covariance matrix can be
constructed using
N
1 X
y[n]y[n]H
(2)
R̂ =
N n=1
where the P × 1 vector y[n] represents the zero mean sampled
output of the receivers. We also assume that the Extended
Factor Analysis [2] model is valid
R = E{R̂} = AAH + Ψ,

(3)

where E{.} represents the expectation of the argument, A is a
P × Q matrix with rank(A) = Q and Ψ is a structured matrix
such that Ψ = M ⊙ Ψ for some mask matrix M consisting of
zeros and ones. We assume to know the structure of Ψ (and
hence M) in advance, and define a selection matrix, S, for
the non-zero elements of M. Using this selection matrix we
can stack the non-zero elements of Ψ as ψ = ST vect(Ψ) and
similarly vect(Ψ) = Sψ. For the classical FA, M = I and
Ψ = D is diagonal. Given a sample covariance matrix we are
interested in finding A and Ψ or only one of them, depending
on the application.
This parametrization of the covariance model is not unique
and for any unitary matrix Q, AQ results in exactly the same
model. As a result, in the complex case, the problem needs
additional Q2 constraints for a unique parametrization [1], [2].
There are several ways to constrain A such as restricting it
to be lower triangular with real (positive) diagonal elements,
or forcing its columns to be orthogonal such that AH A = Γ

is a diagonal matrix. Regardless of the chosen constraints, the
total number of unique parameters for the problem is
n̂ = 2P Q + kMk1 − Q2 .

(4)
H

We show that using the constraints of the form A A = Γ
can be advantageous when solving the EFA model using the
least squares cost function. In the next section we present our
approach.
III. L EAST S QUARES

FOR

EFA

A. Gauss-Newton Method
The aim is to find the FA model for a sample covariance
R̂ matrix by minimizing the following non-linear LS cost
function:
min kR̂ − AAH − Ψk2F .

system of equations in an efficient way. One way to find
A that satisfies the constraints is to use the singular value
decomposition. Let the SVD of A be


 Γ 1/2 H

Q .
(6)
A = U0 Un
0
If A already satisfies the constraints then Q = IQ . Otherwise,
we can always multiply A from right by Q to make sure that
the constraints are satisfied. For simplicity we also define


JA = A∗ ⊗ IP (IP ⊗ A)KP,Q ,
which is the submatrix of J correspoding to derivatives with
respect to A and A∗ .
Theorem 1. Given the SVD of A as defined by (6), the
eigenvalue decomposition of JH
A JA is given by
H
JH
A JA = Ṽ Γ̃Ṽ ,

A,Ψ

Given the non-linearity of the cost function we use the GaussNewton method to find the solution similar to [10], [2].
The GN algorithm for solving this problem consists of the
following updates
θ (k+1) = θ(k) + µk δ
where δ is the direction of descent and the solution to the
following system of linear equations
JH Jδ = JH (r̂ − r),

(5)

where vect(.) is the vectorization operator which stacks the
columns of the argument matrix into a single column vector,
r̂ = vect(R̂), r = vect(AAH + D), and J = ∂r/∂θT is the
Jacobian matrix. Using the relation between the Kronecker
product and vectorization, we find the Jacobian matrix to be


∂r
= A∗ ⊗ IP (IP ⊗ A)KP,Q S ,
J=
T
∂θ
where

T
θ = vect(A)T vect(A)H ψ T ,

KP,Q is a permutation such that vect(XT ) = KP,Q vect(X)
for any P × Q matrix X, ⊗ represents the Kronecker product and ∗ represent the complex conjugate. We also have
KP,Q KQ,P = I.
The most expensive operation during the GN updates is the
solution of (5). We also know that the non-uniqueness of A
causes the Jacobian to be rank deficient by Q2 , which makes
the system singular. A closed form solution of this system
was first introduced in [2] using a block LDU decomposition. While that method is more general and also works for
weighted least squares (WLS), it does not produce a minimum
length solution and also does not clearly show how (5) is
reduced to a non-singular system. In the next section we show
that a unique minimum length solution can be found when the
constraint AH A is enforce at each iteration.
B. Model Reduction
H

In this section we show how the constraints A A = Γ
can be used to transform the problem into a reduced one
which does not suffer from singularities. We also show how
this transformation can be used to solve the Gauss-Newton

(7)

where



Γ̃ = 




Γ⊗I
I⊗Γ
Γ⊗I+I⊗Γ
0


,


is a diagonal matrix containing the eigenvalues and Ṽ is a
unitary matrix containing the eigenvectors. The matrix Ṽ can
be defined as Ṽ = [V|Z] where
V=



IQ ⊗ Un
0
(Γ1/2 ⊗ U0 )G1/2
IP Q
0
,
0
KQ,P
0
(U∗n ⊗ IQ ) (U∗0 ⊗ Γ1/2 )G1/2
and


(IQ ⊗ A)
Z=
G1/2 .
−KQ,P (A∗ ⊗ IQ )
with G = (Γ ⊗ IQ + IQ ⊗ Γ)−1 .
Proof. Carrying out the multiplications on both sides of (7)
verifies the equality, but it is important to show that Ṽ is a
unitary matrix. It is trivial to show that VH V = I using the
H
H
fact that UH
n Un = IP −Q , U0 U0 = IQ , Un U0 = 0 and the
property (A ⊗ B)H (C ⊗ D) = (AH C ⊗ BH D). What is left
to show is that both JA Z and VH Z are zero.
We can easily verify that
JA Z = [(A∗ ⊗ A) − (A∗ ⊗ A)]G = 0,
This shows that Z is in the null-space of J. We also have
ZH Z = G1/2 [(AT A∗ ⊗ IQ ) + (IQ ⊗ AH A)]G1/2
= G1/2 [(Γ ⊗ IQ ) + (IQ ⊗ Γ)]G1/2 = IQ2
which confirms that rank(Z) = Q2 . Now that we have shown
that Z is a basis for the null-space of J, it is sufficient to show
1/2
H
that VH Z = 0. Again, using UH
n A = 0 and U0 A = Γ
H
we can easily verify that V Z = 0.
Using this result we can define a new basis


V 0
V=
,
0 I
which allows us to remove the null-space of J and reformulate
the system of linear equations for solving the direction of

decent δ as
JH Jδ
H H
VV J JVV H δ
(V H JH JV)δ̃

allows us to use the economic-size SVD rather than a full one
(or use EVD of AH A).
To summarize, we use the following updates during the GN
iterations:

= JH (r̂ − r)
= VV H JH (r̂ − r)
= V H JH (r̂ − r),

where in the last equation we changed the variable δ̃ ≡ V H δ.
Because the systems are consistent we have δ = V δ̃. In the
next section we show how to find the solution to this system
of equations without actually constructing the matrices J and
V.
C. Improved Linear Solver
In [2] it was shown that we can split the direction of descent,
δ as
δ = [δ TA , δ Ψ ]T ,
and find each sub-vector sequentially. This can be done
because JA and S must be linearly independent for the (E)FA
problem to be identifiable [2]. First, using the generalized
matrix inversion-lemma [11] we find δ Ψ by solving
(I − T)δ Ψ = SH vect[P⊥ (R̂ − Ψ)P⊥ ],
H

T

(8)

T

where T = S (P ⊗ I + I ⊗ P − P ⊗ P)S, P = U0 UH
0
and P⊥ = I − P. The matrices I − T and T are positive
definite and positive semi-definite respectively. This means
that T is convergent, i.e. its spectral radius satisfies ρ(T) < 1.
This allows us to use simple iterative solvers to find the
solution in an efficient way. Once we have δ Ψ , we find δ A
by solving
H H
(VH JH
A JA V)δ̃ A = V JA (r̂ − r − Sδ Ψ ),

followed by δ A = Vδ̃ A . Using the results from Theorem 1
we have

 −1
Γ ⊗I
−1
,
(VH JH
=
I ⊗ Γ−1
A JA V)
G

which makes finding δ̃ A trivial. Similar to the the definition
of θ, for δ A we have


vect(∆A )
δA =
.
vect(∆∗A )
The matrix form of the solution, ∆A , can be found by
introducing some intermediate results:
E = R̂ − AAH − Ψ − ∆Ψ
1/2
H
∆1 = U0 [G̃ ⊙ (ΓUH
]
0 EU0 + U0 EU0 Γ)Γ
−1/2
∆2 = (I − P)EU0 Γ
∆A = ∆1 + ∆2 ,
(9)
where δ Ψ ≡ vect(∆Ψ ) is the matrix form of the direction of
descent for Ψ and G̃ = unvect(vectdiag(G2 ))1 is a Q × Q
matrix constructed by unvectorizing the diagonal elements of
G2 . By closer inspection we see that ∆1 is the update in the
current subspace of A and ∆2 is the update in its null-space.
It is also clear that we do not need to explicitly find Un which
1 Let γ = vectdiag(Γ) then G̃ = (γ1T + 1 γ T )⊙−2 where ⊙−1 is
Q
Q
the Hadamard or element-wise inversion of a matrix (i.e. for X = [xij ],
X⊙−n = [1/xn
ij ]).

Ψ(k+1) = Ψ(k) + µk ∆Ψ
A(k+1) = A(k) + µk ∆A
where ∆Ψ and ∆A are given by (8) and (9) respectively. Once
the direction of descent is found we can find the optimal µk
using the procedure described in [2] which requires solving
for the roots of a third order polynomial with real coefficients
for which closed form solutions exist.
IV. S IMULATION
We evaluate the convergence speed of various algorithms for
the classical factor analysis model. An array with P = 100
elements is simulated. The matrix A is chosen randomly with
a standard complex Gaussian distribution (i.e. each element
is distributed as CN (0, 1)) and D is chosen randomly with a
uniform distribution between 1 and 5.
For P = 100, the maximum number of sources is
Qmax = 89. We show simulation results for Q = 20, representative for low-rank cases, and for Q = 80 for high-rank
cases. Sources and noise are generated using standard unit
power complex Gaussian distributions.
We compare our algorithm (denoted as “Reduced method”)
with the classical ad-hoc alternating LS (see [12], [8], [2]
for details), the Block LDU based solvers proposed in [2],
Jöreskog [13], and the more recent CM method by [14].
Because the method presented here is based on LS we set the
weighting matrix for the Block LDU based algorithms equal
to identity. Fig. 1 shows the convergence behavior based on
the gradient.
The proposed algorithm has similar convergence properties
as the Block LDU. This is expected because these algorithms
also use the GN updates. CM and the ad-hoc method are
more sensitive to the number of sources, Q, and converge very
slowly for large Q. Also, for large Q, both Jöreskog and CM
might fail to converge for the same model but different noise
realization see Fig. 1c. While it is possible to find models for
which our algorithm also fails to converge as the problem is
non-convex and non-linear it is our experience that in those
scenarios the other methods also fail to converge with high
probability. Also for the same model (i.e. the same A and D)
our method is less sensitive to the noise realization.
V. C ONCLUSION
In this paper we have shown how constraining the columns
of A to be orthogonal in a(n) (E)FA problem allows us to
find a diagonalization of the Gauss-Newton approximation of
the Hessian. Using this diagonalization we have developed
a new non-linear least squares algorithm with competitive
convergence speed. We have also shown that the system of
equations needed for updating the noise covariance matrix, Ψ,
has a convergent matrix which allows fast iterative algorithms
to be applied.

Magnitude of the Gradient P = 100 Q = 20 N = 1000

Magnitude of the Gradient P = 100 Q = 80 N = 1000

Block-LDU ALS
Block-LDU NLLS
CM
Reduced Method
Joreskog 72
Ad hoc

10-2

10-4

10-2

10-4

10-6

10-6

10-8

10-10
100

10-8

101

Block-LDU ALS
Block-LDU NLLS
CM
Reduced Method
Joreskog 72
Ad hoc

10-10
100

102

101

# iterations

102

# iterations

(a) Convergence speed for N = 1000, Q = 20

(b) Convergence speed for N = 1000, Q = 80

Magnitude of the Gradient P = 100 Q = 80 N = 1000
10-2

10-4

10-6

10-8

Block-LDU ALS
Block-LDU NLLS
CM
Reduced Method
Joreskog 72
Ad hoc

10-10
100

101

102

# iterations

(c) Faild Convergence for different noise realization of the same model.

Fig. 1: Convergence for P = 100 sensors and varying number of samples N and sources Q.
Using simulation we demonstrated the newly proposed
algorithm has a competitive convergence properties and is
robust to different noise realizations.
The proposed algorithm does not guarantee the positive
(semi-)definiteness of Ψ. This will be address in future works.
R EFERENCES
[1] A.-J. van der Veen, A. Leshem, and A.-J. Boonstra, “Array
signal processing for radio astronomy,” in The Square Kilometre
Array: An Engineering Perspective, P. Hall, Ed.
Dordrecht:
Springer, 2005, pp. 231–249, iSBN 1-4020-3797-x. Reprinted
from Experimental Astronomy, 17(1-3),2004. [Online]. Available:
http://ens.ewi.tudelft.nl/pubs/aj04ska.pdf
[2] A. M. Sardarabadi and A. J. van der Veen, “Complex factor analysis and
extensions,” IEEE Transactions on Signal Processing, vol. 66, no. 4, pp.
954–967, Feb 2018.
[3] C. Spearman, “The proof and measurement of association between two
things,” The American Journal of Psychology, vol. 15, no. 1, pp. 72–101,
Jan 1904.

[4] D. N. Lawley, “The estimation of factor loadings by the method of
maximum likelihood.” Proceedings of the Royal Society of Edinburgh,
vol. 60, no. 01, pp. 64–82, 1940.
[5] T. W. Anderson and H. Rubin, “Statistical inference in factor analysis,”
In Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability, vol. 5, pp. 111 – 150, 1956. [Online]. Available:
http://www.sciencedirect.com/science/article/pii/0024379573900232
[6] K. G. Jöreskog, “A general approach to confirmatory maximum likelihood factor analysis,” Psychometrika, vol. 34, no. 2, pp. 183–202, 1969.
[7] D. N. Lawley and A. Maxwell, Factor analysis as a statistical method,
ser. 2nd. ed. New York: Am. Elsevier Publ., 1971.
[8] K. Mardia, J. Kent, and J. Bibby, Multivariate Analysis. Academic
Press, 1979.
[9] A. Mouri Sardarabadi and A.-J. van der Veen, “Subspace estimation
using factor analysis,” in 2012 IEEE 7th Sensor Array and Multichannel
Signal Processing Workshop (SAM), June 2012, pp. 477 –480.
[10] S. Y. Lee, “The Gauss-Newton algorithm for the Weighted Least Squares
factor analysis,” Journal of the Royal Statistical Society, vol. 27, no. 2,
June 1978.
[11] C. A. Rohde, “Generalized inverses of partitioned matrices,” Journal of
the Society for Industrial and Applied Mathematics, vol. 13, no. 4, pp.
1033–1035, 1965.

[12] T. Anderson, An Introduction to Multivariate Statistical Analysis, third
edition ed. Wiley, 2003.
[13] K. G. Jöreskog and A. S. Goldberger, “Factor analysis by generalized
least squares,” Psychometrika, vol. 37, no. 3, pp. 243–260, Sep 1972.
[Online]. Available: http://dx.doi.org/10.1007/BF02306782
[14] J.-H. Zhao, P. Yu, and Q. Jiang, “ML estimation for factor
analysis: EM or non-EM?” Statistics and Computing, vol. 18,
pp. 109–123, 2008, 10.1007/s11222-007-9042-y. [Online]. Available:
http://dx.doi.org/10.1007/s11222-007-9042-y

