Linear Map of an ellipsoid in $\mathbb{R}^N$ into another ellipsoid in $\mathbb{R}^n$, with $nStarting from the closed set describing an ellipsoid in $\mathbb{R}^N$:
$$\Omega_x = \{  x \in \mathbb{R}^N : (x-x_0)^T\Sigma_x^{-1}(x-x_0) \leq \varepsilon^2 \}$$
where $\Sigma_x \in \mathbb{R}^{N\times N}$ is a symmetric positive definite square matrix and $x_0 \in \mathbb{R}^N$,
I need to find a way a to proove that its image through the linear map $\varphi: \mathbb{R}^N \rightarrow \mathbb{R}^n,\; n<N$, defined as the product with a full rank rectangular matrix:
$$y = \varphi(x) := Px, \; P\in\mathbb{R}^{n\times N},\;\rho(P)=n$$
is still an ellipsoid, in $\mathbb{R}^n$ of course, defined in this way:
$$\Omega_y = \varphi(\Omega_x) = \{y \in \mathbb{R}^n: y = Px, \;x \in \Omega_x\} = \{y \in \mathbb{R}^n : (y-y_0)^T\Sigma_y^{-1}(y-y_0) \leq \varepsilon^2 \}$$
with:
$$y_0 = Px_0 \in \mathbb{R}^n$$
$$\Sigma_y = P\Sigma_yP^T \in \mathbb{R}^{n\times n}$$
 A: A first change of variable, $x'=x-x_0$, is employed in order to get an ellipsoid centred in the origin:
$$\Omega_x=\{ x\in\mathbb{R}^N: x= x'+ x_0 \;,\;x'^T\Sigma^{-1}x'\le\varepsilon^2 \}$$
Moreover $y=Px=P(x'+x_0)=Px'+Px_0$
Calling $y_0:=Px_0$ then a new variable $y'$ can be defined as $y':=y-y_0 = Px'$
A second change of variable, is needed in order to obtain a sphere of $\mathbb{R}^N$, in place of the original ellipsoid. Considering the factorisation $\Sigma_x$ through the diagonal matrix of the eigenvalues $\Lambda_x$ and the orthogonal matrix with the eigenvectors as columns $\Pi_x $ (this is always possible because $\Sigma_x$ is symmetric positive definite therefore full rank) $\Sigma_x = \Pi_x \Lambda_x\Pi_x^T $:
$$x'^T\Sigma_x^{-1}x'=x'^T\Pi_x\Lambda_x^{-1}\Pi_x^Tx'=x'^T\Sigma_x^{-1}x'=x'^T\Pi_x\Lambda_x^{-\frac{1}{2}}\Lambda_x^{-\frac{1}{2}}\Pi_x^Tx'=\left(\Lambda_x^{-\frac{1}{2}}\Pi_x^Tx'\right)^T\left(\Lambda_x^{-\frac{1}{2}}\Pi_x^Tx'\right)$$ 
because $(\Lambda_x^{-1/2})^T = \Lambda_x^{-1/2}$ since it is a diagonal matrix (it is chosen as the diagonal matrix having in its main diagonal the square root ot the eigenvalues of $\Sigma_x^{-1}$).
By calling $\Lambda_x^{-1/2}\Pi_x^T = B$ and defining $x''=Bx'$ then:
$$x'^T\Sigma_x^{-1}x'= (Bx')^T(Bx') = x''^Tx''$$
$$\Omega_{x''}:=\{ x''\in\mathbb{R}^N:x''=Bx'=B(x-x_0), \; x\in\Omega_x\}$$
Since $\Sigma_x$ is positive definite all its eigenvalues are positive, then $B$ is full rank and then invertible: $x'=B^{-1}x''=\Pi_x\Lambda_x^{1/2}x''$ where $\Lambda_x^{1/2}:=\left(\Lambda_x^{-1/2}\right)^{-1}$ 
Then also $y'=Px'=PB^{-1}x''=P'x''$ where $P':=PB^{-1}$. 
If $\Omega_{y'}:=\{y'\in\mathbb{R}^n:y'=y-y_0,\; y\in\Omega_y\}$ and $\varphi':=x''\longmapsto P'x'' $ then:
$$\Omega_{y'}=\varphi(\Omega_x'')=\{y'\in\mathbb{R}^n:y'=P'x'', \;x''\in\Omega_{x''} \}$$
Now the problem is reduced to finding the linear map of a hypersphere into a lower dimension space. To this aim, the matrix $P'$ is factorised according to the singular value decomposition
$$P'=U\Sigma V^T= \begin{matrix} \end{matrix} = 
 \begin{pmatrix}u_1 & u_2& \cdots & u_n\end{pmatrix}\begin{pmatrix}
  \sigma_1 & 0 & \cdots & 0 & \cdots & 0 \\
  0 & \sigma_2 & \cdots & 0 & \cdots & 0 \\
  \vdots & \vdots   & \ddots & \vdots & & \vdots \\
  0 & 0 & \cdots & \sigma_n & \cdots & 0
 \end{pmatrix}\begin{pmatrix}v_1 \\ v_2 \\ \vdots \\ v_N\end{pmatrix}$$
where $U\in\mathbb{R}^{n\times n}, \; \Sigma\in\mathbb{R}^{n\times N}, \; V\in\mathbb{R}^{N\times N}$, $\sigma_i = \sqrt{\lambda_i(P'P'^T)}$, $u_i$ is the $i$-th normalized eigenvector of $P'P'^T$ whereas $v_i$ is the $i$-th normalized eigenvector of $P'^TP'$. 
The above relationship can be rewritten highlighting the action performed by $P'$ on $v_i$:
$$P'V=U\Sigma V^TV=U\Sigma$$
Now, considering just tbhe $i$-th column of these $n\times N$ matrices:
$ P'v_i = \sigma_iu_i \quad \forall i=1,...,n\;$ whereas $P'v_i = 0 \quad \forall i=n,...,N$.
Then $P'$ maps the base $(v_1,...,v_n)$ into $(\sigma_1 u_1,...,\sigma_n u_n)$ and nulls all the other remaining vectors of $\mathbb{R}^N$ $(v_{n+1},...,v_N)$.
Each vector of $\mathbb{R}^N$ can be represented as linear combination of the orthonormal base $(v_1,...,v_N)$; then $x''$ is mapped onto $\mathbb{R}^n$, deleting its last $N-n$ components w.r.t. such base, stretching the first $n$ components according to different weighting parameters.
Since the linear map of a sphere by a full rank diagonal matrix is an ellipsoid having as semiaxes length the reciprocal of the elements on the diagonal, while the linear map of an ellipsoid through an orthogonal matrix simply rotates without changing the shape the application of $P'$ over $\Omega_{x''}$ has to produce, eventually, an ellipsoid in $\mathbb{R}^n$.
$U$ represents the rotation matrix, aligning the ellipsoid with the eigenvectors of $P'P'^T$ (these eigenvectors are the directions of the ellipsoid main axes) and $\hat \Sigma := \left( \Sigma\Sigma^T\right)^{1/2}$ is the diagonal matrix having the on the diagonal the semiaxes length.
Then $\Omega_{y'}$ is described by the matrix $U\hat \Sigma^2 U^T = U \Sigma\Sigma^T U^T = U \Sigma V^T V \Sigma^T U^T = P'P'^T$:
$$\Omega_{y'}=\{y'\in\mathbb{R}^n: y'^T(P'P'^T)^{-1}y'\le\varepsilon^2\}$$ but $P'=PB^{-1}$ then $(P'P'^T)^{-1}=(PB^{-1}{PB^{-1}}^TP^T)^{-1}=(P\Sigma_xP^T)^{-1}$.
Coming back to the original value $y$:
$$\Omega_{y}=\{y\in\mathbb{R}^n: (y-Px_0)^T(P\Sigma_xP^T)^{-1}(y-Px_0)\le\varepsilon^2\}$$
