线性代数 · 矩阵 · Matlab | Moore-Penrose 伪逆矩阵代码实现

背景 - Moore-Penrose 伪逆矩阵：

• 对任意矩阵 \(A\in\mathbb C^{m\times n}\) ，其 Moore-Penrose 逆矩阵 \(A^+\in\mathbb C^{n\times m}\) 存在且唯一。
  • 定义：若矩阵 G 满足 \(AGA=A,~ GAG=G,~ (AG)^H=AG,~ (GA)^H=GA\) ，则 G 是 Moore-Penrose 逆矩阵，可以记作 \(A^+\) 。
  • 性质：\(A^+\) 满足 \(A^{++}=A,~ A^{H+}=A^{+H},\) \(\mathrm{rank}(A^+)=\mathrm{rank}(A),\) \(\mathrm{Range/Null}(A^+)=\mathrm{Range/Null}(A^H)\) 。
• \(A^+\) 求解方法 1（迹方法）：
  • 令 \(A_{m\times n}\) 的秩为 r。
  • 计算 \(B=A^TA\) 。
  • 令 \(C_1=I_{n\times n}\) 。
  • for i = 2 to r-1，计算 \(C_{i+1}=(1/i)\mathrm{tr}(C_iB)I-C_iB\) 。
  • 得到 \(A^+=\frac{r}{\mathrm{tr}(C_rB)}C^rA^T\) 。
• \(A^+\) 求解方法 2（满秩分解法）：
  • 若 A=FG 是满秩分解，则 \(A^+=G^H(F^HAG^H)^{-1}F^H\) 。

代码 0（matlab 内置函数）：直接使用 pinv() 函数。

代码 1（满秩分解法）：（fullrank_decompose 代码见 代码存档）

function [B] = moore_penrose_pinv(A)
% using full rank decomposition

[m,n] = size(A);

[F,G] = fullrank_decompose(A);
B = G'*(F'*A*G')^(-1)*F';

end

代码 2（迹方法）：

function [X] = moore_penrose_pinv2(A)
% using trace method

[m,n] = size(A);
r = rank(A);

B = A'*A;
C = eye(n);
for i = 1:r-1
    C = 1/i*trace(C*B)*eye(n)-C*B;
end

X = r/trace(C*B)*C*(A');

end

测试代码：

A = [7 12 7 6 9
17 32 18 15 14
14 20 12 14 16
15 16 11 14 18];

B = moore_penrose_pinv2(A)

计算结果：

octave:5> source("my_script.m")
B =

  -0.175926   0.212037  -0.619907   0.474074
   0.112434  -0.043783   0.226257  -0.223280
   0.041005   0.120503  -0.391601   0.233862
  -0.400794  -0.100397   0.627579  -0.279365
   0.333333  -0.133333   0.066667  -0.066667


octave:6> source("my_script.m")
B2 =

  -0.175926   0.212037  -0.619907   0.474074
   0.112434  -0.043783   0.226257  -0.223280
   0.041005   0.120503  -0.391601   0.233862
  -0.400794  -0.100397   0.627579  -0.279365
   0.333333  -0.133333   0.066667  -0.066667

两份代码的计算结果一致，代码正确。