ITK/Release 4/Refactor Numerical Libraries/Inventory/Linear Solvers
It certain ITK applications it is common to need to solve the linear system
A.x = b
- 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.
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.
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:
The initial approach in ITK to solve large sparse linear systems involved the use of the LSQR method.
LSQR is implemented in VXL:
and tested here
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:
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
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.
Unit test for the FEM Registration Filter
As an Image registration problem performing Deformable registration
The two examples