# Generating random symmetric matrix

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)

• UpperTriangularize ? – ybeltukov Sep 14 '14 at 16:15
• Why not just do $M = (A + A^T)/2$? I understand that while the elements are iid, the diagonals and off-diagonals have different variances, but if you're primarily concerned with the eigenvalue statistics, this should be a faster way of generating the matrices. The eigenvalues of this and your matrix should both be Wigner distributed. – rm -rf Sep 14 '14 at 19:47
• @rm-rf For gaussian entries, yes. Otherwise, maybe asymptotically, but for any fixed dimension, the distribution will be different. – Igor Rivin Sep 14 '14 at 19:58
• Ah, yes of course. – rm -rf Sep 14 '14 at 20:00
• Have you seen this? mathematica.stackexchange.com/q/7887/219 – faleichik Dec 25 '15 at 17:53

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 *)

• Duh! For some reason I thought UpperTriangularize[] computed an upper triangular similar matrix. – Igor Rivin Sep 14 '14 at 16:53
• On the other hand, while this clearly solves the egregious problem, this is still over a factor of two slower than the non-symmetric case. Sad. – Igor Rivin Sep 14 '14 at 18:28
• @kguler: shouldn't that be upper+Transpose[upper]-diag since the upper contains the diagonal already? With n=10000, I get 3.2 s for the n-by-n matrix, and 5.1 s for the symmetric one. Not too bad a ratio. – Wouter Sep 14 '14 at 19:33
• @Wouter, upper does not contain the main diagonal; second argument of UpperTriangularize shifts the triangle part off the main digonal. – kglr Sep 14 '14 at 19:41

A slightly simpler code: The presence of the 1 in UpperTriangularize[] grabs the strictly upper triangular components. This way you don't need to define the diag portion for the sum.

    (* Generate symmetric n × m matrix with
component values as random reals in (0,K) *)
RM = RandomReal[K, {n, m}];
SM = UpperTriangularize[RM] + Transpose[UpperTriangularize[RM, 1]]


Of course you can change the RandomReal[] to whatever method you want to generate your entries.

If you have Mathematica 10.3, this is doable by the following code:

randSym[dist_, n_] := Module[{diag, up},
diag = RandomVariate[dist, n];
up = RandomVariate[dist, Binomial[n, 2]];
StatisticsLibraryVectorToSymmetricMatrix[up, diag, n]
]


This method does not sample unnecessary entries hence it is faster. It pastes the sampled values properly into a symmetric matrix.