On reflection, I think what's coupling the different variables together is the backtracking in the line search. I tried a similar experiment with the COIN-OR solver IPOPT
https://projects.coin-or.org/Ipopt
and got the independence I was expecting.
callIpOpt[f[x, y], {}, {{x, -2, 1, 2}, {y, -2, 1, 2}}]
{-0.618508, {x -> 0.994175, y -> 0.993418}, "cons values" -> {},
"var lb \[Lambda]s" -> {1.24674*10^-9, 1.24706*10^-9},
"var ub \[Lambda]s" -> {3.71136*10^-9, 3.70856*10^-9},
"cons \[Lambda]s" -> {}, "Solve_Succeeded"}
callIpOpt[f[x, y], {}, {{x, -2, -1, 2}, {y, -2, -1, 2}}]
{-0.0561414, {x -> -1.00898, y -> -0.957112}, "cons values" -> {},
"var lb \[Lambda]s" -> {2.5286*10^-9, 2.40285*10^-9},
"var ub \[Lambda]s" -> {8.32809*10^-10, 8.47416*10^-10},
"cons \[Lambda]s" -> {}, "Solve_Succeeded"}
callIpOpt[
f[x1, y1] +
f[x2, y2], {}, {{x1, -2, 1, 2}, {y1, -2, 1, 2}, {x2, -2, -1,
2}, {y2, -2, -1, 2}}]
{-0.674649, {x1 -> 0.994175, y1 -> 0.993418, x2 -> -1.00898,
y2 -> -0.957112}, "cons values" -> {},
"var lb \[Lambda]s" -> {1.24674*10^-9, 1.24706*10^-9, 3.76679*10^-9,
3.57946*10^-9},
"var ub \[Lambda]s" -> {3.71136*10^-9, 3.70856*10^-9, 1.24061*10^-9,
1.26237*10^-9}, "cons \[Lambda]s" -> {}, "Solve_Succeeded"}
Thanks for suggesting the formatting tool. I hadn't realized cutting and pasting dropped things, like underscores.