Making a numpy ndarray matrix symmetric

Question

I have a 70x70 numpy ndarray, which is mainly diagonal. The only off-diagonal values are the below the diagonal. I would like to make the matrix symmetric.

As a newcomer from Matlab world, I can't get it working without for loops. In MATLAB it was easy:

W = max(A,A')

where A' is matrix transposition and the max() function takes care to make the W matrix which will be symmetric.

Is there an elegant way to do so in Python as well?

EXAMPLE The sample A matrix is:

1 0 0 0
0 2 0 0
1 0 2 0
0 1 0 3

The desired output matrix W is:

1 0 1 0
0 2 0 1
1 0 2 0
0 1 0 3

Centipede

is answered here – Expel

Answer 1

Found a following solution which works for me:

import numpy as np
W = np.maximum( A, A.transpose() )

Centipede

The problem is that if your matrix A is large, then this wastes a lot of space, because I suppose it creates a second matrix, or does it not? – Rafat

This will not work if the off-diagonal elements you're trying to mirror are negative, no? – Rabbinical

When dealing with negative coefficients, you would simply do W = (A + A.transpose()) / 2 – Haft

Answer 2

Use the NumPy tril and triu functions as follows. It essentially "mirrors" elements in the lower triangle into the upper triangle.

import numpy as np
A = np.array([[1, 0, 0, 0], [0, 2, 0, 0], [1, 0, 2, 0], [0, 1, 0, 3]])
W = np.tril(A) + np.triu(A.T, 1)

tril(m, k=0) gets the lower triangle of a matrix m (returns a copy of the matrix m with all elements above the kth diagonal zeroed). Similarly, triu(m, k=0) gets the upper triangle of a matrix m (all elements below the kth diagonal zeroed).

To prevent the diagonal being added twice, one must exclude the diagonal from one of the triangles, using either np.tril(A) + np.triu(A.T, 1) or np.tril(A, -1) + np.triu(A.T).

Also note that this behaves slightly differently to using maximum. All elements in the upper triangle are overwritten, regardless of whether they are the maximum or not. This means they can be any value (e.g. nan or inf).

Chaparral

Answer 3

For what it is worth, using the MATLAB's numpy equivalent you mentioned is more efficient than the link @plonser added.

In [1]: import numpy as np
In [2]: A = np.zeros((4, 4))
In [3]: np.fill_diagonal(A, np.arange(4)+1)
In [4]: A[2:,:2] = np.eye(2)

# numpy equivalent to MATLAB:
In [5]: %timeit W = np.maximum( A, A.T)
100000 loops, best of 3: 2.95 µs per loop

# method from link
In [6]: %timeit W = A + A.T - np.diag(A.diagonal())
100000 loops, best of 3: 9.88 µs per loop

Timing for larger matrices can be done similarly:

In [1]: import numpy as np
In [2]: N = 100
In [3]: A = np.zeros((N, N))
In [4]: A[2:,:N-2] = np.eye(N-2)
In [5]: np.fill_diagonal(A, np.arange(N)+1)
In [6]: print A
Out[6]: 
array([[   1.,    0.,    0., ...,    0.,    0.,    0.],
       [   0.,    2.,    0., ...,    0.,    0.,    0.],
       [   1.,    0.,    3., ...,    0.,    0.,    0.],
       ..., 
       [   0.,    0.,    0., ...,   98.,    0.,    0.],
       [   0.,    0.,    0., ...,    0.,   99.,    0.],
       [   0.,    0.,    0., ...,    1.,    0.,  100.]])

# numpy equivalent to MATLAB:
In [6]: %timeit W = np.maximum( A, A.T)
10000 loops, best of 3: 28.6 µs per loop

# method from link
In [7]: %timeit W = A + A.T - np.diag(A.diagonal())
10000 loops, best of 3: 49.8 µs per loop

And with N = 1000

# numpy equivalent to MATLAB:
In [6]: %timeit W = np.maximum( A, A.T)
100 loops, best of 3: 5.65 ms per loop

# method from link
In [7]: %timeit W = A + A.T - np.diag(A.diagonal())
100 loops, best of 3: 11.7 ms per loop

Gravure

It's roughly 3x times. Does it hold for larger matrices, say 100x100 or 1000x1000? – Centipede

@xeon I added timings for the sizes you mentioned. These timings are on my machine, YMMV. The difference of ~2x comes from the number of operations performed. Both use numpy objects and numpy methods. Similar to MATLAB, numpy was designed to optimize matrix operations so numpy methods are typically the most efficient. – Gravure

I swapped A.transpose() for A.T to point out they are the same and one requires less typing (docs.scipy.org/doc/numpy/reference/generated/…). – Gravure

Answer 4

can get symmetric positive-definite

import pandas as pd
from sklearn.datasets import make_spd_matrix    # spd - symmetric positive-definite matrix

spd = make_spd_matrix(n_dim=3, random_state=1)
print(pd.DataFrame(spd))

Saddler

Answer 5

All solutions so far use a floating point operation (for every matrix element), either + or max, where you don’t actually need one.

You could just use indexing to copy the transposed bit of the matrix:

idx = np.triu_indices(A.shape[0], 1, A.shape[1])
A[idx] = A.T[idx]

However indexing is a bit slow, so for an alternate (and shorter) way of doing the same, you could use np.where and generate which values to pick with np.tri:

np.where(np.tri(*A.shape, dtype=bool), A, A.T)

Did some benchmarking on a random 1000 x 1000 matrix:

Code run with timeit.timeit('…', globals={'np': np, 'A': A}, number=1000)	time
np.maximum(A, A.T) *	5.29
A + np.triu(A.T, 1)	2.61
A + A.T - np.diag(A.diagonal())	5.71
np.tril(A, -1) + A.T	2.45
idx = np.triu_indices(A.shape[0], 1, A.shape[1]) ; A[idx] = A.T[idx]	21.82
np.where(np.tri(*A.shape, dtype=bool), A, A.T)	2.00

* as noted above this version only works for positive numbers only.

Cerebroside