ST语言实现Moore-Penrose伪逆

最新推荐文章于 2025-07-17 15:42:13 发布

原创最新推荐文章于 2025-07-17 15:42:13 发布 · 335 阅读

2 ·

CC 4.0 BY-SA版权

文章标签：

#人工智能 #算法

摩尔-彭罗斯伪逆是一种特殊的矩阵，它可以在原矩阵不可逆的情况下，作为一种“部分替代”。伪逆矩阵在处理线性方程组时尤为重要，特别是在方程组没有唯一解或存在多重解的情况下。通过伪逆矩阵，我们能够求解那些无法直接用常规逆矩阵求解的线性问题。

摩尔-彭罗斯伪逆的详细数学原理可以参考相关的文献和资料。这里我们不再深入探讨其推导过程，而是着重介绍其在实际应用中的实现方式。伪逆矩阵通常用于求解过定或欠定的方程组，并且在很多科学与工程领域中有着广泛的应用。

同时，LU分解作为求解线性方程组的一种常见方法，也在伪逆的计算中扮演着重要角色。LU分解可以将矩阵分解为下三角矩阵和上三角矩阵的乘积，从而便于求解。

本文将结合LU分解的思想，提出一种实现摩尔-彭罗斯伪逆矩阵的算法。通过LU分解，我们能够有效地求出伪逆矩阵，特别是当矩阵不可逆或不存在唯一解时，伪逆为求解提供了一条可行的路径。

函数Pinv具体定义如下，笔者是求9阶矩阵的伪逆，大家可根据自己所求矩阵的阶数在变量定义处做修改，或将其改为可动态变换的形式。

FUNCTION Pinv : ARRAY[0..8, 0..8] OF LREAL
VAR_INPUT
	
	A: ARRAY[0..8, 0..8] OF LREAL;
	
END_VAR
VAR
	AH: ARRAY[0..8, 0..8] OF LREAL;		// 矩阵A 的转置
	B: ARRAY[0..8, 0..8] OF LREAL;
	E, E1: ARRAY[0..8, 0..8] OF LREAL;
	temp:lreal;
	i,j,k:INT;	
	// Lu 分解
	L,U: ARRAY[0..8, 0..8] OF LREAL;
	
	x, y:ARRAY[0..8] OF LREAL;
	b1:ARRAY[0..8] OF LREAL ;
	invAHA: ARRAY[0..8, 0..8] OF LREAL;
	max_val: LREAL;
	success: BOOL;
	max_row:INT;
	P:ARRAY[0..8, 0..8] OF LREAL;
	
	n1:int:=0;
	findZeros1, findZeros2:ARRAY[0..8] OF LREAL ;
	sum1,sum2:LREAL;	
	LU1:LU;
END_VAR
VAR CONSTANT
	n:INT:=9;
	ep:LREAL:=1E-6;
END_VAR

---------------------------------------------------------------------
// Author:Tong
// Version:2.0
// date : 2024/11
// Funtion:LU分解求矩阵伪逆
success:=TRUE;
// calculate : ep * eye(n)
FOR i :=0 TO n-1 DO 
	E[i,i]:=ep;
	//L[i,i]:=1;			// 创建单位矩阵用于LU分解
	E1[i,i]:=1;
	P[i,i]:=1;
END_FOR
(* Transpose matrix A to matrix AH *)
FOR i := 0 TO n-1 DO
    FOR j := 0 TO n-1 DO
        AH[j, i] := A[i, j];
    END_FOR;
END_FOR;

// B:=A'*A+ep*I;
// 先计算 A'*A
FOR i := 0 TO n - 1 DO
    FOR j := 0 TO n - 1 DO
        temp := 0;
        FOR k := 0 TO n - 1 DO
            temp := temp + AH[i, k] * A[k, j];
        END_FOR;
        B[i, j] := temp;
    END_FOR;
END_FOR;

// 添加 ep * I
FOR i := 0 TO n - 1 DO
    B[i, i] := B[i, i] + ep;
END_FOR;
// 计算 P * B = L * U

LU1(A:= B, L1=> L, U1=> U);
//
FOR i:= 0 TO n-1 DO
	FOR j:=0 TO n-1 DO
		b1[j]:= E1[j,i];
	END_FOR
	y:=forwardSub(L, b1, n);
	x:=backwardSub(U, y, n);
	FOR j:=0 TO n-1 DO
		invAHA[j,i]:=x[j];
	END_FOR
END_FOR

FOR i := 0 TO n-1 DO
    FOR j := 0 TO n-1 DO
        Pinv[i, j] := 0; (* Initialize the element *)
        FOR k := 0 TO n-1 DO
            Pinv[i, j] := Pinv[i, j] + invAHA[i, k] * AH[k, j] ;
        END_FOR;
    END_FOR;
END_FOR;

其中用到了功能块LU，函数backwardSub以及forwardSub，其变量定义即程序部分如下所示：

FUNCTION_BLOCK LU
VAR_INPUT
	
		A: ARRAY[0..8, 0..8] OF LREAL;

END_VAR
VAR_OUTPUT
L1: ARRAY[0..8, 0..8] OF LREAL;	
U1: ARRAY[0..8, 0..8] OF LREAL;
END_VAR
VAR
	//U: ARRAY[0..8, 0..8] OF LREAL;
	L: ARRAY[0..8, 0..8] OF LREAL;	
	U: ARRAY[0..8, 0..8] OF LREAL;
	i, j, k :INT;
	n:INT:=9;		// 矩阵阶数
	ep:LREAL:=1E-6;
END_VAR

---------------------------------------------------------------
// Author:Tong
// Version:2.0
// date : 2024/11
// Funtion:LU分解

// calculate : ep * eye(n)
FOR i :=0 TO n-1 DO 
	L[i,i]:=1;			// 创建单位矩阵用于LU分解
END_FOR
U:=A;
FOR j:=0 TO n-2 DO

	FOR i:= j+1 TO n-1 DO
		IF U[j,j]=0 THEN 
			U[j,j]:=ep;
			IF i=j THEN
				L[i,j]:=1;
			ELSE
				L[i,j]:=0;
			END_IF
		END_IF
		L[i,j]:=U[i,j]/U[j,j];
		FOR k := j TO n-1 DO
			U[i,k]:= U[i,k] - L[i,j] * U[j,k];
		END_FOR	
	END_FOR
END_FOR
L1:=L;U1:=U;

FUNCTION forwardSub : ARRAY[0..8] OF LREAL
VAR_INPUT
	
	L : ARRAY[0..8,0..8] OF LREAL;
	b : ARRAY[0..8] OF LREAL;
	n:INT;
END_VAR
VAR
	i, j:INT;


END_VAR

---------------------------------------------



//  前向替换求解Ly = b
FOR i:=0 TO n-1 DO
	forwardSub[i] := b[i];
	FOR j:=0 TO i-1 DO
		forwardSub[i] := forwardSub[i] - L[i,j] * forwardSub[j];
	END_FOR
	IF L[i,i] = 0 THEN EXIT; END_IF;
	forwardSub[i] := forwardSub[i]/L[i,i];
END_FOR

FUNCTION backwardSub : ARRAY[0..8] OF LREAL
VAR_INPUT
	U : ARRAY[0..8,0..8] OF LREAL;
	y : ARRAY[0..8] OF LREAL;
	n:INT;

END_VAR
VAR
	i, j:INT;

END_VAR

----------------------------------------------------

//  后向替换求解Ux = y
FOR i:=n-1 TO 0 BY -1 DO	
	backwardSub[i]:=y[i];
	FOR j:=i+1 TO n-1 DO
		backwardSub[i]:=backwardSub[i] - U[i,j]*backwardSub[j];
	END_FOR
	IF U[i,i] = 0 THEN EXIT; END_IF;
	backwardSub[i]:=backwardSub[i]/U[i,i];
END_FOR

整体预览如下：