optimize() v = colsum() and starting values

[Demetris Christodoulou]

Dear Statalisters,

I wish to minimize the sum of a Euclidean length between two vectors, but the length is a function of an unknown alpha, where dx and dy are first differences. The step-by-step interpretation of the sum of the Euclidean length is as follows:

mata :
void arcfun(todo, p, dy, dx, v, g, H)
{
alpha = p[1]
sum = 0
for (i=1; i<=rows(dy); i++) {
sum = sum + sqrt( (dx[i] / sqrt(alpha) )^2 + ( dy[i] * sqrt(alpha))^2 )
}
v = sum
}
end
S = optimize_init()
optimize_init_evaluator(S, &arcfun())
optimize_init_evaluatortype(S, "gf0")
optimize_init_which(S,"min")
optimize_init_params(S, 1)
optimize_init_argument(S, 1, dy)
optimize_init_argument(S, 2, dx)
p = optimize(S)

The above interprets v=sum as a scalar. However I have came across two issues that are deeply puzzling to me:

(1)
The function appears to be generally insensitive to other starting values, with the sole exception of the value of 2:

optimize_init_params(S, 2)

When I specify the starting value of 2 then the optimization converges to something entirely different. Why does this happen only with the starting value of 2. By the way, note that alpha can only take values greater than 0.

(2)
In principle, I could reduce the above function using element-by-element operations, as follows:

void arcfun(todo, p, dy, dx, v, g, H)
{
alpha = p[1]
v = colsum( sqrt( ( dx:/sqrt(alpha) ):^2 + ( dy:*sqrt(alpha) ):^2 ) )
}
end

but in this case, Stata issues the warning:

numerical derivatives are approximate
flat or discontinuous region encountered

Conceptually the two implementations should be equivalent. What am I missing here? Also, surprisingly to me, when I exclude the colsum() and leave v as a vector:

v = sqrt( ( dx:/sqrt(alpha) ):^2 + ( dy:*sqrt(alpha) ):^2 )

then optimize() converges to same value as the first implementation at the top of this note, with the literal step-by-step summation. Even more puzzling, this implementation is no longer sensitive to the starting value of 2 !! Here is my string of questions:

(1) Why the first implementation is sensistive only to the starting value of 2? It is possible that it is sensitive to other starting values too that I missed, but for example it is not sensitive to 1.99 or 2.01, just to the integer value of 2.

(2) Why optimize() does not like colsum()?

(3) Does optimize() interpret the optimization of a vector as a summation by default?

(4) Why in the second implementation optimize() is no longer sensitive to the starting value of 2?

[Christophe Kolodziejczyk]

With g evaluators you don't need to compute the sum, it does it for you. If you want to compute the sum yourself use the d0 evaluator. See the help for optimize to learn the differences between evaluators. To answer point 3 of your questions g evaluators expect that you supply a vector with contributions to the objective function and will do the summation for you. d evaluators expect that you supply the value of obejctive function.

If you go for the d0 evaluator try to use quadcolsum() instead of colsum() for numerical stability and more importantly tells quadcolsum or colsum to treat missing values as such and not not as zeros. If there are missing values the sum will return a missing value, see help mata colsum() for more details. I think that it is the main issue and explains why you obtain different results depending on which evaluator you use.
You should also constrain alpha to be greater than 0 by reparametrizing the problem, i.e. estimating log(alpha). If during the optimization of your original problem alpha is close to zero, then you will have problems.

Try to rewrite your evaluator as

void arcfun(todo, p, dy, dx, v, g, H)
{
alpha = exp(p[1])
v = sqrt( ( dx:/sqrt(alpha) ):^2 + ( dy:*sqrt(alpha) ):^2 )
}

[Demetris Christodoulou]

Thanks Christophe for both tips, very useful indeed