R 中奇异矩阵的平方根

Question

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

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

矩阵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)

上面的代码returns FALSE，所以我不得不用ginv求逆。 A的逆运算如下：

A_inv <- ginv(A)

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

library(expm)
sqrtm(A_inv)

函数returns出现如下错误：

Error in solve.default(X[ii, ii] + X[ij, ij], S[ii, ij] - sumU) :
Lapack routine zgesv: system is exactly singular

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

Answer 1

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

逆矩阵
矩阵的平方根

逆

奇异矩阵不存在逆矩阵。在某些应用中，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 的每个对角线条目加上符号 ±。

R 中奇异矩阵的平方根

Square Root of a Singular Matrix in R

r

matrix

inverse

square-root

singular