@@ -264,59 +264,3 @@ trim_env_const_formula <- function(x, keep=NULL){
264264
265265 if (needs_env ) x else trim_env(x , keep )
266266}
267-
268- # ## A patched version of statnet.common::sginv() that uses
269- # ## eigendecomposition rather than the SVD for the case when symmetric
270- # ## non-negative-definite (snnd) is TRUE.
271- # ##
272- # ## TODO: Delete after incorporation into statnet.common.
273- sginv <- function (X , tol = sqrt(.Machine $ double.eps ), ... , snnd = TRUE ){
274- if (! snnd ) statnet.common :: sginv(X , ... , snnd )
275-
276- d <- diag(as.matrix(X ))
277- d <- ifelse(d == 0 , 1 , 1 / d )
278-
279- d <- sqrt(d )
280- if (anyNA(d )) stop(" Matrix a has negative elements on the diagonal, and snnd=TRUE." )
281- dd <- rep(d , each = length(d )) * d
282- X <- X * dd
283-
284- # # Perform inverse via eigendecomposition, removing too-small eigenvalues.
285- e <- eigen(X , symmetric = TRUE )
286- keep <- e $ values > max(tol * e $ values [1L ], 0 )
287- e $ vectors [, keep , drop = FALSE ] %*% ((1 / e $ values [keep ]) * t(e $ vectors [, keep , drop = FALSE ])) * dd
288- }
289-
290- # ' Evaluate a quadratic form with a possibly singular matrix using
291- # ' eigendecomposition after scaling to correlation.
292- # '
293- # ' Let \eqn{A} be the matrix of interest, and let \eqn{D} is a
294- # ' diagonal matrix whose diagonal is same as that of \eqn{A}.
295- # '
296- # ' Let \eqn{R = D^{-1/2} A D^{-1/2}}. Then \eqn{A = D^{1/2} R D^{1/2}} and
297- # ' \eqn{A^{-1} = D^{-1/2} R^{-1} D^{-1/2}}.
298- # '
299- # ' Decompose \eqn{R = P L P'} for \eqn{L} diagonal matrix of eigenvalues
300- # ' and \eqn{P} orthogonal. Then \eqn{R^{-1} = P L^{-1} P'}.
301- # '
302- # ' Substituting,
303- # ' \deqn{x' A^{-1} x = x' D^{-1/2} P L^{-1} P' D^{-1/2} x = h' L^{-1} h}
304- # ' for \eqn{h = P' D^{-1/2} x}.
305- # '
306- # ' @noRd
307- xTAx_seigen <- function (x , A , tol = sqrt(.Machine $ double.eps ), ... ) {
308- d <- diag(as.matrix(A ))
309- d <- ifelse(d < = 0 , 0 , 1 / d )
310-
311- d <- sqrt(d )
312- dd <- rep(d , each = length(d )) * d
313- A <- A * dd
314-
315- e <- eigen(A )
316-
317- keep <- e $ values > max(tol * e $ values [1L ], 0 )
318-
319- h <- drop(crossprod(x * d , e $ vectors [,keep ,drop = FALSE ]))
320-
321- structure(sum(h * h / e $ values [keep ]), rank = sum(keep ), nullity = sum(! keep ))
322- }
0 commit comments