jSimDiag: Simultaneous diagonalization for symmetric matric using... In jacobi: Jacobi Rotations

Usage

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

Arguments

 `a` `eps` `itmax` `vectors` `verbose`

Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61``` ```##---- 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.