jadiag: Joint Approximate Diagonalization of a set of square,...

In jointDiag: Joint Approximate Diagonalization of a Set of Square Matrices

Description

This function performs a Joint Approximate Diagonalization of a set of square, symmetric and real-valued matrices.

Usage

jadiag(M, W_est0 = NULL, eps = .Machine$double.eps, itermax = 200, 
		keepTrace = FALSE)

Arguments

M	DOUBLE ARRAY (KxKxN). Three-dimensional array with dimensions KxKxN representing the set of square, symmetric and real-valued matrices to be jointly diagonalized. N is the number of matrices. Matrices are KxK square matrices.
W_est0	DOUBLE MATRIX (KxK). The initial guess of a joint diagonalizer. If NULL, an initial guess is automatically generated by the algorithm.
eps	DOUBLE. The algorithm stops when the criterium difference between two iterations is less than eps.
itermax	INTEGER. Alternatively, the algorithm stops when itermax sweeps have been performed without reaching convergence. If the maximumu number of iteration is performed, a warning appears.
keepTrace	BOOLEAN. Do we want to keep the successive estimations of the joint diagonalizer.

Details

Given a set C_i of N KxK symmetric and real-valued matrices, the algorithm is looking for a matrix B such that \forall i \in [1,N], B C_i B^T is as close as possible of a diagonal matrix.

Value

B	Estimation of the Joint Diagonalizer.
criter	Successive estimates of the cost function across sweeps.
B_trace	Array of the successive estimates of B across iterations.

Author(s)

Cedric Gouy-Pailler ([e-mail]), from the initial C code by Dinh-Tuan Pham.

References

Pham, D. & Cardoso, J.; Blind separation of instantaneous mixtures of nonstationary sources; IEEE Trans. Signal Process., 2001, 49, 1837-1848

Examples

# generating diagonal matrices
D <- replicate(30, diag(rchisq(df=1,n=10)), simplify=FALSE)
# Mixing and demixing matrices
B <- matrix(rnorm(100),10,10)
A <- solve(B)
C <- array(NA,dim=c(10,10,30))
for (i in 1:30) C[,,i] <- A %*% D[[i]] %*% t(A)
B_est <- jadiag(C)$B
# B_est should be an approximate of B=solve(A)
B_est %*% A
# close to a permutation matrix (with random scales)

Example output

               [,1]          [,2]          [,3]          [,4]          [,5]
 [1,]  1.686871e+00  5.551115e-16 -3.941292e-15  1.776357e-15  6.439294e-15
 [2,]  1.804112e-16 -4.718448e-16 -1.769418e-16  5.238865e-16  1.675743e-15
 [3,] -2.220446e-16  5.107026e-15  3.996803e-15 -1.332268e-15  3.620687e-01
 [4,] -2.775558e-15  9.992007e-16  3.214773e-01  5.995204e-15  1.387779e-15
 [5,] -1.776357e-15  7.482903e-14  1.110223e-16  8.119590e-01  2.220446e-16
 [6,] -8.881784e-16 -3.330669e-16 -1.637579e-15  1.582068e-15  2.609024e-15
 [7,]  0.000000e+00 -4.440892e-16  1.831868e-15 -5.551115e-17 -2.747802e-15
 [8,]  1.110223e-16  1.276756e-15  2.241263e-15  2.437286e-15 -4.343748e-15
 [9,]  1.443290e-15 -1.167015e+00  6.661338e-16  2.142730e-14 -3.885781e-16
[10,] -1.532108e-14 -1.332268e-15 -7.549517e-15  1.998401e-15  6.417089e-14
               [,6]          [,7]          [,8]          [,9]         [,10]
 [1,]  6.938894e-16 -5.828671e-16  5.551115e-17 -1.415534e-15  9.436896e-15
 [2,] -1.439820e-16 -2.810252e-16 -2.220446e-16  4.292811e-02 -3.400058e-16
 [3,]  5.440093e-15  1.110223e-16  6.550316e-15 -1.748601e-15  7.549517e-15
 [4,] -1.082467e-15  3.330669e-16 -2.498002e-15 -2.470246e-15  4.996004e-15
 [5,]  2.231548e-14 -8.881784e-16  1.332268e-15 -2.775558e-15  5.373479e-14
 [6,]  1.013079e-15 -6.356403e-01 -1.582068e-15 -2.720046e-15  2.636780e-15
 [7,]  3.052288e-01 -2.081668e-15  8.132384e-15 -1.498801e-15 -1.332268e-15
 [8,]  1.412238e-14 -1.370432e-15  1.387779e-14 -8.951173e-16 -3.116643e-01
 [9,]  1.360023e-15  1.276756e-15  3.108624e-15  1.332268e-15  2.386980e-15
[10,] -2.836620e-14 -3.608225e-15  4.792559e-01 -1.165734e-15  1.276756e-14

jointDiag documentation built on Jan. 8, 2021, 2:11 a.m.