# GSL linear algebra LU/determinant precision

I am working with symmetric matrices of order $n \times n$ where $n \leq 50$. The diagonal elements of my matrices are a fixed number $d$ and the off diagonal elements are limited to two small numbers $\{a,b\}$.I know $a,b,d$ in advance.

For any such matrix I need to be able to determine whether its determinant is negative or not. I must make sure that I do not mistakenly mark a matrix with non-negative determinant as having a negative determinant. I am not concerned with the converse mistake.

My question is

What would be a sensible way to determine whether the computed determinant is indeed negative?

The way to go would be to determine the precision guarantee of the algorithms implemented by GNU GSL but unfortunately the documentation does not mention that for these specific functions (gsl_linalg_LU_decomp and gsl_linalg_LU_det)

I've already asked this on the relevant GSL mailist but just to get a different perspective, I was wondering, what would be a good way to figure a test that would be 100% precise?

Right now I am just testing whether the determinant is smaller than a rough negative number but I would like to have the smallest possible value $\beta$ such that the determinant is negative whenever it is smaller than $\beta.$

1. I looked up forward LU error bounds in Higham's Accuracy and Stability of Numerical Algorithms, Theorem 9.15 (citing Barrlund and Sun), for the LU decomposition $$A=LU,\quad A+\Delta A=(L+\Delta L)(U+\Delta U)$$ it gives the norm-wise error bound $$\begin{gather} \frac{\|\Delta U\|_F}{\|U\|_2} \leq \frac{\|L^{-1}\|_2\|U^{-1}\|_2\|A\|_2}{1-\|L^{-1}\|_2\|U^{-1}\|_2\|\Delta A\|_2}\frac{\|\Delta A\|_F}{\|A\|_F} =: M, \\ \|\Delta A\|_\infty \leq n^2\gamma_{3n}\rho_n \|A\|_\infty, \end{gather}$$ where $\gamma_k = k\epsilon/(1-k\epsilon)$ ($\epsilon$ is the unit roundoff), and $\rho_n$ is the growth factor, defined as $$\max_{i,j,k} |a^{(k)}_{i,j}|/\max_{i,j} |a_{i,j}|$$ with $a^{(k)}_{i,j}$ being the matrix elements at $k$-th stage of LU factorization by Gaussian elimination.

So, in principle, so long as you can estimate all of the above numbers, the procedure should produce the correct determinant sign so long as the smallest value on the diagonal of $U$ is not so close to zero that its sign could have been influenced by numerical errors: $$\min_k |u_{kk}| \geq M\|U\|_2 \geq \|\Delta U\|_F.$$ Of course, this is not very elegant. Higham also gives a component-wise error bound, which should be stricter.

2. Not likely, but if you know something about how the numbers $a,b$ appear in off-diagonal elements, it is possible there might be a lower bound $m = \min_A |\det A|$. If it so happens that it is easy to calculate, then the determinant's sign can be accurately determined so long as the numerical errors are at most $m$.

3. Eigendecomposition, although more expensive, might also be helpful. For example, GSL (https://www.gnu.org/software/gsl/manual/html_node/Real-Symmetric-Matrices.html#Real-Symmetric-Matrices) promises that

The computed eigenvalues are accurate to an absolute accuracy of \epsilon ||A||_2, where \epsilon is the machine precision.

This would directly address the issue of whether the determinant's sign is correct - that would require all eigenvalues to satisfy $|\lambda| > \epsilon \|A\|_2$.

4. In general, computation of a sign of some quantity is ill-conditioned when that quantity is small - the determinant's sign would be very sensitive to small changes in the input matrix. So the common solution to this problem is just to look for a way to avoid computing the determinant's sign at all, but I don't know how feasible that is.

5. Arbitrary-precision floating-point arithmetic might also help. Although not with GSL, there are libraries (e.g., Eigen) that implement linear algebra in a way that can work with, for example, mpfr.

6. (Edit.) The $LDL^\top$ decomposition is even easier I think: from Barrlund (Eq. 2.1b), ($D$ is the diagonal from $U$, so I don't think it matters if you actually compute $LU$ instead because this is a perturbation analysis) $$\|\Delta D\|_F \leq \frac{\kappa_2(A)}{1-\kappa_2(A)\frac{\|\Delta A\|_2}{\|A\|_2}} \|\Delta A\|_F = \beta,$$ which now only needs estimates of $A$'s 2-norm condition number $\kappa_2(A)$ and the normwise error backward error estimates $\|\Delta A\|_2$, $\|\Delta A\|_F$ from above. The correct sign is now guaranteed by something like $|d_{kk}|\geq \beta$.

• Thank you for your expansive answer. In fact I needed to check that the eigenvalues are non-negative and I wanted to speed up the test by using the fact that my matrices can only have one negative eigenvalue. Considering all the points you mention it seems that the safest and most efficient way is indeed to just compute the eigenvalues. – Jernej Oct 9 '15 at 8:55
• @Jernej Perhaps Cholesky decomposition would be easier since it's really about being psd? If it fails, it's not psd, if it succeeds with $\tilde A = \tilde R^\intercal\tilde R$, then it could only possibly have failed if one of $R$'s diagonal elements is small, in which case you can require that it be larger than (Higham - Theorem 10.8) $$2^{-1/2}\frac{\kappa_2(A)\epsilon}{1-\kappa_2(A)\epsilon}\|R\|_p \geq \|\Delta R\|_F, \qquad \epsilon = \frac{\|\Delta A\|_F}{\|A\|_p}, \qquad p\in\{2,F\},$$ (there are also componentwise bounds) which is similar to the LU case, but rather easier to apply. – Kirill Oct 9 '15 at 22:30
• @Jernej Hm, my last comment answers the converse of your question ("make sure a matrix is p.s.d.", not what you asked). – Kirill Oct 9 '15 at 22:47
• @Jernej See edit. – Kirill Oct 9 '15 at 22:57