Source code for netket.exact

# Copyright 2020, 2021 The NetKet Authors - All rights reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

from typing import Optional

import numpy as _np
from scipy.sparse.linalg import bicgstab as _bicgstab
from scipy.sparse.linalg import LinearOperator as _LinearOperator

from .operator import AbstractOperator as _AbstractOperator
from jax.experimental.sparse import JAXSparse as _JAXSparse


[docs] def lanczos_ed( operator: _AbstractOperator, *, k: int = 1, compute_eigenvectors: bool = False, matrix_free: bool = False, scipy_args: Optional[dict] = None, ): r"""Computes `first_n` smallest eigenvalues and, optionally, eigenvectors of a Hermitian operator using :meth:`scipy.sparse.linalg.eigsh`. Args: operator: NetKet operator to diagonalize. k: The number of eigenvalues to compute. compute_eigenvectors: Whether or not to return the eigenvectors of the operator. With ARPACK, not requiring the eigenvectors has almost no performance benefits. matrix_free: If true, matrix elements are computed on the fly. Otherwise, the operator is first converted to a sparse matrix. scipy_args: Additional keyword arguments passed to :meth:`scipy.sparse.linalg.eigvalsh`. See the Scipy documentation for further information. Returns: Either `w` or the tuple `(w, v)` depending on whether `compute_eigenvectors` is True. - w: Array containing the lowest `first_n` eigenvalues. - v: Array containing the eigenvectors as columns, such that`v[:, i]` corresponds to `w[i]`. Example: Test for 1D Ising chain with 8 sites. >>> import netket as nk >>> hi = nk.hilbert.Spin(s=1/2)**8 >>> hamiltonian = nk.operator.Ising(hi, h=1.0, graph=nk.graph.Chain(8)) >>> w = nk.exact.lanczos_ed(hamiltonian, k=3) >>> w array([-10.25166179, -10.05467898, -8.69093921]) """ from scipy.sparse.linalg import eigsh actual_scipy_args = {} if scipy_args: actual_scipy_args.update(scipy_args) actual_scipy_args["which"] = "SA" actual_scipy_args["k"] = k actual_scipy_args["return_eigenvectors"] = compute_eigenvectors if matrix_free: # wrap the operator.to_linear_operator() in a scipy.sparse.linalg.LinearOperator n = operator.hilbert.n_states A = _LinearOperator( (n, n), operator.to_linear_operator().__matmul__, dtype=operator.dtype, ) else: A = operator.to_sparse() if isinstance(A, _JAXSparse): # jax sparse arrays are not compatible with scipy eigsh. # wrap them in a scipy.sparse.linalg.LinearOperator A = _LinearOperator(A.shape, A.__matmul__, dtype=A.dtype) result = eigsh(A, **actual_scipy_args) if not compute_eigenvectors: # The sort order of eigenvalues returned by scipy changes based on # `return_eigenvalues`. Therefore we invert the order here so that the # smallest eigenvalue is still the first one. return result[::-1] else: return result
[docs] def full_ed(operator: _AbstractOperator, *, compute_eigenvectors: bool = False): """Computes all eigenvalues and, optionally, eigenvectors of a Hermitian operator by full diagonalization. Args: operator: NetKet operator to diagonalize. compute_eigenvectors: Whether or not to return the eigenvectors of the operator. Returns: Either `w` or the tuple `(w, v)` depending on whether `compute_eigenvectors` is True. Example: Test for 1D Ising chain with 8 sites. >>> import netket as nk >>> hi = nk.hilbert.Spin(s=1/2)**8 >>> hamiltonian = nk.operator.Ising(hi, h=1.0, graph=nk.graph.Chain(8)) >>> w = nk.exact.full_ed(hamiltonian) >>> w.shape (256,) """ from numpy.linalg import eigh, eigvalsh dense_op = operator.to_dense() if compute_eigenvectors: return eigh(dense_op) else: return eigvalsh(dense_op)
[docs] def steady_state(lindblad, *, sparse=True, method="ed", rho0=None, **kwargs): r"""Computes the numerically exact steady-state of a lindblad master equation. The computation is performed either through the exact diagonalization of the hermitian :math:`L^\dagger L` matrix, or by means of an iterative solver (bicgstabl) targeting the solution of the non-hermitian system :math:`L\rho = 0` and :math:`\mathrm{Tr}[\rho] = 1`. Note that for systems with 7 or more sites it is usually computationally impossible to build the full lindblad operator and therefore only `iterative` will work. Note that for systems with hilbert spaces with dimensions above 40k, tol should be set to a lower value if the steady state has non-trivial correlations. Args: lindblad: The lindbladian encoding the master equation. sparse: Whether to use sparse matrices (default: False for ed, True for iterative) method: 'ed' (exact diagonalization) or 'iterative' (iterative bicgstabl) rho0: starting density matrix for the iterative diagonalization (default: None) kwargs...: additional kwargs passed to bicgstabl For full docs please consult SciPy documentation at https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.bicgstab.html Keyword Args: maxiter: maximum number of iterations for the iterative solver (default: None) tol: The precision for the calculation (default: 1e-05) callback: User-supplied function to call after each iteration. It is called as callback(xk), where xk is the current solution vector Returns: The steady-state density matrix. """ M = lindblad.hilbert.physical.n_states if method == "ed": if not sparse: from numpy.linalg import eigh from warnings import warn warn( """For reasons unknown to me, using dense diagonalisation on this matrix results in very low precision of the resulting steady-state since the update to numpy 1.9. We suggest using sparse=True, however, if you wish not to, you have been warned. Your digits are your responsibility now. """, stacklevel=2, ) lind_mat = lindblad.to_dense() ldagl = lind_mat.T.conj() * lind_mat w, v = eigh(ldagl) else: from scipy.sparse.linalg import eigsh lind_mat = lindblad.to_sparse() ldagl = lind_mat.T.conj() * lind_mat w, v = eigsh(ldagl, which="SM", k=2) print("Minimum eigenvalue is: ", w[0]) rho = v[:, 0].reshape((M, M)) rho = rho / rho.trace() elif method == "iterative": # An extra row is added at the bottom of the therefore M^2+1 long array, # with the trace of the density matrix. This is needed to enforce the # trace-1 condition. L = lindblad.to_linear_operator(sparse=sparse, append_trace=True) # Initial density matrix ( + trace condition) Lrho_start = _np.zeros((M**2 + 1), dtype=L.dtype) if rho0 is None: Lrho_start[0] = 1.0 Lrho_start[-1] = 1.0 else: Lrho_start[:-1] = rho0.reshape(-1) Lrho_start[-1] = rho0.trace() # Target residual (everything 0 and trace 1) Lrho_target = _np.zeros((M**2 + 1), dtype=L.dtype) Lrho_target[-1] = 1.0 # Iterative solver print("Starting iterative solver...") res, info = _bicgstab(L, Lrho_target, x0=Lrho_start, **kwargs) rho = res[:-1].reshape((M, M)) if info == 0: print("Converged trace is ", rho.trace()) elif info > 0: print("Failed to converge after ", info, " ( trace is ", rho.trace(), " )") elif info < 0: print("An error occurred: ", info) else: raise ValueError("method must be 'ed' or 'iterative'") return rho