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

Question

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

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]])

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

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 не подходит для единичных матриц. Любые предложения с благодарностью. Благодарю.

Answer 1

0 / lu_fac[0][9, 9]

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

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

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

имеет вектор б, который находится в пространстве диапазона А

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

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

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

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

Answer 2

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