Vector valued y

[benneti]

What is the best way to fit a function $f: \mathbb{R} \to \mathbb{R}^n$.
Right now I use something in the line of (best expained by a random example with a function that should be hard to fit)

function _model(x,p)
    eigvals(p[1] * [1 0;0 1] * x +  p[2] * [0 1;1 0] * x)
end
_y(x) = _model(x, [1,2])

function model(xs,p)
    (x->norm(_model(x, p) - _y(x))).(xs)
end
fit = curve_fit(model, 0:100, zero(0:100), [0.4,1.4]) # hard problem to fit as every p with the correct norm could work
coef(fit)

which kinda works, but also seems wrong and not like the intended way.
(Everything that looks better like trying to fit the function directly where every y is a vector fails, either due to isnan or dimension mismath or cannot convert vector to float).

[pkofod]

Do you just want to flatten the vector before the norm-call , or?

[benneti]

I do not follow your suggestion to be honest.
_model are already 1D so as I understand it there is nothing to flatten.
My question was more in the direction whether there is a more elegant/direct way instead of having the wrapper and helper functions.

[pkofod]

I'm just not sure I understand your problem here.

[benneti]

Basically I would like to know whether there is a more elegant way to use the curve_fit for vector functions, as the readme sounds like this could work native.
(On the other hand my example takes the vector function and maps it to a scalar function which works but seems like there would be room for improvement)

Thank you for trying to help!

[pkofod]

Is this the same as here ? #156

[benneti]

It is probably related, but I think it would go a bit further because as a simple check I also tried to fit the fector function but overloaded

Base.isinf(x::Array{Array{Float64,1},1}) = any(isinf.(x))
Base.isinf(x::Array{Float64,1}) = any(isinf.(x))
Base.isnan(x::Array{Array{Float64,1},1}) = any(isnan.(x))
Base.isnan(x::Array{Float64,1}) = any(isnan.(x))

But it still does not work.

Could it be that the cost function f = (p) -> model(xdata, p) - ydata is assumed to be a scalar?
Else it pobably is just that I do not understand how to generalize the broadcasting to multi dimensional output.
Could you provide an example similar to
@. multimodel(x, p) = p[1]*exp(-x[:, 1]*p[2]+x[:, 2]*p[3]) but instead going from R^2 to R from R to R^2?

[pkofod]

I think we need to go back to basics here. What do you mean R to R^2? That's not a nonlinear least squares problem. How do you want to weigh together the two elements of the output? This sounds like either something like a multivalued nonlinear least squares or a system of equations of sorts? Can you describe the problem you're trying to solve, not the problems you're encountering in terms of Julia errors.

[benneti]

Ah OK, I guess I was confused because the Jacobian section of the documentation is kept general for a function f that can go from R^m to R^n (see https://julianlsolvers.github.io/LsqFit.jl/latest/getting_started/) so I just assumed the curve_fit function would be able to handle function of this form. If I know understand you correctly it does only support R^m to R and in this case the issue is non existent and can be closed. The manual wrapper function is then still a working solution.

[pkofod]

I just need to understand what your mathematical problem is, I mean what is the technical or scientific problem you're trying to solve. Without that I feel like we cannot understand each other :) That jacobian is of the vector of model function evaluations, not the objective function itself. The objective function is the least squares function of the residuals (model minus data).

[benneti]

Ah ok, basically I have several energy levels as a funtion of one control parameter and I am trying to fit the eigenvalues of a model Hamiltonian as a function of the same control parameter and several system properties (p[:]) to it.
Or in more mathematical termsI have the energy lines E[i](x) and want to model them by eigenvals(H(x,p))[i] where E[i](x) is real and therefore H(x,p) is a real and symmetric matrix of size length(E)xlength(E).

Now to link to the initial post I was not able to figure out how to fit this directly but I knew how to fit a function x->f(x) where f(x) in R therefore I simply created my own cost function f=norm(E-eigenvals(H))that should take the value 0 everywhere with which the fit worked fine. But reading the docs I assumed that this should be possible without defining the cost function myself.

[pkofod]

I'm not really sure I understand, I'm quite sorry. Maybe bring this up in #math-optimization on the julia slack or on the Julia discourse. You can tag me there, and I'll try to follow the discussion :)