Вычислить журнал определителя в TensorFlow, когда определитель переполняется / теряет значение

Моя функция стоимости включает в себя вычисление log(det(A)) (при условии, что det(A) положительно, поэтому журнал имеет смысл, но A не является эрмитовым, так что разложение Холецкого здесь неприменимо). Когда det(A) очень большой/малый, прямой вызов det(A) вызовет переполнение/недостаточное переполнение. Чтобы обойти это, можно использовать математический факт, что

журнал (det (A)) = Tr (журнал (A)),

где последнее можно оценить с помощью разложения LU (что более эффективно, чем собственное значение/SVD). Этот алгоритм был реализован в numpy как numpy.linalg. slogdet, поэтому проблема в том, как вызвать numpy из TensorFlow.


Вот что я пробовал

import numpy as np
import tensorflow as tf
from tensorflow.python.framework import function

def logdet_np(a):
    _, l = np.linalg.slogdet(a)
    return l

def logdet1(a):
    return tf.py_func(logdet_np, [a], tf.float64)

@function.Defun(tf.float64, func_name='LogDet')
def logdet2(a):
    return tf.py_func(logdet_np, [a], tf.float64)

with tf.Session() as sess:
    a = tf.constant(np.eye(500)*10.)
    #print(sess.run(logdet1(a)))
    print(sess.run(logdet2(a)))

Сначала я определяю функцию python для передачи результата numpy. Затем я определил две функции logdet, используя tf.py_func. Вторая функция украшена function.Defun, которая позже используется для определения функций TensorFlow и их градиентов. Тестируя их, я обнаружил, что первая функция logdet1 работает и дает правильный результат. Но вторая функция logdet2 возвращает KeyError.

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-
packages/tensorflow/python/ops/script_ops.py in __call__(self, token, args)
     77   def __call__(self, token, args):
     78     """Calls the registered function for `token` with args."""
---> 79     func = self._funcs[token]
     80     if func is None:
     81       raise ValueError("callback %s is not found" % token)

KeyError: 'pyfunc_0'

Мой вопрос: что не так с декоратором Defun? Почему он конфликтует с py_func? Как я могу правильно обернуть функции numpy в TensorFlor?


Оставшаяся часть определения градиента для logdet связана с вопросом дифференцирование детерминанта матрицы в тензорном потоке. Согласно решению этого вопроса, можно попытаться написать

@function.Defun(tf.float64, tf.float64, func_name='LogDet_Gradient')
def logdet_grad(a, grad):
    a_adj_inv = tf.matrix_inverse(a, adjoint=True)
    out_shape = tf.concat([tf.shape(a)[:-2], [1, 1]], axis=0)
    return tf.reshape(grad, out_shape) * a_adj_inv
@function.Defun(tf.float64, func_name='LogDet', grad_func=logdet_grad)
def logdet(a):
    return tf.py_func(logdet_np, [a], tf.float64, stateful=False, name='LogDet')

Приведенный выше код будет работать, если удастся решить конфликт между Defun и py_func, что является ключевым вопросом, который я поднял выше.


person Everett You    schedule 26.05.2017    source источник
comment
Вы можете использовать SVD и получить определитель из сингулярных значений   -  person Yaroslav Bulatov    schedule 26.05.2017
comment
Является ли A особенным (например, положительно определенным)?   -  person dmuir    schedule 26.05.2017
comment
@dmuir det(A) положителен, но не все собственные значения положительны (это означает, что отрицательные собственные значения появляются парами).   -  person Everett You    schedule 26.05.2017
comment
@YaroslavBulatov, а вы знаете, каков статус реализации градиента для SVD?   -  person Everett You    schedule 26.05.2017
comment
Разделите матрицу на ее наибольшее абсолютное значение.   -  person P-Gn    schedule 26.05.2017
comment
Это связанный с этим вопрос stackoverflow.com/questions/33714832/   -  person Anton Codes    schedule 27.05.2017
comment
@wontonimo Я знаю об этом вопросе, пожалуйста, посмотрите мой обновленный вопрос. Существует конфликт между Defun и py_func, который не позволяет мне реализовать метод в соответствующем вопросе.   -  person Everett You    schedule 30.05.2017


Ответы (3)


С помощью @MaxB я размещаю здесь код для определения функции logdet для log(abs(det(A))) и ее градиента.

  • logdet вызывает функцию numpy numpy.linalg.slogdet для вычисления логарифма определителя с использованием идеи log(det(A))=Tr(log(A)), которая устойчива к переполнению/недостаточности определителя. Он основан на разложении LU, которое более эффективно по сравнению с методом, основанным на собственном значении/SVD.

  • Функция numpy slogdet возвращает кортеж, содержащий как знак определителя, так и журнал (abs (det (A))). Знаком будем пренебрегать, так как он не будет вносить вклад в сигнал градиента при оптимизации.

  • Градиент logdet вычисляется путем обращения матрицы в соответствии с grad log(det(A)) = inv(A)^T. Он основан на коде TensorFlow на _MatrixDeterminantGrad. с небольшими изменениями.


import numpy as np
import tensorflow as tf
# from https://gist.github.com/harpone/3453185b41d8d985356cbe5e57d67342
# Define custom py_func which takes also a grad op as argument:
def py_func(func, inp, Tout, stateful=True, name=None, grad=None):
    # Need to generate a unique name to avoid duplicates:
    rnd_name = 'PyFuncGrad' + str(np.random.randint(0, 1E+8))
    tf.RegisterGradient(rnd_name)(grad)  # see _MySquareGrad for grad example
    g = tf.get_default_graph()
    with g.gradient_override_map({"PyFunc": rnd_name}):
        return tf.py_func(func, inp, Tout, stateful=stateful, name=name)
# from https://github.com/tensorflow/tensorflow/blob/master/tensorflow/python/ops/linalg_grad.py
# Gradient for logdet
def logdet_grad(op, grad):
    a = op.inputs[0]
    a_adj_inv = tf.matrix_inverse(a, adjoint=True)
    out_shape = tf.concat([tf.shape(a)[:-2], [1, 1]], axis=0)
    return tf.reshape(grad, out_shape) * a_adj_inv
# define logdet by calling numpy.linalg.slogdet
def logdet(a, name = None):
    with tf.name_scope(name, 'LogDet', [a]) as name:
        res = py_func(lambda a: np.linalg.slogdet(a)[1], 
                      [a], 
                      tf.float64, 
                      name=name, 
                      grad=logdet_grad) # set the gradient
        return res

Можно проверить, что logdet работает для очень большого/малого определителя, и его градиент также правильный.

i = tf.constant(np.eye(500))
x = tf.Variable(np.array([10.]))
y = logdet(x*i)
dy = tf.gradients(y, [x])
with tf.Session() as sess:
    sess.run(tf.global_variables_initializer())
    print(sess.run([y, dy]))

Результат: [1151.2925464970251, [array([ 50.])]]

person Everett You    schedule 30.05.2017
comment
В своем вопросе вы никогда не говорили, что только для процессора все в порядке (в противном случае вызов NP может быть плохим). Вы также никогда не говорили, что производительность для этой конкретной операции имеет решающее значение (иначе самое простое решение лучше). - person bobcat; 31.05.2017

Если ваша проблема заключается в переполнении, вы можете бороться с ней с помощью простой математики. введите здесь описание изображения

Таким образом, вам просто нужно получить собственные значения, записать их и суммировать.

person Salvador Dali    schedule 27.05.2017
comment
TF не имеет общего собственного значения op (и они не были бы реальными). Использование SVD - правильный подход (как было предложено). - person bobcat; 27.05.2017
comment
@MaxB тот факт, что собственное значение может быть мнимым, не нарушает математику. Вы по-прежнему можете записывать воображаемые числа и, очевидно, можете суммировать несколько воображаемых чисел. Что касается отсутствия собственных значений, я об этом не знаю, но что мешает человеку реализовать для этого собственные функции C++? Сейчас я вижу комментарий про СВД, и если вы считаете, что этого коммента достаточно, чтобы понять, что делать - объясните нам, пожалуйста, связь между СВД и определителем. - person Salvador Dali; 27.05.2017
comment
Написал как ответ, раз уж вы настояли. - person bobcat; 27.05.2017
comment
@SalvadorDali Спасибо. Я знаю математику, которую вы упомянули. Я обнаружил, что разложение по собственным значениям происходит намного медленнее, чем разложение по LU. На самом деле реализация этого подхода Tr(log(A)) уже есть в numpy. Проблема в том, как вызвать numpy в TensorFlow. Я был бы признателен, если бы вы могли взглянуть на мой обновленный вопрос. - person Everett You; 30.05.2017

Вы можете использовать SVD для разложения A:

A = U S V'

Поскольку определитель произведения является произведением определителей, а определители U и V' равны 1 или -1, а определитель S неотрицательный,

abs(det(A)) = det(S)

Следовательно, журнал (положительного) определителя можно вычислить как

tf.reduce_sum(tf.log(svd(A, compute_uv=False)))


Начиная с TF1.1, в tf.svd отсутствует градиент (в будущих версиях он, вероятно, будет), поэтому я предлагаю взять реализацию из кода kofd:

https://github.com/tensorflow/tensorflow/issues/6503

person bobcat    schedule 27.05.2017
comment
Отсутствие градиента для SVD можно обойти тем фактом, что grad of log(det(A)) является транспонированием inv(A). Я действительно протестировал ваш код, он очень медленный из-за SVD, намного медленнее, чем tf.log(tf.matrix_determinant(A)). Я также обнаружил, что numpy предоставляет функцию slogdet для моей цели, которая основана на более эффективной декомпозиции LU. Теперь я хотел бы знать, как вызвать numpy форму TensorFlow. Пожалуйста, взгляните на мой обновленный вопрос. Буду очень признателен вашим пациентам. - person Everett You; 30.05.2017
comment
@EverettYou Вы можете позвонить в NP из TF: см., например. stackoverflow.com /вопросы/43839431/ - person bobcat; 30.05.2017
comment
@EverettYou или, скорее, gist.github.com/harpone/3453185b41d8d985356cbe5e57d67342 - последний комментарий там. - person bobcat; 30.05.2017
comment
Спасибо за указание на эти ссылки. Они очень полезны. - person Everett You; 30.05.2017