9
$\begingroup$

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)

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ UpperTriangularize ? $\endgroup$
    – ybeltukov
    Sep 14, 2014 at 16:15
  • 1
    $\begingroup$ 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. $\endgroup$
    – rm -rf ♦
    Sep 14, 2014 at 19:47
  • 1
    $\begingroup$ @rm-rf For gaussian entries, yes. Otherwise, maybe asymptotically, but for any fixed dimension, the distribution will be different. $\endgroup$
    – Igor Rivin
    Sep 14, 2014 at 19:58
  • 1
    $\begingroup$ Ah, yes of course. $\endgroup$
    – rm -rf ♦
    Sep 14, 2014 at 20:00
  • $\begingroup$ Have you seen this? mathematica.stackexchange.com/q/7887/219 $\endgroup$
    – faleichik
    Dec 25, 2015 at 17:53

3 Answers 3

Reset to default
10
$\begingroup$ 
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 *)
Improve this answer
$\endgroup$
4
  • 1
    $\begingroup$ Duh! For some reason I thought UpperTriangularize[] computed an upper triangular similar matrix. $\endgroup$
    – Igor Rivin
    Sep 14, 2014 at 16:53
  • 1
    $\begingroup$ 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. $\endgroup$
    – Igor Rivin
    Sep 14, 2014 at 18:28
  • 1
    $\begingroup$ @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. $\endgroup$
    – Wouter
    Sep 14, 2014 at 19:33
  • 2
    $\begingroup$ @Wouter, upper does not contain the main diagonal; second argument of UpperTriangularize shifts the triangle part off the main digonal. $\endgroup$
    – kglr
    Sep 14, 2014 at 19:41
5
$\begingroup$

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]];
    Statistics`Library`VectorToSymmetricMatrix[up, diag, n]
]

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

Improve this answer
$\endgroup$
5
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.