How to use `NonlinearSolve` to perform least squares fitting?

General Usage
1

I am trying to use NonlinearSolve to perform least squares fitting, as shown in this example. One difference from that is in the model function, I have performed some operations on a matrix. And as shown in the figure below, this caused some problem (not sure if this is the crux here):

image
image1506×750 99.9 KB

Is there some errors I have made during this process?

Here is my code:

# --- some helper function ---
function f_gen(k)
	τs = [[1/3, -2/3], [1/3, 1/3], [-2/3, 1/3]]
	fs = map(τs) do τ return cis(2π * dot(k, τ)) end
	return sum(fs)
end

# --- definition of the model ---
function model(pts2d, p)
    γ0, V1 = p
    ys = Array{Float64}(undef, 2, size(pts2d, 1))
    for i in axes(ys, 2)
        f = f_gen(vec(pts2d[i, :]))
        ys[:, i] = eigvals([V1 -γ0*f; -γ0*conj(f) V1])
        # the following line is for testing
        # ys[:, i] = real.(sin.(sum([V1 -γ0*f; -γ0*conj(f) V1], dims=1)))
    end
	# @show ys
    return vec(ys)
end

# --- the fitting not working ---
Random.seed!(1234)
# - params -
γ0 = 2.5
V10 = -0.097611
# - generate testing data -
pts2d = rand(10, 2)
ys0 = model(pts2d, [γ0, V10])
ys = ys0 .+ randn(size(ys0)) .* 0.01
# - fit γ0, V10 -
function residual(p, data)
	# @show data[2]
    pts2d = data[1]
    ys = data[2]
    ys_pred = model(pts2d, p)
    return ys_pred .- ys
end
prob = NonlinearLeastSquaresProblem(residual, [γ0, V10], (pts2d, ys))
res = solve(prob, LevenbergMarquardt())

Any help is warmly welcomed!

2

The problem in the provided snippet (which is different from the stack trace you show, which I guess comes from a different iteration of the script) is that you’re preallocating ys as an Array{Float64}, which is correct when you call the function with Float64 inputs. However, when you pass this function to a NonlinearSolve solver, it attempts to differentiate it, which defaults to a method based on automatic differentiation and is provided by the ForwardDiff.jl library. This library calls the function with arrays of a different element type, a Dual number, containing the derivative information (if you’re not familiar with this topic look up “forward-mode automatic differentiation”). The code then fails when you demand the storage of the Dual number into a Float64 array. If you’re OK with using ForwardDiff, you can bypass this issue by using a more “functional” approach in which you don’t preallocate anything, but rather let the compilation process find out which array is created depending on the input types, e.g.

function model(pts2d, p)
    γ0, V1 = p
    map(eachrow(pts2d)) do row
        f = f_gen(row)
        real.(sin.([V1 - γ0 * conj(f), V1 - γ0 * f]))
    end |> Base.Iterators.flatten |> collect
end

Once you get the script working, there is some stuff that can be done as well to improve performance, if that matters to your project.

3

Hey Daniel, thanks for you reply!

I have tried your method, it does work for the real.sin. function. In fact, I tried this line to test if the eigvals is forbidden in autodiff and mistakenly copied it here (sorry), thus the discrepancy between the code and the stack info.

I have edited in the original question what I am trying to do. Unfortunately, it still doesn’t work even if I have tried your suggenstion. And the stack info stays the same. Would you please take a look on that?

4

The problem of the uncompatability of eigvals and NonlinearSolve is solved by add using GenericLinearAlgebra.