Как векторизовать расширенную индексацию со списком списков в NumPy?

Следующий код выполняется за 45 секунд при использовании чистого Python.

for iteration in range(maxiter):
    for node in range(n):
        for dest in adjacency_list[node]:
            rs[iteration + 1][dest] += beta * rs[iteration][node] / len(adjacency_list[node])

Но, просто инициализируя rs как numpy ndarray вместо списка списков python, код выполняется за 145 секунд. Я действительно не знаю, почему numpy занимает в 3 раза больше времени с индексацией этого массива.

Моя идея заключалась в том, чтобы векторизовать как можно больше вещей, но мне удалось векторизовать только умножение beta/len(adjacency_list[node]). Этот код работает за 77 секунд.

beta_over_out_degree = np.array([beta / len(al) for al in adjacency_list])
for iteration in range(1, maxiter + 1):
    r_next = np.full(shape=n, fill_value=(1 - beta) / n)
    f = beta_over_out_degree * r
    for i in range(n):
        r_next[adjacency_list[i]] += f[i]

    r = np.copy(r_next)
    rs[iteration] = np.copy(r)

Проблема в том, что adjacency_list — это список списков с разным размером столбцов, со 100 000 строк и 1-15 столбцов. Более стандартный подход с матрицей смежности, по крайней мере, как обычный ndarray, не вариант, так как для n=100 000 его форма (n,n) слишком велика для размещения в памяти.

Есть ли способ векторизовать, используя его индексы для расширенного индексирования numpy (возможно, превратив его в numpy ndarray)?

Я также был бы очень признателен за любые другие советы по скорости. Заранее спасибо!

РЕДАКТИРОВАТЬ: благодаря @stevemo мне удалось создать adjacency_matrix с функциональностью csr_matrix и использовать его для итеративного умножения. Теперь программа работает всего за 2 секунды!

for iteration in range(1, 101):
    rs[iteration] += rs[iteration - 1] * adjacency_matrix