# 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

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