MATLAB中巨大矩阵的行列式答案

【问题标题】：Determinants of huge matrices in MATLABMATLAB中巨大矩阵的行列式
【发布时间】：2011-05-19 16:07:15
【问题描述】：

从一个模拟问题，我想在 MATLAB 中计算 1000x1000 量级的复杂方阵。由于这些值是指贝塞尔函数的值，因此矩阵一点也不稀疏。

由于我对行列式相对于某些参数（在我的例子中是搜索的特征函数的能量）的变化感兴趣，我通过首先搜索研究范围的重新缩放因子然后计算来克服目前的问题决定因素，

result(k) = det(pre_factor*Matrix{k});

现在这是一个非常尴尬的解决方案，仅适用于最大 500x500 的矩阵尺寸。

有人知道解决这个问题的好方法吗？与 Mathematica 的接口原则上可能有效，但我对可行性表示怀疑。提前谢谢你

罗伯特

编辑：我没有找到计算问题的方便解决方案，因为这需要更改为更高的精度。相反，我使用了那个

ln det M = trace ln M

也就是说，当我根据 k 导出它时

A = trace(inv(M(k))*dM/dk)

所以我至少有行列式关于k的对数的变化。从问题的物理背景，我可以得出对 A 的约束，这最终给了我一个对我的问题有效的解决方法。不幸的是，我不知道这样的解决方法是否可以推广。

【问题讨论】：

您的矩阵结构有什么特别的地方可能会有所帮助吗？你提到 k = {A1, B1; A2, B2} 其中 A1, A2 是对称的。还要别的吗？ A1 和 A2 是否通勤，或者是否有任何子矩阵很容易反转？
@Jonathan Dursi：感谢您提出深思熟虑的问题。但我担心，总的来说，我看不出为什么 A1 和 A2 应该通勤，以及为什么任何子矩阵都应该很容易可逆。此外，由于有趣的情况接近 det(M) = 0，因此反演效果不佳。
单数矩阵不一定能说明子矩阵的任何有趣之处；如果 B1 和 B2 为零，即使 A1 和 A2 的行为非常好，M 也不会可逆。
@Jonathan Dursi：我想，我明白了。由于所有涉及的矩阵都已完全填充，因此没有理由说明反演在计算上应该很容易。

【解决方案1】：

这不是严格意义上的 matlab 解决方案，但您可能需要考虑使用 Mahout。它专为大规模线性代数而设计。（1000x1000 对于它习惯的比例是没有问题的。）

您可以call into java 向/从 Mahout 传递数据。

【讨论】：

感谢您的建议。 Mahout 项目确实看起来很有希望。当然，这值得一试。但老实说，我对分布式计算机架构（我刚刚遇到的 Apache Hadoop）的了解非常有限，以至于我认为我需要太多时间来设置环境。
@Robert：Hadoop 可以跨多个节点运行，但并非必须如此。我花了大约 5 分钟来设置它以在我的机器上运行。 Matlab 也有一个parallel computing toolbox，如果你想加快 matlab 的速度（虽然我不认为它更精确）。不过我从未尝试过，所以我无法证明它的有效性。
感谢您的鼓励。我将与我们的管理员核实是否可以在某些节点上安装它。不过，您认为 Mahout 可以做的不仅仅是我猜给定问题所需要的双精度吗？我搜索了一下，但找不到有关它的信息（例如 Doc-Page down atm）。真诚的
@Robert：好点，它看起来不像 Mahout 那样。也许试试this？ BigDecimal 是任意精度。
感谢您的链接。我也用谷歌搜索了一下，发现有很多库可以处理 long double （在这里可能就足够了）甚至是任意精度。如果我能找到解决方案，我会发布我的结果。

【解决方案2】：

你说行列式的当前值大约是10^-300。

您是否要使行列式处于某个值，例如 1？如果是这样，重新调整是尴尬的：您正在考虑的矩阵是ill-conditioned，并且考虑到机器的精度，您应该考虑输出行列式为零。换句话说，不可能得到可靠的逆。

我建议修改矩阵的列或行，而不是重新调整它。

我用 R 用随机矩阵（随机正态值）做了一个小测试，看来行列式显然应该是非零的。

> n=100
> M=matrix(rnorm(n**2),n,n)
> det(M)
[1] -1.977380e+77
> kappa(M)
[1] 2318.188

【讨论】：

@wok：谢谢你的回答。实际上，矩阵可能是病态的，因为对于所研究的问题，我正在寻找该矩阵具有消失行列式的参数，因此至少一个特征值应该消失。你认为在这个意义上看决定因素根本不合理吗？
“消失”是指等于零吗？矩阵有什么特定的属性吗？
@wok：是的，我的意思是 det(M) = 0，在我的例子中是系统的共振条件。矩阵的形式为 M = {A1 B1; A2 B2}。 A是对称的，B不是。真诚的。
不确定我能找到更好的方法。我想知道你是否可以碰巧在你考虑的矩阵中有一个明显非零的行列式。
@wok：对于矩阵的低阶，我得到了正确的值（对应于共振的行列式的绝对值比非共振情况小几个数量级）所以我假设通过增加“分辨率”，我不会得到根本不同的行为。

【解决方案3】：

如果速度不是问题，您可能希望使用det(e^A) = e^(tr A) 并将A 一些缩放常数乘以您的矩阵（这样A - I 的光谱半径小于1）。

编辑：在 MatLab 中，矩阵的对数 (logm) 是通过三角化计算的。因此，您最好计算矩阵的特征值并将它们相乘（或者更好的是，添加它们的对数）。您没有指定您的矩阵是否对称：如果是，则查找特征值比不是更容易。

【讨论】：

非常感谢您的 ansatz。我已经尝试过这个，因为我正在搜索局部最小值，所以使用 ln(detM) = tr(lnM) 及其导数的相应表达式，但正如您指出的那样，这在计算上非常慢。
我是用你提供的身份来解决问题的，所以我决定接受你的回答:)

【解决方案4】：

您应该意识到，当您将矩阵乘以常数 k 时，您会将矩阵的行列式缩放 k^n，其中 n 是矩阵的维数。所以对于 n = 1000 和 k = 2，你将行列式缩放

>> 2^1000
ans =
     1.07150860718627e+301

这当然是一个巨大的数字，所以您可能会认为它应该会失败，因为在双精度中，MATLAB 只会表示与 realmax 一样大的浮点数。

>> realmax
ans =
     1.79769313486232e+308

没有必要重新计算该行列式的所有工作，并不是说计算这样一个巨大矩阵的行列式无论如何都是一个非常适定的问题。

【讨论】：

根据定义，行列式是从 Mn(K) 到 K 的唯一交替 多线性 映射，使得 F(Id)=1。
非常感谢您的回答。所研究矩阵的行列式确实非常低，在 N ~ 200 处达到 1.e-300，这使得重新缩放成为必要。您有什么建议可以克服 MATLAB 的双精度限制，比如使用 Mathematica？
为什么要重新缩放？事实上，为什么要使用行列式本身呢？使用行列式的日志。然后，您永远不需要为了将行列式保持在浮点运算的动态范围内而进行任意重新缩放。行列式的对数很容易从 U 的对角线元素的对数中获得，如矩阵的 LU 分解返回的那样。请注意该对角线上的任何负数，但它们很容易处理。
感谢您的建议。我刚刚尝试了这个 ansatz 但事实证明 tr(logm M) 和 det(M) 产生了应有的结果； NaN 用于巨大的矩阵。我将进一步尝试使用 MATLAB 工具箱或 C/C++ 解决方法进行长双精度范围，并在找到解决方案后立即在此处报告。