Assign FindRoot output vector to a vector instead of element by element

Question

Using FindRoot to solve a multivariate equation system, assigning the output to variables one by one works, but assigning the output vector to a vector gives the error (MWE below) SetDelayed::shape: Lists {xsola,ysola} and {x,y}/. FindRoot[{D[(x-1)^2]==0,D[(y-2)^2]==0},{{x,0.3},{y,0.4}}] are not the same shape. The output of FindRoot is {1., 2.} and the list {xsola,ysola} prints as {xsola,ysola}. To me, these lists do look the same shape, but seemingly not to Mathematica.

How to assign a vector output of FindRoot to a vector?

MWE:

Clear[x, y, xsol, ysol, xsola, ysola]
xsol := x /. 
  FindRoot[{D[(x - 1)^2] == 0, 
    D[(y - 2)^2] == 0}, {{x, 0.3}, {y, 0.4}}]
ysol := y /. 
   FindRoot[{D[(x - 1)^2] == 0, 
    D[(y - 2)^2] == 0}, {{x, 0.3}, {y, 0.4}}]
{xsola, ysola} := {x, y} /. 
  FindRoot[{D[(x - 1)^2] == 0, 
    D[(y - 2)^2] == 0}, {{x, 0.3}, {y, 0.4}}]
{x, y} /. 
 FindRoot[{D[(x - 1)^2] == 0, D[(y - 2)^2] == 0}, {{x, 0.3}, {y, 0.4}}]
{xsola, ysola}
{xsol, ysol, xsola, ysola}

The reason to assign the output of FindRoot is to use it in fixed point calculations, similarly to this question: wrong use of _?NumberQ? The equations represented in the MWE as D[(x-1)^2]==0 etc will be solutions of integrals, often not in closed form, depending on parameters. Integrate has trouble with them and NIntegrate does not work because of the non-numeric parameters.

A related question More issues Integrate, NIntegrate, FindRoot suggests using == under FindRoot, which I am already doing.

More background in case anyone is curious and has too much time on their hands: a bigger piece of code (possibly with errors) into which I am trying to fit the FindRoot is

Clear[d, s, vi, vj, vlb, mui, muj, pi, pj, profit, ci, cj, pbr, \
pisol, pjsol, pis, pjs]
$Assumptions = 
  0 < {s, pi, pj, pis, pjs} < 1/2 && 0 < {mui, muj, vlb, vi, vj} < 1 &&
    0 <= {ci, cj} < 1;
vlb[pi_, pj_] = 
 1 + pi - pj - 
  Sqrt[2] Sqrt[s]; ci = 0.; cj = 0.05; mui = 0.5; muj = 0.5; 
d[pi_, pj_, s_, pis_] = 
 mui*(1 - Max[pi, vlb[pi, pj]]) + 
  mui*Integrate[vi - pi + pj, {vi, Max[pi, 0], Max[pi, vlb[pi, pj]]}, 
    Assumptions -> 0 < {s, pi, pj, pis, pjs} < 1/2 && 0 < {vlb} < 1] +
   muj*Integrate[(1 - Max[pi, vj - pj + pi]), {vj, 0, vlb[pj, pis]}, 
    Assumptions -> 0 < {s, pi, pj, pis, pjs} < 1/2 && 0 < {vlb} < 1]; 
profit[pi_, pj_, s_, pis_, c_] = (pi - c)*d[pi, pj, s, pis]; 
pbr[pj_Real?Positive, s_Real?Positive, pis_Real?Positive, 
  c_Real?Positive] = 
 NArgMax[{profit[pi, pj, s, pis, c], c < pi < 1/2}, pi]
pisol[s_Real?Positive] := 
 pi /. FindRoot[{pbr[pj, s, pis, ci] == pi, pbr[pi, s, pjs, cj] == pj,
    pi == pis, 
     pj == pjs}, {{pi, 0.3}, {pj, 0.4}, {pis, 0.35}, {pjs, 0.45}}]
pjsol[s_Real?Positive] := 
 pj /. FindRoot[{pbr[pj, s, pis, ci] == pi, pbr[pi, s, pjs, cj] == pj,
     pi == pis, 
   pj == pjs}, {{pi, 0.3}, {pj, 0.4}, {pis, 0.35}, {pjs, 0.45}}]
pisol[0.25]
Plot[{pisol[s], pjsol[s]}, {s, 0., 0.5}](*Does not solve in 45 min. *)

I have tried to modify this to use derivatives instead of NArgMax:

(*Using derivatives instead of NArgMax*)
Clear[d, s, vi, vj, vlb, \
mui, muj, pi, pj, profit, ci, cj, pfoc, pisol, pjsol, pis, pjs]
$Assumptions = 
  0 < {s, pi, pj, pis, pjs} < 1/2 && 0 < {mui, muj, vlb, vi, vj} < 1 &&
    0 <= {ci, cj} < 1;
vlb[pi_, pj_] = 
 1 + pi - pj - 
  Sqrt[2] Sqrt[s]; ci = 0.; cj = 0.05; mui = 0.5; muj = 0.5; 
d[pi_, pj_, s_, pis_] = 
 mui*(1 - Max[pi, vlb[pi, pj]]) + 
  mui*Integrate[vi - pi + pj, {vi, Max[pi, 0], Max[pi, vlb[pi, pj]]}, 
    Assumptions -> 0 < {s, pi, pj, pis, pjs} < 1/2 && 0 < {vlb} < 1] +
   muj*Integrate[(1 - Max[pi, vj - pj + pi]), {vj, 0, vlb[pj, pis]}, 
    Assumptions -> 0 < {s, pi, pj, pis, pjs} < 1/2 && 0 < {vlb} < 1];
d[0.2, 0.3, 0.1, 0.25]
profit[pi_, pj_, s_, pis_, c_] = (pi - c)*d[pi, pj, s, pis]; 
pfoc[pi_Real?Positive, pj_Real?Positive, s_Real?Positive, 
  pis_Real?Positive, c_Real?Positive] = 
 D[profit[pi, pj, s, pis, c], pi]
pisol[s_Real?Positive] := 
 pi /. FindRoot[{pfoc[pi, pj, s, pis, ci] == 0, 
    pfoc[pj, pi, s, pjs, cj] == 0, pi == pis, 
    pj == pjs}, {{pi, 0.3}, {pj, 0.4}, {pis, 0.35}, {pjs, 0.45}}]
pjsol[s_Real?Positive] := 
 pj /. FindRoot[{pfoc[pi, pj, s, pis, ci] == 0, 
    pfoc[pj, pi, s, pjs, cj] == 0, pi == pis, 
    pj == pjs}, {{pi, 0.3}, {pj, 0.4}, {pis, 0.35}, {pjs, 0.45}}]
pisol[0.25]
Plot[{pisol[s], pjsol[s]}, {s, 0., 0.5}]

Integrate was very slow on the Max of linear functions. The result should be piecewise quadratic and can be done by hand. NIntegrate gives errors about non-numerical boundaries even when replacing the symbolic boundary by a number:

Clear[d, s, vi, vj, vlb, mui, muj, pi, pj, profit, ci, cj, pfoc, \
pisol, pjsol, pis, pjs]
$Assumptions = 
  0 < {s, pi, pj, pis, pjs} < 1/2 && 0 < {mui, muj, vlb, vi, vj} < 1 &&
    0 <= {ci, cj} < 1;
vlb[pi_, pj_] = 
 1 + pi - pj - 
  Sqrt[2] Sqrt[s]; ci = 0.; cj = 0.05; mui = 0.5; muj = 0.5; 
d[pi_Real?Positive, pj_Real?Positive, s_Real?Positive, 
  pis_Real?Positive] = 
 mui*(1 - Max[pi, vlb[pi, pj]]) + 
  mui*NIntegrate[
    vi - pi + pj, {vi, Max[pi, 0], Max[pi, vlb[pi, pj]]}] + 
  muj*NIntegrate[(1 - Max[pi, vj - pj + pi]), {vj, 0, vlb[pj, pis]}];
d[0.2, 0.3, 0.1, 0.25]

yields vi = Max[0.,pi] is not a valid limit of integration.

I have encountered similar problems many times before in Mathematica when trying to find fixed points. I gave up each time after many days of trying, so I do not anticipate a solution now.

As I see it, there is a conflict between Mathematica constraints: Integrate does not work (for no good reason) or takes hours, but NIntegrate requires numeric arguments. For all parameters, numeric values would be fine, but for the candidate fixed point (pis in the above code), the maximization or derivative-taking has to happen before setting pis==pi.

One way may be to solve using all possible guesses of pis before pis==pi, maybe by FunctionInterpolation, but this again is too slow.

Answer 1

The offending line with the SetDelayed::shape error is:

{xsola, ysola} := {x, y} /. 
  FindRoot[{D[(x - 1)^2] == 0, 
    D[(y - 2)^2] == 0}, {{x, 0.3}, {y, 0.4}}]

The FullForm of the offending line looks like SetDelay[List[..], ReplaceAll[..]]. The two arguments to SetDelayed[] should be lists of the same shape, if the first argument is a list. In this case, the second argument is not a list at all. When the first argument is a list, SetDelayed makes a series of definitions, one for each element in the list. The value for each comes from the corresponding element in the second argument. (It does this recursively down nest lists.) With Set instead of SetDelayed as @J.M. recommended, the ReplaceAll argument will evaluate to a list.

Given this tidbit:

I am trying to use this in a bigger piece of code where = creates errors of non-numeric values.

It sounds like the problem is to delay the call of FindRoot until certain conditions are met. It might be better to refactor the code so that those conditions can be passed in functions calls.

Not knowing more, something like the following might be a quick workaround:

ClearAll[getSola];
getSola[] :=
  {xsola, ysola} = {x, y} /. 
    FindRoot[{D[(x - 1)^2] == 0, 
      D[(y - 2)^2] == 0}, {{x, 0.3}, {y, 0.4}}];
getSola[]
{xsol, ysol, xsola, ysola}
(*
  {1., 2.}
  {1., 2., 1., 2.}
*)

Also, one could hold up evaluation as follows depending on a condition, or give an error if the condition was not satisfied:

ClearAll[getSola];
SetAttributes[getSola, HoldAll];
getSola[vars_, cond_] /; TrueQ@cond :=
  vars = {x, y} /. 
    FindRoot[{D[(x - a)^2] == 0, 
      D[(y - b)^2] == 0}, {{x, 0.3}, {y, 0.4}}];

getSola[{xsola, ysola}, VectorQ[{a, b}, NumericQ]]

If a and b are numeric, FindRoot is called and the variables are set; otherwise, nothing is done. getSol[] is similar to ConditionalExpression[], except that ConditionalExpression[] evaluates to Undefined when the condition is False and getSola[] remains unevaluated.

Probably the first alternative is better. At least it fits better with how I imagine the bigger code is organized. But I could easily be wrong.

Answer 2

Long comment about an important side-issue in one of the examples: In one case, the assumptions are defined wrong: a < {x, y, z} < b does not translate into three inequalities. One can refactor with VectorLess or Thread:

Clear[d, s, vi, vj, vlb, mui, muj, pi, pj, profit, ci, cj, pfoc, \
pisol, pjsol, pis, pjs]
($Assumptions = 
   Flatten@{Thread[0 < {s, pi, pj, pis, pjs} < 1/2], 
     Thread[0 < {mui, muj, vlb, vi, vj} < 1], 
     Thread[0 <= {ci, cj} < 1]};
  vlb[pi_, pj_] = 1 + pi - pj - Sqrt[2] Sqrt[s]; ci = 0.; cj = 0.05; 
  mui = 0.5; muj = 0.5;
  d[pi_Real?Positive, pj_Real?Positive, s_Real?Positive, 
    pis_Real?Positive] = 
   mui*(1 - Max[pi, vlb[pi, pj]]) + 
    mui*Integrate[
      vi - pi + pj, {vi, Max[pi, 0], Max[pi, vlb[pi, pj]]}] + 
    muj*Integrate[(1 - Max[pi, vj - pj + pi]), {vj, 0, 
       vlb[pj, pis]}];
  d[0.2, 0.3, 0.1, 0.25]) // AbsoluteTiming

(* {0.162937, 0.545695}  *)

Note that exact (symbolic) solvers sometimes run into trouble with inexact floating-point parameters like ci = 0.1. Here it does not, so we're lucky. However, if it did, then one could do the following (and it's faster by a factor of 1/3000 if you have to call d[] multiple times): If you define d[] before you define the parameters, then the integration will be done with exact symbols just as quickly; and d[] returns the same result. Then you don't have to worry about round-off error causing Integrate to choke.