R/eff.ini.maxeig.general.R

In mxjki/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))
}