1.12. Numerics

The numerical intensive calculations are encapsulated by the OperatorMatrix class, which can encode quantum states or operators in a matrix representation, meaning that each object must have two dimensions. The class was already introduced earlier, but now we would like to illuminate the computational efficiency for advanced users.

A frequently used and computationally complex operation is the calculation of the matrix potential, which is used to solve partial differential equations or first order like Schroedingers equation or a master equation in lindblad form. For example an \(X_{\pi/2}\) rotation on the bloch sphere is given by a unitary:

\begin{equation} U = e^{i \pi \sigma_x /4} \end{equation}

which can be calculated using several numeric methods. The matrix exponential and its frechet derivative are usually calculated together, for example by spectral decomposition:

import numpy as np
from qopt.matrix import DenseOperator

sigma_x = DenseOperator.pauli_x()

print(sigma_x.dexp(tau=.25j * np.pi,
                   method='spectral',
                   compute_expm=True,
                   direction=sigma_x))

(DenseOperator with data:
array([[7.07106781e-01-2.22044605e-16j, 5.55111512e-17+7.07106781e-01j],
       [5.55111512e-17+7.07106781e-01j, 7.07106781e-01+0.00000000e+00j]]), DenseOperator with data:
array([[-5.55360367e-01-1.14514159e-16j,  2.88907583e-16+5.55360367e-01j],
       [ 2.75869449e-16+5.55360367e-01j, -5.55360367e-01-1.88672737e-16j]]))

Aternatively, the scipy method expm_frechet can be used:

print(sigma_x.dexp(tau=.25j * np.pi,
                   method='Frechet',
                   compute_expm=True,
                   direction=sigma_x))

(DenseOperator with data:
array([[0.70710678+0.j        , 0.        +0.70710678j],
       [0.        +0.70710678j, 0.70710678+0.j        ]]), DenseOperator with data:
array([[-0.55536037+0.j        ,  0.        +0.55536037j],
       [ 0.        +0.55536037j, -0.55536037+0.j        ]]))

The derivative can also be calculated by finite differences using the method approx:

print(sigma_x.dexp(tau=.25j * np.pi,
                   method='approx',
                   compute_expm=True,
                   direction=sigma_x))

(DenseOperator with data:
array([[7.07106781e-01-2.22044605e-16j, 5.55111512e-17+7.07106781e-01j],
       [5.55111512e-17+7.07106781e-01j, 7.07106781e-01+0.00000000e+00j]]), DenseOperator with data:
array([[-5.55359092e-01+2.22044605e-06j, -5.55111512e-07+5.55362423e-01j],
       [-5.55111512e-07+5.55362423e-01j, -5.55361312e-01+0.00000000e+00j]]))

And we can also explicitly calculate a number of terms in the expansion serias of the exponential function:

print(sigma_x.dexp(tau=.25j * np.pi,
                   method='first_order',
                   compute_expm=True,
                   direction=sigma_x))

(DenseOperator with data:
array([[7.07106781e-01-2.22044605e-16j, 5.55111512e-17+7.07106781e-01j],
       [5.55111512e-17+7.07106781e-01j, 7.07106781e-01+0.00000000e+00j]]), DenseOperator with data:
array([[0.+0.j        , 0.+0.78539816j],
       [0.+0.78539816j, 0.+0.j        ]]))

print(sigma_x.dexp(tau=.25j * np.pi,
                   method='second_order',
                   compute_expm=True,
                   direction=sigma_x))

(DenseOperator with data:
array([[7.07106781e-01-2.22044605e-16j, 5.55111512e-17+7.07106781e-01j],
       [5.55111512e-17+7.07106781e-01j, 7.07106781e-01+0.00000000e+00j]]), DenseOperator with data:
array([[-0.61685028+0.j        ,  0.        +0.78539816j],
       [ 0.        +0.78539816j, -0.61685028+0.j        ]]))

print(sigma_x.dexp(tau=.25j * np.pi,
                   method='third_order',
                   compute_expm=True,
                   direction=sigma_x))

(DenseOperator with data:
array([[7.07106781e-01-2.22044605e-16j, 5.55111512e-17+7.07106781e-01j],
       [5.55111512e-17+7.07106781e-01j, 7.07106781e-01+0.00000000e+00j]]), DenseOperator with data:
array([[-0.61685028+0.j        ,  0.        +0.54316163j],
       [ 0.        +0.54316163j, -0.61685028+0.j        ]]))