Несогласованные результаты между декомпозицией LU в R и Python

У меня есть следующая матрица A в R:

           # [,1]       [,2]   [,3]       [,4]
# [1,] -1.1527778  0.4444444  0.375  0.3333333
# [2,]  0.5555556 -1.4888889  0.600  0.3333333
# [3,]  0.6250000  0.4000000 -1.825  0.8000000
# [4,]  0.6666667  0.6666667  0.200 -1.5333333

A <- structure(c(-1.15277777777778, 0.555555555555556, 0.625, 0.666666666666667, 
0.444444444444444, -1.48888888888889, 0.4, 0.666666666666667, 
0.375, 0.6, -1.825, 0.2, 0.333333333333333, 0.333333333333333, 
0.8, -1.53333333333333), .Dim = c(4L, 4L), .Dimnames = list(NULL, 
    NULL))

Я вычисляю его LU-разложение следующим образом:

library(Matrix)
ex <- expand(lu(t(A)))
L <- ex$L
P <- ex$P
C <- U <- ex$U
C[lower.tri(U)] <- L[lower.tri(L)]

print(C)

# 4 x 4 Matrix of class "dgeMatrix"
           # [,1]       [,2]       [,3]          [,4]
# [1,] -1.1527778  0.5555556  0.6250000  6.666667e-01
# [2,] -0.3855422 -1.2746988  0.6409639  9.236948e-01
# [3,] -0.3253012 -0.6124764 -1.2291115  9.826087e-01
# [4,] -0.2891566 -0.3875236 -1.0000000 -2.220446e-16

С другой стороны, это код Python для той же работы:

lu, piv = scipy.linalg.lu_factor(A.T, check_finite=False)

print(lu)

# [[ -1.15277778e+00   5.55555556e-01   6.25000000e-01   6.66666667e-01]
 # [ -3.85542169e-01  -1.27469880e+00   6.40963855e-01   9.23694779e-01]
 # [ -2.89156627e-01  -3.87523629e-01   1.22911153e+00  -9.82608696e-01]
 # [ -3.25301205e-01  -6.12476371e-01  -1.00000000e+00   7.69432827e-16]]

Мне интересно, почему две матрицы C и lu в R и Python (соответственно) не совпадают. Дело в том, что я должен получить тот же результат, что и версия Python (т.е. матрица lu). Вы хоть понимаете, что я делаю неправильно?


python r matrix matrix-decomposition
person 989    schedule 10.01.2017    source источник

Ответы (1)


arrow_upward
6
arrow_downward

Стыдно осознавать спустя 1,5 года, что мой первоначальный ответ не совсем верен. Хотя правильно указано, что причиной является недостаток ранга A в вопросе, неправильно приписывать это как основную причину. Настоящая проблема заключается в неуникальном выборе опорной точки, что может произойти (хотя и менее вероятно), даже если A имеет полный ранг. Учитывая, что этот пост был просмотрен более 700 раз и имеет 6 баллов, я мог ввести в заблуждение многих читателей. ИЗВИНИТЕ!

Я отправил Написать отслеживаемую функцию R, которая имитирует dgetrf LAPACK для факторизации LU, и только что ответил на нее. Вопрос содержит функцию LU для факторизации LU без поворота, а ответ содержит две функции LUP и LUP2 для версии с поворотом строки, которые согласуются с dgetrf LAPACK, лежащим в основе плотного метода Matrix::lu и R. базовая функция solve. В частности, функция LUP2 позволяет шаг за шагом отслеживать факторизацию. Я буду использовать эту функцию здесь для расследования.

Таким образом, вы факторизуете транспонирование A.

Из вывода R и Python мы видим, что они дают идентичные 1-й поворот -1.1527778 и 2-й поворот -1.2746988, в то время как 3-й поворот отличается. Таким образом, определенно происходит что-то интересное, когда факторизация выполнила первые два столбца/строки и переходит к третьему столбцу/строке. Давайте приостановим факторизацию LU в R на этом этапе:

oo <- LUP2(t(A), to = 2)
#           [,1]       [,2]       [,3]       [,4]
#[1,] -1.1527778  0.5555556  0.6250000  0.6666667
#[2,] -0.3855422 -1.2746988  0.6409639  0.9236948
#[3,] -0.3253012 -0.6124764 -1.2291115  0.9826087
#[4,] -0.2891566 -0.3875236  1.2291115 -0.9826087
#attr(,"to")
#[1] 2
#attr(,"pivot")
#[1] 1 2 3 4

В этот момент исключение Гаусса уменьшило t(A) до

getU(oo)
#           [,1]       [,2]       [,3]       [,4]
#[1,] -1.1527778  0.5555556  0.6250000  0.6666667
#[2,]  0.0000000 -1.2746988  0.6409639  0.9236948
#[3,]  0.0000000  0.0000000 -1.2291115  0.9826087
#[4,]  0.0000000  0.0000000  1.2291115 -0.9826087
#attr(,"to")
#[1] 2

Ух ты, теперь мы видим кое-что действительно интересное: 3-я и 4-я строки отличаются только изменением знака. Тогда исключение Гаусса не является уникальным, потому что либо -1.2291115, либо 1.2291115 могут быть точкой опоры, поскольку они имеют одно и то же абсолютное значение.

Ясно, что R выбрал -1.2291115 в качестве опорной точки, но Python выбрал в качестве опорной точки 1.2291115. В Python произойдет обмен строками между 3-й и 4-й строками. В своем вопросе вы не сообщили, какой индекс перестановки дает вам Python, но он должен быть 1, 2, 4, 3 вместо 1, 2, 3, 4 в R.

В любом случае множитель U заканчивается строкой нулей внизу, так что t(A) или A не являются полными. Если вы хотите увидеть более согласованное поведение между двумя программами, вам лучше попробовать полноранговую матрицу. В этом случае маловероятно, что у вас будет более одного варианта разворота во время факторизации LU. Вы можете сгенерировать случайную матрицу полного ранга в R с помощью

A <- matrix(runif(16), 4, 4)
person Zheyuan Li    schedule 10.01.2017
comment
Ты прав. Учитывая матрицу полного ранга, я получаю тот же результат. - person 989; 10.01.2017