I'm at my wits end with lscm

[AdamSturge]

I keep getting this error when I run igl::eigs no matter what I do. I'm not sure what counts as too much information for an issue but I've been debugging this for 4h now so I've reached the breaking point.

Assertion failed: a_lhs.cols() == a_rhs.rows() && "invalid matrix product" && "if you wanted a coeff-wise or a dot product use the respective explicit functions", file c:\uot-masters\geometry-processing\geometry-processing-parameterization\shared\libigl\external\nanogui\ext\eigen\eigen\src/Core/ProductBase.h, line 102

I've even removed my vector_area_matrix code and replaced it with igl::vector_area_matrix.

Here's what I do: 1.) Use igl::cotmatrix to compute L 2.) Use igl::massmatrix to compute M (with igl::MASSMATRIX_TYPE_DEFAULT) 3.) Replicate L twice on the diagonal using igl::repdiag 4.) Replicate M twice on the diagonal using igl::repdiag 5.) Compute A using igl::vector_area_matrix 6.) Compute Q as it's defined in the ReadMe 7.) Run igl::eigs with Q and B (using igl::EIGS_TYPE_SM)

If I only ask for 1 eigen pair it runs fine, but more than 1 results in the above error message.

EDIT: I'm on windows

[M-Anwar]

I am also getting this error in my LSCM method, on both windows and linux. I am pretty much doing what you are doing, constructing the matrices as shown in the readme, and sending them to the igl::eigs method.

I tried using the built-in libigl lscm() function using the tutorial 502 code on the camel-head.off example. The libigl lscm method uses a different method to minimize their constrained energy, but that seems to work. Not sure why this method doesn't.

Anyone got any ideas?

[M-Anwar]

Okay I got it working, after spending...3+ hours setting up breakpoints. The problem was not in my code, as far as I could tell. The problem is in eigs.cpp.

On line 140 there is a statement (S.transpose()*B*x).array().abs()>=1e-7).all(). But the dimensions of this matrix operation don't match up S.transpose() is a 1xk matrix but B*x is a nx1 matrix. So this is where the assert kept failing. I also looked up the commits to this file on libigl, and this line was actually added in by Alec last week. I just used an older version of this file, which removes that part of the if-statement.

Running my code on the beetle.obj example, the solver "restarts" a few times, but eventually outputs the correct parameterization. The solver doesn't converge on any of the other examples however. But I think my implementation is correct, so I think I should move on.

If anyone else finds a different solution to this, it would be nice to see.

[ghost]

Yeah, I'm having a lot of trouble getting lscm to run as well. The main problem I get is the "Failed to Converge" error.

[AdamSturge]

I had a look at Eigen and it does have https://eigen.tuxfamily.org/dox-devel/classEigen_1_1GeneralizedEigenSolver.html and https://eigen.tuxfamily.org/dox-devel/classEigen_1_1GeneralizedSelfAdjointEigenSolver.html

but they don't seem to compile for SparseMatrix

[ghost]

I tried Muhammad's solution, and I got to the point where at least it's saying "Restarting," but it's been doing that for so long, I'm not actually sure it solved my issue the way it did for him.

[AdamSturge]

Making the replacement he suggested causes my code to restart a bunch of times then eventually converge for the beetle. I haven't tested it on other meshes. Have you tried using igl::vector_area_matrix() instead of your own code to confirm that it is igl::eigs()

[ghost]

Yeah, that's what I tried too. I'll try again, and make sure the other parts are correctly coded.

[AdamSturge]

After trying out some other meshes it seems that the change @M-Anwar made only works for the beetle, as he said himself

[AdamSturge]

So I think I understand what the line @alecjacobson added is trying to deal with.

if
(
    i==0 || 
    (S.head(i).array()-sigma).abs().maxCoeff()>1e-14 ||
    ((S.transpose()*B*x).array().abs()>=1e-7).all() <-- NEW LINE
)
{

S contains the eigenvalues (sigma) built up so far. The first condition says "if we only have found 1 eigenvalue, then store it". The second line says "if this eigenvalue is sufficiently different from the other eigenvalues we stored, then store it". The commit comment for this new line says "eigs now working even when there are modes with multiplicity". So I believe this line is meant to deal with the case when there are degenerate eigenvectors (eigenvectors that have the same eigenvalue). This would be the case in our assignment as there are 2 eigenvectors that map to 0.

The new line sort of looks like the rhs of $Av = \lambda Bv$, which is a way to formulate the generalized eigenvalue problem

To that end I made the following change

if
(
    i == 0 ||
    (S.head(i).array() - sigma).abs().maxCoeff() > 1e-14 ||
    ((U - x).colwise().norm().array() >= 1e-7).all() <-- MY CHANGE
)
{

All it does is check if the new eigenvector is sufficently different from the eigenvectors found so far. I'm not sure if this makes sense to do in general but it seems to work for the assignment

[alecjacobson]

I believe I have located and fixed a bug in igl::eigs. If you git pull and then git submodule update --recursive you should see a fix.

Awkwardly, this bug was not showing up in Release mode for me and I hadn't tested in Debug mode.

[alecjacobson]

Ack. Maybe wait one sec, there's a git submodule mess now...

[alecjacobson]

Now it's working.

[M-Anwar]

Okay, so I am using the updated eigs function, and my parametrizations turn out quite weird (for the beetle below).

I found that by turning A into a symmetric matrix, by doing A + A^T, fixes this problem. This may be a numerical issue, with the igl::eigs method, but I am not sure.

Any way it seems to be working now. Thanks for the help!

[alecjacobson]

yes, A must be symmetric. This was updated in the README (and we discussed it in class last week).