1
$\begingroup$

I have encountered a problem with Mathematica that I really don't know how to solve. I'm sorry if I am going to be very verbose but I need it to explain the problem properly (and also make the code I'm going to show more comprehensible).

I'm using Mathematica to work for my master's thesis and in these days I'm dealing with the ecosystem model described in this article. Amongst the things I have to do there's the study of the eigenvalue distribution of the community matrix of this ecosystem, i.e. the jacobian matrix of the system of coupled ODEs that describe the population abundances of the various species, computed at the steady state of the system.

In a single sentence my problem is the following: as the size of this matrix increases Mathematica finds that it has also very small positive eigenvalues, while this matrix is always definite negative.

More in detail (now I have to be very verbose): the model describes a system of $m$ animal species competing for $p$ resources. These resources are supplied with constant rates $\vec{s}=(s_1,\dots,s_p)$, and every species can be identified with its metabolic strategy $\vec{\alpha}_{\sigma}=(\alpha_{\sigma 1}, \dots, \alpha_{\sigma p})$ (where $\sigma=\{1,\dots,m\}$ labels the species, and $\alpha_{\sigma i}$ can be interpreted as the ratio at which species $\sigma$ eats resource $i$). Now, assuming $m \geq p$ and that $\vec{\alpha}_\sigma$-s (properly renormalized) satisfy $$\sum_{i=1}^p \alpha_{\sigma i} = 1 \quad \forall \sigma$$ (called "metabolic trade-off condition") it can be shown that the system is stable (i.e. all species survive) if $\vec{s}$ belongs to the convex hull of metabolic strategies, i.e. if: $$\exists \{n_1^*>0, \dots, n_m^*>0\} : \sum_{\sigma=1}^m n_\sigma^* = 1 \quad \text{and} \quad \vec{s} = \sum_{\sigma=1}^m n_\sigma^* \vec{\alpha}_\sigma$$ Now, studying the model analytically it turns out that the jacobian matrix of the system computed at this stable steady state can be written as: $$\mathcal{M}=-DASA^T$$ where: $$D=\operatorname{diag}(n_1^*, \dots, n_m^*) \qquad S=\operatorname{diag}(1/s_1,\dots,1/s_p) \qquad A_{\sigma i} = \alpha_{\sigma i}$$ and it can be shown that for $m \geq p$ it is always definite negative.

Now, with Mathematica I'm trying to study numerically the eigenvalue distribution of $\mathcal{M}$. The code I am using in order to do so is the following:

m = 1000;
p = m;
A = #/Total[#, {2}] &@RandomReal[{0, 1}, {m, p}];       
\[ScriptCapitalN] = (#/Total[#, {2}] &@RandomReal[{0, 1}, {1, m}])[[1]];
S = Sum[\[ScriptCapitalN][[i]]*A[[i]], {i, 1, m}];
\[ScriptCapitalD] = DiagonalMatrix[\[ScriptCapitalN]];
\[ScriptCapitalS] = DiagonalMatrix[1/S];
\[ScriptCapitalM] = -\[ScriptCapitalD].A.\[ScriptCapitalS].Transpose[A];
eigenvalues = Eigenvalues[\[ScriptCapitalM]];
ListPlot[{Re[#], Im[#]} &/@ eigenvalues, PlotRange -> Automatic]
Length[Select[eigenvalues, # > 0 &]]
Max[eigenvalues]

(I added the last two lines to keep track how many eigenvalues are positive and which is the largest eigenvalue).

Now, If I set $m=p$ no problems arise and all eigenvalues are indeed negative. For example, for $m=p=1000$ I get the plot in figure $a)$ (below).

However, problems start popping up if I set $m > p$ (which is the most interesting case for the model). In fact, the results of the analytical study of $\mathcal{M}$ (i.e. the fact that it's negative definite) are valid for $m \geq p$, but using the code I have written it happens that some eigenvalues start "collapsing" around 0, and they turn out to be both positive and negative (even if very very small in magnitude). The number of "collapsed" eigenvalues increases as the ratio $m/p$ increases, until they become the majority of the eigenvalues.

For example, if I use the same code with $m=1000$ and $p=100$ I get figure $b)$ below.

My question is therefore the following: is there something am I doing wrong in the code? Or is this just a problem of numerical precision?

I first thought it could be the latter, so I have tried to increase the numerical precision a bit by changing the third and fourth line of the previous code to:

 A = #/Total[#, {2}] &@
   RandomReal[{0, 1}, {m, p}, WorkingPrecision -> 7];
\[ScriptCapitalN] = (#/Total[#, {2}] &@
     RandomReal[{0, 1}, {1, m}, WorkingPrecision -> 7])[[1]];

However, new problems arise: apart from the fact that computation now requires a lot more time (I had to use much smaller values of $m$ and $p$ in order to get results in minutes) the eigenvalues are now complex, instead of purely real! For example, for $m=500$; $p=80$ I get figure $c)$ below.

Therefore, if this is indeed a problem of numerical precision, how can I fix it? Is there a more efficient way to increase it?

Here are the plots: enter image description here

Sorry again for the verboseness.

$\endgroup$
  • $\begingroup$ WorkingPrecision -> 7 - you're going in the wrong direction. Try cranking it up to, say, WorkingPrecision -> 20, or WorkingPrecision -> $MachinePrecision if you do not wish to be bold. $\endgroup$ – J. M. is away Apr 14 '17 at 18:06
  • $\begingroup$ I have tried with higher numbers like 20 or with $MachinePrecision, but the computation time increases immensely for large matrices... It really becomes infeasible this way. Are there more efficient ways to get more precision without having to wait eons for a result, or do I just have to get over it? $\endgroup$ – LPM Apr 14 '17 at 21:55
2
$\begingroup$

Your first matrix is negative definite. Your second one has rank 100 and dim 1000x1000, so it is negative semidefinite. The smallest in magnitude eigenvalue for the first is around 10^(-11), so it is sufficiently distinct from machine arithmetic fuzz that it is not messed up by numeric error. The 900 zero eigenvalues in the second one are not going to be found as exact zeros by any numerical method, and that's why there is a fairly symmetric range of machine fuzz values appearing around the origin.

$\endgroup$
  • $\begingroup$ So, if I have understood correctly, in the case m>=p I am getting those numbers around zero because (since the matrix is now negative semidefinite) there are some null eigenvalues which numerically are not zero but those very small numbers symmetrically distributed around zero? $\endgroup$ – LPM Apr 14 '17 at 21:52
  • $\begingroup$ Anyway thanks, I had completely missed the fact that in reality the matrix is negative SEMIdefinite. Maybe I have to review some basic algebra hehe $\endgroup$ – LPM Apr 14 '17 at 22:08
  • $\begingroup$ Correct, they are numerically zero and, approximately, show up as symmetrically distributed around the origin. $\endgroup$ – Daniel Lichtblau Apr 15 '17 at 0:21
  • $\begingroup$ Perfect, thank you very much! $\endgroup$ – LPM Apr 15 '17 at 7:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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