Точность вычислений в BLAS.h

[Andrey_Kireew]

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

так

Используется синтаксис C++

*F=*F0+ddot(&N,b0,&inc_x,sm,&inc_x);
 

и

Используется синтаксис C++

val=*F0;
    for (j=0; j<NStr; j++)
    val+=*(b0+j)**(sm+j);
    
    *F=val;
 

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

Может кто знает в чём тут дело, и на какие результаты лучше ориентироваться?

[Dmitriy40]

Очевидной причиной является недостаточная точность вычислений с плавающей точкой, в этих циклах ошибки округления накапливаются немного по разному.
Стандартным советом должен бы стать совет увеличить точность промежуточных переменных, но форматов точнее double вроде бы в стандартных библиотеках нет (и уж тем более они не будут ускоряться векторизацией). А с подключением библиотек работы с числами произольной точности резко упадёт скорость вычислений. Можно подключить и просто перепроверить какой именно результат точнее, но это будет адекватно лишь для конкретных матриц, при любом изменении матриц погрешность может и измениться. Я бы верил первому циклу, т.к. он может (хотя и не обязательно) в промежуточных вычислениях использовать более точный формат чисел. Ещё можно переписать весь вычислительный цикл в командах x87 (сопроцессор с плавающей точкой), у него внутренее представление чисел заметно точнее double, для матриц вменяемых размеров должно хватить.
Ещё, раньше я бы посоветовал провести два (три) раза вычисления с округлением в разные стороны и получить минимум и максимум (и вероятное значение если три раза) для каждой переменной на выходе, но как это сейчас реализовать и пользуется ли вообще кто-то этим методом я не в курсе.
Самым правильным методом для работы с любыми матрицами будет ознакомиться с стандартными методами повышения точности решения СЛАУ и их и использовать (кажется там применяются хитрые перестановки элементов матриц в зависимости от конкретных значений чтобы получить как можно меньшую погрешность каждой операции, например сумму произведений вычислять от меньших чисел к большим). Скорее всего это опять заметно замедлит вычисления.
Ну или осуществлять выборочную проверку точности вычислений для каждого "прогона" программы с разными исходными данными в матрицах использованием медленной библиотеки работы с числами произвольной точности. Выборочную чтобы не замедлять все вычисления и в то же время получить примерную оценку погрешности результата.
Или взять в руку карандаш и руками выписать погрешности каждой операции в программе (по ассемблерному листингу!) и потом долго сворачивать их к одному конечному результату. Величина получится огромной, но зато точной.
В общем проще всего оценить реальную погрешность вычислений и смириться с ней. Если никто не предложит более практичного способа.
PS. А, да, проверить расчёт с произвольной точностью можно в Математике (или других системах алгебры), там кажется можно указывать любую точность чисел с плавающей точкой. Будет медленно, зато сильно точнее.

[Pavia]

Andrey_Kireew
Числа в компьютере имеют ограниченное число разрядов, а реальные числа бесконечные. Приходится числа округлять. Отсюда куча различных эффектов.

В плоть от того что от перестановки слагаемых сумма меняется!

Что касается кода не уверен, что компилятор у вас верно обработает операторы. Я бы скобочки поставил.

val+=*(b0+j)*(*(sm+j));

Ориентироваться лучше на проверенные результаты и повторно использовать код. Доверяй, но и проверяй. Поэтому расчёты ведут разными методами, а потом сравнивают результаты.

Советую прочитать:
- про виды погрешностей.
- про числа обусловленности.
- по сравнение чисел с плавающей точкой.

[Andrey_Kireew]

Разобрался самостоятельно:
основная ошибка накапливается при суммировании чисел double с сильно различающимся порядком (как у меня),
насколько я понял самое точное - это сортировать произведения, а потом складывать в порядке возрастания положительные и в порядке уменьшения отрицательные произведения, потом всё складывается. Разумеется это очень долго и обычно неприемлемо.

От части компенсировать погрешность можно добавляя к слагаемым погрешность округления. Так же можно разбить массив на части и суммировать сначала подмассивы а потом полученные результаты.

Лично я попробовал сложить отдельно положительные и отдельно отрицательные части, результаты заметно изменились, причём в сторону функции ddot()

Видимо в blas.h выполняется что то подобное и поэтому результаты будут точнее.

К стати все намёки на СЛАУ и число обусловленности здесь неуместны, т.к. не имеют ни какого отношения к проблеме.

[mihaild]

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

Еще есть суммирование Кахана - работает за линейное время, и обеспечивает довольно хорошую точность.

Andrey_Kireew в сообщении #1229864 писал(а):

Видимо в blas.h выполняется что то подобное и поэтому результаты будут точнее

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

Andrey_Kireew в сообщении #1229864 писал(а):

намёки на СЛАУ и число обусловленности

Число обусловленности для задачи суммирования определено как $\frac{\sum |x_i|}{|\sum x_i|}$ - и имеет примерно тот же смысл, что и для матриц (суммировать наборы с большим числом обусловленности сложнее).

[Dmitriy40]

Т.к. на СЛАУ намекнул я, поясню причину: при решении СЛАУ матричным методом там возникает точно такая же проблема потеря точности при матричных операциях (если не ошибаюсь тоже именно при сложении слишком разновеликих чисел или вычитании слишком близких). Но у меня осталось впечатление что методы оптимизации точности для решения СЛАУ более разработаны чем для матричного умножения, вот поэтому и порекомендовал почитать про них (методы оптимизации точности).
Фактически Вы и сами нашли причину, и совет ниже дали хороший про оптимизацию сложений.
Добавлю лишь что кроме большой ошибки при сложениях может ещё быть большой ошибка при умножениях чисел вида $2^k \pm \varepsilon, \, \varepsilon < 2^{k-25}$, $25$ - довольно примерно, смысл что при умножении таких чисел произведение станет отличаться от $2^{2k}$ в 50-х битах и вся точность потеряется. Например перемножив числа очень близкие к единице можно получить строго единицу, без хвоста ненулевых цифр, и при дальнейшем вычитании/сравнении получить ошибку (например деление на 0), т.е. равными будут числа, на самом деле не одинаковые. Как от этого надёжно избавиться я не в курсе (а предлагать свои скороспелые варианты не хочу), думаю есть хорошо проработанные теоретические методы.

[Andrey_Kireew]

Действительно, так как Вы написали будет ещё точнее ...
У меня работает библиотека libmwblas.lib кажется она на FORTRAN. Бывает ещё cblas - это видимо на С, но у меня она почему то она не слинковалась. Думаю разница там только в интерфейсе. Какого года точно не знаю, я взял её из папки Matlab2013a.

За ссылку спасибо. Это хороший алгоритм. Я его опробовал - результат отличается от ddot() из blas.h, как мне кажется в худшую сторону.

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

На счёт числа обусловленности задачи суммирования, честно говоря, в первый раз слышу. Для матриц оно определяется через матричные нормы (прямой и обратной матриц), а для вектора через что её выразить? и какой смысл это выражение будет иметь?
Ваша формула выглядит странно - непонятно что там суммируется в числителе, а что в знаменателе.

[mihaild]

Формулу поправил. Идея в том, что если итоговый ответ мал по сравнению с суммой модулей слагаемых - то погрешность будет большая. Я не знаю, есть ли у этого числа какой-то хороший аналитический смысл, но как правило погрешность алгоритмов суммирования зависит от него.

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

[Dmitriy40]

mihaild в сообщении #1229875 писал(а):

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

Не могу полностью с этим согласиться: результат проверки будет точен (адекватен) лишь для конкретного примера всех исходных данных и конкретного алгоритма. При любом их изменении надо или верить в получаемую точность, или снова перепроверять.
Потому если точность нужна гарантированная - использовать методы вычислений с заранее известной точностью и всё равно изредка их перепроверять аналитически или методами ещё лучшей точности.
Если точность достаточна наиболее вероятная - это обычный вопрос о контроле качества многостадийного процесса на производстве, выборочные проверки на некоторых примерах и всё такое. Тут не только алгебра, но и статистика. Ну и немного веры что не нарвёшься на случайный выброс.
Это пишу скорее для Andrey_Kireew.

PS. Делать выводы на основе методов примерно одинаковой (и неизвестной) точности - не слишком разумно, как мне кажется. Если не хочется проверять более точными методами, я бы скажем добавил небольшой шум к исходным данным и посмотрел насколько улетит результат. Или вместо double сделал бы все операции single и тоже посмотрел на итог. Да много чего можно придумать.

[Andrey_Kireew]

спасибо, буду думать