Stochastic results in 2-dimensional biquadratic factorization

[timholy]

If you run the following script multiple times, you'll see you get different results each time, and it only rarely finds the valid real solution:

using HomotopyContinuation
using StaticArrays

function factorbiq(g0, g1, H0, H1; kwargs...)
    @var A0[1:3] b0[1:2] A1[1:3] b1[1:2]
    A0m, A1m = @SMatrix([A0[1] A0[2]; A0[2] A0[3]]), @SMatrix([A1[1] A1[2]; A1[2] A1[3]])
    b0v, b1v = @SVector([b0[1]; b0[2]]), @SVector([b1[1]; b1[2]])
    e = SVector(1, 1)
    q0e = (e' * A0m * e) / 2 + e'*b0v + 1
    q1z = (e' * A1m * e) / 2 - e'*b1v + 1
    H0eq = q1z * A0m + b0v * (-A1m * e + b1v)' + (-A1m * e + b1v) * b0v' + A1m - H0
    H1eq = A0m + b1v * (A0m * e + b0v)' + (A0m * e + b0v) * b1v' + q0e * A1m - H1
    g0eq = q1z * b0v - A1m * e + b1v - g0
    g1eq = q0e * b1v + A0m * e + b0v - g1
    # TODO: elide the steps of the polyalgorithm to make this more precompilable
    s = System([
        g0eq[1]; g0eq[2]; g1eq[1]; g1eq[2];
        H0eq[1,1]; H0eq[1,2]; H0eq[2,2];
        H1eq[1,1]; H1eq[1,2]; H1eq[2,2]
    ])
    result = solve(s; show_progress=false, kwargs...)
    # Return the "most real" solution
    rbest, ratbest = nothing, -Inf
    for r in result
        rat = norm(real(solution(r))) / norm(solution(r))
        if rat > ratbest
            rbest, ratbest = r, rat
        end
    end
    return rbest
end

g0, g1 = [-2.0572295698118172e9, -1.0946872831524246e10], [-2.654355654346107e9, -1.8026431499130545e9]
H0 = [-5.799857157583843e9 2.2112130890346775e9; 2.2112130890346746e9 1.1904739805563393e10]
H1 = [-9.605191673340421e8 2.881331173461073e9; 2.881331173461075e9 1.971553859807829e9]

while true
    global r
    r = factorbiq(g0, g1, H0, H1)
    @show r
    is_nonsingular(r) && is_real(r) && break
end

Note that the equations we're trying to solve are the same each time, and we keep trying until we find a non-singular real solution.

As one example, I got the following:

julia> include("hcbug.jl")
r = PathResult:
 • return_code → :at_infinity
 • solution → ComplexF64[-1.6558468620230283e13 + 2.816040161764201e10im, 1.655824946429631e13 - 2.8401426125692253e10im, -1.6559491629477156e13 + 3.1040601728790916e10im, -948.4051952664054 - 534.7005410844648im, -29260.39143293631 - 104633.62535005403im, 70797.59697366082 + 58920.71165946572im, 2.0911060021426907e8 + 3.897423756070992e8im, 7.452251218640794e8 - 1.9997840668520076e9im, 5667.252738088702 - 75441.9970419194im, -0.1624273916332903 - 2.021261614346973im]
 • accuracy → NaN
 • residual → 628.7
 • condition_jacobian → NaN
 • steps → 283 / 2
 • extended_precision → false
 • path_number → 52

r = PathResult:
 • return_code → :at_infinity
 • solution → ComplexF64[114821.46641677551 + 96919.40948489086im, -114823.24406097757 - 96919.65294483457im, 114825.60314687753 + 96921.59469198558im, -1.6013675947899613e9 - 8.038964125573374e8im, 4.761442621083755e9 + 4.846468403082303e8im, 1.8087067179127766e10 - 1.599135984012342e8im, 1.7418029217087203 + 1.171981216973091im, -2.8955222700932035 - 1.7571167461525765im, 1.1027666635600133e9 - 3.1941767616596586e8im, 1.190150608621988e10 + 3.2216620510423446e8im]
 • accuracy → NaN
 • residual → 1.3605
 • condition_jacobian → NaN
 • steps → 182 / 0
 • extended_precision → false
 • path_number → 68

r = PathResult:
 • return_code → :at_infinity
 • solution → ComplexF64[-7957.989197002466 + 20479.522562566282im, -12231.978437307005 - 7590.600577618955im, 8528.17253971766 + 5743.541243434932im, -3.3479710463098405e15 - 1.1502613927514392e14im, 3.347846492774834e15 + 1.1485966368494695e14im, -3.347572674757794e15 - 1.1473622448925194e14im, 11944.118824901052 - 5520.5739634179945im, 1.6680640660961203 - 0.38723271196744835im, 5.127382665448743e9 + 1.9575881115611458e10im, 6.37125409481006e10 - 2.1914988619238934e10im]
 • accuracy → NaN
 • residual → 27947.0
 • condition_jacobian → NaN
 • steps → 232 / 0
 • extended_precision → false
 • path_number → 185

r = PathResult:
 • return_code → :at_infinity
 • solution → ComplexF64[-797567.3012563197 - 7809.499682550331im, 797550.8981493694 + 7165.001552222063im, -797533.7550200477 - 6575.286166216206im, 4.06657782294473e11 + 2.2480874140239227e10im, -5.545866936260137e11 + 8.05199901720135e8im, 7.152488830169773e11 - 2.5300310836770287e10im, -4.112238701210074 + 636.1333845969035im, 2.684027144599631 - 608.3284962982227im, 1.3168448568716762e9 + 5.023486395936862e9im, 5.019583771816573e9 - 5.865923560223413e9im]
 • accuracy → NaN
 • residual → 130.5
 • condition_jacobian → NaN
 • steps → 269 / 0
 • extended_precision → false
 • path_number → 145

r = PathResult:
 • return_code → :at_infinity
 • solution → ComplexF64[-1.9020709256201603e7 + 2.5088367720778165e6im, 8.893457840818468e6 - 105055.07284417136im, 1.236728454201505e6 - 2.306738474171285e6im, -1.1752130360144535e10 - 1.4410097809620981e9im, 1.8624812241783946e9 + 2.1957191690779352e8im, 8.021838030536382e9 + 1.008087819224245e9im, 263.026407383105 + 3386.2103443089704im, 395.3071470214064 + 235.34350393805923im, -1.2105020403715959e6 - 216872.9676174307im, -826988.4928909729 - 151217.15295182753im]
 • accuracy → NaN
 • residual → 48.766
 • condition_jacobian → NaN
 • steps → 305 / 0
 • extended_precision → false
 • path_number → 235

r = PathResult:
 • return_code → :terminated_max_extended_steps
 • solution → ComplexF64[2.1135364307840996e10 + 1.3615163065341587e9im, -8.791196249456844e9 + 2.0974565663771026e9im, -3.5531959413637304e9 - 5.556657944298898e9im, -6.034221902435255e10 + 1.8673577110093265e9im, 5.1890645506695366e10 + 3.598089720128142e8im, -4.343928263493715e10 - 2.5870326353308268e9im, 113904.89671807428 + 6452.09941490103im, 19367.339269131215 + 29779.036266830986im, -7374.132377167348 - 190232.4640300153im, -7087.312105954894 + 137040.32637638677im]
 • accuracy → NaN
 • residual → 144380.0
 • condition_jacobian → NaN
 • steps → 680 / 0
 • extended_precision → true
 • path_number → 157

r = PathResult:
 • return_code → :success
 • solution → ComplexF64[0.5424262561146389 + 4.484155085839415e-43im, -0.7222287649638645 - 9.865141188846712e-43im, 0.9349262099012551 + 1.2555634240350361e-42im, 3.4025057207910395e9 + 7.703719777548943e-34im, 2.1665621707498975e9 - 2.8888949165808538e-34im, -1.4489885337353756e9 - 2.792598419361492e-33im, 0.7903552508081598 + 1.5694542800437951e-43im, -0.9617825421313116 - 3.1389085600875902e-43im, -3.1411742915498037e9 - 8.546314128218359e-34im, -2.1332545645940182e9 - 6.048623731591163e-34im]
 • accuracy → 5.2966e-15
 • residual → 1.9073e-6
 • condition_jacobian → 2358.5
 • steps → 261 / 0
 • extended_precision → false
 • path_number → 95

It succeeded on the last try (which is why it quit), but if found different non-real solutions on each prior attempt.

In my application, I'm bothered somewhat by the non-reproducibility of the solutions, but I'm bothered even more if it fails to find a genuine solution.

This could partially be a duplicate of #647 (which is much simpler than this example), but I wasn't sure so I thought it would be better to file a separate issue.

[saschatimme]

Hi Tim,

I think there is a better approach to solve your problem using the monodromy_solve function and some (offline) precomputation.

using HomotopyContinuation
using StaticArrays

function factorbiq_system()
  @var A0[1:3] b0[1:2] A1[1:3] b1[1:2]
  @var g0[1:2] g1[1:2] H0[1:2, 1:2] H1[1:2, 1:2]
  A0m, A1m = @SMatrix([A0[1] A0[2]; A0[2] A0[3]]), @SMatrix([A1[1] A1[2]; A1[2] A1[3]])
  b0v, b1v = @SVector([b0[1]; b0[2]]), @SVector([b1[1]; b1[2]])
  e = SVector(1, 1)
  q0e = (e' * A0m * e) / 2 + e' * b0v + 1
  q1z = (e' * A1m * e) / 2 - e' * b1v + 1
  H0eq = q1z * A0m + b0v * (-A1m * e + b1v)' + (-A1m * e + b1v) * b0v' + A1m - H0
  H1eq = A0m + b1v * (A0m * e + b0v)' + (A0m * e + b0v) * b1v' + q0e * A1m - H1
  g0eq = q1z * b0v - A1m * e + b1v - g0
  g1eq = q0e * b1v + A0m * e + b0v - g1
  # TODO: elide the steps of the polyalgorithm to make this more precompilable
  System([
      g0eq[1]; g0eq[2]; g1eq[1]; g1eq[2];
      H0eq[1, 1]; H0eq[1, 2]; H0eq[2, 2];
      H1eq[1, 1]; H1eq[1, 2]; H1eq[2, 2]
    ]; parameters=[g0; g1; H0[:]; H1[:]])
end

F = factorbiq_system()

# First solve using monodromy
R0 = monodromy_solve(F)
# THe system has 6 solutions in general.
# This is non-deterministic, so we need to run it a few times possibly. 
for r in 1:10
  if nsolutions(R0) == 6
    break
  end
  R0 = monodromy_solve(F)
end

# Now we have 6 solutions. Store the good solutions on disc
S0 = write_solutions("S0.txt", solutions(R0))
# Now we have 6 solutions. Store the good solutions on disc
p0 = write_parameters("p0.txt", parameters(R0))

# In a clean session we can just load the solutions and parameters
# and solve the system.
S0 = read_solutions("S0.txt")
p0 = read_parameters("p0.txt")

# Now solve for the specific values we are interested in.
g0, g1 = [-2.0572295698118172e9, -1.0946872831524246e10], [-2.654355654346107e9, -1.8026431499130545e9]
H0 = [-5.799857157583843e9 2.2112130890346775e9; 2.2112130890346746e9 1.1904739805563393e10]
H1 = [-9.605191673340421e8 2.881331173461073e9; 2.881331173461075e9 1.971553859807829e9]

p1 = [g0; g1; H0[:]; H1[:]]

R1 = solve(F, S0; target_parameters=p1, start_parameters=p0)
real_solutions(R1)

This returns 2 solutions

julia> real_solutions(R1)
2-element Vector{Vector{Float64}}:
 [-2.7940017822215576e9, 6.237312709437258e9, -4.353989641624983e9, -0.01053114958353485, -0.19238146498792302, 0.4341175976543101, -1.753008182817548e9, -9.328058435099428e9, 0.5161306281585974, -0.6702613492291367]
 [0.542426256114644, -0.7222287649638696, 0.9349262099012599, 3.402505720791062e9, 2.1665621707498913e9, -1.448988533735391e9, 0.7903552508081622, -0.9617825421313141, -3.1411742915498056e9, -2.1332545645940192e9]

Looking at the scaling of your parameters it could also make sense to find a "generic" solution that is similarly scaled but I don't have to little insight in your specific problem for that.

[PBrdng]

Parameter homotopy with computing the initial solutions by solve also seems to work just fine:

s = System([
      g0eq[1]; g0eq[2]; g1eq[1]; g1eq[2];
      H0eq[1, 1]; H0eq[1, 2]; H0eq[2, 2];
      H1eq[1, 1]; H1eq[1, 2]; H1eq[2, 2]
    ]; parameters=[g0; g1; H0[:]; H1[:]])


g00, g10 = [-2.0572295698118172e9, -1.0946872831524246e10], [-2.654355654346107e9, -1.8026431499130545e9]
H00 = [-5.799857157583843e9 2.2112130890346775e9; 2.2112130890346746e9 1.1904739805563393e10]
H10 = [-9.605191673340421e8 2.881331173461073e9; 2.881331173461075e9 1.971553859807829e9]

interesting_parameters=[g00; g10; H00[:]; H10[:]]
random_parameters = randn(ComplexF64, 12) 

result0 = solve(s; target_parameters = random_parameters)
result = solve(s, result0, start_parameters = random_parameters, target_parameters = interesting_parameters)

Notice that the first solve producing result0 tracks 252 paths, while the second solve only tracks 6 paths. This is why in both Sascha's and my solution we first solve a random instance and then track towards your parameters. This is also the numerically more stable way to solve.

[PBrdng]

I will move this into the discussion section, since it's not really an issue or a bug, but rather the way homotopy continuation works (which can be subtle sometimes).

[PBrdng]

One more comment: If you solve a system more than one time for multiple parameters, you should follow this guide.

[timholy]

Want to add (as thanks, and in case anyone else reads this) that I've now implemented parameter caching. I wish I had done it earlier: I'm particularly impressed at the improved precision, that and the fact that solutions are now 100% repeatable are doing wondrous good for my application.

Many thanks for the advice and demonstrations!

[PBrdng]

I'm glad it works for you now :)

[PBrdng]

There are many other tricks which one can use in various situations. If you find yourself in a similar situation ever again dont hesitate to reach out.