title: Scipy interp1d inversion when using an array input: how to find the root of the output function array?
body: <p>Let's say I have the following 3d-Numpy-array A:</p>
  
  <pre><code>array([[[ 1,  2,  3],
          [ 1,  4,  9]],
  
         [[ 1,  8, 27],
          [ 1, 16, 81]]])
  </code></pre>
  
  <p>I want to interpolate the data using <code>interp1d</code> with <code>axis=2</code>, i.e. interpolate the function values <code>1,2,3</code>, <code>1,4,9</code>, <code>1,8,27</code> and <code>1,16,81</code> with <code>x=np.array([1,2,3])</code> simultaneously (the values represent f(x)=x, x^2, x^3 and x^4, respectively).</p>
  
  <p>Luckily, this is possible since interp1d is able to take ND-arrays as a second argument, as long as the interpolation axis is equal to the length of x:</p>
  
  <pre><code>f = interp1d(x,A,kind='linear',axis=2,fill_value='extrapolate')
  </code></pre>
  
  <p>So far this works pretty well, f(1) produces</p>
  
  <pre><code>array([[1., 1.],
         [1., 1.]]),
  </code></pre>
  
  <p>f(2) produces</p>
  
  <pre><code>array([[ 2.,  4.],
         [ 8., 16.]])
  </code></pre>
  
  <p>and the interpolation f(1.5) produces</p>
  
  <pre><code>array([[1.5, 2.5],
         [4.5, 8.5]])
  </code></pre>
  
  <p>i.e. I get a 2x2 array of interpolation functions f(x), evaluated at x.</p>
  
  <p>Now comes the problem: I want to invert these interpolation functions to get an an x for a specific function value - element wise for each entry of the 2x2 array. </p>
  
  <p>When dealing with a 1D-function g this is usually done by finding the root of the interpolation function minus the requested function value, say a:</p>
  
  <pre><code>g_subtracted = lambda x: f(x) - a
  </code></pre>
  
  <p>and finding the zero, for example using scipy.optimize.newton:</p>
  
  <pre><code>optimize.newton(g_subtracted,1.0)
  </code></pre>
  
  <p>where 1.0 is an initial guess. </p>
  
  <p>Here is my actual question: <strong>how can I find the zero of my 2x2 interp2d function array f element wise?</strong> When I simply do</p>
  
  <pre><code>f_subtracted = lambda x: f(x) - a
  </code></pre>
  
  <p>(this works well for my function array) and then</p>
  
  <pre><code>optimize.newton(f_subtracted,1.0)
  </code></pre>
  
  <p>I get the following error:</p>
  
  <pre><code>File "&lt;ipython-input-187-4cf34581a978&gt;", line 1, in &lt;module&gt;
  newton(f,1.0)
  
  File "/anaconda3/lib/python3.6/site-packages/scipy/optimize/zeros.py", line 201, in newton
  if q1 == q0:
  
  ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
  </code></pre>
  
  <p>Seems to me that optimize.newton is not able to deal with 'function arrays' as input to evaluate them element-wise. Does anyone know a way how to do it?</p>
  
  <p>Thanks a lot in advance!</p>