基于QR分解迭代求解方阵特征值和特征向量( 三 ) _矩阵

template//利用初代变换计算QR分解，并且通过参数表中的引用返回bool Matrix::getQR(Matrix& Q, Matrix& R) {/*if (this->getRank() != m_Columns) {cout << "列不满秩，本算法无法进行QR分解！\n";return false;}*/Matrix a = *this;//调用transpose会使得*this变成原来矩阵的转置，必须要把原来的矩阵保存到aMatrix aTa = (this->transpose()) * a;//aTa为当前矩阵的转置乘以它本身，得到的对称正定的矩阵Matrix identity;identity = identity.cerateAnIdentity(a.m_Columns);//得到与上述aTa等大的单位阵vector tempVec;for (int i = 0; i < a.m_Columns; i++) {tempVec.clear();identity.getOneRow(i, &tempVec);aTa.addOneRowAtBottom(tempVec);}//将单位矩阵拼接到aTa矩阵的下方for (int i = 0; i < a.m_Columns - 1; i++) {for (int j = i + 1; j < a.m_Columns; j++) {double scalingFactor = aTa.m_Matrix[j][i] / aTa.m_Matrix[i][i];for (int k = 0; k < a.m_Columns; k++) {aTa.m_Matrix[j][k] -= scalingFactor * aTa.m_Matrix[i][k];}for (int k = 0; k < 2 * a.m_Columns; k++) {aTa.m_Matrix[k][j] -= scalingFactor * aTa.m_Matrix[k][i];}}}//利用初等行变换和初等列变换，把aTa矩阵的上半部分化成对角矩阵Matrix up, low, upInverse, lowInverse, tempMatrix;for (int i = 0; i < 2 * a.m_Columns; i++) {tempVec.clear();aTa.getOneRow(i, &tempVec);if (i < a.m_Columns)up.addOneRowAtBottom(tempVec);elselow.addOneRowAtBottom(tempVec);}//把aTa的上半部分方阵和下半部分方阵分别提取出来，单独储存for (int i = 0; i < a.m_Columns; i++) {if (up.m_Matrix[i][i] < 0) {cout << "出现复数特征值，本算法无法处理，以下特征值可能出错，请您知晓！\n";}up.m_Matrix[i][i] = sqrt(up.m_Matrix[i][i]);}//我们需要得到上半部分对角矩阵的开方值/*由于此处需要计算up，low两个矩阵的逆矩阵，但是我们的求逆算法在矩阵类的派生类方阵类中，此处无法使用，只能重写*/tempMatrix = up | identity;tempMatrix = tempMatrix.rref();//利用rref求矩阵的逆for (int i = 0; i < a.m_Columns; i++) {tempVec.clear();for (int j = a.m_Columns; j < 2 * a.m_Columns; j++) {tempVec.push_back(tempMatrix.m_Matrix[i][j]);}upInverse.addOneRowAtBottom(tempVec);}tempMatrix = low | identity;tempMatrix = tempMatrix.rref();for (int i = 0; i < a.m_Columns; i++) {tempVec.clear();for (int j = a.m_Columns; j < 2 * a.m_Columns; j++) {tempVec.push_back(tempMatrix.m_Matrix[i][j]);}lowInverse.addOneRowAtBottom(tempVec);}//以上部分只是完成了求逆的运算Q = a * low * upInverse;//Q、R分别的计算公式R = up * lowInverse;return true;}
当然上述代码无法直接编译通过，那是因为这是整个矩阵计算器实现的一个小模块，其中譬如”( )“的函数，重载过的乘法*运算符等等都在其他地方实现了，需要完整的代码请看文尾链接。
下面是例如上述”getQR( )“实现特征值和特征向量计算的代码：
template//计算方阵的特征值和特征向量bool SquareMatrix::getEigenvalue_vector(SquareMatrix& diag, SquareMatrix& Q) {Matrix tempA = *this;//这是一个临时的矩阵，用来保存每一次被QR分解迭代的对象Matrix tempR, tempQ, totalQ;totalQ = totalQ.cerateAnIdentity(size);//totalQ初始化为单位阵，用来累乘tempQconst int iterations = 200;//默认的迭代次数for (int i = 0; i < iterations; i++) {tempA.getQR(tempQ, tempR);tempA = tempR * tempQ;//下一次迭代的矩阵由RQ = Q'AQ给出totalQ = totalQ * tempQ;}for (int i = 0; i < size; i++) {for (int j = 0; j < size; j++) {if (i == j)continue;tempA.setAnElement(i, j, 0);}}//将对角阵的非对角元写成0//totalQ = totalQ.transpose();SquareMatrix reQ(totalQ);SquareMatrix rediag(tempA);Q = reQ;diag = rediag;if (!this->isSymmetric())cout << "由于原矩阵不是对称矩阵，特征向量不唯一，以下只给出一组可能的特征向量！\n";return true;}
同样的，其中部分函数都是实现矩阵计算机过程中另外写的模块，限于篇幅不给出具体的代码。值得注意的是，上述代码QR迭代次数为200次，在大部分情况下，精度和运算时间都能够达到比较令人满意的程度。下面是应用这个原理的一个实际计算例子，对一个五阶元素都在几千的矩阵：