# Improve random generation of matrices

My problem is too plenty of subtleties that a detailed explaination would be boring. I will try to state my problem in a more general way.

I have a 5x5 matrix called $M$ which depends on several parameters that are randomly generated by steps. I use the command RandomReal[-1,1].

Let's say a subset of these parameters is generated and we fix them (because they may satisfy some experimental constraint). Another subset of parameters is still free and we start generating them.

My aim is to find the subset of parameters that, once $M$ is diagonalized (namely $\hat{M} = U_1 M U_2^\dagger$ is diagonal with two different unitary matrices), allow me to satisfy other constraints which depend on the diagonal entries of $\hat{M}$.

Since the number of the free parameters is huge (around $\sim 4$), Mathematica is not able to find the Eigenvalues of $M$ because this would require a to solve high-order equations on these parameters.

Is there an easy way to solve these kind of problems in Mathematica?

EDIT (this is an example)

a = {x, RandomReal[{-1, 1}], RandomReal[{-1, 1}], RandomReal[{-1, 1}],y};
b = {w, RandomReal[{-1, 1}], z, RandomReal[{-1, 1}], RandomReal[{-1,1}]};
c = {RandomReal[{-1, 1}], RandomReal[{-1, 1}], x, w, RandomReal[{-1,1}]};
d = {w, RandomReal[{-1, 1}], x, RandomReal[{-1, 1}], y};
e = {RandomReal[{-1, 1}], RandomReal[{-1, 1}], x, y, RandomReal[{-1, 1}]};

M = {a, b, c, d, e}


My output is

{{x, -0.431882, -0.959419, 0.957675, y}, {w, -0.108609, z, -0.753215, -0.311028}, {-0.858034, 0.397821, x, w, 0.971876}, {w, -0.978155, x, 0.758728, y}, {0.861983, -0.677996, x, y, 0.776895}}


If I compute the Eigenvalues of $M.M^T$ with

Eigenvalues[M.Transpose[M]]


I get an implicit answer in terms of the command $Root[...]$.

Then, I generate the values of $x,y,w,z$ and the matrices $U_1$,$U_2$

x = RandomReal[{-1, 1}];
y = RandomReal[{-1, 1}];
w = RandomReal[{-1, 1}];
z = RandomReal[{-1, 1}];


and the output of M is

M
{{0.620364, -0.431882, -0.959419, 0.957675, -0.462837}, {-0.045156, -0.108609, 0.605672, -0.753215, -0.311028}, {-0.858034, 0.397821, 0.620364, -0.045156, 0.971876}, {-0.045156, -0.978155, 0.620364, 0.758728, -0.462837}, {0.861983, -0.677996, 0.620364, -0.462837, 0.776895}}


Then, I compute the unitary matrices such that $U_1 M U_2^T$ is diagonal

{val1up, vec1up} = Eigensystem[M.Transpose[M]];
U1 = Flatten[Orthogonalize /@ FindClusters[N@val1up -> vec1up], 1] //
Simplify // Chop;
{val2up, vec2up} = Eigensystem[Transpose[M].M];
U2 = Flatten[Orthogonalize /@ FindClusters[N@val2up -> vec2up], 1];


Now, I can compute the eigenvalues of $\hat{M}$

Sqrt[Eigenvalues[Transpose[M].M]]
{2.20432, 1.6556, 1.27769, 1.09486, 0.0393237}


I want to find those values of $x,y,w,z$ such that the entries $\hat{M}_{11}$ and $\hat{M}_{22}$ are $3$ and $4$ respectively.

• RandomReal[{-1,1}] generates a single random number; how many are the simulated parameters? also, it would help if you could provide some reproducible code for others to work with; are the $U$ matrices composed of eigenvectors? – user42582 Jan 7 '18 at 20:02
• I use RandomReal for each parameter. I'll try to write a simplified code – apt45 Jan 7 '18 at 20:23
• @user42582 I edited my question. Thanks – apt45 Jan 7 '18 at 20:48
• What do you try to model? Is there a particular probability distribution that you try to realize? Have you heard of the Bertrand paradox? – Henrik Schumacher Jan 7 '18 at 23:08

Not sure I have this understood correctly but possibly you want to equate coefficients of characteristic polynomials and solve for free parameters? One form of this char poly is just the expansion of (t-30*(t-4)*three factors with unknown roots and the other is the char poly of m.Transpose[m].

a = {x, RandomReal[{-1, 1}], RandomReal[{-1, 1}], RandomReal[{-1, 1}],
y};
b = {w, RandomReal[{-1, 1}], z, RandomReal[{-1, 1}],
RandomReal[{-1, 1}]};
c = {RandomReal[{-1, 1}], RandomReal[{-1, 1}], x, w,
RandomReal[{-1, 1}]};
d = {w, RandomReal[{-1, 1}], x, RandomReal[{-1, 1}], y};
e = {RandomReal[{-1, 1}], RandomReal[{-1, 1}], x, y,
RandomReal[{-1, 1}]};
mat = {a, b, c, d, e};

clist = Module[{ct, cpol1, cpol2},
cpol1 =
CoefficientList[CharacteristicPolynomial[mat.Transpose[mat], t],
t];
cpol2 = (t - 4)*(t - 3)*(t - v3)*(t - v4)*(t - v5);
Most[cpol1/Last[cpol1] - CoefficientList[cpol2, t]]
];


Equating gives five polynomials in seven unknowns. Could make it a square system by giving values for two of the unknown eigenvalues, or two of the matrix parameters. The former gives a challenging system to NSolve but the latter at least is tractable.

solns = NSolve[clist /. {w -> -.8, x -> .33}];
realsolns = Select[{y, z, v3, v4, v5} /. solns, FreeQ[#, Complex] &];
Union[realsolns[[All, 1 ;; 2]]]

(* {{-1.75838653842, 0.366038707986}, {-1.75838653842,
0.366038707986}, {-1.75838653842, 0.366038707986}, {-1.75838653842,
0.366038707986}, {-1.75838653842, 0.366038707986}, {-1.75838653842,
0.366038707986}, {-1.48805372054, 0.96420053861}, {-1.48805372054,
0.96420053861}, {-1.48805372054, 0.96420053861}, {-1.48805372054,
0.96420053861}, {-1.48805372054, 0.96420053861}, {-1.48805372054,
0.96420053861}} *)

• Thank you! How can I extend the code for a 10x10 matrix where I am still interested in solving in order to get 3 and 4 as eigenvalues? In this case I have $v3,v4,v5,v6...v10$ and mathematica is not able to NSolve this system of equations – apt45 Jan 13 '18 at 8:38
• If it is out of reach for NSolve, maybe FindRoot could be used to find a particular solution? – Daniel Lichtblau Jan 13 '18 at 15:23
• It does not work. In the question I gave a very simplified version of my problem. I have a 10x10 matrix $M$ which depends on several parameters, $M = A+\epsilon B$ where $\epsilon\sim 0.25$. The entries of $A$ are just numbers and it is always diagonalizable.The matrix $B$ instead depends on the parameters and it is symmetric.Basically, I want to see for which values of the parameters the first two smallest eigenvalues of $M.M^T$ belong to a given interval $I$. Can you suggest me a method to deal with such kind of problems? – apt45 Jan 13 '18 at 18:04
• I think the problem is non-trivial and I don't pretend you have to answer – apt45 Jan 13 '18 at 18:06
• It is probably difficult. I will suggest creating a separate question for the new full-size problem, and hope for the best. – Daniel Lichtblau Jan 13 '18 at 18:26