用2个方阵模拟matlab的mrdivide
Simulating matlab's mrdivide with 2 square matrices
我有 2 个 19x19 方阵 (a & b)，我正在尝试使用斜杠 (mrdivide) 运算符执行除法，使得
c = a / b
我正在尝试在 OpenCV 中实现它。我发现有几个人建议使用 cv::solve，但到目前为止，我一直无法找到任何能给我提供接近 matlab 的结果的东西。
有人知道我如何用 opencv 实现 mrdivide 吗？
我试过以下代码：
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();
}
这给出了与上面相同的答案，但也是错误的。
在 MATLAB 中，在两个兼容维度的矩阵上使用 mrdivide，使得 a / b 等同于 a * b^{-1}，其中 b^{-1} 是 b 的倒数.因此，您可以做的可能是先反转矩阵 b，然后预乘此 a。
一种方法是在矩阵b上使用cv::invert，然后预乘a。这可以通过以下函数定义来完成（借用上面 post 中的代码）：
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() 中指定方法即可。
这个 Whosebug post 使用三种方法对相对较大的矩阵求逆矩阵进行计时：Fastest method in inverse of matrix
您不应该使用 inv 来求解 Ax=b 或 xA=b 方程。虽然这两种方法在数学上是等效的（x=solve(A,b) 和 x=inv(A)*B），但在处理 floating-point 数字时完全不同！
http://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/
作为一般规则，永远不要乘以逆矩阵。而是对 one-off 系统使用 forward/backward-slash 运算符（或等效的 "solve" 方法），或者在需要时显式执行矩阵分解（想想 LU、QR、Cholesky 等）对多个 b 的
重复使用相同的 A
让我举一个具体的例子来说明反转的问题。我将使用 MATLAB 和 mexopencv，一个允许我们直接从 MATLAB 调用 OpenCV 的库。
（这个例子是蒂姆戴维斯从 this excellent FEX submission 那里借来的，他是 SuiteSparse 背后的同一个人。我展示的是 left-division Ax=b 的情况，但是同样适用于 right-division 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
现在让我们比较cv::invert against cv::solve by finding the solution and computing the residual error using the NORM function (or cv::norm，如果你愿意的话）：
% 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 = A\b;
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 的想法是正确的，但是在 right-division.
的情况下你只是把操作数的顺序弄错了
在 MATLAB 中，运算符 / 和 \（或 mrdivide and mldivide) are related to each other by the equation B/A = (A'\B')' (this is a simple result of transpose properties）。
因此对于 OpenCV 函数，您将编写（注意 A 和 b 的顺序）：
% Ax = b
x = cv.solve(A, b); % A\b or mldivide(A,b)
% xA = b
x = cv.solve(A', b')'; % b/A or mrdivide(b,A)
OpenCV 暴露的 API 在这里有点尴尬，所以我们不得不做所有这些转置。事实上，如果你参考 equivalent LAPACK routines (think DGESV or DGESVX) they actually allow you to specify whether the matrix is transposed TRANS=T or not TRANS=N (at that level, transposition is really just a different memory layout, C or Fortran ordering). MATLAB for instance provides the linsolve 函数，它允许你在选项中指定这些东西......
（顺便说一句，在 C++ OpenCV 中编码时，我更喜欢使用 function-form 操作，例如 cv::transpose as opposed to the matrix-expression variants like Mat::t。前者可以操作 in-place，而后者会创建不必要的临时副本） .
现在，如果您正在寻找性能良好的 linear-algebra C++ 实现，请考虑使用 Eigen (it even integrate nicely with OpenCV)。此外，它是一个纯 template-based 库，因此无需担心链接或二进制文件，只需包含 header 文件。
编辑（回应评论）
@Goz:
look up Return Value Optimisation.
The "unnecessary temporary copies" don't exist
我知道 RVO and move semantics, but it's not of much importance here; The cv::Mat class is copy-friendly 无论如何，有点像 reference-counted 智能指针。这意味着它只在传递 by-value 时进行数据共享的浅拷贝。为新副本创建的唯一部分是垫子 header 中的部分，这些部分在大小方面无关紧要（存储 [=209= 的数量、步长和数据类型等内容）。
我说的是一个明确的深拷贝，而不是你在从函数调用返回时正在考虑的那个...
感谢您的评论，这让我有动力去真正挖掘 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;
}
这个 returns 一个 MatExpr 是一个 class 允许延迟计算 matrix expressions。换句话说，它不会立即执行转置，而是存储对原始矩阵的引用和最终对其执行的操作（转置），但它会推迟评估它直到绝对必要（例如当分配或转换为 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;
}
到目前为止这很简单。 class 应该将输入矩阵存储在 a 中（同样 cv::Mat 实例将共享数据，因此无需复制），以及要执行的操作op，还有其他一些对我们来说不重要的东西。
下面是矩阵运算class MatOp，还有一些子classes（我只展示了转置和逆运算，但还有更多）：
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等各种运算实际上映射到相应的矩阵表达式运算。
谜题的最后一部分是代理 class 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 如何创建和 returns 一个 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 被转置 in-place，并且抽象级别较低。
现在我展示所有这些的原因不是为了争论这个额外的副本，毕竟这没什么大不了的:)
我发现的真正巧妙的事情是以下两个表达式没有做同样的事情并给出不同的结果（而且我不是在谈论逆是否是 in-place 或不是）：
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 函数的包装器）。 IMO 这很聪明！即使您写了 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.
这两种情况
希望这对您有所帮助，抱歉这么长 post ;)
我有 2 个 19x19 方阵 (a & b)，我正在尝试使用斜杠 (mrdivide) 运算符执行除法，使得
c = a / b
我正在尝试在 OpenCV 中实现它。我发现有几个人建议使用 cv::solve，但到目前为止，我一直无法找到任何能给我提供接近 matlab 的结果的东西。
有人知道我如何用 opencv 实现 mrdivide 吗？
我试过以下代码：
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();
}
这给出了与上面相同的答案，但也是错误的。
在 MATLAB 中，在两个兼容维度的矩阵上使用 mrdivide，使得 a / b 等同于 a * b^{-1}，其中 b^{-1} 是 b 的倒数.因此，您可以做的可能是先反转矩阵 b，然后预乘此 a。
一种方法是在矩阵b上使用cv::invert，然后预乘a。这可以通过以下函数定义来完成（借用上面 post 中的代码）：
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() 中指定方法即可。
这个 Whosebug post 使用三种方法对相对较大的矩阵求逆矩阵进行计时：Fastest method in inverse of matrix
您不应该使用 inv 来求解 Ax=b 或 xA=b 方程。虽然这两种方法在数学上是等效的（x=solve(A,b) 和 x=inv(A)*B），但在处理 floating-point 数字时完全不同！
http://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/
作为一般规则，永远不要乘以逆矩阵。而是对 one-off 系统使用 forward/backward-slash 运算符（或等效的 "solve" 方法），或者在需要时显式执行矩阵分解（想想 LU、QR、Cholesky 等）对多个 b 的
A
让我举一个具体的例子来说明反转的问题。我将使用 MATLAB 和 mexopencv，一个允许我们直接从 MATLAB 调用 OpenCV 的库。
（这个例子是蒂姆戴维斯从 this excellent FEX submission 那里借来的，他是 SuiteSparse 背后的同一个人。我展示的是 left-division Ax=b 的情况，但是同样适用于 right-division 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
现在让我们比较cv::invert against cv::solve by finding the solution and computing the residual error using the NORM function (or cv::norm，如果你愿意的话）：
% 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 = A\b;
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 的想法是正确的，但是在 right-division.
在 MATLAB 中，运算符 / 和 \（或 mrdivide and mldivide) are related to each other by the equation B/A = (A'\B')' (this is a simple result of transpose properties）。
因此对于 OpenCV 函数，您将编写（注意 A 和 b 的顺序）：
% Ax = b
x = cv.solve(A, b); % A\b or mldivide(A,b)
% xA = b
x = cv.solve(A', b')'; % b/A or mrdivide(b,A)
OpenCV 暴露的 API 在这里有点尴尬，所以我们不得不做所有这些转置。事实上，如果你参考 equivalent LAPACK routines (think DGESV or DGESVX) they actually allow you to specify whether the matrix is transposed TRANS=T or not TRANS=N (at that level, transposition is really just a different memory layout, C or Fortran ordering). MATLAB for instance provides the linsolve 函数，它允许你在选项中指定这些东西......
（顺便说一句，在 C++ OpenCV 中编码时，我更喜欢使用 function-form 操作，例如 cv::transpose as opposed to the matrix-expression variants like Mat::t。前者可以操作 in-place，而后者会创建不必要的临时副本） .
现在，如果您正在寻找性能良好的 linear-algebra C++ 实现，请考虑使用 Eigen (it even integrate nicely with OpenCV)。此外，它是一个纯 template-based 库，因此无需担心链接或二进制文件，只需包含 header 文件。
编辑（回应评论）
@Goz:
look up Return Value Optimisation. The "unnecessary temporary copies" don't exist
我知道 RVO and move semantics, but it's not of much importance here; The cv::Mat class is copy-friendly 无论如何，有点像 reference-counted 智能指针。这意味着它只在传递 by-value 时进行数据共享的浅拷贝。为新副本创建的唯一部分是垫子 header 中的部分，这些部分在大小方面无关紧要（存储 [=209= 的数量、步长和数据类型等内容）。
我说的是一个明确的深拷贝，而不是你在从函数调用返回时正在考虑的那个...
感谢您的评论，这让我有动力去真正挖掘 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;
}
这个 returns 一个 MatExpr 是一个 class 允许延迟计算 matrix expressions。换句话说，它不会立即执行转置，而是存储对原始矩阵的引用和最终对其执行的操作（转置），但它会推迟评估它直到绝对必要（例如当分配或转换为 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;
}
到目前为止这很简单。 class 应该将输入矩阵存储在 a 中（同样 cv::Mat 实例将共享数据，因此无需复制），以及要执行的操作op，还有其他一些对我们来说不重要的东西。
下面是矩阵运算class MatOp，还有一些子classes（我只展示了转置和逆运算，但还有更多）：
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等各种运算实际上映射到相应的矩阵表达式运算。
谜题的最后一部分是代理 class 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 如何创建和 returns 一个 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 被转置 in-place，并且抽象级别较低。
现在我展示所有这些的原因不是为了争论这个额外的副本，毕竟这没什么大不了的:) 我发现的真正巧妙的事情是以下两个表达式没有做同样的事情并给出不同的结果（而且我不是在谈论逆是否是 in-place 或不是）：
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 函数的包装器）。 IMO 这很聪明！即使您写了 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.
希望这对您有所帮助，抱歉这么长 post ;)