一边看着belisarius关于generation of non-singular integer matrices with uniform distribution of its elements的问题,我正在研究 Dana Randal 的一篇论文,“Efficient generation of random non-singular matrices”。所提出的算法是递归的,涉及生成一个较低维度的矩阵并将其分配给给定的次要。我使用了 Insert
的组合和 Transpose
这样做,但必须有更有效的方法来做到这一点。你会怎么做?
以下是代码:
Clear[Gen];
Gen[p_, 1] := {{{1}}, RandomInteger[{1, p - 1}, {1, 1}]};
Gen[p_, n_] := Module[{v, r, aa, tt, afr, am, tm},
While[True,
v = RandomInteger[{0, p - 1}, n];
r = LengthWhile[v, # == 0 &] + 1;
If[r <= n, Break[]]
];
afr = UnitVector[n, r];
{am, tm} = Gen[p, n - 1];
{Insert[
Transpose[
Insert[Transpose[am], RandomInteger[{0, p - 1}, n - 1], r]], afr,
1], Insert[
Transpose[Insert[Transpose[tm], ConstantArray[0, n - 1], r]], v,
r]}
]
NonSingularRandomMatrix[p_?PrimeQ, n_] := Mod[Dot @@ Gen[p, n], p]
它确实生成了一个非奇异矩阵,并且矩阵元素分布均匀,但要求 p 是素数:
代码也不是每一个都有效,我怀疑是因为我的矩阵构造函数效率低下:
In[10]:= Timing[NonSingularRandomMatrix[101, 300];]
Out[10]= {0.421, Null}
编辑 所以让我浓缩一下我的问题。给定矩阵的次矩阵
m
可以计算如下:MinorMatrix[m_?MatrixQ, {i_, j_}] :=
Drop[Transpose[Drop[Transpose[m], {j}]], {i}]
它是带有
i
的原始矩阵第 - 行和 j
-th 列已删除。我现在需要创建一个大小为
n
的矩阵来自 n
将具有给定的次矩阵 mm
在位置 {i,j}
.我在算法中使用的是:ExpandMinor[minmat_, {i_, j_}, v1_,
v2_] /; {Length[v1] - 1, Length[v2]} == Dimensions[minmat] :=
Insert[Transpose[Insert[Transpose[minmat], v2, j]], v1, i]
例子:
In[31]:= ExpandMinor[
IdentityMatrix[4], {2, 3}, {1, 2, 3, 4, 5}, {2, 3, 4, 4}]
Out[31]= {{1, 0, 2, 0, 0}, {1, 2, 3, 4, 5}, {0, 1, 3, 0, 0}, {0, 0, 4,
1, 0}, {0, 0, 4, 0, 1}}
我希望这可以更有效地完成,这就是我在问题中征求的意见。
根据 blisarius 的建议,我考虑实现
ExpandMinor
通过 ArrayFlatten
.Clear[ExpandMinorAlt];
ExpandMinorAlt[m_, {i_ /; i > 1, j_}, v1_,
v2_] /; {Length[v1] - 1, Length[v2]} == Dimensions[m] :=
ArrayFlatten[{
{Part[m, ;; i - 1, ;; j - 1], Transpose@{v2[[;; i - 1]]},
Part[m, ;; i - 1, j ;;]},
{{v1[[;; j - 1]]}, {{v1[[j]]}}, {v1[[j + 1 ;;]]}},
{Part[m, i ;;, ;; j - 1], Transpose@{v2[[i ;;]]}, Part[m, i ;;, j ;;]}
}]
ExpandMinorAlt[m_, {1, j_}, v1_,
v2_] /; {Length[v1] - 1, Length[v2]} == Dimensions[m] :=
ArrayFlatten[{
{{v1[[;; j - 1]]}, {{v1[[j]]}}, {v1[[j + 1 ;;]]}},
{Part[m, All, ;; j - 1], Transpose@{v2}, Part[m, All, j ;;]}
}]
In[192]:= dim = 5;
mm = RandomInteger[{-5, 5}, {dim, dim}];
v1 = RandomInteger[{-5, 5}, dim + 1];
v2 = RandomInteger[{-5, 5}, dim];
In[196]:=
Table[ExpandMinor[mm, {i, j}, v1, v2] ==
ExpandMinorAlt[mm, {i, j}, v1, v2], {i, dim}, {j, dim}] //
Flatten // DeleteDuplicates
Out[196]= {True}
最佳答案
我花了一段时间才到这里,但由于我花了我博士后的大部分时间来生成随机矩阵,我无法帮助它,所以就这样吧。代码中的主要低效率来自移动矩阵(复制它们)的必要性。如果我们可以重新制定算法,以便我们只修改一个矩阵,我们就可以大获全胜。为此,我们必须计算插入的向量/行结束的位置,因为我们通常会插入较小矩阵的中间并因此移动元素。这个有可能。这是代码:
gen = Compile[{{p, _Integer}, {n, _Integer}},
Module[{vmat = Table[0, {n}, {n}],
rs = Table[0, {n}],(* A vector of r-s*)
amatr = Table[0, {n}, {n}],
tmatr = Table[0, {n}, {n}],
i = 1,
v = Table[0, {n}],
r = n + 1,
rsc = Table[0, {n}], (* recomputed r-s *)
matstarts = Table[0, {n}], (* Horizontal positions of submatrix starts at a given step *)
remainingShifts = Table[0, {n}]
(*
** shifts that will be performed after a given row/vector insertion,
** and can affect the real positions where the elements will end up
*)
},
(*
** Compute the r-s and vectors v all at once. Pad smaller
** vectors v with zeros to fill a rectangular matrix
*)
For[i = 1, i <= n, i++,
While[True,
v = RandomInteger[{0, p - 1}, i];
For[r = 1, r <= i && v[[r]] == 0, r++];
If[r <= i,
vmat[[i]] = PadRight[v, n];
rs[[i]] = r;
Break[]]
]];
(*
** We must recompute the actual r-s, since the elements will
** move due to subsequent column insertions.
** The code below repeatedly adds shifts to the
** r-s on the left, resulting from insertions on the right.
** For example, if vector of r-s
** is {1,2,1,3}, it will become {1,2,1,3}->{2,3,1,3}->{2,4,1,3},
** and the end result shows where
** in the actual matrix the columns (and also rows for the case of
** tmatr) will be inserted
*)
rsc = rs;
For[i = 2, i <= n, i++,
remainingShifts = Take[rsc, i - 1];
For[r = 1, r <= i - 1, r++,
If[remainingShifts[[r]] == rsc[[i]],
Break[]
]
];
If[ r <= n,
rsc[[;; i - 1]] += UnitStep[rsc[[;; i - 1]] - rsc[[i]]]
]
];
(*
** Compute the starting left positions of sub-
** matrices at each step (1x1,2x2,etc)
*)
matstarts = FoldList[Min, First@rsc, Rest@rsc];
(* Initialize matrices - this replaces the recursion base *)
amatr[[n, rsc[[1]]]] = 1;
tmatr[[rsc[[1]], rsc[[1]]]] = RandomInteger[{1, p - 1}];
(* Repeatedly perform insertions - this replaces recursion *)
For[i = 2, i <= n, i++,
amatr[[n - i + 2 ;; n, rsc[[i]]]] = RandomInteger[{0, p - 1}, i - 1];
amatr[[n - i + 1, rsc[[i]]]] = 1;
tmatr[[n - i + 2 ;; n, rsc[[i]]]] = Table[0, {i - 1}];
tmatr[[rsc[[i]],
Fold[# + 1 - Unitize[# - #2] &,
matstarts[[i]] + Range[0, i - 1], Sort[Drop[rsc, i]]]]] =
vmat[[i, 1 ;; i]];
];
{amatr, tmatr}
],
{{FoldList[__], _Integer, 1}}, CompilationTarget -> "C"];
NonSignularRanomMatrix[p_?PrimeQ, n_] := Mod[Dot @@ Gen[p, n],p];
NonSignularRanomMatrixAlt[p_?PrimeQ, n_] := Mod[Dot @@ gen[p, n],p];
这是大矩阵的时间安排:
In[1114]:= gen [101, 300]; // Timing
Out[1114]= {0.078, Null}
对于直方图,我得到了相同的图,并且效率提高了 10 倍:
In[1118]:=
Histogram[Table[NonSignularRanomMatrix[11, 5][[2, 3]], {10^4}]]; // Timing
Out[1118]= {7.75, Null}
In[1119]:=
Histogram[Table[NonSignularRanomMatrixAlt[11, 5][[2, 3]], {10^4}]]; // Timing
Out[1119]= {0.687, Null}
我希望通过仔分割析上述编译代码,可以进一步提高性能。另外,我没有在 Compile 中使用运行时 Listable 属性,虽然这应该是可能的。也可能是执行分配给未成年人的代码部分足够通用,因此可以将逻辑从主要功能中分解出来 - 我还没有调查。
关于wolfram-mathematica - 如何在 Mathematica 中有效地设置矩阵的次要?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/5771217/