• QR算法求矩阵特征值的matlab实现
• 矩阵的特征值计算，用海森伯格QR算法求矩阵全部特征值
• #include #include #include #include #include <cstdlib>#define N 10using namespace std;class MT { public: int x, y; double mt[N][N]; MT() {} MT(int x
#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <cstdlib>

#define N 10

using namespace std;

class MT {
public:
int x, y;
double mt[N][N];

MT() {}
MT(int x, int y) {
this->x = x;
this->y = y;
for (int i = 0; i < this->x; i++)
for (int j = 0; j < this->y; j++)
mt[i][j] = 0.0;
}

MT transpose() {
MT c(this->y, this->x);
for (int i = 0; i < c.x; i++)
for (int j = 0; j < c.y; j++)
c.mt[i][j] = this->mt[j][i];
return c;
}

void print() {
for (int i = 0; i < this->x; i++) {
for (int j = 0; j < this->y; j++)
printf("%.4lf    ", this->mt[i][j]);
puts("");
}
}

void shift(int num) {
for (int i = this->x; i >= 0; i--)
for (int j = this->y; j >= 0; j--)
this->mt[i + num][j + num] = this->mt[i][j];
for (int i = 0; i < num; i++)
for (int j = 0; j < this->y + num; j++)
this->mt[i][j] = 0;
for (int i = 0; i < this->y + num; i++)
for (int j = 0; j < num; j++)
this->mt[i][j] = 0;
for (int i = 0; i < num; i++)
this->mt[i][i] = 1;
this->x += num;
this->y += num;
}
};

MT operator * (MT a, MT b) {
MT c(a.x, b.y);
for (int i = 0; i < c.x; i++)
for (int j = 0; j < c.y; j++) {
c.mt[i][j] = 0;
for (int k = 0; k < a.y; k++)
c.mt[i][j] += a.mt[i][k] * b.mt[k][j];
}
return c;
}

MT operator * (double a, MT b) {
MT c(b.x, b.y);
for (int i = 0; i < c.x; i++)
for (int j = 0; j < c.y; j++) {
c.mt[i][j] = a * b.mt[i][j];
}
return c;
}

MT operator + (MT a, MT b) {
MT c(a.x, a.y);
for (int i = 0; i < c.x; i++)
for (int j = 0; j < c.y; j++) {
c.mt[i][j] += a.mt[i][j] + b.mt[i][j];
}
return c;
}

MT x, H[N], I, u, e1, A;
double sita, beta;

int n;

void init() {
n = 3;
A.x = 3, A.y = 3;
A.mt[0][0] = 2;A.mt[0][1] = -2;A.mt[0][2] = 3;
A.mt[1][0] = 1;A.mt[1][1] = 1;A.mt[1][2] = 1;
A.mt[2][0] = 1;A.mt[2][1] = 3;A.mt[2][2] = -1;
}

void transform(int num) {
// initialize
x.x = n - num + 1, x.y = 1;
for (int i = num - 1; i < n; i++)
x.mt[i - num + 1][0] = A.mt[i][num - 1];

e1.x = x.x, e1.y = 1;
e1.mt[0][0] = 1;
for (int i = 1; i < e1.x; i++) {
e1.mt[i][0] = 0;
}

I.x = I.y = x.x;
for (int i = 0; i < I.x; i++)
for (int j = 0; j < I.y; j++) {
if (i == j) I.mt[i][j] = 1;
else I.mt[i][j] = 0;
}

// run householder
sita = 0;
for (int i = 0; i < x.x; i++)
sita += x.mt[i][0] * x.mt[i][0];
sita = sqrt(sita);
if (x.mt[0][0] < 0) sita = -sita;

MT e1_t = sita * e1;

u = x + e1_t;

beta = sita * (sita + x.mt[0][0]);

MT u_t = u.transpose();
u_t = u * u_t;
u_t = (-1.0 / beta) * u_t;

H[num] = I + u_t;
H[num].shift(num - 1);
printf("H%d:\n", num);
H[num].print();
printf("sita: %.4lf\n\n", sita);

A = H[num] * A;
}

int main() {
init();
int cnt = 0;
while (1) {
cnt++;
for (int i = 1; i < n; i++) {
transform(i);
}
MT R = A;

MT Q = H[n - 1];
for (int i = n - 2; i >= 1; i--)
Q = Q * H[i];
Q = Q.transpose();

puts("Q:");
Q.print();
puts("R:");
R.print();
A = R * Q;
printf("A%d:\n", cnt);
A.print();
system("pause");
}

return 0;
}
展开全文
• 矩阵求特征值特征向量Java代码1.代码介绍2.参考文章 1.代码介绍   代码中实现的模块如下： 矩阵转置 点乘 自动索引生成 矩阵复制 获取给定索引矩阵所在行 矩阵行列式 矩阵正交化 矩阵求逆 矩阵QR分解(注：QR分解...
矩阵求特征值特征向量Java代码1.代码介绍2.参考文章
1.代码介绍
代码中实现的模块如下：

矩阵转置
点乘
自动索引生成
矩阵复制
获取给定索引矩阵所在行
矩阵行列式
矩阵正交化
矩阵求逆
矩阵QR分解(注：QR分解函数只返回一个矩阵，设输入矩阵为m×n，则前m行为Q矩阵，后n行为R矩阵)
矩阵乘法
矩阵相减
求特征值(注：返回值为二维矩阵，对角线上的方为特征值)
求特征向量

所有代码如下：
注意：本文中代码所在包为util，类名为Eig，使用时注意对于及修改；代码传入参数等也请注意阅读函数说明。
package util;

import java.util.Arrays;

/**
* @(#)Eig.java The class of Eig.
* @author: 因吉
* @Email: inki.yinji@qq.com
* @Created: November 23, 2019.
* @Last_modified: November 23, 2019.
*/
public class Eig {

/**
************
* Array transpose.
*
* @param paraArray The transposing array(double[][]).
*
* @return returnArray The transposed array(double[][]).
************
*/
public static double[][] arrayTranspose(double[][] paraArray) {
int tempM = paraArray.length;
int tempN = paraArray[0].length;
double[][] returnArray = new double[tempN][tempM];
for (int i = 0; i < tempM; i++) {
for (int j = 0; j < tempN; j++) {
returnArray[j][i] = paraArray[i][j];
} // Of for j
} // Of for i

return returnArray;
}// Of arrayTranspose

/**
************
* Array corresponding elements multiply and add all elements.点积
*
* @param paraFirstArray  The first array(double[]).
* @param paraSecondArray The second array(double[]).
* @return The result(double).
************
*/
public static double arrayMultiplyAndAdd(double[] paraFirstArray, double[] paraSecondArray) {
int tempM = paraFirstArray.length;
double resultMultipliedArray = 0;

for (int i = 0; i < tempM; i++) {
resultMultipliedArray += paraFirstArray[i] * paraSecondArray[i];
} // Of for i

return resultMultipliedArray;

/**
************
* Get the index of array.
*
* @param paraLen The length of array(int).
* @return The index of array(int[]).
************
*/
public static int[] arrayIndexAuto(int paraLen) {
int[] returnArray = new int[paraLen];
for (int i = 0; i < paraLen; i++) {
returnArray[i] = i;
}
return returnArray;
}// Of arrayIndexAuto

/**
************
* Get the index of array.
*
* @param paraStrat The start index(int).
* @param paraEnd   The end index(int).
* @return The index of array(int[]).
************
*/
public static int[] arrayIndexAuto(int paraStrat, int paraEnd) {
int[] returnArray = new int[paraEnd - paraStrat];
for (int i = 0; i < returnArray.length; i++) {
returnArray[i] = i + paraStrat;
}
return returnArray;
}// Of arrayIndexAuto

/**
************
* Copy array.
*
* @param paraArray The array(double[][]).
* @return The copied array(double[][]).
************
*/
public static double[][] arrayCopy(double[][] paraArray) {
int tempM = paraArray.length;
double[][] resultArray = new double[tempM][];
for (int i = 0; i < tempM; i++) {
resultArray[i] = paraArray[i];
} // Of for i
return resultArray;
}// Of arrayCopy

/**
************
* Get the elements of the corresponding row of the incoming array.
*
* @param paraArray The import array(double[][]).
* @param paraIndex The index i(int[]).
* @return The row elements corresponding incoming array(doubel[][]).
************
*/
public static double[][] arrayRowValue(double[][] paraArray, int[] paraIndex) {
int tempParaIndex = paraIndex.length;
double[][] returnRowValue = new double[tempParaIndex][];

for (int i = 0; i < returnRowValue.length; i++) {
returnRowValue[i] = paraArray[paraIndex[i]];
} // Of for i

return returnRowValue;
}// Of arrayRowValue

/**
************
* The determinant of matrix.求解矩阵行列式
*
* @param paraMatrix The given matrix(double[][]).
*
* @return The answer of determinant.
************
*/
public static double matrixDeterminant(double[][] paraMatrix) {
int tempM = paraMatrix.length;
int tempN = paraMatrix[0].length;

if (tempM == 1) {
return paraMatrix[0][0];
} else if (tempM == 2) {
return paraMatrix[0][0] * paraMatrix[1][1] - paraMatrix[1][0] * paraMatrix[0][1];
} // Of if

double resultValue = 0;
int k;
for (int i = 0; i < tempN; i++) {
double[][] tempMatrix = new double[tempM - 1][tempN - 1];
for (int j = 0; j < tempM - 1; j++) {
k = 0;
while (k < tempN - 1) {
if (k < i) {
tempMatrix[j][k] = paraMatrix[j + 1][k];
} else {
tempMatrix[j][k] = paraMatrix[j + 1][k + 1];
} // Of if
k++;
} // Of while
} // Of for j
resultValue += (Math.pow(-1., i)) * paraMatrix[0][i] * matrixDeterminant(tempMatrix);
} // Of for i

return resultValue;
}// Of matrixDeterminant

/**
************
* The gram schimidt of matrix.求解矩阵正交化
*
* @param paraMatrix The given matrix(double[][]).
*
* @return The gram schimidt of matrix.
************
*/
public static double[][] matrixGramSchimidt(double[][] paraMatrix) {
if (paraMatrix == null) {
return null;
} // Of if

double[][] tempTransposedMatrix = arrayTranspose(paraMatrix);
int tempM = tempTransposedMatrix.length;
int tempN = tempTransposedMatrix[0].length;

double[][] resultMatrix = new double[tempM][tempN];
double tempValue = 0;
double tempFactor = 0;
for (int i = 0; i < tempM; i++) {
for (int j = 0; j < tempN; j++) {
tempValue = tempTransposedMatrix[i][j];
for (int k = 0; k < i; k++) {
tempFactor = (1. * arrayMultiplyAndAdd(tempTransposedMatrix[i], resultMatrix[k]))
tempValue -= tempFactor * resultMatrix[k][j];
} // Of for k
resultMatrix[i][j] = tempValue;
} // Of for j
} // Of for i

return arrayTranspose(resultMatrix);
}// Of matrixGramSchimidt

/**
************
* The inverse of matrix.矩阵求逆
*
* @param paraMatrix The given matrix(double[][]).
*
* @return The inverse of matrix(double[][]).
************
*/
public static double[][] matrixInverse(double[][] paraMatrix) {
if (paraMatrix == null) {
return null;
} // Of if

// Make sure that this matrix is invertible.
if (matrixDeterminant(paraMatrix) == 0) {
System.out.println("Fetal error, the matrix can not be invertible.");
System.exit(0);
} // Of if

int tempM = paraMatrix.length;
int tempN = paraMatrix[0].length * 2;

double[][] tempCopyMatrix = new double[tempM][tempN];

for (int i = 0; i < tempM; i++) {
for (int j = 0; j < tempN; j++) {
if (j < tempM) {
tempCopyMatrix[i][j] = paraMatrix[i][j];
} else {
if (j == i + tempM) {
tempCopyMatrix[i][j] = 1;
} // Of if
} // Of if
} // Of for j
} // Of for i

double tempTimes = 0;
for (int i = 0; i < tempM; i++) {
if (tempCopyMatrix[i][i] == 0) {
int j;
for (j = i + 1; j < tempM; j++) {
if (tempCopyMatrix[j][i] != 0) {
break;
} // Of if
} // Of for j

if (j != i + 1) {
for (int k = 0; k < tempN; k++) {
double temppVlaue = tempCopyMatrix[i][k];
tempCopyMatrix[i][k] = tempCopyMatrix[j][k];
tempCopyMatrix[j][k] = temppVlaue;
} // Of for k
} // Of if
} // Of if

for (int j = i + 1; j < tempM; j++) {
if (tempCopyMatrix[j][i] != 0) {
tempTimes = (tempCopyMatrix[j][i]) / tempCopyMatrix[i][i];
for (int k = i; k < tempN; k++) {
tempCopyMatrix[j][k] /= tempTimes;
tempCopyMatrix[j][k] -= tempCopyMatrix[i][k];
} // Of for k
} // Of if
} // Of for j
} // Of for i

for (int i = 0; i < tempM; i++) {
for (int j = i + 1; j < tempN / 2.; j++) {
if (tempCopyMatrix[i][j] != 0) {
tempTimes = tempCopyMatrix[i][j] / tempCopyMatrix[j][j];
for (int k = j; k < tempN; k++) {
tempCopyMatrix[i][k] -= tempTimes * tempCopyMatrix[j][k];
} // Of for k
} // Of if
} // Of for j
} // Of for i

for (int i = 0; i < tempM; i++) {
tempTimes = tempCopyMatrix[i][i];
for (int j = 0; j < tempN; j++) {
tempCopyMatrix[i][j] /= tempTimes;
} // Of for j
} // Of for i

double[][] resultMatrix = new double[tempM][tempN / 2];
for (int i = 0; i < tempM; i++) {
for (int j = 0; j < tempN / 2; j++) {
resultMatrix[i][j] = tempCopyMatrix[i][j + tempN / 2];
} // Of for j
} // Of for i

return resultMatrix;
}// Of matrixInverse

/**
************
* The QR-decomposition of matrix.矩阵QR分解
*
* @param paraMatrix  The given matrix(double[][]).
* @param paraMatrixQ The given matrix, the size echo paraMatrix(double[][]).
************
*/
public static double[][] matrixQrDecomposition(double[][] paraMatrix) {
double[][] tempOrthogonalMatrix = arrayTranspose(matrixGramSchimidt(paraMatrix));
int tempM = tempOrthogonalMatrix.length;
int tempN = tempOrthogonalMatrix[0].length;

double[][] tempMatrixQ = new double[tempM][tempN];
for (int i = 0; i < tempM; i++) {
double tempMag = magnitude(tempOrthogonalMatrix[i]);
for (int j = 0; j < tempN; j++) {
tempMatrixQ[i][j] = tempOrthogonalMatrix[i][j] / tempMag;
} // Of for j
} // Of for i

double[][] tempMatrixR = matrixMultiply(tempMatrixQ, paraMatrix);
double[][] resultSummary = new double[tempM + tempN][tempM];
for (int i = 0; i < tempN; i++) {
for (int j = 0; j < tempM; j++) {
resultSummary[i][j] = tempMatrixQ[j][i];
} // Of for j
} // Of for i

for (int i = tempN; i < resultSummary.length; i++) {
for (int j = 0; j < tempM; j++) {
resultSummary[i][j] = tempMatrixR[i - tempN][j];
} // Of for j
} // Of for i

return resultSummary;
}// Of matrixQrDecomposition

/**
************
* The magnitude.
*
* @param paraMatrix The given matrix(double[]).
************
*/
private static double magnitude(double[] paraMatrix) {
}// Of magnitude

/**
************
* The matrix multiply.
*
* @param paraFirstMatrix  The given first matrix(double[][]).
* @param paraSecondMatrix The given second matrix(double[][]).
************
*/
public static double[][] matrixMultiply(double[][] paraFirstMatrix, double[][] paraSecondMatrix) {
if (paraFirstMatrix == null || paraSecondMatrix == null) {
return null;
} // Of if

if (paraFirstMatrix[0].length != paraSecondMatrix.length) {
System.out.println("Fetal error: The size of two inputed matrix is illegally.");
System.exit(0);
} // Of if

int tempMa = paraFirstMatrix.length;
int tempNa = paraFirstMatrix[0].length;
int tempNb = paraSecondMatrix[0].length;

double[][] resultMatrix = new double[tempMa][tempNb];
for (int i = 0; i < tempMa; i++) {
for (int j = 0; j < tempNb; j++) {
double tempSum = 0;
for (int k = 0; k < tempNa; k++) {
tempSum += paraFirstMatrix[i][k] * paraSecondMatrix[k][j];
} // Of for k
resultMatrix[i][j] = tempSum;
} // Of for j
} // Of for i
return resultMatrix;
}// Of matrixMultiply

/**
************
* The eig of matrix.
*
* @param paraMatrix The given matrix(double[][]).
* @param paraIter   The maximum iteration.
* @return The eig of matrix.
************
*/
public static double[][] matrixEigValue(double[][] paraMatrix, int paraIter) {
int tempM = paraMatrix.length;
int tempN = paraMatrix[0].length;
int[] tempIndexQ = arrayIndexAuto(tempM);
int[] tempIndexR = arrayIndexAuto(tempM, tempM + tempN);
for (int i = 0; i < paraIter; i++) {
double[][] tempSummary = matrixQrDecomposition(paraMatrix);
double[][] tempMatrixQ = arrayRowValue(tempSummary, tempIndexQ);
double[][] tempMatrixR = arrayRowValue(tempSummary, tempIndexR);
paraMatrix = matrixMultiply(tempMatrixR, tempMatrixQ);
} // Of for i

return paraMatrix;
}// Of matrixEigValue

/**
************
* The eig vector of matrix.求矩阵特征向量
*
* @param paraMatrix The given matrix(double[][]).
* @param paraIter   The maximum iteration.
* @return The eig vector of matrix.
************
*/
public static double[][] matrixEigVector(double[][] paraMatrix, int paraIter) {
int tempM = paraMatrix.length;
int tempN = paraMatrix[0].length;
double[][] tempMatrix = matrixEigValue(arrayCopy(paraMatrix), paraIter);
for (int i = 0; i < tempM; i++) {
for (int j = 0; j < tempN; j++) {
if (i != j) {
tempMatrix[i][j] = 0;
}//Of if
}//Of for j
}//Of for i

return matrixInverse(matrixSubtract(paraMatrix, tempMatrix));
}//Of matrixEigVector

/**
************
* The subtract of matrix.求矩阵相减
*
* @param paraFirstMatrix The first given matrix(double[][]).
* @param paraFirstMatrix The second given matrix(double[][]).
* @return The subtract of two matrix.
************
*/
public static double[][] matrixSubtract(double[][] paraFirstMatrix, double[][] paraSecondMatrix) {
if (paraFirstMatrix == null || paraSecondMatrix == null) {
return null;
}//Of if

int tempMFirst = paraFirstMatrix.length;
int tempNFirst = paraFirstMatrix[0].length;
int tempMSecond = paraSecondMatrix.length;
int tempNSecond = paraSecondMatrix[0].length;

if (tempMFirst != tempMSecond || tempNFirst != tempNSecond) {
System.out.println("Fetal error, the inputed matrix is illegal in matrixSub(double[][], double[][])");
System.exit(0);
}//Of if

double[][] resultMatrix = new double[tempMFirst][tempNFirst];
for (int i = 0; i < tempMFirst; i++) {
for (int j = 0; j < tempNFirst; j++) {
resultMatrix[i][j] = paraFirstMatrix[i][j] - paraSecondMatrix[i][j];
}//Of for j
}//Of for i

return resultMatrix;
}//Of matrixSubtract

/**
************
* The main
************
*/
public static void main(String[] args) {
//		double[][] tempData = { { 1, 0, 0 }, { 1, 1, 0 }, { 1, 1, 1 } };
double[][] tempData = {{2,1}, {1,2}};
double[][] tempEigValue = matrixEigValue(tempData, 1000);
double[][] tempEigVector = matrixEigVector(tempData, 1000);
System.out.println("The eig value is: " + Arrays.deepToString(tempEigValue));
System.out.println("The eig vector is: " +Arrays.deepToString(tempEigVector));
}// Of main
}//Of Eig


2.参考文章

QR算法讲解
矩阵求逆、正交化等


展开全文
• QR分解矩阵特征值、特征向量   最近在看一个高光谱图像压缩...计算特征值整体思路很简单，先使用QR分解出矩阵特征值，然后利用其特征值求解具体特征值对应的特征向量，进而出矩阵的特征向量。下面是其C代码实现
       最近在看一个高光谱图像压缩算法，其中涉及到正交变换，计算正交变换时，需要对普通矩阵求其特征向量。想要在网上找一个现成的程序，可能是我百度的能力不强吧，居然真的没找见。好了废话不多说，下面进入正题。
计算特征值整体思路很简单，先使用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算法求矩阵全部特征值的基本思想是利用矩阵的QR分解通过迭代格式   将A=A1化成相似的上三角阵，从而求出矩阵A的全部特征值。  QR方法的计算步骤如下:    下面就依次进行介绍。  一. ...

QR算法求矩阵全部特征值的基本思想是利用矩阵的QR分解通过迭代格式

将A=A1化成相似的上三角阵，从而求出矩阵A的全部特征值。

QR方法的计算步骤如下:

下面就依次进行介绍。

一. 将一般矩阵化为上Hessenberg阵

1.定义

一个矩阵如果满足i>j+1时aij=0，则将这个矩阵成为上Hessenberg阵。上Hessenberg阵

的形式如下:

2. Householder变换将一般矩阵转化为上Hessnberg阵

首先，选取Householder矩阵H1，使得经过H1相似变换后的矩阵H1AH1的元素a21下面的

元素全部为0，即a31, a41, ....., am1均为0，H1取如下形式

其中 为n-1阶HouseHolder矩阵。然后选取Householder矩阵H2，使得经过H2相似变换

之后的矩阵H2(H1AH1)H2第二列中a32下面的a42, ....am2均为0。如此进行n-2次，可以构造

n-2个householder矩阵H1,H2, Hn-2，使得 Hn-2....H2H2AH1H2....Hm-2 = H(H为上Hessenberg矩阵)。

对于一个n*m的矩阵A,第col次的H可以这样构造求得(col从0开始):

其中，I为n*n的单位矩阵， v'表示矩阵v的转置， sign(x0)表示x0的符号的相反数( 当x0>0时sign=-1,当x<=0时为1),

||x||表示向量x的长度, col等于所求的上hessenberg矩阵的序号，从0开始。

二. 用Givens变换对上hessnberg矩阵作QR分解

此时有  H = R21' * R32' * ... * Rn(n-1)'R = QR。

多次计算H，直到H的变化小于一个较小的阈值时，停止迭代，此时H主对角线上的元素

即为矩阵A的全部特征值。

下面举个例子来说明求解矩阵的全部特征值的过程。

求矩阵的全部特征值

首先将A化成上hessenberg阵，取

x = [0, 6, 4], 则 ||x|| = =

则 w = [0, ,
0] ， v = w + 1 * x = [0, 6+,
4]

则 p = v*v'/v'*v =

于是 H1 = I - 2*p =

所以 H = H1AH1 =

H即为与A相似的上hessenberg矩阵。将H进行QR分解


展开全文
• 矩阵特征值的数值求解方法中，一般有部分特征值和特征向量的幂法和反幂法，以及取全部特征值QR分解方法。在此关注QR分解方法。 QR分解，是把矩阵分解为正交矩阵和三角矩阵的乘积，其中关键点在于正交矩阵的...
• QR算法求矩阵全部特征值的基本思想是利用矩阵的QR分解通过迭代格式   将A=A1化成相似的上三角阵，从而求出矩阵A的全部特征值。  QR方法的计算步骤如下:    下面就依次进行介绍。  一. ...
• java语言编写的一个求特征值和特征向量，算法是：先采用正交分解进行QR分解，然后进行迭代N次后求出对角矩阵和正交矩阵，对角矩阵即为特征值，正交矩阵即为特征向量。
• 矩阵的特征值，主要是用的QR分解，在我的有一次博客里，我已经详细地给出了计算的过程，大家有兴趣可以去看下，经过几天的钻研，终于完成了整个的eig算法。下面我将把我的整个代码附上，有不懂的可以问我，欢迎...
• 算法对图像进行NSCT变换，并在低频系数子图中提取稳定的SIFT特征点，根据特征点的变化出几何变换参数。在对几何失真图像进行校正后，进行分块QR分解。取各块中R矩阵行向量2-范数组成序列，二量化后变换为矩阵...
• 稠密矩阵特征值和特征向量的计算稠密矩阵特征值计算... 对Hessenberg阵H运用QR方法，出H的特征值，H的特征值就是A的特征值。 用逆迭代的算法，对于一个给定的特征值，计算它对应的特征向量。1.将矩阵A转化为Hessenb
• 但当采样数较少，采样协方差矩阵估计值的噪声特征值分散会导致波束形成算法的性能下降问题，QR算法的性能就会下降。针对此缺陷，提出了对角加载奇异值（DSVD）分解的算法，该算法先对采样数据所构成的矩阵进行重构、...
• ## QR分解与最小二乘

千次阅读 2018-08-23 21:51:00
主要内容： 1、QR分解定义 2、QR分解法 3、QR分解与最小二乘 4、Matlab实现 ...QR 分解也是特定特征值算法即QR算法的基础。 定义： 实数矩阵 A 的 QR 分解是把 A 分解为Q、R，这里的 Q ...
• hessqrtz 海森伯格QR算法求矩阵全部特征值 rqrtz 瑞利商位移QR算法求矩阵全部特征值 第7章： 数值微分 第8章： 数值积分 第9章： 方程求根 .............................. .............. 第17章： 数据统计和...
• ## MATLAB常用算法

热门讨论 2010-04-05 10:34:28
hessqrtz 海森伯格QR算法求矩阵全部特征值 rqrtz 瑞利商位移QR算法求矩阵全部特征值 第7章： 数值微分 MidPoint 中点公式求取导数 ThreePoint 三点法求函数的导数 FivePoint 五点法求函数的导数 DiffBSample 三次...
• 使用QR迭代也可以出矩阵的特征值，不过前一步要求矩阵变为上Hessenberg矩阵。将矩阵进行QR分解可以使用Household变换，Schmidt正交化以及Given变换。对于Schmidt不是十分了解。QR迭代的核心是把矩阵...
• hessqrtz 海森伯格QR算法求矩阵全部特征值 rqrtz 瑞利商位移QR算法求矩阵全部特征值 第7章： 数值微分 MidPoint 中点公式求取导数 ThreePoint 三点法求函数的导数 FivePoint 五点法求函数的导数 DiffBSample 三次...
• ## pca算法的实现

千次阅读 2012-12-06 13:02:06
按照pca的处理步骤，终于用matlab把pca的算法写好了，本来想用C++写的，中间遇到了一个问题：就是在求特征值和特征向量的时候比较麻烦，需要QR迭代算法。自己就投机取巧，matlab里只需要调用一个函数就可以了。想着...
• 提示：三阶QR 分解计算实例————保姆级分解过程 文章目录前言一、QR分解理论推导二、案例解析1.对A=(1112−1−12−45)\begin{pmatrix} 1&...QR分解法求解特征值应用范围广泛，是本世界十大算法之一，可以求解任
• 6-1 试验目的计算特征值实现算法 试验容随机产生一个 10 阶整数矩阵各数均在 -5 和 5 之间 1 用 MATLAB 函数 eig 求矩阵全部特征值 2 用幂法求 A 的主特征值及对应的特征向量 3 用基本 QR 算法求全部特征值可用 ...
• 3.17 赫申伯格矩阵全部特征值QR方法 3.18 求实对称矩阵特征值与特征向量的雅可比法 3.19 求实对称矩阵特征值与特征向量的雅可比过关法 第4章 线性代数方程组的求解 4.1 线性方程组类设计 4.2 全选主元高斯消去法...
• 5.实赫申伯格矩阵的QR算法 第7章数据拟合 1.直线拟合 2.线性最小二乘法 3.非线性最小二乘法 4.绝对值偏差最小的直线拟合 第8章方程根和非线性方程组的解法 1.图解法 2.逐步扫描法和二分法 3.割线法和试位...
• 5.实赫申伯格矩阵的QR算法 第7章数据拟合 1.直线拟合 2.线性最小二乘法 3.非线性最小二乘法 4.绝对值偏差最小的直线拟合 第8章方程根和非线性方程组的解法 1.图解法 2.逐步扫描法和二分法 3.割线法和试位法 4....