From d6f579b3cd9b242c608ed19b055d023662f7edc7 Mon Sep 17 00:00:00 2001
From: "Pavel N. Krivitsky" <p.krivitsky@unsw.edu.au>
Date: Mon, 18 Dec 2023 18:22:03 +1100
Subject: [PATCH] Implemented ginv_eigen(), an alternative to MASS:ginv() that
 uses eigendecomposition instead of QR decompsition and its scaled and xTAx
 variants. Bumped minor version to 4.10.

---
 DESCRIPTION      |  4 +--
 NAMESPACE        |  3 ++
 R/matrix.utils.R | 76 ++++++++++++++++++++++++++++++++++++++++--------
 man/ssolve.Rd    | 29 +++++++++++++-----
 man/xTAx.Rd      | 16 ++++++++++
 5 files changed, 106 insertions(+), 22 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 8b9da83..8c5c2c4 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: statnet.common
-Version: 4.9.0-436
-Date: 2023-06-27
+Version: 4.10.0-438
+Date: 2023-12-18
 Title: Common R Scripts and Utilities Used by the Statnet Project Software
 Authors@R: c(
     person(c("Pavel", "N."), "Krivitsky", role=c("aut","cre"), email="pavel@statnet.org", comment=c(ORCID="0000-0002-9101-3362", affiliation="University of New South Wales")),
diff --git a/NAMESPACE b/NAMESPACE
index ce89c57..ffb7833 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -86,6 +86,7 @@ export(eval_lhs.formula)
 export(filter_rhs.formula)
 export(fixed.pval)
 export(forkTimeout)
+export(ginv_eigen)
 export(handle.controls)
 export(is.SPD)
 export(is.linwmatrix)
@@ -140,8 +141,10 @@ export(update_snctrl)
 export(vector.namesmatch)
 export(xAxT)
 export(xTAx)
+export(xTAx_eigen)
 export(xTAx_qrsolve)
 export(xTAx_qrssolve)
+export(xTAx_seigen)
 export(xTAx_solve)
 export(xTAx_ssolve)
 importFrom(coda,as.mcmc)
diff --git a/R/matrix.utils.R b/R/matrix.utils.R
index f02ae27..8afbfcf 100644
--- a/R/matrix.utils.R
+++ b/R/matrix.utils.R
@@ -85,6 +85,27 @@ sandwich_solve <- function(A, B, ...) {
   solve(A, t(solve(A, B, ...)), ...)
 }
 
+#' @describeIn xTAx Evaluate \eqn{x' A^{-1} x} for vector \eqn{x} and
+#'   matrix \eqn{A} (symmetric, nonnegative-definite) via
+#'   eigendecomposition; returns `rank` and `nullity` as attributes
+#'   just in case subsequent calculations (e.g., hypothesis test
+#'   degrees of freedom) are affected.
+#'
+#'   Decompose \eqn{A = P L P'} for \eqn{L} diagonal matrix of
+#'   eigenvalues and \eqn{P} orthogonal. Then \eqn{A^{-1} = P L^{-1}
+#'   P'}.
+#'
+#'   Substituting, \deqn{x' A^{-1} x = x' P L^{-1} P' x
+#'   = h' L^{-1} h} for \eqn{h = P' x}.
+#'
+#' @export
+xTAx_eigen <- function(x, A, tol=sqrt(.Machine$double.eps), ...) {
+  e <- eigen(A, symmetric=TRUE)
+  keep <- e$values > max(tol * e$values[1L], 0)
+  h <- drop(crossprod(x, e$vectors[, keep, drop=FALSE]))
+  structure(sum(h*h/e$values[keep]), rank = sum(keep), nullity = sum(!keep))
+}
+
 
 #' Wrappers around matrix algebra functions that pre-*s*cale their
 #' arguments
@@ -96,16 +117,23 @@ sandwich_solve <- function(A, B, ...) {
 #' functions first scale the matrix's rows and/or columns by its
 #' diagonal elements and then undo the scaling on the result.
 #'
-#' `ssolve()`, `sginv()`, and `snearPD()` wrap [solve()],
-#' [MASS::ginv()], and [Matrix::nearPD()], respectively. `srcond()`
-#' returns the reciprocal condition number of [rcond()] net of the
-#' above scaling. `xTAx_ssolve`, `xTAx_qrssolve`, and
-#' `sandwich_ssolve` wrap the corresponding \pkg{statnet.common}
-#' functions.
+#' `ginv_eigen()` reimplements [MASS::ginv()] but using
+#' eigendecomposition rather than SVD; this means that it is only
+#' suitable for symmetric matrices, but that detection of negative
+#' eigenvalues is more robust.
+#' 
+#' `ssolve()`, `sginv()`, `sginv_eigen()`, and `snearPD()` wrap
+#' [solve()], [MASS::ginv()], `ginv_eigen()`, and [Matrix::nearPD()],
+#' respectively. `srcond()` returns the reciprocal condition number of
+#' [rcond()] net of the above scaling. `xTAx_ssolve()`,
+#' `xTAx_qrssolve()`, `xTAx_seigen()`, and `sandwich_ssolve()` wrap
+#' the corresponding \pkg{statnet.common} functions.
 #'
 #' @param snnd assume that the matrix is symmetric non-negative
-#'   definite (SNND). If it's "obvious" that it's not (e.g., negative
-#'   diagonal elements), an error is raised.
+#'   definite (SNND). This typically entails scaling that converts
+#'   covariance to correlation and use of eigendecomposition rather
+#'   than singular-value decomposition. If it's "obvious" that the matrix is
+#'   not SSND (e.g., negative diagonal elements), an error is raised.
 #'
 #' @param x,a,b,X,A,B,tol,... corresponding arguments of the wrapped functions.
 #'
@@ -135,7 +163,7 @@ ssolve <- function(a, b, ..., snnd = TRUE) {
 #' @rdname ssolve
 #'
 #' @export
-sginv <- function(X, ..., snnd = TRUE) {
+sginv <- function(X, ..., snnd = TRUE){
   d <- diag(as.matrix(X))
   d <- ifelse(d==0, 1, 1/d)
 
@@ -143,14 +171,38 @@ sginv <- function(X, ..., snnd = TRUE) {
     d <- sqrt(d)
     if(anyNA(d)) stop("Matrix a has negative elements on the diagonal, and snnd=TRUE.")
     dd <- rep(d, each = length(d)) * d
-    X <- X * dd
-    MASS::ginv(X, ...) * dd
+    ginv_eigen(X * dd, ...) * dd
   } else {
     dd <- rep(d, each = length(d))
-    MASS::ginv(X*d, ...) * dd
+    MASS::ginv(X * dd, ...) * t(dd)
   }
 }
 
+#' @rdname ssolve
+#' @export
+ginv_eigen <- function(X, tol = sqrt(.Machine$double.eps), ...){
+  e <- eigen(X, symmetric=TRUE)
+  keep <- e$values > max(tol * e$values[1L], 0)
+  tcrossprod(e$vectors[, keep, drop=FALSE] / rep(e$values[keep],each=ncol(X)), e$vectors[, keep, drop=FALSE])
+}
+
+
+#' @rdname ssolve
+#'
+#' @export
+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
+  x <- x * d
+  
+  xTAx_eigen(x, A, tol=tol, ...)
+}
+
 #' @rdname ssolve
 #'
 #' @export
diff --git a/man/ssolve.Rd b/man/ssolve.Rd
index 2bb478b..b14142f 100644
--- a/man/ssolve.Rd
+++ b/man/ssolve.Rd
@@ -3,6 +3,8 @@
 \name{ssolve}
 \alias{ssolve}
 \alias{sginv}
+\alias{ginv_eigen}
+\alias{xTAx_seigen}
 \alias{srcond}
 \alias{snearPD}
 \alias{xTAx_ssolve}
@@ -15,6 +17,10 @@ ssolve(a, b, ..., snnd = TRUE)
 
 sginv(X, ..., snnd = TRUE)
 
+ginv_eigen(X, tol = sqrt(.Machine$double.eps), ...)
+
+xTAx_seigen(x, A, tol = sqrt(.Machine$double.eps), ...)
+
 srcond(x, ..., snnd = TRUE)
 
 snearPD(x, ...)
@@ -27,8 +33,10 @@ sandwich_ssolve(A, B, ...)
 }
 \arguments{
 \item{snnd}{assume that the matrix is symmetric non-negative
-definite (SNND). If it's "obvious" that it's not (e.g., negative
-diagonal elements), an error is raised.}
+definite (SNND). This typically entails scaling that converts
+covariance to correlation and use of eigendecomposition rather
+than singular-value decomposition. If it's "obvious" that the matrix is
+not SSND (e.g., negative diagonal elements), an error is raised.}
 
 \item{x, a, b, X, A, B, tol, ...}{corresponding arguments of the wrapped functions.}
 }
@@ -41,12 +49,17 @@ functions first scale the matrix's rows and/or columns by its
 diagonal elements and then undo the scaling on the result.
 }
 \details{
-\code{ssolve()}, \code{sginv()}, and \code{snearPD()} wrap \code{\link[=solve]{solve()}},
-\code{\link[MASS:ginv]{MASS::ginv()}}, and \code{\link[Matrix:nearPD]{Matrix::nearPD()}}, respectively. \code{srcond()}
-returns the reciprocal condition number of \code{\link[=rcond]{rcond()}} net of the
-above scaling. \code{xTAx_ssolve}, \code{xTAx_qrssolve}, and
-\code{sandwich_ssolve} wrap the corresponding \pkg{statnet.common}
-functions.
+\code{ginv_eigen()} reimplements \code{\link[MASS:ginv]{MASS::ginv()}} but using
+eigendecomposition rather than SVD; this means that it is only
+suitable for symmetric matrices, but that detection of negative
+eigenvalues is more robust.
+
+\code{ssolve()}, \code{sginv()}, \code{sginv_eigen()}, and \code{snearPD()} wrap
+\code{\link[=solve]{solve()}}, \code{\link[MASS:ginv]{MASS::ginv()}}, \code{ginv_eigen()}, and \code{\link[Matrix:nearPD]{Matrix::nearPD()}},
+respectively. \code{srcond()} returns the reciprocal condition number of
+\code{\link[=rcond]{rcond()}} net of the above scaling. \code{xTAx_ssolve()},
+\code{xTAx_qrssolve()}, \code{xTAx_seigen()}, and \code{sandwich_ssolve()} wrap
+the corresponding \pkg{statnet.common} functions.
 }
 \examples{
 x <- rnorm(2, sd=c(1,1e12))
diff --git a/man/xTAx.Rd b/man/xTAx.Rd
index 4bb067e..ceffe12 100644
--- a/man/xTAx.Rd
+++ b/man/xTAx.Rd
@@ -6,6 +6,7 @@
 \alias{xTAx_solve}
 \alias{xTAx_qrsolve}
 \alias{sandwich_solve}
+\alias{xTAx_eigen}
 \title{Common quadratic forms}
 \usage{
 xTAx(x, A)
@@ -17,6 +18,8 @@ xTAx_solve(x, A, ...)
 xTAx_qrsolve(x, A, tol = 1e-07, ...)
 
 sandwich_solve(A, B, ...)
+
+xTAx_eigen(x, A, tol = sqrt(.Machine$double.eps), ...)
 }
 \arguments{
 \item{x}{a vector}
@@ -56,4 +59,17 @@ and \code{nullity} as attributes just in case subsequent calculations
 \item \code{sandwich_solve()}: Evaluate \eqn{A^{-1}B(A')^{-1}} for \eqn{B} a
 square matrix and \eqn{A} invertible.
 
+\item \code{xTAx_eigen()}: Evaluate \eqn{x' A^{-1} x} for vector \eqn{x} and
+matrix \eqn{A} (symmetric, nonnegative-definite) via
+eigendecomposition; returns \code{rank} and \code{nullity} as attributes
+just in case subsequent calculations (e.g., hypothesis test
+degrees of freedom) are affected.
+
+Decompose \eqn{A = P L P'} for \eqn{L} diagonal matrix of
+eigenvalues and \eqn{P} orthogonal. Then \eqn{A^{-1} = P L^{-1}
+  P'}.
+
+Substituting, \deqn{x' A^{-1} x = x' P L^{-1} P' x
+  = h' L^{-1} h} for \eqn{h = P' x}.
+
 }}