# jTucker3Diag: Fits an orthonormal version of the INDSCAL/PARAFAC model. In jacobi: Jacobi Rotations

## Usage

 `1` ```jTucker3Diag(a, eps = 1e-06, itmax = 100, vectors = TRUE, verbose = TRUE) ```

## 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 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112``` ```##---- 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-06, itmax = 100, vectors = TRUE, verbose = TRUE) { n <- dim(a)[1] m <- dim(a)[2] k <- dim(a)[3] nmk <- min(n, m, k) kn <- diag(n) km <- diag(m) kk <- diag(k) ossq <- 0 itel <- 1 repeat { for (i in 1:(n - 1)) for (j in (i + 1):n) { ai <- a[i, , ] aj <- a[j, , ] acc <- ass <- asc <- 0 if (i <= min(m, k)) { acc <- acc + a[i, i, i]^2 ass <- ass + a[j, i, i]^2 asc <- asc + a[i, i, i] * a[j, i, i] } if (j <= min(m, k)) { acc <- acc + a[j, j, j]^2 ass <- ass + a[i, j, j]^2 asc <- asc - a[j, j, j] * a[i, j, j] } u <- eigen(matrix(c(acc, asc, asc, ass), 2, 2))\$vectors[, 1] c <- u[1] s <- u[2] a[i, , ] <- c * ai + s * aj a[j, , ] <- c * aj - s * ai if (vectors) { ki <- kn[i, ] kj <- kn[j, ] kn[i, ] <- c * ki + s * kj kn[j, ] <- c * kj - s * ki } } for (i in 1:(m - 1)) for (j in (i + 1):m) { ai <- a[, i, ] aj <- a[, j, ] acc <- ass <- asc <- 0 if (i <= min(n, k)) { acc <- acc + a[i, i, i]^2 ass <- ass + a[i, j, i]^2 asc <- asc + a[i, i, i] * a[i, j, i] } if (j <= min(n, k)) { acc <- acc + a[j, j, j]^2 ass <- ass + a[j, i, j]^2 asc <- asc - a[j, j, j] * a[j, i, j] } u <- eigen(matrix(c(acc, asc, asc, ass), 2, 2))\$vectors[, 1] c <- u[1] s <- u[2] a[, i, ] <- c * ai + s * aj a[, j, ] <- c * aj - s * ai if (vectors) { ki <- km[i, ] kj <- km[j, ] km[i, ] <- c * ki + s * kj km[j, ] <- c * kj - s * ki } } for (i in 1:(k - 1)) for (j in (i + 1):k) { ai <- a[, , i] aj <- a[, , j] acc <- ass <- asc <- 0 if (i <= min(n, m)) { acc <- acc + a[i, i, i]^2 ass <- ass + a[i, i, j]^2 asc <- asc + a[i, i, i] * a[i, i, j] } if (j <= min(n, m)) { acc <- acc + a[j, j, j]^2 ass <- ass + a[j, j, i]^2 asc <- asc - a[j, j, j] * a[j, j, i] } u <- eigen(matrix(c(acc, asc, asc, ass), 2, 2))\$vectors[, 1] c <- u[1] s <- u[2] a[, , i] <- c * ai + s * aj a[, , j] <- c * aj - s * ai if (vectors) { ki <- kk[i, ] kj <- kk[j, ] kk[i, ] <- c * ki + s * kj kk[j, ] <- c * kj - s * ki } } nssq <- 0 for (v in 1:nmk) nssq <- nssq + a[v, v, v]^2 if (verbose) cat("Iteration ", formatC(itel, digits = 4), "ssq ", formatC(nssq, digits = 10, width = 15), "\n") if (((nssq - ossq) < eps) || (itel == itmax)) break() itel <- itel + 1 ossq <- nssq } d <- rep(0, nmk) for (v in 1:nmk) d[v] <- a[v, v, v] return(list(a = a, d = d, kn = kn, km = km, kk = kk)) } ```

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