# Generating invertible matrix with lines within a given set

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}]

• If you are choosing values from a possible list and you don't want repeats you should look at RandomSample[...]. Also your code is very difficult to follow with all the nested functions, you might consider refactoring it at least for the sake of getting more feedback. – N.J.Evans Oct 11 '17 at 13:01
• @N.J.Evans Thank you for the tip. RandomSample answer the second question. I put the full set options defined in the OP just to work exactly with the set I am interested in. But it could be another set o vectors from $\Bbb R^n$. I only want to find a way to repeat the action of generating a random matrix until it finds an invertible one. – Filburt Oct 11 '17 at 13:11
• The method options2=Table[Reverse[Sort[options]][[i]],{i,1,Length[options]/2}] (which is identical to options2=Take[-options,Length[options]/2]) eliminates half the valid solutions. See my new answer. I hope I clearly explained the symmetry in options. – creidhne Oct 17 '17 at 5:38

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]

• beautiful answer :-) – Filburt Oct 11 '17 at 22:17
• Hey man, could you tell me what is the role of i in your code? What's the difference between your While and While[ Det[matarb] == 0, matarb = genmat[n, m]]; ? It is only to count the number of times the body of While was tested, right? – Filburt Oct 13 '17 at 15:28
• @Filburt Yes, i just a counter. Other than that, it's not needed. If you're looking for only a few of the solution matrices, this Monte Carlo approach is good enough. For more speed, compute the 14 rows just once, outside of the genmat function. It's not necessary to compute them every time genmat runs. If you need every solution matrix (I think there are 13,700 of them), there are better methods. Is this set of 14 (or 2(n^m -1)) rows the only ones you need to test? – creidhne Oct 14 '17 at 7:24
• Yes, once you fix $m$ and $n$, the set options is fixed too and it is the only being tested. Other thing, the same goes for the parentheses, right? I mean, they are not necessary since I use ; to separate the tests withing the body. – Filburt Oct 14 '17 at 12:38
• Sadly, there is no way to give another +1 :-) – Filburt Oct 17 '17 at 11:14

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]}]];


• 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.

• Thank you for the improvements. Happy to see you're interested. There are good knowledge in your answer, thanks for that. I am not a programmer and it's just my first month learning Mathematica. I am happy to see what the forums provided for me. About the algorithm, after all the improvements, the code for generating an arbitrary matrix with those constraints is good enough for me. Now I am looking for a way to generate all these matrices (for m=n=3), and I realize that for my purposes there is no need for taking options completely, only one half (eliminating symmetrics) is sufficient. – Filburt Oct 28 '17 at 14:58
• But this is for another question in the future. – Filburt Oct 28 '17 at 15:00
• @Filburt One more embarrassing fix for the second table -- I didn't see it last night (sigh) ... glad I could help. – creidhne Oct 28 '17 at 17:29