Квадратный корень сингулярной матрицы в R

Question

Мне нужно вычислить матрицу A по степени -1/2, что в основном означает квадратный корень из инверсии исходной матрицы.

Если A является сингулярным, то обобщенное обратное значение Мура-Пенроуза вычисляется с помощью функции ginv из пакета MASS, в противном случае регулярное обратное вычисляется с использованием функции решения.

Матрица А определена ниже:

A <- structure(c(604135780529.807, 0, 58508487574887.2, 67671936726183.9, 
            0, 0, 0, 1, 0, 0, 0, 0, 58508487574887.2, 0, 10663900590720128, 
            10874631465443760, 0, 0, 67671936726183.9, 0, 10874631465443760, 
            11315986615387788, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1), .Dim = c(6L, 
                                                                                   6L))

Я проверяю сингулярность сравнением ранга и размерности.

rankMatrix(A) == nrow(A)

Приведенный выше код возвращает FALSE, поэтому я должен использовать ginv, чтобы получить обратное. Обратное значение A выглядит следующим образом:

A_inv <- ginv(A)

Квадратный корень обратной матрицы вычисляется с помощью функции sqrtm из пакета expm.

library(expm)
sqrtm(A_inv)

Функция возвращает следующую ошибку:

Ошибка в solve.default(X[ii, ii] + X[ij, ij], S[ii, ij] - sumU):
Лапак рутина zgesv: система точно единственная

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

Answer 1

Ваш вопрос касается двух разных проблем:

1. Обратная матрица
2. Квадратный корень из матрицы

обратный

Обратного не существует для особых матриц. В некоторых приложениях, подход Мура-Пенроуза или другое обобщенное обратное может быть взято в качестве подходящей замены обратного. Тем не менее, обратите внимание, что в большинстве случаев компьютерные числа будут вызывать ошибки округления; и эти ошибки могут привести к тому, что единственная матрица будет выглядеть регулярно для компьютера или наоборот.

Если A всегда показывает блочную структуру матрицы, которую вы даете, я предлагаю рассмотреть только ее недиагональный блок

A3 = A[ c( 1, 3, 4 ), c( 1, 3, 4 ) ]

A3
             [,1]         [,2]         [,3]
[1,] 6.041358e+11 5.850849e+13 6.767194e+13
[2,] 5.850849e+13 1.066390e+16 1.087463e+16
[3,] 6.767194e+13 1.087463e+16 1.131599e+16

вместо всего A для лучшей эффективности и меньшего количества вопросов округления. Оставшиеся 1-диагональные записи останутся 1 в обратном корне квадратного, поэтому нет необходимости загромождать вычисления с ними. Чтобы получить представление о влиянии этого упрощения, обратите внимание, что R может рассчитать

A3inv = solve(A3)

пока не смог вычислить

Ainv = solve(A)

Но нам не понадобится A3inverse, как станет ясно ниже.

Квадратный корень

Как правило, квадратный корень из матрицы A будет существовать только в том случае, если матрица имеет диагональную жорданову нормальную форму ( https://en.wikipedia.org/wiki/Jordan_normal_form). Следовательно, по-настоящему общего решения проблемы не существует.

К счастью, подобно тому, как "большинство" (вещественных или комплексных) матриц являются обратимыми, "большинство" (вещественных или комплексных) матриц имеют диагональную комплексную жорданову нормальную форму. В этом случае жорданова нормальная форма

A3 = T·J·T⁻¹

можно рассчитать в R как таковой:

X = eigen(A3)
T = X$vectors
J = Diagonal( x=X$values )

Чтобы проверить этот рецепт, сравните

Tinv = solve(T)
T %*% J %*% Tinv

с A3, Они должны совпадать (до ошибок округления), если A3 имеет диагональную жорданову нормальную форму.

поскольку J является диагональным, его квадратный корень является просто диагональной матрицей квадратных корней

Jsqrt = Diagonal( x=sqrt( X$values ) )

чтобы Jsqrt·Jsqrt = J, Более того, это подразумевает

(T·Jsqrt·T⁻¹)² = T·Jsqrt·T⁻¹·T·Jsqrt·T⁻¹ = T·Jsqrt·Jsqrt·T⁻¹ = T·J·T⁻¹ = A3

так что на самом деле мы получаем

√A3 = T·Jsqrt·T⁻¹

или в коде R

A3sqrt = T %*% Jsqrt %*% Tinv

Чтобы проверить это, рассчитайте

A3sqrt %*% A3sqrt

и сравнить с A3,

Квадратный корень из обратного

Квадратный корень обратного корня (или, в равной степени, обратный корень квадратного) можно легко вычислить после вычисления диагональной жордановой нормальной формы. Вместо J использование

Jinvsqrt = Diagonal( x=1/sqrt( X$values ) )

и рассчитать, аналогично приведенному выше,

A3invsqrt = T %*% Jinvsqrt %*% Tinv

и наблюдать

A3·A3invsqrt² = … = T·(J/√J/√J)·T⁻¹ = 1

единичная матрица, так что A3invsqrt является желаемым результатом.

В случае, если A3 не является обратимым, обобщенный обратный (не обязательно Moore-Penrose) может быть рассчитан путем замены всех неопределенных записей в Jinvsqrt на 0, но, как я уже сказал выше, это следует делать с должной осторожностью в свете общего применения и его устойчивости к ошибкам округления.

В случае, когда A3 не имеет диагональной йордановой нормальной формы, квадратного корня нет, поэтому приведенные выше формулы приведут к другому результату. Чтобы не столкнуться с этим случаем во времена неудачи, лучше всего выполнить тест

A3invsqrt %*% A3 %*% A3invsqrt

достаточно близко к тому, что вы бы посчитали матрицей 1 (это применимо только в том случае, если A3 был обратим в первую очередь).

PS: обратите внимание, что вы можете поставить префикс знака ± для каждой диагональной записи Jinvsqrt на ваш вкус.