I'm currently trying to implement a python script for solving a constrained nonlinear optimization problem with ~800 variables and 2 constraints, one linear and one nonlinear. There already exists a Matlab implementation of this script and I'm trying to use cyipopt to implement an equivalent Python implementation that produces similar results.
The Matlab code has analytic solutions for both the objective and constraint and jacobian, which I've rewritten entirely in Python, and uses L-BFGS for hessian approximation. Calling the Python jacobian functions from fmincon() results in the same solution so I think these functions are correct, and I've also used the derivative checkers in both fmincon() and cyipopt to make sure. However, my Python implementation doesn't converge and fails to match the performance of the Matlab code. In particular, the objective doesn't decrease smoothly. Here's some example output from Ipopt:
This is Ipopt version 3.14.4, running with linear solver MUMPS 5.2.1.
Number of nonzeros in equality constraint Jacobian...: 1602
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 0
Total number of variables............................: 801
variables with only lower bounds: 0
variables with lower and upper bounds: 267
variables with only upper bounds: 0
Total number of equality constraints.................: 2
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.6958872e-01 5.68e-14 1.04e-02 0.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 2.1548097e-01 3.25e+00 7.09e+00 0.6 3.87e+00 - 9.90e-01 4.34e-01f 1
2 1.7668628e-01 3.05e+00 8.15e+00 0.9 2.38e+00 - 1.00e+00 9.72e-01f 1
3 1.7381539e-01 1.78e-01 3.39e+00 0.2 1.09e+00 - 7.37e-01 1.00e+00h 1
4 1.7270930e-01 4.47e-03 2.42e-01 -0.5 2.10e-01 - 1.00e+00 1.00e+00h 1
5 1.7178940e-01 4.69e-03 4.70e-02 -1.2 7.62e-02 - 1.00e+00 1.00e+00f 1
6 1.7055753e-01 5.01e-03 8.89e-03 -1.9 1.21e-01 - 1.00e+00 1.00e+00h 1
7 1.6692096e-01 1.04e-02 1.15e-03 -2.9 7.68e-02 - 1.00e+00 1.00e+00h 1
8 1.6323948e-01 3.18e-01 8.33e-03 -3.3 5.42e-01 - 9.42e-01 1.00e+00h 1
9 1.6005166e-01 1.43e-01 6.64e-03 -9.3 2.50e-01 - 5.61e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 1.5786474e-01 1.22e-01 4.18e-03 -3.7 4.04e-01 - 7.71e-01 1.00e+00h 1
11 1.5672747e-01 6.04e-02 5.11e-03 -4.2 1.61e-01 - 9.61e-01 1.00e+00h 1
12 1.5472463e-01 2.40e-02 2.34e-03 -4.8 1.23e-01 - 9.96e-01 1.00e+00h 1
13 1.5533721e-01 9.66e-04 6.81e-03 -5.0 1.84e-01 - 9.99e-01 1.00e+00H 1
14 1.5377085e-01 3.02e-02 3.08e-03 -5.3 1.65e-01 - 1.00e+00 1.00e+00h 1
15 1.5361137e-01 2.09e-02 4.51e-03 -5.7 1.11e-01 - 1.00e+00 1.00e+00h 1
16 1.5296208e-01 6.95e-03 1.95e-03 -6.0 8.23e-02 - 1.00e+00 1.00e+00h 1
17 1.5273048e-01 1.90e-03 1.01e-03 -7.0 3.83e-02 - 1.00e+00 1.00e+00h 1
18 1.5311695e-01 4.91e-05 5.20e-03 -8.0 2.13e-01 - 1.00e+00 1.00e+00H 1
19 1.5262106e-01 7.92e-03 6.41e-03 -5.9 4.57e-01 - 1.00e+00 2.30e-01h 2
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 1.5261389e-01 1.19e-02 1.27e-02 -3.9 4.90e-01 - 1.00e+00 1.02e-01h 2
21 1.5310130e-01 4.88e-02 4.09e-03 -3.6 2.90e-01 - 4.10e-01 1.00e+00h 1
22 1.5330300e-01 2.75e-02 3.72e-03 -3.9 1.29e-01 - 1.00e+00 1.00e+00h 1
23 1.5201654e-01 1.23e-02 1.40e-03 -3.9 9.69e-02 - 1.00e+00 1.00e+00h 1
24 1.5227903e-01 3.92e-04 6.45e-03 -4.8 2.11e-01 - 9.86e-01 1.00e+00H 1
25 1.5495353e-01 1.86e-04 1.26e-02 -4.5 3.70e-01 - 1.00e+00 1.00e+00H 1
26 1.5557422e-01 3.93e-05 4.07e-03 -4.6 3.53e-01 - 1.18e-01 1.00e+00H 1
27 1.5281467e-01 7.97e-02 4.35e-03 -4.3 3.09e-01 - 1.00e+00 1.00e+00h 1
28 1.6236141e-01 3.83e-02 1.24e-02 -4.4 4.15e-01 - 7.25e-01 8.80e-01H 1
29 1.5299068e-01 8.52e-04 6.17e-03 -4.4 2.48e-01 - 7.17e-02 1.00e+00F 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30 1.5139425e-01 2.53e-02 3.94e-03 -4.4 9.96e-01 - 2.38e-01 1.55e-01h 3
31 1.5138493e-01 5.50e-04 3.76e-03 -4.8 2.40e-01 - 9.94e-01 1.00e+00H 1
32 1.5200809e-01 2.14e-04 5.16e-03 -4.9 1.32e-01 - 9.96e-01 1.00e+00H 1
33 1.5074554e-01 1.86e-02 1.91e-03 -5.0 1.77e-01 - 9.98e-01 1.00e+00h 1
34 1.5130813e-01 1.49e-04 3.53e-03 -5.6 1.56e-01 - 1.00e+00 1.00e+00H 1
35 1.5036412e-01 1.28e-02 1.63e-03 -5.7 1.18e-01 - 1.00e+00 1.00e+00h 1
36 1.5026330e-01 6.65e-04 8.92e-04 -7.2 3.49e-02 - 1.00e+00 1.00e+00h 1
37 1.5076495e-01 1.34e-04 5.36e-03 -7.7 1.20e-01 - 1.00e+00 1.00e+00H 1
38 1.5066012e-01 1.69e-02 6.61e-03 -7.9 1.11e-01 - 1.00e+00 1.00e+00h 1
39 1.5005302e-01 1.25e-02 1.81e-03 -8.0 1.05e-01 - 5.44e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40 1.5001109e-01 3.87e-03 2.08e-03 -8.6 4.98e-02 - 1.00e+00 1.00e+00h 1
41 1.4987818e-01 2.98e-03 2.39e-03 -9.7 5.14e-02 - 1.00e+00 1.00e+00h 1
42 1.5014976e-01 2.70e-05 3.82e-03 -9.8 1.02e-01 - 1.00e+00 1.00e+00H 1
43 1.4997585e-01 1.01e-02 5.16e-03 -10.4 6.66e-02 - 1.00e+00 1.00e+00h 1
44 1.4970903e-01 9.51e-03 2.36e-03 -10.5 8.90e-02 - 8.37e-01 1.00e+00h 1
45 1.4940598e-01 2.37e-03 9.31e-04 -11.0 4.20e-02 - 1.00e+00 1.00e+00h 1
46 1.4968497e-01 2.68e-05 4.63e-03 -11.0 1.21e-01 - 1.00e+00 1.00e+00H 1
47 1.4995823e-01 1.29e-06 7.01e-03 -11.0 3.32e-01 - 1.00e+00 1.00e+00H 1
48 1.4974993e-01 1.46e-02 6.86e-03 -9.0 2.03e+00 - 1.00e+00 1.18e-01f 2
49 1.4948599e-01 1.72e-02 5.95e-03 -9.1 7.57e-01 - 2.53e-01 8.53e-02h 3
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 1.5186145e-01 6.77e-04 8.93e-03 -9.1 2.87e-01 - 7.46e-01 1.00e+00H 1
51 1.5120058e-01 1.35e-03 4.47e-03 -9.1 3.78e-01 - 1.19e-01 1.00e+00H 1
52 1.5550082e-01 1.54e-04 8.87e-03 -9.1 3.78e-01 - 5.42e-01 1.00e+00H 1
53 1.5020757e-01 1.12e-04 4.29e-03 -9.1 4.29e-01 - 1.00e+00 1.00e+00H 1
54 1.4918169e-01 1.28e-02 1.97e-03 -9.1 1.68e+00 - 2.17e-01 6.23e-02h 3
55 1.5133875e-01 4.49e-05 8.71e-03 -9.5 2.28e-01 - 9.87e-01 1.00e+00H 1
56 1.4882545e-01 2.34e-02 2.51e-03 -9.6 1.72e-01 - 4.77e-01 1.00e+00h 1
57 1.4946850e-01 1.48e-02 6.52e-03 -10.4 1.07e-01 - 1.00e+00 1.00e+00h 1
58 1.4859822e-01 5.15e-03 2.33e-03 -10.5 6.48e-02 - 9.27e-01 1.00e+00h 1
59 1.4902635e-01 7.60e-05 4.53e-03 -11.0 1.06e-01 - 1.00e+00 1.00e+00H 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
60 1.4843429e-01 4.65e-03 1.51e-03 -11.0 7.41e-02 - 1.00e+00 1.00e+00h 1
61 1.4839694e-01 1.37e-03 2.49e-03 -11.0 3.13e-02 - 1.00e+00 1.00e+00h 1
62 1.4830507e-01 3.51e-04 8.59e-04 -11.0 2.56e-02 - 1.00e+00 1.00e+00h 1
63 1.4849382e-01 1.36e-05 5.56e-03 -11.0 1.22e-01 - 1.00e+00 1.00e+00H 1
64 1.4823266e-01 3.86e-03 1.66e-03 -11.0 9.19e-02 - 1.00e+00 1.00e+00h 1
65 1.4821480e-01 5.60e-04 1.20e-03 -11.0 1.53e-02 - 1.00e+00 1.00e+00h 1
66 1.4818681e-01 1.50e-04 2.13e-04 -11.0 1.25e-02 - 1.00e+00 1.00e+00h 1
67 1.4816045e-01 1.46e-04 5.02e-04 -11.0 1.61e-02 - 1.00e+00 1.00e+00h 1
68 1.4850582e-01 5.85e-05 5.29e-03 -11.0 1.57e-01 - 1.00e+00 1.00e+00H 1
69 1.5001793e-01 1.58e-06 1.08e-02 -11.0 1.41e-01 - 1.00e+00 1.00e+00H 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
70 1.4861627e-01 2.68e-07 5.43e-03 -11.0 2.13e-01 - 2.10e-01 1.00e+00F 1
71 1.4842174e-01 1.21e-02 4.02e-03 -10.6 3.61e-01 - 1.00e+00 2.50e-01f 3
72 1.4816338e-01 4.71e-03 2.55e-03 -10.7 6.93e-02 - 1.00e+00 1.00e+00h 1
73 1.4841277e-01 1.20e-04 5.26e-03 -10.7 1.25e-01 - 1.00e+00 1.00e+00H 1
74 1.4923205e-01 4.65e-05 7.02e-03 -11.0 9.58e-02 - 1.00e+00 1.00e+00H 1
75 1.4793006e-01 1.50e-02 1.22e-03 -11.0 1.02e-01 - 9.33e-01 1.00e+00h 1
76 1.4790107e-01 1.98e-03 2.30e-03 -11.0 4.36e-02 - 1.00e+00 1.00e+00h 1
77 1.4795403e-01 6.21e-06 2.74e-03 -11.0 3.90e-02 - 1.00e+00 1.00e+00H 1
78 1.4783354e-01 3.00e-03 1.85e-03 -11.0 5.85e-02 - 1.00e+00 1.00e+00h 1
79 1.4775646e-01 1.00e-03 1.66e-03 -11.0 2.36e-02 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
80 1.4772077e-01 2.46e-04 8.16e-04 -11.0 1.37e-02 - 1.00e+00 1.00e+00h 1
81 1.4772337e-01 1.36e-06 1.67e-03 -11.0 2.97e-02 - 1.00e+00 1.00e+00H 1
82 1.4800983e-01 3.43e-07 4.69e-03 -11.0 1.73e-01 - 1.00e+00 1.00e+00H 1
83 1.4785970e-01 3.85e-03 5.48e-03 -10.2 1.53e+00 - 1.00e+00 7.61e-02f 3
84 1.4778198e-01 6.41e-03 6.58e-03 -10.3 3.42e+00 - 6.76e-02 2.92e-02h 3
85 1.4773098e-01 1.02e-02 5.56e-03 -10.3 2.22e+00 - 6.49e-02 4.72e-02h 3
86 1.5775339e-01 1.78e-04 7.17e-03 -10.3 1.27e+00 - 1.35e-01 1.00e+00H 1
87 1.5716495e-01 1.63e-03 5.41e-03 -10.3 1.70e+00 - 1.70e-01 3.17e-02h 4
88 1.5247575e-01 1.33e-04 4.44e-03 -10.3 3.35e-01 - 5.46e-01 1.00e+00H 1
89 1.5188072e-01 1.15e-03 3.96e-03 -10.3 1.80e+00 - 7.48e-01 3.87e-02h 5
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
90 1.5109270e-01 2.99e-03 3.54e-03 -10.3 1.65e+00 - 4.22e-01 5.10e-02h 4
91 1.4947946e-01 5.89e-02 4.20e-03 -10.3 1.50e+01 - 4.46e-02 2.31e-02h 2
92 1.4898332e-01 5.71e-02 3.70e-03 -10.3 6.05e-01 - 3.30e-01 1.30e-01h 3
93 1.4893398e-01 4.01e-02 6.60e-03 -10.3 1.73e-01 - 1.00e+00 1.00e+00h 1
94 1.4813809e-01 2.83e-02 2.88e-03 -10.3 2.64e-01 - 1.00e+00 1.00e+00h 1
95 1.4953181e-01 5.14e-04 7.13e-03 -10.3 1.61e-01 - 1.00e+00 1.00e+00H 1
96 1.4757049e-01 2.19e-02 1.86e-03 -10.3 1.12e-01 - 1.00e+00 1.00e+00h 1
97 1.4838902e-01 1.46e-02 5.65e-03 -10.8 1.05e-01 - 1.00e+00 1.00e+00h 1
98 1.4742559e-01 8.28e-03 5.50e-04 -10.9 7.16e-02 - 1.00e+00 1.00e+00h 1
99 1.4739592e-01 1.63e-04 7.52e-04 -11.0 1.75e-02 - 1.00e+00 1.00e+00h 1
Compared to fmincon():
First-order Norm of
Iter F-count f(x) Feasibility optimality step
0 1 1.695878e-01 5.684e-14 1.102e-02
1 2 1.691645e-01 1.296e-04 4.729e-03 2.229e-02
2 3 1.665545e-01 7.175e-03 9.554e-03 1.857e-01
3 4 1.635913e-01 4.118e-02 6.511e-03 5.518e-01
4 5 1.626062e-01 3.235e-02 1.080e-02 6.491e-01
5 6 1.617645e-01 2.861e-02 3.397e-03 7.418e-01
6 7 1.607960e-01 1.297e-02 2.636e-03 6.037e-01
7 8 1.595868e-01 4.010e-02 4.019e-03 1.089e+00
8 9 1.589008e-01 2.304e-02 4.341e-03 7.246e-01
9 10 1.581489e-01 3.965e-03 1.667e-03 2.636e-01
10 11 1.575713e-01 1.046e-02 1.000e-03 6.139e-01
11 12 1.573578e-01 3.968e-03 1.179e-03 4.697e-01
12 13 1.570832e-01 4.880e-03 1.247e-03 6.278e-01
13 14 1.567794e-01 5.258e-03 1.000e-03 6.424e-01
14 15 1.567468e-01 7.916e-03 2.665e-03 6.705e-01
15 16 1.565624e-01 1.995e-03 1.452e-03 1.778e-01
16 17 1.565058e-01 9.510e-04 1.309e-03 1.582e-01
17 18 1.564559e-01 8.874e-04 1.000e-03 2.408e-01
18 19 1.556438e-01 2.632e-03 7.432e-04 2.010e-01
19 21 1.542575e-01 1.763e-05 1.908e-03 6.447e-01
20 23 1.538128e-01 8.862e-06 2.772e-03 3.140e-01
21 24 1.535477e-01 1.020e-03 1.317e-03 1.241e-01
22 25 1.531595e-01 2.294e-03 7.816e-04 2.747e-01
23 26 1.529536e-01 5.902e-04 8.447e-04 1.609e-01
24 28 1.524555e-01 4.137e-07 1.424e-03 4.367e-01
25 30 1.521461e-01 2.148e-06 1.428e-03 2.705e-01
26 31 1.520693e-01 3.572e-04 9.764e-04 7.670e-02
27 33 1.518785e-01 1.542e-06 1.821e-03 3.113e-01
28 34 1.516523e-01 7.795e-04 1.052e-03 1.813e-01
29 35 1.515454e-01 3.468e-04 6.490e-04 1.431e-01
30 36 1.514033e-01 6.813e-04 7.628e-04 2.445e-01
First-order Norm of
Iter F-count f(x) Feasibility optimality step
31 37 1.511966e-01 1.321e-03 1.065e-03 3.590e-01
32 38 1.509549e-01 2.310e-03 9.453e-04 5.227e-01
33 39 1.507060e-01 2.390e-03 1.893e-03 5.238e-01
34 40 1.505274e-01 1.191e-03 4.680e-04 2.007e-01
35 41 1.504018e-01 3.995e-04 7.486e-04 1.168e-01
36 42 1.501441e-01 8.662e-04 5.371e-04 3.376e-01
37 44 1.498647e-01 3.364e-06 1.694e-03 5.798e-01
38 45 1.496842e-01 9.911e-04 7.416e-04 3.829e-01
39 46 1.496415e-01 3.685e-04 5.740e-04 1.727e-01
40 47 1.495934e-01 9.745e-04 1.368e-03 3.101e-01
41 48 1.494941e-01 6.501e-04 6.743e-04 3.020e-01
42 49 1.493509e-01 9.861e-04 3.953e-04 3.974e-01
43 50 1.492265e-01 7.372e-04 6.141e-04 3.550e-01
44 51 1.490722e-01 1.345e-03 7.354e-04 4.686e-01
45 52 1.489949e-01 9.143e-04 5.345e-04 2.969e-01
46 53 1.489130e-01 8.080e-04 7.258e-04 2.955e-01
47 54 1.488745e-01 1.176e-03 1.445e-03 2.992e-01
48 55 1.488323e-01 4.312e-04 5.706e-04 1.628e-01
49 56 1.487758e-01 3.440e-04 8.374e-04 1.993e-01
50 57 1.486999e-01 1.017e-03 6.874e-04 3.923e-01
51 59 1.486946e-01 8.701e-06 3.005e-03 4.956e-01
52 60 1.485959e-01 3.722e-04 9.674e-04 7.810e-02
53 61 1.485395e-01 3.153e-04 5.199e-04 1.871e-01
54 62 1.484944e-01 3.883e-04 9.311e-04 2.667e-01
55 63 1.484078e-01 6.886e-04 1.208e-03 4.148e-01
56 64 1.482546e-01 2.396e-03 1.258e-03 9.236e-01
57 65 1.481228e-01 2.799e-03 1.541e-03 9.661e-01
58 66 1.479890e-01 7.611e-04 5.626e-04 4.157e-01
59 67 1.478692e-01 1.020e-03 9.555e-04 5.353e-01
60 68 1.477201e-01 1.205e-03 2.336e-03 6.695e-01
First-order Norm of
Iter F-count f(x) Feasibility optimality step
61 69 1.474560e-01 3.556e-03 1.588e-03 1.439e+00
62 70 1.471856e-01 2.773e-03 1.294e-03 1.377e+00
63 71 1.469532e-01 1.536e-03 1.162e-03 8.180e-01
64 72 1.467492e-01 1.390e-03 8.390e-04 7.719e-01
65 73 1.464826e-01 2.369e-03 2.728e-03 1.328e+00
66 74 1.461599e-01 2.225e-03 1.697e-03 1.276e+00
67 75 1.458142e-01 2.839e-03 7.377e-04 1.729e+00
68 76 1.455273e-01 1.941e-03 1.540e-03 1.555e+00
69 77 1.451421e-01 2.727e-03 2.188e-03 2.072e+00
70 78 1.446933e-01 6.136e-03 3.128e-03 2.527e+00
71 79 1.442774e-01 4.200e-03 1.005e-03 1.977e+00
72 80 1.440018e-01 1.392e-03 1.141e-03 4.495e-01
73 81 1.436448e-01 3.002e-03 2.341e-03 1.335e+00
74 82 1.432150e-01 4.996e-03 3.013e-03 1.999e+00
75 83 1.427922e-01 3.273e-03 7.698e-04 1.622e+00
76 84 1.424515e-01 1.755e-03 1.547e-03 1.999e+00
77 85 1.420710e-01 1.985e-03 1.608e-03 2.162e+00
78 86 1.414360e-01 5.327e-03 1.684e-03 3.757e+00
79 88 1.410962e-01 8.909e-05 5.970e-03 7.093e+00
80 89 1.405116e-01 3.316e-03 1.177e-03 8.055e-01
81 90 1.402913e-01 1.138e-03 1.549e-03 3.040e-01
82 91 1.398874e-01 3.111e-03 2.327e-03 1.897e+00
83 93 1.391672e-01 3.531e-06 3.701e-03 5.891e+00
84 95 1.383758e-01 5.642e-05 3.542e-03 6.218e+00
85 96 1.378290e-01 4.166e-03 1.036e-03 2.995e+00
86 97 1.375162e-01 1.780e-03 1.041e-03 1.597e+00
87 98 1.370376e-01 5.208e-03 1.925e-03 3.465e+00
88 99 1.366982e-01 7.099e-03 4.238e-03 3.143e+00
89 100 1.362886e-01 2.651e-03 2.194e-03 2.074e+00
90 101 1.360006e-01 2.295e-03 1.189e-03 1.023e+00
First-order Norm of
Iter F-count f(x) Feasibility optimality step
91 102 1.358196e-01 9.869e-04 1.742e-03 1.240e+00
92 103 1.355212e-01 2.761e-03 1.672e-03 1.891e+00
93 105 1.353705e-01 1.162e-04 4.479e-03 4.375e+00
94 106 1.348992e-01 2.477e-03 1.632e-03 3.747e-01
95 107 1.347110e-01 1.267e-03 1.224e-03 3.268e-01
96 108 1.345712e-01 8.739e-04 1.756e-03 8.004e-01
97 109 1.343092e-01 3.165e-03 2.251e-03 1.691e+00
98 111 1.341921e-01 4.984e-03 2.419e-03 1.037e+00
99 112 1.339379e-01 1.967e-03 9.087e-04 2.052e+00
I don't really have any experience with either the theory or practice of optimization, so I'm wondering if anyone can help spot the issue here? Is there any way I can tweak the options I'm passing to IPOPT to solve the problem or is IPOPT simply less suited to this type of optimization problem?
For reference, here's the python code:
class initial_problem():
def __init__(self, cell_pairs, dC):
self.cell_pairs = cell_pairs
self.dC = dC
def objective(self, x):
x = x.reshape(x0.shape, order='F')
q = x[:,0:2]
p = x[:,2]
b = np.matmul(self.dC, np.multiply(q.T,p).T)
delta_p = np.matmul(self.dC, p)
t1 = np.divide((b - (np.multiply(r1.T,delta_p).T)).T,np.linalg.norm(b - (np.multiply(r1.T,delta_p).T), axis=1)).T
t2 = np.divide((b - (np.multiply(r2.T,delta_p).T)).T,np.linalg.norm(b - (np.multiply(r2.T,delta_p).T), axis=1)).T
E = 0.5 * np.mean(np.power(np.sum(t1 * tau_1, axis=1), 2) +
np.power(np.sum(t2 * tau_2, axis=1), 2))
return E
def gradient(self, x):
# Jacobian of objective
x = x.reshape(x0.shape, order='F')
q = x[:,0:2]
p = x[:,2]
b = np.matmul(self.dC, np.multiply(q.T,p).T)
delta_p = np.matmul(self.dC, p)
t1 = np.divide((b - (np.multiply(r1.T,delta_p).T)).T,np.linalg.norm(b - (np.multiply(r1.T,delta_p).T), axis=1)).T
t2 = np.divide((b - (np.multiply(r2.T,delta_p).T)).T,np.linalg.norm(b - (np.multiply(r2.T,delta_p).T), axis=1)).T
ip1 = np.sum(t1 * tau_1, axis=1)
ip2 = np.sum(t2 * tau_2, axis=1)
drX1 = sx_grad(p[self.cell_pairs[:,0]], p[self.cell_pairs[:,1]], q[self.cell_pairs[:,0],0], q[self.cell_pairs[:,0],1], q[self.cell_pairs[:,1],0], q[self.cell_pairs[:,1],1], r1[:,0], r1[:,1])
drX2 = sx_grad(p[self.cell_pairs[:,0]], p[self.cell_pairs[:,1]], q[self.cell_pairs[:,0],0], q[self.cell_pairs[:,0],1], q[self.cell_pairs[:,1],0], q[self.cell_pairs[:,1],1], r2[:,0], r2[:,1])
drY1 = sy_grad(p[self.cell_pairs[:,0]], p[self.cell_pairs[:,1]], q[self.cell_pairs[:,0],0], q[self.cell_pairs[:,0],1], q[self.cell_pairs[:,1],0], q[self.cell_pairs[:,1],1], r1[:,0], r1[:,1])
drY2 = sy_grad(p[self.cell_pairs[:,0]], p[self.cell_pairs[:,1]], q[self.cell_pairs[:,0],0], q[self.cell_pairs[:,0],1], q[self.cell_pairs[:,1],0], q[self.cell_pairs[:,1],1], r2[:,0], r2[:,1])
dE = np.multiply(ip1 * tau_1[:,0], drX1.T).T + np.multiply(ip1 * tau_1[:,1], drY1.T).T + \
np.multiply(ip2 * tau_2[:,0], drX2.T).T + np.multiply(ip2 * tau_2[:,1], drY2.T).T
dE = dE / self.dC.shape[0]
rows = np.concatenate([self.cell_pairs[:,0],self.cell_pairs[:,0]+self.dC.shape[1],self.cell_pairs[:,0]+2*self.dC.shape[1],
self.cell_pairs[:,1],self.cell_pairs[:,1]+self.dC.shape[1],self.cell_pairs[:,1]+2*self.dC.shape[1]])
dE = np.bincount(rows, weights=np.ravel(dE,order='F'))
return dE
# Define constraints
def constraints(self, x):
# linear constraint
Aeq = np.concatenate((np.zeros(x0.shape[0]*2), np.ones(x0.shape[0])))
lincon = np.dot(Aeq, x)
# nonlinear constraint
x = x.reshape(x0.shape, order='F')
q = x[:,0:2]
p = x[:,2]
b = np.matmul(self.dC, np.multiply(q.T,p).T)
delta_p = np.matmul(self.dC, p)
l1 = np.linalg.norm(b - (np.multiply(r1.T,delta_p).T), axis=1)
l2 = np.linalg.norm(b - (np.multiply(r2.T,delta_p).T), axis=1)
nonlincon = 0.5*(np.mean(l1) + np.mean(l2)) - scale
cons = np.append(nonlincon, lincon)
return cons
def jacobian(self, x):
# Jacobian of nonlinear constraints
x = x.reshape(x0.shape, order='F')
q = x[:,0:2]
p = x[:,2]
b = np.matmul(self.dC, np.multiply(q.T,p).T)
delta_p = np.matmul(self.dC, p)
t1 = np.divide((b - (np.multiply(r1.T,delta_p).T)).T,np.linalg.norm(b - (np.multiply(r1.T,delta_p).T), axis=1)).T
t2 = np.divide((b - (np.multiply(r2.T,delta_p).T)).T,np.linalg.norm(b - (np.multiply(r2.T,delta_p).T), axis=1)).T
drX1 = rx_grad(p[self.cell_pairs[:,0]], p[self.cell_pairs[:,1]], q[self.cell_pairs[:,0],0], q[self.cell_pairs[:,1],0], r1[:,0])
drX2 = rx_grad(p[self.cell_pairs[:,0]], p[self.cell_pairs[:,1]], q[self.cell_pairs[:,0],0], q[self.cell_pairs[:,1],0], r2[:,0])
drY1 = ry_grad(p[self.cell_pairs[:,0]], p[self.cell_pairs[:,1]], q[self.cell_pairs[:,0],1], q[self.cell_pairs[:,1],1], r1[:,1])
drY2 = ry_grad(p[self.cell_pairs[:,0]], p[self.cell_pairs[:,1]], q[self.cell_pairs[:,0],1], q[self.cell_pairs[:,1],1], r2[:,1])
dE = 0.5 * (np.multiply(t1[:,0],drX1.T).T + np.multiply(t1[:,1],drY1.T).T) + \
0.5 * (np.multiply(t2[:,0],drX2.T).T + np.multiply(t2[:,1],drY2.T).T)
rows = np.concatenate([self.cell_pairs[:,0],self.cell_pairs[:,0]+self.dC.shape[1],self.cell_pairs[:,0]+2*self.dC.shape[1],
self.cell_pairs[:,1],self.cell_pairs[:,1]+self.dC.shape[1],self.cell_pairs[:,1]+2*self.dC.shape[1]])
dE = np.bincount(rows, weights=np.ravel(dE/self.dC.shape[0],order='F'))
# Add jacobian of linear constraints
lin_dE = np.concatenate((np.zeros(2*x0.shape[0]), np.ones(x0.shape[0])))
jac = np.concatenate([dE, lin_dE])
return jac
initial_minimization = cyipopt.Problem(
n = x0.size,
m = 2,
problem_obj = initial_problem(self.cell_pairs,self.dC),
lb = np.concatenate((-np.inf*np.ones(x0.shape[0]), -np.inf*np.ones(x0.shape[0]), 0.001*np.ones(x0.shape[0]))),
ub = np.concatenate((np.inf*np.ones(x0.shape[0]), np.inf*np.ones(x0.shape[0]), 1000*np.ones(x0.shape[0]))),
cl = np.array([0, np.mean(p0)*p0.shape[0]]),
cu = np.array([0, np.mean(p0)*p0.shape[0]])
)
initial_minimization.add_option('max_iter',2000)
initial_minimization.add_option('tol',1e-5)
initial_minimization.add_option('nlp_scaling_method','none')
initial_minimization.add_option('print_timing_statistics','yes')
initial_minimization.add_option('hessian_approximation','limited-memory')
res = initial_minimization.solve(x0.ravel(order='F'))
Compared to the original matlab code:
energyFunc = @(x) reducedEnergy( x, sparse(d0), t1, t2, r1, r2, bCells );
nonlinFunc = @(x) mycon( x, sparse(d0), Lscale, r1, r2, bCells );
lb = -inf*ones(size(x0));
lb(:,3) = .001;
ub = inf*ones(size(x0));
ub(:,3) = 1000;
Aeq = [zeros(1,size(d0,2)),zeros(1,size(d0,2)),ones(1,size(d0,2))];
beq = size(d0,2)*mean(p);
optimset = optimoptions('fmincon','Display','iter-detailed','Algorithm','interior-point','MaxFunEvals',2e6,...
'MaxIter',2e3,'TolFun',1e-5,'GradObj','on','GradConstr','on',...
'HessianApproximation','lbfgs');
[x,ERes] = fmincon(energyFunc,x0,[],[],Aeq,beq,lb,ub,nonlinFunc,optimset);
function [ clin, ceq, dlin, deq ] = mycon( x, d0, Lscale, r1, r2, bCells )
NB = size(d0,1);
NC = size(d0,2);
q = x(:,1:2);
p = x(:,3);
clin = [];
B = d0*bsxfun(@times,q,p);
deltaP = d0*p;
B1 = B - bsxfun(@times,deltaP,r1);
B2 = B - bsxfun(@times,deltaP,r2);
L1 = sqrt(sum(B1.^2,2));
B1 = bsxfun(@rdivide,B1,L1);
L2 = sqrt(sum(B2.^2,2));
B2 = bsxfun(@rdivide,B2,L2);
ceq = .5*(mean(L1) + mean(L2)) - Lscale;
if (nargout > 2)
dlin = [];
deq = .5*(bsxfun(@times,B1(:,1),drX1) + bsxfun(@times,B1(:,2),drY1)) + ...
.5*(bsxfun(@times,B2(:,1),drX2) + bsxfun(@times,B2(:,2),drY2));
rows = [bCells(:,1);bCells(:,1)+NC;bCells(:,1)+2*NC;bCells(:,2);bCells(:,2)+NC;bCells(:,2)+2*NC];
deq = accumarray(rows,deq(:)/NB,[3*NC,1]);
end
end
function [ E, dE ] = reducedEnergy( x, d0, t1, t2, r1, r2, bCells )
q = x(:,1:2);
p = x(:,3);
B = d0*bsxfun(@times,q,p);
deltaP = d0*p;
B1 = B - bsxfun(@times,deltaP,r1);
B2 = B - bsxfun(@times,deltaP,r2);
L1 = sqrt(sum(B1.^2,2));
B1 = bsxfun(@rdivide,B1,L1);
L2 = sqrt(sum(B2.^2,2));
B2 = bsxfun(@rdivide,B2,L2);
E = .5*mean( dot(B1,t1,2).^2 + dot(B2,t2,2).^2 );
if (nargout == 2)
NB = size(d0,1);
NC = size(d0,2);
[ drX1, drY1, drX2, drY2 ] = fitDual.AFN.reducedLocalGrads( q, p, bCells, r1 ,r2 );
IP1 = dot(B1,t1,2);
IP2 = dot(B2,t2,2);
dE = bsxfun(@times,IP1.*t1(:,1),drX1) + bsxfun(@times,IP1.*t1(:,2),drY1) + ...
bsxfun(@times,IP2.*t2(:,1),drX2) + bsxfun(@times,IP2.*t2(:,2),drY2);
dE = dE / NB;
rows = [bCells(:,1);bCells(:,1)+NC;bCells(:,1)+2*NC;bCells(:,2);bCells(:,2)+NC;bCells(:,2)+2*NC];
vals = dE(:);
dE = accumarray(rows,vals',[3*NC,1]);
end
end
