# Fit a quadratic form given covariant derivatives on the sphere?

I am trying to solve for a particular vector given covariant first and second derivative for a function on a sphere. If you have a quadratic form restricted to the sphere:

$f(x) = \frac{1}{2}x^T\Gamma x + p^T x + c,\, x\in S^{n-1}$

then the covariant first and second derivatives at $x$ take the form:

$\nabla f(x) = P_x(\Gamma x + p)$

$\nabla^2 f(x) = P_x\Gamma P_x-(p^Tx+x^T\Gamma x)P_x$

where $P_x = (I-xx^T)$ is the projection operator onto the tangent space at $x$. Suppose $\Gamma$ is symmetric, $\Gamma p = 0$ and $\mathrm{Tr}(\Gamma) = 0$. If we know $f, \nabla f$ and $\nabla^2 f$ at a particular point, can we solve for $p$? Is the solution unique?

Edit After some numerical experimentation I realized the expression I gave for $\nabla^2 f$ is wrong. Since the second derivative should map one point in the tangent space to another point in the tangent space, the factor of $I$ in the second term should be replaced with $P_x$. Since $\mathrm{Tr}(P_x) = n-1$, $\mathrm{Tr}(\nabla^2 f) = -(n-1)p^Tx - nx^T\Gamma x$, but now $x^T\nabla^2 f x = 0$ and the argument below fails. Technically I believe this makes $\nabla^2 f$ the Hessian and not the second covariant derivative.

-
If $x$ and $\Gamma x$ are known as well, we have $p = -\nabla^2 f(x)x + \nabla f(x) - \Gamma x$. – user1551 Feb 19 '13 at 7:00
$x$ is known but $\Gamma x$ isn't. – David Pfau Feb 19 '13 at 7:20
Why the complication with the pseudoinverse; why not just use $p=\nabla f-\Gamma x=\nabla f-\nabla^2 f x$? – joriki Feb 19 '13 at 8:39
Good question...because I like making things more complicated than they need to be? I fixed it. – David Pfau Feb 19 '13 at 9:07

Complete revision of the complete rewrite:

Let $P_x = I - x \otimes x$ denote the projection onto the tangent space $T_x S^{n-1}$ to $S^{n-1}$ at $x$. Then $$p = P_x p + (I-P_x)p = P_x p + \left\langle p,x\right\rangle x,$$ so that solving for $p$ is the same as solving for $P_x p$ and for $\left\langle p,x\right\rangle$. Before continuing, let's rewrite what we know: \begin{align*} f(x) &= \tfrac{1}{2}\left\langle x,\Gamma x\right\rangle + \left\langle p,x\right\rangle + c\\ \nabla f(x) &= P_x (\Gamma x + p)\\ \nabla^2 f(x) &= P_x \Gamma P_x - \left(\left\langle p,x \right\rangle + \left\langle x,\Gamma x\right\rangle\right) P_x. \end{align*}

Let us first solve for $\left\langle p,x \right\rangle$ and $\left\langle x,\Gamma x\right\rangle$ in terms of $c$. The condition $\operatorname{Tr}(\Gamma)= 0$ gives us: \begin{align*} 0 &= \operatorname{Tr}(\Gamma)\\ &= \operatorname{Tr}(P_x \Gamma P_x) + \operatorname{Tr}((x \otimes x)\Gamma(x \otimes x))\\ &= \operatorname{Tr}\left(\nabla^2 f(x) + (\left\langle p,x\right\rangle + \left\langle x,\Gamma x \right\rangle)P_x\right) + \operatorname{Tr}\left(\left\langle x,\Gamma x,\right\rangle x \otimes x\right)\\ &=\operatorname{Tr}(\nabla^2 f(x))+(n-1)\left\langle p,x\right\rangle + n\left\langle x,\Gamma x\right\rangle. \end{align*} Combining with our expression for $f(x)$, we find that $$\begin{pmatrix} \tfrac{1}{2} & 1 & 1 \\ n & n-1 & 0 \end{pmatrix} \begin{pmatrix} \left\langle x,\Gamma x\right\rangle \\ \left\langle p,x\right\rangle \\ c \end{pmatrix} = \begin{pmatrix} f(x) \\ -\operatorname{Tr}(\nabla^2 f(x)) \end{pmatrix},$$ which admits the solution, treating $c$ as a free variable, $$\begin{pmatrix} \left\langle x,\Gamma x\right\rangle \\ \left\langle p,x\right\rangle \end{pmatrix} = c \frac{1}{n+1} \begin{pmatrix} 2n-2 \\ -2n \end{pmatrix} + \frac{1}{n+1} \begin{pmatrix} -2n+2 & 2 \\ 2n & -1 \end{pmatrix} \begin{pmatrix} f(x) \\ -\operatorname{Tr}(\nabla^2 f(x)) \end{pmatrix}.$$

Let us now look at $P_x p$ and $P_x \Gamma x$. By the expression for $\nabla f(x)$, we at least have that $$P_x \Gamma x + P_x p = \nabla f(x).$$ Now, since $0 = \Gamma p = \Gamma P_x p + \left\langle p,x\right\rangle \Gamma x$, it follows that \begin{align*} 0 &= P_x \Gamma p \\ &= P_x \Gamma P_x P_x p + \left\langle p,x\right\rangle P_x \Gamma x\\ &= \left(\nabla^2 f(x) + \left(\left\langle p,x\right\rangle + \left\langle x,\Gamma x \right\rangle\right) P_x\right) P_x p + \left\langle p,x\right\rangle P_x \Gamma x\\ &= \left(\nabla^2 f(x) + \left(\left\langle p,x\right\rangle + \left\langle x,\Gamma x \right\rangle\right) P_x \right) P_x p + \left\langle p,x \right\rangle P_x \Gamma x. \end{align*} Hence, $$\begin{pmatrix} \nabla^2 f(x) + \left(\left\langle p,x\right\rangle + \left\langle x,\Gamma x \right\rangle\right) P_x & \left\langle p,x\right\rangle P_x \\ P_x & P_x \end{pmatrix} \begin{pmatrix} P_x p \\ P_x \Gamma x \end{pmatrix} = \begin{pmatrix} 0 \\ \nabla f(x) \end{pmatrix},$$ which is equivalent by an elementary row operation to $$\begin{pmatrix} \nabla^2 f(x) & -\left\langle x,\Gamma x\right\rangle P_x \\ P_x & P_x \end{pmatrix} \begin{pmatrix} P_x p \\ P_x \Gamma x \end{pmatrix} = \begin{pmatrix} -\left(\left\langle p,x\right\rangle + \left\langle x,\Gamma x \right\rangle\right) \nabla f(x) \\ \nabla f(x) \end{pmatrix}.$$ Thus if $\nabla^2 f(x)$ is invertible on $T_x S^{n-1}$ and if $$P_x + \left\langle x,\Gamma x\right\rangle \left(\nabla^2 f(x) \right)^{-1}$$ is invertible on $T_x S^{n-1}$, then you can use the formula for blockwise inversion to solve for $P_x p$ and $P_x \Gamma x$, given your values of $\left\langle p,x\right\rangle$ and $\left\langle x,\Gamma x\right\rangle$. I'm now at a loss, however, as to how pin down $c$, given that every bit of given information really has been used.

-
The particular form of $\Gamma$ given here is motivated by the log likelihood of the Kent distribution, which generalizes the von Mises distribution to the sphere. The conditions that $tr(\Gamma) = 0$ and $||\Gamma||<||p||/2$ I believe are required to make sure that $\nabla^2 f$ is positive definite at $p$ and negative definite at $-p$ and might not actually be needed for solving for $p$. Perhaps either $\Gamma p = 0$ or $\nabla f(x_0)$ could be useful when $c$ isn't known? – David Pfau Feb 19 '13 at 6:36
Are you assuming, then, that $p$ lies on the sphere? Because that might be useful for getting additional equations. – Branimir Ćaćić Feb 19 '13 at 6:37
No, $p$ is unconstrained except that it is in the null space of $\Gamma$. $p/||p||$ is the mean of the Kent distribution, while $||p||$ is the concentration. – David Pfau Feb 19 '13 at 6:47
Alright, I'm pretty sure that the conditions I mentioned in the first comment are superfluous and I removed them. – David Pfau Feb 19 '13 at 8:33
Also I think you have a sign flipped in $\langle x, \nabla^2 f x \rangle$. I get $-\langle p, x \rangle - \langle x, \Gamma x \rangle$. – David Pfau Feb 19 '13 at 8:52

Here is a dirty, inelegant, brute-force approach. Let $e_1=(1,0,\ldots,0)^T$, $Q$ be a real orthogonal matrix with $x$ as its first column (i.e. $Qe_1=x$), $\Lambda=Q^T\Gamma Q$ and $q=Q^Tp$. \begin{align} g := Q^T\nabla f(x) &= (I-e_1e_1^T)(\Lambda e_1 + q),\tag{1}\\ H = Q^T\left(\nabla^2 f(x)\right)Q &= (I-e_1e_1^T)\Lambda(I-e_1e_1^T)-(q^Te_1+e_1^T\Lambda e_1)I,\tag{2}\\ \Lambda q &= 0,\tag{3}\\ \operatorname{tr}\Lambda &= 0.\tag{4} \end{align} Now write $\Lambda=\begin{pmatrix}a&b^T\\b&D\end{pmatrix},\ g=\begin{pmatrix}g_1\\ \tilde{g}\end{pmatrix}$ and $q=\begin{pmatrix}u\\v\end{pmatrix}$ where $a,u\in\mathbb{R}$. Then $(1)-(2)$ give \begin{align} g_1 &= 0,\tag{5}\\ \tilde{g} &= b+v,\tag{6}\\ H &= \begin{pmatrix}-(a+u)\\&D-(a+u)I\end{pmatrix},\tag{7}\\ \end{align} From $(7)$, we obtain $h_{11}=-(a+u)$ and in turn $D$ is the trailing $(n-1)\times(n-1)$ submatrix of $H-h_{11}I$. So, by $(4)$, we get $a=-\operatorname{tr}D$ and hence $$u=\operatorname{tr}D-h_{11}.\tag{8}$$ It remains to solve for $v$. From $(3)$, we get \begin{align} au+b^Tv&=0,\\ bu+Dv&=0. \end{align} Let $w=v-\frac12\tilde{g}$. Rearrange $(6)$ to $b=\frac12\tilde{g}-w$ and substitute the result into the above, we get \begin{align*} au+\left(\frac12\tilde{g}-w\right)^T\left(w+\frac12\tilde{g}\right)&=0,\\ \left(\frac12\tilde{g}-w\right)u+D\left(w+\frac12\tilde{g}\right)&=0. \end{align*} These can be further rewritten as \begin{align} \|w\|^2 &= au + \frac14\|\tilde{g}\|^2,\tag{9}\\ (uI-D)w &= \frac12(uI+D)\tilde{g}.\tag{10} \end{align} Edit: Therefore $w$ (or $v$) is solvable if and only if $g$ satisfies $(5)$, $H$ has the block diagonal form in $(7)$, equation $(10)$ (with $u$ given by $(8)$) has a unique solution $w$ and this solution satisfies $(9)$, or $(10)$ has infinitely many solution and $$\|(uI-D)^\dagger (uI+D)\tilde{g}\|^2 \le -4(h_{11}+u)u + \|\tilde{g}\|^2. \tag{11}$$ The solution is unique iff equality holds.

-
Nice. I may attempt a cleaned-up version of this, but for now a few observations: $a = x^T\Gamma x$, $u = p^Tx$, $\mathrm{tr}D = \mathrm{tr}(I-xx^T)\Gamma(I-xx^T)$. – David Pfau Feb 19 '13 at 10:26

I've managed to derive a near-solution that doesn't require knowledge of $c$ as in Branimir's solution. At least I can reduce it to an equation that depends only on $p^T x$ and $P_x p$ and is easy to solve for one given the other. However I don't know if there is any method to solve for both simultaneously.

From the condition $\mathrm{Tr}(\Gamma) = 0$, we have:

$\mathrm{Tr}\nabla^2 f = -(n-1)(p^T x) - n(x^T\Gamma x)$

If we knew $P_x p$ and $p^Tx$ we could reconstruct $\Gamma$ as follows:

$\Gamma = P_x\Gamma P_x + xx^T\Gamma P_x + P_x\Gamma xx^T + xx^T\Gamma xx^T \\ = (\nabla^2 f+\frac{p^T x-\mathrm{Tr}\nabla^2 f}{n}P_x) + (\nabla f-P_xp)x^T + x(\nabla f - P_x p)^T - \frac{\mathrm{Tr}\nabla^2 f+(n-1)p^T x}{n}xx^T$

and plugging this in to $(\Gamma p)^T(\Gamma p)=0$ gives us:

$||\nabla^2 f P_xp + \frac{p^Tx-\mathrm{Tr}(\nabla^2f)}{n}P_xp + p^Tx(\nabla f - P_xp)||^2 \\ + ||(P_xp)^T(\nabla f - P_xp)x-\frac{(n-1)(p^Tx)+\mathrm{Tr}(\nabla^2 f)}{n}(p^Tx)x||^2 = 0$

Since the norm squared is always positive, both terms must be zero, giving us the conditions:

$\left(\nabla^2 f - \frac{n-1}{n}(p^Tx)P_x-\frac{\mathrm{Tr}(\nabla^2f)}{n}P_x\right)p + (p^Tx)\nabla f = 0$

and

$p^T(P_x p - \nabla f)+\frac{p^Tx}{n}\left(\mathrm{Tr}(\nabla^2 f)+(n-1)p^Tx\right)=0$

Given $p^Tx$, the first equation is a linear function of $P_xp$, and given $P_xp$, the second equation is a quadratic function of $p^Tx$. However, I am at a loss as to how to solve for them simultaneously, and any help is welcome. If we solve for $P_xp$ and plug the solution into the second equation, we get an equation of this form, which doesn't seem to admit an analytic solution. Numerically, Newton's method doesn't seem to converge, so now I'm at a loss.

-