Efficient evaluation of eigenvalues of a unitary operator formed using sparse hermitian matrices

In principle, the Arnoldi method needs only a function that evaluates U.v for a vector v. Such a function can be obtained without computing U as dense matrix first:

{Hz, Hx} = Hkfi[J, g, n];
diag = Exp[-I Normal[Diagonal[Hz]]];
A = -I Hx;
DotU = v |-> diag MatrixExp[A, v];

Note that A is sparse. So MatrixExp[A, v] can be evaluated much more quickly than U = diag MatrixExp[A] and then U. v, while providing essentially the same result. Computing U = diag MatrixExp[A] has probably complexity $O(N^3)$, where $N = 2^n$. A single A.v has probably complexity about O(N) only. I am not sure about the complexity of MatrixExp[A, v]; but experimentally, it seems to be only slightly above $O(N)$. (I guess because the exponential series converges very quickly for the arguments I put in.) In fact, what I suppose what MatrixExp[A, v] does is just evaluating a truncated exponential series of A applied to the vector v as in the following function (see also my answer here):

MatrixPolynomial[DotA_Function, coefficients_?VectorQ, v_?VectorQ] := 
 Module[{degree, Akv, result}, degree = Length[coefficients] - 1;
  Akv = v;
  result = coefficients[[1]] v;
  Do[Akv = DotA[Akv];
   result += coefficients[[k + 1]] Akv;, {k, 1, degree}];
  result
  ]

I expect that MatrixExp[A, v] does not do that with a fixed loop count, but rather by tracking an upper bound for the size of the tail of the series and by stopping when the Taylor remainder is sufficiently small. (One has to estimate the spectral norm, and that can be down with the Power Method, which also needs only powers of A applied to a fixed vector v.)

So, in a perfect world, we would define

coefficients = Table[Exp[-I m \[Phi]], {m, 0, k}];
DotgU = v |-> MatrixPolynomial[DotU, coefficients, v];

and call a matrix-free Arnoldi implemention with the function DotgU as argument. Sadly, Mathematica does not offer an API for that, at least not to my knowledge.

There are several ways around this:

1. Writing a LibaryLink function that calls ARPACK. This is a but difficult, because one has keep the ARPACK library alive while communicating with it from within Mathematica.

2. Rewriting the Arnoldi method from scratch in Mathematica as I did here. It is quite difficult to reach the same level of sophistication than ARPACK.

3. Writing a command line tools in C (or whatever you like) that calls ARPACK directly and also is able to compute MatrixPower[A,v]. Then write the sparse matrices I Hz and I Hx to file, load them into the command line program, do the eigen analysis and write the results back to file.

4. You can also start an Python session from within Mathematica and send the data over there, let the computations be performed there an then get the data back.

The latter would start with

$arnoldiPython = StartExternalSession[Association[
    "System" -> "Python",
    "Evaluator" -> <|"Dependencies" -> {"numpy", "scipy"}, 
      "EnvironmentName" -> "arnoldienv"|>
    ]];

Then you can send strings of Python code over there like this:

ExternalEvaluate[$arnoldiPython,
 "
import scipy as sp; 
import numpy as np; 
"]

You can send data, e.g., a vector v. To set a vector named "v" in Python, you can do

ExternalEvaluate[$arnoldiPython,
 Association[
  "Command" -> "
def set_fun(x) : 
    global v
    v = np.array(x)

set_fun
",
  "Arguments" -> {{1, 2, 3}}
  ]
 ]

And you can get it back with

w = Normal[ExternalEvaluate[$arnoldiPython, "v"]]

{1, 2, 3}

This works, but it is a bit awkward to define the function set_fun and then calling it. Also, I think the vector v is send as a general Python list -- not as a numpy.array. The latter would certainly be more desirable. There is certainly a better way. I just do not know it at the moment. Fortunately, there are many others on this site who are more proficient with the Python interface than I am.