I would like to create an optimization solution for black-box software calculations. Currently, I am using the Levenberg-Marquardt algorithm to update a vector of parameters, $\beta$, with residuals, $r$, determined by outputting the black-box calculation and model parameters to a readable table.

$(J^TJ+\lambda diag(J^TJ))\delta=J^Tr(\beta)$

$\beta_{n}=\beta_{n-1} + \alpha \delta$

$r(\beta) = y_{data} - f_{model}(\beta)$

The Jacobian is estimated using Broyden's method where $\Delta f_n$ and $\Delta x_n$ are the difference in the residual vector and parameter vector between iterations respectively. Since Broyden's method is for root finding, then $\Delta f_n = r_n -r_{n-1} \rightarrow 0 $ should fit the system.

enter image description here

I have little experience in programming, but I think I'm close to having a solution by iterating on the following loop:

deltaf_n = np.reshape(r_n - r_old, (len(r_n), 1))  
deltabeta_n = np.reshape(beta_n - beta_old,  (len(beta_n), 1)) 
broyden_n = broyden_old + np.matmul((deltaf_n - np.matmul(broyden_old,\ 
            deltabeta_n))/np.sum(deltabeta_n**2), deltabeta_n.T)
jTj = np.matmul(broyden_n.T, broyden_n)
delta = np.linalg.solve(jTj + lam*jTj*np.eye(3), -np.matmul(broyden_n.T, r_n))
beta_old = np.copy(beta_n) 
r_old = np.copy(r_n)
broyden_old = np.copy(broyden_n)
beta_n = beta_n + alpha*delta
r_n = residuals(y_data, beta_n)

Unfortunately I'm not getting a solution because the parameters aren't updating as expected: enter image description here

Does anyone have any idea why my parameters aren't getting updated properly? You can see that the quadratic parameter is getting updated, but for some reason the constant parameter is hardly moving.

  • $\begingroup$ Maybe a derivative-free method, such as Nelder-Mead, would be preferable or at least worth trying. en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method $\endgroup$ – Charlie S Feb 11 at 18:43
  • $\begingroup$ @CharlieS, thanks for the reply, that could be an option. I was hoping to avoid having to do the development myself since I'm sure someone out there has already done this. Since responses have been sparse I have started creating a python script which uses Broyden's method to estimate the jacobian for the LM algorithm. It almost works, but isn't updating all the model parameters properly over each iteration. I think my issue is more related to my lack of coding experience than any problem with the method. $\endgroup$ – jborb3663 Feb 11 at 19:08
  • $\begingroup$ Whichever method you end up using (LM may be the best -- who knows?), be sure to scale your variables such that they have similar magnitudes. There is an excellent discussion of it here. alglib.net/optimization/scaling.php $\endgroup$ – Charlie S Feb 11 at 19:31
  • $\begingroup$ I changed the question and added a lot more details to see if someone can spot where I'm making a mistake that prevents the parameters from getting updated properly. $\endgroup$ – jborb3663 Feb 11 at 20:09
  • $\begingroup$ Not an expert on Python, but this looks incorrect: delta = np.linalg.solve(jTj + lamjTjnp.eye(3), -np.matmul(broyden_n.T, r_n)) The matrix should be $J^{T}J+D$ where $D$ is a diagonal matrix (e.g. $(1+\lambda)I$.) However, you seem to to have $J^{T}J+\lambda(J^{T}J)I$, which isn't correct. $\endgroup$ – Brian Borchers Feb 14 at 5:57

The problem ended up being elsewhere in the code with complete working solution posted here:

## Python adaptation of optimization routine as conceptualized by Markus Piro circa 2014.
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

## function to test against
beta_true = np.array([5, 33, -1e4])
def test(x, beta):
    return beta[0]*(x+beta[1])**2 + beta[2]
## function to calculate residuals
def residuals(y_exp, beta):
    return y_exp - test(x_calc, beta)
## The true solution. The thing we are hoping to replicate through LMA
m = 50 #the number of calculated data points
x_calc = np.linspace(-100, 50, num=m)
y_true = test(x_calc, beta_true)

## synthetic experimental data. This will be used in LMA to generate fit
mult_err = np.random.uniform(low=0.8, high=1.2, size=np.shape(x_calc))
add_err = np.random.uniform(low=-1, high=1, size=np.shape(x_calc))
y_exp =  y_true * mult_err + 5e3*add_err

## Start figure where we will add traces
plot, ax1 = plt.subplots()
ax1.plot(x_calc, y_exp, 'b.', label='Simulated Data')
ax1.plot(x_calc, y_true, 'k--', label='Truth', alpha=0.6, linewidth=2.5)
ax1.set_ylim([-15000, 35000])
def trace(y_calc):
    ax1.plot(x_calc, y_calc, 'r-', alpha=0.1)

## begin process outlined in Piro et al.
## 'A Jacobian Free Deterministic Method for Solving Inverse Problems' (unpublished work)
lam = 0.1 #damping parameter
alpha = 0.2 #step length for direction vector
num_its = 100 #how many iterations to do

## initialize parameters to zero
n = np.shape(beta_true)[0] #how many parameters
beta_0 = np.reshape(np.zeros(n), (n, 1)) #initialize guess parameters to zeros
r_0 = np.reshape(residuals(y_exp, beta_0), (len(y_exp), 1)) #calculate initial residuals
beta_1 =beta_0 + 0.1*np.random.rand(n,1)#perturb initial parameter guess to produce beta_1 and calculate r_1
r_1 = np.reshape(residuals(y_exp, beta_1), (len(y_exp), 1))
r_old = np.copy(r_0)
r_k = np.copy(r_1)
beta_old = np.copy(beta_0)
beta_k = np.copy(beta_1)
broyden_old = np.zeros((m,n))
ax1.plot(x_calc, test(x_calc, beta_k), 'r-', alpha=0.3, label='Iteration')

## looping begins here
while j<=num_its:
    t_k = r_k - r_old  #calculate difference in residuals
    s_k = beta_k - beta_old #calculate difference in parameters
    y_calc = test(x_calc, beta_k) #calculate iteration trace for plotting
    trace(y_calc) #add trace to plot

    broyden_k = broyden_old + np.matmul((t_k - \
                np.matmul(broyden_old, s_k))/np.sum(s_k**2), s_k.T) #Approximate the jacobian update
    jTj = np.matmul(broyden_k.T, broyden_k) #convenience variable
    p = np.linalg.solve(jTj + lam*jTj*np.eye(n), -np.matmul(broyden_k.T, r_k)) #compute direction vector
    beta_old = np.copy(beta_k) #set up for iteration
    r_old = np.copy(r_k)
    broyden_old = np.copy(broyden_k)
    beta_k = beta_k + alpha*p
    r_k = np.reshape(residuals(y_exp, beta_k), (len(y_exp), 1))
    #print(np.sum(r_k**2), beta_k)

# show figure with simulated data, true solution, and iteration traces

This code could be adapted to optimize a model produced in black-box software. You would need to add the capability to read tabulated data outputs and figure out some way to return the adjusted parameters back to the black-box software each iteration.enter image description here


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.