QR分解求矩阵特征值、特征向量 C语言
2024-05-15 10:23:55
最近在看一个高光谱图像压缩算法，其中涉及到正交变换，计算正交变换时，需要对普通矩阵求其特征向量。想要在网上找一个现成的程序，可能是我百度的能力不强吧，居然真的没找见。好了废话不多说，下面进入正题。
计算特征值整体思路很简单，先使用QR分解求出矩阵特征值，然后利用其特征值求解具体特征值对应的特征向量，进而求出矩阵的特征向量。下面是其C代码实现：
//------------------------------------头文件---------------------------------
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>//--------------------------这里是一些定义的结构体和数据类型---------
//纯C里面定义的布尔类型
typedef enum{False = 0,True = 1}Bool;
//定义矩阵元素的类型为matrixType
typedef double matrixType;//此结构体用来表示矩阵，其中row为行，column为列，height为高，array用来存放矩阵元素(用一维来模拟二维/三维)
typedef struct
{
unsigned row,column,height;
matrixType *array; //使用时，必须对*array进行初始化
}Matrix;//---------下面是QR分解，求解线性方程所用到的一些函数-----------
/*
matrix为要设置大小并分配内存的矩阵，row、column、height分别为行，列，高。
函数调用成功则则返回true,否则返回false
*/
Bool SetMatrixSize(Matrix *matrix ,const unsigned row,const unsigned column,const unsigned height)
{
unsigned size = row * column * height * sizeof(matrixType);
if(size <= 0 )
{
return False;
}matrix->array = (matrixType*)malloc(size);//如果分配内存成功
if(matrix->array)
{
matrix->row = row;
matrix->column = column;
matrix->height = height;
return True;
}
else
{
matrix->row = matrix->column = matrix->height = 0;
return False;
}
}//设置Matrix矩阵中的所有元素大小为ele
void SetMatrixEle(Matrix *matrix,matrixType ele)
{
unsigned size = matrix->row * matrix->column * matrix->height;
unsigned i;for(i = 0;i < size;++i)
{
matrix->array[i] = ele;
}
}//设置Matrix矩阵中的所有元素大小为0
void SetMatrixZero(Matrix*matrix)
{
SetMatrixEle(matrix,0);
}//判断矩阵是否为空，若为空则返回1，否则返回0
Bool IsNullMatrix(const Matrix* matrix)
{
unsigned size = matrix->row * matrix->column * matrix->column;if(size <= 0 || matrix->array == NULL)
{
return True;
}
return False;
}//销毁矩阵，即释放为矩阵动态分配的内存,并将矩阵的长宽高置0
void DestroyMatrix(Matrix *matrix)
{
if(!IsNullMatrix(matrix))
{
free(matrix->array);
matrix->array = NULL;
}matrix->row = matrix->column = matrix->height = 0;
}//计算矩阵可容纳元素个数，即return row*column*height
unsigned MatrixCapacity(const Matrix*matrix)
{
return matrix->row * matrix->column * matrix->height;
}//||matrix||_2 求A矩阵的2范数
matrixType MatrixNorm2(const Matrix *matrix)
{
unsigned size = matrix->row * matrix->column *matrix->height;
unsigned i;
matrixType norm = 0;for(i = 0;i < size;++i)
{
norm += (matrix->array[i]) *(matrix->array[i]);
}return (matrixType)sqrt(norm);
}//matrixB = matrix(:,:,height)即拷贝三维矩阵的某层，若matrix为二维矩阵，需将height设置为0
Bool CopyMatrix(Matrix* matrixB,Matrix *matrix,unsigned height)
{
unsigned size,i;
Matrix matrixA;//判断height值是否正确
if(height < 0 || height >= matrix->height)
{
printf("ERROR: CopyMatrix() parameter error！\n");
return False;
}//将matrix(:,:,height) 转换为二维矩阵matrixA
matrixA.row = matrix->row;
matrixA.column = matrix->column;
matrixA.height = 1;
matrixA.array = matrix->array + height * matrix->row * matrix->column;//判断两矩阵指向的内存是否相等
if(matrixA.array == matrixB->array)
{
return True;
}//计算matrixA的容量
size = MatrixCapacity(&matrixA);
//判断matrixB的容量与matrixA的容量是否相等
if( MatrixCapacity(matrixB)!= size)
{
DestroyMatrix(matrixB);
SetMatrixSize(matrixB,matrixA.row,matrixA.column,matrixA.height);
}
else
{
matrixB->row = matrixA.row;
matrixB->column = matrixA.column;
matrixB->height = matrixA.height;
}for(i = 0;i < size;++i)
{
matrixB->array[i] = matrixA.array[i];
}return True;
}//matrixC = matrixA * matrixB
Bool MatrixMulMatrix(Matrix *matrixC,const Matrix* matrixA,const Matrix* matrixB)
{
size_t row_i,column_i,i;
size_t indexA,indexB,indexC;
matrixType temp;
Matrix tempC;if(IsNullMatrix(matrixA) || IsNullMatrix(matrixB))
{
return False;
}if(matrixA->column != matrixB->row )
{
return False;
}if(MatrixCapacity(matrixC) != matrixA->row * matrixB->column)
{
SetMatrixSize(&tempC,matrixA->row,matrixB->column,1);
}
else
{
tempC.array = matrixC->array;
tempC.row = matrixA->row;
tempC.column = matrixB->column;
tempC.height = 1;
}for(row_i = 0;row_i < tempC.row;++row_i)
{
for(column_i = 0;column_i < tempC.column;++column_i)
{
temp = 0;
for(i = 0;i < matrixA->column;++i)
{
indexA = row_i * matrixA->column + i;
indexB = i * matrixB->column + column_i;temp += matrixA->array[indexA] * matrixB->array[indexB];
}
indexC = row_i * tempC.column + column_i;tempC.array[indexC] = temp;
}
}if(tempC.array != matrixC->array)
{
DestroyMatrix(matrixC);matrixC->array = tempC.array;
}matrixC->row = tempC.row;
matrixC->column = tempC.column;
matrixC->height = tempC.height;return True;
}//对vector中所有元素排序，sign= 0 时为升序，其余为降序
void SortVector(Matrix *vector,int sign)
{
matrixType mid;
int midIndex;
int size = MatrixCapacity(vector);
int i,j;if(0 == sign)
{
for(i = 0;i < size;++i)
{
mid = vector->array[i];
midIndex = i;
for( j = i + 1; j < size ; ++j)
{
if(mid > vector->array[j])
{
mid = vector->array[j];
midIndex = j;
}
}vector->array[midIndex] = vector->array[i];
vector->array[i] = mid;
}
}
else
{
for(i = 0;i < size;++i)
{
mid = vector->array[i];
midIndex = i;
for( j = i + 1; j < size ; ++j)
{
if(mid < vector->array[j])
{
mid = vector->array[j];
midIndex = j;
}
}vector->array[midIndex] = vector->array[i];
vector->array[i] = mid;
}
}
}//打印矩阵
void PrintMatrix(const Matrix *matrix)
{
size_t row_i,column_i,height_i,index;for(height_i = 0;height_i < matrix->height;++height_i)
{
(matrix->height == 1) ? printf("[:,:] = \n"):printf("[%d,:,:] = \n",height_i);for(row_i = 0;row_i < matrix->row;++row_i)
{
for(column_i = 0;column_i < matrix->column;++column_i)
{
index = height_i * matrix->row * matrix->column + row_i * matrix->column + column_i;
printf("%12.4g",matrix->array[index]);
}
printf("\n");
}
}
}//----------------------QR分解-------------------------------------------//将A分解为Q和R
void QR(Matrix *A,Matrix *Q,Matrix *R)
{
unsigned i,j,k,m;
unsigned size;
const unsigned N = A->row;
matrixType temp;Matrix a,b;//如果A不是一个二维方阵，则提示错误，函数计算结束
if(A->row != A->column || 1 != A->height)
{
printf("ERROE: QR() parameter A is not a square matrix!\n");
return;
}size = MatrixCapacity(A);
if(MatrixCapacity(Q) != size)
{
DestroyMatrix(Q);
SetMatrixSize(Q,A->row,A->column,A->height);
SetMatrixZero(Q);
}
else
{
Q->row = A->row;
Q->column = A->column;
Q->height = A->height;
}if(MatrixCapacity(R) != size)
{
DestroyMatrix(R);
SetMatrixSize(R,A->row,A->column,A->height);
SetMatrixZero(R);
}
else
{
R->row = A->row;
R->column = A->column;
R->height = A->height;
}SetMatrixSize(&a,N,1,1);
SetMatrixSize(&b,N,1,1);for(j = 0 ; j < N;++j)
{
for(i = 0;i < N; ++i)
{
a.array[i] = b.array[i] = A->array[i * A->column + j];
}for(k = 0; k < j; ++k)
{
R->array[k * R->column + j] = 0;for(m = 0;m < N; ++m)
{
R->array[k * R->column + j] += a.array[m] * Q->array[m * Q->column + k];
}for(m = 0;m < N; ++m)
{
b.array[m] -= R->array[k * R->column + j] * Q->array[m * Q->column + k];
}
}temp = MatrixNorm2(&b);
R->array[j * R->column + j] = temp;for(i = 0; i < N; ++i)
{
Q->array[i * Q->column + j] = b.array[i] / temp;
}
}DestroyMatrix(&a);
DestroyMatrix(&b);
}//----------------------使用特征值计算矩阵特征向量-----------------
//eigenVector为计算结果即矩阵A的特征向量
//eigenValue为矩阵A的所有特征值，
//A为要计算特征向量的矩阵
void Eigenvectors(Matrix *eigenVector, Matrix *A,Matrix *eigenValue)
{
unsigned i,j,q;
unsigned count;
int m;
const unsigned NUM = A->column;
matrixType eValue;
matrixType sum,midSum,mid;
Matrix temp;SetMatrixSize(&temp,A->row,A->column,A->height);for(count = 0; count < NUM;++count)
{
//计算特征值为eValue，求解特征向量时的系数矩阵
eValue = eigenValue->array[count] ;
CopyMatrix(&temp,A,0);
for(i = 0 ; i < temp.column ; ++i)
{
temp.array[i * temp.column + i] -= eValue;
}//将temp化为阶梯型矩阵
for(i = 0 ; i < temp.row - 1 ; ++i)
{
mid = temp.array[i * temp.column + i];
for(j = i; j < temp.column; ++j)
{
temp.array[i * temp.column + j] /= mid;
}for(j = i + 1;j < temp.row;++j)
{
mid = temp.array[j * temp.column + i];
for(q = i ; q < temp.column; ++q)
{
temp.array[j * temp.column + q] -= mid * temp.array[i * temp.column + q];
}
}
}
midSum = eigenVector->array[(eigenVector->row - 1) * eigenVector->column + count] = 1;
for(m = temp.row - 2; m >= 0; --m)
{
sum = 0;
for(j = m + 1;j < temp.column; ++j)
{
sum += temp.array[m * temp.column + j] * eigenVector->array[j * eigenVector->column + count];
}
sum= -sum / temp.array[m * temp.column + m];
midSum += sum * sum;
eigenVector->array[m * eigenVector->column + count] = sum;
}midSum = sqrt(midSum);
for(i = 0; i < eigenVector->row ; ++i)
{
eigenVector->array[i * eigenVector->column + count] /= midSum;
}
}
DestroyMatrix(&temp);
}
int main()
{
const unsigned NUM = 50; //最大迭代次数unsigned N = 3;
unsigned k;Matrix A,Q,R,temp;
Matrix eValue;//分配内存
SetMatrixSize(&A,N,N,1);
SetMatrixSize(&Q,A.row,A.column,A.height);
SetMatrixSize(&R,A.row,A.column,A.height);
SetMatrixSize(&temp,A.row,A.column,A.height);
SetMatrixSize(&eValue,A.row,1,1);//A设置为一个简单矩阵
A.array[0] = -1;
A.array[1] = 2;
A.array[2] = 1;
A.array[3] = 2;
A.array[4] = 4;
A.array[5] = -1;
A.array[6] = 1;
A.array[7] = 1;
A.array[8] = -6;//拷贝A矩阵元素至temp
CopyMatrix(&temp,&A,0);//初始化Q、R所有元素为0
SetMatrixZero(&Q);
SetMatrixZero(&R);
//使用QR分解求矩阵特征值
for(k = 0;k < NUM; ++k)
{
QR(&temp,&Q,&R);
MatrixMulMatrix(&temp,&R,&Q);
}
//获取特征值，将之存储于eValue
for(k = 0;k < temp.column;++k)
{
eValue.array[k] = temp.array[k * temp.column + k];
}//对特征值按照降序排序
SortVector(&eValue,1);//根据特征值eValue，原始矩阵求解矩阵特征向量Q
Eigenvectors(&Q,&A,&eValue);//打印特征值
printf("特征值：");
PrintMatrix(&eValue);//打印特征向量
printf("特征向量");PrintMatrix(&Q);
DestroyMatrix(&A);
DestroyMatrix(&R);
DestroyMatrix(&Q);
DestroyMatrix(&eValue);
DestroyMatrix(&temp);return 0;
}//----------------------------代码运行结果--------------------------------
特征值：[:,:] =
4.627
-1.522
-6.105
特征向量[:,:] =
0.3514 0.9389 -0.244
0.9285 -0.3148 0.1432
0.1204 0.1394 0.9591
//---------Matlab 求的特征值，特征向量------------------------------------
D =
4.6272 0 0
0 -1.5222 0
0 0 -6.1051
V =
-0.3514 -0.9389 -0.2440
-0.9285 0.3148 0.1432
-0.1204 -0.1394 0.9591
//-------------------------代码运行软件环境--------------------------
Windows 7 旗舰版
Vs 2012(在其中创建.c文件，以C程序运行)
QR分解求矩阵特征值、特征向量 C语言相关推荐
- QR分解求矩阵全部特征值
QR算法求矩阵全部特征值的基本思想是利用矩阵的QR分解通过迭代格式 将A=A1化成相似的上三角阵,从而求出矩阵A的全部特征值. QR方法的计算步骤如下: 下面就依次进行介绍. 一. 将一般矩阵化为上H ...
- (转)QR分解求矩阵的全部特征值
QR算法求矩阵全部特征值的基本思想是利用矩阵的QR分解通过迭代格式 将A=A1化成相似的上三角阵,从而求出矩阵A的全部特征值. QR方法的计算步骤如下: 下面就依次进行介绍. 一. 将一般矩阵化为上H ...
- 【matlab】 QR分解 求矩阵的特征值
"QR_H.m"function [Q,R] = QR_tao(A) %输入矩阵A %输出正交矩阵Q和上三角矩阵R [n,n]=size(A); E = eye(n); X = z ...
- QR分解求矩阵绝对值-基于HouseHolder变换
思路: 输入矩阵A(mxn)-->HouseHolder变换-->获得矩阵B(Hessenberg矩阵nxn)-->Gievens变换 -->获得Q(标准正交nxn)和R(上三 ...
- C语言通过QR分解计算矩阵的特征值和特征向量
#include <stdio.h> #include <stdlib.h> #include <math.h> #include <time.h>// ...
- 双步位移求解特征值matlab,数值分析——带双步位移的QR分解求特征值算法
C语言实现数值分析中带双步位移的QR分解求特征值算法. 数 值 分 析(B) 大 作 业(二) 1.算法设计: ①矩阵的拟上三角化: 对实矩阵A进行相似变换化为拟上三角矩阵A(n 1),其变换矩阵采用 ...
- 求矩阵特征值和特征向量
求矩阵特征值和特征向量的一个小程序 代码较长,如果不能执行,就是要建立结构体,大家试试吧,希望能用. // // 实对称三对角阵的全部特征值与特征向量的计算 // // 参数: // 1. doubl ...
- qr分解求线性方程组_矩阵分解
矩阵分解 1. 矩阵的三角分解 1.1 高斯消去法解线性方程组 在聊起矩阵分解之前,先看一下我们小学二年级就学过的高斯消去法解线性方程组,其主要思想就是:将方程组写作写作增广矩阵(A|b)的形式,然后 ...
- 【机器学习】用QR分解求最小二乘法的最优闭式解
[机器学习]用QR分解求最小二乘法的最优闭式解 写在前面 QR分解 定义 QR的求解 线性回归模型 用QR分解求解最优闭式解 矩阵的条件数 实验 运行结果 写在前面 今天刷知乎,看到张皓在面试官如何判 ...
最新文章
- IoU、GIoU、DIoU、CIoU损失函数
- js在post后台接口的时候，一行代码完成删除对象中所有值为null、undefined或为空字符串““的属性
- 标准纯C++实现简单的词法分析器(三)
- Web UI设计的关键要素！
- 详解实时查看网卡流量的几款工具
- Android中的AsyncTask异步任务的简单介绍
- c语言中二叉树中总结点,C语言二叉树的三种遍历方式的实现及原理
- 卷积神经网络的几种典型架构
- 可靠性测试竟如此容易
- 凸优化第九章无约束优化 9.2下降方法
- Locust接口压力测试
- Pygame教程（非常详细）
- div+css静态网页设计 web网页设计实例作业 ——中国水墨风的小学学校网站(6页) 专题网页设计作业模板 学校物静态HTML网页模板下载
- pxe kickstart无人值守自动化装机
- 嗅探软件和网络测试,新鲜！山东首条燃气嗅探犬“上岗”，通检测探漏样样精...
- 跟着我干你技术入股，当面临这样的诱惑，我们该怎么办？
- 抖音seo矩阵系统，抖音矩阵系统源码怎么搭建？
- 小程序转发二维码携带参数不生效的问题
- 查看twitter浏览记录_您可以看到谁查看了您的Twitter个人资料吗？
- 如何在SQL Server 2005中修复损坏的数据库