Initial import

README.md

@@ -1,4 +1,6 @@
 randomized_nonlinear_component_analysis
 =======================================
 
-Paper and code for "Randomized Nonlinear Component Analysis", ICML (2014)
+Paper and code for "Randomized Nonlinear Component Analysis" (ICML 2014), in
+collaboration with Suvrit Sra, Alex Smola, Zoubin Ghahramani and Bernhard
+Schoelkopf.

code/rca.r

@@ -0,0 +1,72 @@
+library(kernlab)
+
+geigen <- function (A,B,C,top) {
+  p       <- nrow(B) 
+  q       <- nrow(C) 
+  s       <- min(c(p,q))
+  B       <- (B+t(B))/2
+  C       <- (C+t(C))/2
+  Bfac    <- chol(B)
+  Cfac    <- chol(C)
+  Bfacinv <- solve(Bfac)
+  Cfacinv <- solve(Cfac)
+  D       <- t(Bfacinv)%*%A%*%Cfacinv
+  if (p >= q) {
+      result <- svd(D,nu=top,nv=top)
+      values <- result$d
+      L      <- Bfacinv %*% result$u
+      M      <- Cfacinv %*% result$v
+  } else {
+      result <- svd(t(D),nu=top,nv=top)
+      values <- result$d
+      L      <- Bfacinv %*% result$v
+      M      <- Cfacinv %*% result$u
+  }
+  list(cor=values, xcoef=L, ycoef=M)
+}
+
+rcc <- function (X,Y,l1,l2,top) {
+  geigen(cov(X,Y),var(X)+diag(l1,ncol(X)),
+                  var(Y)+diag(l2,ncol(Y)),top)
+}
+
+aug <- function(x,k,type="fourier") {
+  s <- sigest(x,scaled=NULL)[2]
+
+  if(type == "linear") {
+    return(function(x0) x0)
+  }
+
+  if(type == "nystrom") {
+    w <- x[sample(1:nrow(x),k),]
+    return(function(x0) kernelMatrix(rbfdot(s),x0,w))
+  }
+
+  if(type == "fourier") {
+    w <- matrix(rnorm(ncol(x)*k,sd=sqrt(2*s)),ncol(x))
+    b <- runif(k,0,2*pi)
+    f <- function(x0) x0%*%w+t(matrix(b,k,nrow(x0)))
+    return(function(x0) cos(f(x0)))
+  }
+}
+
+rcca_fit <- function(x,y,kx,ky,type,top) {
+  augx <- aug(x,kx,type)
+  augy <- aug(y,ky,type)
+  C    <- rcc(augx(x),augy(y),1e-10,1e-10,top)
+  list(cor=sum(abs(C$cor[1:top])),a=C$xcoef,b=C$ycoef,augx=augx,augy=augy)
+}
+
+rcca_eval <- function(rcca,x,y){
+  list(x=rcca$augx(x)%*%rcca$a,y=rcca$augy(y)%*%rcca$b)
+}
+
+rpca_fit <- function(x,k,type) {
+  augx <- aug(x,k,type)
+  list(augx=augx,pca=prcomp(augx(x)))
+}
+
+rpca_eval <- function(rpca,x) {
+  predict(rpca$pca,rpca$augx(x))
+}
+

figures/plot_lupi.R

@@ -0,0 +1,28 @@
+library(tikzDevice)
+
+lupi <- rbind(
+c(55,62.45),
+c(56.85,64.55),
+c(56.25,64.55),
+c(53.9,61.65),
+c(53.05,62.15),
+c(54,62.45),
+c(55.8,60.6),
+c(56.65,65.6),
+c(54.85,61.55),
+c(55.7,61.6),
+c(54.5,61.75),
+c(56.6,62.9),
+c(54.2,61.8),
+c(54.8,64.7))
+
+tikz("lupi.tex", height=2, width=5, standAlone=TRUE)
+par(mar=c(2,4,0.5,0.5),cex=1.3)
+
+plot(lupi[,1], t="o", pch=0, ylim=c(50,70), ylab="Classification Acc.",lwd=2)
+points(lupi[,2], t="o", pch=1,lwd=2)
+
+legend("topright", c("SURF", "RCCA"), lty=rep(1,2), pch=0:1,lwd=2,bg="white")
+dev.off()
+
+system("pdflatex lupi.tex; evince lupi.pdf")

figures/plot_norms.R

@@ -0,0 +1,109 @@
+library(tikzDevice)
+
+options( tikzLatexPackages = c(
+  getOption( "tikzLatexPackages" ),
+  "\\usepackage{bm,times}"
+))
+
+my_mar <- c(3.3,3.3,1,1)
+my_mgp <- c(2.1,1,0)
+my_cex <- 1.5
+
+r <- rbind(c(1e-04,5.792990e+04),
+           c(1e-03,6.431905e+03),
+           c(1e-02,6.014578e+02),
+           c(1e-01,6.245001e+01),
+           c(1e+00,6.651322e+00),
+           c(1e+01,7.249647e-01),
+           c(1e+02,8.783448e-02),
+           c(1e+03,1.341658e-02))
+
+tikz("norm_cca_g.tex",standAlone=TRUE,width=3,height=3)
+par(mar=my_mar, mgp=my_mgp,cex=my_cex) 
+plot(log(r),t="o",xlab="$\\log(\\gamma)$",ylab="$\\log(\\|\\hat{\\bm R}^{-1}\\hat{\\bm L}-\\bm R^{-1}\\bm L\\|$)",lwd=3) 
+dev.off()
+
+system("pdflatex norm_cca_g.tex")
+
+r <- rbind(c( 100, 60.574517, 19657.962),
+           c( 250, 38.365135, 11739.667),
+           c( 500, 27.343299,  8340.034),
+           c(1000, 18.891724,  5904.768),
+           c(1500, 14.984350,  4727.440),
+           c(2000, 14.091385,  4141.354),
+           c(2500, 12.474679,  3709.368),
+           c(3000, 11.025231,  3372.574),
+           c(4000,  9.240461,  2910.546),
+           c(5000,  8.677846,  2603.916))
+
+tikz("norm_cca_m.tex",standAlone=TRUE,width=3,height=3);
+par(mar=my_mar, mgp=my_mgp,cex=my_cex) 
+plot(r[,1],r[,3],t="o",xlab="$m$",ylab="$\\|\\hat{\\bm R}^{-1}\\hat{\\bm L}-\\bm R^{-1}\\bm L\\|$",lwd=3) 
+dev.off()
+
+system("pdflatex norm_cca_m.tex")
+
+r <- rbind(c( 100,  1.961450,   393.3971),
+           c( 250,  4.850861,  1295.4472),
+           c( 500,  9.189258,  2782.4135),
+           c(1000, 21.171350,  5834.5632),
+           c(1500, 32.406314,  9154.8163),
+           c(2000, 34.329282, 12056.4939),
+           c(2500, 47.685869, 15105.3356),
+           c(3000, 63.446192, 17977.8319),
+           c(4000, 71.370321, 24444.7370),
+           c(5000, 92.434716, 30033.6785))
+tikz("norm_cca_n.tex",standAlone=TRUE,width=3,height=3);
+par(mar=my_mar, mgp=my_mgp,cex=my_cex) 
+plot(r[,1],r[,3],t="o",xlab="$n$",ylab="$\\|\\hat{\\bm R}^{-1}\\hat{\\bm L}-\\bm R^{-1}\\bm L\\|$",lwd=3)
+dev.off()
+
+system("pdflatex norm_cca_n.tex")
+
+r <- rbind(c( 100, 60.574517, 19657.962),
+           c( 250, 38.365135, 11739.667),
+           c( 500, 27.343299,  8340.034),
+           c(1000, 18.891724,  5904.768),
+           c(1500, 14.984350,  4727.440),
+           c(2000, 13.591385,  4141.354),
+           c(2500, 12.474679,  3709.368),
+           c(3000, 11.025231,  3372.574),
+           c(4000,  9.240461,  2910.546),
+           c(5000,  8.677846,  2603.916))
+fff <- 
+c(60.623741,
+  38.305267,
+  27.056785,
+  19.102907,
+  15.579209,
+  13.478666,
+  12.045186,
+  10.987035,
+   9.501727,
+   8.488104)
+
+tikz("norm_pca_m.tex",standAlone=TRUE,width=3,height=3);
+par(mar=my_mar, mgp=my_mgp,cex=my_cex) 
+plot(r[,1],r[,2],t="o",xlab="$m$",ylab="$\\|\\hat{\\bm K}-\\bm K\\|$",lwd=3) 
+lines(r[,1],fff,col="red",lwd=3)
+dev.off()
+
+system("pdflatex norm_pca_m.tex")
+
+r <- rbind(c( 100, 1.874253),
+           c( 250, 4.923609),
+           c( 500, 9.114774),
+           c(1000,18.907895),
+           c(1500,27.906247),
+           c(2000,39.030703),
+           c(2500,45.969439),
+           c(3000,58.737486),
+           c(4000,75.526354),
+           c(5000,95.037361))
+
+tikz("norm_pca_n.tex",standAlone=TRUE,width=3,height=3);
+par(mar=my_mar, mgp=my_mgp,cex=my_cex) 
+plot(r[,1],r[,2],t="o",xlab="$n$",ylab="$\\|\\hat{\\bm K}-\\bm K\\|$",lwd=3)
+dev.off()
+
+system("pdflatex norm_pca_n.tex")