How to fix some parameters when least-squares-fitting using `NonlinearSolve`?

General Usage
1

I have 2 questions:

  1. I find there are two different practices when performing lsqfit. One defines a function without a sideeffect, as done here; the other defines a function that mutates a predefined variable, examplified here. I think there are some performance differences. Except this, is there some other differences I should know before using any of them?
  2. In the syntax of LsqFit, I can easily fix some of the parameters of the model I have defined, e.g., through a template lambda. Is there similar syntax of NonlinearSolve to avoid redefine the form the model?
2

This example does not mutate a pre-defined variable Sorry, I was looking at the wrong link.

Avoiding allocations (i.e. allocating your temporary arrays once, and using an API where the caller passes in a buffer to mutate) often helps performance.

This syntax is not specific to LsqFit — you can pass a closure (an anonymous function) to any higher-order function, and use it to capture additional parameters or fix some fo the parameters of the model.

However, if you want to use the caller-supplied buffer in-place API, your closure will need to be a little more complicated in order to change the number of parameters. (You will need to pre-allocate a buffer of the correct size, and your closure will need to copy to/from the buffer.)

3

Thanks for you reply!

In fact, what I am trying to do is as following:

# -- define the model --
function expmodel(x, u)
    y0, a, τ, t₀ = u
    return y0 + a * exp(-(x-t₀)/τ)
end
# and unified residual function
function residualfn!(du, u, data)
    (xs, ys) = data   
    du .= expmodel.(xs, Ref(u)) .- ys
    return nothing
end

# -- so that I can fix arbitary parameter at ease --
# -- generate data --
Random.seed!(3465)
y0, a, τ, t₀ = 0.0, 1.0, 2.0, 0
xs = 0:0.1:3;
ys0 = expmodel.(xs, Ref([y0, a, τ, t₀]))  # ideal data
ys = ys0 .+ randn(size(ys0)) .* 0.02      # add noise
# -- fit --
prob1 = NonlinearLeastSquaresProblem(
        NonlinearFunction((du, u, data) -> residualfn!([du; 0.0], [u; t₀], data), resid_prototype = similar(ys)), [y0, a, τ], (xs, ys))

prob2 = NonlinearLeastSquaresProblem(
        NonlinearFunction((du, u, data) -> residualfn!([du[1:2]; 0.0; du[3]], [u[1:2]; t₀; u[3]], data), resid_prototype = similar(ys)), [y0, a, t₀], (xs, ys))

prob3 = NonlinearLeastSquaresProblem(
        NonlinearFunction((du, u, data) -> residualfn!([du; 0.0; 0.0], [u[1:2]; τ; t₀], data), resid_prototype = similar(ys)), [y0, a], (xs, ys))

But as you see, all three trials won’t work. And I couldn’t figure it out. Would you please have a look on that?

4

You go through all of this work to make your residualfn! be non-allocating, and that’s all well and good. But this anonymous function you put on top of it is allocating, specifically you’re making two arrays, [du; 0.0], [u; t₀]. You may want to allocate these vectors once and instead just reuse the cache memory, or rewrite residualfn! so that it uses du and u directly, like du .= expmodel.(xs, Ref(u), t₀) .- ys where you just keep that t₀ separate.

1 Like
5

Not sure if I understand you correctly.
residualfn! only modify du in-place, so if du is preallocated, it is easy to be retrieved from the cache and thus have a performance gain. But I am bit confused, if no parameter in the model is fixed, the problem is defined as prob = NonlinearLeastSquaresProblem(NonlinearFunction(lossfn!, resid_prototype = similar(ydata)), u0, data) (as here, and du is not preallocated but a performance gain is expected.
Here in my case, where some parameters of u are transfered by closure:
prob = NonlinearLeastSquaresProblem(NonlinearFunction((dp, p, data) -> residual!(dp, [p; t₀], data), resid_prototype = similar(ys)), [y0, a, τ], (xs, ys))
I don’t get how I lose the performance gain.
Would you please demonstrate more? I would really appreciate that!