0

我需要计算矩阵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)

上面的代码返回 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 并不总是奇异的,因此我们必须为该问题提供一个通用解决方案。

4

1 回答 1

6

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

  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⁻¹</p>

可以在 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⁻¹</p>

或在 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 不可逆,则可以通过将所有未定义的条目替换为 0 来计算广义逆(不一定是 Moore-Penrose 逆)Jinvsqrt,但正如我上面所说,这应该根据总体情况谨慎完成应用程序及其对舍入误差的稳定性。

如果 A3 没有对角 Jordan 范式,则没有平方根,因此上述公式将产生一些其他结果。为了在运气不好的时候不会遇到这种情况,最好执行一个测试是否

A3invsqrt %*% A3 %*% A3invsqrt

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

JinvsqrtPS:请注意,您可以根据自己的喜好为每个对角线条目添加符号 ± 。

于 2015-08-10T11:44:59.640 回答