jSimDiag: Simultaneous diagonalization for symmetric matric using...

In jacobi: Jacobi Rotations

Usage

jSimDiag(a, eps = 1e-10, itmax = 100, vectors = TRUE, verbose = FALSE)

Arguments

a
eps
itmax
vectors
verbose

Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (a, eps = 1e-10, itmax = 100, vectors = TRUE, verbose = FALSE) 
{
    n <- dim(a)[1]
    kk <- diag(n)
    m <- dim(a)[3]
    itel <- 1
    saa <- sum(a^2)
    fold <- saa - sum(apply(a, 3, function(x) sum(diag(x^2))))
    repeat {
        for (i in 1:(n - 1)) for (j in (i + 1):n) {
            ad <- (a[i, i, ] - a[j, j, ])/2
            av <- a[i, j, ]
            v <- matrix(0, 2, 2)
            v[1, 1] <- sum(ad^2)
            v[1, 2] <- v[2, 1] <- sum(av * ad)
            v[2, 2] <- sum(av^2)
            u <- eigen(v)$vectors[, 2]
            c <- sqrt((1 + u[2])/2)
            s <- sign(u[1]) * sqrt((1 - u[2])/2)
            for (k in 1:m) {
                ss <- s^2
                cc <- c^2
                sc <- s * c
                ai <- a[i, , k]
                aj <- a[j, , k]
                aii <- a[i, i, k]
                ajj <- a[j, j, k]
                aij <- a[i, j, k]
                a[i, , k] <- a[, i, k] <- c * ai - s * aj
                a[j, , k] <- a[, j, k] <- s * ai + c * aj
                a[i, j, k] <- a[j, i, k] <- u[1] * (aii - ajj)/2 + 
                  u[2] * aij
                a[i, i, k] <- aii * cc + ajj * ss - 2 * sc * 
                  aij
                a[j, j, k] <- ajj * cc + aii * ss + 2 * sc * 
                  aij
            }
            if (vectors) {
                ki <- kk[, i]
                kj <- kk[, j]
                kk[, i] <- c * ki - s * kj
                kk[, j] <- s * ki + c * kj
            }
        }
        fnew <- saa - sum(apply(a, 3, function(x) sum(diag(x^2))))
        if (verbose) 
            cat("Iteration ", formatC(itel, digits = 4), "old loss ", 
                formatC(fold, width = 10), "new loss ", formatC(fnew, 
                  width = 10), "\n")
        if (((fold - fnew) < eps) || (itel == itmax)) 
            break()
        itel <- itel + 1
        fold <- fnew
    }
    return(list(a = a, d <- apply(a, 3, diag), k = kk))
  }

jacobi documentation built on May 2, 2019, 4:50 p.m.