I wish to find the numerical eigenvalues of a relatively small (100 dimensional) Hermitian matrix. While the eigenvalue calculation itself is fast, the construction of the matrix is taking longer than I would prefer.

Before I provide the code, I'll describe the matrix I am constructing. The matrix entries are given by

$$M_{i,j} = \begin{cases} \begin{cases} \epsilon_{\frac{i}{2},l+1} & i \in \text{Even}\\ -\epsilon_{\frac{i+1}{2},-l} & i \in \text{Odd} \end{cases} &,\,\,i=j\\ \Delta^{(l)}_{\frac{i}{2},\frac{j+1}{2}} &,\,\,i \in \text{Even}\,\&\,j \in \text{Odd}\\ \Delta^{(l)*}_{\frac{j}{2},\frac{i+1}{2}} &,\,\,i \in \text{Odd}\,\&\,j \in \text{Even} \end{cases} $$

Here $$\begin{align} \Delta_{n,n'}^{(l)} &= i 2\pi \Delta (-1)^l \left(\frac{\alpha_{n',l}}{R} \right)\mathcal{N}_{n,l+1}\mathcal{N}_{n',-l} \int_0^R dr\,r\,J_{l+1}\left(\alpha_{n,l+1}\frac{r}{R}\right)J_{l+1}\left(\alpha_{n',l} \frac{r}{R}\right) \nonumber \\ \Delta_{n',n}^{(l)*}&= -i 2\pi \Delta (-1)^l \left(\frac{\alpha_{n',l+1}}{R} \right)\mathcal{N}_{n,-l}\mathcal{N}_{n',l+1}\int_0^Rdr\,r\,J_l\left(\alpha_{n,l}\frac{r}{R}\right)J_l\left(\alpha_{n',l+1}\frac{r}{R} \right) \end{align}$$ with $$\mathcal{N}_{n,l} = \frac{1}{\sqrt{\pi}R |J_{l+1}(\alpha_{n,l})|}$$ and $$\epsilon_{n,l} = \frac{\alpha_{n,l}^2}{2 R^2} - \mu$$ where $\alpha_{n,l}$ is the $n^{th}$ zero of $J_l(x)$, the Bessel function of the first kind.

Let me also describe what I've already done to reduce the computation time.

Initially, I was performing the integrals numerically and the whole process of populating the matrix took a long time. I then realized that Mathematica could compute these analytically, so I used those simplified expressions instead. This brought the computation time down by half.

I was also populating the matrix using SparseArray initially, but then moved to ParallelTable (see attached code), which brought the computation time down by another factor of about 2.5. However, I am stuck at this point and cannot make it any faster. Currently, a 100 dimensional matrix takes around 40 seconds to construct, but I need around 80 of them (for different $l$ values) and I'll eventually need to scale the dimensionality up to around 400.


Nm[l_, n_] := 1/(Sqrt[π] R Abs[BesselJ[l + 1, BesselJZero[l, n]]])

En[n_, l_, mu_] := (BesselJZero[l, n])^2/(2 R^2) - mu

PairPlus[l_, n1_, n2_, delta_] := (I delta 2 π (BesselJZero[l,n2]/R) (-1)^l Nm[l + 1, n1] Nm[-l, n2] ) (-((R^2 (BesselJ[l, BesselJZero[l, n2]] BesselJ[1 + l,BesselJZero[1 + l, n1]] BesselJZero[l, n2] - BesselJ[l, BesselJZero[1 + l, n1]] BesselJ[1 + l,BesselJZero[l, n2]] BesselJZero[1 + l, n1]))/(BesselJZero[ l, n2]^2 - BesselJZero[1 + l, n1]^2)))

PairNeg[l_, n1_, n2_, delta_] := (-I delta 2 π (BesselJZero[l + 1, n2]/R) (-1)^l Nm[-l, n1] Nm[l + 1,n2]) (-((R^2 (BesselJ[-1 + l, BesselJZero[l, n1]] BesselJ[l,BesselJZero[1 + l, n2]] BesselJZero[l, n1] - BesselJ[-1 + l, BesselJZero[1 + l, n2]] BesselJ[l,BesselJZero[l, n1]] BesselJZero[1 + l, n2]))/(BesselJZero[l, n1]^2 - BesselJZero[1 + l, n2]^2)))

These are just the functions defined above with PairPlus[l_,n1_,n2_,delta_] = $\Delta_{n1,n2}^{(l)}$ and PairNeg[l_,n1_,n2_,delta_] = $\Delta_{n2,n1}^{(l)*}$. I apologize for how ugly they are but I got them from symbolic evaluation of the integrals above.

Now here is where create the matrix element $M_{i,j}$ as described above, using nested if statements

Mij[i_, j_, mu_, delta_, l_] := 
If[i == j, 
   If[EvenQ[i], En[i/2, l + 1, mu], -En[(i + 1)/2, l, mu]
   If[EvenQ[i] && OddQ[j], PairPlus[l, i/2, (j + 1)/2, delta], 
        If[EvenQ[j] && OddQ[i], PairNeg[l, (i + 1)/2, j/2, delta], 0

Basically, I first look at the diagonal entries, which depend on the parity of the row. The Else part of the first If statement then concerns off-diagonal elements, which are 0 if $i$ and $j$ have the same parity, $\Delta$ if $i$ even and $j$ odd, and $\Delta^*$ if $i$ odd and $j$ even.

With these defined, I then populate the $n$-dimensional matrix as follows

M[n_, l_, mu_, delta_] := ParallelTable[Mij[i, j, mu, delta, l], {i, 1, n}, {j, 1, n}]

A further reduction in computation time by another factor of around 10 would be fantastic, if possible!

  • $\begingroup$ Please let me know if I can improve the formatting of this post and/or provide more information. $\endgroup$
    – Aegon
    Commented Jan 8, 2016 at 1:01

1 Answer 1


You are needlessly computing exactly the same BesselJ and BesselJZero function values over and over again. As an example, in the 100x100 case, BesselJZero[1 + l, 1] (l is a constant) is computed 551 times! You just need to compute each one once. The easy way to do that is to memoize:

bj[n_?NumericQ, z_?NumericQ] := bj[n, z] = N[BesselJ[n, z]]

bjz[n_?NumericQ, k_?NumericQ] := bjz[n, k] = N[BesselJZero[n, k]]

Nm[l_, n_] := 1/(Sqrt[π] R Abs[bj[l + 1, bjz[l, n]]])

En[n_, l_, mu_] := (bjz[l, n])^2/(2 R^2) - mu

PairPlus[l_, n1_, n2_, delta_] :=
     (I delta 2 π (bjz[l, n2]/R) (-1)^l Nm[l + 1, n1] Nm[-l, n2]) \
     (-((R^2 (bj[l, bjz[l, n2]] bj[1 + l, bjz[1 + l, n1]] bjz[l, n2] -
              bj[l, bjz[1 + l, n1]] bj[1 + l, bjz[l, n2]] bjz[1 + l, n1]))/
        (bjz[l, n2]^2 - bjz[1 + l, n1]^2)))

PairNeg[l_, n1_, n2_, delta_] :=
     (-I delta 2 π (bjz[l + 1, n2]/R) (-1)^l Nm[-l, n1] Nm[l + 1, n2]) \
     (-((R^2 (bj[-1 + l, bjz[l, n1]] bj[l, bjz[1 + l, n2]] bjz[l, n1] -
              bj[-1 + l, bjz[1 + l, n2]] bj[l, bjz[l, n1]] bjz[1 + l, n2]))/
        (bjz[l, n1]^2 - bjz[1 + l, n2]^2)))

That speeds up the 100x100 case by a factor of about 1250.

  • $\begingroup$ Thanks for that! That is incredible and not something I had ever thought about. $\endgroup$
    – Aegon
    Commented Jan 8, 2016 at 5:28

Your Answer

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

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