Dear all,

I am a Developer of FreeCAD, in which we use Eigen in our Geometric Constraint Solver.

In the past I added support for Sparse QR decomposition of the Jacobian matrix J which has as elements the gradients that the constraints impose on the parameters of the system. We support both Full Pivot Household Dense QR decomposition and Sparse QR decomposition in FreeCAD. Given that the Jacobian is by nature very sparse, the Sparse QR decomposition has been a great success as solving times decrease with complexity hugely. I am personally indebted to your work in Eigen for such a performance improvement.

At this moment I am trying to implement a new function that shows the user the parameters that are not yet constraint (those parameters where a degree of freedom is present). This is of huge importance for mechanical engineers building complex sketches with lots of parameters and constraints, as sometimes, it gets very difficult to find where a constraint is missing in order to get to a fully constraint system.

I have approached this problem and have been successful when using the Dense Full pivoting QR decomposition (P*J*P' = QR). However, I am facing a big problem when using the Sparse QR decomposition. I do not think the problem I am facing is a bug of Eigen, but my limited understanding of Algebra and Algorithms. In any case I would appreciate if somebody with a better understanding could guide me in the problem I describe below.

Dense QR decomposition approach.

1. Let J, the matrix to be decomposed in QR be:
J=
[[ 0, 0, 1, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 1, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, -1, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, -1, 0, 0, 0],
[ 0, 0.0714286, 0, 0, 0, 0, 0,-0.0714286, 0],
[ -0.5, 0, 0, 0, 0, 0, -0.5, 0, 0],
[ 0, 1, 0, 0, 0, 0, 0, 0, -1],
[ 0, 0, 0, 0, 0, 0, 0, 1, 0]]

The columns represent the parameters of the system, which in this case they are 9 parameters. The rows represent the constraints, which in this case they are 8 constraints. The rank of this matrix is 8.

2. This is decomposed into a QR matrix, with the R matrix (the top triangular part of the trapezoidal R matrix to be accurate) being:
R =
[[ 1, 0, 0, 0, 0, 0, 0, 0],
[ 0, 1, 0, 0, 0, 0, 0, 0],
[ 0, 0, -1, 0, 0, 0, 0, 0],
[ 0, 0, 0, -1, 0, 0, 0, 0],
[ 0, 0, 0, 0, -1.41421, 0, 0,-0.0505076],
[ 0, 0, 0, 0, 0, 1, 0,-0.0714286],
[ 0, 0, 0, 0, 0, 0, 0.707107, 0],
[ 0, 0, 0, 0, 0, 0, 0, 0.0505076]]

I appreciate that the diagonal shows decreasing values.

3. In this case, there is a "rowsTransposition" member function that provides the row transpositions. With this I build a permutation matrix which enables me to map the independent parameters (8 out of 9 parameters in this case) corresponding to the diagonal of R to the parameters of the matrix J. I can then show to the user the dependent parameters by highlighting the geometric elements associated to them. The code is this:

Code: Select all

        Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic> rowPermutations;

        rowPermutations.setIdentity(paramsNum);

        if(qrAlgorithm==EigenDenseQR){ // P.J.P' = Q.R see https://eigen.tuxfamily.org/dox/classEigen_1_1FullPivHouseholderQR.html
            const MatrixIndexType rowTranspositions = qrJT.rowsTranspositions();

            for(int k = 0; k < rank; ++k)
                rowPermutations.applyTranspositionOnTheRight(k, rowTranspositions.coeff(k));
        }

        // params (in the order of J) shown as independent from QR
        std::set<int> indepParamCols;
        std::set<int> depParamCols;

        for (int j=0; j < rank; j++) {
             int origRow = rowPermutations.indices()[j];
             indepParamCols.insert(origRow);.
        }

        // If not independent, must be dependent
        for(int j=0; j < paramsNum; j++) {
            auto result = std::find(indepParamCols.begin(), indepParamCols.end(), j);
            if(result == indepParamCols.end()) {
                depParamCols.insert(j);
            }
        }


This detects the parameter of J matrix column index 0, this is the first column, as being dependent, which is fine as the first and seventh cols are the same in the J matrix above:

Code: Select all

 RowTransp =
[
[[2,3,4,5,5,7,6,8]]
]
R row 0 = J col 2
R row 1 = J col 3
R row 2 = J col 4
R row 3 = J col 5
R row 4 = J col 1
R row 5 = J col 7
R row 6 = J col 6
R row 7 = J col 8
Indep params: [1,2,3,4,5,6,7,8,]
Dep params: [0,]


This works ok in my tests. If you are curious about how this look in FreeCAD, you may check this post:
https://forum.freecadweb.org/viewtopic. ... 00#p214852

Sparse QR decomposition approach.

1. The Sparse QR decomposition, same J matrix, provides this R matrix:
R =
[[ -1, 0, 0, 0, 0, 0, 0, 0],
[ 0, -1, 0, 0, 0, 0, 0, 0],
[ 0, 0, -1, 0, 0, 0, 0, 0],
[ 0, 0, 0, -1, 0, 0, 0, 0],
[ 0, 0, 0, 0, -0.707107, 0, 0, 0],
[ 0, 0, 0, 0, 0, 1.41421, 0.0505076, 0],
[ 0, 0, 0, 0, 0, 0,-0.0874818, 0.816497],
[ 0, 0, 0, 0, 0, 0, 0, -0.57735]]

Here I appreciate that the diagonal values of the R matrix (top triangular of the trapezoidal matrix), are no longer ordered in decreasing order.

2. The SparseQR decomposition, does not have a rows transposition member function, which could enable me to map the independent parameters of R into the parameters in J.

I have done quite a reading on the Internet. There is the concept of "left-looking rank-revealing" and "right-looking rank-revealing", which I do not fully understand and I am not sure if it bears any significance to it. It appears to me from what I have read that "left-looking" is generally computationally advantageous over "right-looking", and from the decomposition (A.P=QR), intuitively I think that it may be related to the fact that there is not a row left-multiplication permutation matrix on the left.

However I fail to understand whether there is something I can do on the result of the Sparse QR decomposition to get this "rowsTransposition" matrix or not.

I apologise for this very long post.

Thanks in advance,
abdullah