ssvd.R

n = 500;
m = 700;
d = 10;
X = matrix(data=1*(runif(n*m) < 0.1), nr=n, nc=m);

S = matrix(data=rnorm(m*(d+2)), nr=m, nc=(d+2));

XS = X %*% S
XX = X %*% t(X);
XX2 = XX %*% t(XX);

XXXS = XX %*% XS;
XX2XS = XX2 %*% XS;

Y = cbind(XS, XXXS, XX2XS);

# 3(d+2) x 3(d+2) eigendecomp
YYe = eigen(t(Y) %*% Y);

Q = Y %*% (YYe$vectors %*% diag(1.0 / sqrt(YYe$values)));

B = t(Q) %*% X;

# 3(d+2) x 3(d+2) eigendecomp
EB = eigen(B %*% t(B))
E = diag(sqrt(EB$values[1:d]))
QU = t(t(EB$vectors) %*% t(Q))
U = QU[,1:d]
V = t(B) %*% diag(1.0 / sqrt(EB$values)) %*% EB$vectors;
#V = t(EB$vectors) %*% diag(1.0 / sqrt(EB$values)) %*% B;
V = V[,1:d]

print(sum((X - U %*% E %*% t(V))^2));

A = svd(X)
print(sum((X - A$u[,1:d] %*% diag(A$d[1:d]) %*% t(A$v[,1:d]))^2))
print(sum((X - A$u[,1:(d-1)] %*% diag(A$d[1:(d-1)]) %*% t(A$v[,1:(d-1)]))^2))