如何在Mathematica中高效地设置矩阵的次要元素？

Question

当我看到belisarius对于具有元素均匀分布的非奇异整数矩阵生成的问题时，我正在研究Dana Randal的一篇论文“高效随机非奇异矩阵生成”。所提出的算法是递归的，并涉及生成低维矩阵并将其分配给给定的子矩阵。我使用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列的原始矩阵。

我现在需要创建一个n乘以n的矩阵，在位置{i,j}上放置给定的小矩阵mm。我在算法中使用的是：

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的建议，我研究了通过ArrayFlatten实现ExpandMinor的方法。

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}

- Sasha

---

抱歉，我今天有点懒 : )。Minor 就像这个例子一样，http://en.wikipedia.org/wiki/Minor_(linear_algebra)#Example，是吗？ - Dr. belisarius

@belisarius 那个页面上的次要元素是我所谈论的内容的决定因素。我将编辑我的帖子以做出更详细的解释。 - Sasha

可能相关：https://dev59.com/vm855IYBdhLWcg3wj1IS - Dr. belisarius

你在这里面临的问题是大多数像Insert等这样的Mathematica内置函数是不可变的，它们会创建一个副本。特别是对于较大的矩阵，正是这种复制使得代码效率低下。你唯一的朋友是Part，以向量化的方式用来就地修改矩阵。对于手头的这个案例，我在下面发布了一个解决方案。至于是否可以从中提取出一些通用函数来执行你所要求的一般任务，我还不知道，但似乎是很有可能的。 - Leonid Shifrin

Sasha，恭喜你获得2000的声望值！ - Mr.Wizard

Answer 1

我花了一些时间才来到这里，但由于我在博士后期间花了很多时间生成随机矩阵，所以我忍不住要说一下。代码中的主要低效性来自于需要移动矩阵（复制它们）。如果我们能够重新设计算法，使得我们只在原地修改单个矩阵，那么我们就可以获得很大的优势。为此，我们必须计算插入向量/行将最终出现的位置，考虑到我们通常会在较小的矩阵中间插入并移动元素。这是可能的。以下是代码：

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}

我希望在对上述编译代码进行仔细分析后，可以进一步提高性能。此外，在编译时我没有使用运行时的Listable属性，但这应该是可能的。还有可能代码中执行次要任务的部分是通用的，因此逻辑可以从主函数中分离出来 - 我尚未进行调查。

- Leonid Shifrin

---

@Leonid 注意，ConstantArray 在 Compile 内不受支持并生成对求值器的调用。当用等效的 Table 调用替换时，代码的速度提高了约 5%。我会详细研究代码。请在方便的时候更新代码，并使用我在编辑中添加的 Mod 调用更新定义 NonSignularRanomMatrix 和 NonSignularRanomMatrixAlt。谢谢。 - Sasha

@belisarius 这篇文章讲述了有限域中的生成，但我忘记了最后的模约减。一旦应用它并更新了我的原始代码，矩阵元素确实是均匀分布的。现在有了Leonid的高性能代码，你可以考虑将其用于过滤算法。 - Sasha

@Sasha 是的。下次我重新运行过滤器时，我一定会试试。顺便问一下，你有没有注意到你成功让Leonid写了“嵌套”循环？ :) - Dr. belisarius

@ Sasha 谢谢。我相信可以进行更多的优化（包括你提到的那个） - 我只是发布了第一个正常工作的版本。我很快会用你提到的修改更新代码。 - Leonid Shifrin

@belisarius 很好的发现！一些嵌套循环可以用函数构造来替换，但是昨天我有点懒 :) - Leonid Shifrin

Answer 2

针对你提出的第一个问题（我希望我理解得正确），是否可以将MinorMatrix写成以下形式？

MinorMatrixAlt[m_?MatrixQ, {i_, j_}] := Drop[mat, {i}, {j}]

- 681234