OpenCV源码解析：Jacobi法计算矩阵的特征值和特征向量

（注：CSDN不适合写公式，只好上传图片格式）

 

其中Pkk=Pll=cosθ, Plk=Pkl=sinθ,形式上就是这样，

A*PT

   Aik = Aik×Pkk+Ail×Pkl

   Ail = Aik×Plk+Ail×Pll

P*A

   Aki = Pkk×Aki+ Pkl×Ali

   Ali = Plk×Aki+ Pll×Ali

实际计算时，只计算那些必须计算的量，已经归0的量和不变的量都不用计算，比如,当k=1,l=3时，只计算了下面这6个量A01,A03, A12, A23, A14, A34,其中A13已经归0化(A03=0),

为了达到上述目的，OpenCV充分利用了矩阵的对称性（如Ali=Ail）， 连续采用了三个for循环计算右上角的元素(对右上角的元素， 一定成立），这样既可以避免左下角不必要的运算，同时也免除了向左下角拷贝元素值的必要。

第1个for循环，处理0到k-1行的第k,l个元素，这些元素在对角线的上面；

第2个for循环，k+1到l-1之间的元素，为了确保只用右上角的元素，后面写成了A[astep*i+l]；

第3个for循环，l+1到n-1之间的元素，在对角线的右边；

        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]);

那些对角线的上元素（最后得到的就是特征值），则由下面的代码处理，

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

看到这里，不难明白，为什么OpenCV的eigen/Jacobi函数只适合实对称矩阵的处理，而不适合那些非对称的矩阵。

另外要说明一下源码中的最大值查找，它是冗余的，函数先实现了按行查找最大值，然后又实现了按列查找最大值，对于非对称矩阵，这种算法确保能找到最大值所在的位置；对于对称矩阵，这两种算法就完全重复了，两次查找得到的结果是一样的。

下面的我们看一下OpenCV中Jacob算法的具体实现，源码在lapack.cpp中

Hypot计算的是(a^2 + b^2)^0.5
template<typename _Tp> static inline _Tp hypot(_Tp a, _Tp b)
{
    a = std::abs(a);
    b = std::abs(b);
    if( a > b )
    {
        b /= a;
        return a*std::sqrt(1 + b*b);
    }
    if( b > 0 )
    {
        a /= b;
        return b*std::sqrt(1 + a*a);
    }
    return 0;
}

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; 	// 非对角线元素置0   
            V[i*vstep + i] = (_Tp)1;  		// 所有对角线元素置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];  // 取第k行对角线上的元素（第k行:astep*k,第k列:k）

        // 下面这个if，是对k行右边所有的元素逐个查找，找到绝对值最大的那个，并记录下其所在的位置
        if( k < n - 1 )
        {
		   //第k行中，从第k+1到第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;  // k行中最大值元素所在的位置
        }
		
	    // 下面这个if，是对k列上边所有的元素逐个查找，找到绝对值最大的那个，并记录下其所在的位置
        if( k > 0 ) 
        {
		   //第k列行中，从第0到第k-1个元素，逐个找
            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;  // 第k列中最大值元素所在的位置（行数）
        }
    }

    if( n > 1 ) for( iters = 0; iters < maxIters; iters++ )
    {
        // find index (k,l) of pivot p
        // 从第0行到第n-2行，每行最大的那个值相互比对，逐行找最大值， 
        // 注意，最后一行是第n-1行，右边已经没有元素了，故而不参与计算
        for( k = 0, mv = std::abs(A[indR[0]]), i = 1; i < n-1; i++ )
        {
            _Tp val = std::abs(A[astep*i + indR[i]]); // 第i行的最大值
            if( mv < val )
                mv = val, k = i;
        }
        int l = indR[k];  // 右上角元素最大值所在的列		
	    
        // 从第1列到第n-1列，逐列找 
        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

        // 因为矩阵已经变更，所以第k,l行和第k,l列的最大，最小值可能已经变了，需要更新
        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;
}

static bool Jacobi( float* S, size_t sstep, float* e, float* E, size_t estep, int n, uchar* buf )
{
    return JacobiImpl_(S, sstep, e, E, estep, n, buf);
}

static bool Jacobi( double* S, size_t sstep, double* e, double* E, size_t estep, int n, uchar* buf )
{
    return JacobiImpl_(S, sstep, e, E, estep, n, buf);
}

解析完毕。