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.
A matrix of k stacked pxp matrices with dimension c(kp,p) or an array with dimension c(p,p,k).
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.
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.
Author(s)
Jean-Francois Cardoso. Ported to R by Klaus Nordhausen. C code by Jari Miettinen
References
Cardoso, J.-F. and Souloumiac, A., (1996), Jacobi angles for simultaneous diagonalization, SIAM J. Mat. Anal. Appl., 17, 161–164.
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)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(JADE)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/JADE/rjd.Rd_%03d_medium.png", width=480, height=480)
> ### Name: rjd
> ### Title: Joint Diagonalization of Real Matrices
> ### Aliases: rjd frjd rjd.fortran
> ### Keywords: array
>
> ### ** 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
[,1] [,2] [,3]
[1,] 0.7443232 -0.42436171 -0.5156551
[2,] 0.3228946 0.90457863 -0.2783461
[3,] 0.5845700 0.04067724 0.8103229
> round(U %*% res.matrix$V, 4) # should be a signed permutation
[,1] [,2] [,3]
[1,] 0 1 0
[2,] -1 0 0
[3,] 0 0 1
> # matrix if V is correct.
>
> round(res.matrix$D, 4)
[,1] [,2] [,3]
[1,] 0.4397 0.0000 0.0000
[2,] 0.0000 0.8709 0.0000
[3,] 0.0000 0.0000 0.9855
[4,] 0.7055 0.0000 0.0000
[5,] 0.0000 0.1703 0.0000
[6,] 0.0000 0.0000 0.6199
[7,] 0.8273 0.0000 0.0000
[8,] 0.0000 0.3814 0.0000
[9,] 0.0000 0.0000 0.9080
[10,] 0.6609 0.0000 0.0000
[11,] 0.0000 0.2565 0.0000
[12,] 0.0000 0.0000 0.3720
>
> # 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)
, , 1
[,1] [,2] [,3]
[1,] 0.4397 0.0000 0.0000
[2,] 0.0000 0.8709 0.0000
[3,] 0.0000 0.0000 0.9855
, , 2
[,1] [,2] [,3]
[1,] 0.7055 0.0000 0.0000
[2,] 0.0000 0.1703 0.0000
[3,] 0.0000 0.0000 0.6199
, , 3
[,1] [,2] [,3]
[1,] 0.8273 0.0000 0.000
[2,] 0.0000 0.3814 0.000
[3,] 0.0000 0.0000 0.908
, , 4
[,1] [,2] [,3]
[1,] 0.6609 0.0000 0.000
[2,] 0.0000 0.2565 0.000
[3,] 0.0000 0.0000 0.372
>
> res.array.C <- frjd(X.array)
> round(res.array.C$D, 4)
, , 1
[,1] [,2] [,3]
[1,] 0.4397 0.0000 0.0000
[2,] 0.0000 0.8709 0.0000
[3,] 0.0000 0.0000 0.9855
, , 2
[,1] [,2] [,3]
[1,] 0.7055 0.0000 0.0000
[2,] 0.0000 0.1703 0.0000
[3,] 0.0000 0.0000 0.6199
, , 3
[,1] [,2] [,3]
[1,] 0.8273 0.0000 0.000
[2,] 0.0000 0.3814 0.000
[3,] 0.0000 0.0000 0.908
, , 4
[,1] [,2] [,3]
[1,] 0.6609 0.0000 0.000
[2,] 0.0000 0.2565 0.000
[3,] 0.0000 0.0000 0.372
>
>
>
>
>
> dev.off()
null device
1
>