5
$\begingroup$

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.

EDIT: Thought it might be helpful to explain more about Hamiltonian.

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.

Improve this question
$\endgroup$
1

2 Answers 2

Reset to default
6
$\begingroup$

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.

See also these recent papers:

Improve this answer
$\endgroup$
9
  • $\begingroup$ The "Computing the Action..." algorithm has been implemented in the development version of scipy, as scipy.sparse.linalg.expm_multiply(). $\endgroup$
    – clipper
    May 16, 2013 at 18:33
  • $\begingroup$ Niesen and Wright also have a paper out on Krylov subspace methods for calculating the matrix exponential. $\endgroup$
    – Geoff Oxberry
    May 16, 2013 at 20:00
  • 1
    $\begingroup$ 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. $\endgroup$
    – Jack Poulson
    May 16, 2013 at 20:20
  • 1
    $\begingroup$ @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. $\endgroup$
    – Jed Brown
    May 16, 2013 at 21:06
  • 1
    $\begingroup$ @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. ;-) $\endgroup$
    – Jed Brown
    May 16, 2013 at 21:59
5
$\begingroup$

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.
https://github.com/hadsed/AdiaQC/pull/2

Edit:
Even if the matrix exponent is not sparse, the expm_multiply function can still be accelerated if you have a clever way of quickly computing exponent-vector products, for example if the exponent has been factored or decomposed in certain ways. I don't know enough quantum mechanics to know whether or not this would apply to your situation.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Thanks for that. I do need to do a distributed memory version of my code though. I was going to try and implement scipy.linalg.expm with the Pade approx. in petsc4py, but maybe I will try doing this version instead. I also added some more information about my exponent, so maybe some structure would be apparent to those more knowledgeable than I. By the way, approximately how long does it take for features in scipy-dev to get released? $\endgroup$
    – hadsed
    May 17, 2013 at 1:22
  • $\begingroup$ I still can't quite decipher the notation, maybe I would understand it better if there were a small example, like n=3. Scipy doesn't seem to have a fixed release schedule, and I am not one of the main scipy devs. Based on the past release history of scipy, I would guess that version 0.13 will be released sometime between 4 months to 2 years in the future. $\endgroup$
    – clipper
    May 17, 2013 at 2:24
  • $\begingroup$ 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 $\endgroup$
    – hadsed
    May 17, 2013 at 3:16
  • $\begingroup$ 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. $\endgroup$
    – hadsed
    May 17, 2013 at 3:25
  • $\begingroup$ The way your code is currently written, it seems like you cannot depend on any particular structure because the matrices specified in github.com/hadsed/AdiaQC/tree/master/problems could be arbitrary. I guess that you could also allow the interface to include matrix-vector product functions to be specified in the problems/*.py files. In any case, expm_multiply(A, v) is still faster than expm(A)*v even when A is dense and unstructured. $\endgroup$
    – clipper
    May 17, 2013 at 14:15

Your Answer

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

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