一、前言
- 本章开始,笔者将详细讲解矩阵的QR分解、特征值、特征向量等运算,并给出函数和测试demo的C代码,以及与matlab计算的结果;
- 并且,本章相当于前面几章的大杂烩,前面所有的结构体、宏定义、函数本章基本全部都有用到。
- 此外,由于本章求解QR矩阵/特征值/特征向量,实数和复数的过程原理完全一样,读者可以自行修改参数,简化复数计算的部分即可,或者完全不用修改,直接调用我的函数,虚部设置为0退化为实数运算即可;
- 另,这将也是我此连载文章的系列精品之一,矩阵运算的功能越来越复杂,但是笔者全部分模块进行,化繁为简,化整为零,所以需要调用的小函数越来越多,读者遇到不熟悉的函数可以随时Jump到(一)、(二)、(三)、(四)的内容。
(PS:另外这两个编辑器为啥功能不能优势互补,而且都没有贴matlab code的编辑器?)
二、QR分解
QR分解的数学原理读者可以直接点击链接复习一下,QR分解C语言实现我参考的是此大佬的。笔者这里整理出来了,C code直接给出:
/* QR Decompose */
void QR(const Matrix* A, Matrix* Q, Matrix* R)
{
if (IsNullComplexMatrix(A))
{
puts("ERROE: An invalid matrix!\n");
return;
}
// Not A Square Matrix
if (A->row != A->column)
{
printf("ERROE: Not a square matrix!\n");
return;
}
int i, j, k, m;
int size;
const int N = MatrixRow(A);
ComplexType temp;
// Column Vector Saving Column Vectors of A
Matrix a, b;
InitComplexMatrix(&a, N, 1);
InitComplexMatrix(&b, N, 1);
size = MatrixSize(A);
if (MatrixSize(Q) != size)
{
// free(Q->arrayComplex);
DestroyComplexMatrix(Q); // Free The Initial Matrix
InitComplexMatrix(Q, A->row, A->column); // Reset Size and Initialize
}
if (MatrixSize(R) != size)
{
// free(R->arrayComplex);
DestroyComplexMatrix(R);
InitComplexMatrix(R, A->row, A->column);
}
for (j = 0; j < N; ++j)
{
for (i = 0; i < N; ++i)
{
a.arrayComplex[i] = b.arrayComplex[i] = A->arrayComplex[i * A->column + j]; // Cols Vector of A
}
for (k = 0; k < j; ++k)
{
R->arrayComplex[k * R->column + j]._Val[0] = 0;
R->arrayComplex[k * R->column + j]._Val[1] = 0;
// InitComplex(R->arrayComplex + k * R->column + j); // reset
for (m = 0; m < N; ++m)
{
R->arrayComplex[k * R->column + j] = AddComplex(R->arrayComplex[k * R->column + j], \
_Cmulcc(a.arrayComplex[m], Q->arrayComplex[m * Q->column + k]));
}
for (m = 0; m < N; ++m)
{
b.arrayComplex[m] = SubComplex(b.arrayComplex[m], _Cmulcc(R->arrayComplex[k * R->column + j], \
Q->arrayComplex[m * Q->column + k]));
}
}
temp = MatrixNorm2(&b);
R->arrayComplex[j * R->column + j] = temp;
for (i = 0; i < N; ++i)
{
Q->arrayComplex[i * Q->column + j] = DivComplex(b.arrayComplex[i], temp);
}
}
// Free Local Memory
Matrix ComplexMatrixArray[] = {
a, b };
int numDoubleMatrixArray = sizeof(ComplexMatrixArray) / sizeof(Matrix);
DestroyComplexMatrixArray(ComplexMatrixArray, numDoubleMatrixArray); // Complpex Matrix
// OR
// DestroyComplexMatrix(&a);
// DestroyComplexMatrix(&b);
}
说明:
(1) 因为在后面的特征值的求解函数中需要用到,单单QR分解本身数学意义不大,因此单独的测试就免了;
(2) 后面释放结构体指针(也就是矩阵内存)的函数,调用的是DestroyComplexMatrixArray()而非DestroyComplexMatrix(),前面第一章给了二者实现的代码,读者留意一下,我个人建议读者定义的矩阵不多,最好不要调用这个destroy函数,老老实实一个一个调用DestroyComplexMatrix()/DestroyDoubleMatrix()函数释放内存;
(3) 另外,注意矩阵的个数numDoubleMatrixArray必须在调用DestroyComplexMatrixArray()函数之前求出来赋值给变量然后传给形参,而不能在DestroyComplexMatrixArray()函数里面求解,且ComplexMatrixArray必须定义为数组而不能是指针,因为无论是指针还是数组,实参传给形参时候,系统都会在函数内部将其当作指针处理,sizeof(ptr) ≡ 4Byte(64bit系统下),而不是指向的那段内存的实际大小,但是sizeof(数组名)却是实际的那段内存的大小,这也是数组和指针的区别,C语言里面这一块也是一个比较容易出错的地方(笔者此处也差点搞忘了!!!)。
三、矩阵特征值
特征值计算比较复杂,关于特征值和特征向量的求解,线性代数理论知识大家可以参考此篇文章,关于特征值和特征向量C语言的实现,笔者是根据此位大佬修改的,这里我直接贴出笔者魔改后的代码:
/* eigen values */
void EigenValue(const Matrix* matrix, Matrix* eigenvalue)
{
const int NUM = 100; // Iteration Times
// Local Matrice
Matrix Q, R, temp;
// Initiate
InitComplexMatrix(&Q, matrix->row, matrix->column);
InitComplexMatrix(&R, matrix->row, matrix->column);
InitComplexMatrix(&temp, matrix->row, matrix->column);
if (IsNullComplexMatrix(matrix))
{
puts("ERROE: An invalid matrix!\n");
return;
}
// Copy matrix to temp in order not to change matrix
CopyMatrix(matrix, &temp);
// QR Decompose and
for (int k = 0; k < NUM; ++k)
{
QR(&temp, &Q, &R);
MatrixMulMatrix(&R, &Q, &temp);
}
// Abstract Eigen Values from the Diagonal Elements of temp = Q * R
for (int k = 0; k < temp.row; ++k)
{
eigenvalue->arrayComplex[k] = temp.arrayComplex[k * temp.column + k];
}
// Free Local Memory
Matrix ComplexMatrixArray[] = {
R, Q, temp };
int numComplexMatrixArray = sizeof(ComplexMatrixArray) / sizeof(Matrix);
DestroyComplexMatrixArray(ComplexMatrixArray, numComplexMatrixArray); // Complpex Matrix
// return eigenvalue;
}
说明:
(1) 因为在后面的特征向量的求解中需要用到,因此单独的测试就免了;
(2) 同样这里销毁结构体内部指针(即矩阵)内存调用的也是一次销毁函数DestroyComplexMatrixArray();
(3) 读者注意一下笔者在此特征值求解函数中调用的一些函数,如果忘记了可以跳到前面几章节review。
三、矩阵特征向量
我是一个没有感情的复读机:特征向量的计算更为复杂,关于特征向量的求解,线性代数理论知识大家可以参考此篇文章,关于特征向量C语言的实现,笔者是根据此位大佬修改的,这里我直接贴出笔者修魔改后的代码:
/* Negative Complex : Complex_B = -Complex_A = -creal(Complex_A) - cimag(Complpex_A) */
ComplexType NegativeComplex(const ComplexType Complex_A)
{
ComplexType Complex_B;
Complex_B = _Cmulcr(Complex_A, -1.0);
// OR
/*
Complex_B._Val[0] = -creal(Complex_A);
Complex_B._Val[1] = cimag(Complex_A) * (-1.0);
*/
return Complex_B;
}
/* eigen vectors */
void EigenVector(const Matrix* matrix, const Matrix* eigenvalue, Matrix* eigenvector)
{
if (IsNullComplexMatrix(matrix) || IsNullComplexMatrix(eigenvalue))
{
puts("ERROE: An invalid matrix!\n");
return;
}
int i, j, q;
int m;
// Access to Eigen Values
int count;
int num = MatrixRow(matrix); // = matrix->row: Numbers of Eigen Values or Cols
ComplexType evalue;
// Access to temp
ComplexType sum, midsum, mid;
Matrix temp; // temp = A - λI: (A - λI) * x = 0
InitComplexMatrix(&temp, matrix->row, matrix->column);
for (count = 0; count < num; ++count)
{
// Calculate x: Ax = λ * x
evalue = eigenvalue->arrayComplex[coun