3
$\begingroup$

I have an error function $f : w \rightarrow f(w)$ that I want to minimize, $w$ being a vector of length 211. There are some constraints on $w$.

I managed to compute the jacobian $J$ and even with it the optimization is quite slow. I think this is due to this : when I print the iteration number in the optimization loop and the value of the error function $f$, the Python shell shows this :

number of iteration / value of the error function
1   0.035415193666127553
214 0.3451666379304012
215 0.021196928080386743
428 0.23868912897157107
429 0.015584337890888374
642 0.12928511710763613
643 0.010336203822412956
856 0.1488892149935437
857 0.007432756773027516
1070 0.14502939575869633
1071 0.005535691799374814
... and so on

It seems like the algorithm took 213 iterations to estimate the jacobian and compute the next $w$ at each computation of the error function, whereas since I gave it the jacobian of the error function I expected it to directly be able to compute in one iteration only (with the formula of gradient descent for example). I guess this slows down the algorithm a lot since the constraints will be tested on a lot more vector weights $w$.

To be more precise, I expected this :

number of iteration / value of the error function
1 0.035415193666127553
2 0.3451666379304012
3 0.021196928080386743
4 0.23868912897157107
5 0.015584337890888374
6 0.12928511710763613
7 0.010336203822412956
8 0.1488892149935437
9 0.007432756773027516
10 0.14502939575869633
11 0.005535691799374814
...
 and so on

I am using the minimize function from scipy library. My dataset data has the following dataframe format :

Name of the company | Sector of the company | Weight
       APPLE              TECHNOLOGIES         0.04
      WALMART            CONSUMER GOODS        0.06
     MICROSOFT            TECHNOLOGIES         0.03
       ...                     ...             ...

And the constraints are the following : $$ \sum_{i=1}^{211}w_{i} = 1,$$ $$ \forall i, w_{i} \geq 0,$$ $$ \forall ~sector, ~w_{sector}^{min} \leq w_{sector} \leq w_{sector}^{max}.$$ with $w_{sector} = \sum_{i \in sector}w_{i}$. This quantity is obtained via the following function :

def get_weight_sectors(w):
    weight_sectors=pd.Series(index=list_of_all_sectors)
    for sector in weight_sectors.index:
        weight_sector=data['Weight'].loc[data['Sector of the company']==sector].sum()
        weight_sectors.loc[sector]=weight_sector
    return weight_sectors

So my constraints in the minimize function are :

def constraint_sum_weights(w,k): #k is the number of the iteration
    k+=1                         #that allows to keep track within the loop
    return w.sum()-1.0

def long_only_constraint(w):
    return w

def constraint_sector_min(w):
    weight_sectors=get_weight_sectors(w)
    return weight_sectors-weight_sector_min

def constraint_sector_min(w):
    weight_sectors=get_weight_sectors(w)
    return weight_sector_max-weight_sectors

The function to minimize is :

def function_to_minimize(w,k): #k is the number of the iteration
    #computation of the function
    f = ...
    print(k,f)                 #that is what was shown on the Python shell earlier
    return f

And finally the solving function is :

def find_weight():
    k=0
    init=...
    cons=...
    res=minimize(function_to_minimize,init,args=(k,),method='SLSQP',jac=jacobian,constraints=cons,tol=1E-4)

$\endgroup$
10
  • $\begingroup$ Can you turn this into a minimum, working example? $\endgroup$ – Richard Jul 1 '20 at 14:40
  • $\begingroup$ Thank you for your answer Richard, I edited the problem $\endgroup$ – FredNgu Jul 1 '20 at 14:56
  • $\begingroup$ You still have only a summary of the code. To reproduce your issue, a minimum working example of the code with an appropriate input dataset is necessary. $\endgroup$ – Richard Jul 1 '20 at 14:59
  • $\begingroup$ Oh okay sorry I misunderstood the question. Actually I think that I just observed something that would transform my issue to a general issue, without the need of a reproducible example. I am editing the problem. $\endgroup$ – FredNgu Jul 1 '20 at 15:06
  • 1
    $\begingroup$ You don't indicate if you're written your own optimization routine, or if you're using one of Pythons. If the former, you need to show code. If the latter, you need to clearly identify what library function you're using and how you're calling it. The issue may also depend on your dataset. Either way, a MWE is the best way to get help. $\endgroup$ – Richard Jul 1 '20 at 16:01
2
$\begingroup$

Edit: In light of information provided in the comments by @Gabriel Gouvine , I suggest you forego use of scipy,minimize and SLSQP. Instead, try using Pyomo, which cam call a variety of better, more modern solvers, including all solvers which AMPL and GAMS are capable of calling.

Original answer:

Spoiler alert:

scipy.minimize allow specification of linear constraints and bound constraints. So you should be using that capability ... and leave your Jacobian worries behind.

The constraints are linear. Therefore, the Jacobian of the constraints is the matrix of coefficients in the linear system of constraints.

If the modeling package, solver combination allows it (which any decent such should). linear constraints can be treated specially (as opposed to being treated as general nonlinear constrains), and even more special treatment provided for bound constraints (which are a special case of linear constraints).

Therefore, any decent solver, provided the linear constraints as such, can trivially compute the Jacobian itself, and there is no value in having Jacobian provided as input to the solver.

scipy.minimize allow specification of linear constraints and bound constraints. So you should be using that capability ... and leave your Jacobian worries behind.

$\endgroup$
14
  • $\begingroup$ I agree 100% in general, but AFAIK this is not the case of scipy.minimize, except for the trust-constr method $\endgroup$ – Gabriel Gouvine Jul 4 '20 at 14:08
  • $\begingroup$ @ Gabriel Gouvine According to docs.scipy.org/doc/scipy/reference/generated/… , linear and bound constraints can be entered as such for SLSQP. I have never used this package so my answer relies on the correctness of the documentation. $\endgroup$ – Mark L. Stone Jul 4 '20 at 14:12
  • $\begingroup$ No the documentation is misleading. It refers to the mention of trust-constr just before. It's missing a newline. The case for SLSQP and Cobyla is just below. $\endgroup$ – Gabriel Gouvine Jul 4 '20 at 14:14
  • 1
    $\begingroup$ It's just a typo. Not big of a deal for such a large free project. I just made a pull request to fix it. $\endgroup$ – Gabriel Gouvine Jul 4 '20 at 14:27
  • 1
    $\begingroup$ That's changed on GitHub - just added the newline. I don't know when they update the online doc. Maybe only at release time. $\endgroup$ – Gabriel Gouvine Jul 4 '20 at 17:17
3
$\begingroup$

If you check SLSQP's source code, there is actually no code (as far as I can tell) that makes use of the Jacobian passed by the user. Instead, the code assumes no Jacobian is ever passed from the outside and keeps approximating it using the constraints.

What this means is that even if you do provide a Jacobian, the SLSQP code will ignore it, which is why you do not notice a speed difference.

To confirm the suspicion we can check line 525 of the Fortran SLSQP code:

C   CALL JACOBIAN AT CURRENT X

C   UPDATE CHOLESKY-FACTORS OF HESSIAN MATRIX BY MODIFIED BFGS FORMULA

260 DO 270 i=1,n
       u(i) = g(i) - ddot_sl(m,a(1,i),1,r,1) - v(i)
270 CONTINUE

and see that the code to evaluate the Jacobian is not there.

$\endgroup$
4
  • $\begingroup$ This code does not call the Python function. The Jacobian is evaluated by Scipy beforehand if not provided, and this part just uses it to estimate the Hessian. I think passing the Jacobian of the constraints would indeed prevent the multiple evaluations. $\endgroup$ – Gabriel Gouvine Jul 4 '20 at 12:33
  • $\begingroup$ My first link is the Python function, which doesn't use an external Jacobian, nor does it cache it after the first iteration. $\endgroup$ – Nikos Kazazakis Jul 4 '20 at 17:03
  • $\begingroup$ @Gabriel Gouvine if that's not the case can you cite the relevant line in the source code? $\endgroup$ – Nikos Kazazakis Jul 4 '20 at 17:09
  • $\begingroup$ That's the second link of my answer, line 278 $\endgroup$ – Gabriel Gouvine Jul 4 '20 at 17:15
1
$\begingroup$

Scipy's minimize allows you to pass the Jacobian of the constraints to the minimize function (the documentation says it's used only by SLSQP). Then Scipy won't have to evaluate it by calling the function N times.

The Jacobian of the constraints is evaluated here if it's not present.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.