How can I compute lower-triangle differences for a NumPy array using vectorized subtraction?

A naive approach that does more work than necessary is to compute the entire difference and select the elements you need:

PYTHONCopied
np.where(np.arange(U.size)[:, None] > np.arange(U.size), U[:, None] - U, np.nan)

This is one of the times where np.where is actually useful over a simple mask, although it can be done with a mask as well:

PYTHONCopied
result = np.full((U.size, U.size), np.nan)
index = np.arange(U.size)
mask = index[:, None] > index
result[mask] = [U[:, None] - U][mask]

A more efficient approach might be to use the indices more directly to index into the source:

PYTHONCopied
result = np.full((U.size, U.size), np.nan)
r, c = np.tril_indices(U.size, k=-1)
result[r, c] = U[c] - U[r]

Not most efficient, but only few lines of code:

• Calculate the full matrix via a broadcasted u.
• Set the upper triangle and diagonal to nan.

PYTHONCopied
import numpy as np

u = np.array([3, 5, 7, 9, 11], dtype=float)

result = u - u[:, np.newaxis]
result[np.triu_indices(len(u), k=0)] = np.nan
print(result)
# [[nan nan nan nan nan]
#  [-2. nan nan nan nan]
#  [-4. -2. nan nan nan]
#  [-6. -4. -2. nan nan]
#  [-8. -6. -4. -2. nan]]

Here are timing results using the accepted answer as a baseline – the version using where and the version using masking seem to have an advantage for smaller sizes; differences appear more and more negligible for larger sizes:

Code used for timing (note that I swapped the order of subtraction for baseline_where() and baseline_mask() to be in line with OP's question, also I fixed a typo ((…)[mask] rather than […][mask]) in baseline_mask():

PYTHONCopied
import matplotlib.pyplot as plt
import numpy as np
from timeit import Timer

def proposed(U):
    result = U.astype(float) - U[:, np.newaxis]
    result[np.triu_indices(len(U), k=0)] = np.nan
    return result

def baseline_where(U):
    return np.where(np.arange(U.size)[:, None] > np.arange(U.size),
                    U - U[:, None], np.nan)

def baseline_mask(U):
    result = np.full((U.size, U.size), np.nan)
    index = np.arange(U.size)
    mask = index[:, None] > index
    result[mask] = (U - U[:, None])[mask]
    return result

def baseline_direct_indexing(U):
    result = np.full((U.size, U.size), np.nan)
    r, c = np.tril_indices(U.size, k=-1)
    result[r, c] = U[c] - U[r]
    return result

rand = np.random.default_rng(seed=42)
n = 100  # number of timing runs
sizes = np.round(np.logspace(0, 4, num=9)).astype(int)
timings = {}
for size in sizes:
    u = rand.integers(100, size=size)
    for fct in proposed, baseline_where, baseline_mask, baseline_direct_indexing:
        fct_name = fct.__name__
        fct_timings = timings.get(fct_name, [])
        fct_timings.append(Timer(lambda: fct(u)).timeit(n))
        timings[fct_name] = fct_timings

plt.loglog()
plt.xlabel("len(U)")
plt.ylabel(f"calculation time ({n=} runs) (s)")
for fct_name, fct_timings in timings.items():
    plt.plot(sizes, fct_timings, "x", label=fct_name)
plt.legend()