transform the upper/lower triangular part of a symmetric matrix (2D array) into a 1D array and return it to the 2D format

Question

In this question it is explained how to access the lower and upper triagular parts of a given matrix, say:

m = np.matrix([[11, 12, 13],
               [21, 22, 23],
               [31, 32, 33]])

Here I need to transform the matrix in a 1D array, which can be done doing:

indices = np.triu_indices_from(m)
a = np.asarray( m[indices] )[-1]
#array([11, 12, 13, 22, 23, 33])

After doing a lot of calculations with a, changing its values, it will be used to fill a symmetric 2D array:

new = np.zeros(m.shape)
for i,j in enumerate(zip(*indices)):
    new[j]=a[i]
    new[j[1],j[0]]=a[i]

Returning:

array([[ 11.,  12.,  13.],
       [ 12.,  22.,  23.],
       [ 13.,  23.,  33.]])

Is there a better way to accomplish this? More especifically, avoiding the Python loop to rebuild the 2D array?

Octameter

Answer 1

The fastest and smartest way to put back a vector into a 2D symmetric array is to do this:

---

Case 1: No offset (k=0) i.e. upper triangle part includes the diagonal

import numpy as np

X = np.array([[1,2,3],[4,5,6],[7,8,9]])
#array([[1, 2, 3],
#       [4, 5, 6],
#       [7, 8, 9]])

#get the upper triangular part of this matrix
v = X[np.triu_indices(X.shape[0], k = 0)]
print(v)
# [1 2 3 5 6 9]

# put it back into a 2D symmetric array
size_X = 3
X = np.zeros((size_X,size_X))
X[np.triu_indices(X.shape[0], k = 0)] = v
X = X + X.T - np.diag(np.diag(X))
#array([[1., 2., 3.],
#       [2., 5., 6.],
#       [3., 6., 9.]])

The above will work fine even if instead of numpy.array you use numpy.matrix.

---

Case 2: With offset (k=1) i.e. upper triangle part does NOT include the diagonal

import numpy as np

X = np.array([[1,2,3],[4,5,6],[7,8,9]])
#array([[1, 2, 3],
#       [4, 5, 6],
#       [7, 8, 9]])

#get the upper triangular part of this matrix
v = X[np.triu_indices(X.shape[0], k = 1)] # offset
print(v)
# [2 3 6]

# put it back into a 2D symmetric array
size_X = 3
X = np.zeros((size_X,size_X))
X[np.triu_indices(X.shape[0], k = 1)] = v
X = X + X.T
#array([[0., 2., 3.],
#       [2., 0., 6.],
#       [3., 6., 0.]])

Implosion

Answer 2

Do you just want to form a symmetric array? You can skip the diagonal indices completely.

m=np.array(m)
inds = np.triu_indices_from(m,k=1)
m[(inds[1], inds[0])] = m[inds]

m

array([[11, 12, 13],
       [12, 22, 23],
       [13, 23, 33]])

Creating a symmetric array from a:

new = np.zeros((3,3))
vals = np.array([11, 12, 13, 22, 23, 33])
inds = np.triu_indices_from(new)
new[inds] = vals
new[(inds[1], inds[0])] = vals
new
array([[ 11.,  12.,  13.],
       [ 12.,  22.,  23.],
       [ 13.,  23.,  33.]])

Inherent

I have to process the middle term array a before returning to the symmetric 2D-array – Octameter

You can manipulate m[indup] as long as it returns a 1D numpy array in the correct order. Is there a particular issue with this? – Inherent

a have to come from array([11, 12, 13, 22, 23, 33]) to a corresponding 2D-array array([[11,12,13],[12,22,23],[13,23,33]]) – Octameter

Do you mean you want a symmetric array from a? – Inherent

Works only for 3x3 and smaller – Alda

Answer 3

You can use Array Creation Routines such as numpy.triu, numpy.tril, and numpy.diag to create a symmetric matrix from a triangular. Here's a simple 3x3 example.

a = np.array([[1,2,3],[4,5,6],[7,8,9]])
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])

a_triu = np.triu(a, k=0)
array([[1, 2, 3],
       [0, 5, 6],
       [0, 0, 9]])

a_tril = np.tril(a, k=0)
array([[1, 0, 0],
       [4, 5, 0],
       [7, 8, 9]])

a_diag = np.diag(np.diag(a))
array([[1, 0, 0],
       [0, 5, 0],
       [0, 0, 9]])

Add the transpose and subtract the diagonal:

a_sym_triu = a_triu + a_triu.T - a_diag
array([[1, 2, 3],
       [2, 5, 6],
       [3, 6, 9]])

a_sym_tril = a_tril + a_tril.T - a_diag
array([[1, 4, 7],
       [4, 5, 8],
       [7, 8, 9]])

Stier

I executed same steps but a turns out to be array([[ 22., 12., 13.], [ 12., 44., 23.], [ 13., 23., 66.]]) Diagonal elements are added to themselves – Pomander

@Madhav, numpy may have updated np.tril and np.triu since this post to include diagonals when k=0. Or my original response may have just been wrong. I've corrected it by subtracting the diagonal with np.diag. – Stier

Answer 4

Faster than the accepted solution for large matrices:

import numpy as np
X = np.array([[1,2,3],[4,5,6],[7,8,9]])
values = X[np.triu_indices(X.shape[0], k = 0)]

X2 = np.zeros_like(X)
triu_idx = np.triu_indices_from(X2)
X2[triu_idx], X2[triu_idx[::-1]] = values, values

Wandie

Answer 5

You can use squareform from SciPy. It adds a diagonal of zero which you can avoid as follows:

from scipy.spatial.distance import squareform


m = np.array([[1,2,3],[4,5,6],[7,8,9]])
#array([[1, 2, 3],
#       [4, 5, 6],
#       [7, 8, 9]])

a = m[np.triu_indices(m.shape[0], k = 0)]
# [1 2 3 5 6 9]

m_p = squareform(a) #with zero diagonal
#array([[0, 1, 2, 3],
#      [1, 0, 5, 6],
#      [2, 5, 0, 9],
#      [3, 6, 9, 0]])

#Get indices of diagonal inside a:
diag_ind = np.cumsum(np.insert(np.arange(m.shape[0], 1, -1),0,0))

#Get squareform of array without diagonal elements 
#then replace the zeros with the diagonal elements
m_pp = squareform(np.delete(a, diag_ind))
np.fill_diagonal(m_pp, a[diag_ind])

#m_pp = array([[1, 2, 3],
#      [2, 5, 6],
#      [3, 6, 9]])

The indices are obtained by a cumulative sum of countdown from shape value with zero inserted beforehand to include the zeroth element.

Aitchbone