Reimplementing constraint solver

[DeepSOIC]

Testing if the scaling works. It does!
Running this

Code: Select all

sc = 1e4
import ConstraintSolver as CS
ps = CS.ParameterStore()
p1 = CS.G2D.ParaPoint(ps)
p2 = CS.G2D.ParaPoint(ps)
p2.x.Value = 3 * sc
p2.y.Value = 4 * sc
p1.update()
p2.update()
for p in p1.Parameters + p2.Parameters:
    p.OwnScale = sc

c = CS.G2D.ConstraintDistance(
    p1 = p1,
    p2 = p2,
    store = ps,
    Label = "Constraint1"
)
c.dist.Value = sc
c.Weight = 1.0/sc

c.update()
c.NetError

sys = CS.SubSystem()
sys.addUnknown(p1.Parameters)
sys.addUnknown(p2.Parameters)
sys.addConstraint(c)

vs = CS.ValueSet(CS.ParameterSubset(ps.allFree()))
slv = CS.SketchSolver()
slv.solveDogLeg(sys, vs)

, setting different values to sc doesn't change solver output at all.
Which means, the scaling is working as expected. So hopes are high, those faults related to meter and larger sized sketches can be finally solved. Now, fingers crossed, I ever get that far.

[DeepSOIC]

Now onto LM solver, I immediately stumbled on the reason it has accuracy problems. Its termination condition is:
e1^2 + e2^2 + .. < eps, where default eps = 1e-10
In dogleg, it is:
max(abs(e1), abs(e2), ...) < tolf , where tolf default is again 1e-10
Since DL is comparing linearly and LM is comparing squared, there is a 5 order of magnitude difference in accuracy.

[toralf]

Great progress. Sound interesting. Could you make use of the dual numbers?

[DeepSOIC]

Could you make use of the dual numbers?

I am, see first post of the thread

[toralf]

I know that you wrote it in the first post, but these were your goals. Did they prove to be useful for this reimplementation?

[DeepSOIC]

They already proved useful in sketcher constraints, there is just no question about it. In the implementation so far, they aren't very useful, because I have not written geometry and constraints yet. The two I wrote are too trivial to say "dual numbers were useful", but for more complex stuff, dual numbers are just a lifesaver IMO.

[DeepSOIC]

LM solver is up!
Here's a log where it solves the same super-simple test problem:

Code: Select all

Begin Levenberg-Marquardt solving
  iteration 0
    err = 2.000000
    max(abs(gradient[i])) = 1.600000
    damping factor = 0.000640, computing step...
    sq(step) = 1.998721
    step changed the error by 4.000000 (expected: 4.000000)
    step accepted
    changing damping by factor 0.333300
  iteration 1
    err = 0.000640
    max(abs(gradient[i])) = 0.000512
    damping factor = 0.000213, computing step...
    sq(step) = 0.000000
    step changed the error by 0.000000 (expected: 0.000000)
    step accepted
    changing damping by factor 0.333300
  iteration 2
    err = 0.000000
    max(abs(gradient[i])) = 0.000000
    damping factor = 0.000071, computing step...
    sq(step) = 0.000000
    step changed the error by 0.000000 (expected: 0.000000)
    step accepted
    changing damping by factor 0.333300
  iteration 3
    err = 0.000000
error is small, solved

[ickby]

Nice progress!

It seems there is an error in the implementation.
double dL = err - 0.5*(fx + Jx*h_dl).squaredNorm();
double dF = err - err_new;
double rho = dL/dF;

Hm, back in the days when I implemented dogleg and LM solvers I closely followed a paper (or excerpt of a bock, not sure anymore). I will have a look this evening if I still can find it on my other machine. There is a good chance I switched some variables and made a error, as I experimented a lot with those solvers.

[DeepSOIC]

I closely followed a paper (or excerpt of a bock, not sure anymore). I will have a look this evening if I still can find it on my other machine.

That'd be very nice to have, as I'm struggling to find a good source for dogleg for example. LM is decently explained on Wikipedia (although I'm still struggling to wrap my head around that J^T * J thing).
If I just drop all that "understand the solver algo" part, I feel like I can copy-paste and translate the implementation to the new system a lot faster. But it's interesting , and potentially useful for the future.

[DeepSOIC]

LM's damping factor changer function

Code: Select all

                double dF = e.squaredNorm() - e_new.squaredNorm();
                double dL = h.dot(mu*h+g);

                if (dF>0. && dL>0.) { // reduction in error, increment is accepted
                    double tmp=2*dF/dL-1.;
                    mu *= std::max(1./3., 1.-tmp*tmp*tmp); //   <- this multiplier is plotted