用2平方矩阵模拟matlab的mrdivide

Question

我有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();
}

这给出了与上面相同的答案,但也是错误的.

Answer 1

You should NOT be using inv to solve Ax=b or xA=b equations. While the two methods are mathematically equivalent (x=solve(A,b) and x=inv(A)*B), it's a completely different thing when working with floating-point numbers! http://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/

As a general rule, never multiply by a matrix inverse. Instead use the forward/backward-slash operators (or the equivalent "solve" methods) for a one-off system, or explicitly perform matrix factorization (think LU, QR, Cholesky, etc..) when you want to reuse the same A with multiple b's

让我举一个具体的例子来说明反转的问题.我将使用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

现在,让我们比较cv::invert反对cv::solve通过寻找解决方案,并使用NORM函数计算残差(或者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,但你只是在右分割的情况下得到操作数的错误顺序.

在MATLAB中,运营商/和\(或mrdivide和mldivide)由方程彼此相关的B/A = (A'\B')'(这是一个简单的结果转置属性).

因此,与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在这里有点尴尬,所以我们不得不做所有这些转置.事实上,如果你指的是等效的 LAPACK例程(认为DGESV或DGESVX)他们实际上允许你指定的矩阵是否被调换TRANS=T或不TRANS=N(在这个水平,换位真的只是一个不同的内存布局,C或Fortran语言排序).例如MATLAB提供的linsolve功能允许您在选项中指定这些类型的东西......

(顺便说一下,当在C++ OpenCV中编码时,我更喜欢使用函数形式的操作cv::transpose,而不是像矩阵表达变体那样Mat::t.前者可以就地操作而后者可以创建不必要的临时副本).

现在,如果您正在寻找C++中的良好性能线性代数实现,请考虑使用Eigen(它甚至可以很好地与OpenCV集成).此外,它是一个纯模板库,所以不需要担心链接或二进制文件,只需包含头文件.

编辑(回应评论)

@Goz:

查找返回值优化."不必要的临时副本"不存在

我知道RVO和移动语义,但它在这里并不重要; 无论如何,这个cv::Mat类都是复制友好的,有点像引用计数的智能指针.这意味着它只在传递by-value时执行带有数据共享的浅拷贝.为新副本创建的唯一部分是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();
    ...
}

Now we have to figure out the steps in between. First thing we have is the t member method:

MatExpr Mat::t() const
{
    MatExpr e;
    MatOp_T::makeExpr(e, *this);
    return e;
}

This returns a MatExpr which is a class that allows lazy evaluation of matrix expressions. In other words, it will not perform the transpose right away, instead it stores a reference to the original matrix and the operation to eventually perform on it (transpose), but it will hold off on evaluating it until it is absolutely necessary (for example when assigned or cast to a cv::Mat).

Next let's see the definitions of the relevant parts. Note that in the actual code, these things are split across many files. I only pieced the interesting parts here for easier reading, but it's far from the complete thing:

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;
}

So far this is straightforward. The class is supposed to store the input matrix in a (again cv::Mat instances will share data, so no copying), along with the operation to perform op, and a few other things not important for us.

Here's the matrix operation class MatOp, and some of it subclasses (I'm only showing the transpose and inverse operations, but there are more):

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 is heavily using operator overloading, so all kinds of operations like A+B, A-B, A*B, ... actually map to corresponding matrix expression operations.

The final part of the puzzle is the proxy class InputArray. It basically stores a void* pointer along with info about the thing passed (what kind it is: Mat, MatExpr, Matx, vector<T>, UMat, etc..), that way it knows how to cast the pointer back when requested with something like 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;
}

So now we see how Mat::t creates and returns a MatExpr instance. Which is then received by cv::solve as InputArray. Now when it calls InputArray::getMat to retrieve the matrix, it effectively converts the stored MatExpr to a Mat calling the cast operator:

    MatExpr::operator Mat() const
    {
        Mat m;
        op->assign(*this, m);
        return m;
    }

so it declares a new matrix m, calls MatOp_T::assign with the new matrix as destination. In turn this forces it to evaluate by finally calling cv::transpose. It computes the transposed result into this new matrix as destination.

So we end up having two copies, the original A and the transposed A.t() returned.

Now with all that said, compare it against:

Mat A = ..., b = ..., x;
cv::transpose(A, A);
cv::solve(A, b, x);

In this case, A is transposed in-place, and with less levels of abstraction.

Now the reason I showed all of that is not to argue about this one extra copy, after all it's not that big of a deal :) The really neat thing I found out is that the following two expressions are not doing the same thing and give different results (and I'm not talking about whether the inverse is in-place or not):

Mat A = ..., b = ..., x;
cv::invert(A,A);
x = A*b;

Mat A = ..., b = ..., x;
x = inv(A)*b;

It turns out that the second one is in fact smart enough to call cv::solve(A,b)! If you go back to MatOp_Invert::matmul (which is called when a lazy invert is later chained with another lazy matrix multiplication).

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);
    ...
}

it checks if the first operand in the expression inv(A)*B is an invert operation and the second operand is an identity operation (i.e a plain matrix, not the result of another complicated expression). In that case it changes the stored operation to a lazy solve operation MatOp_Solve (which similarly is a wrapper for the cv::solve function). IMO that's pretty smart! Even though you wrote inv(A)*b, it wont actually compute the inverse, and instead it understands that it is better to rewrite it by solving the system using matrix factorization.

Unfortunately for you, this will only benefit expressions of the form inv(A)*b not the other way around b*inv(A) (that one will end up computing the inverse which is not what we want). So in your case of solving xA=b, you should stick with explicitly calling cv::solve...

Of course this is only applicable when coding in C++ (thanks to the magic of operator overloading and lazy expressions). If you're using OpenCV from another language using some wrappers (like Python, Java, MATLAB), you're probably not getting any of that, and should be explicit in using cv::solve like I did in the previous MATLAB code, for both cases Ax=b and xA=b.

Hope this helps, and sorry for the long post ;)