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

**URL:** <https://discourse.julialang.org/t/how-to-use-nonlinearsolve-to-perform-least-squares-fitting/124449>\
**Category:** General Usage\
**Tags:** question\
**Created:** [January 5, 2025, 3:52pm UTC](https://discourse.julialang.org/t/how-to-use-nonlinearsolve-to-perform-least-squares-fitting/124449 "2025-01-05T15:52:24Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![thudjx](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thudjx/32/15231_2.png) [@thudjx](https://discourse.julialang.org/u/thudjx)\
**Post date:** [January 5, 2025, 3:52pm UTC](https://discourse.julialang.org/t/how-to-use-nonlinearsolve-to-perform-least-squares-fitting/124449/1 "2025-01-05T15:52:24Z")

</div>

I am trying to use `NonlinearSolve` to perform least squares fitting, as shown in this [example](https://discourse.julialang.org/t/curve-fitting-with-nonlinearsolve-jl/119238/10). 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](https://global.discourse-cdn.com/julialang/original/3X/7/a/7a7568cfba1f3c133b7c0b95be2e1c295c633c7f.png)

Is there some errors I have made during this process?

Here is my code:

```julia
# --- 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!

---

<div class="post-metadata">

**Author:** ![Daniel\_Gonzalez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/daniel_gonzalez/32/52412_2.png) [@Daniel\_Gonzalez](https://discourse.julialang.org/u/Daniel_Gonzalez)\
**Post date:** [January 5, 2025, 8:24pm UTC](https://discourse.julialang.org/t/how-to-use-nonlinearsolve-to-perform-least-squares-fitting/124449/2 "2025-01-05T20:24:31Z")

</div>

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.

```julia
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.

---

<div class="post-metadata">

**Author:** ![thudjx](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thudjx/32/15231_2.png) [@thudjx](https://discourse.julialang.org/u/thudjx)\
**Post date:** [January 6, 2025, 1:20am UTC](https://discourse.julialang.org/t/how-to-use-nonlinearsolve-to-perform-least-squares-fitting/124449/3 "2025-01-06T01:20:15Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![thudjx](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thudjx/32/15231_2.png) [@thudjx](https://discourse.julialang.org/u/thudjx)\
**Post date:** [January 6, 2025, 12:50pm UTC](https://discourse.julialang.org/t/how-to-use-nonlinearsolve-to-perform-least-squares-fitting/124449/4 "2025-01-06T12:50:30Z")

</div>

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