[float32] CG/Newton solver terminates early due to improvement overshoots #2313

n3b · 2024-12-27T10:24:04Z

Intro

My setup

What's happening? What did you expect?

I believe an abs() is missing here

Line 1824 in a8a9607

mjtNum improvement = scale * (oldcost - ctx.cost);

There are serious stability issues otherwise that can be seen in cases like throwing 10+ perfectly aligned cubes on top of each other.

Steps for reproduction

Minimal model for reproduction

No response

Code required for reproduction

No response

Confirmations

I searched the latest documentation thoroughly before posting.
I searched previous Issues and Discussions, I am certain this has not been raised before.

yuvaltassa · 2024-12-30T23:27:20Z

An abs is certainly not missing on that line. This is a signed quantity. If it is negative then there is no improvement.

If you can provide a MRE that would be helpful.

yuvaltassa · 2024-12-30T23:31:01Z

Ah but I think I see what you mean. On line 1833, a negative improvement would would lead to termination. I'll investigate. Thanks.

yuvaltassa · 2024-12-30T23:39:57Z

I think the correct thing is not an abs but rather max(0, ):

mjtNum improvement = scale * mju_max(oldcost - ctx.cost, 0);

Right?

Right now any increase in cost leads to termination. Perhaps this is by design? @emotodorov WDYT?

yuvaltassa · 2024-12-30T23:46:50Z

Ya I think a negative improvement is not possible since in this case we will have already terminated here:

mujoco/src/engine/engine_solver.c

Line 1779 in a8a9607

// no improvement: done

I'm closing this issue as you provided no way to reproduce, but I'm quite interested if you were to provide an example.

n3b · 2024-12-31T03:16:30Z

@yuvaltassa sorry, I should have mentioned I'm using single precision. You can see negative improvement in something like the xml below.

I also noticed that sometimes one-sided search can stuck at local minima and exhaust iterations with no improvement, while deriv[0] is above the tolerance. Would it make sense to check for delta in this loop as well?

mujoco/src/engine/engine_solver.c

Line 1292 in a8a9607

while (p1.deriv[0]*dir <= -gtol && ctx->LSiter < m->opt.ls_iterations) {

<mujoco>
    <option timestep="0.015"/>
    <option gravity="0 0 -9.81"/>
    <size memory="100M"/>

    <default>
        <geom solimp=".9 .99 .005"/>
        <default class="box">
            <geom type="box" material="box" size="0.5 0.5 0.5"/>
        </default>
    </default>

    <asset>
        <texture name="texgeom" type="cube" builtin="flat" mark="cross" width="128" height="128"
                 rgb1="0.6 0.6 0.6" rgb2="0.6 0.6 0.6" markrgb="1 1 1"/>
        <material name="box" texture="texgeom" texuniform="true" rgba=".4 .9 .9 1"/>
    </asset>

    <visual>
        <quality shadowsize="4096" offsamples="8"/>
        <map znear="0.1" force="0.05"/>
    </visual>

    <statistic extent="4"/>

    <worldbody>
        <light directional="true" diffuse=".8 .8 .8" pos="0 0 10" dir="0 0 -10"/>
        <geom type="plane" size="10 10 .01"/>

        <body euler="0 0 0" pos="3 0 10">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 10">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 13">
            <freejoint/>
            <geom class="box"/>
            
        </body>

        <body euler="0 0 0" pos="0 0 15">
            <freejoint/>
            <geom class="box"/>
            
        </body>

        <body euler="0 0 0" pos="0 0 17">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 19">
            <freejoint/>
            <geom class="box"/>
            
        </body>

        <body euler="0 0 0" pos="0 0 21">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 23">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 25">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 27">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 29">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 31">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 33">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 35">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 37">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 39">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 41">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 43">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 45">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 47">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 49">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 51">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 53">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 55">
            <freejoint/>
            <geom class="box"/>
        </body>

        <body euler="0 0 0" pos="0 0 57">
            <freejoint/>
            <geom class="box"/>
        </body>

    </worldbody>

</mujoco>

emotodorov · 2024-12-31T03:17:15Z

We are solving a convex optimization problem and have computed a direction which mathematically should be a descent direction. If it is not, that means we are so close to the minimum that floating-point errors are causing an increase in cost; in which case we should indeed terminate.

…

On Mon, Dec 30, 2024 at 3:40 PM Yuval Tassa ***@***.***> wrote: I think the correct thing is not an abs but rather max(0, ): mjtNum improvement = scale * mju_max(oldcost - ctx.cost, 0); Right? Right now any increase in cost leads to termination. Perhaps this is by design? @emotodorov <https://github.com/emotodorov> WDYT? — Reply to this email directly, view it on GitHub <#2313 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AB2V4L7RP7JH7UDIMJLEXF32IHKWJAVCNFSM6AAAAABUINWXSKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDKNRWGAYTANJYGM> . You are receiving this because you were mentioned.Message ID: ***@***.***>

emotodorov · 2024-12-31T03:22:14Z

Ah, like I said, floating-point errors. This is an interesting idea though, when using single-precision. If you do this, set an internal flag to tell you that you are now operating in a different regime where the math and convergence criteria etc. no longer apply. Maybe allow a couple of iterations in that regime and terminate? Not obvious; will have to be fine-tuned or real-world examples.

…

-- Emo

On Mon, Dec 30, 2024 at 7:16 PM Maxim Vorobey ***@***.***> wrote: @yuvaltassa <https://github.com/yuvaltassa> sorry, I should have mentioned I'm using single precision. You can see negative improvement in something like the xml below. I also noticed that sometimes one-sided search can stuck at local minima and exhaust iterations with no improvement, while deriv[0] is above the tolerance. Would it make sense to check for delta in this loop as well? https://github.com/google-deepmind/mujoco/blob/a8a960771234c8d84b085516d9f76a53477af28c/src/engine/engine_solver.c#L1292 <mujoco> <option timestep="0.015"/> <option gravity="0 0 -9.81"/> <size memory="100M"/> <default> <geom solimp=".9 .99 .005"/> <default class="box"> <geom type="box" material="box" size="0.5 0.5 0.5"/> </default> </default> <asset> <texture name="texgeom" type="cube" builtin="flat" mark="cross" width="128" height="128" rgb1="0.6 0.6 0.6" rgb2="0.6 0.6 0.6" markrgb="1 1 1"/> <material name="box" texture="texgeom" texuniform="true" rgba=".4 .9 .9 1"/> </asset> <visual> <quality shadowsize="4096" offsamples="8"/> <map znear="0.1" force="0.05"/> </visual> <statistic extent="4"/> <worldbody> <light directional="true" diffuse=".8 .8 .8" pos="0 0 10" dir="0 0 -10"/> <geom type="plane" size="10 10 .01"/> <body euler="0 0 0" pos="3 0 10"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 10"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 13"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 15"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 17"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 19"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 21"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 23"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 25"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 27"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 29"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 31"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 33"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 35"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 37"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 39"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 41"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 43"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 45"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 47"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 49"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 51"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 53"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 55"> <freejoint/> <geom class="box"/> </body> <body euler="0 0 0" pos="0 0 57"> <freejoint/> <geom class="box"/> </body> </worldbody> </mujoco> — Reply to this email directly, view it on GitHub <#2313 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AB2V4L5V55THTN6Q7XM4JU32IIECJAVCNFSM6AAAAABUINWXSKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDKNRWGA4TONZTGU> . You are receiving this because you were mentioned.Message ID: ***@***.***>

yuvaltassa · 2024-12-31T12:51:04Z

Reopening since we now have an MRE, thanks @n3b !

If you have capacity:

Could you verify if the issues that you're seeing are limited to float32 or if you also see them with float64?
If you've made changes locally that you think are improvements, feel free to submit a PR. Even if you are not sure that what you did is a strict improvement, this would be a good starting point for us to investigate.

Thanks again!

n3b · 2024-12-31T19:38:12Z

@emotodorov thank you, that makes so much sense now. Unfortunately my math background is way below the necessary level to make such assumptions. I never thought of precision issues in that context because I've been seeing unscaled delta cost even below -1e5.

@yuvaltassa thank you,

I tested the above example yesterday and it indeed does not produce negative improvement with doubles (there is some noise but it's very close to zero), but I can't tell for more complex models because:
I made changes indeed if one can call it that. I'm going to use Newton/CG solvers in a game engine and was evaluating their limits for a while now. So even though the solver algorithms are the same with slight differences in explicit vectorization and parallelism, the engine is different and it's rather hard to create an equal reproduction.

I do have some thoughts about improving Hessian updates besides what I already did, and about CG performance as well, but am not sure if it's the right place for it. Neither I'm confident it would make any sense from the math POV, so please let me know if you are up to discuss it.

Oh, and another thing I just remembered. The R dimensions for contacts are being overridden by the first one (using translation diag only). Is that intentional?

mujoco/src/engine/engine_core_constraint.c

Line 1466 in 69c9ac0

R[i+j] = mju_max(mjMINVAL, (1-imp)*d->efc_diagApprox[i+j]/imp);

mujoco/src/engine/engine_core_constraint.c

Line 1528 in 69c9ac0

R[i+j+1] = R[i+1]*friction[0]*friction[0]/(friction[j]*friction[j]);

mujoco/src/engine/engine_core_constraint.c

Line 1543 in 69c9ac0

R[i+j] = Rpy;

n3b added the bug Something isn't working label Dec 27, 2024

yuvaltassa closed this as completed Dec 30, 2024

yuvaltassa reopened this Dec 31, 2024

saran-t changed the title ~~CG/Newton solver terminates early due to improvement overshoots~~ [float32] CG/Newton solver terminates early due to improvement overshoots Jan 3, 2025

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[float32] CG/Newton solver terminates early due to improvement overshoots #2313

[float32] CG/Newton solver terminates early due to improvement overshoots #2313

n3b commented Dec 27, 2024

yuvaltassa commented Dec 30, 2024

yuvaltassa commented Dec 30, 2024

yuvaltassa commented Dec 30, 2024

yuvaltassa commented Dec 30, 2024

n3b commented Dec 31, 2024

emotodorov commented Dec 31, 2024 via email

emotodorov commented Dec 31, 2024 via email

yuvaltassa commented Dec 31, 2024

n3b commented Dec 31, 2024

[float32] CG/Newton solver terminates early due to improvement overshoots #2313

[float32] CG/Newton solver terminates early due to improvement overshoots #2313

Comments

n3b commented Dec 27, 2024

Intro

My setup

What's happening? What did you expect?

Steps for reproduction

Minimal model for reproduction

Code required for reproduction

Confirmations

yuvaltassa commented Dec 30, 2024

yuvaltassa commented Dec 30, 2024

yuvaltassa commented Dec 30, 2024

yuvaltassa commented Dec 30, 2024

n3b commented Dec 31, 2024

emotodorov commented Dec 31, 2024 via email

emotodorov commented Dec 31, 2024 via email

yuvaltassa commented Dec 31, 2024

n3b commented Dec 31, 2024