Introduction

For the non-commutative operators $ X and Y $, the following Lie-Trotter formula holds [^ name].

\exp(h(X+Y)) = \left( \exp(hX/n)\exp(hY/n) \right)^n + O\left(\frac{h^2}{n}\right)

Here, $ h $ is a c-number and is often used as a time step in numerical calculation, so the following is referred to as a time step. $ n $ is the factorization number. The purpose of this article is to confirm this censoring error.

The source is https://github.com/kaityo256/lie-trotter-sample At.

Added on July 14, 2017: The source has been slightly modified so that you can also try quadratic symmetric decomposition. For details, refer to the following article Second-order symmetric decomposition in the Lie-Trotter formula.

Explanation of what the script is doing

Consider a square matrix $ X $, $ Y $ of the appropriate $ d $ dimension. Let's make it a symmetric matrix so that the eigenvalues are real. Strictly calculate $ \ exp (h (X + Y)) $, $ \ exp (hX) $, $ \ exp (hY) $, etc., and calculate the error using the Lie-Trotter formula. For the error, we will use the Frobenius norm.

Make a real symmetric matrix

First, randomly create a symmetric matrix of appropriate $ d $ dimension.

import numpy as np
import math
from scipy import linalg

d = 2
np.random.seed(1)
x1 = np.random.rand(d,d)
x2 = np.random.rand(d,d)
X = x1.dot(x1.T)
Y = x2.dot(x2.T)
Z = X + Y

Make an appropriately random matrix and make the product of the transposed ones $ X $ or $ Y $. In this way, the eigenvalue becomes real because it becomes a symmetric matrix. Let's assume that the sum of $ X $ and $ Y $ is $ Z $.

Calculation of matrix exponential function

Create a function that puts the matrix on the shoulder of the exponent at the specified time interval $ h $. To do so, first diagonalize.

X = U_x \Lambda_x U_x^{-1}

Because it's easy to put the diagonal matrix on the shoulder of the exponential function

\exp(h X) = U_x \exp(h \Lambda_x) U_x^{-1}

Can be calculated. Writing this in Python looks like this.

def calc_exp(X, h):
    rx, Ux = linalg.eigh(X)
    Ux_inv = linalg.inv(Ux)
    tx = np.diag(np.array([math.exp(h * i) for i in rx]))
    eX = Ux.dot(tx.dot(Ux_inv))
    return eX

Explain what you are doing

Call linalg.eigh because it feeds on a symmetric matrix
Put a vector of eigenvalues on the shoulder of an exponential function in a list comprehension
Matrix it with numpy.diag
Undo with $ U_x $ or $ U_x ^ {-1} $

Feeling like that.

Specified number of decompositions Trotter decomposition

It will be easy if the matrix exponential function can be calculated in the specified time interval. For example, in the case of $ h $ in time increments and $ n $ in factorization, first make $ \ exp (h X / n) $ and $ \ exp (h Y / n) $. All you have to do is multiply the identity matrix $ n $ times.

    eX = calc_exp(X,h/n)
    eY = calc_exp(Y,h/n)
    S = np.diag(np.ones(d))
    eXeY = eX.dot(eY)
    for i in range(n):
        S = S.dot(eXeY)

Compare the approximate matrix $ S $ thus created with the exact evaluation $ \ exp (h Z) $. Let's compare with the Frobenius norm here. Putting all the above together makes a function like this.

def trotter(X,Y,Z,h,n):
    eZ = calc_exp(Z,h)
    eX = calc_exp(X,h/n)
    eY = calc_exp(Y,h/n)
    S = np.diag(np.ones(d))
    eXeY = eX.dot(eY)
    for i in range(n):
        S = S.dot(eXeY)
    return linalg.norm(eZ - S)/linalg.norm(eZ)

result

In the two-dimensional case, let's change $ h $ and $ n $ to see what happens to the truncation error.

First, the $ h $ dependency of the censoring error. Let $ n $ be 1,2,4, and evaluate the censoring error with various $ h $ values for each.

It can be seen that the $ h $ dependency of the error is $ (h ^ 2) $ regardless of the value of $ n $. The larger the $ n $, the smaller the error, which is proportional to $ 1 / n $. Fix it at $ h = 1.0 $ and evaluate the censoring error for various $ n $.

It can be seen that the error is certainly $ O (1 / n) $.

Summary

I evaluated the censoring error of the trotter decomposition. Until now, the case $ \ exp (h (A + B)) \ sim \ exp (hA) \ exp (hB) $ of $ n = 1 $ was called "first-order decomposition", so somehow the error is $ O. I thought it was (h) $, but now I know that it was $ O (h ^ 2) $ regardless of the value of $ n $.

[^ name]: Some scary uncles may come to say the right thing about this official name. For the time being, if the operator is bounded, the Lie formula, if it is not bounded, the Trotter formula, and the algorithm that integrates the path using this is called Suzuki-Trotter decomposition. I often call it simply trotter disassembly, but ...

[PYTHON] Examine the Lie-Trotter formula censoring error