基于线性代数的人脸识别初级系统

基于线性代数的人脸识别初级系统

声明：本文BugMaker14与Mid2dog为并列第一作者，同组代码实现与注释标写。
写代码的原因：作者为计算机系学生。某日线性代数课堂学习矩阵的特征值，话题莫名转到“线性代数有什么用”。老师当场解释了线性代数在人脸识别中的运用，并称谁能按照刚才所说的写出C语言代码，就给谁免考。Mid2dog同学由于第一次单元测试分数比较少，有反向被免考的概率。为了阻止挂科，他拉了BugMaker14做一组，一起写了本代码，前后共耗时六周（主要是快期末了，有很多课的论文都在催）。

1、数学基础

一张照片由mn个像素点构成，彩色像素可以通过公式转换为灰度像素，而灰度像素又可以通过[0,2^n)的数值表示，所以一张照片就可以通过mn的矩阵来表达

1. 概念：矩阵特征值

设 A 是n阶方阵，如果存在数m和非零n维列向量 x，使得 Ax=mx 成立，则称 m 是矩阵A的一个特征值（characteristic value)或本征值（eigenvalue)。
——百度百科

我们在计算矩阵特征值的时候，直接地通过数学的方法获得高阶方程并对其求解显然是不能实现的。因为在方程高于六阶时，不像二阶方程的[-b+sqrt(b^2-4ac)]/2a，没有直接的公式来获得答案。所以我们的矩阵应该变形。一个矩阵的特征值最直接的获得就是当这个矩阵是上三角或者下三角矩阵时，这样的矩阵特征值就是对角线上的值。

2. 运算：矩阵的初等相似变换
   不改变矩阵特征值的变换方法是矩阵的初等相似变换。
   

当然这里就有朋友要问了：“那我们是不是可以直接通过矩阵的初等相似变换获得上三角矩阵？”答案显然 是否定的。
高阶方程的一大特点就是基本没有全实数解，多少都会出现成对的虚数解。但是这个方法你细品，毫无虚数解的可能性，所以从数理上肯定是有问题的。至于为什么我不知道。
但是有个东西叫海森堡阵 [滑稽] ，这种方法获得海森堡阵还是没问题的。
数学原理我还是不知道 获得海森堡阵时每次都交换最大的那个行列能最大的减小误差。

3. QR分解

QR算法求矩阵全部特征值的基本思想是利用矩阵的QR分解通过迭代格式将A=A1化成相似的上三角阵，从而求出矩阵A的全部特征值。
求出海森堡阵的形式后，然后接下来有两种方法：一种是given变换，第二种是householder变换，这里我们采用了这个方法。

Mid2dog同学不知道哪里找来的一堆奇奇怪怪的术语，看的我也头大。

____现在有某一个向量 x， 想通过一个标准正交矩阵P将 x 转换为 y，有什么方法可以求出矩阵 P？一种方法是通过上面的givens旋转一步一步完成，P = G1G2G3。且givens变换是值得信赖的正交阵，不会改变向量长度。但还有一个更加快捷的公式，即为 Householder Transformation（豪斯霍尔德变换）。
____这一变换将一个向量变换为由一个超平面反射的镜像，是一种线性变换。其变换矩阵被称作豪斯霍尔德矩阵，在一般内积空间中的类比被称作豪斯霍尔德算子。超平面的法向量被称作豪斯霍尔德向量。
____如果V给出为单位向量而I是单位矩阵，则描述上述线性变换的是豪斯霍尔德矩阵（v*表示向量v的共轭转置）

还配了几张图
下面的图都是引用的，但是时间有点久了都忘了哪里引用了，只是作为数学原理资料，勿喷谢谢。




水印出处都还在，这图只是作为资料，勿喷。

但是我都没看懂，那能怎么办呢代码还是要写的，不然要挂科了， 于是就找了一张非常有爱的图。
他能帮助你理解代码。但这个不是我们写的。

2、上述代码实现

/*finish*/void hsb(double a[][Max])//初等变换成海森堡阵     
{
   
   
    int i,j,k;
    int n = Max;//max是预置的方阵阶数
    double d,t;
	for (k = 1; k <= n - 2; k++) 
    {
   
   
        d=0.0;
		for (j = k; j <= n - 1; j++) //取最大交换
        {
   
   
            t=a[j][k-1];
            if (fabs(t)>fabs(d)) 
			{
   
    
				d = t;
				i = j;
			}
        }
        if (fabs(d) != 0)
        {
   
   
            if (i != k)//初等相似变换
            {
   
   
				for (j = k - 1; j <= n - 1; j++) 
                {
   
   
					t = a[i][j]; 
					a[i][j] = a[k][j]; 
					a[k][j] = t;
                }
				for (j = 0; j <= n - 1; j++) 
                {
   
   
					t = a[j][i]; 
					a[j][i] = a[j][k]; 
					a[j][k] = t;
                }
            }
			for (i = k + 1; i <= n - 1; i++) 
            {
   
   
				t = a[i][k - 1] / d; a[i][k - 1] = 0.0;//消去最大值并完成化零
				for (j = k; j <= n - 1; j++)
					a[i][j] = a[i][j] - t * a[k][j];
				for (j = 0; j <= n - 1; j++)
					a[j][k] = a[j][k] + t * a[j][i];
            }
        }
    }
}

↑这样我们就能获得海森堡阵了。

/*finish*/void Householder(double M[][Max])//豪斯霍尔德法 
{
   
   
	double B, temp;   
	double u[2];     
	double v[2];	   
	double A[3][Max - 1]; 
	for (int h = 0; h < Max - 1; h++) //MAX-1次变换 
	{
   
   
		v[0] = M[h][h];			   //因为已经是海森堡矩阵了，所以每个列向量只有俩是非零，就只用取这俩 
		v[1] = M[h + 1][h];		   
		B = BBB(u, v);             //求出beta，也就是u向量 		
		A[0][h] = u[0];   		   //直接将u向量的元素u0，u1这些元素赋值给A矩阵的相应的位置，然后便于进行之后A矩阵左乘
		A[1][h] = u[1];			   
		A[2][h] = B;   		  		  
		for (int j = 0; j < Max; j++)  //因为只用计算两个，所以实际上只用取两个累加即可    
		{
   
   							        
			temp = (u[0] * M[h][j] + u[1] * M[h + 1][j]) / B;  //相加后除以系数
			M[h][j] -= temp * u[0];//递推过程中一个点加两次，所以不再乘2
			M[h + 1][j] -= temp * u[1];
		}
	}
	for (int j = 0; j < Max - 1; j++) //左乘A,大致与上面一样 
	{
   
   
		u[0] = A[0][j]; 
		u[1] = A[1][j];  
		  B  = A[2][j];  
		for (int h = 0; h < Max; h++)//行化简
		{
   
   
			temp = u[0] * M[h][j] + u[1] * M[h][j + 1];
			temp /= B;
			M[h][j] -= temp*u[