一、前言
本章开始,笔者将详细讲解实数和复数矩阵的模(范数)、协方差等运算,并给出函数和测试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))
return;