@whitequark
I see allowed redundant constraints here, so I've found what 12.7 pt-pt-distances appied twice like as two pt-on-circle. If we are delete all the redundant, we can converge here fine.

@whitequark, Somehow too many redundant constraints can't be converged. I've figured out what 3xPT_ON_CIRCLE can't converge (just try to make it on empty sketch). After disallowing redundant, too many redundants threated as non-converged for no reason. This is probably bug, but I am not sure whether rank test can work fine or not for such amount of redundants.

@Evil-Spirit Yes, I did notice that many redundant constraints tend to degrade convergence. Please figure out why. This makes such a mode much less useful.

@whitequark, I have researched a lot and made some debug functions (visualizing solve way and equations). Finally I've linked Eigen and this has solved the problem. I think the problem is what the Gauss linear system solving is not so precise as it should be.

here is some beta
Evil-Spirit@703bccd

1. What is "solve way"?
2. Are there any downsides to using Eigen?

@whitequark, this is path of points during newton steps like here https://youtu.be/72Dnqvn_uTk

OK. Please prepare a commit that just adds Eigen, and separately these debugging capabilities.

2. I can try to work with Gauss Method to find any ways to improve... But I think it is better to move in direction of Eigen integration because of better stability and convergence and a lots of methods to choose.

@Evil-Spirit I see no problem with Eigen. Just want to know if there is anything it's bad at.

@whitequark
here is eigen https://github.com/Evil-Spirit/solvespace-master/commits/eigen
here is debug https://github.com/Evil-Spirit/solvespace-master/commits/solver-debug

@whitequark, what about to remove System::MAX_UNKNOWNS limit in eigen mode?

@Evil-Spirit Yeah, just rip out all pre-Eigen solver code, no #ifdefs needed, as well as MAX_UNKNOWNS limit. No reason to keep these.

@whitequark ok

@whitequark updated https://github.com/Evil-Spirit/solvespace-master/commits/eigen

@whitequark
the results is:
constarintfreeze.slvs from issue #131
With substitutions:

Equations: 136, Unknows: 146
Jacobian non zeroes: 617/19856 3.1074%
WriteEquations: 1.047 ms
WriteJacobian:  15.663 ms
EvalJacobian:   1.089 ms
Substitution:   36.903 ms
LinearSystem:   29.818 ms
CalculateRank:  36.674 ms
Generate::DIRTY (for bounding box) took 134 ms
Generate::DIRTY took 138 ms

Without substitutions:

Equations: 654, Unknows: 664
Jacobian non zeroes: 1911/434256 0.44006%
WriteEquations: 1.037 ms
WriteJacobian:  70.673 ms
EvalJacobian:   2.617 ms
Substitution:   0.000 ms
LinearSystem:   2795.340 ms
CalculateRank:  2955.157 ms
Generate::ALL (for bounding box) took 5841 ms
Generate::ALL took 5846 ms

@whitequark, @jwesthues

1. We have o(m*n)~o(k^2) (m = number of equations, n = number of unknowns) for SolveBySubstitions and WriteJacobian and now I am working on improving this.
2. The sparsness of matrices is high. I think we have good memory usage reducing, but the speed of linear eq solving (and rank calculating) is 2x slower when it was. I think I will look into finding some better solver than QR decomposion inside eigen libarary. @jwesthues do you have any methods suggestions? SVD? LU? I don't know what does this acronyms means, I am not so good in math.

The most robust way to rank-test is usually an SVD. Take care in how you set the tolerance on the rank test. That's a count of singular values above some threshold; so perhaps set that threshold empirically by looking at the singular values for a variety of sketches with known number of redundant DOF, and putting it halfway (geometric mean) between the biggest "zero" and the smallest "nonzero". Make sure to consider ill-conditioned sketches, like with long skinny triangles and stuff. If that SVD costs significant time, then we could look at a faster rank test.

Then what's the condition number of the resulting matrix, the ratio of the maximum singular value to the minimum non-zero singular value? If it's not too big, then it shouldn't matter much how you solve; so you can probably just do Cholesky or something like before, faster than QR.

@whitequark , @jwesthues
I've found SimplicalLDLT. It working for symmetrical matrices (A * A^T is symmtrical).

Equations: 1280, Unknows: 1552 
Jacobian non zeroes: 7680/1986560 0.3866% 
A*A^T non zeroes: 12880/1638400 0.78613% 
Newton Steps: 5 
WriteEquations: 15.301 ms 
WriteJacobian: 531.120 ms 
EvalJacobian: 47.934 ms 
Substitution: 24.923 ms 
LinearSystem: 65.149 ms

And this is ultrafast. @jwesthues what do you think about https://eigen.tuxfamily.org/dox/classEigen_1_1SimplicialLDLT.html?

@Evil-Spirit Cool. What's the condition number of the matrix? It should intuitively be not too big, in which case Cholesky(-ish) should work robustly, though I don't think I ever actually measured...

@Evil-Spirit Excellent work!

@jwesthues , I will try to calculate the condition number of the matrix and do output for it.
@whitequark Now we have WriteJacobian: 531.120 ms as bottleneck, but this is probably just Expr massive allocations (just need to make faster allocator). Anyway, this is not the problem. The next real problem is rank test. It's slow for rank-revealing QR, but @jwesthues suggesing to try SVD for it... But eigen probably have no any SparseSVD routines.

SVD rank testing is better but unlikely to be faster. I'm not sure about Eigen, but sparse SVD solvers are well-studied, http://www.netlib.org/svdpack/ or whatever.

Iterative methods for the SVD may have a special benefit here if we cache the result, since while dragging, the A matrix doesn't usually change that much from solve to solve. So the SVD shouldn't change that much either, so we might usually be almost converged to start with.

@jwesthues, I have calculated using approach mentioned here http://stackoverflow.com/questions/33575478/how-can-you-find-the-condition-number-in-eigen

JacobiSVD<MatrixXd> svd(A);
double cond = svd.singularValues()(0) 
    / svd.singularValues()(svd.singularValues().size()-1);

the result is 1e+7...1e+14 for A*A^T matrix during newton steps (for a long solution solved in 39 steps )
what does it means?

It's the number of digits of precision that we lose while solving; but we need to define the condition number as the ratio of the maximum singular value to the minimum non-zero singular value. I suspect that your condition number is artificially big because your sketch has unconstrained DOF, and each unconstrained DOF introduces a singular value that's ideally zero.

Ignore all the singular values less than a threshold (e.g., 10^-4 to start), and try to put the threshold roughly halfway (geometric mean) between the smallest singular value that should be non-zero, and the largest singular value that should be zero, over a range of difficult sketches like long skinny triangles. That's probably the most robust way to rank test.

@jwesthues example of singular values list is

sing v = 20558
sing v = 20435
sing v = 20320
sing v = 19977
sing v = 19764
sing v = 19465
sing v = 19080
sing v = 18356
sing v = 17790
sing v = 17736
sing v = 17226
sing v = 16830
sing v = 16279
sing v = 15547
sing v = 150
sing v = 20435
sing v = 20320
sing v = 19977
sing v = 19764
sing v = 19465
sing v = 19080
sing v = 18356
sing v = 17790
sing v = 17736
sing v = 17226
sing v = 16830
sing v = 16279
sing v = 15547
sing v = 15058
sing v = 14446
sing v = 14294
sing v = 13547
sing v = 13190
sing v = 12636
sing v = 12367
sing v = 11415
sing v = 10903
sing v = 10867
sing v = 10418
sing v = 9713.7
sing v = 9456.6
sing v = 9195.2
sing v = 8587.7
sing v = 8362.7
sing v = 8300
sing v = 8253.4
sing v = 6595.2
sing v = 6573
sing v = 6548.5
sing v = 6502.8
sing v = 6142.6
sing v = 6043
sing v = 5838.2
sing v = 5677.5
sing v = 5634.2
sing v = 5488.1
sing v = 5389.9
sing v = 5042.3
sing v = 4844
sing v = 4741.4
sing v = 4661.2
sing v = 4487.3
sing v = 4331.4
sing v = 4279.9
sing v = 4054.2
sing v = 3882.6
sing v = 3556.7
sing v = 3408.2
sing v = 3372.6
sing v = 3255.1
sing v = 3182.6
sing v = 2949.9
sing v = 2821.2
sing v = 2717.1
sing v = 2673.2
sing v = 2569.9
sing v = 2530.4
sing v = 2472.4
sing v = 1.4846
sing v = 1.4643
sing v = 1.4389
sing v = 1.4205
sing v = 1.3988
sing v = 1.3901
sing v = 1.3845
sing v = 1.3757
sing v = 1.3656
sing v = 1.3476
sing v = 1.3302
sing v = 1.3161
sing v = 1.2965
sing v = 1.2926
sing v = 1.2553
sing v = 1.2079
sing v = 1.1996
sing v = 1.1813
sing v = 1.1634
sing v = 1.1313
sing v = 1.1077
sing v = 1.0826
sing v = 1.0512
sing v = 1.0342
sing v = 1.0145
sing v = 0.98174
sing v = 0.94912
sing v = 0.92502
sing v = 0.88509
sing v = 0.85739
sing v = 0.83692
sing v = 0.81474
sing v = 0.78855
sing v = 0.77264
sing v = 0.75487
sing v = 0.73647
sing v = 0.71384
sing v = 0.68019
sing v = 0.66856
sing v = 0.65449
sing v = 0.6515
sing v = 0.63048
sing v = 0.61453
sing v = 0.60073
sing v = 0.58435
sing v = 0.56654
sing v = 0.56132
sing v = 0.54243
sing v = 0.5393
sing v = 0.52999
sing v = 0.51901
sing v = 0.51483
sing v = 0.50115
sing v = 0.49423
sing v = 0.48999
sing v = 0.45135
sing v = 0.43907
sing v = 0.40643
sing v = 0.40353
sing v = 0.38339
sing v = 0.31818
sing v = 0.31435
sing v = 0.31141
sing v = 0.24033
sing v = 0.0014815
sing v = 0.0014584
sing v = 0.0014278
sing v = 0.0013512
sing v = 0.0013218
sing v = 0.0012358
sing v = 0.00121
sing v = 0.0010954
sing v = 0.001019
sing v = 0.0009437
sing v = 0.00090831
sing v = 0.00075787
sing v = 0.00073816
sing v = 0.00070292
sing v = 0.00069507
sing v = 0.00068166
sing v = 0.0006793
sing v = 0.00065714
sing v = 0.00065452
sing v = 0.000614
sing v = 0.00054296
sing v = 0.00051956
sing v = 0.0004122
sing v = 0.00037248
sing v = 0.0003227
sing v = 0.00029575
sing v = 0.00026793
sing v = 0.00019084
sing v = 0.00014426
sing v = 0.0001031
sing v = 8.2979e-005
sing v = 3.1564e-005

Okay, got it. So you've got 32*4*2 rectangle vertex coordinates, minus 31*2 unknowns for the point-coincident at the corner, minus 32*2 for the lengths, minus 32*2 for the parallels, minus 32 for the perpendiculars. That's our 34 DOF.

You write 32*2 + 32*2 + 32 = 160 nontrivial equations, and the A matrix would ideally have full rank. So this isn't a useful test case for what I was talking about, but it's still interesting. We have 159 singular values, for a reason that I don't understand; was a line lost to cut and paste? The reported condition number is slightly smaller than 20435/3.1564e-005, so I'm thinking you lost the first line.

In any case, the singular values fall into three groups, 32 around 10^-4, 64 around 10^0, and 63 (that I think should be 64) around 10^4. I think that corresponds to your three types of constraints. That indicates that we're writing those equations in an unfortunate way--like, if you figure out what kind of equation is giving you the ~10^-4, and multiply all those equations by 10^4, then you'll improve the matrix's condition number. That's easy enough to do experimentally, and then look at the structure of the equations to see if it makes sense.

What are the singular values in some sketches that are actually inconsistent/redundant?

@jwesthues thank you, you are swimming at the deep where I still can understand you. Now I understand a little bit more about how to answer the question what equation is good, and what is bad. You can read this article I have written.

PARALLEL

The restriction[constraint] of parallelism is also based on a pseudoscalar (oblique[skew]) product. Here we use the fact that sin (0) = 0. To be honest, we would have to divide the left side of the equation into[by the] the product of the moduli of vectors, but historically in SolveSpace the equation is written just like this. Since to the zero result of the oblique[skew] product [leads also the case where] the length of the vectors [equals] to zero, the segments tend to decrease all the time. This does not seriously interfere[meaningful, annoying] if the length of the vectors is limited[using length constraints].

I should check the PARALLEL equation because this is not scaled by vectors magnitudes, and scaling it can improve the situation. But there still a question how to make length and angular dimension of the same scale. I have mentioned equations in projective space, so this is also can impove situation.

@jwesthues
As I suspect, the 3k...20k singular values was introduced by parallel constraints. Turning this equation in normalized form removes this big values. But rotation and length still in the different scale.

@Evil-Spirit Cool. One option is just to multiply or divide by a constant scale factor derived from the bounding box of the overall sketch. That's ugly, but it does the right thing as long as the initial guess is drawn with halfway reasonable (within a few orders of magnitude) scale.

What singular values do you get when constraints really are redundant/inconsistent?

@whitequark , ready here
https://github.com/Evil-Spirit/solvespace-master/tree/eigen-master

@whitequark,
https://github.com/Evil-Spirit/solvespace-master/tree/solver-tests

@Evil-Spirit any chance you could rebase these changes on current master? You and JW put a lot of time into it and I'd like to get it merged. The solver stuff is tricky, and I'm reluctant to just make the changes fit in a way that looks right, compiles and seems to work. If you're busy I understand.

@phkahler Sorry, have no time :(

@phkahler made it

Impossible not to put this here
https://youtu.be/-RdOwhmqP5s
for future reference :-)

Merged in #1159