Как решать сингулярные матрицы?

Мне нужно 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)

Есть ли у вас идеи, как «решить» эту проблему, может быть, другую функцию? Или как преобразовать сингулярные матрицы очень-очень незначительно, чтобы не создавать совершенно разные результаты? Спасибо

EDIT: согласно @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)?


r matrix singular
person N.Varela    schedule 06.11.2015    source источник
comment
Из вопроса неясно, является ли это математической задачей, а не программированием. Если инвертируемая сингулярная матрица является промежуточным результатом, то, может быть, другая формула позволит избежать деления на ноль? Потому что, если сингулярности нельзя избежать, она заложена в задаче, и задачу следует переформулировать шире или повысить точность и т. д., зависит от этого.   -  person Roman Susi    schedule 06.11.2015
comment
Когда вы передаете матрицу solve, вы говорите ей найти обратную матрицу, но сингулярные матрицы не имеют обратной. Я думаю, вы просите его сделать невозможное....   -  person Jthorpe    schedule 06.11.2015

Ответы (2)


arrow_upward
3
arrow_downward

Обобщенная инверсия 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

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

person G. Grothendieck    schedule 06.11.2015
comment
Это не обратное в общем смысле - просто удобство/неудобство программирования (я всегда предпочел бы ошибку, чем ответ, который может иметь или не иметь смысла!) - person Jthorpe; 06.11.2015
comment
... или, возможно, нужно будет переопределить, имеет ли это смысл или нет. Извините, не удержался. - person Brian B; 06.11.2015
comment
@Jthorpe, это инверсия подмножества домена, состоящего из ортогонального дополнения к нулевому пространству. Пример, приведенный в вопросе, не особенно иллюстративен, поскольку в этом случае нулевое пространство - это все пространство, но в целом это будет не так. - person G. Grothendieck; 07.11.2015

arrow_upward
0
arrow_downward

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

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

person Paul Harrison    schedule 10.12.2016