Determinants of huge matrices in MATLAB

from a simulation problem, I want to calculate complex square matrices on the order of 1000x1000 in MATLAB. Since the values refer to those of Bessel functions, the matrices are not at all sparse.

Since I am interested in the change of the determinant with respect to some parameter (the energy of a searched eigenfunction in my case), I overcome the problem at the moment by first searching a rescaling factor for the studied range and then calculate the determinants,

result(k) = det(pre_factor*Matrix{k});

Now this is a very awkward solution and only works for matrix dimensions of, say, maximum 500x500.

Does anybody know a nice solution to the problem? Interfacing to Mathematica might work in principle but I have my d开发者_运维技巧oubts concerning feasibility. Thank you in advance

Robert

Edit: I did not find a convient solution to the calculation problem since this would require changing to a higher precision. Instead, I used that

ln det M = trace ln M

which is, when I derive it with respect to k

A = trace(inv(M(k))*dM/dk)

So I at least had the change of the logarithm of the determinant with respect to k. From the physical background of the problem I could derive constraints on A which in the end gave me a workaround valid for my problem. Unfortunately I do not know if such a workaround could be generalized.

You should realize that when you multiply a matrix by a constant k, then you scale the determinant of the matrix by k^n, where n is the dimension of the matrix. So for n = 1000, and k = 2, you scale the determinant by

>> 2^1000
ans =
     1.07150860718627e+301

This is of course a huge number, so you might expect that it should fail, since in double precision, MATLAB will only represent floating point numbers as large as realmax.

>> realmax
ans =
     1.79769313486232e+308

There is no need to do all the work of recomputing that determinant, not that computing the determinant of a huge matrix like that is a terribly well-posed problem anyway.

If speed is not a concern, you may want to use det(e^A) = e^(tr A) and take as A some scaling constant times your matrix (so that A - I has spectral radius less than one).

EDIT: In MatLab, the log of a matrix (logm) is calculated via trigonalization. So it is better for you to compute the eigenvalues of your matrix and multiply them (or better, add their logarithm). You did not specify whether your matrix was symmetric or not: if it is, finding eigenvalues are easier than if it is not.

You said the current value of the determinant is about 10^-300.

Are you trying to get the determinant at a certain value, say 1? If so, rescaling is awkward: the matrix you are considering is ill-conditioned, and, considering the precision of the machine, you should consider the output determinant to be zero. It is impossible to get a reliable inverse in other words.

I would suggest to modify the columns or lines of the matrix rather than rescale it.

I used R to make a small test with a random matrix (random normal values), it seems the determinant should be clearly non-zero.

> n=100
> M=matrix(rnorm(n**2),n,n)
> det(M)
[1] -1.977380e+77
> kappa(M)
[1] 2318.188

This is not strictly a matlab solution, but you might want to consider using Mahout. It's specifically designed for large-scale linear algebra. (1000x1000 is no problem for the scales it's used to.)

You would call into java to pass data to/from Mahout.