作为一种合成噪声生成算法的一部分，我必须即时构建大量的非奇异正方形矩阵。其中 ai,j (i,j:1..n)，对于 ∀ (i,j)，ai,j∈ ℤ 且 0 ≤ ai,j ≤ k 且 Det[a] ≠ 0，但 ai,j 也应该是随机的，遵循 [0,k] 的均匀分布。在当前版本中，问题中 n ≅ 300，k≅ 100。在 Mathematica 中，我可以非常快速地生成随机元素矩阵，但问题在于我还必须检查奇异性。我目前使用行列式值进行检查。问题是对于 300x300 矩阵，这个检查需要近2秒钟，而我无法承受。当然，我可以通过选择随机的第一行并构造连续的正交行来构造行，但我不知道如何保证这些行的元素遵循 [0,k] 的均匀分布。我正在寻找 Mathematica 中的解决方案，但只要能更快地生成矩阵的算法也可以。NB> U[0,k] 条件意味着对于一组矩阵，每个位置 (i,j) 都应该遵循均匀分布。
4个回答
5
在Matlab和Octave中，一个500x500的矩阵的行列式和LU分解基本上是瞬间完成的。 Mathematica有一种模式可以调用LAPACK或类似的库吗？您可能需要注明数组应被视为浮点数而不是符号；这可能会使它更快。作为比较，对于一个5000x5000的矩阵，在我的系统上使用Octave进行LU分解需要8.66秒；500x500应该比那快1000倍左右。
- Jeremiah Willcock
14
2@belisarius：如果矩阵包含实数，那么可能没有必要检查其是否奇异（或近似奇异）。生成奇异矩阵的概率基本上为0（在“实际”实数中确切为0）。 http://mss3.libraries.rutgers.edu/dlr/TMP/rutgers-lib_25959-PDF-1.pdf（使用Google缓存）中的论文指出，所有元素均为-1或1的大小为500x500的矩阵是奇异的概率非常非常小。 - Jeremiah Willcock
@Jeremiah 谢谢你的回答。请注意，由于矩阵元素是整数，奇异性的概率肯定不为零。此外，维度只是当前的实现，我将尝试在未来使用较小的数字。我想我可以调用LAPACK，但从未这样做过。我会尝试研究一下它。 - Dr. belisarius
@belisarius：那么http://techreports.lib.berkeley.edu/accessPages/CSD-91-658.html呢？他们声称可以高效地完成这项任务；请注意，您只需要选择有限域作为Z_{the next prime above k}，我相信在该领域中任何非奇异的东西在有理数中也是非奇异的。 - Jeremiah Willcock
如果你只选择GF(2)，那么在GF(2)中非奇异的元素也是有理数中的非奇异元素。简单来说，你实际上是在检查 det(A) mod 2。如果 det(A) mod 2 = 1，那么 det(A) 是一个奇整数，因此 A 是非奇异的。如果 det(A) mod 2 = 0，那么 det(A) 是一个偶整数，但我们只关心其中的零。det(A) = 2 或 4 或... 都没问题，但我们会通过测试 mod 2 把它们丢弃。这对于任意 N 都成立，但 N 越大，你将丢弃的错误结果越少。也许检查 mod 2 的速度足够快且足够好？ - President James K. Polk
1@belisarius：请参阅http://www.inf.ethz.ch/personal/emo/PublFiles/RankSparseRandMatr_RSA10_97.pdf -- 该研究声称，使用均匀分布的特定变体生成有限域上的矩阵时得到奇异矩阵的概率很小，因此即使您假设所有有限域上的奇异矩阵都是误报，也不会对您的分布产生太大影响。此外，请查看该论文介绍中的引文和脚注以获取更多相关事实。 - Jeremiah Willcock
显示剩余9条评论
5
您可以使用
MatrixRank 替代。在我的机器上，对于大型的nxn整数矩阵，它的速度大约是 n/10 倍快。- Simon
4
在我的电脑上，n/10也是如此。这几乎将我的时间推到了可承受的极限。谢谢！ - Dr. belisarius
1我的机器上没有获得如此戏剧性的加速 - 对于
500x500 矩阵，MatrixRank 大约比 Det 快10倍。如果我们使用 MatrixRank[N@yourMatrix]，这似乎仍然是可靠的，并且可以将速度提高约两倍。人们还可以考虑使用 CUDA - 至少对于某些矩阵操作（例如矩阵乘法），速度提升可能非常显着，特别是对于像CUDA Tesla C1060或更好的强大卡。此外，如果有多个核心可用，则可以使用 ParallelMap[MatrixRank, {m1,m2,...,mn}]，其中 n 是核心数，以换取内存来提高速度。 - Leonid Shifrin@Leonid 不要害羞，把你的深刻评论发布为答案，这样未来的读者就不会错过它们。评论线程中的所有内容都应该是相关的。关于你的N@建议，问题在于将奇异矩阵投入集合中会破坏整个集合（将每个集合视为某些子空间的基础），然后是集合间的排名，因为从此时起我只使用无法检测到病态元素的数值算法。我正在进行原型设计，如果这个方法行得通，那么并行处理技巧将用于实际的扩展问题。 - Dr. belisarius
@belisarius 好的，我明白了。最初它只是 @Simon答案的扩展（
N建议），然后我没有切换模式。并行部分确实起作用-在我的6核机器上给我提供了6个矩阵的0.07秒，这比非并行版本快了正好6倍。如果你有一个多核机器，这样做非常轻松，因此在原型设计时可以尝试一下。考虑到你的情况，使用CUDA可能不适合你，除非你在CUDA上重新实现整数矩阵的Det或MatrixRank函数。 - Leonid Shifrin3
如果您在奇异性测试中使用数值近似矩阵，您将获得更快的速度。
k = 100; n = 500;
mat = RandomInteger[100, {n, n}];
AbsoluteTiming[Det[mat] == 0]
Out[57]= {6.8640000, False}
AbsoluteTiming[Det[N@mat] == 0.] (*warning light!!*)
Out[58]= {0.0230000, False}
AbsoluteTiming[MatrixRank[N@mat] != n]
输出[59] = {0.1710000，False}
不幸的是，最快的测试并不可靠。但是排名测试应该很好用。这里有一个快速示例，在其中我们将最后一行替换为前面行的总和。
mat2 = mat;
mat2[[-1]] = Total[Most[mat]];
AbsoluteTiming[Det[mat2] == 0]
Out[70]= {9.4750000, 真}
AbsoluteTiming[Det[N@mat2] == 0.]
Out[69]= {0.0470000, False}
AbsoluteTiming[MatrixRank[N@mat2] != n]
输出[71] = {0.1440000，真}
原则上，我认为排名测试可能会出现假阴性的小概率情况，比如由于病态条件。由于您的使用将更好地容忍假阳性（即错误的奇异声明），因此您可以改为在素模数上测试奇异性。我认为这是其他人提出的建议之一。
继续上面的例子：
AbsoluteTiming[Det[mat, Modulus -> Prime[1000]]]
输出结果[77]为{0.6320000, 4054}
AbsoluteTiming[Det[mat2, Modulus -> Prime[1000]]]
输出结果[78]= {0.6470000, 0}。
虽然速度慢，但比在有理数上工作要快。就其价值而言，对于大多数用途，我相当有信心通过MatrixRank[N[matrix]]进行更快速的测试得出结果。
Daniel Lichtblau 沃尔夫勒姆研究公司
- Daniel Lichtblau
7
+1 感谢Daniel。至于“由于您的使用将更容忍假阳性”，我仍然需要证明这不会违反过滤矩阵的均匀分布。看起来是对的，但是... - Dr. belisarius
@belisarius 这不会违反分布的均匀性，取决于随机数生成器的质量（非常好）。这个想法是，在这个分布中，你的非奇异矩阵是等可能的选择，因为所有矩阵都是等可能的，只有在奇异时才会丢弃。类似于在立方体中选择不规则随机元素的方法，通过选择并丢弃不在球内的元素来实现。 - Daniel Lichtblau
@belisarius 对不起，我误解了。您的关注点在于丢弃在模p下是奇异的非奇异矩阵。我必须说，这是一个非常低的概率事件，以至于我怀疑您是否能设计出随机测试来检测任何差异。尽管如此，您可以使用几个不同质数模下的测试进一步验证奇异性。现在失败概率为1 /（您一生中不会发生的事情）^ 3。而且这很便宜，因为您不会在一次这样的测试中失败，因此无需进一步测试。 - Daniel Lichtblau
@Daniel，你上一条评论说得很对。我认为正式的演示并不容易，但我可以进行一些实验来“证明”均匀分布。 - Dr. belisarius
@belisarius 如果您需要测试忠实（仅在非奇异时通过，仅在奇异时失败），那么您可以执行其中一个快速测试modp，在mod p表示为奇异时仅在Q上进行测试。这是Mark McClure提出的方案之一。如果k，n的大小都很大，则不应该看到mod p的故障，因此不需要更昂贵的测试而导致减速。但是您也将知道您的代码是防弹的。 - Daniel Lichtblau
显示剩余2条评论
1
这是我之前留言的一个扩展。我同意Dan的观点，数值版本返回假阳性的可能性非常小。尽管如此，您可以通过检查奇异值并可靠地返回False来避免这种情况，如果最小奇异值大于某个误差容限。（不可否认，找到可证明的容限可能有点棘手。）如果最小奇异值过小，则可以将Det应用于整数矩阵。
这是一个函数，对于大多数非奇异矩阵应该快速返回False。如果矩阵接近奇异，则执行更昂贵的整数计算。
singularQ[M_?MatrixQ] := If[
Last[SingularValueList[N[M], Tolerance -> 0]] > 1/10.0^8,
False, Det[M] == 0];
这里有200个符合您描述的矩阵。其中一个被篡改成奇异矩阵。
SeedRandom[1];
matrices = RandomInteger[{0, 100}, {200, 300, 300}];
matrices[[100, -1]] = Sum[RandomInteger[{0, 10}]*matrices[[100, k]],
{k, 1, Length[matrices[[100]]] - 1}];
现在，让我们一起寻找所有奇异矩阵的索引，并且在这个过程中观察。
Flatten@Monitor[Table[
If[singularQ[matrices[[k]]], k, {}],
{k, 1, Length[matrices]}], k]
- Mark McClure
1
1检查近似数值矩阵的秩本质上是通过找到奇异值并计算大于某个公差乘以最大值的奇异值数量来完成的，其中公差是该函数的一个选项。至少我是这样认为的。 - Daniel Lichtblau
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接
random_matrix(ZZ, n, n, algorithm='subspaces', rank=n)，但对于大的n，这样做非常慢...我相信专家们一定能够向你展示更好的方法。 - Simon