7
$\begingroup$

The motivation for my question is the following: in one of Project Euler questions there is a need to count the spanning trees of a rectangular grid graph of dimension 100x500. By the Matrix-Tree theorem this number is equal to the determinant of the matrix obtained by taking the graph Laplacian and deleting its first row and column. The Laplacian is very sparse and positive definite, but it is a 50000x50000 matrix and it seemed to me not feasible to calculate the determinant with the required precision (5 significant digits in exponential notation: may be I am wrong here?).

I have found the following interesting reference

http://arXiv.org/abs/0712.0681v3

which permits to compute the determinant as the determinant of a 100x100 matrix obtained as the left upper corner of the product of 500 200x200 sparse matrices, using the fact that the Laplacian is block tridiagonal (and deriving from the Laplacian a nonsingular 50000x50000 matrix whose determinant is equal to the determinant of the Laplacian with the first row and column deleted), which seemed to be a good thing. However the resulting real 100x100 matrix is nonsymmetric, dense, has some very large entries and some very small (e^(- ...)) eigenvalues (I experimented with smaller grids). Because of the very small eigenvalues I am not sure that any of the numerical methods (Gaussian elimination, SVD etc) will work properly, if I don't keep all the digits of the entries of the matrix. But this means to work with arbitrary precision which is probably very slow.

My question is: does anyone know of a practical method which can be used in such situation to compute the determinant with the required precision (again, 5 significant digits in exponential notation)?

I know the answer to the original Project Euler problem in the form of product of trigonometric terms, so really my interest is not in the solution itself but in these linear algebra questions which come out naturally when I try to apply the Matrix-Tree theorem directly.

Improve this question
$\endgroup$
3
  • 3
    $\begingroup$ There isn't much research in the calculation of determinants because they aren't in general a very useful number numerically. There are some serious scaling issues which prevent it from being used in a useful way. Theoretically of course you can prove some nice things using determinants. The following paper might be interesting to you, though: aimsciences.org/journals/pdfs.jsp?paperID=1224&mode=full $\endgroup$
    – Reid.Atcheson
    Mar 28, 2013 at 18:51
  • $\begingroup$ I tried this on Laplacians of smaller grids (with first row and column deleted), and it actually works, but the error bounds are only asymptotic and based on experiments with smaller grids I don't think one may obtain the required precision in a reasonable amount of time. I did not try it on the 100x100 matrix multiplied by its conjugate, one of the reasons being that the authors say themselves that numerical results on "near singular" matrices are not quite satisfactory. $\endgroup$
    – John Donn
    Apr 4, 2013 at 7:13
  • $\begingroup$ I think it likely that most methods will find issues on nearly singular matrices. Perhaps the best one can hope for really is to get within an order of magnitude. $\endgroup$
    – Reid.Atcheson
    Apr 6, 2013 at 4:18

4 Answers 4

Reset to default
4
$\begingroup$

Wikipedia actually has a nice overview of approaches: http://en.wikipedia.org/wiki/Determinant#Calculation As a general remark, people do not compute determinants of large matrices (large would here be of size >10 or >50) because this is numerically difficult and, likely, not very stable anyway. If you need to do it, I would see if maybe there are randomized algorithms in much the same way as you can use a Monte Carlo approach to approximate the trace.

Improve this answer
$\endgroup$
4
  • 1
    $\begingroup$ If the eigenvalues are reasonably well-behaved, say, between one and two, computing determinants is not that ridiculous. The computational complexity is no worse than computing an LU decomposition of the matrix (and then setting the determinant equal to the product of the diagonal values of the factors). Determinants of dense matrices as large as 5000 x 5000 are reasonable to compute on a laptop. $\endgroup$
    – Jack Poulson
    Mar 30, 2013 at 1:15
  • $\begingroup$ Right, but the matrix eigenvalues for the problem in the original post aren't going to be nicely scaled since they results from the graph Laplacian -- which, if the graph is not too dense, will have the same kind of behavior as the Laplacian on a finite element or finite difference mesh: namely a rather large condition number. $\endgroup$
    – Wolfgang Bangerth
    Mar 30, 2013 at 1:38
  • 2
    $\begingroup$ I managed to compute the 100x100 determinant using the simplest recursion (8) from "Sylvester's Identity and Multistep Integer-Preserving Gaussian Elimination" by Erwin H. Bareiss, linked from the wikipedia page. Using Python longs and implementing the recursion as it is with no particular optimizations, it took 7 hr to compute the 25094 digit integer on my 6 year old laptop. $\endgroup$
    – John Donn
    Apr 4, 2013 at 6:58
  • $\begingroup$ Nice :-) Very interesting to hear that it can actually be done! $\endgroup$
    – Wolfgang Bangerth
    Apr 4, 2013 at 12:42
5
$\begingroup$

The computation of integer-valued matrix determinants has been a subject of considerable research. Using exact arithmetic the Smith normal form can be computed, and from this diagonal form the determinant is easily found.

Saunders and Wan (2004), Smith Normal Form of Dense Integer Matrices, Fast Algorithms into Practice, say "Over the past thirty years, steady progress has been made in algorithms for the Smith forms." They give a reference for an asymptotically fastest known method (but which "may not be pracical"), and then dive into details of improvements on an algorithm by W. Eberly, M. Giesbrecht, and G. Villard.

Essentially one wants to find the greatest common divisor of matrix entries, forming this by a series of row and column operations and using it to eliminate all other nonzero entries in its column.

This unanswered Math.SE Question on "fastest method" for computing integer matrix determinants (in roughly the size you ask about) drew a couple of good Comments, one of which suggests (and reports a successful 100x100 experiment with) Maple's Determinant function using method = integer. Although their algorithm uses (exact) modular arithmetic, it entails a probabilistic element (a small but positive chance of giving a wrong answer). The reason for this is that the computation is carried out for a sequence of primes, reducing the precision needed for each particular run, and in the absence of reaching an apriori bound on the determinant (making the answer deterministic), "convergence" may be decided by obtaining the same result after several such runs (the probabilistic case).

Improve this answer
$\endgroup$
3
$\begingroup$

I have been using pari/gp (confuration displayed in output below) to calculate determinants of large, dense matrices. For example, I use a quad core i7, 16GB RAM, 32GB swap space on internal ssd, linux laptop, to calculate a matrix 1700x1700 to a precision of a little less than 3000 digits. That takes about 17 hours, using all eight threads, 100% of the time, but only a little RAM. Then, it does the determinant in another six hours, one core uses 100% for the determinant, and kswapd uses another core, 100% of the time (and appears to be able to keep up, and use all the swap space). The determinant is not parallelized in pari/gp (as far as I know, except in that kswapd and the gp instance use separate cores). Here is a specific example (that doesn't use parallel code), calculates a random (floating point), 1000-digit precision, 300x300 matrix determinant:

I know it's not matlab, but in my mind, this is a better way.

=========================== cut here
default(nbthreads,2); \\this example sets up two but uses only one thread.
default(threadsize,  500 000 000);
default(parisize,  1 800 000 000); \\ make sure there is enough RAM allocated

initrand()=
{
  cmd = "date +\"%s\"";           \\ a linux command to get the time
  seed = eval(externstr(cmd)[1]); \\ evaluate the external command
  setrand(seed);                  \\ set the random seed in pari/gp
  print ("Random Number Seed = ", seed); \\ this is what it is
  return(seed);                   \\ don't really use this value any more
}
{
  default(realprecision, 1000);  \\ accurate to 1000 digits of precision
  initrand();                    \\ pari/gp seems to need a different seed every time
  M = matrix(300,300,i,j,random(1.0)); \\ some small random matrix
  printf ("%10.7f\n", M[1,1]);   \\ print the first matrix element
  detM = matdet(M);              \\ calculate the determinant
  print (detM);                  \\ print the result of the determinant
  quit(); 
}
========================================= cut here

This is the output, which took a minute or so on this eight-year old laptop (dual core, 2GB ram, not a fast ThinkPad).

                        GP/PARI CALCULATOR Version 2.7.3 (released)
                amd64 running linux (x86-64/GMP-6.0.0 kernel) 64-bit version
             compiled: Mar 20 2015, gcc version 4.8.2 (Ubuntu 4.8.2-19ubuntu1) 
                                 threading engine: pthread
                       (readline v6.3 enabled, extended help enabled)

                           Copyright (C) 2000-2015 The PARI Group

PARI/GP is free software, covered by the GNU General Public License, and comes WITHOUT ANY WARRANTY WHATSOEVER.

Type ? for help, \q to quit.
Type ?12 for how to get moral (and possibly technical) support.

 parisize = 8000000, primelimit = 500000
   ***   Warning: new stack size = 1800000000 (1716.614 Mbytes).
 Random Number Seed = 1442728837
 0.1458753
 -34107869659292788498041296259204073212856192996487539539650250444333175626443030051167504275069295850290060165906721403518710248631573270877318755.16606345550957294257524953917496320048631959234301670865079024142295995123023952749154306602088800645390346232704746782547685913565513791084278569243333238872632785910522536912572116356801097060790810716175151647494785903522152367435398630683566431045727806979724702836066399155098909578586486867501579225868835988817237083450602429253202805656099228929576538193846647490928888795367874033893827126923874835063405546089363260880884557106553366031910609740691344849673406340416580053579023527252144631532815837507760281113528779522639792441184911588021824240226126036961995649658869925314847019638680109421905217088024624421267647700950443089481208575217193882687232501653341331383260930401032570115570957029275695235564734394098177016567882258012726222978777317463869837697774536180341558237759707887932943904115121644125119662351721076540409972736105729
Goodbye!

Of course, in this case, the "-341......05729" is the floating point result. I haven't played with integer matrix elements, but I'm sure it is just as good, if not better.

Improve this answer
$\endgroup$
1
  • $\begingroup$ To do integers, you don't need "default(realprecision,1000);", that is for floats; and change "random(1.0)" to "random(2001)-1000" to make a random matrix with integers in the interval [-1000,1000]; and modify the first print statement to print an integer "print(M[1,1]);" is ok. A 300x300 result takes a few seconds on my old laptop. 1000x1000 takes perhaps a few minutes on a faster computer. $\endgroup$
    – user122986
    Sep 20, 2015 at 6:47
1
$\begingroup$

The Integer Matrix Library provides a C function mDeterminant (const FiniteField p, Double *A, const long n) to compute the determinant det(A) mod p of an integer matrix A (using a row echelon transform), which might be of use.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.