# R/secularRpca.R In onlinePCA: Online Principal Component Analysis

```secularRpca <- function (lambda, U, x, n, f = 1/n, center, tol = 1e-10, reortho = FALSE)
{
# Dimensions check
d <- length(lambda)
stopifnot(length(x) == d)
stopifnot(ncol(U) == d)

# Argument check
stopifnot(f >= 0 && f <= 1)

# Centering
if (!missing(center)) x <- x - center

# Trivial cases
if (f == 0) {
ind <- order(lambda, decreasing=TRUE)
return(list(values=lambda[ind], vectors=U[,ind]))
}
if (f == 1) {
lambda <- c(sum(x^2), rep(0,d-1))
U <- cbind(x, matrix(runif(d^2-d),d,d-1))
U <- qr.Q(qr(U))
return(list(values=lambda, vectors=U))
}

lambda <- (1-f) * lambda
z <- sqrt(f) * crossprod(U,x)
ix <- order(lambda)
U <- U[,ix]
lambda <- lambda[ix]
z <- z[ix]

# Initial deflation
eps <- max(tol, .Machine\$double.eps)
active <- seq_len(d)
zzero <- which(abs(z) < eps)
if (length(zzero))
active <- setdiff(active, zzero)
ind <- which(diff(lambda[active]) < eps)
if (length(ind)) {
ind <- sort(unique(c(ind, ind + 1L)))
ix <- which(diff(ind) > 1)
ub <- ind[c(ix, length(ind))]
lb <- ind[c(1, ix + 1)]
ngroup <- length(ub)
group <- lapply(1:ngroup, function(i) active[lb[i]:ub[i]])
for (i in 1:ngroup) {
g <- group[[i]]
mult <- length(g)
U1 <- U[, g]
z1 <- z[g]
sigma <- sqrt(sum(z1^2))
a <- c(sigma + z1[1], z1[-1])
a <- a/sqrt(sum(a^2))
H <- diag(mult) - 2 * tcrossprod(a)
U[, g] <-U1 %*% H
z[g] <- c(-sigma, rep(0, mult - 1))
active <- setdiff(active, g[-1])
}
rm(g, U1, z1, sigma, a, H)
}

# Resolution of secular equations
pact <- length(active)
dact <- lambda[active]
z2act <- z[active]^2
bounds <- c(dact, dact[pact] + sum(z2act))
amp <- diff(bounds)
secular <- function(lambda) sum(z2act/{
dact - lambda
}) + 1
solver <- function(i) {
delta <- amp[i]/100
repeat {
lb <- bounds[i] + delta
ub <- bounds[i + 1] - delta
flb <- secular(lb)
fub <- secular(ub)
test <- flb * fub
if (test < 0) {
return(uniroot(secular, c(lb, ub), f.lower = flb, f.upper = fub,
tol = tol)\$root)
}
else if (test == 0) {
return(ifelse(flb == 0, lb, ub))
}
else if (flb > 0 && lb - bounds[i] <= tol) {
return((bounds[i] + lb)/2)
}
else if (fub < 0 && bounds[i + 1] - ub <= tol) {
return((ub + bounds[i + 1])/2)
}
else delta <- delta/10
}
}
roots <- sapply(seq_len(pact), solver)
lambda[active] <- roots

# Computation of eigenvectors
if (reortho) {
num <- matrix(roots - rep(dact, each = pact), pact, pact)
den <- matrix(dact - rep(dact, each = pact), pact, pact)
den[den == 0] <- 1
ztilde <- sqrt(apply(num/den, 2, prod))
eigvecs <- matrix(ztilde, pact, pact)/matrix(dact - rep(roots,
each = pact), pact, pact)
norms <- sqrt(.colSums(eigvecs^2, pact, pact))
eigvecs <- eigvecs/rep(norms, each = pact)
}
else {
eigvecs <- matrix(z[active], pact, pact)/matrix(dact -
rep(roots, each = pact), pact, pact)
norms <- sqrt(.colSums(eigvecs^2, pact, pact))
eigvecs <- eigvecs/rep(norms, each = pact)
}
U[, active] <- U[, active] %*% eigvecs
return(list(values = lambda[d:1], vectors = U[, d:1]))
}
```

## Try the onlinePCA package in your browser

Any scripts or data that you put into this service are public.

onlinePCA documentation built on May 2, 2019, 3:28 a.m.