-
Notifications
You must be signed in to change notification settings - Fork 1
/
lanczos.py
120 lines (94 loc) · 3.51 KB
/
lanczos.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
"""
Module containing Lanczos diagonalization and time evolution algorithm
"""
from typing import Union
import numpy as np
from scipy.sparse import csr_matrix
def lanczos_basis(array: Union[csr_matrix, np.ndarray], v_0: np.ndarray, k_dim: int):
"""Tridiagonalises a hermitian array in a krylov subspace of dimension k_dim
using Lanczos algorithm. returns basis of krylov subspace and projection of array
in krylov subspace
Args:
array : Array to tridiagonalise.
v_0 : Initial state.
k_dim : Dimension of the krylov subspace.
Returns:
q_basis : Basis of the krylov subspace.
tridiagonal : Tridiagonal projection of ``array``.
"""
error = np.finfo(np.float64).eps
data_type = np.result_type(array.dtype, v_0.dtype)
v_0 = np.array(v_0).reshape(-1, 1)
v_0 = v_0 / np.linalg.norm(v_0)
array_dim = array.shape[0]
q_basis = np.zeros((k_dim, array_dim), dtype=data_type)
v_p = np.zeros_like(v_0)
projection = np.zeros_like(v_0)
beta = np.zeros((k_dim,), dtype=data_type)
alpha = np.zeros((k_dim,), dtype=data_type)
q_basis[[0], :] = v_0.T
projection = array @ v_0
alpha[0] = v_0.conj().T @ projection
projection = projection - alpha[0] * v_0
beta[0] = np.linalg.norm(projection)
for i in range(1, k_dim, 1):
v_p = q_basis[i - 1, :]
q_basis[[i], :] = projection.T / beta[i - 1]
projection = array @ q_basis[i, :]
alpha[i] = q_basis[i, :].conj().T @ projection
projection = projection - alpha[i] * q_basis[i, :] - beta[i - 1] * v_p
beta[i] = np.linalg.norm(projection)
# additional steps to increase accuracy
delta = q_basis[i, :].conj().T @ projection
projection -= delta * q_basis[i, :]
alpha[i] += delta
# exit if beta is less than numerical precision
if beta[i] < error:
k_dim = i
break
tridiagonal = (
np.diag(alpha[:k_dim], k=0)
+ np.diag(beta[: k_dim - 1], k=-1)
+ np.diag(beta[: k_dim - 1], k=1)
)
q_basis = q_basis[:k_dim]
return q_basis.T, tridiagonal
def lanczos_eig(array: Union[csr_matrix, np.ndarray], v_0: np.ndarray, k_dim: int):
"""
Finds the lowest ``k_dim`` eigenvalues and corresponding eigenvectors of a hermitian array
using Lanczos algorithm.
Args:
array : Array to diagonalize.
v_0 : Vector to initialise Lanczos iteration.
k_dim : Dimension of the krylov subspace.
Returns:
eigen_values : lowest ``k_dim`` Eigenvalues.
eigen_vectors_a : Eigenvectors in hilbert-space.
"""
q_basis, tridiagonal = lanczos_basis(array, v_0, k_dim)
eigen_values, eigen_vectors_t = np.linalg.eigh(tridiagonal)
eigen_vectors_a = q_basis @ eigen_vectors_t
return eigen_values, eigen_vectors_a
def lanczos_expm(
array: Union[csr_matrix, np.ndarray],
v_0: np.ndarray,
k_dim: int,
max_dt: float,
):
"""Calculates action of matrix exponential on the state using Lanczos algorithm.
Args:
array : Array to exponentiate.
v_0 : Initial state.
k_dim : Dimension of the krylov subspace.
max_dt : Maximum step size.
Returns:
v_dt : Action of matrix exponential on state.
"""
q_basis, tridiagonal = lanczos_basis(array, v_0, k_dim)
eigen_values, eigen_vectors_t = np.linalg.eigh(tridiagonal)
v_dt = (
q_basis
@ eigen_vectors_t
@ (np.exp(-1j * max_dt * eigen_values) * eigen_vectors_t[0, :])
)
return v_dt