code and figures.

Author: karlrohe

output/plots/nytScree.pdf

@@ -0,0 +1,43 @@
+library(magrittr)
+n = 10000
+normDat = matrix(rnorm(n*2), ncol = 2) %>% scale()
+expDat = matrix(sample(c(-1,1), 2*n, T)*rexp(n*2)^1.3, ncol = 2) %>% scale()
+s = matrix(c(1, -2,-3,1), ncol = 2) %>% svd()
+R = s$u
+rexpDat = (expDat %*% R)
+sc = 3
+
+
+
+pdf(file = "rotInv.pdf", width = 10, height = 3.5)
+par(mfrow =c(1,3), mar = c(.1,1,4,1), pty="s")
+plot(normDat, pch = ".", xlim = sc*c(-1,1), ylim=sc*c(-1,1), axes=F, ylab = "", xlab="", 
+     main = "Gaussian is\nrotationally invariant.", cex.main = 2)
+arrows(x0 = c(-3,0),y0 = c(0,-3), x1 = c(3,0),y1 = c(0,3), code  =3, col = "grey43", lwd = 2)
+# lines(c(-99,099), c(0,0), col = "grey43", lwd = 2)
+# lines(c(0,0), c(-99,099), col = "grey43", lwd = 2)
+points(normDat, pch = ".")
+points(0,0, col = "red", cex=2, pch  =19)
+
+
+
+plot(rexpDat, pch = ".", xlim = sc*c(-1,1), ylim=sc*c(-1,1), axes=F, ylab = "", xlab="", 
+     main = "This non-Gaussian\nhas radial streaks.", cex.main = 2)
+arrows(x0 = c(-3,0),y0 = c(0,-3), x1 = c(3,0),y1 = c(0,3), code  =3, col = "grey43", lwd = 2)
+# lines(c(-99,099), c(0,0), col = "grey43", lwd = 2)
+# lines(c(0,0), c(-99,099), col = "grey43", lwd = 2)
+points(rexpDat, pch = ".")
+points(0,0, col = "red", cex=2, pch  =19)
+
+
+plot(rexpDat, pch = ".", xlim = sc*c(-1,1), ylim=sc*c(-1,1), axes=F, ylab = "", xlab="", 
+     main = "Varimax correctly\nestimates the basis.", cex.main = 2)
+v = varimax(rexpDat)
+stdBasis = matrix(c(1,0,-1,0,0,1,0,-1), byrow=T, ncol = 2)
+rb = stdBasis%*%R*3.5
+arrows(x0 = rb[c(1,3),1],y0 = rb[c(1,3),2], x1 = rb[-c(1,3),1],y1 = rb[-c(1,3),2], code  =3, col = "grey43", lwd = 2)
+# lines(x = rb[1:2,1], y = rb[1:2,2], col = "grey43", lwd = 2)
+# lines(x = rb[3:4,1], y = rb[3:4,2], col = "grey43", lwd = 2)
+points(rexpDat, pch = ".")
+points(0,0, col = "red", cex=2, pch  =19)
+dev.off()

output/plots/nyt_vsp_centered_factors.png

@@ -0,0 +1,75 @@
+# https://archive.ics.uci.edu/ml/datasets/bag+of+words
+library(data.table)
+library(Matrix)
+library(igraph)
+rm(list=ls())
+source("scripts/functions/vsp_fast.R")
+
+# download the data from https://archive.ics.uci.edu/ml/datasets/Bag+of+Words
+# uncomment this next line to download the 223 MB data file.
+el = read_delim("https://archive.ics.uci.edu/ml/machine-learning-databases/bag-of-words/docword.nytimes.txt.gz", skip = 3, delim = " ", col_names = F)
+el = as.tbl(el) %>% as.matrix
+A = spMatrix(nrow = max(el[,1]), ncol = max(el[,2]), i = el[,1], j = el[,2], x = el[,3]) %>% as("dgCMatrix")
+words = read_delim("https://archive.ics.uci.edu/ml/machine-learning-databases/bag-of-words/vocab.nytimes.txt", delim = "\n", col_names =F)$X1
+colnames(A)= words
+save(A, file = "data/nytA.RData")
+# load(file = "data/nytA.RData")
+
+vsp_factors = vsp(A, k = 50, sym = F, centering = T)
+save(vsp_factors, file = "output/data/vsp_factors_50.RData")
+# load(file = "output/data/vsp_factors_50.RData")
+
+
+# Are any of the singular vectors localized on a few documents?
+U = vsp_factors$pcs$U
+l4 = colSums(U^4)^(1/4)
+plot(l4) #yes
+lines(c(0,10000), .15*c(1,1))
+
+# localized singular vectors only highlight a small number of documents.
+#   they are often an artifact of noise in sparse graphs.
+#   https://papers.nips.cc/paper/8262-understanding-regularized-spectral-clustering-via-graph-conductance.pdf
+# the analysis below removes the singular vectors with 4th moment greater than .15.
+localized = l4>.15
+hist(U[,localized], breaks = 1000)
+hist(U[,!localized], breaks = 1000)
+
+# remove the singular values that correspond to localized
+#   singular vectors.  then, make a screeplot
+pdf(file = "output/plots/nytScree.pdf", height = 4, width = 4)
+vsp_factors$pcs$scree[which(!localized)] %>%
+  plot(main = "There is an eigengap \n at the 8th singular value",
+       xlab = "Leading singular values")
+lines(8.5*c(1,1), c(0,10^9))
+dev.off()
+
+# there is an eigengap at 8
+# so, these are the singular vectors we want to keep and rotate:
+which(!localized)[1:8]
+good = which(!localized)[1:8]
+#  you can pass this to vsp_rotate, with the old vsp object,
+#    setting k = good.
+
+all_factors = vsp_factors
+vsp_factors = vsp_keep(all_factors, keepThese = good, fast = F, recenter=T)
+
+
+png(file = "output/plots/nyt_vsp_centered_pcs.png", height = 7, width = 7, units= "in", res= 200)
+plot(vsp_factors,k=8,whichPlot = "u",nsamp = 5000)
+dev.off()
+
+
+png(file = "output/plots/nyt_vsp_centered_factors.png", height = 7, width = 7, units= "in", res= 200)
+plot(vsp_factors,k=8,whichPlot = "z",nsamp = 5000)
+dev.off()
+
+
+png(file = "output/plots/nyt_vsp_centered_pcs_WORDS.png", height = 7, width = 7, units= "in", res= 200)
+plot(vsp_factors,k=8,whichPlot = "v",nsamp = 5000)
+dev.off()
+
+
+png(file = "output/plots/nyt_vsp_centered_factors_WORDS.png", height = 7, width = 7, units= "in", res= 200)
+plot(vsp_factors,k=8,whichPlot = "y",nsamp = 5000)
+dev.off()
+