Как проверить, совпадают ли две матрицы?

Question

Идея состоит в том, чтобы умножить две матрицы.И сделайте то же умножение, используя Eigen, затем проверьте, является ли результат тем же.

В следующем случае N = 2 возвращает same thing, но N = 1000 возвращает NOT same thing.Почему?

#include <cstdlib>
#include <iostream>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

const int N = 1000;

void mult_matrix(double x[N][N], double y[N][N], double z[N][N]) {
    int rows = N;
    int cols = N;
    for (int i = 0; i < rows; i++)
        for (int j = 0; j < cols; j++)
            for (int k = 0; k < cols; k++)
                z[i][j] += x[i][k] * y[k][j];
}

void check(double *x, double *y, double *z) {

    Matrix<double, Dynamic, Dynamic, RowMajor> m = 
            Matrix<double, Dynamic, Dynamic, RowMajor>::Map(x, N, N) * 
            Matrix<double, Dynamic, Dynamic, RowMajor>::Map(y, N, N);

    cout << m(0, 0) << endl;
    cout << Matrix<double, Dynamic, Dynamic, RowMajor>::Map(z, N, N)(0, 0) << endl;

    if (m == Matrix<double, Dynamic, Dynamic, RowMajor>::Map(z, N, N))
        cout << "same thing" << endl;
    else
        cout << "NOT same thing" << endl;
}

int main() {
    double *a = (double*)malloc(N*N*sizeof(double));
    double *b = (double*)malloc(N*N*sizeof(double));
    double *c = (double*)malloc(N*N*sizeof(double));

    Matrix<double, Dynamic, Dynamic, RowMajor>::Map(a, N, N).setRandom();
    Matrix<double, Dynamic, Dynamic, RowMajor>::Map(b, N, N).setRandom();
    Matrix<double, Dynamic, Dynamic, RowMajor>::Map(c, N, N).setZero();

    mult_matrix((double (*)[N])a, (double (*)[N])b, (double (*)[N])c);
    check(a, b, c);
}

Answer 1

Eigen предоставляет функцию-член isApprox(), которую можно использовать для проверки, равны ли две матрицы в пределах числовой точности.

В вашем коде такое сравнение может бытьпросто достигается путем замены оператора == на isApprox(), например, так:

if (m.isApprox(Matrix<double, Dynamic, Dynamic, RowMajor>::Map(z, N, N)))
  cout << "same thing" << endl;
else
  cout << "NOT same thing" << endl;

Желаемая точность может быть передана как необязательный второй аргумент isApprox().

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

Answer 2

Не ответ, но слишком длинный комментарий.

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

Из-за конечности представления, то есть ограниченного числа значащих цифр, неизбежно возникают ошибки усечения.Например, 0.1 будет представлен не точно как 0.1, а как что-то вроде 0.99999999999993 или 0.100000000000002 (фиктивные значения).Точное значение может зависеть от конкретного базового алгоритма преобразования и политики округления.

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

|a - b| < max(|a|, |b|).ε

, где ε близко к точности станка (около 10^-7 для одинарной точности).

Но в несчастливых случаях, таких как, например, численная оценка производных, возникает явление, известное как катастрофическое аннулирование : при вычитании двух близких значений число точных значащих цифр резко падает.

Например, (sin(1.000001) - sin(1)) / 0.000001 = (0.84147104 - 0.84147100) / 0.0000001 = 0.40000000, в то время как точное значение должно быть 0.5403023.

. И катастрофическое сокращение происходит в матричных произведениях, когда два вектора близки к перпендикуляру.Тогда критерий относительной ошибки больше не работает.

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

Answer 3

От Golub & Van Loan анализ ошибок по точечным произведениям дает следующую оценку:

Let u = 2^-t (t - количество битов в мантиссах) и n - числокомпоненты в строках / столбцах.Тогда с допущением n u < 0.01 (что легко выполняется) ошибка усечения для точечного произведения xTy ограничена

1.01 n u |x|T |y|

(последний фактор - это точечное произведениевекторы, когда вы отбрасываете все отрицательные знаки).

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

---

Конечное примечание: когда xTy = 0, относительная ошибка стремится к бесконечности.

Answer 4

Из-за ошибок округления при сравнении чисел с плавающей запятой с оператором == возможны ошибки.Так что для N=2, это может работать, но для больших N это, скорее всего, не получится.

Попробуйте следующий компаратор вместо ==:

bool double_equals(double a, double b, double epsilon = 0.001)
{
    return std::abs(a - b) < epsilon;
}

после комментариевниже @Jonathan Leffler, вышеупомянутый компаратор не идеален, так как лучше использовать относительную разницу, чем абсолютную разницу.

double reldiff(double a, double b) {
    double divisor = fmax(fabs(a), fabs(b)); /* If divisor is zero, both x and y are zero, so the difference between them is zero */
    if (divisor == 0.0) return 0.0; return fabs(a - b) / divisor; 
}

bool double_equals(double a, double b, double rel_diff)
{
    return reldiff(a, b) < rel_diff;
}