2
\$\begingroup\$

I am trying to optimize writing a script that has 2 minimizations, one dependent on the other. My code is a bit bloated, and I find the global parameters that I am solving for depend quite strongly on the initial guesses for the local parameters (which is quite concerning). So while the code works, I'm looking to optimize it, both for speed and robustness (i.e. the global solutions shouldn't vary so much due to local initial guesses).

Let's say you have multiple systems of linear equations:

#system of equations 1
equation1=dF*pF+dO*pO+dC*pC
equation2=dF*pF2+dO*pO2+dC*pC2
equation3=dF*pF3+dO*pO3+dC*pC3
equation4=dF*pF4+dO*pO4+dC*pC4
#system 2
equation5=dF2*pF+dO2*pO+dC2*pC
equation6=dF2*pF2+dO2*pO2+dC2*pC2
equation7=dF2*pF3+dO2*pO3+dC2*pC3
equation8=dF2*pF4+dO2*pO4+dC2*pC4
...

For each system of equations, there are 4 equations, with 15 adjustable parameters (3 of them shared dF, dO, and dC). However, the other 12 are shared across the multiple various systems. There are given experimental values of "equation1, equation2, equation3...." that the above system of equations is minimized against.

These 12 parameters can be determined from 4 global parameters (k,k1,k2,k3) [io is a constant that is given]

pF=(sqrt(k)*sqrt(k1)*sqrt(kk1+8io(1+k1)-kk1)/(4(1+k1))
pF2=(sqrt(k*k2)*sqrt(k1*k2)*sqrt(kk2k1k2+8io(1+k1*k2)-kk2k1k2)/(4(1+k1k2))
pF3=(sqrt(k*k3)*sqrt(k1*k3)*sqrt(kk3k1k3+8io(1+k1k3)-kk3k1k3)/(4(1+k1*k3))
pF4=(sqrt(k*k2*k3)*sqrt(k1*k2*k3)*sqrt(kk2k3k1k2k3+8io(1+k1k2k3)-kk2k3k1k2k3)/(4(1+k1*k2*k3))
#all other ps can be calculated in the same manner, but they have different equations I won't list here for space reasons

So the idea is the 4 global parameters (k,k1,k2,k3) can be used to calculate the p's. Then with given p's, you go through each system of equation individually and solve for dF, dO, dC. You then minimize the global parameters, using the combined chi2 of the "solved/minimized" local parameters (dF, dO, dC).

Here is the code:

from scipy.optimize import minimize
import numpy as np

experimental_data_list=[[117.77, 117.705, 117.843, 117.597], [110.575, 110.258, 110.167, 110.216], [125.691, 125.006, 125.327, 124.481], [107.491, 108.461, 107.804, 109.383], [128.689, 128.383, 128.668, 128.29], [125.969, 126.326, 126.28, 126.257], [122.439, 122.684, 122.859, 122.194], [125.989, 125.998, 125.985, 125.897], [120.916, 120.18, 120.345, 120.567], [126.772, 126.669, 127.006, 127.592], [120.176, 120.153, 119.864, 120.205]]

def local_calculation(d,pF,pO,pC,pF2,pO2,pC2,pF3,pO3,pC3,pF4,pO4,pC4,experimental_data):
    equation1=(d[0]*pF)+(d[1]*pO)+(d[2]*pC)
    equation2=(d[0]*pF2)+(d[1]*pO2)+(d[2]*pC2)
    equation3=(d[0]*pF3)+(d[1]*pO3)+(d[2]*pC3)
    equation4=(d[0]*pF4)+(d[1]*pO4)+(d[2]*pC4)
    return np.sqrt(np.sum((experimental_data-np.array([equation1,equation2,equation3,equation4]))**2))

def get_populations(k,io):
    pF=(((np.sqrt(k[0])*np.sqrt(k[1]))*(np.sqrt((8*io*(k[1]+1))+(k[0]*k[1])))-(k[0]*k[1]))/(4*(k[1]+1)))/io
    pO=((k[1]*((-np.sqrt(k[0])*np.sqrt(k[1]))*(np.sqrt((8*io*(k[1]+1))+(k[0]*k[1])))+(k[0]*k[1])+(4*io*(k[1]+1))))/(4*((k[1]+1)**2)))/io
    pC=(((-np.sqrt(k[0])*np.sqrt(k[1]))*(np.sqrt((8*io*(k[1]+1))+(k[0]*k[1])))+(k[0]*k[1])+(4*io*(k[1]+1)))/(4*((k[1]+1)**2)))/io
    pF2=(((np.sqrt((k[0]*k[2]))*np.sqrt((k[1]*k[2])))*(np.sqrt((8*io*((k[1]*k[2])+1))+((k[0]*k[2])*(k[1]*k[2]))))-((k[0]*k[2])*(k[1]*k[2])))/(4*((k[1]*k[2])+1)))/io
    pO2=(((k[1]*k[2])*((-np.sqrt((k[0]*k[2]))*np.sqrt((k[1]*k[2])))*(np.sqrt((8*io*((k[1]*k[2])+1))+((k[0]*k[2])*(k[1]*k[2]))))+((k[0]*k[2])*(k[1]*k[2]))+(4*io*((k[1]*k[2])+1))))/(4*(((k[1]*k[2])+1)**2)))/io
    pC2=(((-np.sqrt((k[0]*k[2]))*np.sqrt((k[1]*k[2])))*(np.sqrt((8*io*((k[1]*k[2])+1))+((k[0]*k[2])*(k[1]*k[2]))))+((k[0]*k[2])*(k[1]*k[2]))+(4*io*((k[1]*k[2])+1)))/(4*(((k[1]*k[2])+1)**2)))/io
    pF3=(((np.sqrt((k[0]*k[3]))*np.sqrt((k[1]*k[3])))*(np.sqrt((8*io*((k[1]*k[3])+1))+((k[0]*k[3])*(k[1]*k[3]))))-((k[0]*k[3])*(k[1]*k[3])))/(4*((k[1]*k[3])+1)))/io
    pO3=(((k[1]*k[3])*((-np.sqrt((k[0]*k[3]))*np.sqrt((k[1]*k[3])))*(np.sqrt((8*io*((k[1]*k[3])+1))+((k[0]*k[3])*(k[1]*k[3]))))+((k[0]*k[3])*(k[1]*k[3]))+(4*io*((k[1]*k[3])+1))))/(4*(((k[1]*k[3])+1)**2)))/io
    pC3=(((-np.sqrt((k[0]*k[3]))*np.sqrt((k[1]*k[3])))*(np.sqrt((8*io*((k[1]*k[3])+1))+((k[0]*k[3])*(k[1]*k[3]))))+((k[0]*k[3])*(k[1]*k[3]))+(4*io*((k[1]*k[3])+1)))/(4*(((k[1]*k[3])+1)**2)))/io
    pF4=(((np.sqrt((k[0]*k[2]*k[3]))*np.sqrt((k[1]*k[2]*k[3])))*(np.sqrt((8*io*((k[1]*k[2]*k[3])+1))+((k[0]*k[2]*k[3])*(k[1]*k[2]*k[3]))))-((k[0]*k[2]*k[3])*(k[1]*k[2]*k[3])))/(4*((k[1]*k[2]*k[3])+1)))/io
    pO4=(((k[1]*k[2]*k[3])*((-np.sqrt((k[0]*k[2]*k[3]))*np.sqrt((k[1]*k[2]*k[3])))*(np.sqrt((8*io*((k[1]*k[2]*k[3])+1))+((k[0]*k[2]*k[3])*(k[1]*k[2]*k[3]))))+((k[0]*k[2]*k[3])*(k[1]*k[2]*k[3]))+(4*io*((k[1]*k[2]*k[3])+1))))/(4*(((k[1]*k[2]*k[3])+1)**2)))/io
    pC4=(((-np.sqrt((k[0]*k[2]*k[3]))*np.sqrt((k[1]*k[2]*k[3])))*(np.sqrt((8*io*((k[1]*k[2]*k[3])+1))+((k[0]*k[2]*k[3])*(k[1]*k[2]*k[3]))))+((k[0]*k[2]*k[3])*(k[1]*k[2]*k[3]))+(4*io*((k[1]*k[2]*k[3])+1)))/(4*(((k[1]*k[2]*k[3])+1)**2)))/io
    local_chi2=[]
    for experimental_data in experimental_data_list:
        arguments=(pF,pO,pC,pF2,pO2,pC2,pF3,pO3,pC3,pF4,pO4,pC4,experimental_data)
        local_solution=minimize(local_calculation,args=arguments, x0=np.array([120,110,100]),method = 'Nelder-Mead')
        local_chi2.append(local_solution.fun)
    return sum(local_chi2)


io=280000
global_parameter_solution=minimize(get_populations,args=io, x0=np.array([1000,0.002,8,20]),method = 'Nelder-Mead')
Improve this question
\$\endgroup\$
6
  • \$\begingroup\$ Cite your reference(s) please. The meanings of the experimental observations, and of the global_parameter_solution.x vector, are unclear. I observe either (A.) extremely poor numeric stability, or (B.) valid local minima randomly sampled from a rather mountainous solution surface, not sure which interpretation applies. It would be very nice to work on a simpler subproblem which yields solutions in less than a minute. Are you happy with the condition number of the equations being solved? \$\endgroup\$
    – J_H
    Mar 30 at 19:21
  • \$\begingroup\$ Hello J_H, there is no reference, this is a model I generated myself, so I'm afraid there is no "simpler subproblem" which yields a good solution quickly. I too have observed a mountainous solution surface, albeit I've only looked at the global paramaters (k,k1,k2,k3), I haven't looked at the dependency of the global parameters on fluctuations of the local ones (dF, dO, dC). That being said, as long as the logic flow follows, and the code does represent the logic flow, then the results (as poor as they may be) are what they are (i.e. I would accept them). \$\endgroup\$
    – samman
    Mar 30 at 20:17
  • \$\begingroup\$ I just wanted to add, the main purpose of my post is to see if there a cleaner way to do this type of "dual minimization". But if anything is unclear, please let me know and I can further explain the logic flow and what I am trying to do. \$\endgroup\$
    – samman
    Mar 30 at 20:20
  • \$\begingroup\$ Please do not edit the question, especially the code, after an answer has been posted. Changing the question may cause answer invalidation. Everyone needs to be able to see what the reviewer was referring to. What to do after the question has been answered. If you want a review of updated code please post a follow up question with a link back to this question. \$\endgroup\$
    – pacmaninbw
    Mar 31 at 18:18
  • \$\begingroup\$ @pacmaninbw isn't that redundant though? I have the original code here, so users can see what was approved, and the answer where I applied the improvements to. If I make a follow up post, they'd need to refer back to this answer (plus I don't want to make 10 posts with follow up improvements/answers if I made a mistake in my corrections/improvements). \$\endgroup\$
    – samman
    Mar 31 at 18:31

1 Answer 1

Reset to default
4
\$\begingroup\$

As a start:

Do not write experimental_data_list all in one line; give each row its own line and put it into an np.array.

Combine pF0,1,2,3 into one row of a 3x4 matrix, to also include O and C.

Remove your redundant parens. There are so, so many. The inner expressions in get_populations are formatted in a truly awful way and I am not going to suggest any further simplifications there until that's cleaned up - for a small example,

((np.sqrt(k[0])*np.sqrt(k[1]))

is really just

np.sqrt(k0*k1)

with k0, k1 unpacked from k.

Do not make local_chi2 a list; make it a float starting at 0 and maintain a running total.

Do not call an inner minimize. Recognize that your local_calculation is really a matrix multiplication of d @ pFOC. Ask numpy for a least-squares matrix inversion and use the residual it returns for your chi.

Suggested

Please scrutinize this closely. Add numeric regression tests to make sure that get_populations still works in the same way and step through a few times in a debugger.

from scipy.optimize import minimize
import numpy as np

experimental_data_list = np.array((
    [117.770, 117.705, 117.843, 117.597],
    [110.575, 110.258, 110.167, 110.216],
    [125.691, 125.006, 125.327, 124.481],
    [107.491, 108.461, 107.804, 109.383],
    [128.689, 128.383, 128.668, 128.290],
    [125.969, 126.326, 126.280, 126.257],
    [122.439, 122.684, 122.859, 122.194],
    [125.989, 125.998, 125.985, 125.897],
    [120.916, 120.180, 120.345, 120.567],
    [126.772, 126.669, 127.006, 127.592],
    [120.176, 120.153, 119.864, 120.205],
))


def local_calculation(
    d: np.ndarray,
    pFOC: np.ndarray,
    experimental_data: np.ndarray,
):
    equations = d @ pFOC
    return np.linalg.norm(experimental_data - equations)


def get_populations(k: np.ndarray, io: float) -> float:
    k0, k1, k2, k3 = k

    pF = (
        (((np.sqrt(k0)*np.sqrt(k1))*(np.sqrt((8*io*(k1+1))+(k0*k1)))-(k0*k1))/(4*(k1+1)))/io,
        (((np.sqrt((k0*k2))*np.sqrt((k1*k2)))*(np.sqrt((8*io*((k1*k2)+1))+((k0*k2)*(k1*k2))))-((k0*k2)*(k1*k2)))/(4*((k1*k2)+1)))/io,
        (((np.sqrt((k0*k3))*np.sqrt((k1*k3)))*(np.sqrt((8*io*((k1*k3)+1))+((k0*k3)*(k1*k3))))-((k0*k3)*(k1*k3)))/(4*((k1*k3)+1)))/io,
        (((np.sqrt((k0*k2*k3))*np.sqrt((k1*k2*k3)))*(np.sqrt((8*io*((k1*k2*k3)+1))+((k0*k2*k3)*(k1*k2*k3))))-((k0*k2*k3)*(k1*k2*k3)))/(4*((k1*k2*k3)+1)))/io,
    )

    pO = (
        ((k1*((-np.sqrt(k0)*np.sqrt(k1))*(np.sqrt((8*io*(k1+1))+(k0*k1)))+(k0*k1)+(4*io*(k1+1))))/(4*((k1+1)**2)))/io,
        (((k1*k2)*((-np.sqrt((k0*k2))*np.sqrt((k1*k2)))*(np.sqrt((8*io*((k1*k2)+1))+((k0*k2)*(k1*k2))))+((k0*k2)*(k1*k2))+(4*io*((k1*k2)+1))))/(4*(((k1*k2)+1)**2)))/io,
        (((k1*k3)*((-np.sqrt((k0*k3))*np.sqrt((k1*k3)))*(np.sqrt((8*io*((k1*k3)+1))+((k0*k3)*(k1*k3))))+((k0*k3)*(k1*k3))+(4*io*((k1*k3)+1))))/(4*(((k1*k3)+1)**2)))/io,
        (((k1*k2*k3)*((-np.sqrt((k0*k2*k3))*np.sqrt((k1*k2*k3)))*(np.sqrt((8*io*((k1*k2*k3)+1))+((k0*k2*k3)*(k1*k2*k3))))+((k0*k2*k3)*(k1*k2*k3))+(4*io*((k1*k2*k3)+1))))/(4*(((k1*k2*k3)+1)**2)))/io,
    )

    pC = (
        (((-np.sqrt(k0)*np.sqrt(k1))*(np.sqrt((8*io*(k1+1))+(k0*k1)))+(k0*k1)+(4*io*(k1+1)))/(4*((k1+1)**2)))/io,
        (((-np.sqrt((k0*k2))*np.sqrt((k1*k2)))*(np.sqrt((8*io*((k1*k2)+1))+((k0*k2)*(k1*k2))))+((k0*k2)*(k1*k2))+(4*io*((k1*k2)+1)))/(4*(((k1*k2)+1)**2)))/io,
        (((-np.sqrt((k0*k3))*np.sqrt((k1*k3)))*(np.sqrt((8*io*((k1*k3)+1))+((k0*k3)*(k1*k3))))+((k0*k3)*(k1*k3))+(4*io*((k1*k3)+1)))/(4*(((k1*k3)+1)**2)))/io,
        (((-np.sqrt((k0*k2*k3))*np.sqrt((k1*k2*k3)))*(np.sqrt((8*io*((k1*k2*k3)+1))+((k0*k2*k3)*(k1*k2*k3))))+((k0*k2*k3)*(k1*k2*k3))+(4*io*((k1*k2*k3)+1)))/(4*(((k1*k2*k3)+1)**2)))/io,
    )

    pFOC = np.array((pF, pO, pC))

    local_chi2 = 0

    for experimental_data in experimental_data_list:
        # d @ pFOC ~ experimental_data; what is d?
        d, (residuals,), rank, singular = np.linalg.lstsq(pFOC.T, experimental_data)

        local_chi2 += np.sqrt(residuals)

    return local_chi2


io = 280_000
global_parameter_solution = minimize(
    get_populations,
    args=io,
    x0=(1000, 0.002, 8, 20),
    method='Nelder-Mead',
)

There are other improvements possible, but let's get the basics out of the way first.

Refactor

Pull out common expressions from the pFOC calculation. Further vectorise to a single matrix inverse instead of doing it row-wise. This produces

from typing import Iterable

from scipy.optimize import minimize
import numpy as np

experimental_data_list = np.array((
    [117.770, 117.705, 117.843, 117.597],
    [110.575, 110.258, 110.167, 110.216],
    [125.691, 125.006, 125.327, 124.481],
    [107.491, 108.461, 107.804, 109.383],
    [128.689, 128.383, 128.668, 128.290],
    [125.969, 126.326, 126.280, 126.257],
    [122.439, 122.684, 122.859, 122.194],
    [125.989, 125.998, 125.985, 125.897],
    [120.916, 120.180, 120.345, 120.567],
    [126.772, 126.669, 127.006, 127.592],
    [120.176, 120.153, 119.864, 120.205],
))


def get_pfoc(k: Iterable[float]) -> np.ndarray:
    k0, k1, k2, k3 = k
    ki = np.array((1, k2, k3, k2*k3))
    k01ii = k0*k1*ki*ki
    a = 2*(k1*ki + 1)
    s = np.sqrt(
        k01ii*(4*io*a + k01ii)
    ) - k01ii

    f = s/io/a/2
    c = (2*io*a - s)/io/a/a
    o = k1*ki*c

    return np.vstack((f, o, c))


def get_populations(k: np.ndarray) -> float:
    pFOC = get_pfoc(k)

    # pFOC.T @ d = edl.T
    # 4x3 @ 3x11 = 4x11
    d, *_ = np.linalg.lstsq(pFOC.T, experimental_data_list.T, rcond=None)
    error = (d.T@pFOC - experimental_data_list).ravel()
    return error.dot(error)


io = 280_000

# Regression test for pFOC calculation
assert np.allclose(
    get_pfoc((1000, 0.002, 8, 20)),
    np.array((
        [0.00188615, 0.014887  , 0.03638202, 0.23081748],
        [0.00199224, 0.01551359, 0.03706223, 0.18646849],
        [0.9961216 , 0.96959941, 0.92655575, 0.58271403],
    )),
)

global_parameter_solution = minimize(
    fun=get_populations,
    x0=(1000, 0.002, 8, 20),
    bounds=((0, None),)*4,
)
print(global_parameter_solution)

  message: CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL
  success: True
   status: 0
      fun: 0.23517135040546988
        x: [ 9.999e+02  8.614e-01  4.891e+00  2.441e+00]
      nit: 22
      jac: [ 6.384e-07  1.033e-06  2.839e-06 -9.689e-06]
     nfev: 200
     njev: 40
 hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>
Improve this answer
\$\endgroup\$
5
  • \$\begingroup\$ Hello Reinderien, thank you for your suggestions! The least squared method for overdetermined linear systems is far faster than a default minimization (and is no longer dependent on initial guess issues). I realized there is a lot of reptition in the equations, so I rewrote them, however in doing so I noticed an interesting result. The results of the minimization are now different because the old equation and new equation are identical, but give results with 1e-7 or 1e-8 differences (e.g. 1.00000009 vs. 1.00000008). \$\endgroup\$
    – samman
    Mar 31 at 17:53
  • \$\begingroup\$ I don't quite know which is a better way of doing things? I've noticed that some initial conditions for the global minimization will fail with the old equations, but work with the new equations (again...the equations are identical). But the new code is cleaner and faster now, thank you! \$\endgroup\$
    – samman
    Mar 31 at 17:54
  • \$\begingroup\$ @samman please see edit. lstsq is now called once per iteration. I have also tried this with basinhopping and it really doesn't change, it's a stable result. \$\endgroup\$
    – Reinderien
    Mar 31 at 18:52
  • \$\begingroup\$ 1e-8 differences will be difficult to avoid in practical floating-point math \$\endgroup\$
    – Reinderien
    Mar 31 at 18:54
  • \$\begingroup\$ For posterity, the follow-up chat is here: chat.stackexchange.com/transcript/145060 \$\endgroup\$
    – Reinderien
    Apr 5 at 0:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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