-1
$\begingroup$

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.

$\endgroup$
4
  • 2
    $\begingroup$ 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? $\endgroup$ – Abdullah Ali Sivas Jun 6 at 14:28
  • $\begingroup$ 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. $\endgroup$ – Paddy Jun 6 at 16:29
  • $\begingroup$ I have updated the question with my function definition $\endgroup$ – Paddy Jun 6 at 16:34
  • $\begingroup$ Math checks out, that Jacobian is correct. Why do you think that it is wrong? $\endgroup$ – Abdullah Ali Sivas Jun 6 at 17:40
1
$\begingroup$

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.).

$\endgroup$
2
  • 2
    $\begingroup$ 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). $\endgroup$ – Federico Poloni Jun 6 at 8:18
  • 1
    $\begingroup$ @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). $\endgroup$ – Abdullah Ali Sivas Jun 6 at 13:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.