用C#中的IMSL (Fortran)代替Cholesky分解

Question

我正在将Fortran程序转换为C#。这必须以零碎的方式进行，并在此过程中提供概念的证明。

这些初始步骤之一是复制它使用的IMSL函数。幸运的是，它只使用了选择的少数:一些平凡的随机数生成，一些平凡的正态分布倒置，但也有一个不平凡的: MCHOL。

来自文档：

计算实对称矩阵A+对角矩阵D的上三角分解，其中D是在Cholesky分解过程中顺序确定的，从而使A+D非负定。

常规MCHOL计算A+D的Cholesky因式分解( RTR )，其中A是对称的，D是具有足够大的对角线元素的对角线矩阵，使得A+D是非负定的。这个例程类似吉尔，穆雷和赖特所描述的(1981年，第108页−111)。在这里，我们允许A+D是单数的。

(链接中有更多的细节和示例)。

为了证明我的概念，我需要能够复制MCHOL文档的示例中提供的结果:从示例中传递这个矩阵：

  DATA (A(1,J),J=1,N)/36.0, 12.0, 30.0, 6.0, 18.0/
  DATA (A(2,J),J=1,N)/12.0, 20.0, 2.0, 10.0, 22.0/
  DATA (A(3,J),J=1,N)/30.0, 2.0, 29.0, 1.0, 7.0/
  DATA (A(4,J),J=1,N)/6.0, 10.0, 1.0, 14.0, 20.0/
  DATA (A(5,J),J=1,N)/8.0, 22.0, 7.0, 20.0, 40.0/

并得到以下回报：

  6.000   2.000   5.000   1.000   3.000
  0.000   4.000  -2.000   2.000   4.000
  0.000   0.000   0.000   0.000   0.000
  0.000   0.000   0.000   3.000   3.000
  0.000   0.000   0.000   0.000   2.449

到目前为止，我已经尝试过使用Math.NET，但是它不会在这个示例矩阵上运行，因为它不是正定的。

我也尝试过部分ALGLIB，特别是[医]鸡胆。这似乎是可行的，但只适用于矩阵的一部分：

  6       2       5       1       3
  12      4       -2      2       4
  30      2       0       1       7
  6       10      1       14      20
  8       22      7       20      40

有人知道我在这里做错了什么吗？我需要在这里调用不同的函数吗？

因为一个快速的答案似乎并不在卡片中，如果我理解了基本的数学，那么我至少可以尝试找出这里发生了什么，这可能是最好的。任何理论基础或指针也会受到赞赏。

Answer 1

MCHOL例程不仅仅是做Cholesky分解，它还执行Cholesky的步骤，并跟踪让它继续进行的D对角线。它是“改良的”Cholesky。普通的Cholesky需要一个正定的输入，而这个例子不是。

这是MATLAB中的MCHOL。我将使用这个实现来构建一个.NET版本。

维基百科:Cholesky分解

Answer 2

在你的例子中，a(5,1)=8应该等于a(1,5)=18，所以你的初始矩阵是不对称的。