# How do I calculate the Jacobian of this function?

I have two vectors $$r$$ and $$m$$. Both vectors are $$N\times1$$.

A function is calculated as -

$$F(1:N) = \phi r + (r^3 + rm^2)$$

$$F(N+1:2N) = \phi m + (m^3+mr^2)$$

F1 = (phi*r + r.^3 + r.*m.^2);
F2 = (phi*m + m.^3 + m.*r.^2);

F = [F1;F2];


I am having trouble calculating the Jacboian for this function. Here's what I did -

J11 = phi*eye(N) + diag(3*r.^2 + m.^2);
J12 = diag(2*r.*m);
J21 = diag(2*m.*r);
J22 = phi*eye(N) + diag(3*m.^2 + r.^2);
J = [J11 J12; J21 J22];


where J is the Jacobian. I believe the problem is in the part of the function where $$r$$ and $$m$$ are multiplied with each other, but I am not 100% sure about this and even if I was, I am not sure how I could fix this.

The code above is written in MATLAB. Any help would be really appreciated.

• I see that you edited your question incorporating the modified suggestions below, but didn't accept the answer or upvoted it. Should I interpret this as the suggestion didn't work for you? If that is the case, can you define the function $F$ in the mathematical language? Jun 6 at 14:28
• Oh man! I am so sorry, but I forgot to leave a comment. My code did consist of the diag terms but I had forgot to include it in the question originally. No, what you said does not work. My question's code was meant to consist of your code originally. Jun 6 at 16:29
• I have updated the question with my function definition Jun 6 at 16:34
• Math checks out, that Jacobian is correct. Why do you think that it is wrong? Jun 6 at 17:40

## 1 Answer

If you had used an indexing notation, probably you could have noticed the error yourself, and it was very hard for me to see it too. Take the line J11 = phi*eye(N) + (3*r.^2 + m.^2); for example. You are adding a matrix phi*eye(N) and a vector (3*r.^2 + m.^2) together. This is an allowed operation in MATLAB for a few years now as it has good uses in data science, however, I would have preferred the previous behaviour of throwing an error since it can cause semantic bugs like this one. For example, running [1 0; 0 1] + [2;2] returns

ans =

3     2
2     3


which may be a convenient operation for data science, but a source of bug for you.

Now, let me point out the issue by rewriting your problem in indexed notation. Let $$F = [f_0, f_1, \dots, f_n]^T$$ where $$f_i = \phi r_i + r_i^3 + r_i m_i^2$$ and $$G = [g_0, g_1, \dots, g_n]^T$$ where $$g_i = \phi m_i + m_i^3 + m_i r_i^2$$. We want to compute $$\nabla \begin{bmatrix} F \\ G \end{bmatrix}$$ which is going to be a matrix. Let's start by writing the first row, which consists of the partial derivatives of $$f_0 = \phi r_0 + r_0^3 + r_0 m_0^2$$: $$\partial_{r_0} f_0 = \phi + 3r_0^2 + m_0^2,$$ $$\partial_{r_i} f_0 = 0 \quad \text{for all } i>0,$$ $$\partial_{m_0} f_0 = 2r_0m_0,$$ and finally $$\partial_{m_i} f_0 = 0 \quad \text{for all } i>0.$$ Hence, all Jacobian computations are wrong, given that I interpreted your notation correctly. (But I don't know how to make sense of the definition of $$F$$ otherwise)

As a result, the Jacobian should be computed as

J11 = phi*eye(N) + diag(3*r.^2 + m.^2);
J12 = diag(2*r.*m);
J21 = diag(2*m.*r);
J22 = phi*eye(N) + diag(3*m.^2 + r.^2);

J = [J11 J12; J21 J22];


Depending on the number of unknowns, you may want to employ the sparse data structures MATLAB offers (speye, spdiags, etc.).

• Note that Julia has an elegant fix for this language issue: M + v throws an error, M .+ v works as in your example (singleton dimension expansion). Jun 6 at 8:18
• @FedericoPoloni That indeed is an elegant fix. MATLAB could benefit from this kind of operator overloading, but apparently there are reasons why .+ cannot be a valid combination of . and + (though I can not find any information beyond a vague statement that there are technical issues related). Jun 6 at 13:11