# 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 ? Sep 14, 2014 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, 2014 at 19:47
• @rm-rf For gaussian entries, yes. Otherwise, maybe asymptotically, but for any fixed dimension, the distribution will be different. Sep 14, 2014 at 19:58
• Ah, yes of course.
– rm -rf
Sep 14, 2014 at 20:00
• Have you seen this? mathematica.stackexchange.com/q/7887/219 Dec 25, 2015 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. Sep 14, 2014 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. Sep 14, 2014 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. Sep 14, 2014 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, 2014 at 19:41

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.

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.