ITK/Release 4/Refactor Numerical Libraries/Inventory/Linear Solvers

The Problem

It certain ITK applications it is common to need to solve the linear system

           A.x = b

Where

• A is a large sparse matrix
• x is a vector
• b is a vector

The common approach to solve this problem is to use linear solvers.

Use Cases

FEM

The most significant case where large sparse matrix appear is in the computation of Finite Element Methods (FEM). In ITK, this is done for the purpose of computing Deformable Registration between images.

Exploration

In the past, ITK developers have explored alternative solutions for solving large sparse matrices.

The findings were codified at the time in the following Wiki page:

• Proposals:Sparse_Linear_Solvers

LSQR

The initial approach in ITK to solve large sparse linear systems involved the use of the LSQR method.

• http://www.stanford.edu/group/SOL/software/lsqr.html

In VXL

LSQR is implemented in VXL:

  ITK/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_lsqr.h
  ITK/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_lsqr.cxx

and tested here

  ITK/Modules/ThirdParty/VNL/src/vxl/core/vnl/examples/vnl_lsqr_test.cxx

Unfortunately the VXL implementation is done using source code taken from netlib, and that code happens to be copyrighted by ACM. The license that ACM provides for algorithms preclude their use in the ITK ecosystem.

 //
 // vnl_lsqr implements an algorithm for large, sparse linear systems and
 // sparse, linear least squares. It is a wrapper for the LSQR algorithm
 // of Paige and Saunders (ACM TOMS 583). The sparse system is encapsulated
 // by a vnl_linear_system.
 //
 // The LSQR code from netlib was copyrighted by ACM, therefore it had to 
 // be removed, and has been replaced with code taken from 
 // http://www.stanford.edu/group/SOL/software.html
 // Thas was distributed under a BSD license. The Fortran90 files were
 // manually translated to C++ to create the solverBase class.
 //

The solverBase class can be found in:

  ThirdParty/VNL/src/vxl/v3p/netlib/linalg/lsqrBase.h
  ThirdParty/VNL/src/vxl/v3p/netlib/linalg/lsqrBase.cxx

This new code is cleanly written under a BSD license.

At some point we must send this code back to Professor Saunders along with the thanks for helping us navigate the problem by pointing us to the newer BSD versions of the FORTRAN code.

Use From ITK

The lsqr code is used in ITK in the classes

  Numerics/FEM/include/itkFEMLinearSystemWrapperVNL.h
  Numerics/FEM/src/itkFEMLinearSystemWrapperVNL.cxx

Uncertainty exist on whether the ITK code is properly using the LSQR Solver.

In particular, ITK code does not make a particular effort on

• Verifying the pre-conditions of the system,
• Setting the many parameters of the LSQR solver to anything different than the default values

It may be that a quick Wiki-based tutorial on the practical ways of properly using LSQR is the solution to this problem.

It may also be that multiple predefined settings could be provided, in order to address specific numerical problems.

Use Cases

Unit test

Unit test for the FEM Registration Filter

 /home/ibanez/src/ITK/Modules/Registration/FEM/test/itkFEMRegistrationFilterTest.cxx

Examples

As an Image registration problem performing Deformable registration

 ITK/Examples/Registration/

The two examples

• DeformableRegistration1.cxx
• DeformableRegistration11.cxx