qr_solve(): unexpected failure on 4x4 system #983
Description
Attempting to solve the following system of equations using qr_solve(A, b) results in an error, although it appears the solution should exist:
Reproducing code example, cpython version = 3.12.1 mpmath version = 1.3.0:
from mpmath import mp
A = mp.matrix([[mp.one, -mp.pi / 20, (-mp.pi / 20) ** 2, (-mp.pi / 20) ** 3],
[mp.one, mp.zero, mp.zero, mp.zero],
[mp.one, mp.pi / 20, (mp.pi / 20) ** 2, (mp.pi / 20) ** 3],
[mp.one, mp.pi / 10, (mp.pi / 10) ** 2, (mp.pi / 10) ** 3]])
b = mp.matrix([[mp.sin(-mp.pi / 20)],
[mp.zero],
[mp.sin(mp.pi / 20)],
[mp.sin(mp.pi / 20)]])
# Uncommenting the line below allows the QR householder method to work. I
# think this is because it allows a slight non-zero result to appear so
# execution can continue. Any 'dps' value greater than 10 causes
# failure.
#
# mp.dps = 10
x_qr, residual_qr = mp.qr_solve(A, b) # <- Gives error.
print(f"\nSolution using QR decomposition x_qr:\n{x_qr}")
print(f"Residual Ax_qr - b:\n{residual_qr}")
Error message:
Traceback (most recent call last):
File "...\mpmath_issue.py", line 48, in <module>
x_qr, residual_qr = mp.qr_solve(A, b)
^^^^^^^^^^^^^^^^^
File "...\mpmath\matrices\linalg.py", line 407, in qr_solve
H, p, x, r = ctx.householder(ctx.extend(A, b))
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "...\mpmath\matrices\linalg.py", line 341, in householder
raise ValueError('matrix is numerically singular')
ValueError: matrix is numerically singular
By comparison, lu_solve(A, b) has no problem:
print(f"Condition number of A: {mp.cond(A)}")
# Solution using LU decomposition -> Works fine.
x_lu = mp.lu_solve(A, b)
residual_lu = A * x_lu - b
print(f"\nSolution using LU decomposition x_lu:\n{x_lu}")
print(f"\nResidual Ax_lu - b:\n{residual_lu}")
Gives result:
Condition number of A: 694.870840206284
Solution using LU decomposition x_lu:
[ 0.0]
[1.16187485778415]
[ 0.0]
[-6.7270020477122]
Residual Ax_lu - b:
[0.0]
[0.0]
[0.0]
[0.0]
Cross-checking with Wolfram|alpha confirms the system is invertible / solvable, and lines up with the lu_solve() solution:
I believe the problem is likely this line in householder(ctx, A) in mpmath/matrices/linalg.py:
if not abs(s) > ctx.eps: # <- Problem area?
raise ValueError('matrix is numerically singular')
I suspect this check might be in the wrong location... it may not be required here at all? From memory this should check that if the sum is zero, then nothing further needs to be done with that particular column?