Оптимизация математического расчета - PullRequest
Новая
       16

Оптимизация математического расчета

0 голосов
/ 08 ноября 2018

У меня есть модель для четырех возможностей покупки пары предметов (покупка обоих, ни одного или только одного), и мне нужно оптимизировать (псевдо) функцию логарифмического правдоподобия.Частично это, конечно, является вычислением / определением функции правдоподобия псевдо-лога.

Ниже приведен мой код, где бета представляет собой двумерный вектор для каждого клиента (есть U клиентов иU различных бета-векторов), X - это двумерный вектор для каждого элемента (различный для каждого из N элементов), а гамма - симметричная матрица со скалярным значением гамма (i, j) для каждой пары элементов.А df - это кадр данных покупок - одна строка для каждого клиента и N столбцов для товаров.

Мне кажется, что все эти циклы неэффективны и занимают слишком много времени, но я неуверен, как ускорить этот расчет и был бы признателен за любую помощь, улучшая его.Заранее спасибо! 

def pseudo_likelihood(Args):
    Beta = np.reshape(Args[0:2*U], (U, 2))
    Gamma = np.reshape(Args[2*U:], (N,N))
    L = 0
    for u in range(0,U,1):
        print datetime.datetime.today(), " for user {}".format(u)
        y = df.loc[u][1:]
        beta_u = Beta[u,:]
        for l in range(N):
            print datetime.datetime.today(), " for item {}".format(l)
            for i in range(N-1):
                if i == l:
                    continue
                for j in range(i+1,N):
                    if (y[i] == y[j]):
                        if (y[i] == 1):
                            L += np.dot(beta_u,(x_vals.iloc[i,1:]+x_vals.iloc[j,1:])) + Gamma[i,j] #Log of the exponent of this expression
                        else:
                            L += np.log(
                                1 - np.exp(np.dot(beta_u, (x_vals.iloc[i, 1:] + x_vals.iloc[j, 1:])) + Gamma[i, j])
                                - np.exp(np.dot(beta_u, x_vals.iloc[i, 1:])) * (
                                            1 - np.exp(np.dot(beta_u, x_vals.iloc[j, 1:])))
                                - np.exp(np.dot(beta_u, x_vals.iloc[j, 1:])) * (
                                            1 - np.exp(np.dot(beta_u, x_vals.iloc[i, 1:]))))
                    else:
                        if (y[i] == 1):
                            L += np.dot(beta_u,x_vals.iloc[i,1:]) + np.log(1 - np.exp(np.dot(beta_u,x_vals.iloc[j,1:])))
                        else:
                            L += (np.dot(beta_u, x_vals.iloc[j,1:])) + np.log(1 - np.exp(np.dot(beta_u, x_vals.iloc[i,1:])))

            L -= (N-2)*np.dot(beta_u,x_vals.iloc[l,1:])
            for k in range(N):
                if k != l:
                    L -= np.dot(beta_u, x_vals.iloc[k,1:])

    return -L

Для добавления / уточнения - я использую этот расчет, чтобы оптимизировать и найти параметры бета и гаммы, которые сгенерировали данные для этой функции псевдослучайного вероятности.
Я используюscipy optimize.minimize с помощью метода «Пауэлл».

1 Ответ

0 голосов
/ 15 ноября 2018

Обновление для всех, кто заинтересован- Я нашел numpy.einsum, чтобы ускорить вычисления здесь более чем на 90%.

np.einsum выполняет матричные / векторные операции с использованием нотации Эйнштейна. Напомним, что для двух матриц A, B их произведение можно представить в виде суммы 

a_ij*b_jk

т.е. элемент ik матрицы AB является суммой по j для a_ij * b_jk

Используя функцию einsum, я мог заранее рассчитать все значения, необходимые для итеративного расчета, сэкономив драгоценное время и сотни, если не тысячи, ненужных вычислений. Я переписал код следующим образом: 

def pseudo_likelihood(Args):
    Beta = np.reshape(Args[0:2*U], (U,2))
    Gamma = np.reshape(Args[2*U:], (N,N))
    exp_gamma = np.exp(Gamma)
    L = 0
    for u in xrange(U):
        y = df.loc[u][1:]
        beta_u = Beta[u,:]
        beta_dot_x = np.einsum('ij,j',x_vals[['V1','V2']],beta_u)
        exp_beta_dot_x = np.exp(beta_dot_x)
        log_one_minus_exp = np.log(1 - exp_beta_dot_x)
        for l in xrange(N):
            for i in xrange(N-1):
                if i == l:
                    continue
                for j in xrange(i+1,N):
                    if (y[i] == y[j]):
                        if (y[i] == 1):
                            L += beta_dot_x[i] + beta_dot_x[j] + Gamma[i,j] #Log of the exponent of this expression
                        else:
                            L += math.log(
                                1 - exp_beta_dot_x[i]*exp_beta_dot_x[j]*exp_gamma[i,j]
                                - exp_beta_dot_x[i] * (1 - exp_beta_dot_x[j])
                                - exp_beta_dot_x[j] * (1 - exp_beta_dot_x[i]))
                    else:
                        if (y[i] == 1):
                            L += beta_dot_x[i] + log_one_minus_exp[j]
                        else:
                            L += (beta_dot_x[j]) + log_one_minus_exp[i]

            L -= (N-2)*beta_dot_x[l]
            for k in xrange(N):
                if k != l:
                    L -= sum(beta_dot_x) + beta_dot_x[l]

    return -L
...