Let $E$ be the $n$-dimensional ellipsoid defined by $$E:=\{x \in \mathbb{R}^n: (x-c)^T A (x-c) \le 1\},$$ where $c \in \mathbb{R}^n$ is the center of the ellipsoid, and $A \in \mathbb{R}^{n \times n}$ is a symmetric positive definite matrix.
Question: How can one efficiently compute the coordinate-axis-aligned bounding box that just barely contains the ellipsoid?
For a 2D example, see the following picture:
Note: I ask this question, and answer it myself, because this question (in general form) is surprisingly absent from math.stackexchange even after 10+ years. Good answers to this question are hard to find on the internet more generally. After googling around, I eventually had to figure this out myself, and am posting here to spare future people the same trouble. Many websites discuss the question in the special case of $2D$ and $3D$, but the the format of the ellipse is given in terms of axes and angles rather than SPD matrices, and the formulas do not generalize to n-dimensions. The good answer is given by achilles hui in the comments to the following closed question: Bounding box of Ellipsoid but no proof is provided there, and the question is closed so I cannot provide the answer with proof there. Even if that question were reopened, it is focused on the 3D case with axes and angles rather than the n-dimensional case with SPD matrices.

Given vector $\rm{c} \in \Bbb R^n$ and matrix $\rm{Q} \succ \rm{O}_n$, let
$$\mathcal E := \left\{ \rm{x} \in \Bbb R^n \mid \left( \rm{x} - \rm{c} \right)^\top \rm{Q}^{-1} \left( \rm{x} - \rm{c} \right) \leq 1 \right\}$$
Let $g (\rm{x}) := \left( \rm{x} - \rm{c} \right)^\top \rm{Q}^{-1} \left( \rm{x} - \rm{c} \right)$. The vector field orthogonal to the boundary of ellipsoid $\mathcal E$ is
$$\nabla g (\rm{x}) = 2 \, \rm{Q}^{-1} \left( \rm{x} - \rm{c} \right)$$
Let us pick $i \in [n]$ and focus on the $i$-th axis. Let $\rm{P}_i := \rm{e}_i \rm{e}_i^\top$ be the projection matrix that projects onto the $i$-th axis. At the two points where ellipsoid $\mathcal E$ touches the (smallest) bounding box, we have $\rm{P}_i \nabla g (\rm{x}) = \nabla g (\rm{x})$, i.e.,
$$\left( \rm{I}_n - \rm{P}_i \right) \underbrace{ {\rm Q}^{-1} \left( \rm{x} - \rm{c} \right)}_{=: {\rm y}} = 0_n$$
Hence, $y_i$ is free and all other entries of $\rm y$ are zero, i.e., ${\rm y} = t \, {\rm e}_i$, or, ${\rm x} = {\rm c} + t \, {\rm Q} \, {\rm e}_i$. Intersecting this line with the boundary of ellipsoid $\mathcal E$, we obtain
$$t^2 = \left( {\rm e}_i^\top {\rm Q} \, {\rm e}_i \right)^{-1} = q_{ii}^{-1}$$ or, $t = \pm \frac{1}{\sqrt{q_{ii}}}$. Thus, ellipsoid $\mathcal E$ touches the (smallest) bounding box at points
$${\rm x} = {\rm c} + t \, {\rm Q} \, {\rm e}_i = {\rm c} \pm \frac{1}{\sqrt{q_{ii}}} \, {\rm Q} \, {\rm e}_i$$
and, projecting onto the $i$-th axis,
$$x_i = c_i \pm \frac{1}{\sqrt{q_{ii}}} \, {\rm e}_i^\top {\rm Q} \, {\rm e}_i = c_i \pm \frac{q_{ii}}{\sqrt{q_{ii}}} = c_i \pm \sqrt{q_{ii}}$$
Hence, the bounding box is
$$\color{blue}{\left[ c_1 - \sqrt{q_{11}}, c_1 + \sqrt{q_{11}} \right] \times \left[ c_2 - \sqrt{q_{22}}, c_2 + \sqrt{q_{22}} \right] \times \cdots \times \left[ c_n - \sqrt{q_{nn}}, c_n + \sqrt{q_{nn}} \right]}$$
Related