jSVD: Singular-value decomposition of a matrix using Jacobi plane...

In jacobi: Jacobi Rotations

Usage

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

Arguments

x
eps1
eps2
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 (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.