Let be $\vec v \in \mathbb{R^2}$ a vector and consider two different basis:
$$B_v:{\vec v_1,\vec v_2}$$ and $$B_w:{\vec w_1,\vec w_2}$$
Vector $\vec v $ can be represented in two ways:
$$\vec v =a_1 \cdot \vec v_1+a_2 \cdot \vec v_2 =x_1 \cdot \vec w_1+x_2 \cdot \vec w_2$$
or in matrix form:
$$\vec v =V \cdot a=W\cdot x$$
NOTE matrices V and W have the corresponding vectors of the basis as columns
Now suppose we know the components of $\vec v$ with respect to basis $B_v$ and we want to find the $x$ components of $\vec v$ with respect to basis $B_w$, thus:
$$V \cdot a=W\cdot x \implies x=W^{-1}Va=Ma$$
matrix $M=W^{-1}V$ is the matrix of change of basis from $B_v$ to $B_w$.
NOTE
This concept can be extended to vectors $\vec v \in \mathbb{R^n}$.
In your case, assuming $B$ as the canonical basis:
W= \begin{bmatrix}
1 & 0 \\
0 & 1
\end{bmatrix}
and
V=\begin{bmatrix}
-4 & 1 \\
1 & 0
\end{bmatrix}
thus
$$M=W^{-1}V=IV=V=\begin{bmatrix}
-4 & 1 \\
1 & 0
\end{bmatrix}$$
Ah! you are all perfectly correct till the last few equations. One of the three probe points is wrong, it's not on the ellipsoid, so it got a different value.
The correct coordinate for the first point shall be
$$
\hat{x} = [1+\sqrt 2, 1, 1, 1]\rightarrow \hat{x}^{\top}Q\hat{x}= 0.5
$$
Key is to notice, $\Sigma$ is the covariance, not standard deviation, the extant of the axis shall be sqrt ed.
In your code, this line is wrong
ellipsoid_m1 = (Sigma @ unit_sphere_ijk) + mu[:, np.newaxis]
Recall that to generate a Gaussian random variable of covariance $\Sigma$, you should not use $\Sigma$ to transform standard Gaussian. Instead you need to use its Cholesky factor to transform the standard normal variable.
$$
\Sigma=AA^T\\
z\sim N(0,I)\\
Az\sim N(0,AA^T)
$$
Really nice that you provide that well made colab notebook! it makes verifying your work super easy.
Best Answer
Notation: I will use the subscript c meaning that the object is considered in the canonical basis, while b if considered in basis b, and $x_c=(X,Y,Z)$.
The quadratic form is indipendent from basis, so we can write
$x_c^tA_cx_c=x_b^tA_bx_b$.
Using the relation $Px_c=x_b$, we get:
$(Px_c)^tA_b(Px_c)=x_c^t \left(P^t A_b P\right)x_c=x_c^tA_cx_c$
Thus, $A_b=\left(P^{-1}\right)^tA_cP^{-1}$
Note: while in this case your change of basis matrix was symmetric and orthogonal,so there was no difference in $P, P^{-1}, P^t$, you have in general to pay close attention to this differences. Those difficulties lead to the concept of contravariant and covariant transformation