Description Usage Arguments Details Value Note Author(s) References Examples
This function performs a Joint Approximate Diagonalization of a set of square, symmetric and real-valued matrices.
1 2 |
C |
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. |
W0 |
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. |
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.
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. |
Two versions of the quadratic optimization are present in the paper
referenced below. These two
versions have different complexities, O(N K^3)
and O(K^5)
.
Currently only the version with O(N K^3)
is implemented.
Cedric Gouy-Pailler (cedric.gouypailler@gmail.com), from the initial matlab code by R. Vollgraf.
R. Vollgraf and K. Obermayer; Quadratic Optimization for Approximate Matrix Diagonalization; IEEE Transaction on Signal Processing, 2006
1 2 3 4 5 6 7 8 9 10 11 | # 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 <- qdiag(C)$B
# B_est should be an approximate of B=solve(A)
B_est %*% A
# close to a permutation matrix (with random scales)
|
[,1] [,2] [,3] [,4] [,5]
[1,] -1.898481e-14 2.512797e-11 -8.356649e-13 -2.249234e-11 -3.543247e-01
[2,] -2.181394e-11 -7.976356e-12 2.936032e-09 1.256645e-11 3.500783e-12
[3,] 3.831213e-11 2.490159e-11 -8.256923e-12 -3.415123e-01 -7.230647e-12
[4,] -1.467854e-09 2.470076e-08 -2.275125e-12 -3.068029e-11 -6.040557e-12
[5,] 3.138977e-11 2.024279e-11 -1.118830e-11 -1.066362e-11 -1.056424e-11
[6,] -1.279137e-10 -1.458761e-10 -8.254508e-14 2.775946e-12 -2.283618e-12
[7,] 3.626671e-10 -2.778966e-01 3.797934e-12 -9.948101e-12 1.594190e-12
[8,] 3.178246e-01 2.163110e-10 -5.295764e-12 -1.106515e-11 2.614436e-13
[9,] 3.233705e-11 1.237959e-11 -8.823442e-12 7.136708e-12 -3.754039e-12
[10,] -2.961187e-12 2.457309e-11 -3.378768e-01 -3.480641e-11 -3.255660e-12
[,6] [,7] [,8] [,9] [,10]
[1,] -1.083014e-11 5.106560e-11 1.574815e-11 1.524891e-12 1.624278e-11
[2,] -1.719180e-12 2.930432e-01 -6.818268e-12 5.409007e-13 -1.211953e-11
[3,] -5.819539e-12 -2.195275e-11 5.602609e-12 2.050027e-13 2.310527e-11
[4,] 4.723666e-12 -6.453182e-11 3.528640e-01 8.331114e-13 -6.612388e-09
[5,] -3.623661e-01 -4.324163e-11 -1.238102e-11 -8.759382e-13 8.679526e-12
[6,] -1.767642e-12 -9.156731e-12 2.934949e-09 1.503321e-09 3.159958e-01
[7,] 3.227488e-13 5.677538e-12 1.113340e-08 1.797874e-12 -6.973426e-11
[8,] -4.145934e-12 5.980427e-11 6.652543e-10 -2.706835e-12 9.776964e-11
[9,] 1.439515e-12 -5.245709e-11 1.238898e-12 4.422307e-01 -3.198898e-09
[10,] -4.466150e-12 2.611318e-09 4.812747e-12 1.426637e-12 1.038947e-11
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.