Generating invertible matrix with lines within a given set

Question

Consider the set options given as below

n=2;m=3;
options = 
 Complement[
  DeleteDuplicates[
   Table[Table[
     Product[Tuples[
          Table[Tuples[{1, -1}, n], {m}]][[i]][[k]][[Tuples[
           Table[i, {i, 1, n}], m][[j]][[k]]]], {k, 1, m}], {j, 1, 
      Length[Tuples[Table[i, {i, 1, n}], m]]}], {i, 1, 
     Length[Tuples[Table[Tuples[{1, -1}, n], {m}]]]}]], 
  Union[Tuples[{1}, n^m], Tuples[{-1}, n^m]]]

I am looking for a way to generate a random $n^m\times n^m$($=8\times8$) invertible matrix with the vector $(1,...,1)$ as the first line and with vectors from options as the other lines.

With

prematarb = RandomChoice[options, n^m - 1];
matarb = Prepend[prematarb, Table[1, {n^m}]];
MatrixForm[matarb]
Det[matarb]

we can generate a matrix like I want and calculate its determinant. But to find one with non-null determinant, I had to run the code about 50 times. Is there a way to make the code automatically repeat until it finds one invertible matrix?

One observation is that the matrices generated through the code above, in general, have repeated lines. Is there a way to minimize these repetitions? Could I change the SeedRandom to achieve the aim in a faster way? I see that in the case $n=m=3$, it is far too hard to find the matrix I want generating randomly. The case above has 14=Length[options] options to fill each line of a $7\times 8$ matrix. But in the case $n=m=3$, we have 126 options to fill each line of a $26\times27$ matrix, and still it is common to see repeated lines within the generated matrices. So if there is a way to make things faster (like modifying SeedRandom or using some way for don't taking lines that already was taken to fill the lines), it could be really handy.

EDIT: We can reduce the set options to the half of the size by removing the symmetric elements, this way:

options2=Table[Reverse[Sort[options]][[i]], {i, 1, Length[options]/2}]

Answer 1

Using your method, here's a way to repeatedly generate a random matrix until the matrix is invertible. First, turn your matrix generation into a function so it's easy to call from two places.

Edit: Remove generation of options from genmat since options doesn't change unless n or m are changed. ConstantArray is faster than Table for making a list of 1s. These changes find solution matrices more quickly (see comments).

genmat[n_, m_] := Prepend[RandomChoice[options,n^m - 1],ConstantArray[1,n^m]]

Next, initialize options. Do this once for new values of n or m.

n = 2; m = 3;
options = 
  Complement[
   DeleteDuplicates[
    Table[Table[
      Product[Tuples[Table[Tuples[{1, -1}, n], {m}]][[i]][[k]][[
         Tuples[Table[i, {i, 1, n}], m][[j]][[k]]]], {k, 1, m}], {j, 
       1, Length[Tuples[Table[i, {i, 1, n}], m]]}],
       {i, 1, Length[Tuples[Table[Tuples[{1, -1}, n], {m}]]]}]], 
   Union[Tuples[{1}, n^m], Tuples[{-1}, n^m]]];

The following While command will run until the determinant of matarb is not zero. Run again to find another matrix.

matarb = genmat[n, m]; (*initialize matrix*)
i = 0; (*counter*)
While[
 Det[matarb] == 0, {matarb = genmat[n, m], i += 1}];
TableForm[{MatrixForm[matarb], Det[matarb], i}, 
 TableDirections -> Row]

The result is a matrix, its determinant, and the number of matrices tested before finding a solution.

Eliminate singular matrices

Finding solution matrices is quicker if there are fewer singular cases to test. As Filburt says in his question, the options matrix is symmetrical in this sense: the top half is the reverse of the negative of the bottom half. The symmetry is easy to see for the n=2, m=3 case:

choices = Length[options]/2 (*half the number of rows in the options matrix*);

TableForm[{MatrixPlot[Take[options, choices]],
  MatrixPlot[Reverse[-Take[options, -choices]]]}, 
 TableHeadings -> {{"Top half", 
    "Negative of bottom half,\nReversed"}}, TableDirections -> Row]

There are 17,297,280 possible test matrices using the original method where 14 rows of options are randomly selected 7 rows at a time. More than 99% of them are singular.

To test symmetry for any n and m, simply:

Take[options,choices] == Reverse[-Take[options,-choices]]

True

- Case where n = 2, m = 3

Because the bottom half of options is the negative of the top half, every corresponding row from the top half and the reverse of the bottom half of options results in a singular matrix. For example, if a test run selects rows 1, 3, 5, and 7 from the top half, only rows 2, 4, 6 from the bottom half will not result in a singular matrix.

Because of the symmetry in the options matrix, only the top half is needed.

cols = n^m;
top = Take[options, choices];

It's unnecessary to compute the bottom half of options, because options == Join[top, Reverse[-top]].

We can entirely eliminate singular test matrices by randomly choosing rows from top and the complement of the rows from -top. Note that top and -top are not singular.

Here's how to make a list of random rows to be selected from top:

tr = RandomSample[Range[1, choices], RandomInteger[choices]]

{6, 1, 5, 4}

Given tr rows, these are the rows to select from -top:

Complement[Range[choices], tr]

{2, 3, 7}

Then every random matarb matrix is not singular (all singular test cases are eliminated):

tr = RandomSample[Range[1,choices], RandomInteger[choices]];
Det[matarb = Prepend[
  Join[top[[tr, All]], -top[[Complement[Range[choices],tr]]]],
    ConstantArray[1, cols]]]

The number of non-singular matrices (13,700 for n=2, m=3) is:

$$ \sum _{k=0}^{n^m-1} p\left(n^m-1,k\right) $$

where $p({n},{k})\text{=}\frac{n!}{(n-k)!}$.

- Case where n = 2, m = any

Everything about the n=2, m=3 case applies when m is any value. Every randomly chosen matrix is not singular.

- Case where n is not 2, m = any

When n is not 2, the top and bottom halves of options are symmetrical in the same way as the n=2 case. Test this as before with:

Take[options, choices] == Reverse[-Take[options, -choices]]

However, although choosing complementary rows reduces the number of singular test matrices, Monte Carlo testing is needed to eliminate the ones that remain. I haven't coded a way to exploit the symmetry in the options matrix when n is not 2, but the same complementary-rows technique will help. Meanwhile, this method using While will find solutions.

choices = Length[options]/2;
cols = n^m;
matarb = Prepend[options[[RandomSample[Range[1,choices],cols-1], All]], 
  ConstantArray[1, cols]];
While[Det[matarb] == 0, 
  matarb = Prepend[
    options[[RandomSample[Range[1,choices],cols-1], All]], 
    ConstantArray[1,cols]]];
Det[matarb]

Answer 2

Filburt's question in a 2-parter... The first part was how to select rows from options without selecting the same row twice for each matarb. That's done by selecting rows using RandomSample and by not selecting corresponding rows from the top and bottom halves of options.

The second part is how to improve the code. This replacement code for computing the options matrix doesn't change the results, but it uses half the memory and it's significantly faster for larger values of n and m.

n = 2; m = 3;
cols = n^m;
(* n == 1 and m == 1 are special cases *)
ii = If[n == 1 || m == 1, 2^(n*m), 2^(n*m - 1)];
tbla = Tuples[ 
   Table[Tuples[{1, -1}, n], m]];
tblb = Tuples[
   Table[i, {i, 1, n}], m];
options = Complement[
   Table[Table[
     Product[tbla[[i]][[k]][[tblb[[j]][[k]]]], {k, 1, m}],
     {j, 1, cols}], {i, 1, ii}],
   Union[{ConstantArray[-1, cols]}, {ConstantArray[1, cols]}]];

Comments about changes:

• For readability, compute cols = n^m. The value is used several times.

• Replace Length with computed values.

It's clear that:

Length[Tuples[Table[i, {i, 1, n}], m]]}] == cols

for all n and m, so we can use {j, 1, cols}. Also:

Length[Tuples[Table[Tuples[{1, -1}, n], {m}]]] == 2^(n*m)

for all n and m, so we can use {i, 1, ii}, where ii is computed with If[m == 1 || n == 1, 2^(n*m), 2^(n*m - 1)].

The 2^(n*m) term grows quickly, so reducing it is useful. Notice that ii is 2^(n*m - 1) for most values of n and m. This change uses half the memory, and makes the code run almost twice as fast. The 2^(n*m) value produces duplicate rows that are removed by Complement, so the result is unchanged.

• Tuples[Table[Tuples[{1, -1}, n], m]] doesn't change inside Product. Compute its value once as tbla. Do the same for tblb = Tuples[Table[i, {i, 1, n}], m]. For the n=3, m=4 case, computing options runs more than 50 times faster.

• Union[{ConstantArray[-1, cols]}, {ConstantArray[1, cols]}] is a bit faster than Union[Tuples[{1}, n^m], Tuples[{-1}, n^m]].

With these changes, the code is faster and uses less memory.