如何在Mathematica中有效地設置矩陣的次要？

Question

雖然看着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列對應的原始矩陣。

我現在需要創建一個大小爲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} 

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} 

我希望仔細分析上面的編譯代碼，可以進一步提高性能。另外，我沒有在Compile中使用運行時Listable屬性，但這應該是可能的。也可能是執行賦值給未成年人的代碼部分是足夠通用的，這樣邏輯可以被分解出主函數 - 我還沒有調查過。

Answer 2

對於你的問題的第一部分（我希望我的理解正確），可以 MinorMatrix可以寫成？

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