Numerical tolerance for spectral decompositions of random matrices

Department of Applied Mathematics and Statistics
Center for Imaging Science
Johns Hopkins University
and
Paradigm, Inc.
Yale University
and
University of Maryland
and
North Carolina State University

---

Avanti Athreya, Michael Kane, Bryan Lewis, Zachary Lubberts, Vince Lyzinski, Youngser Park, Carey E. Priebe, Minh Tang, “Numerical tolerance for spectral decompositions of random matrices,” submitted, 2020.

Abstract

We precisely quantify the impact of statistical error in the quality of a numerical approximation to a random matrix eigendecomposition, and under mild conditions, we use this to introduce an optimal numerical tolerance for residual error in spectral decompositions of random matrices. We demonstrate that terminating an eigendecomposition algorithm when the numerical error and statistical error are of the same order results in computational savings with no loss of accuracy. We also repair a flaw in a ubiquitous termination condition, one in wide employ in several computational linear algebra implementations. We illustrate the practical consequences of our stopping criterion with an analysis of simulated and real networks. Our theoretical results and real-data examples establish that the tradeoff between statistical and numerical error is of significant import for data science.

Codes and Experiments

Simulation (run this on cis machine!)

procrustes <- function(X,Y, type = "I"){
    if(type == "C"){
        X <- X/norm(X, type = "F")*sqrt(nrow(X))
        Y <- Y*norm(X, type = "F")*sqrt(nrow(Y))
    }
    if(type == "D"){
        tX <- rowSums(X^2)
        tX[tX <= 1e-15] <- 1
        tY <- rowSums(Y^2)
        tY[tY <= 1e-15] <- 1
        X <- X/sqrt(tX)
        Y <- Y/sqrt(tY)
    }

    tmp <- t(X)%*%Y
    tmp.svd <- svd(tmp)
    W <- tmp.svd$u %*% t(tmp.svd$v)
    err <- norm(X%*%W - Y, type = "F")
    return(list(W=W,err=err))
}

bb <- 10
#bb <- seq(200,250,by=10) 
B <- cbind( c(.5, .2,.2), c(.2,.5,.2),c(.2,.2, .5) ) / bb
m <- nrow(B)
bvec <- rep(300,m)
N <- sum(bvec); #rho <- bvec/N

rR <- rep(1:m,times=bvec)
P <- B[rR,rR]
U <- irlba(P,nu=m,nv=m,tol=1e-06)$u # default is 1e-05, use 1e-06 instead.

nmc <- 50
kvec <- seq(1,20)
err <- iter <- matrix(0,nmc,length(kvec))
for (mc in 1:nmc) {
    set.seed(124+mc)
    g <- sbm.game(N,B,bvec)#; summary(g)

    ind <- 1
    for (k in 1:length(kvec)) {
        out <- irlba(g[],nu=m,nv=m,tol=2^-kvec[k])
        iter[mc,ind] <- out$iter
        Uhat <- out$u
        err[mc,ind] <- procrustes(U,Uhat)$err
        ind <- ind+1
    }
}
err.mean <- colMeans(err)
err.sd <- apply(err,2,sd) / sqrt(nmc)

df1 <- data.frame(k=kvec, mean=colMeans(err), err=apply(err,2,sd)/sqrt(nmc), type="error")
df2 <- data.frame(k=kvec, mean=colMeans(iter), err=apply(iter,2,sd)/sqrt(nmc), type="iteration")
df <- rbind(df1,df2)

bi <- 3000
b <- seq(150,250,by=10)

out <- foreach (i=b) %dopar% {
    B <- cbind( c(.5,.2,.2), c(.2,.5,.2),c(.2,.2,.5) ) / i
    m <- nrow(B)
    bvec <- rep(bi,m)
    N <- sum(bvec); #rho <- bvec/N

    rR <- rep(1:m,times=bvec)
    P <- B[rR,rR]
    U <- irlba(P,nu=m,nv=m,tol=1e-06)$u

    nmc <- 50
    kvec <- seq(1,20)
    err <- iter <- matrix(0,nmc,length(kvec))
    for (mc in 1:nmc) {
        set.seed(124+mc)
        g <- sbm.game(N,B,bvec)#; summary(g)

        ind <- 1
        for (k in 1:length(kvec)) {
            out <- irlba(g[],nu=m,nv=m,tol=2^-kvec[k])
            iter[mc,ind] <- out$iter
            Uhat <- out$u
            err[mc,ind] <- procrustes(U,Uhat)$err
            ind <- ind+1
        }
    }
    err.mean <- colMeans(err)
    err.sd <- apply(err,2,sd) / sqrt(nmc)

    df1 <- data.frame(k=kvec, mean=colMeans(err), err=apply(err,2,sd)/sqrt(nmc), type="err")
    df2 <- data.frame(k=kvec, mean=colMeans(iter), err=apply(iter,2,sd)/sqrt(nmc), type="iter")
    df <- rbind(df1,df2)
#    save(df, file=paste0("~/Dropbox/Tolerance/df-n9k-b",i,".Rbin"))
}