# R/eff.ini.maxeig.general.R In EfficientMaxEigenpair: Efficient Initials for Computing the Maximal Eigenpair

#' @title General matrix maximal eigenpair
#'
#' @description Calculate the maximal eigenpair for the general matrix.
#'
#' @param A The input general matrix.
#' @param v0_tilde The unnormalized initial vector \eqn{\tilde{v0}}.
#' @param z0 The type of initial \eqn{z_0} used to calculate the approximation of \eqn{\rho(Q)}.
#'        There are three types: 'fixed', 'Auto' and 'numeric' corresponding to three choices
#'        of \eqn{z_0} in paper.
#' @param z0numeric The numerical value assigned to initial \eqn{z_0} as an approximation of
#'        \eqn{\rho(Q)} when z_0='numeric'.
#' @param xi The coefficient used to form the convex combination of \eqn{\delta_1^{-1}} and
#'        \eqn{(v_0,-Q*v_0)_\mu}, it should between 0 and 1.
#' @param digit.thresh The precise level of output results.
#'
#' @return A list of eigenpair object are returned, with components \eqn{z}, \eqn{v} and \eqn{iter}.
#' \item{z}{The approximating sequence of the maximal eigenvalue.}
#' \item{v}{The approximating eigenfunction of the corresponding eigenvector.}
#' \item{iter}{The number of iterations.}
#'
#' @seealso \code{\link{eff.ini.maxeig.tri}} for the tridiagonal matrix maximal
#' eigenpair by rayleigh quotient iteration algorithm.
#' \code{\link{eff.ini.maxeig.shift.inv.tri}} for the tridiagonal matrix
#' maximal eigenpair by shifted inverse iteration algorithm.

#' @examples A = matrix(c(1, 1, 3, 2, 2, 2, 3, 1, 1), 3, 3)
#' eff.ini.maxeig.general(A, v0_tilde = rep(1, dim(A)[1]), z0 = 'fixed')
#'
#' A = matrix(c(1, 1, 3, 2, 2, 2, 3, 1, 1), 3, 3)
#' eff.ini.maxeig.general(A, v0_tilde = rep(1, dim(A)[1]), z0 = 'Auto')
#'
#' ##Symmetrizing A converge to second largest eigenvalue
#' A = matrix(c(1, 3, 9, 5, 2, 14, 10, 6, 0, 11, 11, 7, 0, 0, 1, 8), 4, 4)
#' S = (t(A) + A)/2
#' N = dim(S)[1]
#' a = diag(S[-1, -N])
#' b = diag(S[-N, -1])
#' c = rep(NA, N)
#' c[1] = -diag(S)[1] - b[1]
#' c[2:(N - 1)] = -diag(S)[2:(N - 1)] - b[2:(N - 1)] - a[1:(N - 2)]
#' c[N] = -diag(S)[N] - a[N - 1]
#'
#' z0ini = eff.ini.maxeig.tri(a, b, c, xi = 7/8)$z[1] #' eff.ini.maxeig.general(A, v0_tilde = rep(1, dim(A)[1]), z0 = 'numeric', #' z0numeric = 28 - z0ini) #' @export eff.ini.maxeig.general = function(A, v0_tilde = NULL, z0 = NULL, z0numeric, xi = 1, digit.thresh = 6) { if (xi < 0 | xi > 1) stop("The coefficient xi should between 0 and 1!") N = dim(A)[1] h = rep(NA, N) m = 0 # check input matrix Q if (max(rowSums(A)) <= 1e-10) { Q = A } else { # print('Adjust the diagonal of matrix by Q=A-mI') m = max(rowSums(A)) print(paste0("m=", max(rowSums(A)))) Q = A - m * diag(1, N) } if ((sum(rowSums(Q)[1:(N - 1)] < 0) == 0) & (rowSums(Q)[N] < 0)) { h = rep(1, N) } else { h[1] = 1 h[2:N] = solve(Q[-N, -1], -Q[1:(N - 1), 1]) } Q=diag(1/h)%*%Q%*%diag(h) q = -diag(Q) P = diag(1/q) %*% Q + diag(rep(1, N)) if (-Q[N, N] - sum(Q[N, 1:(N - 1)] * h[1:(N - 1)])/h[N] < sum(Q[N, 1:(N - 1)] * h[1:(N - 1)])/h[N]/100) { x = rep(1, N) } else { x = rep(NA, N) x[1] = 1 x[2:N] = solve(diag(1, N - 1) - P[-1, -1], P[2:N, 1]) } if (is.null(v0_tilde)) { v0_tilde = sqrt(x) } if (is.null(z0)) { Q1= t(Q) mu = rep(1, N) mu[2:N] = solve(Q1[-N, -1], -Q1[1:(N - 1), 1]) deltapart1 = mu * sqrt(x) deltapart1sum = cumsum(deltapart1) deltapart2 = mu * x^(3/2) deltapart2sum = rev(cumsum(rev(deltapart2))) delta1 = max(c(sqrt(x[1:(N - 1)]) * deltapart1sum[1:(N - 1)] + 1/sqrt(x[1:(N - 1)]) * deltapart2sum[2:N], sqrt(x[N]) * deltapart1sum[N]))/(1 - x[2]) v0 = v0_tilde/sqrt(sum(v0_tilde^2 * mu)) zstart = xi * 1/delta1 + (1 - xi) * sum(v0 * (-Q %*% v0) * mu) ray = ray.quot.general(A = Q, mu = mu, v0_tilde = v0_tilde, zstart = zstart, digit.thresh = digit.thresh) } else { v0 = v0_tilde/sqrt(sum(v0_tilde^2)) zstart = switch(z0, fixed = 0, Auto = sum(v0 * (-Q %*% v0)), numeric = z0numeric) ray = ray.quot.general(A = Q, mu = rep(1, N), v0_tilde = v0_tilde, zstart = zstart, digit.thresh = digit.thresh) } return(list(z = m - unlist(ray$z), v = ray$v, iter = ray$iter))
}


## Try the EfficientMaxEigenpair package in your browser

Any scripts or data that you put into this service are public.

EfficientMaxEigenpair documentation built on May 2, 2019, 2:17 a.m.