2

I have the upper triangle entries (including diagonal) of a symmetric matrix in a flat list (concatenated rows), and I want to use them to fill in the full matrix, including the lower triangle. What's the fastest way to do this?

Here is my current approach. Seems like a lot of work for such a simple operation.

import numpy as np
def utri2mat(utri,ntotal):
    iu1 = np.triu_indices(ntotal)
    ret = np.zeros([ntotal,ntotal])
    ret[iu1] = utri
    ret = ret + ret.transpose() - np.diag(ret.diagonal())
    return ret
Divakar
  • 218,885
  • 19
  • 262
  • 358
cgreen
  • 351
  • 1
  • 17
  • I don't see how it can be made any simpler. Your symmetric matrix as 3 distinct parts, upper and lower triangle and diagonal. In one way or other you have to set them individually with your flat list. – hpaulj Dec 12 '15 at 01:13

4 Answers4

1

This version is a slight variation of yours:

import numpy as np

def utri2mat(utri):
    n = int(-1 + np.sqrt(1 + 8*len(utri))) // 2
    iu1 = np.triu_indices(n)
    ret = np.empty((n, n))
    ret[iu1] = utri
    ret.T[iu1] = utri
    return ret

I replaced

    ret = ret + ret.transpose() - np.diag(ret.diagonal())

with an assignment of utri directly to the transpose of ret:

    ret.T[iu1] = utri

I also removed the argument ntotal and instead computed what n has to be based on the length of utri.

Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214
1

Inspired by this solution, you can use boolean indexing to set elements and as such might be pretty efficient. Here's one way to implement it -

def mask_based_utri2mat(utri,ntotal):
    # Setup output array
    out = np.empty((ntotal,ntotal))

    # Create upper triang. mask
    mask = np.triu(np.ones((ntotal,ntotal),dtype=bool))

    # Set upper triang. elements with mask
    out[mask] = utri

    # Set lower triang. elements with transposed mask
    out.T[mask] = utri
    return out

Runtime tests -

In [52]: # Inputs
    ...: ntotal = 100
    ...: utri = np.random.rand(np.triu_indices(ntotal)[0].size)
    ...: 

In [53]: np.allclose(mask_based_utri2mat(utri,ntotal),utri2mat(utri,ntotal))
Out[53]: True

In [54]: %timeit utri2mat(utri,ntotal)
1000 loops, best of 3: 270 µs per loop

In [55]: %timeit mask_based_utri2mat(utri,ntotal)
10000 loops, best of 3: 127 µs per loop

In [56]: # Inputs
    ...: ntotal = 1000
    ...: utri = np.random.rand(np.triu_indices(ntotal)[0].size)
    ...: 

In [57]: np.allclose(mask_based_utri2mat(utri,ntotal),utri2mat(utri,ntotal))
Out[57]: True

In [58]: %timeit utri2mat(utri,ntotal)
10 loops, best of 3: 53.9 ms per loop

In [59]: %timeit mask_based_utri2mat(utri,ntotal)
100 loops, best of 3: 15.1 ms per loop
Community
  • 1
  • 1
Divakar
  • 218,885
  • 19
  • 262
  • 358
  • In my tests `utri2mat` works with a flat array, your `mask_based_utri2mat` requires a `(ntotal,ntotal)` array. The problem is with `out[mask] = utri`. `out[np.where(mask)]=utri` does it correctly. `np.triu_indices` does what `where`. – hpaulj Dec 12 '15 at 20:37
  • @hpaulj Not sure if I follow correctly - Even OP's `utri2mat` generates that array with `ret = np.zeros([ntotal,ntotal])`, right? – Divakar Dec 12 '15 at 20:42
  • There was a July 2014 patch to `np.triu`: https://github.com/numpy/numpy/commit/1392888ec4ddb37d4fa7bbb9231712af0dda4ea6 . Older versions returned an `int` array rather than a `bool`. Without that patch Warren's answer is better. – hpaulj Dec 12 '15 at 21:17
  • @hpaulj Better runtime wise? I am on the latest `1.10.1` NumPy version and tried with `ntotal = 100` and the runtimes with `mask_based_utri2mat` are `2x+` better than both OP and Warren's answer. For `n = 1000` it's `3x+`. – Divakar Dec 13 '15 at 08:27
  • I wonder if ` mask=np.tri(n).astype(bool)` is any faster. Warren's solution may slower because it uses `where` indexing instead of boolean. `np.tri` is the base 'tri' function. – hpaulj Dec 13 '15 at 11:13
1

Here's my nomination for a faster, and possibly better, way to make a symmetric matrix from flat values:

def make_sym(val, n):
    # uses boolean mask
    # uses the same lower tri as np.triu
    mask = ~np.tri(5,k=-1,dtype=bool)
    out = np.zeros((n,n),dtype=val.dtype)
    out[mask] = val
    out.T[mask] = val
    return out

testing:

In [939]: val=np.arange(1,16)
In [940]: make_sym(val, 5)
Out[940]: 
array([[ 1,  2,  3,  4,  5],
       [ 2,  6,  7,  8,  9],
       [ 3,  7, 10, 11, 12],
       [ 4,  8, 11, 13, 14],
       [ 5,  9, 12, 14, 15]])

Like the other answers it uses out.T[] to assign the lower triangle.

Warren's answer uses np.triu_indices, which are the where values. This type of indexing is a bit slower than boolean masking.

But as I noted the np.triu that Divakar uses does not return a boolean mask in earlier numpy versions (e.g. 1.9). This is what prompted me to dig into the issue.

In 1.10 this function was rewritten as:

mask = np.tri(*m.shape[-2:], k=k-1, dtype=bool)
return np.where(mask, np.zeros(1, m.dtype), m)

I gain a bit of speed by replacing the where with ~mask. Same result, but just cutting out an intermediate step.

hpaulj
  • 221,503
  • 14
  • 230
  • 353
0

Assuming that you have a vector containing the upper triangular values of a symmetric matrix (n x n) then you can re-build the full matrix as follows: 

import numpy as np

# dimension of the full matrix
n = 80

# artificial upper triangle entries n(n-1) / 2 if matrix is symmetric
entries = np.array(range((80*79) / 2))

full_matrix = np.zeros((n,n))
inds = np.triu_indices_from(full_matrix, k = 1)
full_matrix[inds] = entries
full_matrix[(inds[1], inds[0])] = entries

print(full_matrix)

Verify:

np.allclose(full_matrix, full_matrix.T)
seralouk
  • 30,938
  • 9
  • 118
  • 133