5. Eigen-Solving & 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 eigen-decompositions as each block can be decomposed individually. For example, imagine we want to supply the eigen-decomposition of the Heisenberg Hamiltonian to an Evolution object to perform long-time, 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 Z-spin, 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 in syncro mode.

  • 'slepc-nompi': like 'slepc', but performs computation in the single, main process.

  • 'AUTO' - choose a good backend, the default.

The possible partical decompositions are:

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 post-fix), get k=1 eigenstate, and target the ‘(s)mallest (a)lgebraic’ eigenvalue].

5.2.1. Interior eigen-solving

SLEPc is highly recommended for performing these using ‘shift-invert’. See the following functions:

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 rank-adaptively, which allows the efficient estimation of an operator’s rank, see estimate_rank().