Is there a way to use scipy.optimize.minimize with 2D np array as x0?

Question

I am new to optimization techniques that python offers and I couldn't find a way on how to adjust optimize.minimize from scipy for my needs.

What I would like to do is to minimize the sum of square differences between the market price and the modeled price for the given set of bonds by changing the rate (so the rate would be my x0). This sum is calculated via price_eval, which takes few lists or nested lists as arguments, for instance the nested list of cash flow amounts per every bond. To make the latter more clear - if I have 3 bonds and they have 3, 2 and 4 future cash flows of some amounts like a,b,c; d,e; h;i;j;k respectively, one of the arguments would be the list CF_list = [[a,b,c],[d,e],[h;i;j;kl]].

Here's a function:

def price_eval(rates_list, ttm_list, price_list, CF_list, CF_time_list, N_cf_list):

    price_model_list = []
    sum_squares = 0

    for i in range(0, len(price_list)):
        
        price = 0
        number_cf_per_bond = N_cf_list[i]
        
        for j in range(0, number_cf_per_bond):
            
            price += (CF_list[i][j]) / math.pow(1+rates_list[i][j],CF_time_list[i][j])
   
        diff = price_model_list[i] - price_list[i]
        sum_squares += math.pow(diff,2)

    return sum_squares

I was trying to define the solver as

res = scipy.optimize.minimize(price_eval, rates_guess, args=(ttm_list, price_list, CF_list, CF_time_list, N_cf_list))

where I provide the nested list rates_guess as the initial guess for rates_list[i][j] for solver (it doesn't really matter if I change it to a numpy array or not, gives the same result). The whole problem probably comes from the fact that scipy.optimize.minimize only takes 1-D arrays as x0, while in my case x0 needs to be two-dimensional due to data stored per every bond per every cash flow. I get an error

...
grad[k] = (f(*((xk + d,) + args)) - f0) / d[k]

TypeError: can only concatenate list (not "float") to list

So, my question is whether anyone has an idea on how to tackle this? Or whether it's possible to somehow adjust solver in a way I can provide a nested list/2D np array as x0 which should be given back with new values in a same format? If not, what else can I do?

Thanks a lot for any help!

Answer 1

Both rates_list and rates_guess need to be 1d arrays/lists. The way I usually deal with this is by having functions that flatten your preferred structure into a list, and 'inflate' it back to the desired shape. So you can have, for example, something like this

def to_flat(nested_list):
    '''
    flatten a list of lists
    '''
    flat_list = []
    for l in nested_list:
        for e in l:
            flat_list.append(e)
    return flat_list

# test
rates_nested_list = [[1,1,1],[2,2],[3,3,3,3]]
rates_flat_list = to_flat(rates_nested_list)
print('flat: ',rates_flat_list)

def from_flat(flat_list, list_to_copy_shape_from):
    ''' 
    here the second argument is a nested list in the desired shape -- values not relevant only shape
    '''
    nested_list = []
    idx = 0
    for l in list_to_copy_shape_from:
        nested_list.append([])
        for e in l:
            nested_list[-1].append(flat_list[idx])
            idx += 1
    return nested_list

# test
template = [[0,0,0],[0,0],[0,0,0,0]]
rates_nested_list2 = from_flat(rates_flat_list, template)
print('nested:',rates_nested_list2)

check the output:

flat:  [1, 1, 1, 2, 2, 3, 3, 3, 3]
nested: [[1, 1, 1], [2, 2], [3, 3, 3, 3]]

Now in your price_eval function you will be given, by the optimizer, a flattened rates list and you need to 'inflate' it, so your code would look like this

def price_eval(rates_flat_list, ttm_list, price_list, CF_list, CF_time_list, N_cf_list):

    # create a nested list of rates from a flat one, using CF_list as a template
    rates_flat_list = from_flat(rates_flat_list,CF_list)

    price_model_list = []
    sum_squares = 0

    for i in range(0, len(price_list)):
        
        price = 0
        number_cf_per_bond = N_cf_list[i]
        
        for j in range(0, number_cf_per_bond):
            
            price += (CF_list[i][j]) / math.pow(1+rates_list[i][j],CF_time_list[i][j])

        # !!! I think your code is missing the following line !!!
        price_model_list.append(price)
        diff = price_model_list[i] - price_list[i]
        sum_squares += math.pow(diff,2)

    return sum_squares

Now when you call the optimizer, as x0 you need to pass a flattened list, so assuming rates_guess is a nested list you would call

res = scipy.optimize.minimize(price_eval, x0 = to_flat(rates_guess), args=(ttm_list, price_list, CF_list, CF_time_list, N_cf_list))

finally, when you get res it will have a solution as a flat list so you may want to call from_flat on that again

Finally, since you are doing a least squares fit anyway, you may want to use scipy.optimize.least_squares instead of minimize -- it has pretty much the same inputs but should be better in fitting