Add per-iteration condition number estimate

[nickmccleery]

Since we wound up rolling our own QR based solve method on the back of newton_faer, that actually buys us the opportunity to spit out some nice diagnostic stats at each inner (linear) system solve.

In our python version we only have condition number for the R matrix, and not A, but I think we should be able to look for and flag numerical issues all the same. For reference, there we have:

    def factor(self, matrix):
        """
        Factor the matrix using QR decomposition.
        """
        if isinstance(matrix, sparse.csr_matrix):
            matrix = matrix.toarray()
        try:
            result = qr(matrix, mode="economic", pivoting=True)
            if len(result) == 3:
                # Cast to the specific return type we expect.
                q, r, p = cast(Tuple[np.ndarray, np.ndarray, np.ndarray], result)
                self.q = q
                self.r = r
                self.p = p
                self.matrix_original = matrix

                # Build stats.
                diag = np.diag(r)  # Diagonal of R.
                diag_abs = np.abs(diag)  # Absolute values of diagonal.
                s_max = diag_abs.max() if diag_abs.size else 0.0
                s_min = diag_abs.min() if diag_abs.size else 0.0
                condition_number = (s_max / s_min) if (s_min > 0) else np.inf
                rank = np.linalg.matrix_rank(r)

                # Check reconstruction error: A[:,p] = QR
                matrix_permutation = matrix[:, p]
                recon_err = np.linalg.norm(matrix_permutation - np.matmul(q, r)) / (
                    np.linalg.norm(matrix_permutation) + 1e-16
                )

                logger.debug(f"System shape: {matrix.shape}")
                logger.debug(f"Rank estimate (np.linalg.matrix_rank(R)): {rank}")
                logger.debug(f"Frobenius norm, Q: {np.linalg.norm(q)}")
                logger.debug(f"Frobenius norm, R: {np.linalg.norm(r)}")
                logger.debug(f"Reconstruction error: {recon_err}")
                logger.debug(f"Condition number: {condition_number}")
            else:
                raise SolverError("QR factorization did not return expected results")
        except Exception as e:
            raise SolverError(f"Dense QR factorization failed: {e}")

We might benefit from implementing something similar here.

Related: #6