minimizing a function with an array input

Question

I want to use scipy.optimize.minimize on a function with an array type variable input. This is what I hope to do.

I have the following signal,

import numpy as np
time = np.linspace(0, 1, 501)
data = np.cos(2 * np.pi * 4 * time) + np.cos(2 * np.pi * 9 * time) + np.cos(2 * np.pi * 20 * time)
noise = np.sqrt(1 / 25) * np.random.randn(501)
signal = data + noise 

and I am hoping to find a curve fit for this signal. Since I created the data myself, I know that a sum of cosine functions will work for this. So the function that I hope to optimize is the following:

def cos_sum(x, P):

    assert isinstance(P, np.ndarray)
    assert P.shape[0] == P.shape[1]

    sums = []
    for param in P:
        a, b, c = param
        sums.append(a * np.cos(b * (x - c)))
    sums = np.array(sums)

    return np.sum(sums, axis=0) 

In order to use minimize to find the correct parameters, I create this residual function.

def resid(params, x):
    assert isinstance(params, np.ndarray)
    fit = cos_sum(x, params)
    residual = np.sqrt(np.mean(np.abs(fit - signal)) ** 2)
    return residual

I now need to create a guess. Since I already know how the signal was created, I made the following guess:

guess_A = np.random.normal(1, .2, size=3)
guess_B = 2 * np.pi * np.array([4, 9, 20], dtype=float)
guess_C = np.random.normal(0, .2, size=3)
guess = np.array([guess_A, guess_B, guess_C]).T

However, I am unable to run the following:

from scipy.optimize import minimize
optimization = minimize(resid, guess, args=(time))

and I receive the following error:

Traceback (most recent call last):
  File "/Users/nickeisenberg/GitRepos/Python_Misc/Misc/minimize_curvefit_vector_variables.py", line 70, in <module>
    optimization = minimize(resid, guess, args=(time))
  File "/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/scipy/optimize/_minimize.py", line 676, in minimize
    res = _minimize_bfgs(fun, x0, args, jac, callback, **options)
  File "/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/scipy/optimize/_optimize.py", line 1296, in _minimize_bfgs
    sf = _prepare_scalar_function(fun, x0, jac, args=args, epsilon=eps,
  File "/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/scipy/optimize/_optimize.py", line 263, in _prepare_scalar_function
    sf = ScalarFunction(fun, x0, args, grad, hess,
  File "/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py", line 158, in __init__
    self._update_fun()
  File "/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py", line 251, in _update_fun
    self._update_fun_impl()
  File "/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py", line 155, in update_fun
    self.f = fun_wrapped(self.x)
  File "/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py", line 137, in fun_wrapped
    fx = fun(np.copy(x), *args)
  File "/Users/nickeisenberg/GitRepos/Python_Misc/Misc/minimize_curvefit_vector_variables.py", line 53, in resid
    fit = cos_sum(x, params)
  File "/Users/nickeisenberg/GitRepos/Python_Misc/Misc/minimize_curvefit_vector_variables.py", line 30, in cos_sum
    assert P.shape[0] == P.shape[1]
IndexError: tuple index out of range

Is this possible to do?

Answer 1

The problem is that minimize treats parameters as 1D array. Adding print just before failing assert shows, that it reshaped guess array to 1D array.

Changing cos_sum so that it accepts 1D array of params fixes this problem. Here is the code.

def cos_sum(x, P):
    assert isinstance(P, np.ndarray)

    sums = []
    for i in range(0, P.shape[0], 3):
        a, b, c = P[i], P[i+1], P[i+2]
        sums.append(a * np.cos(b * (x - c)))
    sums = np.array(sums)

    return np.sum(sums, axis=0) 

Result:

optimization = minimize(resid, guess, args=(time))
print(optimization.x)

[  0.96805816  25.1679919    0.25020317   1.00261543  56.44511497
   0.77872223   1.00966167 125.71622787   0.55004217]

plt.figure(figsize=(8, 4))
plt.plot(data)
plt.plot([cos_sum(t, optimization.x) for t in time], 'C1--')
plt.show()