动手实现矩阵广义逆求解（2022.5.4）

引言
1、逆矩阵和广义逆矩阵简介
- 1.1 逆矩阵
- 1.2 广义逆矩阵
- 1.3 广义逆矩阵求解流程
2、广义逆矩阵求解代码
- 2.1 Matlab代码及运行结果
- - 2.1.1 使用标准函数 pinv
- 2.2 C#代码及运行结果
- - 2.2.1 使用MathNet库
  - 2.2.2 自己写
- 2.3 C++代码及运行结果
- - 2.3.1 使用Eigen库
  - 2.3.2 自己写
- 2.4 Java代码及运行结果
- - 2.4.1 使用ujmp库
  - 2.4.2 自己写
- 2.5 Python代码及运行结果
- - 2.5.1 使用numpy库
  - 2.5.2 自己写
- 2.6 R代码及运行结果
3、总结

引言

在矩阵论和线性代数等理论当中，作为基础部分的逆矩阵和广义逆矩阵都是十分重要的研究分支，广义逆矩阵起源于线性方程组的求解问题，实数域R和复数域C不断扩充矩阵元素范围和线性空间，计算方式的进步不断推动着广义逆的研究进展。这里利用初等行变换法和满秩分解法来分别求解方阵A的逆矩阵与广义逆矩阵，愿为数学爱好者和工程师提供参考和帮助。

1、逆矩阵和广义逆矩阵简介

1.1 逆矩阵

逆矩阵定义：对于一个n阶方阵A而言，若存在矩阵B，使得AB=BA=E（E为单位阵）(式1-1)AB=BA=E（E为单位阵） \tag{式1-1}AB=BA=E（E为单位阵）(式1-1)A=B−1(式1-2)A=B^{-1} \tag{式1-2}A=B−1(式1-2)B=A−1(式1-3)B=A^{-1} \tag{式1-3}B=A−1(式1-3) 那么将B称为方阵A的逆矩阵，B也用A−1A^{-1}A−1表示，这与倒数的定义有点类似，倒数的定义为乘积为1的两个数互为倒数。n阶方阵A可逆的充分必要条件是A满秩，即rank(A)=n(式1-4)rank(A)=n\tag{式1-4}rank(A)=n(式1-4) 也称A为满秩矩阵或非奇异矩阵，因此可逆矩阵也称为非奇异矩阵。满秩矩阵都能通过有限次初等行变换 / 列变换化为单位矩阵，因此满秩矩阵A的逆矩阵可表示成有限个初等矩阵的乘积。
逆矩阵具有如下七个性质：

方阵A与B具有相同的地位，A是B的逆矩阵，B也是A的逆矩阵，二者互为逆矩阵；
零矩阵不可逆，即不存在方阵B，使得OB=BO=EOB=BO=EOB=BO=E；
单位方阵E可逆，且逆矩阵为E^-1，E=E−1E=E^{-1}E=E−1；
若方阵A存在可逆矩阵，则逆矩阵必唯一；
若方阵A可逆，则其逆矩阵A−1A^{-1}A−1也可逆，有(A−1)−1=A{(A^{-1})}^{-1}=A(A−1)−1=A；
若方阵A可逆，则其转置矩阵ATA^{T}AT也可逆，有(AT)−1=(A−1)T{(A^{T})}^{-1}={(A^{-1})}^{T}(AT)−1=(A−1)T；
若A、B均为n阶可逆方阵，则相乘矩阵AB也可逆，有(AB)−1=B−1A−1{(AB)}^{-1}={B}^{-1}{A}^{-1}(AB)−1=B−1A−1。

利用初等行变换或初等列变换都可直接计算逆矩阵，初等行变换就是指下列3种变换：（1）以一个非零的数乘矩阵的某一行；（2）把矩阵的某一行的c倍加到另一行，这里c是实数；（3）互换矩阵中两行的位置。任意一个矩阵经过一系列初等行变换总能变成阶梯型矩阵。而矩阵秩的求解可通过将其变换为阶梯型矩阵后来统计得到。
(AE)→basicrowtransform(EA−1)\begin{pmatrix} &A &E \\ \end{pmatrix}\xrightarrow[basic row transform]{}\begin{pmatrix} &E &A^{-1} \\ \end{pmatrix}(AE)basicrowtransform(EA−1)
或者
(AE)→basiccoltransform(EA−1)\begin{pmatrix} &A \\ &E \end{pmatrix}\xrightarrow[basic col transform]{}\begin{pmatrix} &E \\ &A^{-1} \end{pmatrix}(AE)basiccoltransform(EA−1)

1.2 广义逆矩阵

在实际问题中，遇到的矩阵不一定是方阵，即使是方阵也不一定是非奇异的（即不可逆），因此需要将逆矩阵的概念进一步推广，因而学者们逐渐提出了广义逆矩阵。
        广义逆矩阵定义：设A∈C^mxn，若存在矩阵A⁺∈C^nxm，且A⁺能够同时满足以下4个公式（被称为Penrose-Moore方程），则称A⁺是A的广义逆矩阵，A⁺与A互为广义逆矩阵。（其中，A^H代表A的共轭转置矩阵）
（1）AA+A=A（1）AA^{+}A=A（1）AA+A=A （2）A+AA+=A+（2）A^{+}AA^{+}=A^{+}（2）A+AA+=A+ （3）(AA+)H=AA+（3）{(AA^{+})}^{H}=AA^{+}（3）(AA+)H=AA+ （4）(A+A)H=A+A（4）{(A^{+}A)}^{H}=A^{+}A（4）(A+A)H=A+A
        满秩分解法作为求解广义逆矩阵的一种方法，该方法认为：若A的秩为r，则A能被分解为两个秩为r的矩阵B∈C^mxr和C∈C^rxn，使得
A=BC(式1-5)A=BC\tag{式1-5}A=BC(式1-5)        其中B称为列满秩矩阵，C称为行满秩矩阵，H称为拟Hermite标准形，H∈RmxnH∈R^{mxn}H∈Rmxn，H的后m-r行上的元素全为零，H存在j₁、j₂、…、j_r列可构成m阶单位阵的前r列。存在可逆矩阵P可将矩阵A化为拟Hermite标准形H，可令前r行为C，P−1P^{-1}P−1的前r列为B，后m-r列为S，
P−1=(B,S)，B=（p1,p2,...,pr）,S=（pr+1,...,pm）P^{-1}=(B , S)，B=（p_{1},p_{2},...,p_{r}）,S=（p_{r+1},...,p_{m}）P−1=(B,S)，B=（p1,p2,...,pr）,S=（pr+1,...,pm）PA=H，即A=P−1HPA=H，即A=P^{-1}HPA=H，即A=P−1HH=(CO(m−r)×n)H=\begin{pmatrix} C\\ O_{(m-r)\times n}\\ \end{pmatrix}H=(CO(m−r)×n)A=P−1H=(B,S)(CO(m−r)×n)=BCA =P^{-1}H= (B , S)\begin{pmatrix} C\\ O_{(m-r)\times n}\\ \end{pmatrix}=BCA=P−1H=(B,S)(CO(m−r)×n)=BCP−1H=(B,S)(ErDOO)=(B,BD)=AP^{-1}H=(B , S)\begin{pmatrix} &E_{r} &D\\ &O &O\\ \end{pmatrix}=(B,BD)=AP−1H=(B,S)(ErODO)=(B,BD)=A        矩阵A的j₁、j₂、…、j_r列分别等于P−1HP^{-1}HP−1H的第j₁、j₂、…、j_r列。因此，满秩分解的具体步骤为：
    （1）用初等行变换将A化为拟Hermite标准形H，其j₁、j₂、…、j_r列构成m阶单位矩阵的前r列；
    （2）写出A有满秩分解A=BC，其中B是由A的j₁、j₂、…、j_r列构成的矩阵，C是由H的前r行构成的矩阵。

则A的广义逆矩阵可由下式计算得到：
A+=CH(CCH)−1(BHB)−1BH(式1-6)A^{+}=C^{H}{(CC^{H})}^{-1}{(B^{H}B)}^{-1}B^{H}\tag{式1-6}A+=CH(CCH)−1(BHB)−1BH(式1-6)
        将上式代入定义中的4个公式中可进行验证：
AA+A=BC⋅CH(CCH)−1(BHB)−1BH⋅BC=AAA^{+}A=BC\cdot C^{H}{(CC^{H})}^{-1}{(B^{H}B)}^{-1}B^{H}\cdot BC=AAA+A=BC⋅CH(CCH)−1(BHB)−1BH⋅BC=AA+AA+=CH(CCH)−1(BHB)−1BH⋅BC⋅CH(CCH)−1(BHB)−1BH=CH(CCH)−1(BHB)−1BH=A+A^{+}AA^{+}=C^{H}{(CC^{H})}^{-1}{(B^{H}B)}^{-1}B^{H}\cdot BC\cdot C^{H}{(CC^{H})}^{-1}{(B^{H}B)}^{-1}B^{H}=C^{H}{(CC^{H})}^{-1}{(B^{H}B)}^{-1}B^{H}=A^{+}A+AA+=CH(CCH)−1(BHB)−1BH⋅BC⋅CH(CCH)−1(BHB)−1BH=CH(CCH)−1(BHB)−1BH=A+(AA+)H=[BC⋅CH(CCH)−1(BHB)−1BH]H=]B(BHB)−1BH]H=B(BHB)−1BH=AA+{(AA^{+})}^{H}={[BC\cdot C^{H}{(CC^{H})}^{-1}{(B^{H}B)}^{-1}B^{H}]}^{H}={]B{(B^{H}B)}^{-1}B^{H}]}^{H}=B{(B^{H}B)}^{-1}B^{H}=AA^{+}(AA+)H=[BC⋅CH(CCH)−1(BHB)−1BH]H=]B(BHB)−1BH]H=B(BHB)−1BH=AA+(A+A)H=[CH(CCH)−1(BHB)−1BH⋅BC]H=[CH(CCH)−1C]H=CH(CCH)−1C=A+A{(A^{+}A)}^{H}={[C^{H}{(CC^{H})}^{-1}{(B^{H}B)}^{-1}B^{H}\cdot BC]}^{H}={[C^{H}{(CC^{H})}^{-1}C]}^{H}=C^{H}{{(CC^{H})}^{-1}}C=A^{+}A(A+A)H=[CH(CCH)−1(BHB)−1BH⋅BC]H=[CH(CCH)−1C]H=CH(CCH)−1C=A+A        广义逆矩阵存在且唯一，证明过程如下：设X、Y都是A的广义逆矩阵，则X和Y均满足定义中的4个等式，因此有
X=XAX=XAYAX=X(AY)H(AX)H=X(AXAY)H=X(AY)H=XAY=XAYAY=(XA)H(YA)HY=(YAXA)HY=(YA)HY=YAY=YX = XAX = X AYA X = X (AY)^{H} (AX)^{H} = X(AXAY)^{H} =X (AY)^{H} =X AY = XA YAY = (XA)^{H} (YA)^{H} Y = (YAXA)^{H} Y=(YA)^{H} Y = YAY = YX=XAX=XAYAX=X(AY)H(AX)H=X(AXAY)H=X(AY)H=XAY=XAYAY=(XA)H(YA)HY=(YAXA)HY=(YA)HY=YAY=Y        同时由于广义逆矩阵唯一，可知若A∈C^nxn且满秩Rank(A)=n，则
A+=A−1(式1-7)A^{+}=A^{-1}\tag{式1-7}A+=A−1(式1-7)        可将A−1A^{-1}A−1代入广义逆矩阵定义的四个公式中进行验证。
（1）AA−1A=A（1）AA^{-1}A=A（1）AA−1A=A （2）A−1AA−1=A−1（2）A^{-1}AA^{-1}=A^{-1}（2）A−1AA−1=A−1 （3）(AA−1)H=AA−1（3）{(AA^{-1})}^{H}=AA^{-1}（3）(AA−1)H=AA−1 （4）(A−1A)H=A−1A（4）{(A^{-1}A)}^{H}=A^{-1}A（4）(A−1A)H=A−1A

1.3 广义逆矩阵求解流程

图1 广义逆矩阵计算流程

矩阵的满秩分解结果并不唯一，因为设A∈Rrm×nA∈R_{r}^{m\times n}A∈Rrm×n,BCBCBC和B1C1B_{1}C_{1}B1C1都是A的满秩分解，则存在r阶可逆矩阵P，使得
A=BC,A=B1C1A=BC,A=B_{1}C_{1}A=BC,A=B1C1A=(BP)⋅(P−1C)A=(BP)\cdot (P^{-1}C)A=(BP)⋅(P−1C)令B1=BP，C1=P−1C令B_{1}=BP，C_{1}=P^{-1}C令B1=BP，C1=P−1C

2、广义逆矩阵求解代码

下面利用Matlab、C#、C++、Java、Python和R这六种高级编程语言分别实现矩阵广义逆的求解。

2.1 Matlab代码及运行结果

作为强大的矩阵运算神器，这里使用标准函数pinv测试结果，帮助中有关于该函数的详细介绍，自己就不在此献丑了。

2.1.1 使用标准函数 pinv

 A1=[1,2,3,4,5,6,7;8,9,10,11,12,13,14;15,16,17,18,19,20,21;22,23,24,25,26,27,28;29,30,31,32,33,34,35]pinv(A1)A2=[11,12,13,14;15,16,17,18;19,20,21,22;23,24,25,26;27,28,29,30;31,32,33,34]pinv(A2)

2.2 C#代码及运行结果

C#语言十分灵活，经验丰富，还是比较自信的，于是自己动手实现与MathNet库进行对比，效果还好!

2.2.1 使用MathNet库

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;namespace pinvofmatrix
{class Program{static void pinvOfMatrix(double[,] a){Matrix<double> A = DenseMatrix.OfArray(a);Console.WriteLine("原矩阵：");Console.WriteLine(A);Matrix<double> pinvA = A.PseudoInverse();Console.WriteLine("广义逆矩阵：");Console.WriteLine(pinvA);}static void Main(string[] args){double[,] a1 = new double[5, 7] { { 1, 2, 3, 4, 5, 6, 7 }, { 8, 9, 10, 11, 12, 13, 14 }, { 15, 16, 17, 18, 19, 20, 21 }, { 22, 23, 24, 25, 26, 27, 28 }, { 29, 30, 31, 32, 33, 34, 35 } }; ;pinvOfMatrix(a1);double[,] a2 = new double[6, 4] { { 11, 12, 13, 14 }, { 15, 16, 17, 18 }, { 19, 20, 21, 22 }, { 23, 24, 25, 26 }, { 27, 28, 29, 30 }, { 31, 32, 33, 34 } };pinvOfMatrix(a2);}}
}

2.2.2 自己写

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;namespace MypinvOfMatrix
{class MatrixOperator{public static double[,] Trans(double[,] a0, int m, int n)//矩阵转置函数{double[,] c = new double[n, m];for (int i = 0; i < n; i++)for (int j = 0; j < m; j++)c[i, j] = a0[j, i];return c;}public static double[,] Multi(double[,] a0, int m, int n, double[,] b0, int m2, int n2)//矩阵相乘函数{if (n != m2){Exception myException = new Exception("数组维数不匹配");throw myException;}double[,] c = new double[m, n2];for (int i = 0; i < m; i++)for (int j = 0; j < n2; j++){c[i, j] = 0;for (int k = 0; k < n; k++)c[i, j] += a0[i, k] * b0[k, j];}return c;}public static double[,] ni(double[,] a0, int m, int n)//矩阵求逆（方阵){double[,] a = new double[m, n],c = new double[m, n];if (m != n){Console.WriteLine("非方阵，不可逆！");}else{if (rank(a0, m, n) != m)Console.WriteLine("该方阵为奇异矩阵，不可逆！");else{for (int i = 0; i < m; i++){for (int j = 0; j < n; j++){a[i, j] = a0[i, j];}}for (int i = 0; i < m; i++){for (int j = 0; j < n; j++){if (i == j) c[i, j] = 1;else c[i, j] = 0;}}for (int i = 0; i < m; i++){double aii = a[i, i];//记录每行主对角线元素if (Math.Abs(aii) <= (1e-6))//如果为0{int k;for (k = i + 1; k < m; k++){double aki = a[k, i];//记录每行主对角线元素if (Math.Abs(aki) > (1e-6))//如果不为0{for (int j = 0; j < n; j++){a[i, j] = a[i, j] - a[k, j];c[i, j] = c[i, j] - c[k, j];}aii = a[i, i];for (int j = 0; j < n; j++){a[i, j] = a[i, j] / aii;c[i, j] = c[i, j] / aii;}break;}}/*if (k == m){Exception myException = new Exception("没有逆矩阵");throw myException;}* */}else//如果不为0{for (int j = 0; j < n; j++){a[i, j] = a[i, j] / aii;c[i, j] = c[i, j] / aii;}}for (int l = i + 1; l < m; l++){double ali = a[l, i];if (Math.Abs(ali) > (1e-6)){for (int j = 0; j < n; j++){a[l, j] -= ali * a[i, j];c[l, j] -= ali * c[i, j];}}}}for (int i = m - 1; i >= 0; i--){for (int k = i - 1; k >= 0; k--){double aki = a[k, i];if (Math.Abs(aki) > (1e-6)){for (int j = n - 1; j >= 0; j--){a[k, j] -= aki * a[i, j];c[k, j] -= aki * c[i, j];}}}}}}return c;}public static int rank(double[,] a0, int m, int n)//求矩阵的秩的函数{int r = 0, m1 = 0;double[,] a = new double[m, n];for (int i = 0; i < m; i++)for (int j = 0; j < n; j++){a[i, j] = a0[i, j];}if (m > n) m1 = n;else m1 = m;for (int i = 0; i < m1; i++){double aii = a[i, i];//记录每行主对角线元素if (Math.Abs(aii) <= (1e-6))//如果为0{int k;for (k = i + 1; k < m; k++){double aki = a[k, i];//记录每行对应第i行主对角线元素的元素if (Math.Abs(aki) > (1e-6))//如果不为0{for (int j = 0; j < n; j++){a[i, j] = a[i, j] - a[k, j];}double ai = a[i, i];for (int j = 0; j < n; j++)a[i, j] /= ai;break;}}}else  //主元素不为0{for (int j = 0; j < n; j++)a[i, j] /= aii;}for (int l = i + 1; l < m; l++){double ali = a[l, i];if (Math.Abs(ali) > (1e-6)){for (int j = 0; j < n; j++){a[l, j] -= ali * a[i, j];}}}}for (int i = m1 - 1; i >= 0; i--){for (int k = i - 1; k >= 0; k--){double aki = a[k, i];if (Math.Abs(aki) > (1e-6)){for (int j = m1 - 1; j >= 0; j--){a[k, j] -= aki * a[i, j];}}}}for (int i = 0; i < m; i++){for (int j = 0; j < n; j++){if (Math.Abs(a[i, j]) > (1e-6)){r++;break;}}}return r;}public static double[,] pinv(double[,] A, int m, int n)//求矩阵广义逆函数{int M;if (m > n) M = n;else M = m;double[,] E = new double[m, m];for (int i = 0; i < m; i++){for (int j = 0; j < m; j++){if (j == i) E[i, j] = 1;else E[i, j] = 0;}}for (int i = 0; i < M; i++){double aii = A[i, i];//记录每行主对角线元素if (Math.Abs(aii) <= (1e-6))//如果为0{int k;for (k = i + 1; k < m; k++){double aki = A[k, i];//记录每行对应第i行主对角线元素的元素if (Math.Abs(aki) > (1e-6))//如果不为0{for (int j = 0; j < n; j++){A[i, j] = A[i, j] - A[k, j];}for (int j = 0; j < m; j++){E[i, j] = E[i, j] - E[k, j];}double ai = A[i, i];for (int j = 0; j < n; j++){A[i, j] /= ai;}for (int j = 0; j < m; j++){E[i, j] /= ai;}break;}}}else  //主元素不为0{for (int j = 0; j < n; j++){A[i, j] /= aii;}for (int j = 0; j < m; j++){E[i, j] /= aii;}}for (int l = i + 1; l < m; l++){double ali = A[l, i];if (Math.Abs(ali) > (1e-6)){for (int j = 0; j < n; j++){A[l, j] -= ali * A[i, j];}for (int j = 0; j < m; j++){E[l, j] -= ali * E[i, j];}}else continue;//为零，则下一行}}for (int i = M - 1; i >= 0; i--){for (int k = i - 1; k >= 0; k--){double aki = A[k, i];if (Math.Abs(aki) > (1e-6)){for (int j = n - 1; j >= 0; j--){A[k, j] -= aki * A[i, j];}for (int j = m - 1; j >= 0; j--){E[k, j] -= aki * E[i, j];}}}}int rA = MatrixOperator.rank(A, m, n);for (int q = 0; q < rA; q++){if (Math.Abs(A[q, q]) <= (1e-6))//主对角线元素若为0，需要交换行或交换列{for (int i = q + 1; i < M; i++){int biaoji = -1;if (A[i, i] == 1){biaoji = i;}else if (Math.Abs(A[i, i]) <= (1e-6))//主对角线元素若为0{for (int j = 0; j < n; j++){if (A[i, j] == 1){biaoji = i;break;}}}if (biaoji > 0)//该行有1，已经找到其行号和列号{//将i行和q行交换for (int w = 0; w < n; w++){double temp = A[q, w];A[q, w] = A[i, w];A[i, w] = temp;}for (int w = 0; w < m; w++){double temp1 = E[q, w];E[q, w] = E[i, w];E[i, w] = temp1;}}}}}int z = rA;double[,] G = new double[m, n], P = new double[m, m], P1 = new double[m, m], B = new double[m, z], C = new double[z, n];for (int i = 0; i < m; i++){for (int j = 0; j < n; j++){G[i, j] = A[i, j];}}for (int i = 0; i < m; i++){for (int j = 0; j < m; j++){P[i, j] = E[i, j];}}P1 = MatrixOperator.ni(P, m, m);int rB = 0, rC = 0;for (int i = 0; i < m; i++)for (int j = 0; j < z; j++)B[i, j] = P1[i, j];rB = MatrixOperator.rank(B, m, z);for (int i = 0; i < z; i++)for (int j = 0; j < n; j++)C[i, j] = G[i, j];rC = MatrixOperator.rank(C, z, n);double[,] BT = new double[z, m], CT = new double[n, z];BT = MatrixOperator.Trans(B, m, z);CT = MatrixOperator.Trans(C, z, n);double[,] A1 = new double[n, m];A1 = MatrixOperator.Multi(MatrixOperator.Multi(CT, n, z, MatrixOperator.ni(MatrixOperator.Multi(C, z, n, CT, n, z), z, z), z, z), n, z,MatrixOperator.Multi(MatrixOperator.ni(MatrixOperator.Multi(BT, z, m, B, m, z), z, z), z, z, BT, z, m), z, m);return A1;}public static int chudengtransform(double[,] a, int m, int n)//对矩阵进行初等变换后再进行初等列变换化为行阶梯型矩阵{int zhi = 0;for (int i = 0; i < m; i++){double aii = a[i, i];//记录每行主对角线元素if (Math.Abs(aii) <= (1e-6))//如果为0{int k;for (k = i + 1; k < m; k++){double aki = a[k, i];//记录每行对应第i行主对角线元素的元素if (Math.Abs(aki) > (1e-6))//如果不为0{for (int j = 0; j < n; j++){a[i, j] = a[i, j] - a[k, j];}double ai = a[i, i];for (int j = 0; j < n; j++)a[i, j] /= ai;break;}}}else  //主元素不为0{for (int j = 0; j < n; j++)a[i, j] /= aii;}for (int l = i + 1; l < m; l++){double ali = a[l, i];if (Math.Abs(ali) > (1e-6)){for (int j = 0; j < n; j++){a[l, j] -= ali * a[i, j];}}else continue;//为零，则下一行}}for (int i = 0; i < n; i++){for (int k = i + 1; k < n; k++){double aik = a[i, k];if (Math.Abs(aik) > (1e-6)){for (int j = 0; j < m; j++){a[j, k] -= a[j, i] * aik;}}}}for (int i = 0; i < m; i++){for (int j = 0; j < n; j++){if (Math.Abs(a[i, j]) > (1e-6)){zhi++;break;}}}for (int q = 0; q < zhi; q++){if (Math.Abs(a[q, q]) <= (1e-6))//主对角线元素若为0，需要交换行或交换列{for (int i = q + 1; i < m; i++){int biaoji = -1;if (a[i, i] == 1){biaoji = i;}if (biaoji > 0)//该行有1，已经找到其行号和列号{//将i行和q行交换for (int w = 0; w < n; w++){double temp = a[q, w];a[q, w] = a[i, w];a[i, w] = temp;}//将j列和q列交换for (int w = 0; w < m; w++){double temp = a[w, q];a[w, q] = a[w, biaoji];a[w, biaoji] = temp;}}}}}return zhi;}public static void chudenghangbianhuan(double[,] A, double[,] E, int a, int b)//初等行变换化为标准阶梯型{int m = a, m1 = 0;int n = b;if (m > n) m1 = n;else m1 = m;for (int i = 0; i < m1; i++){double aii = A[i, i];//记录每行主对角线元素if (Math.Abs(aii) <= (1e-6))//如果为0{int k;for (k = i + 1; k < m; k++){double aki = A[k, i];//记录每行对应第i行主对角线元素的元素if (Math.Abs(aki) > (1e-6))//如果不为0{for (int j = 0; j < n; j++){A[i, j] = A[i, j] - A[k, j];}for (int j = 0; j < m; j++){E[i, j] = E[i, j] - E[k, j];}double ai = A[i, i];for (int j = 0; j < n; j++){A[i, j] /= ai;}for (int j = 0; j < m; j++){E[i, j] /= ai;}break;}}}else  //主元素不为0{for (int j = 0; j < n; j++){A[i, j] /= aii;}for (int j = 0; j < m; j++){E[i, j] /= aii;}}for (int l = i + 1; l < m; l++){double ali = A[l, i];if (Math.Abs(ali) > (1e-6)){for (int j = 0; j < n; j++){A[l, j] -= ali * A[i, j];}for (int j = 0; j < m; j++){E[l, j] -= ali * E[i, j];}}else continue;//为零，则下一行}}for (int i = m1 - 1; i >= 0; i--){for (int k = i - 1; k >= 0; k--){double aki = A[k, i];if (Math.Abs(aki) > (1e-6)){for (int j = n - 1; j >= 0; j--){A[k, j] -= aki * A[i, j];}for (int j = m - 1; j >= 0; j--){E[k, j] -= aki * E[i, j];}}}}int rA = MatrixOperator.rank(A, m, n);for (int q = 0; q < rA; q++){if (Math.Abs(A[q, q]) <= (1e-6))//主对角线元素若为0，需要交换行或交换列{for (int i = q + 1; i < m1; i++){int biaoji = -1;if (A[i, i] == 1){biaoji = i;}else if (Math.Abs(A[i, i]) <= (1e-6))//主对角线元素若为0{for (int j = 0; j < n; j++){if (A[i, j] == 1){biaoji = i;break;}}}if (biaoji > 0)//该行有1，已经找到其行号和列号{//将i行和q行交换for (int w = 0; w < n; w++){double temp = A[q, w];A[q, w] = A[i, w];A[i, w] = temp;}for (int w = 0; w < m; w++){double temp1 = E[q, w];E[q, w] = E[i, w];E[i, w] = temp1;}}}}}}}class Program{static void Main(string[] args){int m, n;Console.Write("矩阵的行数:");m = Convert.ToInt16(Console.ReadLine());Console.Write("矩阵的列数:");n = Convert.ToInt16(Console.ReadLine());double[,] a = new double[m, n], A = new double[m, n];Console.WriteLine("请输入矩阵元素");for (int i = 0; i < m; i++){for (int j = 0; j < n; j++)a[i, j] = double.Parse(Console.ReadLine());}Console.WriteLine("原矩阵A:");for (int i = 0; i < m; i++){for (int j = 0; j < n; j++){A[i, j] = a[i, j];Console.Write(a[i, j] + " ");if (j == (n - 1))Console.WriteLine();}}double[,] A11 = MatrixOperator.pinv(A, m, n);Console.WriteLine("广义逆矩阵pinv_A：");for (int i = 0; i < n; i++){for (int j = 0; j < m; j++){Console.Write(A11[i, j] + "  ");if (j == m - 1) Console.WriteLine();}}//接下来计算矩阵的满秩分解 求a=B*CConsole.WriteLine("满秩分解法求解广义逆矩阵pinv_A过程如下：");int z = MatrixOperator.rank(a, m, n);Console.WriteLine("A矩阵的秩为" + z);double[,] E = new double[m, m], G = new double[m, n], P = new double[m, m], P1 = new double[m, m], B = new double[m, z], C = new double[z, n];for (int i = 0; i < m; i++){for (int j = 0; j < m; j++){if (j == i) E[i, j] = 1;else E[i, j] = 0;}}MatrixOperator.chudenghangbianhuan(a, E, m, n);for (int i = 0; i < m; i++){for (int j = 0; j < n; j++){G[i, j] = a[i, j];}}for (int i = 0; i < m; i++){for (int j = 0; j < m; j++){P[i, j] = E[i, j];}}P1 = MatrixOperator.ni(P, m, m);Console.WriteLine("初等行变换后的矩阵如下：");Console.WriteLine("A----G：");for (int i = 0; i < m; i++){for (int j = 0; j < n; j++){Console.Write(G[i, j] + " ");if (j == (n - 1))Console.WriteLine();}}Console.WriteLine("E----P：");for (int i = 0; i < m; i++){for (int j = 0; j < m; j++){Console.Write(P[i, j] + " ");if (j == (m - 1))Console.WriteLine();}}Console.WriteLine("P逆:");for (int i = 0; i < m; i++){for (int j = 0; j < m; j++){Console.Write(P1[i, j] + " ");if (j == (m - 1))Console.WriteLine();}}Console.WriteLine("满秩分解（A=BC)后的矩阵如下：");int rB = 0, rC = 0;for (int i = 0; i < m; i++)for (int j = 0; j < z; j++)B[i, j] = P1[i, j];rB = MatrixOperator.rank(B, m, z);Console.WriteLine("B矩阵的秩为" + rB);for (int i = 0; i < m; i++){for (int j = 0; j < z; j++){Console.Write(B[i, j] + "  ");if (j == z - 1) Console.WriteLine();}}for (int i = 0; i < z; i++)for (int j = 0; j < n; j++)C[i, j] = G[i, j];rC = MatrixOperator.rank(C, z, n);Console.WriteLine("C矩阵的秩为" + rC);for (int i = 0; i < z; i++){for (int j = 0; j < n; j++){Console.Write(C[i, j] + "  ");if (j == n - 1) Console.WriteLine();}}double[,] BT = new double[z, m], CT = new double[n, z];BT = MatrixOperator.Trans(B, m, z);CT = MatrixOperator.Trans(C, z, n);double[,] A1 = new double[n, m];A1 = MatrixOperator.Multi(MatrixOperator.Multi(CT, n, z, MatrixOperator.ni(MatrixOperator.Multi(C, z, n, CT, n, z), z, z), z, z), n, z,MatrixOperator.Multi(MatrixOperator.ni(MatrixOperator.Multi(BT, z, m, B, m, z), z, z), z, z, BT, z, m), z, m);Console.WriteLine("矩阵A的广义逆矩阵如下：");for (int i = 0; i < n; i++){for (int j = 0; j < m; j++){Console.Write(A1[i, j] + "  ");if (j == m - 1) Console.WriteLine();}}}}
}

2.3 C++代码及运行结果

个人对C++还是比较肯定的，这门语言十分高效，为了学以致用，自己还是动手实现了一下，与Eigen库对比，结果不错，尽管这里使用了指针动态分配内存来存储二维数组，但感觉不是很好，希望给广大开发者以借鉴，最后还是推荐使用vector库存储二维矩阵。

2.3.1 使用Eigen库

#include<iostream>
#include<Eigen/Core>
#include<Eigen/SVD>   template<typename _Matrix_Type_> _Matrix_Type_ pseudoInverse(const _Matrix_Type_ &a, double epsilon = std::numeric_limits<double>::epsilon())
{Eigen::JacobiSVD< _Matrix_Type_ > svd(a, Eigen::ComputeThinU | Eigen::ComputeThinV);double tolerance = epsilon * std::max(a.cols(), a.rows()) *svd.singularValues().array().abs()(0);return svd.matrixV() *  (svd.singularValues().array().abs() > tolerance).select(svd.singularValues().array().inverse(), 0).matrix().asDiagonal() * svd.matrixU().adjoint();
}void pinvofMatrix(Eigen::MatrixXd A)
{std::cout <<"原矩阵：" << std::endl << A << std::endl;std::cout <<"广义逆矩阵：" << std::endl << pseudoInverse(A) << std::endl;
}
int main()
{Eigen::MatrixXd A1(5, 7);A1 << 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35;pinvofMatrix(A1);Eigen::MatrixXd A2(6, 4);A2 << 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34;pinvofMatrix(A2);return 0;
}

2.3.2 自己写

#include <iostream>
#include <vector>
using namespace std;struct Matrix
{public:int m = 0;int n = 0;double **value;void Init(double *a, int M, int N){if (M > 0 && N > 0){m = M;n = N;value = (double **)malloc(m * sizeof(double *));for (int i = 0; i < m; i++){value[i] = (double *)malloc(n * sizeof(double));}for (int i = 0; i < m; i++){for (int j = 0; j < n; j++){value[i][j] = *(a + i * n + j);}}}}void Init(int N){if ( N > 0){m = N;n = N;value = (double **)malloc(m * sizeof(double *));for (int i = 0; i < m; i++){value[i] = (double *)malloc(n * sizeof(double));}for (int i = 0; i < m; i++){for (int j = 0; j < n; j++){if (i == j)value[i][j] = 1;elsevalue[i][j] = 0;}}}}void Init(int M,int N){if (M > 0 && N > 0){m = M;n = N;value = (double **)malloc(m * sizeof(double *));for (int i = 0; i < m; i++){value[i] = (double *)malloc(n * sizeof(double));}for (int i = 0; i < m; i++){for (int j = 0; j < n; j++){value[i][j] = 0;}}}}void print(){if (m > 0 && n > 0){for (int i = 0; i < m; i++){for (int j = 0; j < n; j++){cout << value[i][j] << " ";}cout << endl;}}}
};
class MatrixOperator
{public:static Matrix Trans(Matrix a)//矩阵转置{Matrix mat;if (a.m > 0 && a.n > 0){mat.m = a.n;mat.n = a.m;mat.value = (double **)malloc(mat.m * sizeof(double *));for (int i = 0; i < mat.m; i++){mat.value[i] = (double *)malloc(mat.n * sizeof(double));}for (int i = 0; i < mat.m; i++){for (int j = 0; j < mat.n; j++){mat.value[i][j] = a.value[j][i];}}}return mat;}static Matrix Multiply(Matrix a, Matrix b)//两个矩阵相乘{Matrix c;if (a.n != b.m){cout << "两个矩阵无法相乘！" << endl;}else{c.m = a.m;c.n = b.n;c.value = (double **)malloc(c.m * sizeof(double *));for (int i = 0; i < c.m; i++){c.value[i] = (double *)malloc(c.n * sizeof(double));}for (int i = 0; i < c.m; i++){for (int j = 0; j < c.n; j++){c.value[i][j] = 0;for (int k = 0; k < a.n; k++)c.value[i][j] += a.value[i][k] * b.value[k][j];}}}return c;}static int Rank(Matrix a0)//求矩阵的秩{int r = 0, m1 = 0;Matrix a;a.Init(a0.m, a0.n);for (int i = 0; i < a.m; i++)for (int j = 0; j < a.n; j++)a.value[i][j] = a0.value[i][j];if (a.m > a.n)m1 = a.n;elsem1 = a.m;for (int i = 0; i < m1; i++)//从上到下，将主对角线元素化为1，左下角矩阵化为0{double aii = a.value[i][i];//记录每行主对角线元素if (abs(aii) <= (1e-6))//如果为0{int k;for (k = i + 1; k < a.m; k++){double aki = a.value[k][i];//记录每行对应第i行主对角线元素的元素if (abs(aki) >(1e-6))//如果不为0{for (int j = 0; j < a.n; j++){a.value[i][j] -= a.value[k][j];}double ai = a.value[i][i];for (int j = 0; j < a.n; j++)a.value[i][j] /= ai;break;}}}else  //主元素不为0{for (int j = 0; j < a.n; j++)a.value[i][j] /= aii;}for (int l = i + 1; l < a.m; l++){double ali = a.value[l][i];if (abs(ali) >(1e-6)){for (int j = 0; j < a.n; j++){a.value[l][j] -= ali * a.value[i][j];}}}}for (int i = m1 - 1; i >= 0; i--)//从下到上，将右上角矩阵化为0{for (int k = i - 1; k >= 0; k--){double aki = a.value[k][i];if (abs(aki) > (1e-6)){for (int j = a.n-1; j >= 0; j--){a.value[k][j] -= aki * a.value[i][j];}}}}for (int i = 0; i < a.m; i++){for (int j = 0; j < a.n; j++){if (a.value[i][i] >(1e-6)){r++; break;}}}return r;}static Matrix inv(Matrix a0)//方阵求逆{Matrix c;if (a0.m != a0.n){cout << "非方阵，不可逆！" << endl;}else{if (Rank(a0) != a0.m)cout << "该方阵为奇异矩阵，不可逆！" << endl;else{Matrix a;a.Init(a0.m, a0.n);for (int i = 0; i < a.m; i++)for (int j = 0; j < a.n; j++)a.value[i][j] = a0.value[i][j];c.m = a.m;c.n = a.n;c.value = (double **)malloc(c.m * sizeof(double *));for (int i = 0; i < c.m; i++){c.value[i] = (double *)malloc(c.n * sizeof(double));}for (int i = 0; i < c.m; i++){for (int j = 0; j < c.n; j++){if (i == j) c.value[i][j] = 1;else c.value[i][j] = 0;}}for (int i = 0; i < a.m; i++){double aii = a.value[i][i];//记录每行主对角线元素if (abs(aii) <= (1e-6))//如果为0{int k;for (k = i + 1; k < a.m; k++){double aki = a.value[k][i];//记录每行主对角线元素if (abs(aki) >(1e-6))//如果不为0{for (int j = 0; j < a.n; j++){a.value[i][j] -= a.value[k][j];c.value[i][j] -= c.value[k][j];}aii = a.value[i][i];for (int j = 0; j < a.n; j++){a.value[i][j] /= aii;c.value[i][j] /= aii;}break;}}}else//如果不为0{for (int j = 0; j < a.n; j++){a.value[i][j] /= aii;c.value[i][j] /= aii;}}for (int l = i + 1; l < a.m; l++){double ali = a.value[l][i];if (abs(ali) >(1e-6)){for (int j = 0; j < a.n; j++){a.value[l][j] -= ali * a.value[i][j];c.value[l][j] -= ali * c.value[i][j];}}}}for (int i = a.m - 1; i >= 0; i--){for (int k = i - 1; k >= 0; k--){double aki = a.value[k][i];if (abs(aki) > (1e-6)){for (int j = a.n - 1; j >= 0; j--){a.value[k][j] -= aki * a.value[i][j];c.value[k][j] -= aki * c.value[i][j];}}}}}}return c;}static Matrix pinv(Matrix a0)//矩阵求广义逆{Matrix A;A.Init(a0.m, a0.n);for (int i = 0; i < A.m; i++)for (int j = 0; j < A.n; j++)A.value[i][j] = a0.value[i][j];int M = 0;if (A.m > A.n) M = A.n;else M = A.m;Matrix E;E.Init(A.m);for (int i = 0; i < M; i++){double aii = A.value[i][i];//记录每行主对角线元素if (abs(aii) <= (1e-6))//如果为0{int k;for (k = i + 1; k < A.m; k++){double aki = A.value[k][i];//记录每行对应第i行主对角线元素的元素if (abs(aki) >(1e-6))//如果不为0{for (int j = 0; j < A.n; j++){A.value[i][j] -= A.value[k][j];}for (int j = 0; j < A.m; j++){E.value[i][j] -= E.value[k][j];}double ai = A.value[i][i];for (int j = 0; j < A.n; j++){A.value[i][j] /= ai;}for (int j = 0; j < A.m; j++){E.value[i][j] /= ai;}break;}}}else  //主元素不为0{for (int j = 0; j < A.n; j++){A.value[i][j] /= aii;}for (int j = 0; j < A.m; j++){E.value[i][j] /= aii;}}for (int l = i + 1; l < A.m; l++){double ali = A.value[l][i];if (abs(ali) >(1e-6)){for (int j = 0; j < A.n; j++){A.value[l][j] -= ali * A.value[i][j];}for (int j = 0; j < A.m; j++){E.value[l][j] -= ali * E.value[i][j];}}else continue;//为零，则下一行}}for (int i = M - 1; i >= 0; i--){for (int k = i - 1; k >= 0; k--){double aki = A.value[k][i];if (abs(aki) > (1e-6)){for (int j = A.n - 1; j >= 0; j--){A.value[k][j] -= aki * A.value[i][j];}for (int j = A.m - 1; j >= 0; j--){E.value[k][j] -= aki * E.value[i][j];}}}}int rA = MatrixOperator::Rank(A);for (int q = 0; q < rA; q++){if (abs(A.value[q][q]) <= (1e-6))//主对角线元素若为0，需要交换行或交换列{for (int i = q + 1; i < M; i++){int biaoji = -1;if (A.value[i][i] == 1){biaoji = i;}else if (abs(A.value[i][i]) <= (1e-6))//主对角线元素若为0{for (int j = 0; j < A.n; j++){if (A.value[i][j] == 1){biaoji = i;break;}}}if (biaoji > 0)//该行有1，已经找到其行号和列号{//将i行和q行交换for (int w = 0; w < A.n; w++){double temp = A.value[q][w];A.value[q][w] = A.value[i][w];A.value[i][w] = temp;}for (int w = 0; w < A.m; w++){double temp1 = E.value[q][w];E.value[q][w] = E.value[i][w];E.value[i][w] = temp1;}}}}}int z = rA, rB = 0, rC = 0;Matrix G, P, P1, B, C;G.Init(A.m, A.n);P.Init(E.m, E.n);for(int i=0;i<A.m;i++)for (int j = 0; j < A.n; j++){G.value[i][j] = A.value[i][j];}for (int i = 0; i<E.m; i++)for (int j = 0; j < E.n; j++){P.value[i][j] = E.value[i][j];}P1 = inv(P);B.Init(A.m, z);for (int i = 0; i < A.m; i++)for (int j = 0; j < z; j++)B.value[i][j] = P1.value[i][j];rB = Rank(B);C.Init(z, A.n);for (int i = 0; i < z; i++)for (int j = 0; j < A.n; j++)C.value[i][j] = G.value[i][j];rC = Rank(C);Matrix Bt = Trans(B), Ct = Trans(C),A1 = Multiply(Multiply(Ct, inv(Multiply(C, Ct))),Multiply(inv(Multiply(Bt, B)), Bt));return A1;}
};int main()
{cout << "*********************" << endl;//double mb[5][7] ={ {1, 2, 3, 4, 5, 6, 7}, { 8, 9, 10, 11, 12, 13, 14 }, {15, 16, 17, 18, 19, 20, 21}, { 22, 23, 24, 25, 26, 27, 28 }, { 29, 30, 31, 32, 33, 34, 35 } };//Matrix mB;//mB.Init((double *)mb, 5, 7);//cout << "原矩阵A1:" << endl;//mB.print();//Matrix pinv_mB = MatrixOperator::pinv(mB);//cout << "矩阵A1的广义逆矩阵:" << endl;//pinv_mB.print();cout << "*********************" << endl;double ma[6][4] = { { 11,12,13,14 },{ 15,16,17,18 },{ 19,20,21,22 },{ 23,24,25,26 },{ 27,28,29,30 },{ 31,32,33,34 } };Matrix mA;mA.Init((double *)ma, 6, 4);cout << "原矩阵A2:" << endl;mA.print();Matrix pinv_mA = MatrixOperator::pinv(mA);cout << "矩阵A2的广义逆矩阵:" << endl;pinv_mA.print();//也可考虑使用vector作为二维数组传参，更加稳定且较为方便//vector< vector<double> > arry = { { 1,2,3 },{4,5,6 } };//for (int i = 0; i < arry.size(); i++)//{//   for (int j = 0; j < arry[0].size(); j++)//    {//     cout << arry[i][j] << " ";//  }// cout << endl;//}return 0;
}

2.4 Java代码及运行结果

Java语言简单灵活，相对C++而言更易上手，实现较快，自己依靠多年实践经验动手实现，可与ujmp库对比结果，计算较为准确。

2.4.1 使用ujmp库

package com.test;import org.ujmp.core.DenseMatrix;
import org.ujmp.core.Matrix;public class pinvOfMatrix
{public static void pinvMatrix(double[][] a){Matrix A = DenseMatrix.Factory.importFromArray(a);System.out.print("A:");System.out.println(A);     Matrix pinvA = A.pinv();System.out.print("A的广义逆矩阵为:");System.out.println(pinvA);}public static void main(String[] args){double[][] a1 = {{1, 2, 3, 4, 5, 6, 7},{8, 9, 10, 11, 12, 13, 14},{15, 16, 17, 18, 19, 20, 21},{22, 23, 24, 25, 26, 27, 28},{29, 30, 31, 32, 33, 34, 35}};pinvMatrix(a1);double[][] a2={{11, 12, 13, 14},{15, 16, 17, 18},{19, 20, 21, 22},{23, 24, 25, 26},{27, 28, 29, 30},{31, 32, 33, 34}};pinvMatrix(a2);}
}

2.4.2 自己写

package com.test;public class MatrixOperator
{private static double[][] Trans(double[][] a)//矩阵转置{double[][] mat = null;if (a.length > 0 && a[0].length > 0){mat = new double[a[0].length][a.length];for (int i = 0; i < a[0].length; i++){for (int j = 0; j < a.length; j++){mat[i][j] = a[j][i];}}}return mat;}private static double[][] Multiply(double[][] a, double[][] b)//两个矩阵相乘{double[][] c = null;if (a[0].length != b.length)System.out.println("两个矩阵无法相乘！");else{c = new double[a.length][b[0].length];for (int i = 0; i < c.length; i++){for (int j = 0; j < c[0].length; j++){c[i][j] = 0;for (int k = 0; k < b.length; k++)c[i][j] += a[i][k] * b[k][j];}}}return c;}private static int Rank(double[][] a0)//求矩阵的秩{int r = 0, m1 = 0,m=a0.length,n=a0[0].length;double[][] a = new double[m][n];for (int i = 0; i < m; i++)for (int j = 0; j < n; j++)a[i][j] = a0[i][j];if (m > n)m1 = n;elsem1 = m;for (int i = 0; i < m1; i++)//从上到下，将主对角线元素化为1，左下角矩阵化为0{double aii = a[i][i];//记录每行主对角线元素if (Math.abs(aii) <= (1e-6))//如果为0{int k;for (k = i + 1; k < m; k++){double aki = a[k][i];//记录每行对应第i行主对角线元素的元素if (Math.abs(aki) >(1e-6))//如果不为0{for (int j = 0; j < n; j++){a[i][j] -= a[k][j];}double ai = a[i][i];for (int j = 0; j < n; j++)a[i][j] /= ai;break;}}}else  //主元素不为0{for (int j = 0; j < n; j++)a[i][j] /= aii;}for (int l = i + 1; l < m; l++){double ali = a[l][i];if (Math.abs(ali) >(1e-6)){for (int j = 0; j < n; j++){a[l][j] -= ali * a[i][j];}}}}for (int i = m1 - 1; i >= 0; i--)//从下到上，将右上角矩阵化为0{for (int k = i - 1; k >= 0; k--){double aki = a[k][i];if (Math.abs(aki) > (1e-6)){for (int j = m1-1; j >= 0; j--){a[k][j] -= aki * a[i][j];}}}}for (int i = 0; i < m; i++){for (int j = 0; j < n; j++){if (Math.abs(a[i][j]) >(1e-6)){r++; break;}}}return r;}private static double[][] inv(double[][] a0)//方阵求逆{double[][] c = null;int m = a0.length, n = a0[0].length;if (m != n){System.out.println("非方阵，不可逆！");}else{if (Rank(a0) != m)System.out.println("该方阵为奇异矩阵，不可逆！");else{double[][] a = new double[m][n];for (int i = 0; i < m; i++)for (int j = 0; j < n; j++)a[i][j] = a0[i][j];c = new double[m][n];for (int i = 0; i < m; i++){for (int j = 0; j < n; j++){if (i == j) c[i][j] = 1;else c[i][j] = 0;}}for (int i = 0; i < m; i++){double aii = a[i][i];//记录每行主对角线元素if (Math.abs(aii) <= (1e-6))//如果为0{int k;for (k = i + 1; k < m; k++){double aki = a[k][i];//记录每行主对角线元素if (Math.abs(aki) >(1e-6))//如果不为0{for (int j = 0; j < n; j++){a[i][j] -= a[k][j];c[i][j] -= c[k][j];}aii = a[i][i];for (int j = 0; j < n; j++){a[i][j] /= aii;c[i][j] /= aii;}break;}}}else//如果不为0{for (int j = 0; j < n; j++){a[i][j] /= aii;c[i][j] /= aii;}}for (int l = i + 1; l < m; l++){double ali = a[l][i];if (Math.abs(ali) >(1e-6)){for (int j = 0; j < n; j++){a[l][j] -= ali * a[i][j];c[l][j] -= ali * c[i][j];}}}}for (int i = m - 1; i >= 0; i--){for (int k = i - 1; k >= 0; k--){double aki = a[k][i];if (Math.abs(aki) > (1e-6)){for (int j = n - 1; j >= 0; j--){a[k][j] -= aki * a[i][j];c[k][j] -= aki * c[i][j];}}}}}}return c;}private static double[][] pinv(double[][] a0)//矩阵求广义逆{int m = a0.length,n = a0[0].length;double[][] A = new double[m][n];for (int i = 0; i < m; i++)for (int j = 0; j < n; j++)A[i][j] = a0[i][j];int M = 0;if (m > n) M = n;else M = m;double[][] E = new double[m][m];for(int i=0;i<m;i++)for(int j=0;j<m;j++)if(i==j)E[i][j] = 1;elseE[i][j] = 0;for (int i = 0; i < M; i++){double aii = A[i][i];//记录每行主对角线元素if (Math.abs(aii) <= (1e-6))//如果为0{int k;for (k = i + 1; k < m; k++){double aki = A[k][i];//记录每行对应第i行主对角线元素的元素if (Math.abs(aki) >(1e-6))//如果不为0{for (int j = 0; j < n; j++){A[i][j] -= A[k][j];}for (int j = 0; j < m; j++){E[i][j] -= E[k][j];}double ai = A[i][i];for (int j = 0; j < n; j++){A[i][j] /= ai;}for (int j = 0; j < m; j++){E[i][j] /= ai;}break;}}}else  //主元素不为0{for (int j = 0; j < n; j++){A[i][j] /= aii;}for (int j = 0; j < m; j++){E[i][j] /= aii;}}for (int l = i + 1; l < m; l++){double ali = A[l][i];if (Math.abs(ali) >(1e-6)){for (int j = 0; j < n; j++){A[l][j] -= ali * A[i][j];}for (int j = 0; j < m; j++){E[l][j] -= ali * E[i][j];}}else continue;//为零，则下一行}}for (int i = M - 1; i >= 0; i--){for (int k = i - 1; k >= 0; k--){double aki = A[k][i];if (Math.abs(aki) > (1e-6)){for (int j = n - 1; j >= 0; j--){A[k][j] -= aki * A[i][j];}for (int j = m - 1; j >= 0; j--){E[k][j] -= aki * E[i][j];}}}}int rA = Rank(A);for (int q = 0; q < rA; q++){if (Math.abs(A[q][q]) <= (1e-6))//主对角线元素若为0，需要交换行或交换列{for (int i = q + 1; i < M; i++){int biaoji = -1;if (A[i][i] == 1){biaoji = i;}else if (Math.abs(A[i][i]) <= (1e-6))//主对角线元素若为0{for (int j = 0; j < n; j++){if (A[i][j] == 1){biaoji = i;break;}}}if (biaoji > 0)//该行有1，已经找到其行号和列号{//将i行和q行交换for (int w = 0; w < n; w++){double temp = A[q][w];A[q][w] = A[i][w];A[i][w] = temp;}for (int w = 0; w < m; w++){double temp1 = E[q][w];E[q][w] = E[i][w];E[i][w] = temp1;}}}}}int z = rA;double[][] G = new double[m][n], P=new double[m][m], P1=new double[m][m], B=new double[m][z], C=new double[z][n];for(int i=0;i<m;i++)for (int j = 0; j < n; j++){G[i][j] = A[i][j];}for (int i = 0; i<m; i++)for (int j = 0; j < m; j++){P[i][j] = E[i][j];}P1 = inv(P);for (int i = 0; i < m; i++)for (int j = 0; j < z; j++)B[i][j] = P1[i][j];
//      int rB = Rank(B);for (int i = 0; i < z; i++)for (int j = 0; j < n; j++)C[i][j] = G[i][j];
//      int rC = Rank(C);double[][] Bt = Trans(B), Ct = Trans(C),A1 = Multiply(Multiply(Ct, inv(Multiply(C, Ct))),Multiply(inv(Multiply(Bt, B)), Bt));return A1;}static void Print(double[][] a)//打印矩阵函数{for(int i=0;i<a.length;i++){for(int j=0;j<a[i].length;j++)System.out.print(a[i][j]+" ");System.out.println();}}public static void main(String[] args){double[][] A1= new double[][]{{1,2,3,4,5,6,7},{8,9,10,11,12,13,14}, {15,16,17,18,19,20,21},{22,23,24,25,26,27,28},{29,30,31,32,33,34,35}};double[][] pinv_A1 = pinv(A1);double[][] A2= new double[][]{{11,12,13,14},{15,16,17,18},{19,20,21,22},{23,24,25,26},{27,28,29,30},{31,32,33,34}};double[][] pinv_A2 = pinv(A2);System.out.println("原矩阵A：");Print(A1);System.out.println("矩阵A对应的广义逆矩阵A+：");Print(pinv_A1);System.out.println("原矩阵A：");Print(A2);System.out.println("矩阵A对应的广义逆矩阵A+：");Print(pinv_A2);}
}

2.5 Python代码及运行结果

由于自己也跑过很多深度学习代码，感觉Python语言的API强大，可利用各种依赖包和库，无论是爬虫分析、人工智能、自然语言处理，这门语言灵活，易上手，因此出于好奇和自信，自己决定动手实现，学以致用，与numpy库对比，计算结果精确。

2.5.1 使用numpy库

import numpy as npdef Shucu(X,row,col): #打印二维矩阵函数for i in range(0, row):for j in range(0, col):print(" %f " % (X[i][j]), end='')print()def pinvofMatrix(X,row,col): #计算广义逆并输出pinv_X = np.linalg.pinv(X)print('原矩阵A：')Shucu(X, row, col)print('A的广义逆矩阵：')Shucu(pinv_X, col, row)if __name__ == '__main__':m ,n = 5, 7A1 = np.array([[1, 2, 3, 4, 5, 6, 7],[8, 9, 10, 11, 12, 13, 14],[15, 16, 17, 18, 19, 20, 21],[22, 23, 24, 25, 26, 27, 28],[29, 30, 31, 32, 33, 34, 35]])pinvofMatrix(A1,m,n)A2 = np.array([[11, 12, 13, 14],[15, 16, 17, 18],[19, 20, 21, 22],[23, 24, 25, 26],[27, 28, 29, 30],[31, 32, 33, 34]])m, n = 6, 4pinvofMatrix(A2, m, n)

2.5.2 自己写

import numpy as npdef Trans(A):# 矩阵转置m, n = len(A),len(A[0])At = [[0 for j in range(m)] for i in range(n)]for i in range(0, n):for j in range(0, m):At[i][j] = A[j][i]return np.array(At)def Multiply(A, B):# 矩阵相乘m1, n1 = len(A),len(A[0])m2, n2 = len(B), len(B[0])C = []if n1 != m2:print('两个矩阵无法相乘！')else:C = [[0 for j in range(n2)] for i in range(m1)]for i in range(0,m1):for j in range(0,n2):C[i][j] = 0for k in range(0,m2):C[i][j] += A[i][k] * B[k][j]return Cdef Rank(a0):# 求矩阵的秩r = 0m1 = 0m, n = len(a0), len(a0[0])a = [[0 for j in range(n)] for i in range(m)]for i in range(0,m):for j in range(0,n):a[i][j] = a0[i][j]if m > n:m1 = nelse:m1 = mfor i in range(0,m1):# 从上到下，将主对角线元素化为1，左下角矩阵化为0aii = a[i][i]; # 记录每行主对角线元素if abs(aii) <= 1e-6:# 如果为0for k in range(i+1,m):aki = a[k][i]; # 记录每行对应第i行主对角线元素的元素if abs(aii) > 1e-6:  # 如果不为0for j in range(0,n):a[i][j] -= a[k][j]ai = a[i][i]for j in range(0,n):a[i][j] /= aibreakelse:# 主元素不为0for j in range(0,n):a[i][j] /= aiifor l in range(i+1, m):ali = a[l][i]if abs(ali) > 1e-6:for j in range(0,n):a[l][j] -= ali * a[i][j]i = m1 -1while i >=0:k = i -1while k >= 0:aki = a[k][i]if abs(aki) > 1e-6:j = m1 -1while j >= 0:a[k][j] -= aki * a[i][j]j = j-1k = k -1i = i-1for i in range(0,m):for j in range(0,n):if abs(a[i][j]) > 1e-6:r = r+1breakreturn rdef inv(a0):c = []m,n = len(a0), len(a0[0])if m != n:print('非方阵，不可逆！')else:if Rank(a0) != m:print('该方阵为奇异矩阵，不可逆！')else:# 满秩a = [[0 for j in range(n)] for i in range(m)]c = [[0 for j in range(n)] for i in range(m)]for i in range(0, m):for j in range(0, n):a[i][j] = a0[i][j]if i == j:c[i][j] = 1else:c[i][j] = 0for i in range(0,m):aii = a[i][i]; # 记录每行主对角线元素if abs(aii) <= 1e-6:# 如果为0for k in range(i+1,m):aki = a[k][i]# 记录每行主对角线元素if abs(aki) > 1e-6:# 如果不为0for j in range(0,n):a[i][j] -= a[k][j]c[i][j] -= c[k][j]aii = a[i][i]for j in range(0,n):a[i][j] /= aiic[i][j] /= aiibreakelse:# 如果不为0for j in range(0,n):a[i][j] /= aiic[i][j] /= aiifor l in range(i+1,m):ali = a[l][i]if abs(ali) > 1e-6:for j in range(0,n):a[l][j] -= ali * a[i][j]c[l][j] -= ali * c[i][j]i = m -1while i >= 0:k = i -1while k >= 0:aki = a[k][i]if abs(aki) > 1e-6:j = n -1while j >= 0:a[k][j] -= aki * a[i][j]c[k][j] -= aki * c[i][j]j = j-1k = k-1i = i - 1return cdef pinv(a0):m, n = len(a0), len(a0[0])A = [[0 for j in range(n)] for i in range(m)]for i in range(0, m):for j in range(0, n):A[i][j] = a0[i][j]M = 0if m > n:M = nelse:M = mE = [[0 for j in range(m)] for i in range(m)]for i in range(0, m):for j in range(0, m):if i == j:E[i][j] = 1;else:E[i][j] = 0for i in range(0,M):aii = A[i][i]# 记录每行主对角线元素if abs(aii) <= 1e-6:# 如果为0for k in range(i+1,m):aki = A[k][i]# 记录每行对应第i行主对角线元素的元素if abs(aki) > 1e-6:# 如果不为0for j in range(0,n):A[i][j] -= A[k][j]for j in range(0,m):E[i][j] -= E[k][j]ai = A[i][i]for j in range(0,n):A[i][j] /= ai;for j in range(0,m):E[i][j] /= aibreakelse:# 主元素不为0for j in range(0,n):A[i][j] /= aiifor j in range(0,m):E[i][j] /= aiifor l in range(i+1,m):ali = A[l][i]if abs(ali) > 1e-6:for j in range(0,n):A[l][j] -= ali * A[i][j]for j in range(0,m):E[l][j] -= ali * E[i][j]else:continue# 为零，则下一行i = M -1while i >= 0:k = i-1while k >= 0:aki = A[k][i]if abs(aki) > 1e-6:j = n - 1while j >= 0:A[k][j] -= aki * A[i][j]j = j-1j = m-1while j >= 0:E[k][j] -= aki * E[i][j]j = j-1k = k-1i = i -1rA = Rank(A)for q in range(0,rA):if abs(A[q][q]) <= 1e-6:# 主对角线元素若为0，需要交换行或交换列for i in range(q+1,M):biaoji = -1if A[i][i] == 1:biaoji = 1elif abs(A[i][i]) <= 1e-6:# 主对角线元素若为0for j in range(0,n):if A[i][j] == 1:biaoji = ibreakif biaoji > 0:# 该行有1，已经找到其行号和列号# 将i行和q行交换for w in range(0,n):temp = A[q][w]A[q][w] = A[i][w]A[i][w] = tempfor w in range(0,m):temp1 = E[q][w]E[q][w] = E[i][w]E[i][w] = temp1z = rAG = [[0 for j in range(n)] for i in range(m)]P = [[0 for j in range(m)] for i in range(m)]P1 = [[0 for j in range(m)] for i in range(m)]B = [[0 for j in range(z)] for i in range(m)]C = [[0 for j in range(n)] for i in range(z)]for i in range(0,m):for j in range(0,n):G[i][j] = A[i][j]for i in range(0,m):for j in range(0,m):P[i][j] = E[i][j]P1 = inv(P)for i in range(0,m):for j in range(0,z):B[i][j] = P1[i][j]for i in range(0,z):for j in range(0,n):C[i][j] = G[i][j]Bt = Trans(B)Ct = Trans(C)A1 = Multiply(Multiply(Ct, inv(Multiply(C, Ct))), Multiply(inv(Multiply(Bt, B)), Bt))return A1def Shucu(X): #打印二维矩阵函数row,col = len(X), len(X[0])for i in range(0, row):for j in range(0, col):print(" %f " % (X[i][j]), end='')print()def pinvofMatrix(A): #计算广义逆并输出pinv_A = pinv(A)print('原矩阵A：')Shucu(A)print('A的广义逆矩阵：')Shucu(pinv_A)if __name__ == '__main__':A1 = np.array([[1, 2, 3, 4, 5, 6, 7],[8, 9, 10, 11, 12, 13, 14],[15, 16, 17, 18, 19, 20, 21],[22, 23, 24, 25, 26, 27, 28],[29, 30, 31, 32, 33, 34, 35]])A2 = np.array([[11, 12, 13, 14],[15, 16, 17, 18],[19, 20, 21, 22],[23, 24, 25, 26],[27, 28, 29, 30],[31, 32, 33, 34]])pinvofMatrix(A1)pinvofMatrix(A2)

2.6 R代码及运行结果

R语言在统计分析方面具有极大优势，自己接触较少，这里就不动手实现了，使用安装好的MASS库，调用ginv函数查看结果。

a1<- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35)
A1<- matrix(a1,5,7)
MASS::ginv(A1)
a2<- c(11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34)
A2<- matrix(a2,6,4)
MASS::ginv(A2)

3、总结

若A∈Rm×nA∈R^{m\times n}A∈Rm×n，则Ax=bAx=bAx=b的极小范数最小二乘解x唯一存在，且x=A+bx=A^{+}bx=A+b 无论是科学研究还是社会生产实践，很多问题都会涉及到矩阵，实践出真知，理论联系实际，理论指导实践。如今，广义逆矩阵已被广泛应用到优化理论、误差理论与测量平差、控制理论和系统识别等领域，尽管各种编程语言都可以调用相关的矩阵运算库来实现线性代数功能，但如果想要自己实现矩阵运算，不拘泥于现有封装库，仍需深刻把握理论数学知识，希望自己通过实践加强对理论的理解，从而为解决矩阵运算问题提供思路。（注：C语言实现可用GNU Scientific Library (GSL) 库）。
最后，向广大劳动者致敬，祝大家青年节快乐，愿不忘初心，希望对感兴趣的事情要勇于实践和探索哦！

自己动手实现广义逆矩阵求解（2022.5.4）相关推荐

高等工程数学 —— 第四章（1）线性方程组的直接解法与广义逆矩阵求解矛盾方程组
高等工程数学 -- 第四章 (1)线性方程组的直接解法与广义逆矩阵求解矛盾方程组文章目录高等工程数学 -- 第四章 (1)线性方程组的直接解法与广义逆矩阵求解矛盾方程组线性方程组的直接解法 Ga ...
五、广义逆矩阵–求解线性方程组
五.广义逆矩阵–求解线性方程组 1. 广义逆矩阵A+ 为解决各种线性方程组(系数矩阵是非方阵和方阵为奇异),将逆矩阵的概念推广到"不可逆方阵"和"长方形矩阵"上 ...
求解线性最小二乘问题的奇异值分解及广义逆法的C++实现
求解线性最小二乘问题的广义逆法的C++实现 1,功能 2,方法说明 3,函数语句与形参说明第一步,求对系数矩阵进行奇异值分解(muav函数) #include "stdlib.h" ...
2022河南萌新联赛
2022河南萌新联赛第(一)场:河南工业大学 A - Alice and Bob B - 打对子 C - 割竿榄 D - 纪念品领取 E - 聚会 F - 买车 G - 热身小游戏 H - 兴奋值 I ...
按照日期：蓝桥杯真题、洛谷题单、力扣题单汇总
2020年 2022.03.23绝世武功 2020.12.26框子求循环数组的m个最大和 2020.12.28暴力三阶幻方 2020.12.29未名湖的烦恼 2021年 2021.01.25包子凑数 ...
[更新]Win11自带邮件添加Gmail
Win11自带邮件添加Gmail 描述常见问题问题原因动手实操描述 =========== 2022.03.25 更新 ============ 补充了一些常见问题,感谢评论区的各位热心大佬, ...
暑期及9月份CSP-J1 CSP-S1初赛培训计划及学习要点
信息学奥赛一本通初赛真题解析(2022) 2022版信息学奥赛一本通初赛篇C++版一.系统讲解学习计算机基础知识.程序设计基础知识.数学问题.阅读程序.完善程序这5个模块的内容,打好坚实的基础 ...
2022年认证杯SPSSPRO杯数学建模C题(第一阶段)污水流行病学原理在新冠疫情防控方面的作用求解全过程文档及程序
2022年认证杯SPSSPRO杯数学建模 C题污水流行病学原理在新冠疫情防控方面的作用原题再现: 2019 年新型冠状病毒肺炎疫情暴发至今已过两年,新型冠状病毒历经多次变异,目前已有 11 种 ...
2022年数维杯数学建模C题电动汽车充电站的部署优化策略求解全过程文档及程序
2022年数维杯数学建模 C题电动汽车充电站的部署优化策略原题再现: 近年来,随着化石能源的逐渐枯竭和环境污染的不断加剧,电动汽车(EV)作为传统燃油车的主要替代品之一,得到了快速的发展.据国 ...

自己动手实现广义逆矩阵求解（2022.5.4）