# Joint Diagonalization of Real Matrices 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. ```r 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) ``` ## 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. ## Details Denote the square matrices as $A_i$, $i=1,...,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. ## Returns 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. ## References Cardoso, J.-F. and Souloumiac, A., (1996), Jacobi angles forsimultaneous diagonalization, _SIAM J. Mat. Anal. Appl._, *17*,161-164. Miettinen, J., Nordhausen, K. and Taskinen, S. (2017), Blind SourceSeparation Based on Joint Diagonalization in R: The Packages JADE andBSSasymp, _Journal of Statistical Software_, *76*, 1-31,\. ## Author(s) Jean-Francois Cardoso. Ported to `R` by Klaus Nordhausen. C code by Jari Miettinen ## Examples ```r 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) ```

rjd function