c++ - 用OpenCV模拟matlab的mldivide

Question

我昨天问了这个问题:Simulating matlab's mrdivide with 2 square matrices

这就是 mrdivide 工作的原因。但是现在我遇到了 mldivide 的问题，目前实现如下:

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 正在运行的事实应该意味着 mldivide 正在运行，但我无法得到与 matlab 相同的结果。同样，结果完全不同。

值得注意的是我正在尝试做一个 [19x19]\[19x200] 所以这次不是方阵。

Answer 1

就像我之前在您的 other question 中提到的那样, 我正在使用 MATLAB 和 mexopencv ，这样我就可以轻松比较 MATLAB 和 OpenCV 的输出。

也就是说，我无法重现您的问题:我随机生成矩阵，并重复比较 N=100 次。我正在使用针对 OpenCV 3.0.0 编译的 mexopencv 运行 MATLAB R2015a:

N = 100;
r = zeros(N,2);
d = zeros(N,1);
for i=1:N
    % double precision, i.e CV_64F
    A = randn(19,19);
    B = randn(19,200);

    x1 = A\B;
    x2 = cv.solve(A,B);   % this a MEX function that calls cv::solve

    r(i,:) = [norm(A*x1-B), norm(A*x2-B)];
    d(i) = norm(x1-x2);
end

所有结果一致，误差非常小，大约为 1e-11:

>> mean(r)
ans =
   1.0e-12 *
    0.2282    0.2698

>> mean(d)
ans =
   6.5457e-12

(顺便说一句，我还尝试了 x2 = cv.solve(A,B, 'IsNormal',true); 设置了 cv::DECOMP_NORMAL 标志，并且结果也没有那么不同)。

这让我相信要么您的矩阵恰好突出了 OpenCV 求解器中的某些边缘情况，在这种情况下它无法给出正确的解决方案，要么更有可能是您的代码中的其他地方存在错误。

我首先要仔细检查您如何加载数据，尤其要注意矩阵的布局方式(显然 MATLAB 是列优先的，而 OpenCV 是行优先的)...

此外，您从未告诉我们有关您的矩阵的任何信息；它们是否表现出某种特征，是否存在任何对称性，它们是否主要为零，它们的秩等。

在 OpenCV 中，默认求解器方法是 LU 分解，如果合适，您必须自己明确更改它。 MATLAB手上会automatically选择最适合矩阵 A 的方法，LU 只是可能的分解之一。

---

编辑(回复评论)

在MATLAB中使用SVD分解时，左右特征向量U和V的符号为arbitrary (这真的来自 DGESVD LAPACK 例程)。为了得到一致的结果，一个约定是要求每个特征向量的第一个元素是某个符号，然后将每个 vector 乘以+1或-1以适本地翻转符号。我还建议查看 eigenshuffle .

再一次，这是我所做的一项测试，以确认我在 MATLAB 和 OpenCV 中获得了类似的 SVD 结果:

N = 100;
r = zeros(N,2);
d = zeros(N,3);
for i=1:N
    % double precision, i.e CV_64F
    A = rand(19);

    % compute SVD in MATLAB, and apply sign convention
    [U1,S1,V1] = svd(A);
    sn = sign(U1(1,:));
    U1 = bsxfun(@times, sn, U1);
    V1 = bsxfun(@times, sn, V1);
    r(i,1) = norm(U1*S1*V1' - A);

    % compute SVD in OpenCV, and apply sign convention
    [S2,U2,V2] = cv.SVD.Compute(A);
    S2 = diag(S2);
    sn = sign(U2(1,:));
    U2 = bsxfun(@times, sn, U2);
    V2 = bsxfun(@times, sn', V2)';  % Note: V2 was transposed w.r.t V1
    r(i,2) = norm(U2*S2*V2' - A);

    % compare
    d(i,:) = [norm(V1-V2), norm(U1-U2), norm(S1-S2)];
end

同样，所有结果都非常相似，误差接近机器 epsilon 并且可以忽略不计:

>> mean(r)
ans =
   1.0e-13 *
    0.3381    0.1215

>> mean(d)
ans =
   1.0e-13 *
    0.3113    0.3009    0.0578

有一件事我在 OpenCV 中不确定，但是 MATLAB 的 svd 函数返回按降序排序的奇异值(与 eig 函数不同)，其中列相应顺序的特征向量。

现在，如果由于某种原因不能保证对 OpenCV 中的奇异值进行排序，如果您想将结果与 MATLAB 进行比较，则也必须手动进行排序，如:

% not needed in MATLAB
[U,S,V] = svd(A);
[S, ord] = sort(diag(S), 'descend');
S = diag(S);
U = U(:,ord)
V = V(:,ord);