From 3628ac1d4be44726b31adf601c13bfcf04952de1 Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Mon, 11 Dec 2023 20:45:39 +1100 Subject: [PATCH] Try to minimize matrix inversions in calculating quadratic forms, including using eigendecomposition to compute most xTAx. --- R/approx.hotelling.diff.test.R | 4 ++-- R/ergm.MCMLE.R | 32 ++++++++++++-------------------- R/ergm.estimate.R | 2 +- R/ergm.utility.R | 22 ++++++++++++++++++++++ 4 files changed, 37 insertions(+), 23 deletions(-) diff --git a/R/approx.hotelling.diff.test.R b/R/approx.hotelling.diff.test.R index 0c215846e..ce10d4de1 100644 --- a/R/approx.hotelling.diff.test.R +++ b/R/approx.hotelling.diff.test.R @@ -166,7 +166,7 @@ approx.hotelling.diff.test<-function(x,y=NULL, mu0=0, assume.indep=FALSE, var.eq if(p==0) stop("data are essentially constant") - ivcov.d <-sginv(vcov.d[!novar,!novar,drop=FALSE]) + ivcov.d <- sginv(vcov.d[!novar,!novar,drop=FALSE]) method <- paste0("Hotelling's ", NVL2(y, "Two", "One"), @@ -189,7 +189,7 @@ approx.hotelling.diff.test<-function(x,y=NULL, mu0=0, assume.indep=FALSE, var.eq " do not vary in x or in y but have differences equal to mu0"), "; they have been ignored for the purposes of testing.") } - T2 <- xTAx((d-mu0)[!novar], ivcov.d) + T2 <- xTAx_qrssolve((d-mu0)[!novar], vcov.d[!novar,!novar,drop=FALSE]) } NANVL <- function(z, ifNAN) ifelse(is.nan(z),ifNAN,z) diff --git a/R/ergm.MCMLE.R b/R/ergm.MCMLE.R index 32dfc3ed5..a205f1094 100644 --- a/R/ergm.MCMLE.R +++ b/R/ergm.MCMLE.R @@ -315,15 +315,7 @@ ergm.MCMLE <- function(init, s, s.obs, estdiff <- NVL3(esteq.obs, colMeans(.), 0) - colMeans(esteq) pprec <- diag(sqrt(control$MCMLE.MCMC.precision), nrow=length(estdiff)) Vm <- pprec%*%(cov(esteq) - NVL3(esteq.obs, cov(.), 0))%*%pprec - - # Ensure tolerance hyperellipsoid is PSD. (If it's not PD, the - # ellipsoid is workable, if flat.) Special handling is required - # if some statistic has a variance of exactly 0. - novar <- diag(Vm) == 0 - Vm[!novar,!novar] <- snearPD(Vm[!novar,!novar,drop=FALSE], posd.tol=0, base.matrix=TRUE)$mat - iVm <- sginv(Vm, tol=.Machine$double.eps^(3/4)) - diag(Vm)[novar] <- sqrt(.Machine$double.xmax) # Virtually any nonzero difference in estimating functions will map to a very large number. - d2 <- xTAx(estdiff, iVm) + d2 <- xTAx_seigen(estdiff, Vm) if(d2<2) last.adequate <- TRUE } @@ -411,7 +403,7 @@ ergm.MCMLE <- function(init, s, s.obs, } }else if(control$MCMLE.termination=='confidence'){ if(!is.null(estdiff.prev)){ - d2.prev <- xTAx(estdiff.prev, iVm) + d2.prev <- xTAx_seigen(estdiff.prev, Vm) if(verbose) message("Distance from origin on tolerance region scale: ", format(d2), " (previously ", format(d2.prev), ").") d2.not.improved <- d2.not.improved[-1] if(d2 >= d2.prev){ @@ -448,9 +440,9 @@ ergm.MCMLE <- function(init, s, s.obs, estdiff <- estdiff[!hotel$novar] estcov <- estcov[!hotel$novar, !hotel$novar] - d2e <- xTAx(estdiff, iVm[!hotel$novar, !hotel$novar]) + d2e <- xTAx_seigen(estdiff, Vm[!hotel$novar, !hotel$novar]) if(d2e<1){ # Update ends within tolerance ellipsoid. - T2 <- try(.ellipsoid_mahalanobis(estdiff, estcov, iVm[!hotel$novar, !hotel$novar]), silent=TRUE) # Distance to the nearest point on the tolerance region boundary. + T2 <- try(.ellipsoid_mahalanobis(estdiff, estcov, Vm[!hotel$novar, !hotel$novar]), silent=TRUE) # Distance to the nearest point on the tolerance region boundary. if(inherits(T2, "try-error")){ # Within tolerance ellipsoid, but cannot be tested. message("Unable to test for convergence; increasing sample size.") .boost_samplesize(control$MCMLE.confidence.boost) @@ -591,26 +583,26 @@ ergm.MCMLE <- function(init, s, s.obs, } #' Find the shortest squared Mahalanobis distance (with covariance W) -#' from a point `y` to an ellipsoid defined by `x'U x = 1`, provided +#' from a point `y` to an ellipsoid defined by `x'U^-1 x = 1`, provided #' that `y` is in the interior of the ellipsoid. #' #' @param y a vector #' @param W,U a square matrix #' #' @noRd -.ellipsoid_mahalanobis <- function(y, W, U){ +.ellipsoid_mahalanobis <- function(y, W, U, tol=sqrt(.Machine$double.eps)){ y <- c(y) - if(xTAx(y,U)>=1) stop("Point is not in the interior of the ellipsoid.") + if(xTAx_seigen(y,U,tol=tol)>=1) stop("Point is not in the interior of the ellipsoid.") I <- diag(length(y)) - WU <- W%*%U - x <- function(l) c(solve(I+l*WU, y)) # Singluar for negative reciprocals of eigenvalues of WiU. - zerofn <- function(l) ERRVL(try({x <- x(l); xTAx(x,U)-1}, silent=TRUE), +Inf) + WUi <- t(qr.solve(qr(U, tol=tol), W)) + x <- function(l) c(solve(I+l*WUi, y)) # Singluar for negative reciprocals of eigenvalues of WiU. + zerofn <- function(l) ERRVL(try({x <- x(l); xTAx_seigen(x,U,tol=tol)-1}, silent=TRUE), +Inf) # For some reason, WU sometimes has 0i element in its eigenvalues. - eig <- Re(eigen(WU, only.values=TRUE)$values) + eig <- Re(eigen(WUi, only.values=TRUE)$values) lmin <- -1/max(eig) l <- uniroot(zerofn, lower=lmin, upper=0, tol=sqrt(.Machine$double.xmin))$root x <- x(l) - xTAx_qrssolve(y-x, W) + xTAx_seigen(y-x, W, tol=tol) } diff --git a/R/ergm.estimate.R b/R/ergm.estimate.R index 6bb13bf27..30eff0735 100644 --- a/R/ergm.estimate.R +++ b/R/ergm.estimate.R @@ -225,7 +225,7 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL, # or more statistics. if(inherits(Lout$par,"try-error")){ Lout$par <- try(eta0 - + sginv(-Lout$hessian, tol=.Machine$double.eps^(3/4)) %*% + + sginv(-Lout$hessian) %*% xobs, silent=TRUE) } diff --git a/R/ergm.utility.R b/R/ergm.utility.R index 698a24758..1627adfda 100644 --- a/R/ergm.utility.R +++ b/R/ergm.utility.R @@ -286,3 +286,25 @@ sginv <- function(X, tol = sqrt(.Machine$double.eps), ..., snnd = TRUE){ keep <- e$values > max(tol * e$values[1L], 0) e$vectors[, keep, drop=FALSE] %*% ((1/e$values[keep]) * t(e$vectors[, keep, drop=FALSE])) * dd } + +#' Compute \eqn{x'A^{*}x'} for a vector \eqn{x} and a matrix \eqn{A} +#' assuming that \eqn{A} is PSD and \eqn{A^*} is the Moore--Penrose +#' pseudoinverse, and after scaling \eqn{A} into a correlation matrix +#' (and scaling back the result). +#' +#' @noRd +xTAx_seigen <- function(x, A, tol=sqrt(.Machine$double.eps), ...) { + d <- diag(as.matrix(A)) + d <- ifelse(d<=0, 0, 1/d) + + d <- sqrt(d) + dd <- rep(d, each = length(d)) * d + A <- A * dd + + e <- eigen(A) + + keep <- e$values > max(tol * e$values[1L], 0) + + h <- drop(crossprod(x*d, e$vectors[,keep,drop=FALSE])) + sum(h*h/e$values[keep]) +}