Misunderstanding of the gradient calculation I'm trying to understand how loss_metric class in dlib calculates the gradient. Here is the code(full version):
// It should be noted that the derivative of length(x-y) with respect
// to the x vector is the unit vector (x-y)/length(x-y).  If you stare
// at the code below long enough you will see that it's just an
// application of this formula.
if (x_label == y_label)
{
    // Things with the same label should have distances < dist_thresh between
    // them.  If not then we experience non-zero loss.
    if (d2 < dist_thresh-margin)//d2 - distance between x and y samples.
    {
        gm[r*temp.num_samples() + c] = 0;
    }
    else
    {
   // The whole objective function is multiplied by this to scale the loss
   // relative to the number of things in the mini-batch.
   // scale = 0.5/num_pos_samps;
        loss += scale*(d2 - (dist_thresh-margin));
        //r - x sample index, c - y sample index
        gm[r*temp.num_samples() + r] += scale/d2;
        gm[r*temp.num_samples() + c] = -scale/d2;
    }
}
else
{
    // Things with different labels should have distances > dist_thresh between
    // them.  If not then we experience non-zero loss.
    if (d2 > dist_thresh+margin || d2 > neg_thresh)
    {
        gm[r*temp.num_samples() + c] = 0;
    }
    else
    {
        loss += scale*((dist_thresh+margin) - d2);
        // don't divide by zero (or a really small number)
        d2 = std::max(d2, 0.001f);
        gm[r*temp.num_samples() + r] -= scale/d2;
        gm[r*temp.num_samples() + c] = scale/d2;
    }
}

//...
// gemm - matrix multiplication
// grad - final gradient
// grad_mul - gm
// output_tensor - output tensor of the last layer
tt::gemm(0, grad, 1, grad_mul, false, output_tensor, false); 

Let's look at the loss for the same classes (1030 line), I think the gradient of scale*(d2 - (dist_thresh-margin)) should be equal to the gradient of C1*(||X1 - X2|| - (C2-C3)), where CN - constant, and XN - output vector, so their gradient should be = C1 for X1 and -C1 for X2, but instead of this, we get some another calculations on the lines 1031, 1032, 1056. 
Same for the gradient for the different classes (starting from 1048 line).
Unfortunately, even hint in the first comment does not make it clearer. I don't have enough experience for understanding this, but I think some slightly more experienced person can show where I made a mistake. 
So, what exactly gradient formula is used here? How did we get it?
 A: I didn't read the entire code but it seems the formula used is related to the following
$$\|x-y\|=\sqrt{\sum (x_n-y_n)^2}$$
$$\frac{\partial c_1(\|x-y\|-(c_2-c_3))}{\partial x_n}=\frac{c_1}{\|x-y\|}\frac{\partial \sum (x_n-y_n)^2}{\partial x_n}=\frac{c_1}{\|x-y\|}(x_n-y_n)$$
${c_1}/{\|x-y\|}$ on the right hand side should be the scale/d2 in the code.
Update 
We have a couple of positive and negative pairs after computing the distance between each pair, we get a list of distances. The purpose of this loss is to let positive pairs have smaller distances (below some threshold) and negative pair to have larger distances (above some threshold). So it becomes rather similar to a binary classification problem (below some threshold is positive, above is negative). 
We can give the misclassified pairs a loss then that could be ||X1 - X2|| - dist_thresh. In such case, if all the positive distances are dist_thresh-0.0001 and all the negative distances are dist_thresh+0.0001, there won't be any loss anymore but the distances between positive and negative pairs are still pretty close (differed by only 0.0002). To address this we can add a margin term to the loss, which turns out to be ||X1 - X2|| - (dist_thresh +/- margin).
The scale parameter is just the reciprocal of the number of pairs, so that multiplying it gives the average loss of the batch.
