R中奇异矩阵的平方根

Question

我需要计算矩阵 A 的 -1/2 次方，这基本上意味着初始矩阵逆的平方根。

如果 A 是奇异的，则使用 MASS 包中的 ginv 函数计算 Moore-Penrose 广义逆，否则使用 solve 函数计算正则逆。

矩阵A定义如下：

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)

逆矩阵的平方根是使用 expm 包中的 sqrtm 函数计算的。

library(expm)
sqrtm(A_inv)

函数返回如下错误：

solve.default(X[ii, ii] + X[ij, ij], S[ii, ij] - sumU) 中的错误：
Lapack 例程 zgesv：系统完全是奇异的

那么在这种情况下我们如何计算平方根呢？请注意，矩阵 A 并不总是奇异的，因此我们必须为该问题提供一个通用解决方案。

Answer 1

您的问题涉及两个不同的问题：

1. 矩阵的逆
2. 矩阵的平方根

逆向

奇异矩阵不存在逆。在某些应用中，Moore-Penrose 或一些其他广义逆可以被视为逆的合适替代。但是，请注意，在大多数情况下，计算机数字会产生舍入误差；这些错误可能会使奇异矩阵在计算机看来是规则的，反之亦然。

如果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 的平方根仅在矩阵具有对角 Jordan 范式 (https://en.wikipedia.org/wiki/Jordan_normal_form) 时才存在。因此，没有您需要的问题的真正通用解决方案。

幸运的是，就像“大多数”（实数或复数）矩阵是可逆的一样，“大多数”（实数或复数）矩阵具有对角复数 Jordan 范式。在这种情况下，乔丹范式

A3 = T·J·T⁻¹

可以在 R 中这样计算：

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

要测试这个食谱，请比较

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

A3。如果 A3 具有对角 Jordan 范式，则它们应该匹配（直至舍入误差）。

由于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比较。

逆的平方根

一旦计算了对角 Jordan 范式，就可以很容易地计算出逆的平方根（或者，同样地，平方根的倒数）。而不是J 使用

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

并计算，类似于上面，

A3invsqrt = T %*% Jinvsqrt %*% Tinv

观察

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

单位矩阵，因此 A3invsqrt 是所需的结果。

如果 A3 不可逆，则可以通过将 Jinvsqrt 中的所有未定义条目替换为 0 来计算广义逆（不一定是 Moore-Penrose 逆），但正如我上面所说，这应该小心谨慎考虑到整个应用程序及其对舍入误差的稳定性。

如果 A3 不存在对角 Jordan 范式，则不存在平方根，因此上述公式会产生其他结果。为了在运气不好的时候不会遇到这种情况，最好进行测试

A3invsqrt %*% A3 %*% A3invsqrt

足够接近您认为的 1 矩阵（这仅适用于 A3 首先是可逆的）。

PS：请注意，您可以根据自己的喜好为 Jinvsqrt 的每个对角线条目添加一个符号 ±。