朱金灿
因为要对一个矩阵求特征向量矩阵和特征值，本想使用MTL库，因为MTL从成熟程度和运算效率都有保证，没想到使用MTL库求特征向量矩阵和特征值还要依赖其它库，只好另辟蹊径了。幸好我找到了OpenCV，一个由intel资助的开源库，全称叫Open Source Computer Vision Library(大意就是开源计算机视觉库吧)。
要下载OpenCV库，可以从OpenCV中文站上下载，网址是：http://www.opencv.org.cn/index.php/Download。下载下来，我发现OpenCV库做得真的是很人性化——可以通过安装包进行安装。安装好后已经有支持VS 2005的解决方案文件了，你基本不用配置编译选项了。因为我只想里面的矩阵运算，因此这次我只编译了cxcore工程（呵呵，更深的学习还有待以后）。运行VS 2005的生成命令后生成库文件，debug模式下生成：cxcored.lib和cxcore100d.dll，release模式下生成cxcore100.dll和cxcore.lib（可以看到debug文件比release文件多了一个字母d）。
现在正式开始使用OpenCV库，要使用它，当然得把相关头文件拷出来，这里使用它的矩阵运算就是把安装目录下的cv文件夹拷出来，然后在你的程序中include。
使用库的第一个问题就是实现标准类型向库类型的转化（有时我就搞不懂，为何库不使用标准C的类型呢？比如char、int之类的，也许这个是出于效率方面的考虑吧）。OpenCV库求特征值和特征向量的函数是：
void cvEigenVV( CvArr* mat, CvArr* evects, CvArr* evals, double eps=0 );
mat
输入对称方阵。在处理过程中将被改变。
evects
特征向量输出矩阵， 连续按行存储
evals
特征值输出矩阵，按降序存储(当然特征值和特征向量的排序是同步的)。
eps
对角化的精确度 (典型地， DBL_EPSILON=≈10-15 就足够了)。
测试开始，建一个控制台程序，写一个测试函数：
void TestOpenCV(void)
{
// a为输入对称方阵
double a[5][5] =
{
{10.0, 1.0, 2.0, 3.0, 4.0},
{ 1.0, 9.0, -1.0, 2.0, -3.0},
{ 2.0, -1.0, 7.0, 3.0, -5.0},
{ 3.0, 2.0, 3.0, 12.0, -1.0},
{ 4.0, -3.0, -5.0, -1.0, 15.0}
};
// 构造输入方阵
CvMat SrcMatrix =cvMat(5,5, CV_64FC1,a);
double b[5][5] =
{
{0.0, 0.0, 0.0, 0.0, 0.0},
{ 0.0, 0.0, 0.0, 0.0, 0.0},
{ 0.0, 0.0, 0.0, 0.0, 0.0},
{ 0.0, 0.0, 0.0, 0.0, 0.0},
{ 0.0, 0.0, 0.0, 0.0, 0.0}
};
// 构造输出特征向量矩阵
CvMat ProVector =cvMat(5, 5, CV_64FC1, b);
// 构造输出特征值矩阵
double c[5] = {0.0, 0.0, 0.0, 0.0, 0.0};
CvMat ProValue = cvMat(5,1, CV_64FC1, c);
cvEigenVV(&SrcMatrix,&ProVector,&ProValue,1.0e-6F);
printf("/n");
int i =0;
int j =0;
for(i=0;i<5;i++)
{
for(j=0;j<5;j++)
printf("%f/t",cvmGet(&ProVector,i,j));
printf("/n");
}
for(i=0;i<5;i++)
{
printf("%f/n",cvmGet(&ProValue,i,0));
}
}
使用cvEigenVV函数要注意以下三点：
1. 构造对称方阵、特征向量输出矩阵和特征值输出矩阵必须分静态数组和动态开辟内存的数组，如果是静态数组如上面测试程序就可以使用上面的构造矩阵的方式，如果是动态开辟内存的数组，那么在构造矩阵类时就要使用下面的方式：
double** pfMatrix =new double*[nInBandNum];
double** pfVector=new double*[nInBandNum];
double* pProValue =new double[nInBandNum];
for(i=0;i<nInBandNum;i++)
{
pfMatrix[i]=new double[nInBandNum];
pfVector[i]=new double[nInBandNum];
pProValue[i]= 0.0;
}
// 初始化参数
for(i=0;i<nInBandNum;i++)
{
for(j=0;j<nInBandNum;j++)
{
pfMatrix[i][j]=0.0;
pfVector[i][j] = 0.0;
}
}
// 对输入对称方阵进行赋值，这一步放过
// 下面是构造输入方阵、特征向量矩阵和特征值矩阵
CvMat *pSrcMatrix = cvCreateMat(nBandNum, nBandNum, CV_64FC1);
CvMat *pProVector = cvCreateMat(nBandNum, nBandNum, CV_64FC1);
CvMat *pProValue = cvCreateMat(nBandNum,1, CV_64FC1);
for(i=0;i<nBandNum;i++)
{
cvmSet(pProValue,i,0,pfProVector[i]);
for(j=0;j<nBandNum;j++)
{
cvmSet(pSrcMatrix,i,j,pfMatrix[i][j]);
cvmSet(pProVector,i,j,pfVector[i][j]);
}
}
cvEigenVV(pSrcMatrix,pProVector,pProValue);
// 释放矩阵所占内存空间
cvReleaseMat(&pSrcMatrix);
cvReleaseMat(&pProVector);
cvReleaseMat(&pProValue);
呵呵，是不是感觉有点麻烦。我想究其原因，是因为OpenCV的矩阵类是一个一维数组，静态数组的地址是连续的，因为直接使用数组的首地址就可以访问其它数组元素，而多维动态数组的行与行的地址是不连续的，因此必须通过赋值进行构造。
2. 对称方阵、特征向量输出矩阵和特征值矩阵的行列必须保证匹配，比如对称方阵是5行5列，那么特征向量输出矩阵就是5行5列，特征值矩阵就是1行5列。
3. 特征值矩阵是降序排列的（从大到小）。
有兴趣的朋友可以看看cvEigenVV的源码是如何求特征向量矩阵和特征值矩阵，我大致看了一下，是使用雅可比法的。
参考网站及文献：
1. http://www.opencv.org.cn/index.php
2. OpenCV参考手册，中文翻译者:于仕琪，中科院自动化所自由软件协会，HUNNISH，阿须数码