[Insight-developers] Problems with the vnl library

[Paul Hughett]

I've encountered several problems using the vnl library, and wanted
to make them known so other people can comment, and perhaps provide
solutions.

1. vnl_matrix_inverse does not check for singular matrices.  As far as
I can tell, this function will not provide any error signal or exception
when given a singular matrix.  In fact, trying the invert the singular
matrix [1 2; 3 6] yields the "inverse" [ 2.78842e+15 -9.29473e+14 ; 
-1.39421e+15 4.64737e+14 ] without any indication of a problem.

2. The eigenvalues computed by vnl_symmetric_eigensystem are not very
accurate.  A good symmetric eigensystem solver should produce
eigenvalues within a few eps of the correct values.  However, the
eigenvalues computed for the principal moments in the test program
itkImageMomentsTest are incorrect by about 6e-8 or 3e-7, depending on
the platform; since type double was used, these errors should be
closer to 2e-16.  In fact, the same eigenvalue computation using
octave (a Matlab clone) yields eigenvalues exact to the full machine
precision.

3. The diagonalizing transformation computed by the function
vnl_symmetric_eigensystem is improper.  That is, its determinant is -1
and it defines a mirror reflection when used as a coordinate
transformation.  This is an issue in the itkImageMoments class, where
the computation of the principal axes should yield a proper
transformation.  I am currently working around this by computing the
determinant and multiplying the last row of the transformation matrix
by it, but this is a rather ugly implementation.

Paul Hughett