Как решить особые матрицы?

Question

Мне необходимо solve более тысячи матриц в списке. Однако я получаю ошибку Lapack routine dgesv: system is exactly singular, Моя проблема в том, что мои входные данные не являются единичными матрицами, но после вычислений, которые мне нужно сделать на матрицах, некоторые из них становятся единственными. Воспроизводимый пример с подмножеством моих данных, однако, не возможен, поскольку это было бы слишком долго (я уже пробовал). Вот базовый пример моей проблемы (A будет матрица после некоторых вычислений и R следующий расчет, который мне нужно сделать):

A=matrix(c(1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1), nrow=4)
R = solve(diag(4)-A)

У вас есть идеи, как "решить" эту проблему, может быть, другую функцию? Или как очень мало преобразовать особые матрицы, чтобы не создавать совершенно разных результатов? Спасибо

РЕДАКТИРОВАТЬ: в соответствии с @Roman Susi я включаю функцию (со всеми расчетами), я должен сделать:

function(Z, p) {
  imp <- as.vector(cbind(imp=rowSums(Z)))
  exp <- as.vector(t(cbind(exp=colSums(Z))))
  x = p + imp
  ac = p + imp - exp  
  einsdurchx = 1/as.vector(x)    
  einsdurchx[is.infinite(einsdurchx)] <- 0
  A = Z %*% 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)
}

Проблема в строке 8 функции расчета solve(diag(length(p))-A), Здесь я могу предоставить новые примеры данных для Z а также pОднако в этом примере он работает нормально, так как я не смог воссоздать пример, который приводит к ошибке:

p = c(200, 1000, 100, 10)
Z = matrix(
  c(0,0,100,200,0,0,0,0,50,350,0,50,50,200,200,0),
  nrow = 4, 
  ncol = 4,
  byrow = T) 

Итак, вопрос в соответствии с @Roman Susi: есть ли способ изменить расчеты до того, чтобы det(diag(length(p))-A) никогда не получает 0 для того, чтобы solve уравнение? Я надеюсь, что вы можете понять, что я хочу:) Идеи, спасибо. Edit2: Может быть, спросить проще: как избежать сингулярности в этой функции (по крайней мере, до строки 8)?

Answer 1

Обобщенный обратный ginv в пакете MASS можно обрабатывать единичные матрицы, но нужно будет определить, имеет ли это смысл для вашей проблемы.

 library(MASS)
 ginv(diag(4) - A)

давая:

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0
[4,]    0    0    0    0

Также есть Ginv функция в пакете ibdreg.

Answer 2

Функции QR-разложения R могут иметь ваш ответ. Они обеспечивают надежное решение линейных уравнений. QR-разложение не обеспечивает обратное, а скорее матричное разложение, которое часто можно использовать там, где будет использоваться обратное.

Для прямоугольных матриц QR-разложение может использоваться для нахождения наименьших квадратов. Для квадратных ((близких) сингулярных матриц, qr() обнаруживает эту (почти) особенность, и qr.coef() затем можно использовать для получения решения без каких-либо ошибок, но, возможно, с некоторыми NA (которые можно преобразовать в нули).