C语言实现实数和复数矩阵及其各种运算（四）

一、前言

本章开始，笔者将详细讲解实数和复数矩阵的模（范数）、协方差等运算，并给出函数和测试demo的C代码，以及与matlab计算的结果；

并且，本章也出现了我在编写函数、后期与matlab联调跑数据时候踩到的一个坑，读者稍微留意一下。

另，这将是我此连载文章的系列精品之一，矩阵运算的功能越来越复杂，但是笔者全部分模块进行，化繁为简，化整为零，所以需要调用的小函数越来越多，读者遇到不熟悉的函数可以随时Jump到（一）、（二）、（三）的内容。

二、矩阵的模（范数）

这里笔者提到的矩阵的模（范数），数学上严格来讲，指的是矩阵的F范数，即矩阵的各个元素平方之和再开平方根，它通常也叫做矩阵的L2范数。关于向量和矩阵的其他的模（范数）的数学定义和求解，读者可以check此链接。笔者单独引入L2范数概念和求解函数，是因为后面协方差矩阵求解过程有一个步骤需要求解矩阵L2范数，所以我直接单独写一个函数实现。
由于实数矩阵和复数矩阵的L2范数求解逻辑完全一样，这里我只给出复数矩阵L2范数的函数code，话不多说，直接把C code贴出来：

/* 2-norm of a Matrix */
ComplexType MatrixNorm2(const Matrix* matrix)
{// OR//    int size = matrix->row * matrix->column;ComplexType norm;norm._Val[0] = 0; norm._Val[1] = 0;if (IsNullComplexMatrix(matrix)){puts("ERROE: An invalid matrix!\n");return (csqrt(norm));}else{for (int i = 0; i < matrix->row; ++i){for(int j = 0; j < matrix->column; ++j)norm = AddComplex(norm, _Cmulcc(matrix->arrayComplex[i * matrix->column + j], matrix->arrayComplex[i * matrix->column + j]));}// OR
//      for (row_i = 0; row_i < MatrixSize(matrix); ++row_i)
//      {//          for (column_j = 0; column_j < MatrixSize(matrix); ++column_j)
//              norm = *(AddComplex(norm, _Cmulcc(matrix->arrayComplex[row_i * matrix->row + column_j], matrix->arrayComplex[row_i * matrix->row + column_j])));
//      }return (csqrt(norm));   // csqrt w.r.t a Complex is different frm sqrt w.r.t a Double}
}

说明：
(1) matlab里面求解向量和矩阵的各种模（范数）都是直接调用norm()函数然后修改参数即可。这里我贴上上面链接里给的matlab code以作补充：

clear all;clc;
%% 求向量的范数
X = [2, 3, -5, -7];    % 初始化向量X
XLfs1 = norm(X, 1);    % 向量的1-范数
XLfs2 = norm(X, 2);    % 向量的2-范数
XLfsz = norm(X, inf);  % 向量的正无穷范数
XLfsf = norm(X, -inf); % 向量的负无穷范数%% 求矩阵的范数
A=[2, 3, -5, -7;4, 6,  8, -4;6, -11, -3, 16];     %初始化矩阵AJZfs1 = norm(A, 1);        % 矩阵的1-范数
JZfs2 = norm(A, 2);        % 矩阵的2-范数
JZfswq = norm(A, inf);     % 矩阵的无穷范数
JZhfs = sum(svd(A));       % 矩阵的核范数
JZL1fs = sum(sum(abs(A))); % 矩阵的L1范数
JZFfs=norm(A, 'fro');      % 矩阵的F范数
JZL21fs=norm(A(:, 1), 2) + norm(A(:, 2), 2) + norm(A(:, 3), 2)++ norm(A(:, 4), 2); % 矩阵的L21范数

(2) L2范数求解结果是一个简单的实数/复数，因此这个函数直接用return返回一个值，就不用指针多此一举了；
(3) 读者可以自行测试L2范数函数，由于后续协方差矩阵求解函数调用了此函数，所以笔者不单独写demo测试了。

三、矩阵的协方差

(1) 由于实数矩阵和复数矩阵的协方差矩阵求解会有区别，因此此章笔者会着重讲解一下；
(2) 关于矩阵协方差矩阵的数学定义和求解公式，读者可以check此链接，笔者不做数学推导的工作；
(3) 我在项目中用到的是按列排列求的协方差矩阵，如果原矩阵是M X N，那么按列求协方差矩阵应该是N X N，读者一定要注意矩阵的size关系，尤其不是方阵时候，笔者在下面的代码中设定的每个矩阵的size都是严格按照数学原理的，所以对于方阵/非方阵都是通用的；
(3) 思路：先写一个协方差函数求解协方差矩阵的每一个协方差cell，然后再写一个函数调用协方差函数依次填充协方差矩阵。

1. 实数矩阵的协方差矩阵

直接贴出C code：

/* Covariance Matrix Cell */
DoubleType covDoubleMatrixCell(const DoubleType x[], const DoubleType y[], const int size )  // size refers to row number of matrix
{DoubleType cov; // = NULL;
//  cov = (DoubleType *)malloc(sizeof(DoubleType ));int i;DoubleType averX =0.0, averY = 0.0, sumX = 0.0, sumY = 0.0, sumXY = 0.0, subx =0.0, suby = 0.0;// Sum Values of Row/Columnfor (i = 0; i < size; i++){sumX += x[i];   // double + doublesumY += y[i];   // double + double}// Average Values of Row/ColumnaverX = sumX / size;averY = sumY / size;// Covariance Valuesfor (i = 0; i < size; i++){subx -= averX;suby -= averY;sumXY = sumXY + subx * suby;}cov = sumXY / (size - 1.0);return cov;
}/* Covarianxe Matrix of Complex Matrix */
void CovarianceDoubleMatrix(const Matrix2Double* matrix, const Matrix2Double* TransMat, Matrix2Double * CovMat)
{if (IsNullDoubleMatrix(matrix))return;else if (CovMat->row != DoubleMatrixColumn(matrix) || CovMat->column != DoubleMatrixColumn(matrix))return;else{for (int i = 0; i < CovMat->row; i++){for (int j = 0; j < CovMat->column; ++j){CovMat->arrayDouble[CovMat->column * i + j] = \covMatrixCell(TransMat->arrayDouble + i * TransMat->column, \TransMat->arrayDouble + j * TransMat->column, TransMat->column);}   }}
}

说明：
(1) 第一个函数covDoubleMatrixCell()是求取协方差矩阵的每个元胞，第二个函数CovarianceDoubleMatrix()是将元胞填入对应的协方差矩阵的位置；
(2) 关于第一个函数，由于是求取两个列向量之间的协方差，所以形参是两个数组（或者指针变量也可以）。由于指针依次遍历不太方便，因此我都是将矩阵转置，这样变成行向量，使得指针遍历访问角标直接依次+1即可。读者也可以不作转置变换，这样的话每次+i*column也可以，但是笔者建议最好不要这样，因为会很容易出错（亲身经历）。因为是实数矩阵，前面讲过，与复数转置之间的区别，这里只是位置对称变换，当然我给的函数实数复数完全通用，注意变量类型即可；
(3) 关于第二个函数，此函数只是填充协方差矩阵操作。注意角标即可。
(4) 本来准备写一个demo测试的，但是实数说到底也是复数，我写的所有的复数域的函数其实实数域完全也可用，我单独拿出来讲实数，只不过想给读者把二者一些原理性的差距单列出来。

2. 复数矩阵的协方差矩阵

直接贴出C code：

/* Covariance Matrix Cell: void */
ComplexType covMatrixCell(const ComplexType x[], const ComplexType y[], const int size )  // size refers to row number of matrix
{ComplexType cov; // = NULL;
//  cov = (ComplexType*)malloc(sizeof(ComplexType));int i;ComplexType averX, averY, sumX, sumY, sumXY, subx, suby;// Initiate Complex VariablesInitComplex(&averX);InitComplex(&averY);InitComplex(&sumX);InitComplex(&sumY);InitComplex(&sumXY);InitComplex(&subx);InitComplex(&suby);// Sum Values of Row/Columnfor (i = 0; i < size; i++){// Solution_1sumX = AddComplex(x[i], sumX);   // Complex + ComplexsumY = AddComplex(y[i], sumY);   // Complex + Complex// Solution_2
//      sumX = AddComplex(x[i * size], sumX);   // Complex + Complex
//      sumY = AddComplex(y[i * size], sumY);   // Complex + Complex// Solution_3
//      sumX._Val[0] = creal(x[i]) + sumX._Val[0];
//      sumX._Val[1] = cimag(x[i]) + sumX._Val[1];
//      sumY._Val[0] = creal(y[i]) + sumY._Val[0];
//      sumY._Val[1] = cimag(y[i]) + sumY._Val[1];}// Average Values of Row/ColumnaverX = _Cmulcr( sumX, pow(size * 1.0, -1.0) );   // Tip: 1/size = pow(size, -1)averY = _Cmulcr( sumY, pow(size * 1.0, -1.0) );// Covariance Valuesfor (i = 0; i < size; i++){// Solution_1subx = SubComplex(x[i], averX);suby = SubComplex(y[i], averY);// Solution_2
//      subx = SubComplex(x[i * size], averX);
//      suby = SubComplex(y[i * size], averY);// Solution_3
//      subx._Val[0] = creal(x[i]) - creal(averX);
//      subx._Val[1] = cimag(x[i]) - cimag(averX);
//      suby._Val[0] = creal(y[i]) - creal(averY);
//      suby._Val[1] = cimag(y[i]) - cimag(averY);sumXY = AddComplex(sumXY, _Cmulcc(subx, suby));
//      sumXY._Val[0] = sumXY._Val[0] + _Cmulcc(subx, suby)._Val[0];
//      sumXY._Val[1] = sumXY._Val[1] + _Cmulcc(subx, suby)._Val[1];}//  *cov = _Cmulcr(sumXY, pow(size - 1, -1));  // OR:cov = _Cmulcr( sumXY, pow(size - 1.0, -1.0) );return cov;
}/* Covarianxe Matrix of Complex Matrix */
void CovarianceComplexMatrix(const Matrix* matrix, const Matrix* TransMat, const Matrix *TransPoseMat, Matrix* CovMat)
{if (IsNullComplexMatrix(matrix))return;else if (CovMat->row != MatrixColumn(matrix) || CovMat->column != MatrixColumn(matrix))return;else{for (int i = 0; i < CovMat->row; i++){for (int j = 0; j < CovMat->column; ++j){CovMat->arrayComplex[CovMat->column * i + j] = \covMatrixCell(TransMat->arrayComplex + i * TransMat->column, \TransPoseMat->arrayComplex + j * TransPoseMat->column, TransMat->column);}   }}
}

说明：
(1) 注释也是我写算法时候，没有自己单独定义复数加减运算的函数后来加上的，读者也可以直接忽略；
(2) covMatrixCell()函数与实数矩阵在编程逻辑上无明显区别，只是需要自己定义复数的运算函数；
(3) CovarianceComplexMatri()函数形参依次是：原矩阵、转置矩阵、转位置矩阵、协方差矩阵，因此这一点，要比实数域矩阵多一个转位置矩阵，事实上，狭义上来讲，实数域的转置矩阵与复数域的转位置矩阵的数学逻辑，即计算结果是一样的，广义上，实数域的转置矩阵的本质还是复数域的转置矩阵。笔者书写的复数转置矩阵函数TransMatrix()，复数域/实数域转置矩阵通用；
(4) 复数和实数矩阵协方差矩阵求解的区别在第二个函数，即协方差矩阵函数，对于实数矩阵来讲，转置和转位置结果是一样的，因此求协方差时候调用任意一种转置矩阵均可。其区别部分我将其代码单列出来：

// Complex Matrix
CovMat->arrayComplex[CovMat->column * i + j] = \covMatrixCell(TransMat->arrayComplex + i * TransMat->column, \TransPoseMat->arrayComplex + j * TransPoseMat->column, TransMat->column);
// Double Matrix:
// Solution_1
CovMat->arrayDouble[CovMat->column * i + j] = \covDoubleMatrixCell(TransMat->arrayDouble + i * TransMat->column, \TransMat->arrayDouble + j * TransMat->column, TransMat->column);
// Solution_2
CovMat->arrayDouble[CovMat->column * i + j] = \covDoubleMatrixCell(TransPoseMat->arrayDouble + i * TransPoseMat->column, \TransPoseMat->arrayDouble + j * TransPoseMat->column, TransMat->column);

(5) 测试demo比较多，调用了很多前面章节用到的函数：

#include<stdio.h>
#include<stdlib.h>
#include<complex.h>
#include <math.h>  // pow()typedef _Dcomplex ComplexType;
typedef double DoubleType;typedef struct
{int row, column;ComplexType* arrayComplex;
}Matrix;typedef struct
{int row, column;DoubleType* arrayDouble;
}Matrix2Double;typedef enum
{False = 0, True = 1
}Bool;/*Initiation of Complex Number*/
void InitComplex(ComplexType* Complex)
{Complex->_Val[0] = 0.0;Complex->_Val[1] = 0.0;
}/*Validity of Complex Matrix
*/
Bool IsNullComplexMatrix(const Matrix* matrix)
{int size = matrix->row * matrix->column;if (size <= 0 || matrix->arrayComplex == NULL){return True;}return False;
}/*Initiation of Complex Matrix*/
void InitComplexMatrix(Matrix* matrix, int row, int column)
{int size = row * column * sizeof(ComplexType);if (size <= 0){puts("ERROE: An invalid matrix!\n");return;}matrix->arrayComplex = (ComplexType*)malloc(size);if (matrix->arrayComplex){matrix->row = row;matrix->column = column;for (int row_i = 0; row_i < row; ++row_i){for (int column_j = 0; column_j < column; ++column_j){InitComplex(matrix->arrayComplex + row_i * matrix->column + column_j);}}}
}/*Initiation of Double Matrix
*/
void InitDoubleMatrix(Matrix2Double* matrix, int row, int column)
{int size = row * column * sizeof(DoubleType);if (size <= 0){puts("ERROE: An invalid matrix!\n");return;}matrix->arrayDouble = (DoubleType*)malloc(size); // initiate pointer/rowif (matrix->arrayDouble){matrix->row = row;matrix->column = column;for (int row_i = 0; row_i < row; ++row_i){for (int column_j = 0; column_j < column; ++column_j){matrix->arrayDouble[row_i * matrix->column + column_j] = 0.0;}}}
}/*Free Memory of Complex Matrix
*/
void DestroyComplexMatrix(Matrix* matrix)
{if (!IsNullComplexMatrix(matrix)){free(matrix->arrayComplex);matrix->arrayComplex = NULL;}matrix->row = matrix->column = 0;
}/*Free Memory of Complex Matrice Array
*/
void DestroyComplexMatrixArray(Matrix matrixArray[], int num)  // Array Transfer--->Pointer Transfer
{if (num)      // if no cell{for (int i = 0; i < num; i++){DestroyComplexMatrix(&matrixArray[i]);  // Nested Call of DestroyComplexMatrix()}}
}/*return matrix row size
*/
int MatrixRow(const Matrix* matrix)
{return matrix->row;
}/*return matrix column size
*/
int MatrixColumn(const Matrix* matrix)
{return matrix->column;
}/*Add Complex: Complex_C = Complex_A + Complex_B
*/
ComplexType AddComplex(const ComplexType Complex_A, const ComplexType Complex_B)
{ComplexType Complex_C;Complex_C._Val[0] = creal(Complex_A) + creal(Complex_B);Complex_C._Val[1] = cimag(Complex_A) + cimag(Complex_B);return Complex_C;
}/*Subvision Complex: Complex_C = Complex_A - Complex_B
*/
ComplexType SubComplex(const ComplexType Complex_A, const ComplexType Complex_B)
{ComplexType Complex_C;Complex_C._Val[0] = creal(Complex_A) - creal(Complex_B);Complex_C._Val[1] = cimag(Complex_A) - cimag(Complex_B);return Complex_C;
}/*Conjugate/Transposal Complex: Complex_B = creal(Complex_A) - cimag(Complpex_A)
*/
ComplexType ConjugateComplex(const ComplexType Complex_A)
{ComplexType Complex_B;Complex_B._Val[0] = creal(Complex_A);Complex_B._Val[1] = cimag(Complex_A) * (-1.0);return Complex_B;
}/*Transposition: matrixB = matrixA': Position and Conjugate Complex of Matrix Cell are all changed: Coressponding to transpose of double matrix !!!
*/
void TransMatrix(const Matrix* matrixA, Matrix* matrixB)
{if (IsNullComplexMatrix(matrixA)){puts("ERROE: An invalid matrix!\n");return;}else if ((matrixA->row != matrixB->column) || (matrixA->column != matrixB->row)){puts("ERROE: An incompatable matrix!\n");return;}else{for (int row_i = 0; row_i < matrixA->row; row_i++){for (int column_j = 0; column_j < matrixA->column; column_j++){matrixB->arrayComplex[column_j * matrixB->column + row_i] = \ConjugateComplex(matrixA->arrayComplex[row_i * matrixA->column + column_j]);}}}
}/*Change Postion: [i, j] -- > [j, i]: Only Postion of Matrix Cell is changed: In order to get access a column vector pointer easily !!!
*/
void TransPoseMatrix(const Matrix* matrixA, Matrix* matrixB)
{                                                   if (IsNullComplexMatrix(matrixA)){puts("ERROE: An invalid matrix!\n");return;}else if ((matrixA->row != matrixB->column) || (matrixA->column != matrixB->row)){puts("ERROE: An incompatable matrix!\n");return;}else{for (int row_i = 0; row_i < matrixA->row; row_i++){for (int column_j = 0; column_j < matrixA->column; column_j++){matrixB->arrayComplex[column_j * matrixB->column + row_i] = \matrixA->arrayComplex[row_i * matrixA->column + column_j]; }}}
}/*Covariance Matrix Cell: void
*/
ComplexType covMatrixCell(const ComplexType x[], const ComplexType y[], const int size)
{ComplexType cov;int i;ComplexType averX, averY, sumX, sumY, sumXY, subx, suby;InitComplex(&averX);InitComplex(&averY);InitComplex(&sumX);InitComplex(&sumY);InitComplex(&sumXY);InitComplex(&subx);InitComplex(&suby);for (i = 0; i < size; i++){sumX = AddComplex(x[i], sumX);sumY = AddComplex(y[i], sumY);averX = _Cmulcr(sumX, pow(size * 1.0, -1.0));averY = _Cmulcr(sumY, pow(size * 1.0, -1.0));}for (i = 0; i < size; i++){subx = SubComplex(x[i], averX);suby = SubComplex(y[i], averY);sumXY = AddComplex(sumXY, _Cmulcc(subx, suby));}cov = _Cmulcr(sumXY, pow(size - 1.0, -1.0));return cov;
}/* Covarianxe Matrix of Complex Matrix
*/
void CovarianceComplexMatrix(const Matrix* matrix, const Matrix* TransMat, const Matrix *TransPoseMat, Matrix* CovMat)
{if (IsNullComplexMatrix(matrix))return;else if (CovMat->row != MatrixColumn(matrix) || CovMat->column != MatrixColumn(matrix))return;else{for (int i = 0; i < CovMat->row; i++){for (int j = 0; j < CovMat->column; ++j){CovMat->arrayComplex[CovMat->column * i + j] = \covMatrixCell(TransMat->arrayComplex + i * TransMat->column, \TransPoseMat->arrayComplex + j * TransPoseMat->column, TransMat->column);}   }}
}
int main(void)
{Matrix C;InitComplexMatrix(&C, 3, 3);  // 3 * 3    C.arrayComplex[0]._Val[0] = 2;  C.arrayComplex[0]._Val[1] = 2;C.arrayComplex[1]._Val[0] = 3;  C.arrayComplex[1]._Val[1] = 7;C.arrayComplex[2]._Val[0] = 3;  C.arrayComplex[2]._Val[1] = -6;C.arrayComplex[3]._Val[0] = 6;  C.arrayComplex[3]._Val[1] = 2;C.arrayComplex[4]._Val[0] = -7; C.arrayComplex[4]._Val[1] = 7;C.arrayComplex[5]._Val[0] = 0;  C.arrayComplex[5]._Val[1] = 6;C.arrayComplex[6]._Val[0] = 0;  C.arrayComplex[6]._Val[1] = 6;C.arrayComplex[7]._Val[0] = 5;  C.arrayComplex[7]._Val[1] = -4;C.arrayComplex[8]._Val[0] = -2;  C.arrayComplex[8]._Val[1] = -1;Matrix C_Tran;InitComplexMatrix(&C_Tran, MatrixColumn(&C), MatrixRow(&C));  // 3 * 3TransMatrix(&C, &C_Tran);Matrix C_TransPose;InitComplexMatrix(&C_TransPose, MatrixColumn(&C), MatrixRow(&C));  // 3 * 3TransPoseMatrix(&C, &C_TransPose);Matrix C_Cov;InitComplexMatrix(&C_Cov, MatrixColumn(&C), MatrixColumn(&C));  // 3 * 3CovarianceComplexMatrix(&C, &C_Tran, &C_TransPose, &C_Cov);for (int i = 0; i < C_Cov.row; i++)for (int j = 0; j < C_Cov.column; j++)printf("%lf, %lf\n", creal(C_Cov.arrayComplex[i * C_Cov.column + j]), cimag(C_Cov.arrayComplex[i * C_Cov.column + j]));// Free Array of Memories of Complex Matrice by Calling ComplexMatrixArray() Matrix ComplexMatrixArray[] = { C, C_Tran, C_TransPose, C_Cov };int numComplexMatrixArray = sizeof(ComplexMatrixArray) / sizeof(Matrix);DestroyComplexMatrixArray(ComplexMatrixArray, numComplexMatrixArray);return 0;
}

测试结果我直接通过截图的形式给出，VS C（上）和matlab（下）：

PS：
(1) matlab调用cov()函数，第二个参数默认和为0时候都是以列排列求协方差矩阵，而为1时候表示以行排列求协方差矩阵，二者结果完全不同，我这里数据用的是3 X 3的方阵，如果不是方阵，协方差矩阵的维度都会不同。因此读者请一定要注意我写协方差矩阵时候矩阵大小的书写，请严格按照数学原理模板来设置，否则会出现笔者编程时候出现的致命编译错误（内存相关的exception）和算法解算错误，这也是本人写程序后期调试时候遇到的大坑之一！！！
(2) cov(C, 0)和cov(C, 1)的结果我也贴一下matlab的截图：

四、总结

本章涉及到的矩阵运算本身数学公式就比较多，建议读者先把对应的数学原理搞清楚再对应我的函数看，才能更好消化。此外，对比较复杂的算法和功能，我们可以将其拆分为多个步骤，化繁为简，编写几个小函数，不仅有利于理清编程思路，也方便后续compile出错时分步debug，快速定位到出错的环节–小function。在C里面我们提倡结构化编程，在C++里，我们提倡模板化编程，在本章中，我充分融入了以上两种编程思想，笔者希望与读者一起体会和学习！

下一章，我将继续详细讲解复数矩阵的另外几个最为难啃的运算，包括矩阵QR分解、特征值、特征向量，敬请期待…

传送门：（五）…