Scipy: Find root of function with vector input and scalar output

Question

Here is my minimal working example. It's quite simple, I basically am given a density value z and I would like to find any point on that contour. I do this by root finding.

import numpy as np
from scipy.stats import multivariate_normal
from scipy.optimize import root

# Define the density
mu = np.zeros(2)
Sigma = np.eye(2)
target = multivariate_normal(mu, Sigma)

# Contour value
z = target.pdf(target.rvs())

# Find a point on its contour
def func(x):
    return target.pdf(x) - z
root(func, mu)

This throws the following error

TypeError: fsolve: there is a mismatch between the input and output shape of the 'func' argument 'func'.Shape should be (2,) but it is (1,).

Answer 1

The documentation of root says that, fun is

A vector function to find a root of.

I guess they mean by that the fun is a mapping from nDim -> nDim and not from nDim -> 1.

So one option is two blow up your function artificially to fulfill this requirement.

def func(x):
    y = target.pdf(x) - z
    return [y] * len(x)

Another option is to use minimize instead of root.

import numpy as np
from scipy.stats import multivariate_normal
from scipy.optimize import root, minimize

target = multivariate_normal(np.zeros(2), np.eye(2))
z = target.pdf(target.rvs())


def fun_root(x):
    return [target.pdf(x) - z]*2


def fun_minimize(x):
    return (target.pdf(x) - z)**2


x0 = np.random.uniform(low=-1, high=+1, size=2)

res_root = root(fun=fun_root, x0=x0)
print('root')
print(res_root.x)
print(res_root.fun)

print('minimize')
res_min = minimize(fun=fun_minimize, x0=x0, tol=1e-10)
print(res_min.x)
print(res_min.fun)

In addition, I would suggest using another initial guess. And it seems to me that there are a lot of possible solutions for your example, so be aware of that.