Точность в цифрах: проблемы при сравнении чисел

Question

Сначала немного фона.Я нахожу собственные значения и собственные векторы реальной симметричной матрицы, в которой строки суммируются с 0. Более конкретно, когда я нахожу собственный вектор, я использую $ argsort $, чтобы найти перестановку, которая сортирует одно из собственных значений, и применить перестановку кисходная матрица.

Теперь я реализовал код на python, используя пакет numpy.Сам код является рекурсивным, и если он находит набор значений в собственном векторе, которые равны, он извлекает симметричную подматрицу, соответствующую индексам, для которых у нас равные значения, и снова применяет алгоритм к этой матрице.

Хотя все это очень хорошо, и в основном это грубая работа, я был удивлен, когда группа индексов, которые должны были бы соответствовать равным элементам в собственном векторе, не были признаны имеющими равные значения.Проблема заключалась в том, что значения вычислялись с точностью до машины по некоторому алгоритму (возможно, Ланцошу, но я не совсем знаком с numpy).Это пример вывода, в котором я явно проверяю разницу между двумя записями в собственном векторе:

    >>> T=spectral.seriation(A,index)

    columns [ 0  1  2  3  4  5  6  7  8  9 10 11]

    [  3.30289130e-01  -2.75240941e-01  -2.75240941e-01   3.30289130e-01
    -2.75240941e-01   3.30289130e-01  -2.75240941e-01   3.30289130e-01
    3.30289130e-01  -2.75240941e-01  -1.69794463e-16  -2.75240941e-01]

    [ 4  6  9  1  2 11 10  0  5  7  8  3]

    difference   -5.55111512313e-17

Подпрограмма seriation () является рекурсивной функцией.Массив чисел с плавающей запятой - это рассматриваемый собственный вектор, а массив ниже, который дает отсортированный порядок столбцов.Обратите внимание, что столбцы [4,6,9,1,2,11] имеют одинаковое значение.Тем не менее, вычисления собственного вектора и собственного значения всегда являются приблизительными, и, действительно, когда я вывожу разницу между записью в столбце 9 и столбце 2, она не равна нулю.Там, где алгоритм должен группировать [4,6,9,1,2,11], он группирует только [4,6,9] и помещает остальное в другую группу, бросая гаечный ключ в работы.

Таким образом, вопрос заключается в следующем: существует ли метод для выполнения вычислений произвольной точности в numpy?Если это не удастся, что будет «хорошим» решением этой проблемы?

Кроме того, я, вероятно, должен упомянуть, что математически доказано, что эти записи должны быть равны.Это свойство матрицы, но, надеюсь, это не относится к вопросу.

Answer 1

Двойные числа не являются точно действительными числами [даже не рациональными].Существует бесконечное число рациональных чисел в каждом диапазоне [ну, каждый диапазон, по крайней мере, с двумя элементами, если быть точным], но только конечное число бит для их представления.
Таким образом, вы должны ожидать некоторых ошибок округления для «точных» вычислений.

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

Answer 2

Проверьте numpy.allclose и numpy.isclose функции для проверки равенства в пределах допуска.

Answer 3

При выполнении вычитания двух чисел с плавающей запятой сопоставимого размера точность не должна быть проблемой, т. Е. Если [2] и [9] действительно одинаковы, тогда разница будет равна нулю.

Я подозреваю, что на самом деле дело в том, что по умолчанию выходные данные отображают числа с 8 десятичными разрядами, но кроме того, числа отличаются, как правило, двойное число имеет около 16 десятичных знаков точности (чтобы выяснить, как работает numpy.finfo(numpy.float).eps для получите эпсилон машины, который дает наименьшее возможное число выше 1)

Попробуйте проверить числа, используя выходной формат "%.16f\n%.16f" % myarray[[2, 9]].

Если они отличаются, но вы довольны 7d.p сходства, то вы можете усечь результаты, используя что-то вроде numpy.around(differences, 7).

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

sigcnd, expn = numpy.frexp(myarray)
sigcnd = numpy.around(sigcnd, 7)
truncated_myarray = numpy.ldexp(sigcnd, expn)

Answer 4

Если вы хотите, чтобы индексы почти равных элементов соответствовали заданному допуску, вы можете сделать что-то вроде:

def almost_matches(x, array, rtol=1e-05, atol=1e-08):
    answer = []
    for y in xrange(len(array)):
        if abs(x-array[y]) <= (atol + rtol * abs(array[y])):
            answer.append(y)
    return answer

(используя то же приблизительное сравнение, что и numpy.allclose () использует)

>>> a = [3.30289130e-01,  -2.75240941e-01,  -2.75240941e-01,   3.30289130e-01, -2.75240941e-01, 3.30289130e-01,  -2.75240941e-01,   3.30289130e-01, 3.30289130e-01,  -2.75240941e-01,  -1.69794463e-16,  -2.75240941e-01]
>>> almost_matches(min(a), a)
[1, 2, 4, 6, 9, 11]