opencv CvSolve函数深度解析

Opencv CvSolve函数主要是用来求解线性系统Ax=b的方程，X的解。solve函数跟它的算法是一样的，也是用来求解线性系统。

        设方程Ax = b.根据有效的方程个数和未知数的个数，可以分为以下3种情况：

1）rank（A） < n，也就是说方程个数小于未知数的个数，约束不够，方程存在无数组解，

2)  rank（A） =  n  方程个数等于未知数的个数， 方程存在唯一的精确解，解法通常有我们熟悉的消元法，LU分解法

3)  rank(A) > n，方程个数多于未知数个数，这个时候约束过于严格，没有精确解，这种方程又称之为超定方程。通常工程应用都会遇到这种情况，找不到精确解的情况下，我们选取最优解。这个最优解，又称之为最小二乘解。

前面2种情况是比较好理解的，我们在这里就不多说了，我们重点研究的是第3种情况，也是我们应用中碰到最多最常见的情况。

这个超定方程的最优解是怎么求的呢？常规的消元法是无解的，于是有人提出并证明了这样一个理论：

以下是定理3的证明：

我们把Ax=b的问题直接转化成A'A x = A'b的求解，那么求超定线性方程组的方法，应该就转换成（A'*A)x = (A'*b)咯。而这个新的方程又具有什么样的特性呢？

那么，这个应该是常规线性方程求解啦，为什么不能直接LU分解或是用消元法呢，这个文档里说了，也确实是可以的
http://wenku.baidu.com/link?url=bTW3JfQzwzXtMWA2JvF4fCS_sXmWYelI2Yx7jKyquo-P8vyKwBb2qxYYJyIgnMFV_rc-PF6yQQaDEmNOR7jGSowVpShTMDg9W_1S_t3Uv7G

update(2017.5.8):这是有前提的，必须是R(A) = n时，AtAx = Atb才有唯一解，否则，我们求的还是最小二乘解。

而在solve函数的代码中，用特征值分解（Jacobi）呢？

Jacobi方法用于求解对称矩阵的特征值，对称矩阵不一定是满秩矩阵，不一定能够顺利的LU分解，但是可以作SVD分解！如果是调用cvSolve的话，求最小二乘解，

那么最后一项要传入CV_SVD.不管是SVD还是jacobi都是跟特征值相关。Jacobi求奇异值分解的方法，

详情参考http://blog.csdn.net/k531623594/article/details/50628163.

代码片段加注解：

template<typename _Tp> bool
JacobiImpl_( _Tp* A, size_t astep, _Tp* W, _Tp* V, size_t vstep, int n, uchar* buf )
{
    const _Tp eps = std::numeric_limits<_Tp>::epsilon();
    int i, j, k, m;

    astep /= sizeof(A[0]);
    if( V )
    {
        vstep /= sizeof(V[0]);
        for( i = 0; i < n; i++ )
        {
            for( j = 0; j < n; j++ )
                V[i*vstep + j] = (_Tp)0;
            V[i*vstep + i] = (_Tp)1;
        }
    }

    int iters, maxIters = n*n*30;

    int* indR = (int*)alignPtr(buf, sizeof(int));
    int* indC = indR + n;
    _Tp mv = (_Tp)0;

    for( k = 0; k < n; k++ )
    {
        W[k] = A[(astep + 1)*k];
        if( k < n - 1 )
        {
            for( m = k+1, mv = std::abs(A[astep*k + m]), i = k+2; i < n; i++ )
            {
                _Tp val = std::abs(A[astep*k+i]);
                if( mv < val )
                    mv = val, m = i;
            }
            indR[k] = m; //第a[k][k]右边最大的非对角元素的列下标
        }
        if( k > 0 )
        {
            for( m = 0, mv = std::abs(A[k]), i = 1; i < k; i++ )
            {
                _Tp val = std::abs(A[astep*i+k]);
                if( mv < val )
                    mv = val, m = i;
            }
            indC[k] = m;//第a[k][k]上边最大的非对角元素的行下标
        }
    }

    if( n > 1 ) for( iters = 0; iters < maxIters; iters++ )
    {
        // find index (k,l) of pivot p
        for( k = 0, mv = std::abs(A[indR[0]]), i = 1; i < n-1; i++ )
        {
            _Tp val = std::abs(A[astep*i + indR[i]]);
            if( mv < val )
                mv = val, k = i;
        }
        int l = indR[k];
        for( i = 1; i < n; i++ )
        {
            _Tp val = std::abs(A[astep*indC[i] + i]);
            if( mv < val )
                mv = val, k = indC[i], l = i;
        }

        _Tp p = A[astep*k + l];
        if( std::abs(p) <= eps )
            break;
        _Tp y = (_Tp)((W[l] - W[k])*0.5);
        _Tp t = std::abs(y) + hypot(p, y);
        _Tp s = hypot(p, t);
        _Tp c = t/s;
        s = p/s; t = (p/t)*p;
        if( y < 0 )
            s = -s, t = -t;
        A[astep*k + l] = 0;

        W[k] -= t;
        W[l] += t;

        _Tp a0, b0;

#undef rotate
#define rotate(v0, v1) a0 = v0, b0 = v1, v0 = a0*c - b0*s, v1 = a0*s + b0*c

        // rotate rows and columns k and l
        for( i = 0; i < k; i++ )
            rotate(A[astep*i+k], A[astep*i+l]);
        for( i = k+1; i < l; i++ )
            rotate(A[astep*k+i], A[astep*i+l]);
        for( i = l+1; i < n; i++ )
            rotate(A[astep*k+i], A[astep*l+i]);

        // rotate eigenvectors
        if( V )
            for( i = 0; i < n; i++ )
                rotate(V[vstep*k+i], V[vstep*l+i]);

#undef rotate

        for( j = 0; j < 2; j++ )
        {
            int idx = j == 0 ? k : l;
            if( idx < n - 1 )
            {
                for( m = idx+1, mv = std::abs(A[astep*idx + m]), i = idx+2; i < n; i++ )
                {
                    _Tp val = std::abs(A[astep*idx+i]);
                    if( mv < val )
                        mv = val, m = i;
                }
                indR[idx] = m;
            }
            if( idx > 0 )
            {
                for( m = 0, mv = std::abs(A[idx]), i = 1; i < idx; i++ )
                {
                    _Tp val = std::abs(A[astep*i+idx]);
                    if( mv < val )
                        mv = val, m = i;
                }
                indC[idx] = m;
            }
        }
    }
	//降序排列特征值和相应的特征向量
    // sort eigenvalues & eigenvectors
    for( k = 0; k < n-1; k++ )
    {
        m = k;
        for( i = k+1; i < n; i++ )
        {
            if( W[m] < W[i] )
                m = i;
        }
        if( k != m )
        {
            std::swap(W[m], W[k]);
            if( V )
                for( i = 0; i < n; i++ )
                    std::swap(V[vstep*m + i], V[vstep*k + i]);
        }
    }

    return true;
}

      值得一提的是，当b = 0时，超定齐次线性方程 A*X =0;求解最小二乘解，在||X||=1的约束下，其最小二乘解为矩阵A'A最小特征值所对应的特征向量。

参考http://blog.csdn.net/ccwwff/article/details/45912671

简单证明：

min ||Ax|| 
s.t    ||x||=1
目标函数：||Ax|| = x'A'Ax = x'lamda x=lamda||x||=lamda,其中lamda是A'A的特征值。
于是可知，得到了A'A的最小特征值，就得到了最优值，而其最小特征值对应的特征向量就是最优解.