为什么我不能在犰狳中得到这个对称正定矩阵的平方根?

Why can't I get the square root of this symmetric positive-definite matrix in armadillo?

我有矩阵 A:

0 0 0 0 1

0 0 0 0 0 1

0 0 0 0 0 1

0 0 0 0 0 1

1 0 0 0 0 0

1 1 1 1 0 0.


这个矩阵显然是对称的,所以A*A显然是正定矩阵。我试图计算它在犰狳中的平方根但失败了。我做错了什么?

代码如下:

#include <iostream>
#include <armadillo>

using namespace arma;
using namespace std;

int main(){
    mat A(6,6);
    for(int i=0;i<6;i++){ // this part just reads the matrix,
        for(int j=0;j<6;j++){
            int x;
            scanf("%d",&x);
            A(i,j)=x;
        }
    }
    cout << sqrtmat_sympd(A*A);
}

我编译了它,它适用于一些矩阵,比如单位矩阵,但是当我尝试使用这个矩阵时,它说:

error: sqrtmat_sympd(): transformation failed
terminate called after throwing an instance of 'std::runtime_error'
  what():  sqrtmat_sympd(): transformation failed
Aborted (core dumped)

它适用于某些确定的 semi-positive 矩阵而不适用于其他矩阵的原因似乎是有时本应为零的特征值结果是范数较小的负数。这会导致问题,因为该算法很可能会在某个地方取特征值的平方根。

为了能够得到对称定semi-positive对称矩阵的平方根,我们可以使用特征值分解特征,然后取对角线部分的平方根。(注意改变"negative zeros" 到实际零。

下面是执行此操作的一些代码:

mat raiz(const mat& A){
  vec D;
  mat B;
  eig_sym(D,B,A);

  unsigned int n = D.n_elem;

  mat G(n,n,fill::zeros);

  for(unsigned int i=0;i<n;i++){
      if(D(i)<0 ) D(i)=0;
      G(i,i)=sqrt(D(i));
  }

  return B*G*B.t();
}

注意sqrtmat() 函数可能更合适,因为它已经尝试近似矩阵平方根。