Numpy make the product between all elemens and then insert into a triangular 2d array

Question

Suppose we got a 1D array below

arr = np.array([a,b,c])

The first thing I need to do is the make the product of all of the elments, i.e

[ab,ac,bc]

Then construct a 2d triangular array with this element

[
[a,ab,ac],
[0,b,bc],
[0,0,c]
]

score 5 · Accepted Answer · answered Jul 10 '20 at 02:09

5

Create a diagonal with your 1-D array and fill the upper triangle of it with upper triangle of outer:

out = np.diag(arr)
#upper triangle indices
uidx = np.triu_indices(arr.size,k=1)
#replacing upper triangle with outer
out[uidx]=np.outer(arr,arr)[uidx]

answered Jul 10 '20 at 02:09

Ehsan

12,072
2
20
33

1

This will be slower at scale, because you're allocating memory for an `(arr.size, arr.size)` array twice. If this calculation is being repeated on large arrays in a loop, this slowdown could become materially relevant. – Nick Becker Jul 10 '20 at 03:22
@NickBecker I disagree. Computationally speaking this answer and the post below should be of the same order. Memory-wise it uses double. They both are essentially the same. It is up to the reader which one they find more readable. – Ehsan Jul 10 '20 at 03:29
1

Agreed, they are the same in terms of complexity. But the actual elapsed time is likely to be slower at scale/repeated due to the cost of the allocations. Timing a 20 iteration for loop with `arr = np.array(range(1,5000))`, I see a 30% difference in elapsed time. To be fair, this probably isn't relevant. Both are good answers :) – Nick Becker Jul 10 '20 at 03:39
@NickBecker Thank you for the test run. When we look at improving performance, usually 30% will not be sth a user worries about (mostly in terms of x times is considered significant). However, of course the faster the better and your point is valid. Upvoting your answer. – Ehsan Jul 10 '20 at 03:44
1

Your answer is also good, as it's more clear in terms of readability. Upvoting as well – Nick Becker Jul 10 '20 at 03:50
@NickBecker Don't want to sound snobbish in any way, but just to share my experience thus far working with upper/lower indices, we are better off masking if possible, as it seems you guys care about efficiency. Here's a related post - https://stackoverflow.com/questions/47314754/how-to-get-triangle-upper-matrix-without-the-diagonal-using-numpy. – Divakar Jul 10 '20 at 06:08
Look at that. Even better. These kinds of improvements can absolutely matter in some workloads, so don't worry about sounding snobbish at all – Nick Becker Jul 10 '20 at 14:53

score 3 · Answer 2 · answered Jul 10 '20 at 00:49

One way to do this is to calculate the outer product of your 1d array and then use masking informed by the knowledge that you only want the upper triangle of the 2d triangular matrix.

import numpy as np

a = np.array([5,4,3])
n = len(a)

outer = np.outer(a, a)
outer[np.tril_indices(n)] = 0
outer[np.diag_indices(n)] = a

outer
array([[ 5, 20, 15],
       [ 0,  4, 12],
       [ 0,  0,  3]])

Divakar · Answer 3 · 2020-07-10T06:24:12.450

We can use masking to achieve our desired result, like so -

def upper_outer(a):
    out = a[:,None]*a
    out[np.tri(len(a), k=-1, dtype=bool)] = 0
    np.fill_diagonal(out,a)
    return out

Sample run -

In [84]: a = np.array([3,6,2])

In [86]: upper_outer(a)
Out[86]: 
array([[ 3, 18,  6],
       [ 0,  6, 12],
       [ 0,  0,  2]])

Benchmarking

Other approaches :

# @Nick Becker's soln
def tril_diag(a):
    n = len(a)
    outer = np.outer(a, a)
    outer[np.tril_indices(n)] = 0
    outer[np.diag_indices(n)] = a
    return outer

# @Ehsan's soln
def triu_outer(arr):
    out = np.diag(arr)
    uidx = np.triu_indices(arr.size,k=1)
    out[uidx]=np.outer(arr,arr)[uidx]
    return out

Using benchit package (few benchmarking tools packaged together; disclaimer: I am its author) to benchmark proposed solutions.

import benchit
in_ = [np.random.rand(n) for n in [10,100,200,500,1000,5000]]
funcs = [upper_outer, tril_diag, triu_outer]
t = benchit.timings(funcs, in_)
t.rank()
t.plot(logx=True, save='timings.png')

For large datasets, we can also use numexpr to leverage multi-cores -

import numexpr as ne

def upper_outer_v2(a):
    mask = ~np.tri(len(a), dtype=bool)
    out = ne.evaluate('a2D*a*mask',{'a2D':a[:,None], 'a':a, 'mask':mask})
    np.fill_diagonal(out,a)
    return out

New timings plot :

I like your answer and analysis. Upvoted. Interesting note on masks vs. indexing. — Ehsan, Jul 10 '20 at 06:55

Paul Panzer · Answer 4 · 2020-07-10T08:05:49.030

3

There is a blas function for (almost) that:

# example
a = np.array([1.,2.,5.])

from scipy.linalg.blas import dsyr

# apply blas function; transpose since blas uses FORTRAN order
out = dsyr(1,a,1).T
# fix diagonal
out.reshape(-1)[::a.size+1] = a
out
# array([[ 1.,  2.,  5.],
#        [ 0.,  2., 10.],
#        [ 0.,  0.,  5.]])

benchit (thanks @Divakar)

edited Jul 10 '20 at 08:05

answered Jul 10 '20 at 07:53

Paul Panzer

51,835
3
54
99

You got that working huh! Gotta lookup that new thing - `dsyr`. – Divakar Jul 10 '20 at 08:20
1

@Divakar Yep, neat! Love the progress bars. – Paul Panzer Jul 10 '20 at 08:25
Lol, beats all (gotta be careful and say so far) :-) – Ehsan Jul 10 '20 at 08:26
The future is now :) – Nick Becker Jul 10 '20 at 15:00

Numpy make the product between all elemens and then insert into a triangular 2d array

4 Answers4

Benchmarking