4

Consider a quadratic $n\times n$ Matrix $A$ and the general question "how to find the determinant $\det(A)$ when too lazy for a Laplace Expansion but lucky enough to get a singular value decomposition for free".

Why? Within the C++ library Lapacke, which I use, for all its power there is no determinant function. However, there is a multitude of state-of-the-art singular value decomposition tools xgesvd(..). Hence we start out at

$$A=USV'$$

where $U,V$ are orthonormal and $S$ is a diagonal matrix which's diagonal is comprised of the singular values in descending order. Already the aim comes tantalizingly close:

$$\det(A)=\det(USV')=\det(U)\det(S)\det(V')=\pm\prod\limits_{j\in n}s_{j,j}$$

the "$\pm$" coming from the fact that it is unknown whether $\det(U)$ or $\det(V)$ are $+1$ or $-1$.

Hence the question: Is there a lazier way than a complete expansion to simply determine the sign of a determinant of ideally a general but at least an orthonormal matrix? Or of a broader frame of mind: Is there another simple Lapack-compatible way than my SVD approach?

Markus-Hermann
  • 143
  • Relevant thread: http://math.stackexchange.com/questions/205427/determinant-of-a-big-matrix – littleO Feb 28 '14 at 12:23
  • @little(): Thanks for your link. One of the answers uses a LU decomposition. My general lazyness would allow me to calculate the det(T) where T is triangular (being just the product of the diagonal elements). Why don't you put it in a proper answer? I would accept it. – Markus-Hermann Feb 28 '14 at 12:36

2 Answers2

3

Use the QR decomposition. LAPACK produces it using Householder reflectors, which each have determinant $-1$, so you only need the number of reflectors, typically $n-1$ for an $n\times n$ matrix, and the non-unit part of the determinant is the product over the diagonal entries of $R$.

Since the QR decomposition is one of the first steps in the SVD computation, this will also be faster. And usually, QR decomposition is numerically more stable than the LU decomposition, even with pivoting.

Lutz Lehmann
  • 126,666
  • Thank you, LutzL. This sounds good. But I will have to think about it first, before I accept your answer (which will happen sometime next week I am afraid). Two small issues: 1.) If memory serves there was (like in LU with pivoting) a permutation matrix involved (maybe it remains a unity matrix if all diagonal elements of the triangular matrix are !=0 which is the only non-trivial case for det(A)), else I guess I would have to calculate sign(P) first, as of yet no idea how, but it sounds solvable). [..] – Markus-Hermann Feb 28 '14 at 16:53
  • [..] 2.) I wasn't aware yet, that a Householder transform has determinant -1 :-) So your answer really stands and falls with the question "does Lapack truly do (n-1) HH transforms or is there at least a way to determine how many they do". This last one will require some minimal research. Independent of (1): If the answer to (2) is "yes" I will accept your answer. I would upvote you already, but I am afraid my newby-rep isn't sufficient so far. – Markus-Hermann Feb 28 '14 at 16:53
  • The simple Householder reflector $S_v=I-2\frac{vv^T}{v^Tv}$ flips two half spaces and leaves the hyperplane of reflection orthogonal to $v$ invariant. So one dimension gets multiplied by $-1$, all other by $1$ in an adapted orthogonal basis, giving determinant $-1$. – Lutz Lehmann Feb 28 '14 at 16:58
  • Learned something new here! And by now I also found http://www.netlib.no/netlib/lapack/double/dgeqrf.f with a doc that fulfilled the criterion for accepting your answer. Thanks again! – Markus-Hermann Feb 28 '14 at 17:03
  • No pivots or permutations, but the LAPACK protocol for QR needs a double tau array for the coefficients in $S_k=I-\tau_kv_kv_k^T$. The original matrix is overwritten with the $R$ matrix, including the diagonal, and the vectors $v_k$ in the index range $(k+1:n)$ that fit below the diagonal under the assumption that the scaling is such that $v_k[k]=1$ (and $v_k[1:k-1]=0$). – Lutz Lehmann Feb 28 '14 at 17:03
  • 1
    For complex matrices it is a bit more complicated. The determinant of a LAPACK Householder reflection is not simply -1, but $\det(I-\tau·v·v^H) = 1-\tau·v^H·v$. CGEQRF chooses $\tau$ and $v$ such that the determinant has absolute value 1. – Lemming Apr 06 '18 at 08:04
  • @Lemming : But then again the question for the sign, in the sense of positive or negative, of the determinant of a complex matrix is usually less meaningful. – Lutz Lehmann Apr 06 '18 at 08:33
0

For complex matrices it is a bit more complicated. The determinant of a LAPACK Householder reflection is not simply -1, but $\det(I-\tau·v·v^H) = 1-\tau·v^H·v$. (This is Silvester's determinant identity applied to single-row and single-column matrices.) CGEQRF chooses $\tau$ and $v$ such that the determinant has absolute value 1. More precisely, $1-\tau·v^H·v = -\frac{\tau^2}{|\tau|}$

Lemming
  • 101