5. EigenSolving & Other Linear Algebra¶
5.1. Dense / full decomposition¶
Currently full decompositions use numpy. They are as follows:
5.1.1. Automatic Exploitation of Symmetries via Blocking¶
Sometimes symmetries present themselves in an operator as, after
a certain permutation, diagonal blocks. This fact can be exploited
to speed up full eigendecompositions as each block can be
decomposed individually. For example, imagine we want to supply
the eigendecomposition of the Heisenberg Hamiltonian to an
Evolution
object to perform longtime, exact,
dynamics:
[1]:
import quimb as qu
H = qu.ham_heis(12)
[2]:
%%timeit
el, ev = qu.eigh(H)
8.81 s ± 566 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
This is already quite slow, and will not be sustainable for larger lengths. The Heisenberg Hamiltonian however has an abelian \(Z_2\) subgroup symmetry, conserved Zspin, that manifests itself in the computational basis. Thus if instead we specify autoblock=True
:
[3]:
%%timeit
el, ev = qu.eigh(H, autoblock=True)
372 ms ± 192 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
We can see a drastic speedup (helped by the fact that the whole algorithm uses numba
).
To work more fundamentally with the indentified blocks see the
compute_blocks()
function.
5.2. Partial decomposition¶
Partial decompositions are mostly just specified by supplying the k
kwarg to the above functions. These also take a backend
argument which can be one of:
'scipy'
: is generally reliable'numpy'
: can be faster for small or dense problems'lobpcg'
: useful for fast, low accruacy generalized eigenproblems (like periodic DMRG)'slepc'
: Usually the fastest for large problems, with many options. Will either spawn MPI workers or should be used insyncro
mode.'slepcnompi'
: like'slepc'
, but performs computation in the single, main process.'AUTO'
 choose a good backend, the default.
The possible partical decompositions are:
eig()
with(k > 0)
eigh()
with(k > 0)
eigvals()
with(k > 0)
eigvalsh()
with(k > 0)
eigenvectors()
with(k > 0)
eigvecs()
with(k > 0)
eigvecsh()
with(k > 0)
So for example the groundstate()
function
for a Hamiltonian H
is an alias to:
[4]:
psi = qu.eigvecsh(H, k=1, which='sa')
[find eigenvectors, Hermitian operator (h
postfix), get k=1
eigenstate,
and target the ‘(s)mallest (a)lgebraic’ eigenvalue].
5.2.1. Interior eigensolving¶
SLEPc is highly recommended for performing these using ‘shiftinvert’. See the following functions:
eigh(..., k=k, sigma=x)
withk > 0
etc., or
With the last three allowing the specification of a window relative to the total spectrum of the operator.
5.2.2. Fast Randomized Linear Algebra¶
quimb
has an implementation of a fast randomized SVD  rsvd()

that can be significantly quicker than svd()
or svds()
,
especially for large k
. This might be useful for e.g. tensor network linear operator decompositions.
It can perform the SVD rankadaptively, which allows the efficient estimation of an operator’s rank,
see estimate_rank()
.