用 2 个平方矩阵模拟 matlab 的 mrdivide
Simulating matlab's mrdivide with 2 square matrices
我有2个19x19的方阵(a &b)和我试图使用斜杠(mrdivide)运算符执行除法，使
c = a / b
我正试图在OpenCV中实现这一点。我发现一些人建议使用cv::solve，但到目前为止，我一直无法找到任何能给我一个接近matlab的结果。
有没有人知道我如何实现mrdivide与opencv?
我试过下面的代码:
cv::Mat mldivide(const cv::Mat& A, const cv::Mat& B )
{
//return b * A.inv();
cv::Mat a;
cv::Mat b;
A.convertTo( a, CV_64FC1 );
B.convertTo( b, CV_64FC1 );
cv::Mat ret;
cv::solve( a, b, ret, cv::DECOMP_NORMAL );
cv::Mat ret2;
ret.convertTo( ret2, A.type() );
return ret2;
}
然后我实现了mrdivide如下:
cv::Mat mrdivide(const cv::Mat& A, const cv::Mat& B )
{
return mldivide( A.t(), B.t() ).t();
}
(编辑:根据答案，当我正确使用它时，它实际上给了我正确的答案!)
这给了我一个错误的答案，即不像matlab。根据评论，我也试过
cv::Mat mrdivide(const cv::Mat& A, const cv::Mat& B )
{
return A * B.inv();
}
这个答案和上面的一样，但也是错误的。
你不应该用inv来解Ax=b或xA=b方程。虽然这两个方法在数学上是等价的(x=solve(A,b)和x=inv(A)*B)，但在处理浮点数时，这是完全不同的事情!http://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/
作为一般规则，永远不会乘以矩阵逆． 相反，对于一次性系统使用正/反斜线运算符(或等效的"solve"方法)，或者当您希望将相同的A与多个b重用时显式执行矩阵分解(例如LU, QR, Cholesky等)
让我举一个具体的例子来说明反转的问题。我将使用MATLAB和mexopencv，一个允许我们直接从MATLAB调用OpenCV的库。
(这个例子是从Tim Davis提交的优秀FEX中借用的，他也是开发SuiteSparse的人。我展示的是左除Ax=b的情况，但右除xA=b也是如此。
让我们首先为Ax=b系统构建一些矩阵:
% Ax = b
N = 16; % square matrix dimensions
x0 = ones(N,1); % true solution
A = gallery('frank',N); % matrix with ill-conditioned eigenvalues
b = A*x0; % Ax=b system
下面是16x16矩阵A和16x1向量b的样子(请注意，真正的解x0只是一个1的向量):
A = b =
16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1 136
15 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1 135
0 14 14 13 12 11 10 9 8 7 6 5 4 3 2 1 119
0 0 13 13 12 11 10 9 8 7 6 5 4 3 2 1 104
0 0 0 12 12 11 10 9 8 7 6 5 4 3 2 1 90
0 0 0 0 11 11 10 9 8 7 6 5 4 3 2 1 77
0 0 0 0 0 10 10 9 8 7 6 5 4 3 2 1 65
0 0 0 0 0 0 9 9 8 7 6 5 4 3 2 1 54
0 0 0 0 0 0 0 8 8 7 6 5 4 3 2 1 44
0 0 0 0 0 0 0 0 7 7 6 5 4 3 2 1 35
0 0 0 0 0 0 0 0 0 6 6 5 4 3 2 1 27
0 0 0 0 0 0 0 0 0 0 5 5 4 3 2 1 20
0 0 0 0 0 0 0 0 0 0 0 4 4 3 2 1 14
0 0 0 0 0 0 0 0 0 0 0 0 3 3 2 1 9
0 0 0 0 0 0 0 0 0 0 0 0 0 2 2 1 5
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 2
现在让我们通过找到解决方案并使用NORM函数(或者cv::norm，如果你愿意)计算残差来比较cv::invert和cv::solve:
% inverting (OpenCV)
x1 = cv.invert(A)*b;
r1 = norm(A*x1-b)
% inverting (MATLAB)
x2 = inv(A)*b;
r2 = norm(A*x2-b)
% solve using matrix factorization (OpenCV)
x3 = cv.solve(A,b);
r3 = norm(A*x3-b)
% solve using matrix factorization (MATLAB)
x4 = Ab;
r4 = norm(A*x4-b)
下面是找到的解决方案(我减去1，这样你就可以看到它们离真正的解决方案x0有多远):
>> format short g
>> [x1 x2 x3 x4] - 1
ans =
9.0258e-06 3.1086e-15 -1.1102e-16 2.2204e-16
-0.0011101 -1.0181e-13 -2.2204e-15 -2.3315e-15
-0.0016212 -2.5123e-12 3.3751e-14 3.3307e-14
0.0037279 4.1745e-11 -4.3476e-13 -4.3487e-13
-0.0022119 4.6216e-10 5.2165e-12 5.216e-12
-0.0010476 1.3224e-09 -5.7384e-11 -5.7384e-11
0.0035461 2.2614e-08 5.7384e-10 5.7384e-10
-0.0040074 -4.1533e-07 -5.1646e-09 -5.1645e-09
0.0036477 -4.772e-06 4.1316e-08 4.1316e-08
-0.0033358 4.7499e-06 -2.8922e-07 -2.8921e-07
0.0059112 -0.00010352 1.7353e-06 1.7353e-06
-0.0043586 0.00044539 -8.6765e-06 -8.6764e-06
0.0069238 -0.0024718 3.4706e-05 3.4706e-05
-0.0019642 -0.0079952 -0.00010412 -0.00010412
0.0039284 0.01599 0.00020824 0.00020823
-0.0039284 -0.01599 -0.00020824 -0.00020823
最重要的是，以下是每种方法的错误:
r1 =
0.1064
r2 =
0.060614
r3 =
1.4321e-14
r4 =
1.7764e-15
后两个是更精确的数量级，甚至不接近!这只是一个16变量的方程组。反演在数值上的可靠性较差，特别是当矩阵很大且稀疏时…
现在回答你的问题，你使用cv::solve的想法是正确的，但你只是在右除法的情况下得到了错误的操作数顺序。
在MATLAB中，运算符/和(或mrdivide和mldivide)通过等式B/A = (A'B')'相互关联(这是转置性质的简单结果)。
所以对于OpenCV函数，你会写(注意A和b的顺序):
% Ax = b
x = cv.solve(A, b); % Ab or mldivide(A,b)
% xA = b
x = cv.solve(A', b')'; % b/A or mrdivide(b,A)
OpenCV公开的API在这里有点尴尬，所以我们必须做所有这些转置。事实上，如果您参考等效的LAPACK例程(例如DGESV或DGESVX)，它们实际上允许您指定矩阵是否在TRANS=T或TRANS=N中进行转置(在该级别上，转置实际上只是不同的内存布局，C或Fortran排序)。例如，MATLAB提供了linsolve函数，它允许您在选项中指定这些类型的东西…
(顺便说一句，在c++ OpenCV中编码时，我更喜欢使用函数形式的操作，如cv::transpose，而不是矩阵表达式变体，如Mat::t。前者可以就地操作，而后者会创建不必要的临时副本)。
现在，如果你正在寻找一个性能良好的c++线性代数实现，可以考虑使用Eigen(它甚至可以很好地与OpenCV集成)。另外，它是一个纯基于模板的库，所以不需要担心链接或二进制文件，只需要包含头文件。
编辑(回复评论)
@Goz:
查找返回值优化。"不必要的临时副本"不存在
我知道RVO和move语义，但它在这里不是很重要;cv::Mat类是复制友好的，有点像引用计数的智能指针。这意味着当按值传递时，它只执行具有数据共享的浅拷贝。为新副本创建的唯一部分是mat头中的部分，这些部分在大小方面无关紧要(存储诸如维度/通道数，步长和数据类型等内容)。
我说的是一个显式的深度复制，而不是你从函数调用返回时考虑的那个…
感谢你的评论，它让我有动力去真正挖掘OpenCV源代码，这不是最容易阅读的…代码几乎没有注释，有时很难理解。看到OpenCV真正关心性能，复杂性是可以理解的，实际上令人印象深刻的是，许多功能以各种方式实现(常规CPU实现，循环展开版本，SIMD矢量化版本(SSE, AVX, NEON等)，使用各种后端的并行和线程版本，来自英特尔IPP的优化实现，使用OpenCL或CUDA的GPU加速版本，Tegra, OpenVX的移动加速版本等)
让我们以下面的例子为例，跟踪我们的步骤:
Mat A = ..., b = ..., x;
cv::solve(A.t(), b, x);
，函数定义如下:
bool cv::solve(InputArray _src, InputArray _src2arg, OutputArray _dst, int method)
{
Mat src = _src.getMat(), _src2 = _src2arg.getMat();
_dst.create( src.cols, _src2.cols, src.type() );
Mat dst = _dst.getMat();
...
}
现在我们必须找出中间的步骤。首先是t成员方法:
MatExpr Mat::t() const
{
MatExpr e;
MatOp_T::makeExpr(e, *this);
return e;
}
返回MatExpr，这是一个允许矩阵表达式延迟求值的类。换句话说，它不会立即执行转置，而是存储对原始矩阵的引用和最终要对其执行的操作(转置)，但是它会在绝对必要时(例如，当赋值或强制转换为cv::Mat时)才对其进行求值。
接下来让我们看看相关部分的定义。注意，在实际的代码中，这些东西被拆分到许多文件中。为了方便阅读，我只把有趣的部分拼凑在这里，但这远远不是完整的:
class MatExpr
{
public:
MatExpr()
: op(0), flags(0), a(Mat()), b(Mat()), c(Mat()), alpha(0), beta(0), s()
{}
explicit MatExpr(const Mat& m)
: op(&g_MatOp_Identity), flags(0), a(m), b(Mat()), c(Mat()),
alpha(1), beta(0), s(Scalar())
{}
MatExpr(const MatOp* _op, int _flags, const Mat& _a = Mat(),
const Mat& _b = Mat(), const Mat& _c = Mat(),
double _alpha = 1, double _beta = 1, const Scalar& _s = Scalar())
: op(_op), flags(_flags), a(_a), b(_b), c(_c), alpha(_alpha), beta(_beta), s(_s)
{}
MatExpr t() const
{
MatExpr e;
op->transpose(*this, e);
return e;
}
MatExpr inv(int method) const
{
MatExpr e;
op->invert(*this, method, e);
return e;
}
operator Mat() const
{
Mat m;
op->assign(*this, m);
return m;
}
public:
const MatOp* op;
int flags;
Mat a, b, c;
double alpha, beta;
Scalar s;
}
Mat& Mat::operator = (const MatExpr& e)
{
e.op->assign(e, *this);
return *this;
}
MatExpr operator * (const MatExpr& e1, const MatExpr& e2)
{
MatExpr en;
e1.op->matmul(e1, e2, en);
return en;
}
a中(同样，cv::Mat实例将共享数据，因此不会复制)，以及执行op的操作，以及其他一些对我们不重要的事情。
这是矩阵操作类MatOp，以及它的一些子类(我只展示转置和逆操作，但还有更多):
class MatOp
{
public:
MatOp();
virtual ~MatOp();
virtual void assign(const MatExpr& expr, Mat& m, int type=-1) const = 0;
virtual void transpose(const MatExpr& expr, MatExpr& res) const
{
Mat m;
expr.op->assign(expr, m);
MatOp_T::makeExpr(res, m, 1);
}
virtual void invert(const MatExpr& expr, int method, MatExpr& res) const
{
Mat m;
expr.op->assign(expr, m);
MatOp_Invert::makeExpr(res, method, m);
}
}
class MatOp_T : public MatOp
{
public:
MatOp_T() {}
virtual ~MatOp_T() {}
void assign(const MatExpr& expr, Mat& m, int type=-1) const
{
Mat temp, &dst = _type == -1 || _type == e.a.type() ? m : temp;
cv::transpose(e.a, dst);
if( dst.data != m.data || e.alpha != 1 ) dst.convertTo(m, _type, e.alpha);
}
void transpose(const MatExpr& e, MatExpr& res) const
{
if( e.alpha == 1 )
MatOp_Identity::makeExpr(res, e.a);
else
MatOp_AddEx::makeExpr(res, e.a, Mat(), e.alpha, 0);
}
static void makeExpr(MatExpr& res, const Mat& a, double alpha=1)
{
res = MatExpr(&g_MatOp_T, 0, a, Mat(), Mat(), alpha, 0);
}
};
class MatOp_Invert : public MatOp
{
public:
MatOp_Invert() {}
virtual ~MatOp_Invert() {}
void assign(const MatExpr& e, Mat& m, int _type=-1) const
{
Mat temp, &dst = _type == -1 || _type == e.a.type() ? m : temp;
cv::invert(e.a, dst, e.flags);
if( dst.data != m.data ) dst.convertTo(m, _type);
}
void matmul(const MatExpr& e1, const MatExpr& e2, MatExpr& res) const
{
if( isInv(e1) && isIdentity(e2) )
MatOp_Solve::makeExpr(res, e1.flags, e1.a, e2.a);
else if( this == e2.op )
MatOp::matmul(e1, e2, res);
else
e2.op->matmul(e1, e2, res);
}
static void makeExpr(MatExpr& res, int method, const Mat& m)
{
res = MatExpr(&g_MatOp_Invert, method, m, Mat(), Mat(), 1, 0);
}
};
static MatOp_Identity g_MatOp_Identity;
static MatOp_T g_MatOp_T;
static MatOp_Invert g_MatOp_Invert;
OpenCV大量使用操作符重载，所以各种操作，如A+B,A-B,A*B，…实际映射到相应的矩阵表达式操作
InputArray。它基本上存储了一个void*指针以及关于传递的东西的信息(它是什么类型:Mat,MatExpr,Matx,vector<T>,UMat等)，这样它就知道如何在请求时将指针转换回InputArray::getMat:
typedef const _InputArray& InputArray;
class _InputArray
{
public:
_InputArray(const MatExpr& expr)
{ init(FIXED_TYPE + FIXED_SIZE + EXPR + ACCESS_READ, &expr); }
void init(int _flags, const void* _obj)
{ flags = _flags; obj = (void*)_obj; }
Mat getMat_(int i) const
{
int k = kind();
int accessFlags = flags & ACCESS_MASK;
...
if( k == EXPR ) {
CV_Assert( i < 0 );
return (Mat)*((const MatExpr*)obj);
}
...
return Mat();
}
protected:
int flags;
void* obj;
Size sz;
}
现在我们看到Mat::t如何创建并返回MatExpr实例。然后被cv::solve接收为InputArray。现在，当它调用InputArray::getMat来检索矩阵时，它有效地将存储的MatExpr转换为Mat，调用强制转换操作符:
MatExpr::operator Mat() const
{
Mat m;
op->assign(*this, m);
return m;
}
，所以它声明了一个新的矩阵m，调用MatOp_T::assign与新的矩阵作为目标。反过来，这迫使它通过最终调用cv::transpose来求值。它计算转置结果到这个新矩阵作为目标。
所以我们最终有两个副本，原始的A和转置的A.t()返回。
现在，将它与:
进行比较Mat A = ..., b = ..., x;
cv::transpose(A, A);
cv::solve(A, b, x);
在这种情况下，A在位置上调换了，并且具有更少的抽象级别。
现在我展示所有这些的原因并不是为了争论这一个额外的副本，毕竟这不是什么大不了的事情:)我发现的真正整洁的事情是，以下两个表达式不做同样的事情，并给出不同的结果(我不是在谈论逆是否到位):
Mat A = ..., b = ..., x;
cv::invert(A,A);
x = A*b;
Mat A = ..., b = ..., x;
x = inv(A)*b;
事实证明，第二个实际上足够聪明，可以调用cv::solve(A,b)!如果你回到MatOp_Invert::matmul(当一个惰性反转稍后与另一个惰性矩阵乘法链接时调用)。
void MatOp_Invert::matmul(const MatExpr& e1, const MatExpr& e2, MatExpr& res) const
{
if( isInv(e1) && isIdentity(e2) )
MatOp_Solve::makeExpr(res, e1.flags, e1.a, e2.a);
...
}
检查表达式inv(A)*B中的第一个操作数是否为反转操作，第二个操作数是否为单位操作(即普通矩阵，而不是另一个复杂表达式的结果)。在这种情况下，它将存储操作更改为延迟求解操作MatOp_Solve(类似地，它是cv::solve函数的包装器)。在我看来，这很聪明!即使你写了inv(A)*b，它也不会真正计算逆，相反，它理解通过使用矩阵分解来解决这个系统更好地重写它。
不幸的是，这将只有利于形式为inv(A)*b的表达式，而不是b*inv(A)的表达式(这将最终计算逆，这不是我们想要的)。因此，在您解决xA=b的情况下，您应该坚持显式调用cv::solve…
当然，这只适用于在c++中编码时(多亏了操作符重载和惰性表达式的魔力)。如果你使用OpenCV从另一种语言使用一些包装器(如Python, Java, MATLAB)，你可能不会得到任何，应该明确使用cv::solve，就像我在以前的MATLAB代码中所做的那样，对于Ax=b和xA=b两种情况。
希望这对你有帮助，很抱歉写了这么长一篇文章;)
在MATLAB中，对两个维数相容的矩阵使用mrdivide，使得a / b等价于a * b^{-1}，其中b^{-1}是b的逆。因此，您可以做的也许是先对矩阵b求反，然后再对这个a进行预乘。
一种方法是在矩阵b上使用cv::invert，然后与a预乘。这可以通过下面的函数定义来完成(借用你在上面的帖子中的代码):
cv::Mat mrdivide(const cv::Mat& A, const cv::Mat& B)
{
cv::Mat bInvert;
cv::invert(B, bInvert);
return A * bInvert;
}
另一种方法是使用cv::Mat接口内置的inv()方法并使用该方法将矩阵本身相乘:
cv::Mat mrdivide(const cv::Mat& A, const cv::Mat& B)
{
return A * B.inv();
}
我不确定哪一个更快，所以你可能要做一些测试，但这两种方法中的任何一种都应该工作。然而，为了提供一些关于可能时间的见解，有三种方法可以在OpenCV中反转矩阵。您只需将第三个参数重写为cv::invert或在cv::Mat.inv()中指定该方法即可。
这篇StackOverflow文章介绍了使用三种方法对一个相对较大的矩阵进行逆运算的时间:
- 在C++STL中是否有Polyval(Matlab函数)等价物？
- 有可能在Armadillo中复制MATLAB circshift方法吗
- 使用 MATLAB 编码器生成C++代码：编译错误"undefined reference to `rgb2gray_tbb_real64'"
- MATLAB to C++: csvread() not supported by MATLAB Coder
- 加载由 MATLAB Coder 生成的带有函数的 DLL，该函数调用外部函数
- MATLAB:跟踪imufilter对象中的状态变化
- 如何将 c++ 中的客户端 TCP 的替身列表发送到 Matlab 中的 TCP 服务器?
- 将三角函数的正确值与matlab进行比较
- 连接 MATLAB 和 Visual Studios 的问题
- 如何将C++连接到 Matlab
- 如何发送 Mat H=findHomography 返回 Matlab
- 如何在C++中编写 MATLAB fix(X） 函数?
- 提供变量作为 MATLAB 系统命令的输入参数，以便C++可执行文件
- 有没有更好的方法可以使用特征/C++实现 matlab 的逻辑索引？
- 当我运行MEX文件时，MATLAB崩溃
- 从 MATLAB 到 C++：相当于带有选项 'remove' 的 bwmorph
- MATLAB的imfinfo in openCV
- 将C++中的多个参数传递给MatLab共享库函数
- 将 matlab 数组(MAT-文件）转换为C++数组
- 用 2 个平方矩阵模拟 matlab 的 mrdivide
