在 C# 中替换 IMSL MCHOL (Fortran) 中的 Cholesky 分解

发布于 2024-12-11 00:53:37 字数 1850 浏览 0 评论 0原文

我正在将 Fortran 程序转换为 C#。这必须一点一点地完成，并一路上进行概念验证。

这些初始步骤之一是复制它使用的 IMSL 函数。幸运的是，它只使用了选定的几个：一些平凡的随机数生成，一些平凡的正态分布反转，但也有一个不那么平凡的：MCHOL。

来自文档：

计算实数对称矩阵 A 加对角矩阵 D 的上三角分解，其中 D 在 Cholesky 过程中按顺序确定因式分解以使 A + D 非负定。

例程 MCHOL 计算 A + D 的 Cholesky 因式分解，RTR，其中 A 是对称的，D 是具有足够大对角线元素的对角矩阵，使得 A + D 是非负定的。该例程与 Gill、Murray 和 Wright（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 的某些部分，特别是 spdmatrixcholesky。它似乎有效，但仅适用于矩阵的一部分：

  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

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

由于似乎不可能快速给出答案，因此如果我了解基础数学可能是最好的，这样我至少可以尝试弄清楚这里发生了什么。任何理论基础或指针也值得赞赏。

原文

I am converting a Fortran program to C#. This has to be done in bits and pieces, with proof of concepts along the way.

One of these initial steps is to replicate the IMSL functions it uses. Luckily, it only uses a chosen few: some trivial random number generation, some trivial normal distribution inverts but also one not-so-trivial one: MCHOL.

From the documentation:

Computes an upper triangular factorization of a real symmetric matrix A plus a diagonal matrix D, where D is determined sequentially during the Cholesky factorization in order to make A + D nonnegative definite.

Routine MCHOL computes a Cholesky factorization, RTR, of A + D where A is symmetric and D is a diagonal matrix with sufficiently large diagonal elements such that A + D is nonnegative definite. The routine is similar to one described by Gill, Murray, and Wright (1981, pages 108−111). Here, though, we allow A + D to be singular.

(There are more details and a sample at the link).

For my proof of concept, I need to able to replicate the results as provided in the sample for the MCHOL documentation: passing in this matrix from the sample:

  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/

And getting the following in return:

  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

So far, I've tried using Math.NET, but it will not run on this example matrix because it is not positive definite.

I've also tried parts of ALGLIB, specifically spdmatrixcholesky. It seems to work, but only for part of the matrix:

  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

Does anyone have any idea what I am doing wrong here? Do I need to call a different function here?

Since a quick answer doesn't seem to be in the cards, it might be best if I understood the underlying mathematics so I can at least try to figure out what is going on here. Any theoretical underpinning or pointers appreciated as well.

分享到QQ

分享到微博