Suppose first that I want to generate a matrix whose elements are independent and identically distributed (i.i.d.) with distribution dist
. This is easy:
randMat[dist_, n_, m_] := RandomVariate[dist, {n, m}]
Now, suppose I wanted to generate a random symmetric matrix whose (non-lower-triangular) elements are i.i.d. Also easy, you say:
randSym[dist_, n_] :=
Module[{uptri =
Table[If[i > j, RandomVariate[dist], 0], {i, 1, n}, {j, 1, n}],
diag = DiagonalMatrix[RandomVariate[dist, n]]},
diag + uptri + Transpose[uptri]]
Now, the problem is that for $n=10000=m,$ the first command takes 2.5 seconds, while the second takes 250(!) seconds. In fact, it takes much longer to generate the matrix than to compute its eigenvalues! The question, then, is: what is the most elegant way to solve the second problem?
(by the way, in MATLAB this is quite efficient because there is a primitive to extract the upper triangular part of the matrix, so something like
uptri(foo) + uptri(foo', 1)
does the right thing)
Answer
ClearAll[randomSymMat];
randomSymMat = Module[{mat = RandomVariate[#, {#2, #2}], upper, diag},
upper = UpperTriangularize[mat, 1];
diag = DiagonalMatrix[Diagonal@mat];
diag + upper + Transpose[upper]] &;
dist = NormalDistribution[3, 5];
First[AbsoluteTiming[res = randomSymMat[dist, 1000]]]
(* 0.065046 *)
res == Transpose[res]
(* True *)
Comments
Post a Comment