R: Как добавить дрожание только на сингулярные матрицы внутри функции?

Question

У меня есть следующая функция, которую мне нужно (m) применить к списку из более чем 1500 больших матриц (Z) и списка векторов (p) той же длины. Тем не менее, я получаю ошибку, что некоторые матрицы являются единственными, как я уже опубликовал here. Вот моя функция:

kastner <- function(item, p) { print(item) 
           imp <- rowSums(Z[[item]]) 
           exp <- colSums(Z[[item]]) 
           x = p + imp 
           ac = p + imp - exp 
           einsdurchx = 1/as.vector(x)  
           einsdurchx[is.infinite(einsdurchx)] <- 0 
           A = Z[[item]] %*% diag(einsdurchx) 
           R = solve(diag(length(p))-A) %*% diag(p)  
           C = ac * einsdurchx 
           R_bar = diag(as.vector(C)) %*% R 
           rR_bar = round(R_bar) 
           return(rR_bar) 
} 

и моя команда mapply, которая также печатает имена бегущего матрицы:

KASTNER <- mapply(kastner, names(Z), p, SIMPLIFY = FALSE) 

Для того, чтобы преодолеть проблему сингулярности, я хочу добавить небольшое количество jitter в сингулярных матриц. Проблема начинается в строке 9 функции R = solve(diag(length(p))-A) %*% diag(p), так как этот термин (diag(length(p))-A) становится сингулярным и не может быть solve d. Я попытался добавить дрожание ко всем Z-матрицам в первой строке функции, используя: Z <- lapply(Z,function(x) jitter(x, factor = 0.0001, amount = NULL)), но это очень мало и создает все еще ошибки.

Поэтому моя идея состоит в том, чтобы проверить с if/else или что-то подобное if эта матрица diag(length(p))-A сингулярен (возможно с использованием собственных векторов для проверки коллинеарности) и добавить на этих матриц джиттера, else (если нет) команда solve должна выполняться, как это. Идеи, как реализовать это в функции? Благодаря

Вот некоторые примеры данных, хотя нет никаких проблем с особенностью, как я был не в состоянии восстановить эту ошибку на линии 9:

Z <- list("111.2012"= matrix(c(0,0,100,200,0,0,0,0,50,350,0,50,50,200,200,0), 
           nrow = 4, ncol = 4, byrow = T), 
      "112.2012"= matrix(c(10,90,0,30,10,90,0,10,200,50,10,350,150,100,200,10), 
            nrow = 4, ncol = 4, byrow = T)) 
p <- list("111.2012"=c(200, 1000, 100, 10), "112.2012"=c(300, 900, 50, 100)) 

Edit: небольшое количество о дрожании не должно быть проблематичным мои данные, поскольку у меня есть, вероятно, более 80% нулей в моих матрицах и большие значения. И меня интересуют только эти большие значения, но большое количество 0s, вероятно, является причиной сингулярности, но необходимо.

Answer 1

Поскольку вы не представили рабочий пример, я не смог проверить это легко, поэтому бремя доказывания на вас. :) В любом случае, это должна быть отправная точка для дальнейшего возиться. Комментарии в коде.

kastner <- function(item, p) { print(item) 
    imp <- rowSums(Z[[item]]) 
    exp <- colSums(Z[[item]]) 
    x = p + imp 
    ac = p + imp - exp 
    einsdurchx = 1/as.vector(x)  
    einsdurchx[is.infinite(einsdurchx)] <- 0 

    # start a chunk that repeats until you get a valid result 
    do.jitter <- TRUE # bureaucracy 
    while (do.jitter == TRUE) { 
    # run the code as usual 
    A = Z[[item]] %*% diag(einsdurchx) 
    # catch any possible errors, you can even catch "singularity" error here by 
    # specifying error = function(e) e 
    R <- tryCatch(solve(diag(length(p))-A) %*% diag(p), error = function(e) "jitterme") 
    # if you were able to solve(), and the result is a matrix (carefuly if it's a vector!)... 
    if (is.matrix(R)) { 
     # ... turn the while loop off 
     do.jitter <- FALSE 
    } else { 
     #... else apply some jitter and repeat by construcing A from a jittered Z[[item]] 
     Z[[item]] <- jitter(Z[[item]]) 
    } 
    } 
    C = ac * einsdurchx 
    R_bar = diag(as.vector(C)) %*% R 
    rR_bar = round(R_bar) 
    return(rR_bar) 
}