Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

Chebyshev.cpp File Reference

Source code for Chebyshev polynomial approximations. More...

#include "linalg.h"
#include "quantum.h"
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_errno.h>

Include dependency graph for Chebyshev.cpp:


void GenerateChebyshevMatrix (itype N, DenseMatrix< complex > &M)
 Computes an NxN matrix, where $ M_{m,n} $ is the coefficient of $ x^m $ in the $ n $ 'th Chebyshev polynomial.
void ResizeChebyshevMatrix (itype N, DenseMatrix< complex > &M)
 Efficiently resizes an existing Chebyshev polynomial matrix to NxN, using the coefficients already computed.
void GenerateChebyshevEvolutionVector (double time, itype N, DenseVector< complex > &V)
 Computes an N-element vector, where $ v_n $ is the coefficient of the $ n $ 'th Chebyshev polynomial in the expansion of $ e^{-ix} $ .
int ChebyshevPrep (double timestep, double prec, DenseVector< complex > &PowerCoeffs, DenseVector< double > &CoeffsBoundSeq)
 Precomputes the coefficients needed to evolve a state using ChebyshevStep(), for time timestep, with accuracy prec.

Detailed Description

Source code for Chebyshev polynomial approximations.

A standard and nontrivial task in quantum simulations is to compute the unitary evolution of a state vector:

\[ \left|\psi\right\rangle \longrightarrow \left|\psi(t)\right\rangle = e^{-i\mathbf{H}t}\left|\psi\right\rangle \]

Matrix exponentiation is nontrivial, especially if the matrix in question is very large. However, computing $ \left|\psi(t)\right\rangle $ does not actually require computing $ e^{-i\mathbf{H}t} $ , or even $ \mathbf{H} $ . It can be written as a power series in $ \mathbf{H} $ :

\[ e^{-i\mathbf{H}t}\left|\psi\right\rangle \approx \left(\mathbf{1} -i\mathbf{H}t - \frac{t^2}{2}\mathbf{H}^2\ldots\right) \left|\psi\right\rangle \]

This requires only that we compute $ \mathbf{H}^n\left|\psi\right\rangle $ , which in turn can be done recursively if we can compute $ \left|\psi\right\rangle \longrightarrow \mathbf{H}\left|\psi\right\rangle$ .

The Taylor expansion given above turns out to converge rather slowly; the Taylor series is optimized for local convergence at $ t=0 $ , at the cost of all larger $ t $ . Other power series expansions exist, and have better global convergence properties. The power series with the best global convergence properties over a fixed intervals is the Chebyshev series:

This makes Chebyshev polynomials ideally suited for quantum state evolution. When they are generalized to matrices, the argument $ x $ becomes a matrix, which typically has many eigenvalues. We want the approximation of $ f(x) $ to be accurate for all the eigenvalues. This motivates the uniform convergence of the Chebyshev polynomial expansion. It is necessary to bound the highest and lowest eigenvalues of $ \mathbf{H} $ , and then rescale energy and time so that $ ||\mathbf{\tilde{H}}||_\infty < 1 $ , but once this is done, computing the correct power series is straightforward. From the rescaled time $ \tilde{t} $ , it is possible to identify how many coefficients will be needed to approximate the matrix exponential well; beyond a critical value, the error declines exponentially.

Representing $ e^{-i\mathbf{H}t} $ as a power series in $ \mathbf{H} $ requires first computing the Chebyshev polynomials themselves. The method GenerateChebyshevMatrix() does this via a simple recursion relation, storing the coefficient of $ x^m $ in $ T_n(x) $ as an element $ M_{mn} $ in a (triangular) matrix. If more coefficients are needed later, ResizeChebyshevMatrix() can compute just the additional ones needed.

The second step involves computing the expansion of $ e^{-i\tilde{t}x} $ in Chebyshev polynomials. The method GenerateChebyshevEvolutionVector() does this using Bessel functions (borrowed from the Gnu Scientific Library), and stores the expansion in a complex vector. The coefficients of $ \mathbf{\tilde{H}}^n $ in the approximated matrix exponential can then be computed by matrix-vector multiplication.

For a particular Hamiltonian and time, ChebyshevPrep() first determines what order approximation will be required to bound the error by some very small number, then computes the coefficients. A single step of duration $ t=1 $ can then be accomplished by ChebyshevStep(). The entire process of evolving for time $ t $ (including both ChebyshevPrep() and ChebyshevStep()), can be done in one step using ChebyshevEvolve().

However, taking a single step of large duration can lead to substantially increased error -- not because of error in the Chebyshev polynomials, but because the coefficients become so large for high-order approximations that they amplify rounding errors. This could very likely be taken care of by a more sophisticated method of evaluating polynomials, but this would require keeping additional temporary vectors (e.g., $ \mathbf{\tilde{H}}^n\left|\psi\right\rangle $ for several values of $ n $ ), which could be memory-intensive for extremely large vectors. As an alternative, ChebyshevStepEvolve() will automatically break up large timesteps into several smaller jumps of bounded size.

While the Chebyshev polynomial method is (unlike split-operator or Cayley methods) not explicitly norm-preserving, its accuracy is so absurdly good (errors are typically $ O(10^{-13}) $ in double precision) that unitarity is not a concern.

Generated on Wed Jun 14 22:25:26 2006 for linalg by  doxygen 1.4.4