# Computing the Hessian from the Jacobian of the gradient

I'm trying to compute the Hessian of the function,

$$f(X)=a^tX^2b$$

I used the perturbation approach to expand $$f(X+H)$$ as

$$f(X+H)=a^t(X+H)b=\underbrace{a^tX^2b}_{f} + \underbrace{a^t(XH+HX)b}_{\nabla_Hf}+a^tH^2b$$

and then took the second term and rewrote it as a matrix inner product,

\begin{align} a^t(XH+HX)b&=\text{tr}(XHba^t)+\text{tr}(HXba^t)\\ &=\text{tr}(ba^tXH)+\text{tr}(Xba^tH)\\ &=\text{tr}((ba^tX+Xba^t)H)\\ &=\langle X^tab^t+ab^tX^t,H\rangle \end{align}

thus obtaining $$\nabla f=X^tab^t+ab^tX^t$$.

Next, I used the perturbation approach again to compute the Jacobian of $$\nabla f$$,

$$(X+H)^tab^t+ab^t(X+H)^t=\underbrace{X^tab^t+ab^tX^t}_{\nabla f} + \underbrace{H^tab^t+ab^tH^t}_{J_H\nabla f}$$

Using the property for Kronecker products, $$\text{vec}(ABC)=(C^t\otimes A)\;\text{vec}(B)$$, we have

\begin{align} H^tab^t+ab^tH^t&=(ba^t\otimes I)K^{(n,n)}\text{vec}(H)+(I\otimes ab^t)K^{(n,n)}\text{vec}(H)\\ &=\big((ba^t\otimes I)+(I\otimes ab^t)\big)K^{(n,n)}\text{vec}(H) \end{align}

with $$K^{(n,n)}$$ the commutation matrix.

Therefore, the Jacobian of $$\nabla f$$, and thus presumably the Hessian of $$f$$, is

$$J\nabla f=\big((ba^t\otimes I)+(I\otimes ab^t)\big)K^{(n,n)}$$

However this matrix is not symmetric, so clearly I've either made a mistake in my calculations, or there's something about the theory that I'm not understanding.

A: With the inclusion of the commutation matrix, the Jacobian becomes symmetric.

• Your derivation of the Jacobian is almost correct, you just forgot to include the Commutation Matrix. Note that the non-vectorized Jacobian is actually a fourth-order tensor, which when flattened into a second-order tensor does not produce a symmetric matrix.
– greg
Commented Jul 21 at 18:46
• @greg could you point out what step is incorrect in my derivation of the Jacobian, I'm having trouble seeing where exactly the mistake is and where the commutation Matrix should be introduced.
– Set
Commented Jul 21 at 19:25
• In the second derivative, you forgot the transpose on the term $(X+H)^T$
– greg
Commented Jul 21 at 19:58
• @greg oh shoot you're right, and it looks like with the introduction of the commutation matrix that the Jacobian does indeed become symmetric.. is this not at odds with your comment about the flattened Jacobian not being symmetric?
– Set
Commented Jul 21 at 20:49
• i stand corrected
– greg
Commented Jul 22 at 2:10

$$\def\qty#1{\left(#1\right)}$$ $$\def\d{{\sf d}}$$

An alternative approach is to write $$f = a^j {X^i}_j {X^k}_i b_k$$ Then take the differential $$\d f = a^j \qty{\d{X^i}_j} {X^k}_i b_k + a^j {X^i}_j \qty{\d{X^k}_i} b_k$$ Write $$\d{X^i}_j =\delta^i_p \delta^q_j \d{X^p}_q \\ \d{X^k}_i = \delta^k_p \delta^q_i \d{X^p}_q$$ so that $$\d f = \qty{a^q {X^k}_p b_k + a^j {X^q}_j b_p} \d{X^p}_q$$ and $$\frac{\partial f}{\partial {X_p}^q}=a^q {X^k}_p b_k + a^j {X^q}_j b_p$$ Now write the next differential $$\d\qty{\frac{\partial f}{\partial {X_p}^q}} = a^q \qty{\d{X^k}_p} b_k + a^j \qty{\d{X^q}_j} b_p$$ As before $$\d{X^k}_p = \delta^k_l \delta_p^m \d{X^l}_m \\ \d{X^q}_j = \delta_l^q \delta_j^m \d{X^l}_m$$ so that \begin{align} \d\qty{\frac{\partial f}{\partial {X_p}^q}} & = \qty{a^q \delta^k_l \delta_p^m b_k + a^j \delta_l^q \delta_j^m b_p} \d{X^l}_m \\ & =\qty{a^q \delta_p^m b_l + a^m \delta_l^q b_p} \d{X^l}_m \end{align} Therefore the second partial derivative is $$\frac{\partial^2 f}{\partial {X_p}^q \partial {X_l}^m} = a^q \delta_p^m b_l + a^m \delta_l^q b_p$$ An example for $$a, b \in \mathbb{R}^{3 \times 1}$$ and $$X \in \mathbb{R}^{3 \times 3}$$ is shown below

• I'm not familiar with some of your notation, what do the sub and superscripts denote? Is this Einstein notation?
– Set
Commented Jul 21 at 16:43
• A matrix $\bf X$ can be written as ${\bf X }= {X^j}_i \; {\bf e}_j \otimes {\bf e}^i$ so ${X^j}_i$ are the components of the matrix given a basis. The Einstein summation convention is used to simplify the expressions. Commented Jul 21 at 16:47

You can solve these kinds of derivatives pretty easily using tensorgrad:

    i = symbols("i")
a = Variable("a", i)
b = Variable("b", i)
X = Variable("X", i, j=i)
graph = F.graph('a -i- X1 -j-i- X2 -j-i- b', a=a, X1=X, X2=X, b=b)
expr = Derivative(Derivative(graph, X), X)
save_steps(expr)


This uses Penrose graphical notation to take the derivative, and even prints out the step-by-step to do so (see below.)

The only tricky part is that the final output is an order 4 tensor, so it's not so easy to "turn it back" into algebraic notation, which usually only allows vectors and matrices.

But we can write it in index notation as $$\delta_{j',i''} a_{i'} b_{j''} + \delta_{i',j''} a_{i''} b_{j'}$$, which seems to be the same as what Ted gets.

• It is not clear how to translate your set of images into a usable mathematical expression Commented Jul 21 at 22:12
• I know. We need a better reference for how it all works. The best I currently know is Penrose's own figures: en.wikipedia.org/wiki/… . But it's basically the index method just drawn as a graph. Commented Jul 21 at 22:15