# Computing matrix exponential with PETSc/SLEPc

My question is basically the best way to compute the matrix exponential in parallel. I've got a sequential code written in Python (https://github.com/hadsed/AdiaQC) and I need to extend it to run in a high performance environment (cluster, supercomputer, maybe GPU).

The code simulates adiabatic quantum computation, and the main things to worry about is a complex-valued vector representing the quantum state, and a Hamiltonian matrix which I do an eigendecomposition on, as well as an exponent for the matrix exponential operator that I apply to the state vector. Hamiltonian and exponent are dense, and it grows as $2^n\times 2^n$, and state vector of course is $2^n$ as well. Here $n$ is the number of qubits.

There seems to be some development in SLEPc for a matrix exponential module, but I don't think I can wait for it.

I'm thinking that the eigendecomposition method for the matrix exponential probably isn't the most efficient, but I'm not exactly sure how to efficiently do the other methods (and if they really are better). Any suggestions? I have gotten petsc4py and slepc4py working finally, so I'd like to exploit that as much as I can.

Let me know if I need to be more specific about the actual problem.

This is my Ising Hamiltonian: $H = \sum\limits_{ij} \beta_{ij} z_i z_j + \sum\limits_{i} \alpha_i z_i + \sum\limits_i \delta_i x_i$ where

$$z_i = I \otimes ... \otimes ~ \sigma^z_i ~ \otimes ... \otimes I$$

$\sigma^z$ is the Pauli-z matrix, 2x2, and basically the identity but the lower diagonal element is -1. $I$ is just the identity, 2x2. So $\sigma^z_i$ is at the $i$th position of this bunch of Kronecker products. We cycle through i := {0, 1, ... , n} where n is the number of qubits. You do the same thing for $x_i$ (this time with the Pauli-x matrix, which is a 2x2 matrix, ones on the antidiagonal. With the two $z$s you just shift the position of two Pauli-z matrices instead of one in this big Kronecker product thing and sum those.

But this Hamiltonian is time-dependent. $\alpha$, $\beta$, and $\delta$ are the only time-dependent components. Actually, they are coefficient matrices/vectors, but they are multiplied by something that makes it time-dependent. So,

$$H_{exponent} = \frac{t}{T} H_{\alpha \beta} + (1 - \frac{t}{T}) H_{\delta}$$

where $H_{\alpha \beta}$ is just the first two sums in the Ising Hamiltonian, and same for the $H_{\delta}$. This is what I take the exponential of to apply to my state vector $\psi$ (I'm slightly lying here to simplify, it's not as simple as $t/T$ but it doesn't matter, they're just numbers that change with time).

Something to notice is that $H_{\alpha \beta}$ is completely diagonal. Really convenient, but unfortunately $H_{\delta}$ is completely the opposite of diagonal (read: dense as hell). In the beginning, $H_{exponent} = H_{\delta}$, and at the end it's $H_{\alpha \beta}$.

Anyway, explaining all this in the hopes that maybe someone can spot some exploitable structure.

Does your dense Hamiltonian have structure (e.g., sparse plus low-rank) or is it unstructured dense (no compressed representation available)? Is it well-conditioned? You can use Elemental, either on its own or through PETSc's interface, for the parallel dense linear algebra. I'm not aware of a suitable library implementation of the exponential, but you can build one out of the components provided in these libraries. If you end up developing a general-purpose parallel implementation, please consider contributing it to PETSc.

• The "Computing the Action..." algorithm has been implemented in the development version of scipy, as scipy.sparse.linalg.expm_multiply(). May 16, 2013 at 18:33
• Niesen and Wright also have a paper out on Krylov subspace methods for calculating the matrix exponential. May 16, 2013 at 20:00
• Elemental actually contains routines for real and complex functions of Hermitian matrices, e.g., see this driver, but, if only the action of the matrix exponential is required, it is likely best to avoid explicitly forming it. May 16, 2013 at 20:20
• @JackPoulson He wants $e^{i H} \Psi$ (not Hermitian). Anyway, I thought full eigendecomposition had fallen out of favor, but considering that he says he has an eigen-decomposition, it is a trivial amount of work. May 16, 2013 at 21:06
• @JackPoulson Cool, I didn't know you had an example of that. I didn't mean to imply that it was deeply different. If he wants to stick to Python, I guess there will be a couple hoops to jump through to get support for Elemental eigensolvers. One of us should contribute an interface to SLEPc so that he can use slepc4py. ;-) May 16, 2013 at 21:59

The development version of scipy has recently added expm_multiply which computes the action of the matrix exponential without explicitly computing the matrix exponential itself. This uses the algorithm in "Computing the Action of the Matrix Exponential..."
http://eprints.ma.man.ac.uk/1591/
https://github.com/scipy/scipy/pull/2456

The implementation has been recently patched to work with complex numbers, but it is less tested for this case.
https://github.com/scipy/scipy/pull/2495

Its application to your problem would look something like this.
• OK. First, imagine that $\alpha$, $\beta$, and $\delta$ are a vector (or matrix in case of $\beta$) of ones. For n = 3, we can calculate $\sum\nolimits_i \alpha \sigma^z_i$ by doing $\alpha_1(\sigma^z \otimes I \otimes I) + \alpha_2 (I \otimes \sigma^z \otimes I) + \alpha_3 (I \otimes I \otimes \sigma^z)$. Same concept for $\beta$ but with two $\sigma^z$s rotating in that sum of kron-products. For the $\sum\delta_i x_i$ part, you just do the same thing but now you use $\sigma^x$ (en.wikipedia.org/wiki/Pauli_matrices). So once I get that matrix, I look at $H_{exponent}$, because to May 17, 2013 at 3:16
• I need to approximate the time-evolution operator, and that is $exp(H_{exponent})$ and apply it to state vector $\psi$. $H_{exponent}$ changes in time though. Here is the exponent: github.com/hadsed/AdiaQC/blob/master/solve.py#L47 And what I call $\alpha$ and such there is the matrices that involve the $\alpha$ coefficients that I was mentioning before (i.e. $\alpha_{code} = \sum\limits_i \alpha_i z_i$). Same with the others. May 17, 2013 at 3:25