Solve in 32 bits

[oameye]

Is their a way solve the system for 32 Floats?

using HomotopyContinuation
@var u1, v1, ω, α, γ, λ, ω0

eqs = [
    -u1*ω^2 + u1*ω0^2 + (3/4)*u1^3*α + (3/4)*u1*v1^2*α + (-1/2)*u1*λ*ω0^2 + v1*γ*ω,
    -v1*ω^2 + v1*ω0^2 + (3/4)*v1^3*α - u1*γ*ω + (3/4)*u1^2*v1*α + (1/2)*v1*λ*ω0^2,
]

F = System(eqs, parameters = [ω, α, γ, λ, ω0], variables = [u1, v1])

input_array = Float32[1.0, 1.0, 0.005, 0.05, 1.0]

sol = HomotopyContinuation.solve(
    F;
    target_parameters=input_array,
    threading=true
)
typeof(solutions(sol))

Vector{Vector{ComplexF64}}

[saschatimme]

I don't really understand. Are you looking for solutions with 32 bit precision? In this case just convert the result. We compute internally if needed even with higher precision than 64bit floating point.

[oameye]

I was hoping too gain some performance benefits from computing with single precision (Float32 type). However, looking at internals you guys always compute with 64bit floating point (double precision). Do you think it would be a lot of work to rewrite the internals to be bit agnostic?

[saschatimme]

This is quite a lot of work and I don't think worth it since you are risking lots of path failures ( = possibly missing solutions). I would rather look into high level algorithmic improvements like picking better start systems

[oameye]

Could you elude on why it would lead to path failures? Is it because the homotopy "step size" is too small?

[saschatimme]

https://github.com/NonlinearOscillations/HarmonicBalance.jl/blob/44349d505de8d537bba8b4fd409949fa9ce8afb4/src/solve_homotopy.jl#L241-L247

Is this the main part where you call HC and where your performance bottleneck is? I only glanced over your intro, but if you really do parameter sweeps then this is incredibly inefficient. I would pick a random complex point close to your sweep path, compute all solutions for this and then perform a parameter homotopy from this random point to your parameter point. (And this can even be improved further).

[saschatimme]

Could you elude on why it would lead to path failures? Is it because the homotopy "step size" is too small?

If a solution path passes a singularity closely then the condition number of the system increases, resulting in the need of high precision to still move forward with small steps. If precision is too low, you would not be able to move forward anymore since all computational results from the evaluation and/or linear algebra are just noise

[oameye]

Is this the main part where you call HC and where your performance bottleneck is? I only glanced over your intro, but if you really do parameter sweeps then this is incredibly inefficient. I would pick a random complex point close to your sweep path, compute all solutions for this and then perform a parameter homotopy from this random point to your parameter point. (And this can even be improved further).

What you describe is our default method, we dubbed WarmUp() (see here). However, this does not always find all solutions. Sometimes it misses some which leads to "noise" in for example the first plot of this tutorial.

Indeed doing a total_degree at every parameter point is significantly slower, but the end result always yields all solutions in our sweeps.

(And this can even be improved further).

Any suggestions? Maybe utilizing solutions nearby parameters points?

[oameye]

If a solution path passes a singularity closely then the condition number of the system increases, resulting in the need of high precision to still move forward with small steps. If precision is too low, you would not be able to move forward anymore since all computational results from the evaluation and/or linear algebra are just noise

Thanks for the clarification.

[saschatimme]

Interesting that there are so many failures. My guess would be that you are actually operating intentionally close to singularities in the parameter space and thus maybe more tricky. Definitely an interesting thing to look closer.

Couple ideas that I would try before giving up on the "warmup" method. These can be applied together / independently / partly:

1. After the first generic parameter solve perform a round of monodromy to possibly recover failed solutions.
2. Instead of a single set of generic start parameters generate a hand full and check the consensus on the number of solutions.
3. Don't sample completely random parameters but instead take the average of your parameters and add some coordinate wise random complex noise (relative to the magnitude of each coordinate). If the system has a low condition number around a parameter value then the solutions don't move much. So short distance between parameters = short distance between solutions.
4. When you do your sweep, as long as you find as many solutions as the generic root count is you can just reuse the homotopy result for the next parameter
5. In general it is probably faster to perform the "warmup" logic multiple times and take the highest result instead of doing the total degree

[saschatimme]

Also not 100% sure this is applicable, but for me it looks like you actually want to sample the discriminant of the steady state problem not the steady state space. I wrote once an article where we computed very similar drawings to your tutorial https://arxiv.org/pdf/2009.13408

[oameye]

Thank you! We will try some things out. We noticed it seems to explicitly not find the zero solution, so we could always easily test for the zero solution and it to the collection.

We indeed are interested in the discriminant (signaling a phase transition/bifurcation if happens for real solutions), but also the solutions themselves.