C语言实现线性代数中基的标准正交化

    线性代数中一个很重要的概念就是基，一组基是构成线性空间“最高效的方式“，即若多余这个向量个数会出现冗余，若少于这个向量个数会出现无法构成该向量空间。由一般向量构成的一组基往往不是标准正交的，所以需要进行标准正交化！

    线性代数“基“的概念在很多领域都有应用，包括机器人学，图像处理，数据降维等，具体可以参考相关线性代数教材。

    本文采用C语言实现了基的标准正交化，正交化方法为格兰姆-施密特正交化（参见相关线性代数教材或维基百科）。程序功能为：询问用户输入向量的维数，用户输入一组基，构成n*n的矩阵（为了符合输入习惯，采用行向量的方式进行输入），然后对矩阵进行转置，将行向量变为列向量（符合线性代数的习惯），在进行判断之前，先判断用户输入的向量能否构成一组基，判断方式在isInverse（）这个函数中，采用的是消元法判断矩阵是否可逆，若可逆则进行正交化，若不可逆，输出该矩阵不可逆，退出程序，具体程序如下：

/*================================================== # Author:      Joker@HIT    [e-mail] # Filetype:    C source code # Environment: Linux & Ubuntu 14.04 # Tool:        Vim & Gcc  # Date:        Sat Aug 20 2016 # Descprition: Input a couple of basis as rows of  a matrix. the program tansposes it,and examine it's inverse whether exists or not! If exists, orthogonizeit standardly!==================================================*/#include<stdio.h>#include<stdlib.h>#include<stdbool.h>#include<math.h>void transpose(int, double [*][*]);void show_matrix(int, double [*][*]);void orthogonize(int, double [*][*]);bool isInverse(int, double [*][*]);void exchange(int, int, double [*][*]);int main(void){int dimention, i, j;bool judge;printf("Enter the dimention of the basis: ");scanf("%d", &dimention);double basis_matrix[dimention][dimention];printf("Input the basis (as rows of the matrix):\n");for(i=0;i<dimention;i++)for(j=0;j<dimention;j++)scanf("%lf", *(basis_matrix+i)+j);transpose(dimention, basis_matrix);printf("The matrix you input as follow:\n");show_matrix(dimention, basis_matrix);judge = isInverse(dimention, basis_matrix);if(!judge){printf("The column vector of matrix can't consists of the basis!\n");exit(1);}printf("The standard orthogonal basis matrix as follow:\n");orthogonize(dimention, basis_matrix);show_matrix(dimention, basis_matrix);return 0;}void transpose(int dim, double matrix[dim][dim]){int i, j;double replica[dim][dim];for(i=0;i<dim;i++)for(j=0;j<dim;j++)*(*(replica+i)+j) = *(*(matrix+i)+j);for(i=0;i<dim;i++)for(j=0;j<dim;j++)*(*(matrix+i)+j) = *(*(replica+j)+i);return;}void show_matrix(int dim, double matrix[dim][dim]){int i, j;for(i=0;i<dim;i++){for(j=0;j<dim;j++)printf("%-9.3f", *(*(matrix+i)+j));putchar('\n');}return;}void orthogonize(int dim, double matrix[dim][dim]){int i, j, k;double numerator, demoninator;double proj_vector[dim];for(i=1;i<dim;i++){for(k=0;k<dim;k++)proj_vector[k] = 0;for(j=0;j<i;j++){numerator = demoninator = 0;for(k=0;k<dim;k++){numerator += matrix[k][i] * matrix[k][j];demoninator += matrix[k][j] * matrix[k][j];}for(k=0;k<dim;k++)proj_vector[k] += numerator / demoninator * matrix[k][j];}for(k=0;k<dim;k++)matrix[k][i] = matrix[k][i] - proj_vector[k];}for(i=0;i<dim;i++){demoninator = 0;for(j=0;j<dim;j++)demoninator += matrix[j][i] * matrix[j][i];for(j=0;j<dim;j++)matrix[j][i] = matrix[j][i] / sqrt(demoninator);}return;}bool isInverse(int dim, double matrix[dim][dim]){double replica[dim][dim];double coef;int i, j, k;bool result = true;for(i=0;i<dim;i++)for(j=0;j<dim;j++)replica[i][j] = matrix[i][j];for(i=0;i<dim-1;i++){exchange(dim, i, replica);if(fabs(replica[i][i])<0.00001){result = false;break;}for(j=i+1;j<dim;j++){coef = replica[j][i] / replica[i][i];for(k=0;k<dim;k++)replica[j][k] = replica[j][k] - coef * replica[i][k];}}if(fabs(replica[dim-1][dim-1])<0.00001)result = false;return result;}void exchange(int dim, int flag, double matrix[dim][dim]){double temp[dim];int i, j;for(i=flag+1;i<dim;i++)if(fabs(matrix[flag][flag])<fabs(matrix[i][flag])){for(j=0;j<dim;j++){temp[j] = matrix[flag][j];matrix[flag][j] = matrix[i][j];matrix[i][j] = temp[j];}}return;}

运行结果如下图所示：

其实上述代码在exchange（）函数中可以进行一定的改进，采用指针数组的方式，然后排序操作是对指针进行排序，而不是对数组进行排序，优化程序运行时间和空间！

若是读者对代码有什么意见，欢迎留言或者邮箱探讨