Try to minimize matrix inversions in calculating quadratic forms.

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)

R/ergm.MCMLE.R

@@ -321,9 +321,7 @@ ergm.MCMLE <- function(init, s, s.obs,
       # 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_qrssolve(estdiff, Vm)
       if(d2<2) last.adequate <- TRUE
     }
 
@@ -411,7 +409,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_qrssolve(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 +446,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_qrssolve(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 +589,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=1e-07){
   y <- c(y)
-  if(xTAx(y,U)>=1) stop("Point is not in the interior of the ellipsoid.")
+  if(xTAx_qrssolve(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_qrssolve(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_qrssolve(y-x, W, tol=tol)
 }

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)
     }