Hi!
I'm working on a 3-D surface mesh and trying to obtain gaussian curvature values by fitting a quadric polynomial to each vertex neighborhood. In order to obtain the polynomial coefficients I'm using the standard cholesky decomposition (LL^T), for solving my system Ax = b, where A is a symmetric, positive definite matrix.

When solving the system I'm getting the strangest values only for one of the solution (among about 2000 vertices!). Here's a sample of the what would be the inverse of my symmetric, positive definite matrix:

Vertex 85:
288.413 6.45563 3.95006 -48.7131 -17.8428 -30.1609
6.45563 0.199845 0.103744 -1.09986 -0.409136 -0.673235
3.95006 0.103744 0.0979785 -0.666895 -0.241934 -0.41527
-48.7131 -1.09986 -0.666895 8.24633 3.03219 5.09083
-17.8428 -0.409136 -0.241934 3.03219 1.14272 1.8648
-30.1609 -0.673235 -0.41527 5.09083 1.8648 3.16057
Vertex 86:
496411 2072.94 -8090.48 -48680.5 23267.4 -88125.4
2072.94 8.69566 -33.8001 -203.285 97.1605 -367.998
-8090.48 -33.8001 131.916 793.391 -379.213 1436.27
-48680.5 -203.285 793.391 4773.85 -2281.71 8641.99
23267.4 97.1605 -379.213 -2281.71 1090.6 -4130.56
-88125.4 -367.998 1436.27 8641.99 -4130.56 15644.5
Vertex 87:
523.131 -1.91535 -7.666 -49.6631 17.1292 -96.191
-1.91535 0.0430769 0.0135821 0.181351 -0.061675 0.352979
-7.666 0.0135821 0.174082 0.72835 -0.249803 1.41043
-49.6631 0.181351 0.72835 4.72212 -1.62466 9.12386
17.1292 -0.061675 -0.249803 -1.62466 0.58246 -3.16276
-96.191 0.352979 1.41043 9.12386 -3.16276 17.7163

Just by comparing the first element of each matrix you can see the differences: 288 (vertex 85), 496411 (vertex 86) and back to normal values on vertex 87: 523.131...

Here's part of the code I'm using (B is my symmetric, positive definite matrix):
MatrixXd B = A.adjoint()*A;
MatrixXd x = B.llt().solve(A.adjoint()*b)

I don't know if I'm missing something... How to get rid of these unexpected results?
Thanks,
Miguel

No idea. Perhaps you cold try with .ldlt() and .lu() to see if you get a different answer. If not you could also try to avoid squaring:

x = A.jacobiSvd(ComputeThinU|ComputeThinV).solve(b);

If you still get the same result, then blame your matrix A or right hand side b!

btw, if you want to handle larger meshes you might consider using a SparseMatrix with a sparse solver.

Thanks for your reply...

Unfortunately I tried your suggestions (.ldlt, .lu and .jacobisvd) but in all three cases I get the exact same result for the solution (vector x), it's the strangest thing since my matrices A and b do not look that different:

- For matrix A


- For matrix b


I've also solved the system using Matlab, and the results are exactly the same!!! So I'm sure is not an Eigen or Matlab issue... For sure, as you said, it must be something weird with my matrices A and/or b, but I cannot see it!

Hope anybody can give me a clue!
Thanks,
Miguel

msotaquira wrote:Thanks for your reply...

Unfortunately I tried your suggestions (.ldlt, .lu and .jacobisvd) but in all three cases I get the exact same result for the solution (vector x), it's the strangest thing since my matrices A and b do not look that different:


I've also solved the system using Matlab, and the results are exactly the same!!! So I'm sure is not an Eigen or Matlab issue... For sure, as you said, it must be something weird with my matrices A and/or b, but I cannot see it!

Hope anybody can give me a clue!

Check the condition number of the matrix. Save the result of

Code: Select all

A.jacobiSvd(ComputeThinU|ComputeThinV)

and apply the .singularValues() method then form the ratio of the last to the first (which is actually the reciprocal condition number but it avoids problems with division by zero in the case of a computationally singular matrix). If the reciprocal condition number is close to zero then you have an ill-conditioned matrix and calculation of the inverse is unstable.

dmbates wrote...
Check the condition number of the matrix. Save the result of

Code: Select all

A.jacobiSvd(ComputeThinU|ComputeThinV)

and apply the .singularValues() method then form the ratio of the last to the first (which is actually the reciprocal condition number but it avoids problems with division by zero in the case of a computationally singular matrix). If the reciprocal condition number is close to zero then you have an ill-conditioned matrix and calculation of the inverse is unstable

Hi dmbates, I tried your suggestion and regarding the reciprocal condition here's what I found:
[code]- vertex 85: 3.79459153 x 10-3
- vertex 86: 8.99388902 × 10-5
- vertex 87: 2.72287945 x 10-3[\code]

so the reciprocal condition for vertex 86 is two orders of magnitude less than the other two vertices! Finally I can see that there is something wrong with this matrix!

In my case I'm obtaining matrix A by projecting the coordinates of vertex's neighbors into a local coordinate system centered at current vertex. For example: if current vertex is 86, then I find three orthonormal vectors at this vertex and then transform each neighbor (vertices 85, 83, 87, etc.) to this local system. Elements of A are formed by the modified x,y,z coordinates of each of these transformed locations.

Now the big question: how can I avoid having this ill-conditioned matrix? Is it possible to avoid it?

Any suggestions are welcome!