jointDiag (version 0.4)

ffdiag: Joint Approximate Diagonalization of a set of square, symmetric and real-valued matrices

Description

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

Usage

ffdiag(C0, V0 = NULL, eps = .Machine$double.eps, itermax = 200, 
		keepTrace = FALSE)

Arguments

C0

DOUBLE ARRAY (KxKxN). Three-dimensional array with dimensions KxKxN representing the set of square and real-valued matrices to be jointly diagonalized. N is the number of matrices. Matrices are KxK square matrices.

V0

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 maximum number of iteration is performed, a warning appears.

keepTrace

BOOLEAN. Do we want to keep the successive estimations of the joint diagonalizer.

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.

Details

Given a set \(C_i\) of N KxK 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.

References

A. Ziehe, P. Laskov, G. Nolte and K.-R. Mueller; A Fast Algorithm for Joint Diagonalization with Non-orthogonal Transformations and its Application to Blind Source Separation; Journal of Machine Learning Research vol 5, pages 777-800, 2004

Examples

# NOT RUN {
# 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 <- ffdiag(C)$B
# B_est should be an approximate of B=solve(A)
B_est %*% A
# close to a permutation matrix (with random scales)
# }