jSimSVD: SVD of Simultaneous Diagonalized matrices In jacobi: Jacobi Rotations

Usage

 `1` ```jSimSVD(x, eps = 1e-06, itmax = 1000, vectors = TRUE, verbose = FALSE) ```

Arguments

 `x` `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``` ```##---- 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 (x, eps = 1e-06, itmax = 1000, vectors = TRUE, verbose = FALSE) { n <- dim(x)[1] m <- dim(x)[2] nmat <- dim(x)[3] itel <- 1 kkk <- diag(n) lll <- diag(m) sxx <- sum(x^2) fold <- Inf print(dim(x)) repeat { for (i in 1:(n - 1)) { if (i > m) next() for (j in (i + 1):n) { v <- matrix(0, 2, 2) for (imat in 1:nmat) { xij <- ifelse(j > m, 0, x[i, j, imat]) xjj <- ifelse(j > m, 0, x[j, j, imat]) xii <- x[i, i, imat] xji <- x[j, i, imat] v[1, 1] <- v[1, 1] + (xij^2 + xji^2) v[1, 2] <- v[2, 1] <- v[1, 2] + (xii * xji - xjj * xij) v[2, 2] <- v[2, 2] + (xii^2 + xjj^2) } u <- eigen(v)\$vectors[, 1] for (imat in 1:nmat) { xi <- x[i, , imat] xj <- x[j, , imat] x[i, , imat] <- u[2] * xi + u[1] * xj x[j, , imat] <- u[2] * xj - u[1] * xi } if (vectors) { ki <- kkk[i, ] kj <- kkk[j, ] kkk[i, ] <- u[2] * ki + u[1] * kj kkk[j, ] <- u[2] * kj - u[1] * ki } } } ss <- sum(apply(x, 3, diag)^2) fnew <- sqrt((sxx - ss)/sxx) if (verbose) cat(" Left iteration ", formatC(itel, digits = 4), "loss ", formatC(fnew, digits = 6, width = 10), "\n") for (k in 1:(m - 1)) { if (k > n) next() for (l in (k + 1):m) { v <- matrix(0, 2, 2) for (imat in 1:nmat) { xlk <- ifelse(l > n, 0, x[l, k, imat]) xll <- ifelse(l > n, 0, x[l, l, imat]) xkl <- x[k, l, imat] xkk <- x[k, k, imat] v[1, 1] <- v[1, 1] + (xkl^2 + xlk^2) v[1, 2] <- v[2, 1] <- v[1, 2] + (xll * xlk - xkk * xkl) v[2, 2] <- v[2, 2] + (xkk^2 + xll^2) } u <- eigen(v)\$vectors[, 1] for (imat in 1:nmat) { xk <- x[, k, imat] xl <- x[, l, imat] x[, k, imat] <- u[2] * xk - u[1] * xl x[, l, imat] <- u[1] * xk + u[2] * xl } if (vectors) { lk <- lll[, k] ll <- lll[, l] lll[, k] <- u[2] * lk - u[1] * ll lll[, l] <- u[1] * lk + u[2] * ll } } } ss <- sum(apply(x, 3, diag)^2) fnew <- sqrt((sxx - ss)/sxx) if (verbose) cat("Right iteration ", formatC(itel, digits = 4), "loss ", formatC(fnew, digits = 6, width = 10), "\n") if (((fold - fnew) < eps) || (itel == itmax)) break() itel <- itel + 1 fold <- fnew } return(list(d = apply(x, 3, diag), u = t(kkk), v = lll)) } ```

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