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);
}
解析完毕。
