Proof of the affine property of normal distribution for a landscape matrix The widely used/mentioned/assumed affine property of multivariate normal distributions says that: 
Given a random vector $x \in R^N$ with a multivariate normal distribution -- $x \sim N_x(\mu_x, \Sigma_x)$ -- then the random vector $y = Ax + b$ obtained by applying an affine/linear transformation to $x$ also has a normal distribution --> $y \sim N_y(A\mu_x+b, A\Sigma_x A^T)$
The above property is easy to prove if $A$ is an $N \times N$ matrix by writing $x = A^{-1}(y-b)$ and substituting it into $N_x(\mu, \Sigma)$ as shown below:
\begin{aligned}
p_y(y) & \propto p_x(A^{-1}(y-b)) \\
& \propto exp\{-0.5 \times (A^{-1}(y-b)-\mu_x)^T\Sigma_x^{-1}(A^{-1}(y-b)-\mu_x)\}\\
& = exp\{-0.5 \times (y - (A\mu_x + b))^T A^{-T}\Sigma_x ^{-1}A(y - (A\mu_x + b))\}\\
& = exp\{-0.5 \times (y - (A\mu_x + b))^T (A \Sigma_x A^T)^{-1}(y - (A\mu_x + b))\}\\
&\sim N_y(A\mu_x+b, A \Sigma_x A^T)
\end{aligned}
My questions are the following:


*

*Does the affine propert hold true even if  A is a landscape $M \times N$ matrix with $M < N$ ? (most textbooks/lecture-notes say so and many papers assume this before deriving other things)

*If the affine property is true, how do you prove it? because when A is a landscape $M \times N$ matrix with $M < N$ you cannot compute A^{-1} and hence you cannot express the random vector $x$ as  $x = A^{-1}(y-b)$

 A: One should approach this through characteristic functions. Recall that $X$ is normal $N(\mu,\Sigma)$ if and only if, for every deterministic vector $t$ of size $N\times1$,
$$
E(\mathrm e^{\mathrm it'X})=\mathrm e^{\mathrm it'\mu-t'\Sigma t/2},
$$
where $t'$ denotes the transpose of $t$. For every $(A,b)$ of compatible sizes, if $Y=AX+b$, one gets
$$
E(\mathrm e^{\mathrm it'Y})=\mathrm e^{\mathrm it'b}E(\mathrm e^{\mathrm it'AX})=\mathrm e^{\mathrm it'b}E(\mathrm e^{\mathrm is'X}),
$$
where $s'=t'A$. Since $s=(t'A)'=A't$, one sees that
$$
E(\mathrm e^{\mathrm it'Y})=\mathrm e^{\mathrm it'b}\mathrm e^{\mathrm is'\mu-s'\Sigma s/2}=\mathrm e^{\mathrm it'b+\mathrm it'A\mu-t'A\Sigma A't/2}.
$$
The RHS is the characteristic function of the normal distribution $N(b+A\mu,A\Sigma A')$, and this is enough to identify the distribution of $Y$ as such. Note that no inverse or pseudo-inverse is involved and that this applies to $Y$ of any dimension.
A: Yes it does hold for a $M \times N$ matrix
You can use the generalized inverse in place of the normal inverse, where the original inverse is now the left pseudo-inverse, such that
$A^{-1}_{left}A=I$ (Property 1: Not that this inverse is not commutative)
Applying $A_{left}^{-1}$ to both sides of $y=Ax+b$, we get
$A^{-1}_{left}y = A^{-1}_{left}(Ax+b) =  A^{-1}_{left}Ax + A^{-1}_{left}b = x + A^{-1}_{left}b$
Isolating x then gives us
$x = A^{-1}_{left}(y-b)$
which is similar to the $A^{-1}(y-b)$ which you have used in your proof above. 
If you were to trace through your proof (left as an exercise), now using the $A^{-1}_{left}$ instead of $A^{-1}$, along with property 1 above, you should arrive at the final form as you have derived.
