jTucker3Block {jacobi}R Documentation

Fits a TUCKER-3 model using Jacobi plane rotations.

Usage

jTucker3Block(a, dims, eps = 1e-06, itmax = 100, vectors = TRUE, verbose = TRUE)

Arguments

a
dims
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, dims, 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)
    p <- dims[1]
    q <- dims[2]
    r <- dims[3]
    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 <= p) 
                for (u in 1:q) for (v in 1:r) {
                  acc <- acc + a[i, u, v]^2
                  ass <- ass + a[j, u, v]^2
                  asc <- asc + a[i, u, v] * a[j, u, v]
                }
            if (j <= p) 
                for (u in 1:q) for (v in 1:r) {
                  acc <- acc + a[j, u, v]^2
                  ass <- ass + a[i, u, v]^2
                  asc <- asc - a[j, u, v] * a[i, u, v]
                }
            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 <= q) 
                for (u in 1:p) for (v in 1:r) {
                  acc <- acc + a[u, i, v]^2
                  ass <- ass + a[u, j, v]^2
                  asc <- asc + a[u, i, v] * a[u, j, v]
                }
            if (j <= q) 
                for (u in 1:p) for (v in 1:r) {
                  acc <- acc + a[u, j, v]^2
                  ass <- ass + a[u, i, v]^2
                  asc <- asc - a[u, i, v] * a[u, j, v]
                }
            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 <= r) 
                for (u in 1:p) for (v in 1:q) {
                  acc <- acc + a[u, v, i]^2
                  ass <- ass + a[u, v, j]^2
                  asc <- asc + a[u, v, i] * a[u, v, j]
                }
            if (j <= r) 
                for (u in 1:p) for (v in 1:q) {
                  acc <- acc + a[u, v, j]^2
                  ass <- ass + a[u, v, i]^2
                  asc <- asc - a[u, v, i] * a[u, v, 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 <- kk[i, ]
                kj <- kk[j, ]
                kk[i, ] <- c * ki + s * kj
                kk[j, ] <- c * kj - s * ki
            }
        }
        nssq <- 0
        for (i in 1:p) for (j in 1:q) for (l in 1:r) nssq <- nssq + 
            a[i, j, l]^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 <- a[1:p, 1:q, 1:r]
    return(list(a = a, d = d, kn = kn, km = km, kk = kk))
  }

[Package jacobi version 1.0 Index]