Finding matrices , given

Finding the matrices & finding the smallest probability for interval

Initialization--Setup

A0={{1,0,0,0},{1,1,1,0},{1,0,1,1},{0,0,1,0}}
A1={{1, 1, 0, 0},{0, 1, 1, 1},{0, 1, 0, 2},{0, 0, 0, 1}}
AllMatrix={A0, A1};

U={{1, 0,1,0},{0,1,0,1}};

level=13;

Length[U]

U

MatrixForm[{{1,0,1,0},{0,1,0,1}}]

U

Clear[a]

a/.LinearOptimization[-Total[a],{a.U\[VectorLessEqual]{11,12,11,10}, a\[VectorGreaterEqual]ConstantArray[0,Length[U]]}, {a}]

Functions

vecDivide: Input: x(\in \R^{Length[U]}) Given a matrix U the output says whats the argmax (Sum a_i), such that Sum a_iU_i\leq x
indexing(x, y) gives back the y-th element of the list x; in particular: indexing[AllMatrix, #] & /@ {1, 3, 2, 1} gives the matrices from AllMatrix of index 1, 3, 2, 1 in this order.
Productifier: Given a set of indices, productifier gives back the product of matrices (from AllMatrix) wiht the given indices

vecDivide[x_]:=a/.LinearOptimization[-Total[a],{a.U\[VectorLessEqual]x, a\[VectorGreaterEqual]ConstantArray[0,Length[U]]}, {a}]
h[x_]:=U.x
indexing[x_, y_]:=x[[y]]
productifier[x_]:=Dot@@@indexing[AllMatrix,#]&/@x

MaxEigo: Finds the max modulus eigenvalue
smallestSpec: finds the worst spectral radius from the list x

MaxEigo[x_] := Max[Abs /@ Eigenvalues[x]]
smallestSpec[x_] := Min[MaxEigo /@x]

Init-- calculations

nProd: calulates all the products of level given in level;
allIndex: a list of all indices

Uanized: Gives back the list of ROWS of the matrices we get by left multiplying all matrices of nProd with U
Divided: ~ Gives back the rows of matrices , where  is the best matrix we can get in our theorem. --- Gives back the , such that Sum a_iU_i\leq the rows of Uanized.
AMatrices: The above but shaped correctly into matrices,
worstcase: smallest rowsums of the matrices A.
indices: the index of matrices corresponding to the worst rowsums
WorstMatrixCombos: all list of indices where the product corresponds to the worst row sums.

nProd=Dot@@@Tuples[AllMatrix,level ];
allIndex=Tuples[{1, 2},level ];

nProd

Uanized=Flatten[h/@nProd, 1];

Uanized

Divided=vecDivide/@Uanized;

AMatrices=MatrixForm/@Partition[Divided,2];

worstcase=Min[Total/@Divided];
indices=Quotient[#,3]&/@Position[Total/@Divided, worstcase];
WorstMatrixCombos=Flatten[indexing[allIndex, #]&/@indices, 1];

AMatrices

SmallestSpectralRadius=smallestSpec[nProd]

(1/2) (123+55 (Sqrt[5]))

N[(1/2) (123+55 (Sqrt[5]))]

Output

probability=worstcase^{-1/level}

N[probability]

N[SmallestSpectralRadius^{-1/level}]