Scipy optimization with matrix multiplication

Question

I've tried to use spicy.optimize.minimize to solve a matrix multiplication optimization problem, however, the result gives me a dimension error, can someone help me with it?

import numpy as np
from scipy.optimize import minimize

# define known variables, mu, sigma, rf
mu = np.matrix([[0.12], 
                [0.08], 
                [0.05]])

sigma = np.matrix([[0.5, 0.05, 0.03],
                   [0.05, 0.4, 0.01],
                   [0.03, 0.01, 0.2]])

rf = 0.02

def objective_fun(x):
'''
This is the objective function
'''
    s = np.sqrt(x.T * sigma * x)/(mu.T * x - rf)
    return s

def constraint(x):
    con = 1 
    for i in np.arange(0,3):
        con = con - x[i] 
    return con

# set up the boundaries for x
bound_i = (0, np.Inf)
bnds = (bound_i, bound_i, bound_i)

#set up the constraints for x
con = {'type':'eq', 'fun':constraint}

# initial guess for variable x
x = np.matrix([[0.5],
               [0.3],
               [0.2]])

sol = minimize(objective_fun, x, method = 'SLSQP', bounds = bnds, constraints = con)

The error gives me:

ValueError                                Traceback (most recent call last)
<ipython-input-31-b8901077b164> in <module>
----> 1 sol = minimize(objective_fun, x, method = 'SLSQP', bounds = bnds, constraints = con)

e:\Anaconda3\lib\site-packages\scipy\optimize\_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
    606     elif meth == 'slsqp':
    607         return _minimize_slsqp(fun, x0, args, jac, bounds,
--> 608                                constraints, callback=callback, **options)
    609     elif meth == 'trust-constr':
    610         return _minimize_trustregion_constr(fun, x0, args, jac, hess, hessp,

e:\Anaconda3\lib\site-packages\scipy\optimize\slsqp.py in _minimize_slsqp(func, x0, args, jac, bounds, constraints, maxiter, ftol, iprint, disp, eps, callback, **unknown_options)
    397 
    398             # Compute objective function
--> 399             fx = func(x)
    400             try:
    401                 fx = float(np.asarray(fx))

e:\Anaconda3\lib\site-packages\scipy\optimize\optimize.py in function_wrapper(*wrapper_args)
    324     def function_wrapper(*wrapper_args):
    325         ncalls[0] += 1
--> 326         return function(*(wrapper_args + args))
    327 
    328     return ncalls, function_wrapper

<ipython-input-28-b1fb2386a380> in objective_fun(x)
      3     This is the objective function
      4     '''
----> 5     s = np.sqrt(x.T * sigma * x)/(mu.T * x - rf)
      6     return s

e:\Anaconda3\lib\site-packages\numpy\matrixlib\defmatrix.py in __mul__(self, other)
    218         if isinstance(other, (N.ndarray, list, tuple)) :
    219             # This promotes 1-D vectors to row vectors
--> 220             return N.dot(self, asmatrix(other))
    221         if isscalar(other) or not hasattr(other, '__rmul__') :
    222             return N.dot(self, other)

ValueError: shapes (1,3) and (1,3) not aligned: 3 (dim 1) != 1 (dim 0)

However, I tried individually every function I wrote, they all have no errors in the end, like, if after defining the x matrix as shown in the code, I simply run objective_fun(x) in the console, and I immediately get an answer:

optimize_fun(x)
matrix([[5.90897598]])

That means that my function can do the matrix multiplication correctly, so what is wrong with the code here?

Answer 1

The docs for minimize() says that x0 should be an (n,) shaped array, but you are trying to treat it like a (3,1) array. I'm not sure on the inner workings of minimize() but I suspect when it steps over different values of the fit parameters it converts to the format that it thinks it wants. Anyways, the following minor corrections make it so the code works.

import numpy as np
from scipy.optimize import minimize

# define known variables, mu, sigma, rf
mu = np.matrix([[0.12], 
                [0.08], 
                [0.05]])

sigma = np.matrix([[0.5, 0.05, 0.03],
                   [0.05, 0.4, 0.01],
                   [0.03, 0.01, 0.2]])

rf = 0.02

def objective_fun(x):
  '''
  This is the objective function
  '''
  x = np.expand_dims(x, 1) # convert the (3,) shape to (3,1). Then we can do our normal matrix math on it
  s = np.sqrt(x.T * sigma * x)/(mu.T * x - rf) # Transposes so the shapes are correct
  return s

def constraint(x):
  con = 1 
  for i in np.arange(0,3):
      con = con - x[i] 
  return con

# set up the boundaries for x
bound_i = (0, np.Inf)
bnds = (bound_i, bound_i, bound_i)

#set up the constraints for x
con = {'type':'eq', 'fun':constraint}

# initial guess for variable x

x = np.array([0.5, 0.3, 0.2]) # Defining the initial guess as an (3,) array)

sol = minimize(objective_fun, x, method = 'SLSQP', bounds = bnds, constraints = con)
print(sol) # and the solution looks reasonable

Output

     fun: 5.86953830952583
     jac: array([-1.70555401, -1.70578796, -1.70573896])
 message: 'Optimization terminated successfully.'
    nfev: 32
     nit: 6
    njev: 6
  status: 0
 success: True
       x: array([0.42809911, 0.29522438, 0.27667651])

Take a look at the comments I put in for an explanation on what you need to do.