MATLAB 中大矩阵的行列式
从一个模拟问题中,我想在 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 的约束,这最终为我提供了对我的问题有效的解决方法。不幸的是,我不知道这样的解决方法是否可以推广。
from a simulation problem, I want to calculate complex square matrices on the order of 1000x1000 in MATLAB. Since the values refer to those of Bessel functions, the matrices are not at all sparse.
Since I am interested in the change of the determinant with respect to some parameter (the energy of a searched eigenfunction in my case), I overcome the problem at the moment by first searching a rescaling factor for the studied range and then calculate the determinants,
result(k) = det(pre_factor*Matrix{k});
Now this is a very awkward solution and only works for matrix dimensions of, say, maximum 500x500.
Does anybody know a nice solution to the problem? Interfacing to Mathematica might work in principle but I have my doubts concerning feasibility.
Thank you in advance
Robert
Edit: I did not find a convient solution to the calculation problem since this would require changing to a higher precision. Instead, I used that
ln det M = trace ln M
which is, when I derive it with respect to k
A = trace(inv(M(k))*dM/dk)
So I at least had the change of the logarithm of the determinant with respect to k. From the physical background of the problem I could derive constraints on A which in the end gave me a workaround valid for my problem. Unfortunately I do not know if such a workaround could be generalized.
如果你对这篇内容有疑问,欢迎到本站社区发帖提问 参与讨论,获取更多帮助,或者扫码二维码加入 Web 技术交流群。
绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱,如果其他用户或者作者回复了您的评论,将不能在第一时间通知您!
发布评论
评论(4)
您应该意识到,当您将矩阵乘以常数 k 时,您会将矩阵的行列式缩放 k^n,其中 n 是矩阵的维度。因此,对于 n = 1000 和 k = 2,您可以将行列式缩放为
这当然是一个巨大的数字,因此您可能会认为它会失败,因为在双精度中,MATLAB 将仅表示与 realmax 一样大的浮点数。
不需要做所有重新计算该行列式的工作,而且计算这样一个巨大矩阵的行列式无论如何都是一个非常适定的问题。
You should realize that when you multiply a matrix by a constant k, then you scale the determinant of the matrix by k^n, where n is the dimension of the matrix. So for n = 1000, and k = 2, you scale the determinant by
This is of course a huge number, so you might expect that it should fail, since in double precision, MATLAB will only represent floating point numbers as large as realmax.
There is no need to do all the work of recomputing that determinant, not that computing the determinant of a huge matrix like that is a terribly well-posed problem anyway.
如果速度不是问题,您可能需要使用 det(e^A) = e^(tr A) 并将一些缩放常数乘以矩阵作为
A
(这样A - I
的谱半径小于 1)。编辑:在 MatLab 中,矩阵的对数 (
logm
) 是通过三角化计算的。因此,最好计算矩阵的特征值并将它们相乘(或者更好的是添加它们的对数)。您没有指定矩阵是否对称:如果是,则查找特征值比不对称更容易。If speed is not a concern, you may want to use
det(e^A) = e^(tr A)
and take asA
some scaling constant times your matrix (so thatA - I
has spectral radius less than one).EDIT: In MatLab, the log of a matrix (
logm
) is calculated via trigonalization. So it is better for you to compute the eigenvalues of your matrix and multiply them (or better, add their logarithm). You did not specify whether your matrix was symmetric or not: if it is, finding eigenvalues are easier than if it is not.你说行列式的当前值大约是10^-300。
您是否试图获得某个值(例如 1)的行列式?如果是这样,重新缩放会很尴尬:您正在考虑的矩阵是病态的,并且,考虑到机器的精度,您应该将输出行列式视为零。换句话说,不可能得到可靠的逆。
我建议修改矩阵的列或行,而不是重新缩放它。
我用 R 用随机矩阵(随机正态值)进行了一个小测试,看来行列式应该明显非零。
You said the current value of the determinant is about 10^-300.
Are you trying to get the determinant at a certain value, say 1? If so, rescaling is awkward: the matrix you are considering is ill-conditioned, and, considering the precision of the machine, you should consider the output determinant to be zero. It is impossible to get a reliable inverse in other words.
I would suggest to modify the columns or lines of the matrix rather than rescale it.
I used R to make a small test with a random matrix (random normal values), it seems the determinant should be clearly non-zero.
严格来说,这不是一个 matlab 解决方案,但您可能需要考虑使用 Mahout。它专为大规模线性代数而设计。 (1000x1000 对于它所使用的比例来说没有问题。)
您会 调用 java 将数据传递给 Mahout 或从 Mahout 传递数据。
This is not strictly a matlab solution, but you might want to consider using Mahout. It's specifically designed for large-scale linear algebra. (1000x1000 is no problem for the scales it's used to.)
You would call into java to pass data to/from Mahout.