矩阵妙法—从基础运算到进阶使用

引言

矩阵运算(非MATLAB下)

     ~~~~          ~~~    学过线性代数或数值计算等数学课程的人应该都能意识到，矩阵是一种非常重要的数据表示形式，它的由来主要是数学家在研究解方程组时所衍生出的一种求解工具。
     ~~~~          ~~~    而今天笔者想谈的是如何在C/C++中，使用二维数组来表示并模拟矩阵的一些运算情况。虽然矩阵包含一系列运算（如加减，求逆，转置、稀疏矩阵等），但本文还想谈一些非数值计算下的一些使用心得。

矩阵加法、乘法模拟

1.矩阵加法

加减原理概要

     ~~~~          ~~~    矩阵加法运算相对来说较为简单，除了注意它的一些数学前提以外（如运算的两矩阵的行和列应该是一致的），还应注意编程语言使用上的些细节（这里以C++为例）。
　　下面这个例子简单地演示了矩阵Ａ和矩阵Ｂ两者相加，其中row和col分别代表行数和列数，矩阵C保存相加的结果。

实例分析

     ~~~~          ~~~    其中，需要注意的是函数声明时采用的是一维数组的形式，这样可以很方便的采用线性遍历的方式(起点->终点)完成对应相加操作。而若采用二维数组作为形参，此时则需指明数组的行数与列数，使得该函数不具有通用性。

#include<iostream>

using namespace std;

int matrix_plus(double A[], double B[], int row, int col, double C[]);

int main()
{
	double A[3][3]={{1.0,2.0,3.0},
					{4.0,5.0,6.0},
					{7.0,8.0,9.0}};
	double B[3][3]={{2.0,-2.0,1.0},
					{1.0,3.0,9.0},
					{17.0,-3.0,7.0}};
	double C[3][3];
	
	double *pA=&A[0][0], *pB=&B[0][0], *pC=&C[0][0];
	
	int row = 3, col = 3;
	
	printf("矩阵A和矩阵B相加的结果为:\n");
	matrix_plus(pA, pB, row, col, pC);
	for(int i=0; i<row; i++){
		for(int j=0; j<col; j++){
			printf("%10.3f", C[i][j]);
		}
		cout<<endl;
	}
	
return 0;
}

int matrix_plus(double A[], double B[], int row, int col, double C[])
{
	for(int i=0; i<row*col; i++){ //一维线性遍历(因二维数组本质上也是一维数组,按行排列的)
		C[i] = A[i] + B[i];  //对应项相加
	}
	return 0;
}

     ~~~~          ~~~    看完这个例子，回想我们之前听过的一些话语：“对于算法来说，重要的是思想，而不是使用何种编程语言”。诚然，这种说法也不算错误，但对于具体的实现来说，语言使用上的细节也是值得注意的。

2.矩阵乘法

乘法原理概要

     ~~~~          ~~~    矩阵的乘法相对于加减来说会复杂一些，如该矩阵乘法$C_{m\times k}=A_{m\times n}{B}_{n\times k}$，其中乘积矩阵C中各个元素的值为：
$$c_{ij}=\sum _{t=0}^{n-1}{a}_{it}{b}_{tj}\text{，其中}\left(i=\text{0,1,}\cdot \cdot \cdot ,m-\text{1;}j=\text{0,1,}\cdot \cdot \cdot ,k-1\right)$$
　　简单说来，就是矩阵A的第i行的所有元素与矩阵B的第j列的所有元素对应相乘，并把所相乘的结果相加。最终得到的值就是矩阵C的第i行第j列的值。如果还有些不明白可以继续看下面这张运算图。$$\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\end{array}\right]\times \left|\begin{array}{cc}{b}_{11}& b_{12}\\ b_{21}& b_{22}\\ b_{31}& b_{32}\end{array}\right|=\left[\begin{array}{cc}{a}_{11}{b}_{11}+a_{12}{b}_{21}+a_{13}{b}_{31}& a_{11}{b}_{12}+a_{12}{b}_{22}+a_{13}{b}_{32}\\ a_{21}{b}_{11}+a_{22}{b}_{21}+a_{23}{b}_{31}& a_{21}{b}_{12}+a_{22}{b}_{22}+a_{23}{b}_{32}\end{array}\right]$$

实例分析

     ~~~~          ~~~    这里依然采用的是一维遍历方向，具体实现方法已集成于函数体matrix_mul中。整体的思想上是用两层大的for循环来控制矩阵C的计算(m*k阶)，故函数体中循环变量i变化范围为0 ~ m-1，循环变量j变化范围为0 ~ k-1。

#include<iostream>

using namespace std;

int matrix_mul(double A[], double B[], int m, int n, int k, double C[]);

int main()
{
	double A[3][3]={{1.0,2.0,3.0},
					{4.0,5.0,6.0},
					{7.0,8.0,9.0}};
	double B[3][3]={{2.0,-2.0,1.0},
					{1.0,3.0,9.0},
					{17.0,-3.0,7.0}};
	double C[3][3];
	
	double *pA=&A[0][0],*pB=&B[0][0],*pC=&C[0][0];
	
	int row = 3, col = 3, new_col = 3;   //new_col为乘积矩阵C所对应的列数
	
	printf("矩阵A和矩阵B乘的结果为:\n");
	matrix_mul(pA, pB, row, col, new_col, pC); 
	for(int i=0; i<row; i++){
		for(int j=0; j<new_col; j++){
			printf("%10.3f ", C[i][j]);
		}
		printf("\n");
	}

return 0;
}

int matrix_mul(double A[], double B[], int m, int n, int k, double C[])
{
	int u;
	
	/*计算m*k阶C矩阵的结果*/
	for (int i=0; i<m; i++){               //下标i: C矩阵共m行
		for (int j=0; j<k; j++){          //下标j: C矩阵共k列
			u = i * k + j;               //下标u: 采用的是一维遍历方式
			C[u] = 0.0;                 //先清零待计算的元素
			
			for(int l=0; l<n; l++){           //下标l: B矩阵共n行(或A矩阵共n列)
				C[u] += A[i*n+l]*B[l*k+j];   //相乘累加: 均采用一维遍历方式, A走行方向, B走列方向
			}
		}
	}
return 0;
}

【测试结果】

     ~~~~          ~~~    看完该例子后，可能有些读者会不大理解程序中出现的这三个表达式：
$$i\ast k+j(将i替换为m或许有助于你的理解,i取不同值对应矩阵C不同的行,靠j来选择不同的行中第几个元素)$$$$i\ast n+l(原理同上,i取不同值对应矩阵A不同的行,靠l来选择不同的行中第几个元素)$$$$l\ast k+j(略...)$$

矩阵旋转、螺旋矩阵

矩阵旋转概要

     ~~~~          ~~~    矩阵旋转在代数学中应称之为转置，但落实到二维数组中进行模拟时也还是以数组下标的变化规律作为问题的切入点。首先，咱们先通过一个4*4的矩阵来分析顺时针旋转以及逆时针旋转的规律。

　　图中已将两种旋转情况合并在一张图中进行了分析，旋转元素分别为图中的①~④。

实例分析

#include<iostream>

using namespace std;

int mat_rot(int (*mat)[5], int row, int col);

int main()
{
	int mat[3][5]={0};

	cout<<"请输入3*5矩阵中的元素, 待会给出旋转结构:\n";
	for(int i=0; i<3; i++){
		for(int j=0; j<5; j++){
			scanf("%d", &mat[i][j]);
		}
	}
	mat_rot(mat,3,5);
	
return 0;
}

int mat_rot(int (*mat)[5], int row, int col)
{
	int dir = 1;  //1:顺时针, 0:逆时针
	
	printf("请输入旋转方向. 1:顺时针 0:逆时针.\n");
	
	scanf("%d",&dir);
	if(dir==1){
		for(int i=0; i<row; i++){
			for(int j=0; j<col; j++){
				printf("%2d ", mat[row-j-1][i]);
			}
			cout<<endl;
		}
	}else{
		for(int i=0; i<5; i++){
			for(int j=0; j<3; j++){
				printf("%2d ", mat[j][col-i-1]);
			}
			cout<<endl;
		}
	}
	cout<<endl;
		
	return 0;
}

【测试结果】

螺旋旋转概要

     ~~~~          ~~~    当数据值在矩阵中按一定的行走规律相应变化时，螺旋矩阵也就相应显现。所以依据此种描述思想，我们也应该用编程来模拟这种数据值的行走过程。这里我想以两个较为典型的例子来说明，具体请观察下图。

     ~~~~          ~~~    从图中可观察得出，蛇形线(S型填数)的基本步法为：向右移动->向下移动->向左移动->向下移动，这四个动作(依次循环)。
     ~~~~          ~~~    而螺旋线的基本步法为：向下->向左->向上->向右，这四个动作。只不过需要处理每个方向上的行走范围问题。
     ~~~~          ~~~    下面给出一份示例参考程序，读者可根据此图好好体会程序中是如何分解处理不同的动作。

实例分析

#include<iostream>

using namespace std;

#define MAX_S 4
#define MAX_H 20

int mat_snake(int (*mat)[MAX_S]);      //蛇形线
int mat_helix(int (*mat)[MAX_H]);     //螺旋线

int main()
{
	int mat_s[MAX_S][MAX_S]={0};
	int mat_h[MAX_H][MAX_H]={0};
	
	mat_helix(mat_h);	
	cout<<endl;
	mat_snake(mat_s);
			
	return 0;
}

/*宏MAX_S可改变矩阵大小, 变量n可改变起点数值*/
int mat_snake(int (*mat)[MAX_S])
{
	int n,x=0,y=0,tot=0;   //x,y用于指示坐标, tot指示变化的数值
	
	cout<<"snake_begin!\n";
	cout<<"请输入蛇形坐标起点数值："; 
	cin>>n;
	tot = mat[0][0] = n;  //蛇形走位的起点数值
	
	//蛇形基本步法
	while(tot < n + MAX_S*MAX_S){  
		
		/*向右移动,++y*/
		while(y<MAX_S-1)  //矩阵右边界
			mat[x][++y]=++tot;
		
		/*向下,++x*/
		mat[++x][y]=++tot;
		
		/*向左移动,--y*/
		while(y>0)  //矩阵左边界
			mat[x][--y]=++tot;
		
		/*向下,++x*/
		mat[++x][y]=++tot;
	}
	
	for(int i=0; i<MAX_S; i++){
		for(int j=0; j<MAX_S; j++){
			printf("%2d ", mat[i][j]);
		}
		cout<<endl;
	}
	cout<<"snake_end!\n";
	
return 0;
}

/*起点默认为右上角(数值为1), 宏MAX_H用于控制矩阵的最大范围, 变量n控制当前矩阵大小*/
int mat_helix(int (*mat)[MAX_H])
{
	int n,x,y,tot=0;   //x,y用于指示坐标, tot指示变化的数值
	
	cout<<"helix_begin!\n";
	cout<<"请输入螺旋线方阵大小："; 
	cin>>n;
	tot = mat[x=0][y=n-1] = 1;  //螺旋线走位的起点(右上角)
	
	//螺旋线基本步法
	while(tot < n*n){   
		
		/*向下：不超出给定矩阵的下边界, 以及防范填充非空的格子*/
		while(x<n-1 && !mat[x+1][y]){
			mat[++x][y]=++tot;
		}
				
		/*向左：*/
		while(y>0 && !mat[x][y-1]){
			mat[x][--y]=++tot;
		}
				
		/*向上：*/
		while(x>0 && !mat[x-1][y]){
			mat[--x][y]=++tot;
		}
				
		/*向右：*/
		while(y<n-1 && !mat[x][y+1]){
			mat[x][++y]=++tot;
		}			
	}
	
	for(int i=0; i<n; i++){
		for(int j=0; j<n; j++){
			printf("%2d ", mat[i][j]);
		}
		cout<<endl;
	}
	cout<<"helix_begin!\n";
	
return 0;
}

【测试结果】

参阅书籍
C/C++常用算法手册
狄泰软件C进阶剖析教程
同济版-线性代数