import numpy as np
from sympy import Matrix
from scipy.optimize import minimize

def Testfunction(x):
    a = np.array([[1],[3],[6]])  #3x1 np.ndarray
    b = np.array([[2]])          #1x1 np.ndarray
    c = Matrix(x@x.T).vech()     #3x1 mutable dense matrix
    y = np.transpose(np.subtract(a,c)) * b * np.subtract(a,c)
    return y                     #y should be a 1x1 np.ndarray

x0= np.array([[1],[3]])       # 2x1 np.ndarray
print(Matrix(x0@x0.T).vech()) # works (3x1 mutable dense matrix)
print(Testfunction(x0))       # works (3x3 instead of 1x1 np.ndarray though)

result = minimize(Testfunction, x0, method='nelder-mead', options={'xtol': 1e-8, 'disp': True}) #does not work as explained above