Speed up construction of a dense matrix

Question

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.

---

Code:

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!

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.