PARAFAC WITH ORTHOGONALITY IN ONE MODE AND APPLICATIONS IN
DS-CDMA SYSTEMS
Mikael Sørensen
1 ∗, Lieven De Lathauwer
2†, Luc Deneire
1 1Laboratoire I3S, CNRS/UNSA, Les Algorithmes - Euclide-B 06903 Sophia Antipolis, France 2K.U.Leuven - E.E. Dept. (ESAT) - SCD-SISTA, Kasteelpark Arenberg 10, B-3001 Leuven-Heverlee, Belgiumsorensen@i3s.unice.fr , Lieven.DeLathauwer@kuleuven-kortrijk.be , deneire@i3s.unice.fr ABSTRACT
Blind deterministic receivers for DS-CDMA systems based on the PARAFAC model have been proposed in several papers since their conception in [1]. In many cases, the transmitted signals can be considered uncor-related. Hence, we develop PARAFAC receivers for uncorrelated signals. We introduce several numerical algorithms for orthogonality constrained PARAFAC on which receivers for uncorrelated signals can be based. Simulation results show an increase in performance when the PARAFAC receiver takes the uncorrelated-ness of the transmitted signals into account.
Index Terms— PARAFAC, DS-CDMA, blind signal separation.
1. INTRODUCTION
A blind receiver for DS-CDMA systems was proposed in [1]. The receiver takes a deterministic approach based on the multilinear PARAFAC model [2]. By exploiting the uniqueness properties of the PARAFAC model blind equalization of DS-CDMA systems was made possible. Later on, a blind PARAFAC receiver that allows for more users was proposed in [3]. It was based on a link between the PARAFAC model and simultaneous matrix diagonalization, as explained in [4].
In many cases the symbol data will be uncorrelated. Hence, the first purpose of this paper is to incorpo-rate the data uncorrelatedness assumption into existing blind PARAFAC receivers for DS-CDMA signals.
Second, in order to implement blind PARAFAC re-ceivers for uncorrelated signals we will also propose orthogonality constrained versions of the Simultane-ous matrix Diagonalization PARAFAC (SD-PARAFAC) method [4] and the Simultaneous Schur Decomposition
∗The work of Mikael Sørensen is supported by the EU by a Marie-Curie Fellowship (EST-SIGNALprogram : http://est-signal.i3s.unice.fr) under contract No MEST-CT-2005-021175.
†Research supported by: (1) Research Council K.U.Leuven: GOA-Ambiorics, CoE EF/05/006 Optimization in Engineering (OPTEC),CIF1, STRT1/08/23, (2) F.W.O.: (a) project G.0321.06, (b) Re-search Communities ICCoS, ANMMM and MLDM, (3) the Belgian Federal Science Policy Office: IUAP P6/04 (DYSCO, “Dynamical sys-tems, control and optimization”, 2007–2011), (4) EU: ERNSI.
PARAFAC (SSD-PARAFAC) method [5]. We will also discuss an orthogonality constrained Alternating Least Squares PARAFAC (ALS-PARAFAC) method [2, 1].
The paper is organized as follows. In the rest of the introduction the notations will be presented, followed by a quick review of the DS-CDMA signal model and PARAFAC receiver. Section 2 will propose orthogonal-ity constrained PARAFAC methods capable of incorpo-rating the uncorrelatedness assumption of the data. In section 3 some computer results will be reported and we will conclude in section 4.
1.1. Notations
Vectors, matrices and tensors are denoted by lower case boldface, upper case boldface and upper case calligraphic letters, respectively. Let ◦, ⊗ and de-note the outer, Kronecker and Khatri-Rao product, respectively, and let (·)T, (·)∗, (·)H
, Re{·}, · F, Col (·)
denote the transpose, conjugate, conjugate-transpose, real part, Frobenius norm and the column space of a matrix, respectively. Let triu (·) and diag (·) denote the operator that sets the strictly lower triangular el-ements and off-diagonal elel-ements of a matrix equal to zero, respectively. Moreover, let A ∈ CI×J, then
Vec (A) ∈ CIJ will denote a column vector with the property (Vec (A))i+(j−1)I = (A)ij. Let a ∈ CIJ, then the
reverse operation is Unvec (a) = A ∈ CI×J such that
(a)i+(j−1)I = (A)ij. Let A ∈ CI×J, then Dk(A) ∈ CJ×J
de-notes the diagonal matrix holding row k of A on its diagonal. Given a matrix A∈ CI×J, then a
jdenotes the jth column vector of A. Finally, IR ∈ CR×Rdenotes the
identity matrix.
1.2. DS-CDMA signal model and PARAFAC receiver In DS-CDMA systems the transmitted encoded symbol
skr from user r at the kth symbol period is given by
yjkr = hjrskr, where hjrdenotes the jth chip of the applied
spreading code. Assume there are I receive antennas and R users transmitting over a flat fading channel with a gain airbetween the ith receiver and the rth transmitter.
4142
Then the output of the ith antenna is xijk= R r=1 airyjkr = R r=1 airhjrskr.
Let us stack the channel gains, spreading codes and the transmitted symbols in the matrices A∈ CI×R, H ∈
CJ×Rand S∈ CK×R. The received data can then be stored
in the tensorX ∈ CI×J×Ksuch that
X =
R
r=1
ar◦ hr◦ sr, (1)
where ar, hrand srdenote the rth column vector of A, H
and S, respectively. Using guard chips, the model still holds when reflections take place in the far field of the ar-ray; hrthen consists of the convolution of the rth
spread-ing sequence and the rth channel impulse response [1]. The decomposition (1) ofX is sometimes referred to as the PARAFAC decomposition. If R is the minimal value for which (1) holds, then it is called the rank ofX. It can can be shown that under certain conditions the decom-position (1) is unique up to a permutation and scaling of the columns of the matrices A, H and S [4, 6]. We say that the decomposition (1) is essentially unique.
Due to the essential uniqueness properties of the PARAFAC model it is possible to blindly estimate the model parameters, i.e., solely based on the observation data, up to the mentioned ambiguities. This is the basis of the PARAFAC receiver proposed in [1].
2. PARAFAC WITH ORTHOGONALITY CONSTRAINT IN ONE MODE
The transmitted signals can usually be considered un-correlated and for sufficiently large K we have SHS ≈
αIR,α ∈ R. Hence, in this paper we will incorporate this
knowledge in the PARAFAC receiver.
When R≤ min (IJ, K) and R(R − 1) ≤ 12I(I− 1)J(J − 1),
then in subsection 2.1 we will show how to incorpo-rate an orthogonality constraint in one mode in the SD-PARAFAC method. Next, when I, K ≥ R or J, K ≥ R, then in subsection 2.2 we will show how to incorporate an orthogonality constraint in one mode in the SSD-PARAFAC method. The orthogonality constrained ALS procedure will be discussed in subsection 2.3. The lat-ter can be used to refine the results obtained with the SD-PARAFAC or SSD-PARAFAC procedures.
2.1. SD-PARAFAC with orthogonality constraint in one mode
ConsiderX ∈ CI×J×K in (1). Let X(i) ∈ CJ×K denote the
matrix such that (X(i))jk = (X)ijk, then X(i)= HDi(A) ST.
Stack the matrices{X(i)}I
i=1into the matrix
CIJ×K X (1) ⎡ ⎢⎢⎢⎢ ⎢⎢⎢⎢ ⎢⎣ X(1) ... X(I) ⎤ ⎥⎥⎥⎥ ⎥⎥⎥⎥ ⎥⎦= ⎡ ⎢⎢⎢⎢ ⎢⎢⎢⎢ ⎣ HD1(A) ... HDI(A) ⎤ ⎥⎥⎥⎥ ⎥⎥⎥⎥ ⎦S T= (A H) ST.
Let X(1) = UΣVHdenote the Singular Value
Decom-position (SVD), where U∈ CIJ×R, V∈ CK×Rare
column-wise orthonormal matrices andΣ ∈ RR×Ris a positive
diagonal matrix. Assume that the tensorX ∈ CI×J×Kof
rank R is constructed from the column-wise orthonormal symbol matrix S, and the two other matrices A and H. Assume also that (A H) has full column rank, which is a necessary condition for the essential uniqueness of the PARAFAC decomposition [7]. This assumption is true under mild conditions when R≤ min (IJ, K) [4].
Since S is of full column rank by definition, we have that Col (UΣ) = Col (A H) ST= Col (A H). Hence, there exists a nonsingular matrix F ∈ CR×R such that A H = UΣF. This result together with the relation X(1) = (A H) ST = UΣVH implies that ST = F−1VH.
Since S and V are column-wise orthonormal, we get VHV= FSTS∗FH = IR⇔ FFH = IRand hence F must be
a unitary matrix.
In order to estimate F, the SD-PARAFAC method was proposed in [4]. The problem amounts to simulta-neously diagonalizing a set matrices by congruence. In the case that F is unitary, the problem consists of maxi-mizing the cost function
f (F)= R
k=1
diag FHM(k)F2F, M(k)T= M(k),
see [4] for the construction of{M(k)} ⊂ CR×R. A complex
symmetric variant of JADE [8], known as simultaneous Takagi factorization [9], will be applied to find F.
Once F has been estimated, the unitary matrix S fol-lows immediately from S = V∗F∗. The matrices A and Hfollow from a set of decoupled SVD problems since Unvec ((UΣF)r)= Unvec (ar⊗ hr)= hraTr is a rank-1
ma-trix. In other words, hrand arcan be estimated from the
best rank-1 approximation of Unvec ((UΣF)r), r∈ [1, R].
2.2. SSD-PARAFAC with Orthogonality Constraint in One Mode
When I, K ≥ R or J, K ≥ R, we can consider the orthog-onality constrained PARAFAC estimation problem as a Simultaneous Schur Decomposition (SSD) problem [5]. We will restrict the discussion to the case I, K ≥ R, but similar results can be obtained when J, K ≥ R.
Let the tensor X ∈ CI×J×K of rank R be constructed
from a column-wise orthonormal symbol matrix S with
R≤ K, and the two other matrices A and H. Let X(j) ∈
CI×K denote the matrix such that (X(j))
ik = (X)ijk, then
X(j) = ADj(H) ST. Stack the matrices {X(j)}Jj=1 into the
matrix CIJ×K X (2) ⎡ ⎢⎢⎢⎢ ⎢⎢⎢⎢ ⎢⎣ X(1) ... X(J) ⎤ ⎥⎥⎥⎥ ⎥⎥⎥⎥ ⎥⎦= ⎡ ⎢⎢⎢⎢ ⎢⎢⎢⎢ ⎣ AD1(H) ... ADJ(H) ⎤ ⎥⎥⎥⎥ ⎥⎥⎥⎥ ⎦S T = (H A) ST. 4143
Let X(2) = (H A) ST = UΣVH denote the SVD.
Again, assuming PARAFAC essential uniqueness, we can assume that (H A) has full column rank. Hence Col (UΣ) = Col (H A) which implies the existence of a nonsingular matrix F∈ CR×Rsuch that H A = UΣF.
Due to the relations SHS= IR, X(2)= (H A) ST= UΣVH
and H A = UΣF we obtain the factorization S = V∗F∗ with F being a unitary matrix.
The first step of the SSD-PARAFAC method is to find the matrix V from the SVD of X(2). Let X(2) = X(2)V =
(H A) FH and let A = QR be the QR factorization of A, where Q ∈ CI×I is a unitary matrix and R ∈ CI×R
is an upper triangular matrix. Then an estimate of the unknowns will be obtained from the cost function
f (Q, R, H, F) = J j=1 X(j)− QRDj(H) FH2F = J j=1 QH X(j)F− RD j(H)2F.
To simplify the problem, we will in the first step of the SSD-PARAFAC method estimate the unitary matrices Q and F that make the matrices{ X(j)} as upper triangular as possible. Formally, we maximize
g (Q, F) = J j=1 triu QH X(j)F 2 F.
This approach can be seen as an orthogonality con-strained version of the SSD-PARAFAC method pre-sented in [5]. For the estimation of Q and F we will apply the extended QZ method [10].
Once the matrices Q and F have been found, the upper triangular matrix R and the diagonal matrices {Dj(H)} can be found from a set of decoupled rank-1
approximation problems as discussed next. Let R(j) = triu
QH X(j)F
, then the problem of estimating R and {Dj(H)} is equivalent to the minimization of
h (R, H) = J j=1 R(j)− RD j(H)2F= R r=1 R(r)− rrd (r)T 2 F, where R(r) =r(1)r , · · · , r(J)r , d(r)=d(1)r , · · · , d(J)r T , r(j)r and
rr are the rth column of the upper triangular part of
R(j)and R, respectively, and d(j)r is the rth entry of the
diagonal part of Dj(H).
Let R(r) = UΣVH denote the SVD of R(r), then rr =
σ1u1and d (r)
= v∗
1, whereσ1denotes the largest singular
value of R(r)and u1and v1denote its corresponding left
and right singular vector, respectively.
2.3. ALS-PARAFAC with Orthogonality Constraint in One Mode
The conventional ALS method attempts to minimize the cost function
f (A, H, S) = X(1)− (A H) ST2F.
The main idea behind the ALS method is to update a subset of the parameters{A, H, S} conditioned on previ-ously obtained estimates of the remaining parameters. When H and S are fixed, then the least squares estimate of A will be used as the conditional update of A. Due to multilinearity similar conditional updates of H and Sare applied. When the symbol matrix S is column-wise orthonormal, the ALS method reduces to a cyclic two-step procedure as explained next. A least squares estimate of S can be found from the expression
g (S) = X(1)− (A H) ST2F = X(1)2F+ A H2F− 2Re Tr S∗(A H)HX(1) . Let UΣVH = (A H)H
X(1) denote the SVD, then the
optimal S is S= V∗UT. This is known as the unitary Procrustes problem [11]. Now, let S be fixed, then an estimate of A and H can be obtained from a series of best rank-1 approximations, as explained in subsection 2.1. This two-step conditional updating procedure is repeated until convergence. Similar techniques have been proposed in [12], [13], [14].
3. SIMULATION RESULTS
In order to determine if imposing orthogonality on the symbol matrix will reduce the Bit Error Rate (BER), some simulations have been conducted. In all the simulations the real and imaginary elements of A ∈ CI×R and H∈
CJ×Rare drawn from a Gaussian distribution with zero
mean and unit variance. The symbols are drawn from a QPSK sequence and Additive White Gaussian Noise (AWGN) is added to the transmitted signal.
When the Frobenius norm of the difference between the estimated matrix F at iteration k and k+ 1 is less than = 10−5, then we decide that the extended QZ method has converged. Moreover, we decide that the ALS method and the simultaneous Takagi sweeping pro-cedure have converged when the value of their respec-tive cost function at iteration k and k+ 1 has changed less than = 10−5.
The first simulation will compare the SD-PARAFAC receiver [3], with the extended QZ iteration as its core iteration, against the orthogonality constrained SD-PARAFAC receiver (called SD-SD-PARAFACO). The model parameters were set to I= 4, J = 4, K = 1000, R = 8 and the SNR value is varying from 0 to 20 dB with an incre-ment of 5 dB. The mean BER over 100 trials can be seen in figure 1. A gain between 1 and 2 dB is obtained by
imposing column-wise orthonormality on the symbol matrix S.
The second simulation will compare the ALS-PARAFAC receiver [1] with the orthogonality con-strained SSD-PARAFAC receiver followed by two or-thogonality constrained ALS-PARAFAC refinement iterations (called SSD-ALS-PARAFACO). The model parameters were set to I = 4, J = 4, K = 120, R = 4 and the SNR value is varying from 0 to 20 dB with an increment of 5 dB. The QPSK sequence is constructed such that SHS= 120 · I4. The mean BER over 100 trials
can be seen in figure 2. A gain between 0.5 and 2 dB is obtained by imposing column-wise orthonormality on the symbol matrix.
0 5 10 15 20 10−3 10−2 10−1 100 SNR (dB) BER SD−PARAFAC SD−PARAFACO
Fig. 1.The mean BER when SNR is varying from 0 to 20 dB.
0 5 10 15 20 10−4 10−3 10−2 10−1 100 SNR (dB) BER ALS−PARAFAC SSD−ALS−PARAFACO
Fig. 2.The mean BER when SNR is varying from 0 to 20 dB. 4. CONCLUSION
We have presented a family of blind PARAFAC receivers for DS-CDMA capable of dealing with uncorrelated transmitted signals. It is based on orthogonality con-strained versions of the SD-PARAFAC, SSD-PARAFAC and ALS-PARAFAC methods. Computer simulations showed that a gain between 0.5 and 2 dB was ob-tained when the uncorrelatedness of the transmitted signals was taken into account by the receiver. The pro-posed techniques for orthogonality constrained tensor decomposition could also be generalized to more chal-lenging propagation schemes, starting from [15], [16]. A
technique that combines the PARAFAC structure with statistical independence, rather than uncorrelatedness, is reported in [17]. Variants of [17] can be based on results derived in the present paper.
5. REFERENCES
[1] N. D. Sidiropoulos, G. B. Giannakis, and R. Bro, “Blind PARAFAC receivers for DS-CDMA systems,” IEEE Trans. Signal Process., vol. 48, no. 3, pp. 810–823, 2000.
[2] R. A. Harshman, “Foundations of the Parafac procedure: Mod-els and conditions for an explanatory multimodal factor analy-sis,” UCLA Working Papers in Phonetics, vol. 16, pp. 1–84, 1970, http://publish.uwo.ca/ harshman.
[3] L. De Lathauwer and J. Castaing, “Tensor-based techniques for the blind separation of DS-CDMA signals ,” Signal Processing, vol. 87, pp. 322–336, 2007.
[4] L. De Lathauwer, “A link between the canonical decomposition in multilinear algebra and simultaneous matrix diagonalization,”
SIAM J. Matrix Anal. Appl., vol. 28, no. 3, pp. 642–666, 2006.
[5] L. De Lathauwer, B. De Moor, and J. Vandewalle, “Computation of the canonical decomposition by means of a simultaneous gen-eralized Schur decomposition,” SIAM J. Matrix Anal. Appl., vol. 26, no. 2, pp. 295–327, 2004.
[6] A. Stegeman and N. D. Sidiropoulos, “On Kruskal’s unique-ness condition for the CANDECOMP/PARAFAC decomposi-tion,” Linear Algebra and its Applications, vol. 420, 2007.
[7] X. Liu and N. D. Sidiropoulos, “Cram´er-Rao lower bounds for low-rank decomposition of multidimensional arrays,” IEEE
Trans. Signal Process., vol. 49, no. 9, pp. 2074–2086, 2001.
[8] J.-F. Cardoso and A. Souloumiac, “Blind Beamforming for non-Gaussian Signals,” IEEE Proceedings-F, vol. 140, no. 6, pp. 362– 370, 1993.
[9] L. De Lathauwer and B. De Moor, “On the blind separation of non-circular sources,” in Proc. EUSIPCO 2002, September 3-6 2002.
[10] A.-J. Van Der Veen and A. Paulraj, “An Analytical Constant Modulus Algorithm,” IEEE Trans. Signal Process., vol. 44, no. 5, pp. 1136–1155, 1996.
[11] R. A. Horn and C. Johnson, Matrix Analysis, Cambridge Univer-sity Press, 1985.
[12] J.M.F. ten Berge, Least squares optimization in
multivari-ate analysis, DSWO Press, Leiden, The Netherlands, 1993,
http://www.rug.nl/psy/onderzoek/onderzoeksprogrammas/pdf/ leastsquaresbook.pdf.
[13] J. Chen and Y. Saad, “On the tensor SVD and the optimal low rank orthogonal approximation of tensors,” SIAM J. Mat. Anal.
Appl., vol. 30, no. 4, pp. 1709–1734, 2009.
[14] A. Karfoul, L. Albera, and L. De Lathauwer, “Canonical de-composition of even higher order cumulant arrays for blind un-derdetermined mixture identification,” in SAM 2008, July 21-23 2008.
[15] L. De Lathauwer and A. de Baynast, “Blind deconvolution of DS-CDMA signals by means of decomposition in rank-(1,L,L) terms,” IEEE Trans. Signal Process., vol. 56, no. 4, pp. 1562–1571, 2008.
[16] D. Nion and L. De Lathauwer, “Block component model based blind DS-CDMA receivers,” IEEE Trans. Signal Process., vol. 56, no. 11, pp. 5567–5579, 2008.
[17] M. De Vos, D. Nion, S. Van Huffel, and L. De Lathauwer, “A com-bination of parallel factor and independent component analysis,” Tech. Rep. 04-09, ESAT-SISTA, K.U. Leuven, Belgium, 2009.
