Получаете удовольствие от использования lu-факторизации для решения сингулярной квадратной матрицы?

У меня есть сингулярная матрица A (10 * 10) с недостаточным рангом (ранг = 9), и у меня есть вектор b, который находится в пространстве диапазонов A. Теперь меня интересует какое-то решение для Ax = b. Для конкретности вот мой A

array([[ 0.        ,  0.        ,  0.        ,  0.86826141,  0.        ,
             0.        ,  0.88788426,  0.        ,  0.4089203 ,  0.88134901],
           [ 0.        ,  0.        ,  0.46416372,  0.        ,  0.        ,
             0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
           [ 0.        ,  0.        ,  0.        ,  0.        ,  0.31303966,
             0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
           [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
             0.        ,  0.3155742 ,  0.        ,  0.64059294,  0.        ],
           [ 0.        ,  0.        ,  0.        ,  0.        ,  0.51349938,
             0.        ,  0.        ,  0.        ,  0.53593509,  0.        ],
           [ 0.        ,  0.01252787,  0.        ,  0.6870415 ,  0.        ,
             0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
           [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
             0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
           [ 0.        ,  0.        ,  0.        ,  0.16643105,  0.        ,
             0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
           [ 0.08626592,  0.        ,  0.        ,  0.        ,  0.        ,
             0.        ,  0.        ,  0.        ,  0.        ,  0.66939531],
           [ 0.43694586,  0.        ,  0.        ,  0.        ,  0.        ,
             0.95941661,  0.        ,  0.52936733,  0.79687149,  0.81463887]])

b генерируется с использованием A.dot(np.ones(10)). Теперь я хотел решить эту проблему с помощью факторизации lu, и для этого я сделал следующее.

lu_fac=scipy.linalg.lu_factor(X)
scipy.linalg.lu_solve(lu_fac,b)

Который дает

array([ nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan])

Кроме того, lu_factor в этом случае работает нормально (иногда он выдает предупреждение во время выполнения, говорящее "Диагональное число %d равно нулю. Сингулярная матрица"). Для полноты вот код для проверки того, что PLU из lu_factor такой же, как A :

L=np.tril(lu_fac[0])
np.fill_diagonal(L,1)
U=np.triu(lu_fac[0])
perm=np.arange(10)
ipiv=lu_factor[1]
for i in range(10):
  temp=perm[i]
  perm[i]=perm[ipiv[i]]
  perm[ipiv[i]]=temp
np.allclose(X[perm,:],L.dot(U))

Теперь я знаю, что моя матрица сингулярна и у моей проблемы бесконечно много решений. Но меня интересует любое решение, и я просто запутался, почему lu-факторизация терпит неудачу, не может ли она установить свободные переменные в 0 и найти какое-то решение, как нас учат? Также в чем дело с предупреждением во время выполнения "Диагональное число %d равно нулю. Сингулярная матрица". Примечание. Меня не интересует подход svd/qr для решения этой проблемы, мне просто любопытно узнать, почему lu не работает для сингулярных матриц. Любые предложения приветствуются. Спасибо.


python scipy linear-algebra matrix-factorization
person user1131274    schedule 03.10.2017    source источник

Ответы (2)


arrow_upward
1
arrow_downward 

0 / lu_fac[0][9, 9]

возвращает nan, потому что эта запись — последняя диагональная запись U — равна нулю. Таким образом, это nan становится значением 9-й переменной. Затем оно подставляется в приведенные выше уравнения, и, естественно, остальное тоже получается как nan. Код SciPy LU, или, скорее, код Fortran, который он обертывает, не предназначен для матриц с недостаточным рангом, поэтому он не собирается создавать значения для переменных, которые не могут быть определены.

Также в чем дело с предупреждением во время выполнения «Диагональное число% d равно нулю. Сингулярная матрица».

Предупреждение ясно: алгоритм обнаружил сингулярную матрицу, чего не ожидалось. Это также говорит вам, что реализация не предназначена для использования с сингулярными матрицами.

есть вектор b, который находится в пространстве диапазонов A

Это теоретически. На практике нельзя быть уверенным в том, что что-либо находится в пространстве диапазонов матрицы с недостаточным рангом из-за ошибок, присущих арифметике с плавающей запятой. Вы можете вычислить b = A.dot(...), а затем попытаться решить Ax=b, и решения не будет из-за ошибок, возникающих при манипулировании числами с плавающей запятой.

Кстати: вы упомянули, что факторизация PLU существует для каждой квадратной матрицы, но SciPy не обязательно предназначена для ее вычисления. Например,

scipy.linalg.lu_factor(np.array([[0, 1], [0, 0]]))

возвращает матрицу с NaN. В вашем случае NaN появляются позже, при попытке найти решение и встречаются с нулевым диагональным элементом фактора U.

person Community    schedule 04.10.2017
comment
Значит ли это, что LU не рекомендуется для сингулярных матриц? И если я правильно понимаю вопросы не зная что-то равно 0 или близко к нулю? Большое спасибо. - person user1131274; 04.10.2017
comment
Да. LU предназначен для точного решения системы, что неправильно пытаться делать с сингулярными матрицами. Это контрастирует с методами, включающими ортогональные разложения, которые представляют собой надежный подход, основанный на минимизации суммы квадратов. - person ; 04.10.2017

arrow_upward
0
arrow_downward

Как упоминалось здесь, матрица имеет Факторизация LU тогда и только тогда, когда rank(A11) + k >= rank([A11 A12]) + rank([A11 A21]). В вашем случае rank(A11) = 3, k = 5 и rank([A11 A12]) + rank([A11 A21]) = 9. Таким образом, ваша матрица не удовлетворяет условиям и не имеет факторизации LU.

person Kevin K.    schedule 03.10.2017
comment
Спасибо за ваш ответ. Но lu_factor выполняет декомпозицию PLU, которая существует всегда. Также я проверил, что PLU, который я получаю от lu_factor, такой же, как A. Так что это не так. При выполнении lu_solve ничего не получается. Спасибо еще раз. - person user1131274; 04.10.2017