由于最近利用PCA统计图像主方向，用到了特征值求解算法，网上搜罗看普通解法是jacobi算法，也正好复习一下，在此做一分享，希望对有需要的人有所帮助：
###原理
jacobi算法用于求实对称矩阵的全部特征值及特征向量，利用平面旋转矩阵对矩阵A做相似变换，化A为对角阵，进而求出特征值及特征向量
对于p!=q,下面定义n阶矩阵Uqp是平面旋转矩阵

变换过程

在保证相似条件下，使主对角线外元素趋于0，jacobi算法是一种迭代过程，每一步迭代，需要求出旋转后的新的矩阵，其过程如下：


###算法过程

  • 初始化特征向量Vec为对角阵
  • 在A的非主对角线元素中，找到一个绝对值最大的元素A[p,q]
  • 利用公式1.3计算theta值
  • 利用公式1、2、3计算特征向量及特征值
  • 当迭代前非对角元素最大值小于阈值或迭代次数小于给定次数、停止计算
  • 对特征值及特征向量进行排序

###code

void Jacobi(double *A, int dims, double *eigenvalue, double *vec, int times, double errorValue)
{
    // initial vector
    memset(vec, 0, sizeof(double) * dims * dims);
    for (int i = 0; i < dims; i++)
    {
        vec[i * dims + i] = 1;
    }
    int nCount = 0;
    while (1)
    {
        // 在非对角元素上找最大元素
        double maxValue = abs(A[1]);
        int row = 0;
        int col = 1;
        for (int i = 0; i < dims; i++)
        {
            for (int j = 0; j < dims; j++)
            {
                if (i != j && abs(A[i * dims + j]) > maxValue)
                {
                    maxValue = abs(A[i * dims + j]);
                    row = i;
                    col = j;
                }
            }
        }
        if (maxValue < errorValue || nCount > times)
        {
            break;
        }
        double Apq = A[row * dims + col];
        double App = A[row * dims + row];
        double Aqq = A[col * dims + col];
        double theta = 0.5 * atan2(-2 * Apq, Aqq - App);
        double cost = cos(theta);
        double sint = sin(theta);
        A[row * dims + row] = App * cost * cost + Aqq * sint * sint + 2 * Apq * cost * sint;
        A[col * dims + col] = App * sint * sint + Aqq * cost * cost - 2 * Apq * cost * sint;
        A[row * dims + col] = 0.5 * (Aqq - App) * sin(2 * theta) + Apq * cos(2 * theta);
        A[col * dims + row] = A[row * dims + col];

        for (int i = 0; i < dims; i++)
        {
            if ( i != row && i != col)
            {
                maxValue = A[i * dims + row];
                A[i * dims + row] = cost * maxValue + sint * A[i * dims + col];
                A[i * dims + col] = cost * A[i * dims + col] - sint * maxValue;    
            }
        }
        for (int j = 0; j < dims; j++)
        {
            if (j != row && j != col)
            {
                maxValue = A[row * dims + j];
                A[row * dims + j] = cost * maxValue + sint * A[col * dims + j];
                A[col * dims + j] = cost * A[col * dims + j] - sint * maxValue;
            }
        }
        // 计算特征向量
        for (int i = 0; i < dims; i++)
        {
            maxValue = vec[i * dims + row];
            vec[i * dims + row] = cost * maxValue + sint * vec[i * dims + col];
            vec[i * dims + col] = cost * vec[i * dims + col] - sint * maxValue;
        }
        nCount++;
    }
    // 对特征值进行排序，特征值既A主对角线上的元素
    // map 在插入元素时，会自动按照从大道小排序
    map<double, int> mapEigen;
    for (int i = 0; i < dims; i++)
    {
        eigenvalue[i] = A[i * dims + i];
        mapEigen.insert(make_pair(eigenvalue[i], i));
    }
    double *vecTemp = new double[dims * dims];
    map<double, int>::reverse_iterator iter = mapEigen.rbegin();
    for (int j = 0; iter != mapEigen.rend(), j < dims; iter++, j++)
    {
        for (int i = 0; i < dims; i++)
        {
            vecTemp[i * dims + j] = vec[i * dims + iter->second];
        }
        // 特征值重排序
        eigenvalue[j] = iter->first; // key value 
    }
    memcpy(vec, vecTemp, sizeof(double) * dims * dims);
    delete[] vecTemp;
}

文中用到了stl中的标准c++库map，其在添加元素时，自动对特征值进行排序，用时还是特别方便的。

###Mark
对于一个实例如下图所示

其中a为输入矩阵，v和c分别为其特征向量及特征值，黑色背景为上面代码所演算特征值及特征向量，可以看出，有一个特征向量符号正好相反，可能原因是对称矩阵在对角化过程中结果不唯一性所致，其结果仍然是正确的。