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]