I have an ellipsoid centered at $0$ (the contour of a Gaussian distribution centered at $0$ with covariance matrix $\Sigma=\Lambda^{-1}$) $$ x^\top \Lambda x = \gamma $$ and I know that the gradient at a point is given by $$ g(x) = -\Lambda x $$ Is there an expression for the tangent at a point? All I know is that $t(x)^\top g(x) = 0$.
In practice the tangent plane will be a hyperplane so there will be many vectors to choose from. However, I am looking for a tangent basis
On
Question: "Is there an expression for the tangent at a point?"
Answer: Let $k$ be the real numbers. If your ellipsoid $E$ is given as the zero set of the following polynomial
$$F(x_1,..,x_n):=a_1x_1^2+a_2x_2^2+\cdots + a_nx_n^2-a_{n+1}$$
with $a_i>0$ real numbers, it follows the tangent hyperplane $T_p(E)$ to $E$ at a point $p:=(p_1,..,p_n)$ is the zero set of the following linear polynomial:
$$H_p(E)(x_1,..,x_n):=a_1p_1(x_1-p_1)+\cdots + a_np_n(x_n-p_n).$$
Hence $T_p(E)$ is the set of points $q:=(q_1,..,q_n)\in k^n$ with $H_p(E)(q_1,..,q_n)=0$.
Note: A general formula for the (embedded) tangent plane $T_p(X)$ at the (smooth) hypersurface $X:=V(F)$ where $F(x_1,..,x_n)$ is a polynomial, is given by
$$h_p(E):=\frac{\partial F}{\partial_{x_1}}(p)(x_1-p_1)+ \cdots \frac{\partial F}{\partial_{x_n}}(p)(x_n-p_n)=0.$$
By
$$\frac{\partial F}{\partial_{x_i}}(p)$$
I mean the partial derivative of $F$ wrto the $x_i$-variable, evaluated at the point $p$. The embedded tangent space $T_p(E)$ is not a vector subspace of $k^n$.
From what I understand is that you want to obtain a basis for the hyperplane tangent to the ellipse in systematic way. The hyperplane equation is $$ p^Tg(x) = 0 $$ as you already obtained, where any vector $p\in\mathbb{R}^n$ complying with this equation will be tangent at the ellipse for that $x$. This reduces to the question:
If I have a plane equation $p^Ta=0$ with $\|a\|=1$, how can I find a basis for this subspace in the most general case?
This can be applied to your case by picking the $x\in\mathbb{R}^n$ of interest and using the procedure I provide in the following with $a=g(x)/\|g(x)\|$.
The normal of the plane is $a\in\mathbb{R}^n$. All vectors in the plane are orthogonal to $a$ since $p^Ta = 0$. In the following we will obtain an orthogonal basis for the plane. First, denote with $e_1,\dots,e_n$ the standard orthonormal basis for $\mathbb{R}^n$ (i.e. the columns of the $n\times n$ identity matrix). These comply $e_i^Te_j = 0$ if $i\neq j$ and $e_i^Te_i = 1$.
Now, lets study the following transformation: $$ A = I - 2\frac{(a - e_n)(a-e_n)^T}{(a-e_n)^T(a-e_n)} $$ where $I$ is the identity matrix. A straightforward computation shows that $A = A^T$ and $A^TA=AA = I$. Hence, $A$ is an orthogonal matrix. Now, lets compute $Aa$: $$ \begin{aligned} Aa&=a- 2\frac{(a - e_n)(a-e_n)^T}{(a-e_n)^T(a-e_n)}a\\ &= a-(a-e_n)\frac{2(a^Ta-e_n^Ta)}{a^Ta-2a^Te_n+e_n^Te_n}\\ & = a - (a-e_n)\frac{2(1-e_n^Ta)}{2-2e_n^Ta} = e_n \end{aligned} $$ Thus, $Aa = e_n$ means $a = Ae_n$ since $AA = I$. Hence, $A$ is a matrix that aligns the plane of interest to the "ground" plane orthogonal to $e_n$. This "ground" plane has basis $e_1,\dots,e_{n-1}$ and since $A$ is non-singular we can move between both planes with no problem. If you are interested in how I computed this matrix, it is a reflection over the bisecting plane between $e_n$ and $a$ as explained here.
Now, lets look at how the basis $e_1,\dots,e_{n-1}$ looks in the original space: $$ \begin{aligned} b_1 &= Ae_1\\ &\vdots\\ b_{n-1}&=Ae_{n-1} \end{aligned} $$
and verify the desired properties for the vectors $b_1,\dots,b_{n-1}$:
All $b_1,\dots,b_{n-1}$ belong to the plane $p^Ta = 0$: We check this by computing $b_i^Ta=0$ for any $i=1,\dots,n-1$. $$ b_i^Ta = \left(Ae_i\right)^T\left(Ae_n\right) = e_i^T(AA)e_n = e_i^Te_n = 0 $$ since $i\neq n$ and $AA=I, A=A^T$.
$\{b_1,\dots,b_{n-1}\}$ is a basis: We check this by showing that these are $n-1$ independent vectors. This follows simply by noticing that $b_i,b_j$ are orthogonal for $i\neq j$ by a similar argument as before: $$ b_i^Tb_j = \left(Ae_i\right)^T\left(Ab_j\right) = e_i^T(AA)b_j = e_i^Tb_j = 0 $$
Small technical detail: this cannot work if $a = e_n$ since $A$ would not be defined. However, you can pick any other $e_i$ with $i\neq n$ without loss of generality for that matter.
Summary
You have $g(x)$ as the tangent vector of the hyperellipse, and define $\hat{g}(x) = g(x)/\|g(x)\|$.
You define the matrix $$ A(x) = I - 2\frac{(\hat{g}(x) - e_n)(\hat{g}(x)-e_n)^T}{(\hat{g}(x)-e_n)^T(\hat{g}(x)-e_n)} $$
You compute the vectors $t_1(x),\dots,t_{n-1}(x)$ as: $$ \begin{aligned} t_1(x) &= A(x)e_1\\ &\vdots\\ t_{n-1}(x)&=A(x)e_{n-1} \end{aligned} $$
Then, $\{t_1(x),\dots,t_{n-1}(x)\}$ is an orthogonal basis for the tangent plane of your hyper-ellipse at the point $x$.
Last comment: note that this is not the only basis for the tangent plane.
I hope that I did not messed up anywhere and this is correct, please check. Good luck!