0

Hi I'm trying to fit a curve with gsl (Gnu Scientific Library). For the curve fit I need a Jacobian matrix, something I've never heard of before. I'm trying to wrap my head around it, but I simply don't understand how to translate the many tutorials on the topic to my curve.

A simplified version of the function is as follows:

$$ f(x, a, b, n) = \begin{cases} 0, & \text{for } x<a \newline (nx-na)/(b-a), & \text{for } a \leq x < b \newline n , & \text{for } b \leq x \newline \end{cases} $$

This function is 0 at the start, until $x$ is greater or equal to $a$, then it rises to $n$ and arrives there at $b$. If someone could explain to me how to compute the Jacobian matrix for that, it would already help a lot.

The full function is as follows: $$ f(x, a, b, c, d, n, m) = \begin{cases} 0, & \text{for } x<a \newline (nx-ca)/(b-a), & \text{for } a \leq x < b \newline (n - m)(x-b)/(c-b)+n , & \text{for } b \leq x < c \newline -m(x-c)/(d-c)+m, & \text{for } c \leq x < d \newline 0, & \text{for } d \leq x \end{cases} $$

Which would basically look like a trapezoid. See $n$ and $m$ as amplitudes of a "flat" top. and $a$ as the start, $b$ the start of the flat top, $c$ the end of the flat top, and $d$ as the end. For this I compute the following matrix with c. I'm not really sure how to translate that to math notation. I think I'm making a mistake somewhere

for (size_t i = 0; i < X; ++i) {
        double x = X[i];
        for(size_t j = 0; j < 6; j++) {
            gsl_matrix_set(J, i, j, 0);
        }
        if (x < a) {
            // do nothing, is already initialized to zero
        } else if (x < b) {
             gsl_matrix_set(J, i, 0, n*(x-b)/(pow(b-a,2)));  
             gsl_matrix_set(J, i, 1, n*(a-x)/(pow(b-a,2))); 
             gsl_matrix_set(J, i, 4, (x-start)/(b- a)); 
        } else if (x< c) {
            gsl_matrix_set(J, i, 1, (c-x)*(n-m)/(pow(c-b, 2))); 
            gsl_matrix_set(J, i, 2, (a-x)*(n-m)/(pow(c-b, 2))); 
            gsl_matrix_set(J, i, 4, (c-x)/(c- a)); 
            gsl_matrix_set(J, i, 5, (x-start)/(c- a)); 
        } else if (x < d) {
            gsl_matrix_set(J, i, 2, m*(d- x)/(pow(end-c,2))); 
            gsl_matrix_set(J, i, 3, m*(x-c)/(pow(d-c,2))); 
            gsl_matrix_set(J, i, 5, (d-x)/(d-c));
        } else {
            // do nothing, is already initialized to zero
        }
    }
}

```
Typhaon
  • 121
  • The short answer for the Jacobian matrix is: Suppose you want to change from coordinates $(x,y)$ to coordinates $(u,v)$ where the new coordinates are certain functions of the old ones. (For instance, you could have $(u,v)=(x+y,y-x)$.) Then each of $u,v$ can be differentiated with respect to each of $x,y$ which yields four derivatives in total. In the context of the multivariable chain rule, these derivatives is best arranged as a 2-by-2 matrix. (If you have more input/output variables, then the size of the matrix changes.) But how that applies to gsl is not entirely obvious... – Semiclassical Mar 14 '23 at 15:45
  • @Semiclassical I'm using gsl 1.15, and I'm simply trying to translate some code that I wrote in python. The type of fit function that I'm guessing is correct is gsl_multifit_fdfsolver. Honestly I'm really not sure if I'm using the right function for this. To be fair, the last derivative I had to compute by hand is 15 years ago, so I'm a bit lost. How would you compute that matrix specifically for the first function if x contains lets say 10 values? – Typhaon Mar 14 '23 at 15:54
  • If it's just one function $f$ which depends on $x=(x_1,x_2,\cdots, x_{10})$ then the Jacobian matrix is just a list of the 10 partial derivatives of $x$ with respect to its inputs. (This is also known as the gradient of $f$.) For each element, you differentiate with respect to one of the coordinates while leaving the rest fixed. If you had more than one output variable, then you'd need such a list of derivatives for each of the variables and this would give you a matrix. But if you don't need that, then gsl_vector is probably the better approach. – Semiclassical Mar 14 '23 at 16:01
  • Ok, I get what you mean, but I still don't really know how to compute it. I'm working on a timeseries and $a, b, c, d$ are essentially timestamps. So the y value changes either by nothing or really drastically, depending on with how much you change the parameter. How would I compute the partial derivatives? Can you show an example? – Typhaon Mar 16 '23 at 10:12

0 Answers0