1 | jTucker3Block(a, dims, eps = 1e-06, itmax = 100, vectors = TRUE, verbose = TRUE)
|
a |
|
dims |
|
eps |
|
itmax |
|
vectors |
|
verbose |
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 113 114 115 116 117 118 119 120 121 122 | ##---- 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))
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.