In neurodata/lol: Linear Optimal Low-Rank Projection

require(lolR)
require(ggplot2)
require(MASS)
require(abind)
n=1000
d=30
r=3
#rot=-pi/4  # rotation in radians
#R = rbind(c(cos(rot), -sin(rot), 0, 0, 0), c(sin(rot), cos(rot), 0, 0, 0),
#          c(0,0,1,0,0), c(0,0,0,1,0), c(0,0,0,0,1))  # rotation matrix for rotation by 45 degrees in first 2 dimensions

Generate Data

mu0 = array(.25, dim=c(d))
Sigma0 = diag(c(array(1, dim=c(d/3)), array(.25, dim=c(2*d/3))))
Sigma1 = diag(c(array(.25, dim=c(d/3)), array(1, dim=c(d/3)), array(.25, dim=c(d/3))))
mu1 = array(.25, dim=c(d))
mu0[1:d/3] = 1
mu1[(d/3+1):(2*d/3)] = 1
mu0 = array(0, dim=c(d));  mu1=array(0, dim=c(d))
#Sigma1 = t(R) %*% Sigma0 %*% R

mus = abind(mu0, mu1, along=2)
Sigmas = abind(Sigma0, Sigma1, along=3)
data = lolR:::lol.sims.sim_gmm(mus, Sigmas, n=n, priors = c(0.5, 0.5))
X = data$X; Y = data$Y

Visualize Data

Below, we visualize the data rotated 45 degrees in the first and second dimensions:

X.dat = data.frame(x1=X[,10], x2=X[,11], class=Y)
X.dat$class <- factor(X.dat$class)
ggplot(X.dat, aes(x=x1, y=x2, color=class, group=class, fill=class)) +
  geom_point() + 
  stat_ellipse(type = "norm", level=.6825)

Projection vectors using LO

We look at the second 2 vectors the data projects onto with LOL (the top 2 non-mean vectors). These should be "between" the two ellipses:

lol.result <- lol.project.lol(X, Y, r=20, orthogonalize=FALSE)
lol.result$A <- lol.result$A

project.X <- lol.embed(X, lol.result$A)
proj.dat <- data.frame(x1=project.X[,2], x2=project.X[,3], class=factor(Y))
proj.dat$class <- factor(proj.dat$class)
ggplot(proj.dat, aes(x=x1, y=x2, color=class, group=class)) +
  geom_point() + 
  stat_ellipse(type = "norm", level=.6825)

As we can see, the directions are orthogonal to the point clouds due to the class-conditional covariance structure.

Projection vectors using QO

lol.result <- lol.project.lol(X, Y, r=20, second.moment='quadratic', orthogonalize=FALSE)
lol.result$A <- lol.result$A

project.X <- lol.embed(X, lol.result$A)
proj.dat <- data.frame(x1=project.X[,2], x2=project.X[,3], class=factor(Y))
proj.dat$class <- factor(proj.dat$class)
ggplot(proj.dat, aes(x=x1, y=x2, color=class, group=class)) +
  geom_point() + 
  stat_ellipse(type = "norm", level=.6825)

as we can see, the directions are exactly in-line with the point clouds.

lol.result <- lol.project.pls(X, Y, r=20)
lol.result$A <- lol.result$A

project.X <- lol.embed(X, lol.result$A)
proj.dat <- data.frame(x1=project.X[,1], x2=project.X[,2], class=factor(Y))
proj.dat$class <- factor(proj.dat$class)
ggplot(proj.dat, aes(x=x1, y=x2, color=class, group=class)) +
  geom_point() + 
  stat_ellipse(type = "norm", level=.6825)