Integrate penalty functions into minimization

[Cong Wang]

Dear all,

I wonder if there is a recommended way to integrate penalty functions into minimization using lmfit. 

A picture below summarizes the situation of minimizing a penalized objective.  (Sorry for the terrible handwriting if you find it unrecognizable)

It seems lmfit.minimize with least_squares method only accepts an array of residuals, that is the y_i - yhat_i term in the objective function.  So as a workaround, my implementation penalizes every residual by doing the following.  As a result, it does its job in the most part. 

However, I believe this implementation is nowhere close to a good one, especially from a math perspective (the square operation would lead to three terms in the end).  Therefore, I would really appreciate your opinions on how to integrate penalty functions into minimization using lmfit. 

Sincerely,

Cong

[Matt Newville]

Hi Cong, 

I think you're close to the right way to implement a penalty to the fit.   This is sometimes called "regularization".  I think some other uses call this a "restraint" (as opposed to a "constraint"), or even just a lambda multiplier. 

The minimization routines do the sum-of-squares given the residual array.   So you just need to append the appropriate value to the residual arrays so that it will give the value you want when squared.   That will also avoid the cross-terms that you worry (rightly so) about. 

In code, I would suggest something like in your objective function for `minimize` (this won't work with Models) 

     residual_array = ymodel - ydata

  

     # penalty as a 1 element nd.array:

     # notes:   

     #      1. "lambda" is a builtin keyword in Python - you'll need a different name

     #      2.  I *think* you mean that 'a' and 'b' will be nd-arrays

     #      3.  If this penalty was a single, scalar value, do 

     #                  penalty_array = np.array([penalty_scalar_value])

     #           to turn it into a 1-element array so that 'concatenate' below will work

     penalty_array = lam_factor * (a*a - b*b -1)    # assumes 'a' and 'b' are arrays       

     residual_with_penalty = np.concatenate((residual_array, penalty_array))

     return residual_with_penalty

For sure, tuning the value of 'lambda` can be challeging.  But if you think of it as a weighting factor for your penalty versus the fit to the data, it can make it interpretable as "who are you going to believe: the data or the expectation that your penalty demands"?

I say it won't work with Models.  If you set on using that, you *could* write a Model function that created a model ndarray and then concatenate any penalty values onto the array returned by the model.  But then you would also have to create a "data_for_penalized_model" array that concatenated the right number of zeros onto the actual data.  I think that is not any easier and since it means your fitting script has to account for the penalty in two different places, it's probably more error-prone.

Hope that helps,

--Matt

      

--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/4a103009-7736-4f51-848b-5cb3945448d2n%40googlegroups.com.

[Cong Wang]

Dear Matt,

Thank you very much!  It really makes a lot of sense to append results of penalty terms to residuals.  So the final residual array would look like [residual1, residual2, ..., residualN, penalty1, penalty2, ..., penaltyN]. 

I will give it a try, and perhaps struggle with tuning the lam_factor.  Nevertheless, your idea is absolutely a better way (perhaps the most elegant) to integrate penalty functions into minimization. 

It surely helps!  Thank you again!

Cong