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]
]

Answer 1

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]

Answer 2

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]])

Answer 3

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 :

Answer 4

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)