O.rka - 1 year ago 185
Python Question

# How to convert triangle matrix to square in NumPy?

I'm doing some computations on a full matrix that is redundant (i.e. can be a triangle matrix without losing info). I realized I can compute only the lower portion of the triangle for faster results. How can I project the lower triangle into the upper once I'm done?

In other words, how can I reverse the

`np.tril`
method?

``````print DF_var.as_matrix()
# [[1 1 0 1 1 1 0 1 0 0 0]
#  [1 1 1 1 1 0 1 0 1 1 1]
#  [0 1 1 0 0 0 0 0 0 0 0]
#  [1 1 0 1 0 0 0 0 0 0 0]
#  [1 1 0 0 1 0 0 0 0 0 0]
#  [1 0 0 0 0 1 1 0 0 0 0]
#  [0 1 0 0 0 1 1 0 0 0 0]
#  [1 0 0 0 0 0 0 1 1 0 0]
#  [0 1 0 0 0 0 0 1 1 0 0]
#  [0 1 0 0 0 0 0 0 0 1 0]
#  [0 1 0 0 0 0 0 0 0 0 1]]
print np.tril(DF_var.as_matrix())
# [[1 0 0 0 0 0 0 0 0 0 0]
#  [1 1 0 0 0 0 0 0 0 0 0]
#  [0 1 1 0 0 0 0 0 0 0 0]
#  [1 1 0 1 0 0 0 0 0 0 0]
#  [1 1 0 0 1 0 0 0 0 0 0]
#  [1 0 0 0 0 1 0 0 0 0 0]
#  [0 1 0 0 0 1 1 0 0 0 0]
#  [1 0 0 0 0 0 0 1 0 0 0]
#  [0 1 0 0 0 0 0 1 1 0 0]
#  [0 1 0 0 0 0 0 0 0 1 0]
#  [0 1 0 0 0 0 0 0 0 0 1]]
``````

How to convert it back to a full matrix?

Assuming `A` as the input array, few methods are listed below.

Approach #1 : Using `np.triu` on a transposed version of `A` -

``````np.triu(A.T,1) + A
``````

Approach #2 : Avoid `np.triu` with direct summation between A.T and A and then indexing to set diagonal elements -

``````out = A.T + A
idx = np.arange(A.shape[0])
out[idx,idx] = A[idx,idx]
``````

Approach #3 : Same as previous one, but compact using in-builts for indexing -

``````out = A.T + A
np.fill_diagonal(out,np.diag(A))
``````

Approach #4 : Same as previous one, but with boolean indexing to set diagonal elements -

``````out = A.T + A
``````

Approach #5 : Using mask based selection for diagonal elements with `np.where` -

``````np.where(np.eye(A.shape[0],dtype=bool),A,A.T+A)
``````

Approach #6 : Using mask based selection for all elements with `np.where` -

``````np.where(np.triu(np.ones(A.shape[0],dtype=bool),1),A.T,A)
``````

Runtime tests

Functions -

``````def func1(A):
return np.triu(A.T,1) + A

def func2(A):
out = A.T + A
idx = np.arange(A.shape[0])
out[idx,idx] = A[idx,idx]
return out

def func3(A):
out = A.T + A
np.fill_diagonal(out,np.diag(A))
return out

def func4(A):
out = A.T + A
return out

def func5(A):
return np.where(np.eye(A.shape[0],dtype=bool),A,A.T+A)

def func6(A):
return np.where(np.triu(np.ones(A.shape[0],dtype=bool),1),A.T,A)
``````

Timings -

``````In [140]: # Input array
...: N = 5000
...: A = np.tril(np.random.randint(0,9,(N,N)))
...:

In [141]: %timeit func1(A)
...: %timeit func2(A)
...: %timeit func3(A)
...: %timeit func4(A)
...: %timeit func5(A)
...: %timeit func6(A)
...:
1 loops, best of 3: 617 ms per loop
1 loops, best of 3: 354 ms per loop
1 loops, best of 3: 354 ms per loop
1 loops, best of 3: 395 ms per loop
1 loops, best of 3: 597 ms per loop
1 loops, best of 3: 440 ms per loop
``````

Looks like the approaches # 2 & #3 are pretty efficient!

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download