I am using gekko to solve for 14 variables by minimizing around 10,000 equations with IMODE=3.
Each equation is the squared error between a response y and the output of a polynomial model at row i in the training data.
eq[i] = (y[i] - model[i]) ** 2
In each row, the polynomial model has around 10 to 100 terms, where the 14 optimized variables are found. The variables are very interactive in the model, meaning that multiple variables are multiplied together multiple times.
Question: What strategies can I employ to speed up the solving time?
Here is a much simpler reproducible example where the model tries to fit a straight line:
from gekko import GEKKO
import numpy as np
m = GEKKO() # instantiate gekko model
# instantiate free variables
a = m.FV(lb=0, ub=2)
a.STATUS = 1
b = m.FV(lb=0, ub=2)
b.STATUS = 1
c = m.FV(lb=0, ub=2)
c.STATUS = 1
n_eqs1 = 1000 # number of equations in dataset1
n_eqs2 = 500 # number of equations in dataset2
n_terms = 12 # number of terms in each equation
noise_scl = 1 # amount of noise represented as the std of the normal distributions
# training datasets
x = {
"dataset1": np.arange(n_eqs1)[:, np.newaxis]
+ np.random.normal(loc=0, scale=noise_scl, size=(n_eqs1, n_terms)),
"dataset2": np.arange(n_eqs2)[:, np.newaxis]
+ np.random.normal(loc=0, scale=noise_scl, size=(n_eqs2, n_terms)),
}
# response
y = np.arange(n_eqs1)
for x_ds in x.values():
for i in range(x_ds.shape[0]):
# minimize equations
m.Minimize(
(
y[i]
- (
x_ds[i, 0] * a
+ x_ds[i, 1] * a**2
+ x_ds[i, 2] * a * b
+ x_ds[i, 3] * a * (b**2)
+ x_ds[i, 4] * (a**2) * b
+ x_ds[i, 5] * (a**2) * (b**2)
+ x_ds[i, 6] * c
+ x_ds[i, 7] * (c**2)
+ x_ds[i, 8] * c * b
+ x_ds[i, 9] * c * (b**2)
+ x_ds[i, 10] * (c**2) * b
+ x_ds[i, 11] * (c**2) * (b**2)
)
/ n_terms
)
** 2
)
m.options.IMODE = 3
m.solve(disp=True)
# depending on the amount of noise, the optimized values should tend towards 1
print(f"a = {a.value[0]:3f}\n" f"b = {b.value[0]:3f}\n" f"c = {c.value[0]:3f}")
Try using IMODE=2 that is designed for data regression. Total build time will decrease but the solve time should be similar. There isn't much that can be done to decrease the solve time except try better initial guess values, simplify the model, scale the data, remove outliers, or reduce the number of data points.
Total Build time: 0.04206514358520508
Total Solve time: 7.088646650314331
Solver time: 4.3853
a = 0.976090
b = 0.980885
c = 1.048744
Modified for Regression (IMODE=2) with Timing
from gekko import GEKKO
import numpy as np
import time
start = time.time()
m = GEKKO() # instantiate gekko model
# instantiate free variables
a = m.FV(lb=0, ub=2)
a.STATUS = 1
b = m.FV(lb=0, ub=2)
b.STATUS = 1
c = m.FV(lb=0, ub=2)
c.STATUS = 1
n_eqs1 = 1000 # number of equations in dataset1
n_eqs2 = 500 # number of equations in dataset2
n_terms = 12 # number of terms in each equation
noise_scl = 1 # amount of noise represented as the std of the normal distributions
# training datasets
x = {
"dataset1": np.arange(n_eqs1)[:, np.newaxis]
+ np.random.normal(loc=0, scale=noise_scl, size=(n_eqs1, n_terms)),
"dataset2": np.arange(n_eqs2)[:, np.newaxis]
+ np.random.normal(loc=0, scale=noise_scl, size=(n_eqs2, n_terms)),
}
# response
y = np.arange(n_eqs1)
# create combined dataset
xd = np.vstack((x["dataset1"],x["dataset2"]))
yd = np.append(y,y[:n_eqs2])
yp = m.Param(yd)
x_ds = m.Array(m.Param,n_terms)
for i in range(n_terms):
x_ds[i].value = xd[:,i]
# minimize equations
m.Minimize((yp - ( x_ds[0] * a
+ x_ds[1] * a**2
+ x_ds[2] * a * b
+ x_ds[3] * a * (b**2)
+ x_ds[4] * (a**2) * b
+ x_ds[5] * (a**2) * (b**2)
+ x_ds[6] * c
+ x_ds[7] * (c**2)
+ x_ds[8] * c * b
+ x_ds[9] * c * (b**2)
+ x_ds[10] * (c**2) * b
+ x_ds[11] * (c**2) * (b**2)
) / n_terms)**2)
m.options.IMODE = 2
print('Total Build time:',time.time()-start)
start = time.time()
m.solve(disp=False)
print('Total Solve time:',time.time()-start)
print('Solver time:',m.options.SOLVETIME)
# depending on the amount of noise, the optimized
# values should tend towards 1
print(f"a = {a.value[0]:3f}\n" \
f"b = {b.value[0]:3f}\n" \
f"c = {c.value[0]:3f}")
The original model with IMODE=3 solves with 0.32 sec to build the model in Gekko, 6.88 sec once Gekko compiles and solves the problem with 2.92 sec of that time spent by the solver.
Total Build time: 0.3288764953613281
Total Solve time: 6.876295566558838
Solver time: 2.9223
a = 1.021506
b = 0.978350
c = 1.007642
Original (IMODE=3) with Timing
from gekko import GEKKO
import numpy as np
import time
start = time.time()
m = GEKKO() # instantiate gekko model
# instantiate free variables
a = m.FV(lb=0, ub=2)
a.STATUS = 1
b = m.FV(lb=0, ub=2)
b.STATUS = 1
c = m.FV(lb=0, ub=2)
c.STATUS = 1
n_eqs1 = 1000 # number of equations in dataset1
n_eqs2 = 500 # number of equations in dataset2
n_terms = 12 # number of terms in each equation
noise_scl = 1 # amount of noise represented as the std of the nrml distributions
# training datasets
x = {
"dataset1": np.arange(n_eqs1)[:, np.newaxis]
+ np.random.normal(loc=0, scale=noise_scl, \
size=(n_eqs1, n_terms)),
"dataset2": np.arange(n_eqs2)[:, np.newaxis]
+ np.random.normal(loc=0, scale=noise_scl, \
size=(n_eqs2, n_terms)),
}
# response
y = np.arange(n_eqs1)
for x_ds in x.values():
for i in range(x_ds.shape[0]):
# minimize equations
m.Minimize(
(
y[i]
- (
x_ds[i, 0] * a
+ x_ds[i, 1] * a**2
+ x_ds[i, 2] * a * b
+ x_ds[i, 3] * a * (b**2)
+ x_ds[i, 4] * (a**2) * b
+ x_ds[i, 5] * (a**2) * (b**2)
+ x_ds[i, 6] * c
+ x_ds[i, 7] * (c**2)
+ x_ds[i, 8] * c * b
+ x_ds[i, 9] * c * (b**2)
+ x_ds[i, 10] * (c**2) * b
+ x_ds[i, 11] * (c**2) * (b**2)
)
/ n_terms
)
** 2
)
m.options.IMODE = 3
print('Total Build time:',time.time()-start)
start = time.time()
m.solve(disp=False)
print('Total Solve time:',time.time()-start)
print('Solver time:',m.options.SOLVETIME)
# depending on the amount of noise,
# the optimized values should tend towards 1
print(f"a = {a.value[0]:3f}\n" \
f"b = {b.value[0]:3f}\n" \
f"c = {c.value[0]:3f}")