Pythonic способ векторизации двойного суммирования - PullRequest
Новая
       11

Pythonic способ векторизации двойного суммирования

3 голосов
/ 28 марта 2019

Я пытаюсь преобразовать формулу двойного суммирования в код, но не могу выяснить правильное ее матричное / векторное представление.

Первое суммирование - от i до n, а второе - от j> i до n.

Я предполагаю, что есть гораздо более эффективный и питонский способ написания этого?

Я прибегнул к циклам nested, чтобы просто заставить его работать, но, как и ожидалось, он работает очень медленно сбольшой набор данных: 

def wapc_denom(weights, vols):
    x = []
    y = []

    for i, wi in enumerate(weights):
        for j, wj in enumerate(weights):
            if j > i:
                x.append(wi * wj * vols[i] * vols[j])
        y.append(np.sum(x))

    return np.sum(y)

Редактировать:

Используя руководство из ответа smci, я думаю, что у меня есть потенциальное решение: 

def wapc_denom2(weights, vols):
    return np.sum(np.tril(np.outer(weights, vols.T)**2, k=-1))

Ответы [ 3 ]

3 голосов
/ 29 марта 2019

Предполагая, что вы хотите посчитать каждый член только один раз (для этого вам нужно переместить x = [] во внешний цикл), один дешевый способ вычисления суммы будет

Создать фиктивные данные 

weights = np.random.random(10)
vols = np.random.random(10)

Выполните расчет 

wv = weights * vols
result = (wv.sum()**2 - wv@wv) / 2

Убедитесь, что это то же самое 

def wapc_denom(weights, vols):
    y = []

    for i, wi in enumerate(weights):
        x = []
        for j, wj in enumerate(weights):
            if j > i:
                x.append(wi * wj * vols[i] * vols[j])
        y.append(np.sum(x))

    return np.sum(y)

assert np.allclose(result, wapc_denom(weights, vols))

Почему это работает?

Что мы делаем, это вычисляем суммуиз полной матрицы вычтите диагональ и разделите на два.Это дешево, потому что легко проверить, что сумма внешнего произведения является просто произведением суммированных факторов.

1 голос
/ 28 марта 2019
  • wi * wj * vols[i] * vols[j] является контрольным.vols - это еще один вектор, поэтому сначала вы хотите вычислить вектор wv = w * vols
  • , а затем (wj * vols[j]) * (wi * vols[i]) = wv^T * wv - это ваше выражение (матричный внешний продукт);это вектор столбца * вектор строки.Но на самом деле вы хотите только сумму.Так что я не вижу необходимости в построении вектора y.append(np.sum(x)), вы все равно будете его только суммировать np.sum(y)
  • , а часть if j > i означает, что вам нужна только сумма нижнего треугольникачасть и исключить диагональ.
    • РЕДАКТИРОВАТЬ: результат полностью определяется только из wv, я не думал, что нам нужна матрица, чтобы получить сумму, и нам не нужна диагональ;@PaulPanzer нашел самое компактное выражение.
0 голосов
/ 28 марта 2019

Вы можете использовать триангуляции в numpy, отметьте np.triu и np.meshgrid. Есть: 

np.product(np.triu(np.meshgrid(weights,weights), 1) * np.triu(np.meshgrid(vols,vols), 1),0).sum(1).cumsum().sum()

Пример: 

w = np.arange(4) +1
v = np.array([1,3,2,2])
print(np.triu(np.meshgrid(w,w), k=1))
>>array([[[0, 2, 3, 4],
        [0, 0, 3, 4],
        [0, 0, 0, 4],
        [0, 0, 0, 0]],

       [[0, 1, 1, 1],
        [0, 0, 2, 2],
        [0, 0, 0, 3],
        [0, 0, 0, 0]]])

# example of product + triu + meshgrid (your x values):
print(np.product(np.triu(np.meshgrid(w,w), 1) * np.triu(np.meshgrid(v,v), 1),0))
>>array([[ 0,  6,  6,  8],
   [ 0,  0, 36, 48],
   [ 0,  0,  0, 48],
   [ 0,  0,  0,  0]])

print(np.product(np.triu(np.meshgrid(w,w), 1) * np.triu(np.meshgrid(v,v), 1),0).sum(1).cumsum().sum())
>> 428
print(wapc_denom(w, v))
>> 428
...