I'm working on a robotics problem where I'm keeping state in 3D. The state has an error bound described by a covariance matrix I'm keeping internally. A cool property of matrices is that you can take their dot product with (in 2D)
$$[cos(t) \ sin(t)]^T, \ t \in [0, 2\pi)$$
or (in 3D)
$$ [sin(\phi)cos(\theta) \ \ sin(\phi)sin(\theta) \ \ cos(\phi)]^T, \ \phi \in [0, \pi), \ \theta \in [0, 2\pi)$$
to plot an ellipse (or ellipsoid) which neatly shows how they warp a circle (or sphere). For covariance matrices this corresponds to a notion about how uncertain your state estimate is in various directions.
Right now my ellipses are showing variance, because they're from the covariance matrix, but I actually would like to plot ellipses to show standard deviation. My question is: How can I calculate the "co-standard-deviations" matrix from the covariance matrix?
I've tried just taking a matrix square root (and also an element-wise square root), but the ellipse ends up skewed.
I know this can't be right, because the spatial directions along which my sensors have maximum uncertainty should be the same, whether I'm showing standard deviation or variance.

The ellipsoid I want can be expressed as
$$ \{x | x^TP^{-1}x = 1 \} $$
where $P$ is a covariance matrix and $x$ is a vector/point. This is all the points on the boundary of an ellipsoid with principal axes along the eigenvectors of $P$, with length equal to the square roots of the corresponding eigenvalues. source
But what I'm really looking for to do my visualization is a matrix $A$, such that $A \cdot spheroid$ = this ellipsoid, whereas before I was doing $P \cdot spheroid$ = variance ellipsoid.
The answer turns out to be that you want the matrix $A \ s.t. AA^T = P$. I was missing the transpose!
Here are a few different proofs to solidify the intuition:
If $P$ is positive semidefinite and symmetric, as it must be if it's a covariance matrix, then it is diagonalizable as $P = VDV^{-1}$, where $V$ holds the eigenvectors and $D$ is a diagonal matrix holding eigenvalues.
If we choose unit length eigenvectors, noting that eigenvectors are orthogonal for a symmetric matrix, then $V$ is an orthonormal matrix, and $V^{-1} = V^T$.
Next, we note that because $D$ is diagonal, we can split it apart really easily as $D = D_{\sqrt{}}D_{\sqrt{}}$, where $D_{\sqrt{}}$ holds the square root of the eigenvalues on its diagonal. Once again, this will work because $P$ doesn't have any negative eigenvalues by property of being symmetric positive semidefinite.
Putting it all together we get
$$ P = VD_{\sqrt{}}D_{\sqrt{}}V^T = (VD_{\sqrt{}}) (D_{\sqrt{}}V^T) = (VD_{\sqrt{}})(VD_{\sqrt{}}^T)^T = (VD_{\sqrt{}})(VD_{\sqrt{}})^T = AA^T$$
So I can find $A$ with a simple eigenvalue-and-vector decomposition of $P$. Note $A$ is no longer symmetric, so I couldn't repeat this process recursively.
Next, find $P^{-1}$:
$$ P = VDV^T \rightarrow P^{-1}P = P^{-1}VDV^T \rightarrow I = P^{-1}VDV^T$$ and now start right multiplying by inverses $$ \rightarrow V = P^{-1}VD \rightarrow VD^{-1}V^T = P^{-1} $$
Next, use this $A$ in the constraint of the ellipse we want:
We're going to take a whole bunch of $VD_{\sqrt{}} \cdot spheroid\ point$ to get our $x$s that satisfy the constraint. Let's call an spheroid point $s$. The constraint now becomes:
$$ (VD_{\sqrt{}}s)^T P^{-1} (VD_{\sqrt{}}s) = 1 $$
Is this true?
$$ \rightarrow (VD_{\sqrt{}}s)^T VD^{-1}V^T (VD_{\sqrt{}}s) $$ $$ = s^T(VD_{\sqrt{}})^T VD^{-1}V^T (VD_{\sqrt{}}s) $$ $$ = s^TD_{\sqrt{}}^T V^T VD^{-1}V^T VD_{\sqrt{}}s $$ $$ = s^TD_{\sqrt{}} V^{-1} VD^{-1}V^{-1} VD_{\sqrt{}}s $$ $$ = s^TD_{\sqrt{}} D^{-1} D_{\sqrt{}}s $$ $$ = s^Ts = 1$$
the last because any point on the spheroid has magnitude 1, so its dot product with itself will be 1.
Plotting all $x = VD_{\sqrt{}}s$, I see ellipses now align!
One last proof to really solidify why this is:
$V$ is orthonormal, so it's a rotation matrix. So for the variance case we first rotate our spheroid with $V^T$ to recover the same spheroid, then scale by $D$, then rotate that ellipsoid by $V$. And in the standard deviation case, we scale the spheroid by $D_{\sqrt{}}$ and the rotate that sqrts-scaled ellipsoid by $V$.
And now we see the answer to this puzzle to make $A$ symmetric again: Use $Q = V^T$, so $A = VD_{\sqrt{}}V^T$. This extra rotation won't affect the rotation of the final ellipsoid, because, as in the variance case, it's applied to a spheroid and just returns the same spheroid.