想法是乘以两个矩阵。并使用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);
}
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()
。
正如评论中所讨论的那样,可能总会存在这样的比较可能无法可靠地工作的情况。但使用像Qazxswpoi或isApprox()
这样的Eigen函数比任何简单的手工制作解决方案都更有效。
由于舍入误差,将浮点数与isMuchSmallerThan()
运算符进行比较会出错。所以对于==
来说,它可能有用,但对于大型的N=2
,它很可能会失败。
尝试使用以下比较器而不是N
:
==
以下评论由@Jonathan Leffler提出,上述比较器并不理想,因为使用相对差异比绝对差异更好。
bool double_equals(double a, double b, double epsilon = 0.001)
{
return std::abs(a - b) < epsilon;
}
评论不是答案,但答案太长。
坏消息是没有办法比较两个浮点值的相等性。
由于表示的有限性,即有限数量的有效数字,不可避免地发生截断误差。例如,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;
}
不会像0.1
那样完全表示,而是像0.1
或0.99999999999993
(虚构值)。确切的值可以取决于特定的基本转换算法和舍入策略。
在幸运的情况下,当你进行计算时,trunction错误会轻微累积,因此有效数字的数量也会轻微减少。因此,使用有界相对误差测试相等性是有意义的,例如:
0.100000000000002
其中|a - b| < max(|a|, |b|).ε
接近机器精度(单精度约为ε
)。
但是在不幸的情况下,例如衍生物的数值评估,会发生称为灾难性消除的现象:当您减去两个附近的值时,精确有效数字的数量会急剧下降。
例如10^-7
,而精确值应该是(sin(1.000001) - sin(1)) / 0.000001 = (0.84147104 - 0.84147100) / 0.0000001 = 0.40000000
。
当两个向量接近垂直时,矩阵产品确实会发生灾难性的消除。然后相对误差标准不再起作用。
最糟糕的情况是当你想要检查一个数量为零时,例如在查找函数的根时:“接近零”值可以具有任意数量级(想想当变量是相同问题的解决方案时)用米表示,然后以毫米表示。没有相对误差标准可以工作,但即使没有绝对误差也可以工作,除非你有关于幅度的额外信息。没有通用的比较器可以工作。
从Golub和Van Loan,对点积的误差分析给出了以下估计:
设0.5403023
(u = 2^-t
是尾数中的位数),t
是行/列中的组件数。然后假设n
(容易持有),点积n u < 0.01
上的截断误差受限于
xTy
(当你放弃所有负号时,最后一个因素是向量的点积)。
这为您提供了一种可靠的方法来为产品矩阵的元素设置精度标准。
最后注意:当1.01 n u |x|T |y|
时,相对误差趋于无穷大。