Optimizing a function where one of the parameters is an array

Question

I want to optimize a function by varying the parameters where two of the parameters are actually arrays. I've tried to do

...
# initial parameters
params0 = np.array([p1, p2, ... , p_array1, p_array2])
p_min = minimize(myfunc, params0, args)
...

where the pj's are scalars and p_array1 and p_array2 are arrays of the same length, but this gave me an error saying

ValueError: setting an array element with a sequence.

I've also tried passing p_array1 and p_array2 as scalars into myfunc and then create predetermined arrays from those two inside myfunc (e.g. setting p_array1 = p_array1*np.arange(6) and similarly for p_array2), eliminating the error, but I don't want them to be predetermined -- instead I want 'minimize' to figure out what they should be.

Is there any way that I can utilize one of Scipy's optimization functions without getting this error while still keeping p_array1 and p_array2 as arrays and not scalars?

EDIT

Sorry for being very broad but here is my code:

NOTE: 'myfunc' here is actually norm_residual .

import pandas as pd
import numpy as np

def f(yvec, t, a, b, c, d, M, theta):
    # the system of ODEs to be solved
    x, y = yvec
    dydt = [ a*x - b*y**2 + 1, -c*x - d*x*y + np.sum(M * np.cos(theta*t)) ]
    return dydt

ni = 3 # the number of periodic forcing functions to add to the DE system
M = 0.56*np.random.rand(ni) # the initial amplitudes of forcing functions
theta = np.pi/6*np.arange(ni) # the initial coefficients of the forcing functions

# initialize the parameters
params0 = [0.75, 0.23, 1.0, 0.2, M, theta]

# grabbing the data to be used later
data = pd.read_csv('data.csv')
y_data = data['Y']

N = y_data.shape[0] #20
t = np.linspace(0, N, N) # array of t values to integrate over
yvec0 = [0.3, 0.34] # initial conditions for x and y respectively

def norm_residual(params, *args):
    """
    Computes the L^2 norm of the residual of y and the data (y as defined above).
    Input:    params = array of parameters (scalars or arrays) for the DE system
              args = other arguments to pass into the function f or to use
                   to compute the residual.
    Output: err = L^2 error of the solution vector (scalar).
    """
    data, yvec0, t = args
    a, b, c, d, M, theta = params
    sol = odeint(f, yvec0, t, args=(a, b, c, d, M, theta))
    x = sol[:, 0]; y = sol[:, 1]
    res = data - y
    err = np.linalg.norm(res, 2)
    return err

from scipy.optimize import minimize

p_min = minimize(norm_residual, params0, args=(y_data, yvec0, t))
print(p_min)

And the traceback

Traceback (most recent call last):
  File "model_ex_1.py", line 62, in <module>
    p_min = minimize(norm_residual, params0, args=(y_anom, yvec0, t))
  File "/usr/lib/python2.7/dist-packages/scipy/optimize/_minimize.py", line 354, in minimize
    x0 = np.asarray(x0)
  File "/usr/lib/python2.7/dist-packages/numpy/core/numeric.py", line 482, in asarray
    return array(a, dtype, copy=False, order=order)
ValueError: setting an array element with a sequence.

Whatever you do, you have to use a flattened-view in regards to minimize-API! If there are params: ```p1, p2, p3, p_array1, p_array2``` where ```len(p_array1)``` = N and ```len(p_array2)``` = M, you will want to provide an x0 of size ```3 + N + M``` and you could unpack those in your function (probably first lines) like ```p1, p2, p3 = x[:3]```, ```p_array1 = x[3:3+N]``` and ```p_array2 = x[3+3+N:]``` and do whatever you want to do on those mixed types. ```params0``` could probably be created by ```np.hstack((p1, p2, p3, p_array1, p_array2))``` or in some similar way (under some assumptions). — sascha, May 21 '18 at 15:27
@YakymPirozhenko I have edited my post to include more of the code and traceback. See the note above. — BRz, May 21 '18 at 20:26
Why was my post down-voted? I have put in some research effort as indicated in my trial and errors as shown in my post. If it was due to clarity or usefulness, I've also included the function definition of my system of ODEs. — BRz, May 22 '18 at 13:53

score 0 · Answer 1 · answered May 21 '18 at 15:15

You cannot put a list in a numpy array if the other elements are scalars.

>>> import numpy as np
>>> foo_array = np.array([1,2,3,[5,6,7]])
Traceback (most recent call last):
  File "<pyshell#1>", line 1, in <module>
    foo_array = np.array([1,2,3,[5,6,7]])
ValueError: setting an array element with a sequence.

score 0 · Answer 2 · answered May 21 '18 at 16:21

It would be helpful if you post myfunc but you can do this -

def foo():
    return [p0,p1,p2..pn]

params0 = numpy.array([foo(), p_array1, p_array2])
p_min = minimize(myfunc, params0, args)

OR from Multiple variables in SciPy's optimize.minimize

import scipy.optimize as optimize

def f(params):
    # print(params)  # <-- you'll see that params is a NumPy array
    a, b, c = params # <-- for readability you may wish to assign names to the component variables
    return a**2 + b**2 + c**2

initial_guess = [1, 1, 1]
result = optimize.minimize(f, initial_guess)
if result.success:
    fitted_params = result.x
    print(fitted_params)
else:
    raise ValueError(result.message)

I have edited my post and included more code if it helps. I will try to work on this idea. — BRz, May 21 '18 at 20:23

BRz · Accepted Answer · 2018-05-22T13:36:50.147

I figured it out! The solution that I found to work was to change

params0 = [0.75, 0.23, 1.0, 0.2, M, theta]

in line 6 to

params0 = np.array([ 0.75, 0.23, 1.0, 0.2, *M, *theta], dtype=np.float64)

and in my function definition of my system of ODEs to be solved, instead of having

def f(yvec, t, a, b, c, d, M, theta):
    x, y = yvec
    dydt = [ a*x - b*y**2 + 1, -c*x - d*x*y + np.sum(M * np.cos(theta*t)) ]
    return dydt

I now have

def f(yvec, t, myparams):
    x, y = yvec
    a, b, c, d = myparams[:4]
    ni = (myparams[4:].shape[0])//2 # halved b/c M and theta are of the same shape
    M = myparams[4:ni+4]
    theta = myparams[ni+4:]
    dydt = [ a*x - b*y**2 + 1, -c*x - d*x*y + np.sum(M * np.cos(theta*t)) ]
    return dydt

NOTE: I had to add "dtype=np.float64" for 'params0' because I was getting the error

AttributeError: 'numpy.float64' object has no attribute 'cos'

when I did not have it there and it appears that 'cos' does not know how to handle 'ndarray' objects. The workaround can be found here.

Thanks everyone for the suggestions!

Optimizing a function where one of the parameters is an array

3 Answers3