When we have an inhomogeneous differential equation $\hat{L}f(x) = g(x)$ where $L$ is a linear operator, it makes sense to use the method of Green's functions, i.e. $\hat{L}G(x,y) = \delta(x-y)$, so that we can write $$f(x) = \int G(x,y) g(x)dy$$. This to me looks rather like having a matrix equation $Ax = b$ and writing $x = A^{-1}b$, so that $G$ is the integral kernel of the inverse operator of $\hat{L}$. The only wrinkle is that the above is not a general solution and so we need to add a homogeneous term $$f(x) = f_0(x) + \int G(x,y) g(x)dy$$ where $\hat{L} f_0(x) = 0$. This makes sense because $\hat{L}$ is not really invertible and $G$ is merely a right-inverse but not a left-inverse.
My question is this: what is the analogous relationship for matrix-vector equations with non-invertible matrices? By analogy I would expect that the general solution to $Ax=b$ would look like $$x = x_0 + Gb$$ where $AG = I$ and $Ax_0 = 0$, however I haven't been able to find any discussion of this anywhere.
In the context of elliptic differential operators, the "Green's Operator" I know of is an operator $G$ such that $$ AG=GA=\text{Id}-S $$ Where $S$ denotes the projection onto $\ker(A)$. Then if $x$ were a solution to $$ Ax=b $$ then $$ x-S(x) = GAx=Gb $$ so $$ x=Gb+S(x) $$ Of course in the context of differential operators, the underlying spaces are infinite dimensional, and there are more technicalities to concern, such as appropriate norms on the domains and properties of $A$. Nevertheless, I think the above captures the essential part if one only considers an analogy in the finite dimensional case.