4

I am new to Matlab. I would like to check the so call "logarithmic law" for determinant of random matrices with Matlab, but still do not know how.

Logarithmic law:

Let A be a random Bernoulli matrix (entries are iid, taking value +-1 with prob. 1/2) of size n by n. We may want to compare the probability density function of (log(det(A^2))-log(factorial(n-1)))/sqrt(2n) with the pdf of Gaussian distribution. The logarithmic law says that the pdf of the first will approach to that of the second when n approaches infinity.

My Matlab task is very simple: check the comparison for, say n=100. Anyone knows how to do so?

Thanks.

Amro
  • 123,847
  • 25
  • 243
  • 454
H. H.
  • 41
  • 1
  • [pdf](http://www.mathworks.com.au/help/toolbox/stats/pdf.html) is here in Matlab docs. –  Aug 29 '11 at 02:02
  • @H.H.: I did a quick a search on google for "Logarithmic Law for Random Determinants" and found [two](http://books.google.com/books?id=fsMamZXZwIsC&pg=PA60) [book references](http://books.google.com/books?id=LeGMfJwUK3cC&pg=PA47).. I think the last term should be: `sqrt(2*log(n))` – Amro Aug 29 '11 at 03:06

1 Answers1

5

Consider the following experiment:

n = 100;                           %# matrix size
num = 1000;                        %# number of matrices to generate

detA2ln = zeros(num,1);
for i=1:num
    A = randi([0 1],[n n])*2 - 1;  %# -1,+1
    detA2ln(i) = log(det(A^2));
end

%# `gammaln(n)` is more accurate than `log(factorial(n-1))`
myPDF = ( detA2ln - gammaln(n) ) ./ sqrt(2*log(n));
normplot(myPDF)

enter image description here

Note that for large matrices, the determinant of A*A will be too large to represent in double numbers and will return Inf. However we only require the log of the determinant, and there exist other approachs to find this result that keeps the computation in log-scale.

In the comments, @yoda suggested using the eigenvalues detA2(i) = real(sum(log(eig(A^2))));, I also found a submission on FEX that have a similar implementation (using LU or Cholesky decomposition)

Amro
  • 123,847
  • 25
  • 243
  • 454
  • Thanks a lot, it works for me. (btw, there was a typo in my formula, the denominator should be sqrt(2*log(n)). – H. H. Aug 29 '11 at 11:39
  • Also, is there anyway to check for matrices of larger size? My Matlab works for n=100, but not for n=200, and how about n=1000? The comparison for n=100 seems still unconvincing. – H. H. Aug 29 '11 at 11:41
  • 1
    @H.H. Replace the last line in the `for` loop with this: `detA2(i) = real(sum(log(eig(A^2))));` This will be slower, but will let you run it for matrices as large as your system/MATLAB will allow. – abcd Aug 29 '11 at 16:41
  • @yoda: thanks I added a note how to safely compute the log of the determinant as you suggested – Amro Aug 29 '11 at 17:30