tmp-tests/test-downdating2.R

In privefl/bigutilsr: Utility Functions for Large-scale Data

set.seed(1)
n <- 2e3
p <- 20e3
U <- matrix(0, n, 10); U[] <- rnorm(length(U))
V <- matrix(0, p, 10); V[] <- rnorm(length(V))
X <- tcrossprod(U, V) + 5 * matrix(rnorm(n * p), n, p)
system.time(svd2 <- PRIMME::svds(X, 10, isreal = TRUE))

ind <- -(1:1000)
A <- -X[, -ind]
UT_A <- crossprod(svd2$u, A)
M_A <- A - svd2$u %*% UT_A
svd_A <- svd(M_A)
U_A <- svd_A$u
# eig_A <- eigen(crossprod(A) - crossprod(UT_A), symmetric = TRUE)
# VD_A <- sweep(eig_A$vectors, 2, sqrt(eig_A$values), '*')
# VDI_A <- sweep(VD_A, 2, eig_A$values, '/')
# U_A <- A %*% VD_A - svd2$u %*% UT_A %*% VDI_A

BT_V <- svd2$v[-ind, ]
# B <- matrix(0, ncol(X), length(ind)); B[cbind(-ind, seq_along(ind))] <- 1
# all.equal(BT_V, crossprod(B, svd2$v))
# all.equal(crossprod(B), diag(100))
# M_B <- -tcrossprod(svd2$v, BT_V)
# M_B[cbind(-ind, seq_along(ind))] <- M_B[cbind(-ind, seq_along(ind))] + 1
# svd_B <- svd(M_B)
eig_B <- eigen(diag(nrow(BT_V)) - tcrossprod(BT_V), symmetric = TRUE)
VD_B <- sweep(eig_B$vectors, 2, sqrt(eig_B$values), '*')
VDI_B <- sweep(VD_B, 2, eig_B$values, '/')
# U_B <- -svd2$v %*% crossprod(BT_V, VDI_B)
# U_B[-ind, ] <- U_B[-ind, ] + VDI_B


Q1 <- UT_A %*% BT_V + diag(svd2$d)
VD_A <- sweep(svd_A$v, 2, svd_A$d, '*')
Q2 <- crossprod(VD_A, BT_V)
# VD_B <- sweep(svd_B$v, 2, svd_B$d, '*')
Q3 <- UT_A %*% VD_B
Q4 <- crossprod(VD_A, VD_B)
Q <- cbind(rbind(Q1, Q2), rbind(Q3, Q4))
svd_Q <- svd(Q, nu = 10, nv = 10)
U2 <- cbind(svd2$u, U_A) %*% svd_Q$u
BT_V2 <- svd2$v[ind, ]
V2 <- cbind(BT_V2, BT_V2 %*% crossprod(-BT_V, VDI_B)) %*% svd_Q$v

# MTM <- crossprod(A) - crossprod(crossprod(svd2$u, A))

X2 <- X[, ind]
u0 <- svd2$u
system.time(
  svd2.1 <- PRIMME::svds(X2, 10, isreal = TRUE))
system.time(
  svd2.2 <- PRIMME::svds(X2, 10, isreal = TRUE, u0 = u0, v0 = V[ind, ]))
system.time(
  svd2.3 <- PRIMME::svds(X2, 10, isreal = TRUE, u0 = u0, v0 = svd(V[ind, ])$u))
system.time(
  svd2.4 <- PRIMME::svds(X2, 10, isreal = TRUE, u0 = u0,
                         v0 = crossprod(X2, sweep(u0, 2, svd2$d, '/'))))
system.time(
  svd2.5 <- PRIMME::svds(X2, 10, isreal = TRUE, u0 = NULL, v0 = V2))

plot(U2, svd2.1$u); abline(0, 1, col = "red"); abline(0, -1, col = "red")
plot(V2, svd2.1$v); abline(0, 1, col = "red"); abline(0, -1, col = "red")
all.equal(svd_Q$d[1:10], svd2.1$d)