CODEWITHSUNDEEP
pythonnumpymatrixrandom
luw
luw

Reputation: 217

Is there an efficient way to generate a symmetric random matrix?

Is there an efficient way (using numpy) to generate a symmetric random matrix with entries uniformly distributed in [0,1)?

Upvotes: 2

Views: 3407

Answers (2)

Paul Panzer
Paul Panzer

Reputation: 53029

Here is a method using scipy.spatial.distance.squareform:

squareform switches back and forth between the full and "compressed" forms of a symmetric matrix:

>>> full = squareform(np.arange(1,11))
>>> full
array([[ 0,  1,  2,  3,  4],
       [ 1,  0,  5,  6,  7],
       [ 2,  5,  0,  8,  9],
       [ 3,  6,  8,  0, 10],
       [ 4,  7,  9, 10,  0]])
>>> squareform(full)
array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10])

As it was designed with distance matrices in mind it leaves the diagonal at zero, so we have to manually fill that in. For that we use einsum which used the way we do returns a writable view of the diagonal,

>>> from scipy.spatial.distance import squareform
>>> 
>>> N = 5
>>> a = squareform(np.random.random(N*(N-1)//2))
>>> np.einsum('ii->i', a)[:] = np.random.random(N)
>>> a
array([[0.29946651, 0.3636706 , 0.00708741, 0.87536594, 0.62197293],
       [0.3636706 , 0.31774527, 0.05597852, 0.10800514, 0.99871399],
       [0.00708741, 0.05597852, 0.83912235, 0.86241008, 0.01806965],
       [0.87536594, 0.10800514, 0.86241008, 0.11039534, 0.64213608],
       [0.62197293, 0.99871399, 0.01806965, 0.64213608, 0.84755054]])

Upvotes: 2

Jakob Drachmann Havtorn
Jakob Drachmann Havtorn

Reputation: 128

Let U be a square matrix of uniformly distributed random numbers. You can then add the lower triangular part of U with itself transposed (including the diagonal only once) to get a symmetric matrix with random numbers from the same distribution as U.

import numpy as np 

U = np.random.uniform(low=0, high=1.0, size=(1000, 1000))
S = np.tril(U) + np.tril(U, -1).T

print(np.histogram(S.flatten()))
print(np.histogram(S[0,:]))
print(np.histogram(S[:,0]))

The matrix as a whole as well as any row or column will be uniformly distributed in [0,1) by the documentation for np.random.uniform

Speed-wise I get

%timeit U = np.random.uniform(low=0, high=1.0, size=(1000, 1000))
10.6 ms ± 46.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%timeit S = np.tril(U) + np.tril(U, -1).T
5.76 ms ± 75.2 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

As other people have noted, you can also do

S = (U + U.T) / 2

to get symmetry but this will give you triangular distributed random numbers in the off-diagonal since you are summing two uniform random variables.

Upvotes: 5

Related Questions