# jSVD: Singular-value decomposition of a matrix using Jacobi plane... In jacobi: Jacobi Rotations

## Usage

 `1` ```jSVD(x, eps1 = 1e-06, eps2 = 1e-06, itmax = 1000, vectors = TRUE, verbose = FALSE) ```

## Arguments

 `x` `eps1` `eps2` `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``` ```##---- 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, eps1 = 1e-06, eps2 = 1e-06, itmax = 1000, vectors = TRUE, verbose = FALSE) { n <- nrow(x) m <- ncol(x) itel <- 1 mx <- 0 kkk <- diag(n) lll <- diag(m) sxx <- sum(x^2) sxm <- sqrt(sxx/(n * m)) repeat { for (i in 1:(n - 1)) { if (i > m) next() for (j in (i + 1):n) { xi <- x[i, ] xj <- x[j, ] xij <- ifelse(j > m, 0, x[i, j]) xjj <- ifelse(j > m, 0, x[j, j]) xii <- x[i, i] xji <- x[j, i] mx <- max(mx, abs(xij)/sxm, abs(xji)/sxm) v <- matrix(0, 2, 2) v[1, 1] <- xij^2 + xji^2 v[1, 2] <- v[2, 1] <- xii * xji - xjj * xij v[2, 2] <- xii^2 + xjj^2 u <- eigen(v)\$vectors[, 1] x[i, ] <- u[2] * xi + u[1] * xj x[j, ] <- 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 } } } ff <- sqrt((sxx - sum(diag(x)^2))/sxx) if (verbose) cat(" Left iteration ", formatC(itel, digits = 4), "maxel ", formatC(mx, width = 10), "loss ", formatC(ff, width = 10), "\n") for (k in 1:(m - 1)) { if (k > n) next() for (l in (k + 1):m) { xk <- x[, k] xl <- x[, l] xlk <- ifelse(l > n, 0, x[l, k]) xll <- ifelse(l > n, 0, x[l, l]) xkk <- x[k, k] xkl <- x[k, l] mx <- max(mx, abs(xkl)/sxm, abs(xlk)/sxm) v <- matrix(0, 2, 2) v[1, 1] <- xkl^2 + xlk^2 v[1, 2] <- v[2, 1] <- xll * xlk - xkk * xkl v[2, 2] <- xkk^2 + xll^2 u <- eigen(v)\$vectors[, 1] x[, k] <- u[2] * xk - u[1] * xl x[, l] <- 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 } } } ff <- sqrt((sxx - sum(diag(x)^2))/sxx) if (verbose) cat("Right iteration ", formatC(itel, digits = 4), "maxel ", formatC(mx, width = 10), "loss ", formatC(ff, width = 10), "\n") if ((mx < eps1) || (ff < eps2) || (itel == itmax)) break() itel <- itel + 1 mx <- 0 } return(list(d = diag(x), u = t(kkk), v = lll)) } ```

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