R devel - Nov 2003 - bug in det using method="qr" (PR#1244) (PR#4450)

I just detected, that det() is not working on complex matrices any more,
due to the fix to the bug reports noted above. I am not happy with this,
as determinants are perfectly usable on complex matrices.

AFAIUI the bugs resulted from less than optimal behaviour of qr() in
certain cases. IMHO this is due to the unhappy decision to use a default for
parameter tol to decide whether the the decomposition is rank deficient.

A better fix for (PR#1244) should be considered. I propose, using qr method
with tol=0 as default, for det(). Even more preferrable, tol=0 can be made the
default for qr(), forcing all applications to set a reasonable tol for their
own.

wbk
-- 
Dipl.-Math. Wilhelm Bernhard Kloke
Institut fuer Arbeitsphysiologie an der Universitaet Dortmund
Ardeystrasse 67, D-44139 Dortmund, Tel. 0231-1084-257

Which version of R are you running?  If it is 1.8.0 then the bug
report doesn't make sense because the NEWS file for 1.8.0 has an
entry
     o	det() uses an LU decomposition and LAPACK.  The `method'
	argument to det() no longer has any effect.
so method="qr" is redundant.  This is documented in ?det as well:

Note:

     Often, computing the determinant is _not_ what you should be doing
     to solve a given problem.

     Prior to version 1.8.0 the 'det' function had a 'method'
argument
     to allow use of either a QR decomposition or an
     eigenvalue-eigenvector decomposition.  The 'determinant' function
     now uses an LU decomposition and the 'det' function is simply a
     wrapper around a call to 'determinant'.

The default method for det() now calls determinant() and it is
determinant() that does not have a method for complex matrices.

If you can propose a stable calculation for the determinant of a
complex matrix (including what should be done for logarithm = TRUE) we
can add it to the list of things to implement.  If you provide an
implementation then it will be available even faster.

The reason that we switched from det() to determinant() with the
optional logarithm parameter is because calculation of a determinant
does not scale well and often gives nonsense answers.  (Also, as
stated above, often computing a determinant is not what you should be
doing to solve a given problem.)  Philippe Grosjean used a det
calculation of a large (600 by 600 IIRC) matrix with random standard
normal entries as part of his benchmark suite.  It is a good idea (and
we appreciate Philippe doing the benchmarking) except if you check
what the result is; it is almost always +/- Inf so the answer is
meaningless.


wb@arb-phys.uni-dortmund.de writes:
> I just detected, that det() is not working on complex matrices any more,
> due to the fix to the bug reports noted above. I am not happy with this,
> as determinants are perfectly usable on complex matrices.
> 
> AFAIUI the bugs resulted from less than optimal behaviour of qr() in
> certain cases. IMHO this is due to the unhappy decision to use a default
for
> parameter tol to decide whether the the decomposition is rank deficient.
> 
> A better fix for (PR#1244) should be considered. I propose, using qr
> method with tol=0 as default, for det(). Even more preferrable,
> tol=0 can be made the default for qr(), forcing all applications to
> set a reasonable tol for their own.

>>>>> "Wilhelm" == Wilhelm B Kloke
<wb@arb-phys.uni-dortmund.de>
>>>>>     on Wed, 12 Nov 2003 13:37:13 +0100 (CET) writes:
    Wilhelm> I just detected, that det() is not working on
    Wilhelm> complex matrices any more, due to the fix to the
    Wilhelm> bug reports noted above. I am not happy with this,
    Wilhelm> as determinants are perfectly usable on complex
    Wilhelm> matrices.

Btw, this is not a bug in the strict sense, 
but a feature request;  OTOH, I know that det() did work for
complex matrices.

    Wilhelm> AFAIUI the bugs resulted from less than optimal
    Wilhelm> behaviour of qr() in certain cases. IMHO this is
    Wilhelm> due to the unhappy decision to use a default for
    Wilhelm> parameter tol to decide whether the the
    Wilhelm> decomposition is rank deficient.

    Wilhelm> A better fix for (PR#1244) should be considered. I
    Wilhelm> propose, using qr method with tol=0 as default, for
    Wilhelm> det(). Even more preferrable, tol=0 can be made the
    Wilhelm> default for qr(), forcing all applications to set a
    Wilhelm> reasonable tol for their own.

I'm not entirely sure about the correctness of your
interpretation of bug fixes here.
Nevertheless, I agree that we should consider allowing
determinant() and hence det() to work for complex matrices --
and change determinant.matrix()'s default argument
`logarithm = TRUE' to 
`logarithm = !is.complex(x)'

Then replace

    if (is.complex(x)) 
        stop("determinant not currently defined for complex matrices")

by something like (we had in R 1.6)

    if (is.complex(x)) {
        qx <- qr(x, tol = 0)
        if (qx$rank < n) 
            return(0)
        x <- prod(diag(qx$qr))
        return(if (n%%2 == 1) x else -x)
    }

---

Martin Maechler <maechler@stat.math.ethz.ch>
http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum  LEO C16	Leonhardstr. 27
ETH (Federal Inst. Technology)	8092 Zurich	SWITZERLAND
phone: x-41-1-632-3408		fax: ...-1228			<><