User-defined Diagonalization#
Often the most time consumuing operation that scqubits
performs is a diagonalization of the system Hamiltonian during calls to eigensys
or eigenvals
. scqubits
allows users to specify what library or procedure is used, if something other than the default method (typically provided via the scipy
library) is desired.
Such customization is done by setting esys_method
and evals_method
(and potentially esys_method_options
and evals_method_options
to further customize diagonalization options) of the different quantum system classes (which can be done at object initialization time, or alternatively, after the object has already been created).
Furthermore, esys_method
and evals_method
can be either:
A string that
scqubits
maps to a predefined diagonalization function, with some particular set of preset options.An arbitrary, user-defined callable object (e.g., a function) that can perform the diagonalization.
Parameters esys_method_options
and evals_method_options
can be also set, if the user wishes to pass custom options to the diagonalization routines defined by esys_method
and evals_method
respectively. Setting these is a powerful way to easily customize how the diagonalization is performed.
Warning
At the moment custom diagonalization is not supported by the Circuit module. This functionality will be added in a future version of scqubits
, however.
Predefined procedures#
In the case of scenario (1) above, users can check all the predefined diagonalization procedures that have a string-representation by listing the keys of the DIAG_METHODS
dictionary:
[2]:
print(*scq.DIAG_METHODS.keys(), sep="\n")
evals_scipy_dense
esys_scipy_dense
evals_scipy_sparse
esys_scipy_sparse
evals_scipy_sparse_SM
esys_scipy_sparse_SM
evals_scipy_sparse_LA_shift-inverse
esys_scipy_sparse_LA_shift-inverse
evals_scipy_sparse_LM_shift-inverse
esys_scipy_sparse_LM_shift-inverse
evals_primme_sparse
esys_primme_sparse
evals_primme_sparse_SM
esys_primme_sparse_SM
evals_primme_sparse_LA_shift-inverse
esys_primme_sparse_LA_shift-inverse
evals_primme_sparse_LM_shift-inverse
esys_primme_sparse_LM_shift-inverse
evals_cupy_dense
esys_cupy_dense
evals_cupy_sparse
esys_cupy_sparse
evals_jax_dense
esys_jax_dense
In the above, the first part of the string represents whether just the eigenvalues or both eigenvalues and eigenvectors are returned. Naturally, procedures that return the eigenvalues (eigenvalues and eigenvectors) are used by the method eigenvals
(eigensys
) of the different quantum objects that scqubits
implements. The second word labels what library is used (e.g., scipy
or cuda
, etc.), and the third, whether sparse or dense diagonalization should be performed.
Everything that follows, is simply some description of specific diagonalization options that are used.
For example, a method esys_cupy_sparse
will invoke the diagonalization function eigsh from the sparse implemetnation of the cupy library (see here), which will return the lowest set of requested eigenvalues and eigenvectors (i.e., will have the options which='SA'
passed in).
On the other hand, if the user sets esys_method
to "esys_scipy_sparse_LA_shift-inverse"
when initializing a, for example, Transmon
object, will results in scqubits
internally invoking the scipy
’s sparse diagonalization routine (here eighs
), whenever the Transmon.eigensys
method is called. Furthermore, the shift-inverse method will be used, with algebraic (and not by magnitude) eienvalue ordering (i.e., which='LA'
will be passed to eighs
- see scipy
’s
documnetation for more details).
Warning
It is crucial to note that in some cases, a particular method may be the wrong choice. For example in quantum systems where the eigenenergies are not guaranteed to be positive, using the “SM” method (that sorts eigenvalues by magnitude), may lead to incorrect results.
scqubits
supports a number of backend libraries such as, scipy, primme, cupy (which allows for GPU-based diagonalization), as well as jax (which also, optionally, allows for GPU-based operations).
Note
Some of the libraries that can be used for diagonalization (e.g., cupy or primme) are only optional depedenecies of scqubits
. Naturally, however, they do need to be installed before a digonalization method that relies on them is selected and used.
As an example, let us consider a fluxonium qubit with a custom diagonalization method (and options set by providing a dictionary to evals_method_options
), that take advante of the sparase diagonalization routine provided by the primme library (which is assumed to be installed on the system).
[3]:
fluxonium= scq.Fluxonium(EJ=8.9, EC=2.5, EL=0.5, flux = 0.5, cutoff = 120,
evals_method='evals_primme_sparse',
evals_method_options=dict(tol=1e-4) # custom tolerance
)
evals = fluxonium.eigenvals()
We can compare the eigenvalues obtained useing the primme library, to the default scqubits
method (by resettting evals_method
and evals_method_options
to None
) before calling eigenvals
for the second time.
[4]:
fluxonium.evals_method = None
fluxonium.evals_method_options = None
evals_default = fluxonium.eigenvals()
print(evals_default)
[-0.63001146 -0.26659009 8.75633751 11.69881054 16.67252211 17.42908833]
Comparing the two approaches shows the same results
[5]:
np.allclose(evals, evals_default)
[5]:
True
Custom procedures#
As outlined above, the diagonalization procedures can be even further customized. This is done by, instead of passing a string-based description as esys_method
or evals_method
, providing a callable object (e.g., function). In principle this approach could be used to take advantage of 3rd party libraries (e.g., that use GPUs) that may not have yet been directly utilized in scqubits
. In the following trivial example, a simple user-defined function custom_esys
(that internally calls
scipy.linalg.eigh
) is used for the diagonalization.
[6]:
def custom_esys(matrix, evals_count, **kwargs):
evals, evecs = sp.linalg.eigh(
matrix, subset_by_index=(0, evals_count - 1),
eigvals_only=False,
**kwargs
)
return evals, evecs
tmon = scq.Transmon(EJ=30.02, EC=1.2, ng=0.0, ncut=501,
esys_method=custom_esys,
esys_method_options=dict(driver='evr') # here we set a custom driver
)
We can always check what all the diagonalization related variables in a given quantum system (here a tmon
) object are set to
[7]:
print(tmon.esys_method,
tmon.esys_method_options,
tmon.evals_method,
tmon.evals_method_options, sep="\n")
<function custom_esys at 0x7fea5477d5e0>
{'driver': 'evr'}
None
None
We can next diagonalize the system with a standard call to eigensys
, which will return both the eigenvalues as well as the eigenvectors, but using the custom_esys
method defined above
[8]:
evals, evecs = tmon.eigensys()
Let us compare the above results with the default diagonalization method (which in the case of the transmon qubit is done using scipy.linalg.eigvalsh_tridiagonal
):
[9]:
tmon.esys_method = None
tmon.esys_method_options = None
evals_default, evecs_default = tmon.eigensys()
Comparing the eigenvalues obtained using the two procedures gives
[10]:
np.allclose(evals, evals_default)
[10]:
True
and for the eigenvecotrs:
[11]:
print(*[np.allclose(evecs[:,i], evecs_default[:,i]) for i,_ in enumerate(evals)])
True True True True True True
Diagonalization of composite systems#
Custom diagonaliation procedure can set by the user for both, the predefined qubit classes (as shown above), but also for systems consisting of joint Hilbert spaces. Since often one is intersted in a low-lying energy landscape of larger systems, scqubits
uses the concept of hierarchical diagonalization to minimize the size of the relevant Hilber space that has to be considered: Each system is first diagonalized separately, the coupling terms
are then rewritten in the (often some subset of) repsective eigenbasis, and finally the combained sytem is diagonalized at the end. This gives a lot of control to the user, and scqubits
allows one to set custom (and even different) diagonalization procedure for each step of this process.
Here as an example of a joint system consisting of a Fluxonium and a Transmon qubit, where each qubit, as well as the combined Hamiltonian are all diagonalized using a different method.
[12]:
tmon = scq.Transmon(EJ=30.02, EC=1.2, ng=0.0, ncut=101, ) # default diaongalization method
fluxonium= scq.Fluxonium(EJ=8.9, EC=2.5, EL=0.5, flux = 0.5, cutoff = 110,
evals_method='evals_primme_sparse' # use the primme library
)
hs = scq.HilbertSpace([tmon, fluxonium],
evals_method='evals_scipy_sparse' # use scipy
)
hs.add_interaction(
g_strength = 0.250,
op1 = tmon.n_operator,
op2 = fluxonium.n_operator,
add_hc = True
);
[13]:
evals = hs.eigenvals(evals_count=10)
Let us reset the diagonalization methods to default, for both the fluxonium
as well as hs
:
[14]:
fluxonium.evals_method=None
hs.evals_method=None
and compare the results results with those obtained above
[15]:
evals_default = hs.eigenvals(evals_count=10)
np.allclose(evals, evals_default)
[15]:
True
Note that regardless of the (predefined) diagonalization procedure that has been used, the eigenvectors that are returened, are objects of the type qutip.Qobj
, and hence remember their Hilbert space (i.e. tensor product) structure.
[16]:
hs.esys_method='esys_primme_sparse'
evals, evecs = hs.eigensys()
[17]:
evecs[3]
[17]:
GPU support#
Diagonalization on GPUs is currently supported through the cupy, and optionally, the jax libraries. As already outlined above, to use user any of these packages for diagonalization, they have to be installed on the user’s machine.
It is worth pointing out, that depending on one’s hardware, and the quantum system studied, the performance of GPU vs CPU diagonalization may vary dramatically from setup to setup.
Warning
Using multiple GPU-reliant libraries within the same python
session may not always be supported by either the libraries, or the hardware vendors (e.g. GPUs may need resetting and/or reinitialization between uses). Users may need to consult documentation on details of how that is done in their setup.
GPU use with cupy
#
We can check that the cupy
library is installed on the current machine, by importing it
[7]:
import cupy
cupy.__version__
[7]:
'12.1.0'
To diagionalize using cupy
, we can set the esys_method
and/or evals_method
as in the following
[4]:
zero_pi = scq.ZeroPi(
grid = scq.Grid1d(-6*np.pi, 6*np.pi, 200),
EJ = 10.0,
EL = 0.04,
ECJ = 20.0,
EC = 0.04,
ng = 0.1,
flux = 0.23,
ncut = 5, #smaller than needed, for testing purposes
esys_method='esys_cupy_dense',
evals_method='evals_cupy_dense',
)
[5]:
fig, ax = zero_pi.plot_wavefunction(which=1, mode='real', zero_calibrate=True);
[6]:
flux_list = np.linspace(0, 1, 27)
zero_pi.plot_evals_vs_paramvals('flux', flux_list, evals_count=8);
GPU use with jax
#
The jax
library supports both CPU, and well as GPU based backends. The default on this machine can be checked with
[23]:
import jax
jax.devices()
[23]:
[gpu(id=0)]
Setting esys_method
and/or evals_method
to one of the jax
-based routines, will use the backed that is currently set, which here is GPU-based
[29]:
fluxonium= scq.Fluxonium(EJ=8.9, EC=2.5, EL=0.5, flux = 0.5, cutoff = 120,
evals_method='evals_jax_dense'
)
evals = fluxonium.eigenvals()