Nested indexing in NumPy einsum?

Question

I'm trying to write the following expression using the einsum function in NumPy:

for j in range(100):
    p[j] = 0
    for i in range(100):
        if i!=j:
            p[j] += S[i,j]*B[T[i,j], i]

p.shape = (1,100) S.shape = (100,100) B.shape = (N,100) T.shape = (100,100)

N is chosen such that it exceeds the maximum value in the matrix T.

Is this possible to implement using the einsum function? If not, is there another way to optimise this?

I've searched the documentation, and it looks like it'll need some kind of broadcasting operation before being inserted into the einsum expression. But I can't get the expression right.

Chrysophylaxs · Accepted Answer · 2024-03-02 14:12:15Z

1

I don't think einsum has any nested indexing. But you can use vectorization:

import numpy as np

# Preliminaries:
S = np.random.rand(100, 100)
B = np.random.rand(10, 100)
T = np.random.randint(10, size=(100, 100))

# Solution:
idx = np.arange(100)[:, None]
arr = S * B[T, idx]

p = arr.sum(axis=0) - arr.diagonal()

I subtract the diagonal to account for the fact that you don't want to sum when i == j. You may choose to use keepdims=True in the sum to get p as shape (1, 100) instead.

answered Mar 2, 2024 at 14:12

Chrysophylaxs

6,5933 gold badges13 silver badges25 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

hpaulj Over a year ago

np.einsum('ij,ij->j', S, B[T, idx]) performs the sum/multiply, but probably without much of an improvement in speed.

Collectives™ on Stack Overflow

Nested indexing in NumPy einsum?

1 Answer 1

1 Comment

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

1 Comment

Your Answer

Sign up or log in

Post as a guest

Related