As an exploration, I disabled the check for a singular matrix in
Matrix::fsolve. Attached is a small patch which just comments out the check on the reciprocal condition number.
disable_rcond.diff (778 Bytes)
MatrixType code identifies
A as a singular matrix, but that can be overridden with the
matrix_type function. With the patch above, I get an LU decomposition and the “correct” answer.
AFull = matrix_type (A, 'full');
Q = AFull \ B
The code uses a test for singularity which is
rcond + 1.0 == 1.0
This will occur if
rcond is smaller than
eps (1.0). Alternatively, this will occur for a condition number >
1 / (2*eps) or 2.25e15. Given that the original reporter is working with a matrix that has a condition number of 2e73, this will clearly be tagged as singular.
Maybe we should do what Matlab does. The wording of the warning is
“Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.275234e-067.”
This correctly says that trouble may result, but it goes ahead and treats the matrix as an ordinary matrix and attempts an LU decomposition.
Octave, on the other hand, is definitive that this is a singular matrix:
“matrix singular to machine precision, rcond = 8.27523e-67”
It then goes on to choose an entirely different algorithm (minimum norm) to solve the problem.
If we wanted too, we could use the
INFO output from
DGETRF to detect truly singular matrices, but otherwise process matrices with large condition numbers.
INFO (output) INTEGER
= 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal
> 0: if INFO = i, U(i,i) is exactly zero. The fac-
torization has been completed, but the factor U is
exactly singular, and division by zero will occur if
it is used to solve a system of equations.
Of course, this would be a huge paradigm change for Octave and would need some robust debate about the pros and cons of such a move.