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];     %初始化矩阵A
 
JZfs1 = 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/Column
	for (i = 0; i < size; i++)
	{
   
		sumX += x[i];   // double + double
		sumY += y[i];   // double + double
	}
	// Average Values of Row/Column
	averX = sumX / size;
	averY = sumY / size;

	// Covariance Values
	for (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))