如何测试矩阵是否为对角矩阵？

Question

我需要测试一个方差矩阵是否为对角矩阵。如果不是，那么我将进行Cholesky LDL分解。但我想知道最可靠和最快的方法来测试矩阵是否为对角矩阵？我正在使用Fortran。

我首先想到的是将矩阵的所有元素求和，并从该总和中减去对角线元素。如果答案为0，则矩阵为对角矩阵。您有更好的建议吗？

在Fortran中，我会编写以下代码：

!A is my matrix
k=0.0d0
do i in 1:n #n is the number of rows/colums
k = k + A(i,i)
end do

if(abs(sum(A)-k) < epsilon(k)*sum(A)) then
#do cholesky LDL, which I have to write myself, haven't found any subroutines for that  in Lapack or anywhere else
end if

- Jouni

---

仅仅挑剔一下：你的意思是LDL分解，而不是LDL。;-) - Stobor

另外：LAPACK LDL'分解：http://www.netlib.org/lapack/single/ssptrf.f LAPACK Cholesky LL'分解：http://www.netlib.org/lapack/single/spotrf.f - Stobor

非常感谢你提出的所有观点！ :D - Jouni

那个LAPACK LDL'分解只有单精度，我需要双精度... - Jouni

顺便说一下，那个反例在这里不起作用，A是协方差矩阵，因此其中没有负值。但正如所说，Sharptooth的算法更好。这只是我想到的第一个想法。很高兴我问了。 :) - Jouni

Answer 1

更好的方法是遍历所有的非对角线元素，并测试它们是否接近于零（直接比较浮点数是否相等容易受到舍入误差影响，导致错误的结果）。

首先，一旦发现任何违反条件的元素，你可以立即停止遍历，如果存在违反条件的矩阵，则这可以显著减少时间。

其次，这可能使得编译器更好地展开循环（Fortran 编译器以良好的优化策略而闻名），并且由于指令之间的依赖关系更少，因此可以在芯片上更快地执行。

此外，你提出的算法容易出现溢出和误差累积问题，而“遍历和测试”算法则没有这些问题。

- sharptooth

---

非常感谢，我会按照那种方式去做。 - Jouni

A是对称的，所以我只需要遍历矩阵的一半... 再次感谢，我不是一个“真正”的程序员，所以不能说关于多处理器、指令间依赖等方面的事情。 :D 只是一个统计学家。 ;) - Jouni

最好的方法是遍历所有非对角线元素并测试它们是否为零，将其更改为“如果它们小于epsilon”。很少会出现零值，通常是0.0002332等。 - Rook

好的观点。我是这样比较元素的： 如果(abs(Ht(k,i)) > epsilon(0.0d0)*abs(Ht(k,i)))，那么 不确定我是否正确理解了epsilon函数？我需要在那里使用abs(Ht(k,i))吗？我只是读了一些例子，我真的不理解那个乘法的真正含义。 - Jouni

@Jouni - 不，我不是指Epsilon作为最小的处理器相关值（根据我的了解，这似乎是你在尝试做的）。我的意思更多的是像“如果(element_value.GT.greatest_element_value_in_matrixeps) THEN”，其中eps=10*-5。... 意思是，如果某个元素比主对角线上的元素小5个数量级，则可以视为零（因此不会使矩阵任何更少对角线化——这只是数值过程剩余的微小差异）。 - Rook

这些注释应该有某种格式，我个人认为。就是为了这种东西。 - Rook

Answer 2

在矩阵中搜索非零值

logical :: not_diag
integer :: i, j

not_diag = .false.

outer: do i = 2, size(A,1)
  do j = i, size(A, 2)
    if (A(i,j) > PRECISION) then
      not_diag = .true.
      exit outer
    end if
  end
end outer

if (not_diag) then
  ! DO LDL' decomposition
end if

要使用双精度LAPACK例程，请将第一个's'替换为'd'。 因此，spotrf变成了dpotrf。

http://www.netlib.org/lapack/double/

- Ian

---

是的，但没有dsptrf.f，只有ssptrf.f可以进行LDL'分解。dpotrf可以进行LL'分解。 - Jouni