Matrix free arnoldi

I want to do a "matrix-free eigenvalue evaluation" of a unitary matrix using ARPACK package. Unitary matrix comes from a sparse hermitian matrix, something that happens in many-body problem.

I used the following code to construct a random sparse matrix ( here it is small size for reason that the matrix free version takes exponentially long time as compared to this version for evaluation)

import time
import numpy as np
import scipy as sp
from scipy.sparse.linalg import LinearOperator
from scipy.sparse.linalg import expm, expm_multiply  # For matrix exponentials
from scipy.sparse.linalg import eigs, eigsh
from scipy.stats import unitary_group
# Generate a random unitary matrix
N = 20  # Size of the matrix
den = 0.01

# Define Hermitian matrices H1 and H2 (only stored implicitly)
H1 = sp.sparse.random_array((N,N), density=den)
H1 = (H1 + H1.conj().T) / 2 
H2 = sp.sparse.random_array((N,N), density=den)
H2 =(H2 + H2.conj().T) / 2 

U= expm(-1j * H2) @ expm(-1j * H1)

start = time.time()
# Compute eigenvalues using scipy
eigenvalues,_= eigs(U,5,sigma= np.exp(1j * np.pi / 4), which = 'SR')
end = time.time()
print("The time of execution of dense matrix program is :",(end-start) * 10**3, "ms")
# Print the results
print("Eigenvalues of the unitary matrx:")
print(eigenvalues)

Now I have the unitary matrix as $U= e^{-i H_{2}} e^{-i H_{1}}$. To have it in a matrix free form I wrote the following code and used ARPACK to find the eigenvalues

def apply_U(v):
    v= expm_multiply(-1j * H2, expm_multiply(-1j* H1,v))
    return v   
A = LinearOperator((N, N), matvec=apply_U)

start = time.time()
# Compute eigenvalues using scipy
eigenvalues,_= eigs(A,5,sigma= np.exp(1j * np.pi / 4), which = 'SR')
end = time.time()
print("The time of execution of matrix-free is :",(end-start) * 10**3, "ms")
# Print the results
print("Eigenvalues of the unitary matrx:")
print(eigenvalues)

I have checked that both of them give the same result except the matrix-free version is very slow. My expectations were that it would be faster. Please help me optimize the code, as I am new to python ( I coded this all with my intuition from Mathematica in which I am fluent). I would like to know if it is an expected behavior from the matrix free evaluations or is there any problem with my code specifically, thus making it slow

Listen to your warnings. The default from your random generation is CSR but the warnings tell you you should convert to CSC:

import functools
import time

import numpy as np
import scipy.sparse as sps
import scipy.sparse.linalg as spl


def generate(
    rand: np.random.Generator, den: float = 0.01,
    N: int = 20,  # Size of the matrix
) -> tuple[sps.csc_array, sps.csc_array]:
    # Generate a random unitary matrix
    # Define Hermitian matrices H1 and H2 (only stored implicitly)
    H1 = sps.random_array(shape=(N, N), density=den, rng=rand, format='csc')
    H1 += H1.conj().T
    H2 = sps.random_array(shape=(N, N), density=den, rng=rand, format='csc')
    H2 += H2.conj().T
    return H1*0.5, H2*0.5


def eig_method(H1: sps.csc_array, H2: sps.csc_array):
    # Compute eigenvalues using scipy
    U = spl.expm(-1j * H2) @ spl.expm(-1j * H1)
    eigenvalues, _ = spl.eigs(U, k=5, sigma=np.exp(0.25j * np.pi), which='SR')
    return eigenvalues


def apply_U(
    v: np.ndarray,
    nH1: sps.csc_array, nH2: sps.csc_array,
    trace1: np.complex128, trace2: np.complex128,
) -> np.ndarray:
    return spl.expm_multiply(
        nH2,
        spl.expm_multiply(nH1, v, traceA=trace2),
        traceA=trace1,
    )


def matrix_free(H1: sps.csc_array, H2: sps.csc_array) -> np.ndarray:
    N, N = H1.shape
    nH1 = -1j*H1
    nH2 = -1j*H2
    trace1 = nH1.trace()
    trace2 = nH2.trace()
    A = spl.LinearOperator(
        shape=(N, N),
        matvec=functools.partial(apply_U, nH1=nH1, nH2=nH2, trace1=trace1, trace2=trace2),
    )
    # Compute eigenvalues using scipy
    eigenvalues, _ = spl.eigs(A, k=5, sigma=np.exp(0.25j * np.pi), which='SR')
    return eigenvalues


def main() -> None:
    rand = np.random.default_rng(seed=0)
    H1, H2 = generate(rand=rand)

    t0 = time.perf_counter()
    ea = eig_method(H1, H2)
    t1 = time.perf_counter()
    eb = matrix_free(H1, H2)
    t2 = time.perf_counter()
    print(t1 - t0)
    print(t2 - t1)

    ea.sort()
    eb.sort()
    assert np.allclose(ea, eb, rtol=0, atol=1e-14)
    assert np.allclose(
        ea,
        [
            0.76873299 - 0.63956985j, 0.91509182 - 0.40324554j,
            0.91845690 - 0.39552108j, 0.95435119 - 0.29868681j,
            0.99985900 - 0.01679200j,
        ],
        rtol=0, atol=1e-8,
    )


if __name__ == '__main__':
    main()