Ошибка умножения собственных разреженных матриц с малым числом

Question

Я использую библиотеку C ++ Eigen 3 в своей программе.В частности, мне нужно умножить две собственные разреженные матрицы и сохранить результат в другую собственную разреженную матрицу.Тем не менее, я заметил, что если некоторые записи в собственной разреженной матрице меньше 1e-13, соответствующая запись в результате будет 0 вместо небольшого числа.Скажем, я умножаю разреженную единичную матрицу a и другую разреженную матрицу b.Если максимальная запись b, то есть b (0,0) меньше, чем 1e-13, например 9e-14, верхняя запись результата c = a * b, то есть c (0,0),0 вместо 9e-14.

Вот код, который я тестирую,

#include <iostream>
#include <math.h>
#include <Eigen/Sparse>
#include <Eigen/Dense>
#include <Eigen/Eigen>
#include <Eigen/Dense>
#include <Eigen/LU>

using namespace std;
using namespace Eigen;

int main() {

DynamicSparseMatrix<double, RowMajor> a(2,2);
DynamicSparseMatrix<double, ColMajor> b(2,2);
DynamicSparseMatrix<double, RowMajor> c(2,2);
a.coeffRef(0,0) = 1;
a.coeffRef(1,1) = 1;
b.coeffRef(0,0) = 9e-14;
b.coeffRef(1,1) = 1;
cout << "a" << endl;
cout << a << endl;
cout << "b" << endl;
cout << b << endl;
c = a * b;
cout << "c" << endl;
cout << c << endl;
Matrix<double, Dynamic, Dynamic> a2(2,2);
Matrix<double, Dynamic, Dynamic> b2(2,2);
Matrix<double, Dynamic, Dynamic> c2(2,2);
a2(0,0) = 1;
a2(1,1) = 1;
b2(0,0) = 9e-14;
b2(1,1) = 1;
cout << "a2" << endl;
cout << a2 << endl;
cout << "b2" << endl;
cout << b2 << endl;
c2 = a2 * b2;
cout << "c2" << endl;
cout << c2 << endl;

return 1;
}

Вот странный вывод

a

1 0

0 1

b

Ненулевые записи:

(9e-14,0) (1,1)

Указатели на столбцы:

0 1 $

9e-14 0

0 1

c

0 0

01

a2

1 0

0 1

b2

9e-14 0

0 1

c2

9e-14 0

0 1

Вы можете видеть, что умножение плотных матриц хорошо, но результат разреженных матрицнеправильно в верхней левой записи, и b имеет странный выходной формат.

Я отладил исходный код Eigen, но не смог найти, где два числа умножены в матрице.Есть идеи?

Answer 1

Есть несколько причин обрезать небольшие значения до нуля в библиотеке линейной алгебры.

Когда вы приблизитесь к нулю, формат с плавающей запятой не сможет хорошо представить вычисления. Например, 9e-14 + 1.0 может равняться 1.0, потому что недостаточно точности для представления истинного результата.

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

Наконец, это помогает сохранить разреженность матриц. Если при вычислении создаются очень маленькие ненулевые значения, вы можете их обрезать и оставить матрицу разреженной. Во многих физических приложениях, таких как анализ методом конечных элементов, небольшие значения не являются физически значимыми, и их можно безопасно удалить.

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

Общая картина такова: если ваш расчет основан на значениях с плавающей запятой, очень близких к нулю, вам следует выбрать другой численный метод. Выход никогда не будет надежным. Например, вы можете заменить трехмерные вращения матриц кватернионными вращениями. Эти методы называются «численно устойчивыми», и они являются основным направлением численного анализа.

Answer 2

Я не уверен точно, где и как выполняется усечение в Eigen, но значение эпсилона, используемое для усечения, определено в NumTraits< Scalar >::dummy_precision() в Eigen/src/Core/NumTraits.h.

Answer 3

Я знаю, что этот вопрос сейчас очень старый.Но для дальнейшего использования: DynamicSparseMatrix устарело, вместо этого следует использовать обычный SparseMatrix (в несжатом режиме, если это необходимо).

Кроме того, продукты с разреженной матрицей больше не будут автоматически сокращать результаты [ 1 ].Если кто-то на самом деле хочет сократить результат произведения с разреженной матрицей, нужно написать это в явном виде, например:

C = A*B;                     // don't suppress numerical zeros
C = (A*B).pruned();          // suppress numerical zeros (exact)
C = (A*B).pruned(ref, eps);  // suppress anything smaller than ref*eps

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