Ваш вопрос касается двух разных проблем:
- Обратная матрица
- Квадратный корень из матрицы
Обратный
Обратного не существует для сингулярных матриц. В некоторых приложениях подход Мура-Пенроуза или другой обобщенный обратный ход можно использовать в качестве подходящей замены обратного. Однако учтите, что компьютерные числа в большинстве случаев приводят к ошибкам округления; и эти ошибки могут сделать сингулярную матрицу регулярной для компьютера или наоборот.
Если 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.
Квадратный корень из обратного
Квадратный корень из обратного (или, что то же самое, из обратного корня sqare) может быть легко вычислен после вычисления диагональной жордановой нормальной формы. Вместо J используйте
Jinvsqrt = Diagonal( x=1/sqrt( X$values ) )
и вычислим, аналогично предыдущему,
A3invsqrt = T %*% Jinvsqrt %*% Tinv
и наблюдать
A3 · A3invsqrt² =… = T · (Дж / √J / √J) · T⁻¹ = 1
единичная матрица, так что A3invsqrt является желаемым результатом.
В случае, если A3 не является обратимым, можно вычислить обобщенное обратное (не обязательно Мура-Пенроуза), заменив все неопределенные записи в Jinvsqrt на 0, но, как я сказал выше, это следует делать с должной осторожностью в свете общее приложение и его устойчивость к ошибкам округления.
В случае, если A3 не имеет диагональной жордановой нормальной формы, квадратный корень отсутствует, поэтому приведенные выше формулы дадут другой результат. Чтобы не столкнуться с этим случаем во время неудач, лучше всего выполнить проверку,
A3invsqrt %*% A3 %*% A3invsqrt
достаточно близко к тому, что вы бы считали матрицей 1 (это применимо только в том случае, если A3 изначально была обратимой).
PS: Обратите внимание, что вы можете поставить знак ± для каждой диагональной записи Jinvsqrt по своему вкусу.
person
Bernhard Bodenstorfer
schedule
10.08.2015
