The optimal $u$ and $v$ computed with `mbound()` are plotted in figure 2.

The `wPCA()` function has a `bnd` parameter, which can be equal to `all`, or `col`, or `row`, or `opt`. If `bnd="all"` then $c_{ij}=\max w_{ij}$, as in @kiers_97, if `bnd="row"` then $c_{ij}=\max_j w_{ij}$, as in @groenen_giaquinto_kiers_03. If `bnd="col"` then $c_{ij}=\max_i w_{ij}$ and if `bnd="opt"` then we compute the optimal $u_iv_j$. First look at $p=1$. For all four values of the `bnd` parameter we find a chi-square equal to 709.9526292976 with $168-((24+7)-1)=138$ degrees of freedom. The iteration count for `bnd="all","col","row","opt"` is, respectively, 208, 151, 21, 17. Because in this example the between-row/within-column variation is so much larger than the within-row/between-column variation, it makes sense that `bnd="row"` performs well, almost as bad as `bnd="opt"`, while `bnd="col"` is bad, almost as bad as `bnd="all"`. The one-dimensional plots of $A$ and $B$ from $z=ab'$ are in figure 3.

We also computed the spectral norms of the derivative at the solution, using the functions in `lowrank.R`. For `bnd="all","col","row","opt"` the convergence rate is, respectively, 0.9710924907, 0.955475149, 0.660381091, 0.6152936489. For $p=2$ for all four values of the `bnd` parameter we find chi-square equal to 215.349822881 with $168-((24+7)*2-4)=110$ degrees of freedom. The iteration count for `bnd="all","col","row","opt"` is, respectively, 164, 99, 46, 35.

For `bnd="all","col","row","opt"` the convergence rate is, respectively, 0.9715807406, 0.961724624, 0.9042846128, 0.8856193743. #Discussion The main question, which we did not answer, is if the decrease in computer time by lowering the number of iterations compensates for the extra computation of the optimal rank-one bounds. Computing the optimal bounds is a quadratic programming problem of size $nm$, which can be substantial and rapidly becomes impractical for really large data sets. In that case, alternative strategies are needed, for example starting with row or column bounds and performing just a few iterations of an iterative quadratic programming method. Of course we also have not exploited $\ell_1$ and $\ell_\infty$ solutions, which have linear programming solutions that are less expensive than the quadratic ones. Also note that in each majorization step we have solved a low-rank approximation problem. We would also have a convergent algorithm if we replaced the singular value decomposition in each majorization step by one or more alternating least squares iterations. This will generally lead to more majorization iterations, but each majorization iteration becomes less expensive. In addition we can use $A$ and $B$ from $X\approx AB'$ for the previous majorization iteration as starting points for the next one. There is an interesting recent paper by @szlam_tulloch_tygert_17, who show that just a few alternating least squares iterations will generally provide a very good low-rank approximation. This implies that it may not pay off to optimize the number of iterations needed for convergence, if convergence is not really of much practical value anyway. It also indicates that performing only a small number of "inner" alternating least squares iterations in each majorization step is probably a good strategy. Again, that is something to be explored. #Appendix: Code ##mbound.R ```r mbound <- function (w) { n <- nrow (w) m <- ncol (w) ww <- log (w) g1 <- matrix(0, n * m, n) g2 <- matrix(0, n * m, m) k <- 1 for (j in 1:m) { for (i in 1:n) { g1[k, i] <- 1 g2[k, j] <- 1 k <- k + 1 } } g <- cbind (g1, g2) f <- log (as.vector (w)) h <- lsi (g, f, g, f) u <- exp (h[1:n]) v <- exp (h[n + (1:m)]) return (list (u = u, v = v)) } ``` ##wPCA.R ```r gfunc <- function (z, x, u, v, w) { return ((z + (w / outer (u, v)) * (x - z)) * sqrt(outer(u,v))) } sfunc <- function (z, u, v) { return (z / sqrt (outer (u, v))) } lfunc <- function (x, p) { l <- as.matrix(eigen (crossprod(x))$vectors[,1:p]) return (x %*% tcrossprod (l)) } wPCA <- function (x, w, p = 2, bnd = "opt", itmax = 1000, eps = 1e-6, verbose = FALSE) { n <- nrow (x) m <- ncol (x) if (bnd == "all") { u <- rep(1, n) v <- rep (max (w), m) } if (bnd == "row") { u <- apply (w, 1, max) v <- rep (1, m) } if (bnd == "col") { u <- rep (1, n) v <- apply (w, 2, max) } if (bnd == "opt") { uv <- mbound (w) u <- uv$u v <- uv$v } zold <- lfunc (x, p) fold <- sum (w * (x - zold) ^ 2) itel <- 1 repeat { z1 <- gfunc (zold, x, u, v, w) z2 <- lfunc (z1, p) znew <- sfunc (z2, u, v) fnew <- sum (w * (x - znew) ^ 2) if (verbose) { cat ( "iteration ", formatC ( itel, digits = 0, width = 3, format = "d" ), " fold", formatC ( fold, digits = 6, width = 10, format = "f" ), " fnew", formatC ( fnew, digits = 6, width = 10, format = "f" ), "\n" ) } if (((fold - fnew) < eps) || (itel == itmax)) break itel <- itel + 1 zold <- znew fold <- fnew } return (list (z = znew, f = fnew, itel = itel, u = u, v = v)) } ``` ##lowrank.R ```r gfunc <- function (z, x, u, v, w) { return ((z + (w / outer (u, v)) * (x - z)) * sqrt(outer(u,v))) } sfunc <- function (z, u, v) { return (z / sqrt (outer (u, v))) } lfunc <- function (x, p) { l <- as.matrix(eigen (crossprod(x))$vectors[,1:p]) yy <- x %*% tcrossprod (l) return (yy) } ffunc <- function (z, x, p, u, v, w) { z1 <- gfunc (z, x, u, v, w) z2 <- lfunc (z1, p) z3 <- sfunc (z2, u, v) return (as.vector(z3)) } cfunc <- function (z, x, p, u, v, w) { n <- length (u) m <- length (v) z <- matrix (z, n, m) return (ffunc (z, x, p, u, v, w)) } hfunc <- function (delta, u, v, w) { ((1 - w / outer (u, v)) * delta) * sqrt (outer (u, v)) } lrPerturb <- function(z, delta, x, p, u, v, w) { m <- ncol (x) gz <- gfunc (z, x, u, v, w) hd <- hfunc (delta, u, v, w) gh <- crossprod(gz, hd) + crossprod(hd, gz) hp <- matrix(0, m, m) e <- eigen(crossprod(gz)) lt <- e$vectors vt <- e$values lp <- as.matrix(lt[, 1:p]) for (s in 1:p) { vi <- ifelse (vt == vt[s], 0, 1 / (vt - vt[s])) cinv <- lt %*% diag (vi) %*% t(lt) hp <- hp + cinv %*% gh %*% outer (lp[, s], lp[, s]) } dz <- hd %*% tcrossprod(lp) - gz %*% (hp + t(hp)) return (sfunc (dz, u, v)) } lrDerive <- function(z, x, p, u, v, w) { n <- nrow(z) m <- ncol(z) k <- 1 dz <- matrix(0, n * m, n * m) for (j in 1:m) for (i in 1:n) { dx <- matrix(0, n, m) dx[i, j] <- 1 dz[, k] <- as.vector(lrPerturb(z, dx, x, p, u, v, w)) k <- k + 1 } return(dz) } ``` #References