Difficulty with 4 doors example, expecting roughly equal modes

[itsDionix]

I am struggling to reproduce the results of the 4 doors example from "Multi-modal and Inertial Sensor Solutions for Navigation-type Factor Graphs", Fourie 2017. This might be related to #1570.

Included is an adapted Minimal Working Example (MWE) and two example figures of the output of the MWE with N=5000 and N=1000. The example is simplified to only 3 doors.

Specifically, I am trying to reproduce the two roughly equal modes in Figure 6-16 from Fourie's 2017 thesis, cited above. The belief for x3 I am getting from the automatic initialization in the MWE (orange dotted line in the attached figures) presents two modes that are similar, which is good. However, the issue is that proceeding to solve the graph makes the result worse, most of the time. The solved belief (solid orange line in the attached figures) places almost all probability on the left mode, when both should be equal.

At the moment, is there a way to get better performance (where both modes are roughly equal) when solving the full graph?

I am seeking to report the results in my Master's thesis, so I want to make sure that I am getting representative results.

Thank you for the time and any help you can give.

I used Julia 1.9 and IncrementalInference 0.34.0.
Next follow example results and the MWE code.

Example output with N=5000:

Example output with N=1000:

MWE:

using IncrementalInference
using Plots
cv = 3.0
components = [Normal(-100, cv); Normal(0, cv); Normal(100, cv)]
doorPrior = Mixture(Prior,
    components,
    [1 / 3; 1 / 3; 1 / 3])

## Build the factor graph object
fg = initfg()
fg.solverParams.N = 5000   # 500 is enough to see the phenomenon, used more just in case
# fg.solverParams.gibbsIters = 30 # Increasing from the default doesn't seem to help

gr(dpi=300)
plt = Plots.plot()
xlabel!("\$x\$")
ylabel!("Likelihood")
box = (-130, 230)
x = range(box[1], box[2], 2001)
xlims!(box[1], box[2])

# first pose location
v1 = addVariable!(fg, :x1, ContinuousScalar)
# see a door for the first time
addFactor!(fg, [:x1], doorPrior)
## drive to second pose location
addVariable!(fg, :x2, ContinuousScalar)
addFactor!(fg, [:x1; :x2], LinearRelative(Normal(50.0, 2.0)))
# drive to third pose location
v3 = addVariable!(fg, :x3, ContinuousScalar)
addFactor!(fg, [:x2; :x3], LinearRelative(Normal(50.0, 4.0)))
# see a door for the second time
addFactor!(fg, [:x3], doorPrior)

initAll!(fg)

bel = getBelief(fg, :x1).belief
plot!(x, (x) -> bel([x])[1], label="x1 init", ylims=(0, Inf), linestyle=:dot, c=1)
bel = getBelief(fg, :x3).belief
plot!(x, (x) -> bel([x])[1], label="x3 init", ylims=(0, Inf), linestyle=:dot, c=2)
display(plt)

solveTree!(fg)

bel = getBelief(fg, :x1).belief
plot!(x, (x) -> bel([x])[1], label="x1 solved", ylims=(0, Inf), c=1)
bel = getBelief(fg, :x3).belief
plot!(x, (x) -> bel([x])[1], label="x3 solved", ylims=(0, Inf), c=2)
savefig("MWE.png")
display(plt)

[dehann]

Hi @itsDionix , apologies for the late reply. This example is in the CI tests. There are two ways to induce the multidoor example, either through a mixture prior (as you have) or via the addFactor!(..., multihypo=[1;0.25;0.25;0.25;0.25]).

automatic initialization
proceeding to solve the graph makes the result worse

The initialization step is akin to the prediction step used by filtering techniques which will produce all modes quite easily. The solve is supposed to retain those modes in roughly equal modes, but prune out less likely ones. The purpose of parameter N is how many likely modes should be kept. The current stable branch code requires around 30 kernels per mode, so roughly 3 modes in 100 samples. The future of the solver code here will greatly improve on this.

Since you are getting less than ideal results it is probably a good idea on our end to review and improve the mixture belief performance. I'll add this into the roadmap, but is will only happen after two other known issues are resolved.

The best way for us to confirm is to make this test more stringent (the multihypo method):

IncrementalInference.jl/test/testMultiHypo3Door.jl

Line 57 in 4bf5d1a

f1 = addFactor!(fg, [:x0; :l0; :l1; :l2; :l3], LinearRelative(Normal(0, meas_noise)), multihypo=[1.0; (1/4 for _=1:4)...])

At the moment, is there a way to get better performance (where both modes are roughly equal) when solving the full graph?

Not in the next few weeks, however, we have been reworking a lot of fundamental in the solver code over the past months to improve other issues. These upgrades are about halfway and may improve the performance as they get merged in.

In the mean time, here are a few suggestions of this to check and try for your example:

• Multiple runs amplify this or that mode, and not always asymmetry towards the same mode.
• getSolverParams(fg).inflationCycles = 5 # normal range 1 - 5
• getSolverparams(fg).inflation = 10.0 # normal range 1 - 1000.0
• getSolverParams(fg).gibbsIters = 1 # normal range 1-5
• getSolverParams(fg).nullSurplusAdd = 0.5 # normal range 0.0 - 0.7
• getSolverParams(fg).N = 200 # normal range 75-500

I'd like to keep this issue open until this is better resolved, please do post any additional findings which might help find a bug. I'll also add that the code related to multihypo has been greatly refactored in the past months and we are relying on CI testing to catch bugs, but some bugs have slipped through in the past. Overall though the code quality and maintainability has been getting much better, but we have more work to do.

[dehann]

This might be related to #1570.

Maybe similar in issue name, but that problem was resolved during some of the earlier work. I'd recommend the focus of this issue thread look into "expecting near equal modes".

[itsDionix]

Thank you for the help!

My experiments with the suggested settings did not result in an improvement. In case it helps, I will leave documented that:

With default parameters and N=300, most times the left mode in solved x3 is the biggest of the two - this was the case in 10 out of 10 runs.

With default parameters, N=300 and gibbsIters = 1, the asymmetry in the modes of solved x3 is reduced, but the solved x1 has 3 modes.

I also tried defining the mixture with MixtureModel from Distributions, via:

doorMix = Distributions.MixtureModel(components, [1 / 3; 1 / 3; 1 / 3])
doorPrior = Prior(doorMix)

This did not result in a noticeable difference, though.

[dehann]

Thanks for the feedback.

most times the left mode in solved x3 is the biggest of the two - this was the case in 10 out of 10 runs.

This is definitely a bug then. The results should have been better (and did not have this behavior in the past).

With default parameters, N=300 and gibbsIters = 1, the asymmetry in the modes of solved x3 is reduced, but the solved x1 has 3 modes.

This might be due to a parasitic effect and was previously resolved with more iterations. This is moot until the bug above of "biased to one side" is resolved.

I also tried defining the mixture with MixtureModel... This did not result in a noticeable difference, though.

More confirmation that the issue is in the inference somewhere.

---

Not directly relevant but worth noting for general, IncrementalInference.Mixture should be refactored / removed / consolidated with Distirbutions.MixtureModel:

• Refactor IIF.Mixture as direct Density #1420

[dehann]

Also see:

• "Basic 3 door, trinary multihypothesis tests" often fails numerically #1570