JADE (version 2.0-4)

rjd: Joint Diagonalization of Real Matrices

Description

This is an R version of Cardoso's rjd matlab function for joint diagonalization of k real-valued square matrices. A version written in C is also available and preferrable.

Usage

rjd(X, eps = 1e-06, maxiter = 100, na.action = na.fail)
frjd(X, weight = NULL, maxiter = 100, eps = 1e-06, na.action = na.fail)
frjd.int(X, maxiter = 100, eps = 1e-06)
rjd.fortran(X, weight = NULL, maxiter = 100, eps = 1e-06, na.action = na.fail)

Value

A list with the components

V

An orthogonal matrix.

D

A stacked matrix with the diagonal matrices or an array with the diagonal matrices. The form of the output depends on the form of the input.

iter

The frjd function returns also the number of iterations.

Arguments

X

A matrix of k stacked pxp matrices with dimension c(kp,p) or an array with dimension c(p,p,k). In case of frjd_int it has to be an array.

weight

A vector of length k to give weight to the different matrices, if NULL, all matrices have equal weight

eps

Convergence tolerance.

maxiter

Maximum number of iterations.

na.action

A function which indicates what should happen when the data contain 'NA's. Default is to fail.

Author

Jean-Francois Cardoso. Ported to R by Klaus Nordhausen. C code by Jari Miettinen

Details

Denote the square matrices as \(A_i\), \(i=1,\ldots,k\). The algorithm searches an orthogonal matrix \(V\) so that \(D_i=V'A_iV\) is diagonal for all \(i\). If the \(A_i\) commute then there is an exact solution. Otherwise, the function will perform an approximate joint diagonalization by trying to make the \(D_i\) as diagonal as possible.

Cardoso points out that notion of approximate joint diagonalization is ad hoc and very small values of eps make in that case not much sense since the diagonality criterion is ad hoc itself.

rjd, frjd and rjd.fortran terminate with an error in case maxiter is reach without convergence whereas frjd_int returns the current state at when maxiter is reached and does not warn about convergence problems.

References

Cardoso, J.-F. and Souloumiac, A., (1996), Jacobi angles for simultaneous diagonalization, SIAM J. Mat. Anal. Appl., 17, 161--164.

Miettinen, J., Nordhausen, K. and Taskinen, S. (2017), Blind Source Separation Based on Joint Diagonalization in R: The Packages JADE and BSSasymp, Journal of Statistical Software, 76, 1--31, <doi:10.18637/jss.v076.i02>.

Examples

Z <- matrix(runif(9), ncol = 3)
U <- eigen(Z %*% t(Z))$vectors
D1 <- diag(runif(3))
D2 <- diag(runif(3))
D3 <- diag(runif(3))
D4 <- diag(runif(3))

X.matrix <- rbind(t(U) %*% D1 %*% U, t(U) %*% D2 %*% U,
                  t(U) %*% D3 %*% U, t(U) %*% D4 %*% U)
res.matrix <- rjd(X.matrix)
res.matrix$V
round(U %*% res.matrix$V, 4) # should be a signed permutation 
                             # matrix if V is correct.

round(res.matrix$D, 4)

# compare to C version

#res.matrix.C <- frjd(X.matrix)
#res.matrix.C$V
#round(U %*% res.matrix.C$V, 4)
#round(res.matrix.C$D, 4)

X.array <- aperm(array(t(X.matrix), dim = c(3,3,4)), c(2,1,3))

res.array <- rjd(X.array)
round(res.array$D, 4)

res.array.C <- frjd(X.array)
round(res.array.C$D, 4)

res.array.C2 <- frjd.int(X.array)
round(res.array.C2$D, 4)